The fragility of thin discs in galaxies – I. Building tailored N-body galaxy models

ABSTRACT

1 INTRODUCTION

About 70 per cent of the observed galaxies in the local Universe host (either thin or thick) kinematically cold disc-like structures in their stellar and gas distributions, regardless of the classification (Ilbert et al. 2006; Weinmann et al. 2006; Choi, Park & Vogeley 2007; Park et al. 2007; van den Bosch et al. 2007; Ledo et al. 2010). The high occurrence of disc galaxies and the paucity of past interaction signatures – such as gas bridges or distinct rotating discs (e.g. Corsini 2014; Mazzilli Ciraulo et al. 2021) or extended stellar shells (e.g. Hernquist & Quinn 1988; Romanowsky et al. 2012; Onaka et al. 2018; Pop et al. 2018; Bílek et al. 2022) – may appear at odds with the current hierarchical model of structure growth, in which minor and major mergers contribute significantly to the mass assembly of galaxies. Indeed, the effect of past interactions can dim in time, making necessary in-depth studies of the galaxy stellar populations to unveil their assembly history (Davison et al. 2021; Mazzilli Ciraulo et al. 2021).

Kinematically cold thin discs are fragile structures subject to morphological transformations during interactions with their environment (Genel et al. 2014; Vogelsberger et al. 2014a, b; Sijacki et al. 2015; Joshi et al. 2020; Galán-de Anta et al. 2022), including galaxy mergers. For this reason, discs have been proposed as natural clocks of their last merger event (Toomre 1977; Barnes & Hernquist 1992; Hammer et al. 2009; Taranu, Dubinski & Yee 2013; Deeley et al. 2017).

It has, however, been recently observed that thin and kinematically cold structures could survive merger events spanning a wide range of mass ratios (from 1:3 to <1:10; Abadi et al. 2003; Robertson et al. 2006; Purcell, Kazantzidis & Bullock 2009; Lotz et al. 2010; Moster et al. 2010). In some cases, disc galaxies might even survive major encounters, leaving a prominent disc component in the merger remnant (e.g. Springel 2005; Naab, Jesseit & Burkert 2006; Robertson et al. 2006; Governato et al. 2007; Capelo et al. 2015). Athanassoula et al. (2016), Sparre & Springel (2017), and Peschken, Łokas & Athanassoula (2020) show that merger remnants of wet major mergers can induce the formation of galactic discs.

Alternative kinematically cold tracers of the last merger event are nuclear stellar discs (NSDs), originally unveiled in Hubble Space Telescope images (Jaffe et al. 1994; van den Bosch et al. 1994), where they appeared as razor-thin disc structures of a few hundred pc across lying at the centre of galaxies (Pizzella et al. 2002). Ledo et al. (2010) presented a catalogue of NSDs in a wide variety of early-type galaxies, getting a rough estimation for the number of NSDs residing in galactic nuclei to be about 20 per cent, making them a structure commonly present in the Universe. However, the properties of the stellar populations (i.e. age, metallicity, and star formation time-scale) have been derived only for a few NSDs (Corsini et al. 2016; Sarzi et al. 2016).

Ledo et al. (2010) and Sarzi, Ledo & Dotti (2015) tested for the first time the fragility of NSDs against mergers, by performing a set of pure N-body simulations consisting of a NSD, a stellar halo, and a supermassive black hole (SMBH) in interaction with a secondary SMBH. In particular, Sarzi et al. (2015) explored a broad region of the merger parameter space [circular orbits with different inclinations, impact parameters, and masses, including major (1:1), intermediate (1:5), and minor (1:10) mergers] and showed that these discs cannot survive any major encounter but can withstand minor ones.

These first studies relied on rather idealized representations of the inner regions of the NSD host galaxies. More realistic studies modelling the whole stellar distribution of NSD hosts are needed in order to confirm the conclusions presented in Sarzi et al. (2015), allowing to consistently gauge the effect on the large-scale galactic disc and on the NSD.

To isolate the role of mergers and other physical processes in galaxy evolution, it is common to construct the initial conditions for the simulations as stationary equilibrium configurations. The methods for constructing these models can be grouped into a few categories. Jeans equations offer the fastest and least demanding approach, usually relying only on the first two moments of the stellar distribution function (DF; e.g. Cappellari 2008; Mamon, Biviano & Boué 2013), although the use of higher order moments (e.g. Łokas & Mamon 2003; Richardson & Fairbairn 2013; Read & Steger 2017) may help to lift the so-called mass–anisotropy degeneracy inherent to this approach (Dejonghe & Merritt 1992). It should be stressed that, while being fast and simple to produce, N-body realizations of galaxy models based on Jeans equations are not guaranteed to be in full equilibrium (Kazantzidis, Magorrian & Moore 2004).

At the other end of the spectrum are made-to-measure N-body codes capable of guiding the model to meet specific observational constraints (Syer & Tremaine 1996; De Lorenzi et al. 2007; Yurin & Springel 2014), which are very flexible but also expensive to run. The Schwarzschild (1979) orbit-superposition method (in numerous implementations) has also been used to construct flexible models of observed galaxies (e.g. Cretton et al. 1999; Gebhardt et al. 2003; van den Bosch et al. 2008) and generate the initial conditions for N-body simulations (e.g. Vasiliev & Athanassoula 2015). Finally, self-consistent models based on DFs have been used both in the theoretical and simulation context (e.g. Kuijken & Dubinski 1995; Debattista & Sellwood 2000) and observational applications (e.g. Widrow & Dubinski 2005; Piffl, Penoyre & Binney 2015; Taranu et al. 2017; Bienaymé, Leca & Robin 2018). Nevertheless, even the sophisticated N-body and DF-based approaches often do not produce systems in exact equilibrium, and need to be followed by an initial ‘relaxation’ stage before running merger simulations (e.g. fig. 2 and section 3.2.1 of Garavito-Camargo et al. 2019). On the other hand, N-body simulations are also subjected to numerical heating that leads to an artificial expansion of the flattened components and to a randomization of the circular orbits that also need to be taken in consideration (Ludlow et al. 2019, 2021; Wilkinson et al. 2023, hereafter L19, L21, and W22, respectively).

In this paper, we present a new method for constructing N-body galaxy models based on DFs and tailored to observed kinematics. We apply our procedure to the NSD host FCC 170 (NGC 1381), for which Pinna et al. (2019, hereafter P19) constructed detailed two-dimensional (2D) maps of stellar ages, abundances, and stellar kinematics along the line of sight (LOS). We demonstrate that these models reproduce well the observed data and are in near-perfect equilibrium, eliminating the need to perform an initial relaxation stage. We also show how numerical heating affects the flattening of both the thin disc and the NSD of our N-body FCC 170 model. With the mass resolution achieved in our tests, this numerical artefact is well under control for both the kpc-scale and nuclear discs of FCC 170, introducing an artificial 2 and 20 per cent thickening of such discs, respectively, after 10 Gyr of isolated evolution. The current study is the first step towards the exploration of the fragility of thin-disc structures (kpc scale and nuclear scale) against different galactic mergers, which will be the topic of a follow-up paper currently in preparation (hereafter Paper II).

The paper is structured as follows: In Section 2, we give a brief description of the FCC 170 data obtained from the Multi Unit Spectroscopic Explorer (MUSE) and its observed kinematics. Section 3 explains how we build the pure N-body models based on DFs. In Section 4, we describe the set-up of the code used to evolve the N-body model, along with the discussion on the stability and evolution of the N-body galaxy in isolation. In Section 5, we give our conclusions.

2 SPECTROSCOPIC OBSERVATIONS AND STELLAR KINEMATICS OF FCC 170

We selected the edge-on galaxy FCC 170 (NGC 1381), hosted in the Fornax cluster, making use of the data obtained by the Fornax3D survey (Sarzi et al. 2018; Iodice et al. 2019) with the MUSE integral-field unit installed at the Very Large Telescope. The MUSE datacubes were taken using the wide-field mode, providing a spatial sampling of 0.2 arcsec × 0.2 arcsec on a 1 arcmin × 1 arcmin field of view. The wavelength range of the MUSE datacubes is enclosed between 4650 and 9300 Å, with a spectral sampling of 1.25 Å pixel⁻¹ and an average spectral resolution FWHM_int = 2.8 Å. The spatial scaling of the MUSE images is |$0.2{\,\mathrm{ arcsec}}\, {\rm px}^{-1}$| and we assume a distance towards our target galaxy FCC 170 of 21.9 Mpc, in concordance to P19. The actual field of view consists of a mosaic of two pointings: a central pointing that covers the inner regions of the galaxy and an offset pointing that covers the outer disc and halo region of the galaxy. The central and offset pointings have integration times of 60 and 90 min, respectively, due to different signal-to-noise ratios (SNRs), in order to reach the same limiting surface brightness of |$\mu _{\rm B}=25 \rm \, mag\, arcsec^{-2}$|⁠. The acquisition and reduction of the datacubes are extensively described in Sarzi et al. (2018) and Iodice et al. (2019).

The MUSE pointings have been reduced using their own dedicated pipeline (Weilbacher et al. 2012; Weilbacher, Streicher & Palsa 2016) within the environment ESOREFLEX (Freudling et al. 2013), as described in Sarzi et al. (2018) and Iodice et al. (2019). In the reduction, they took care of sky subtraction, telluric correction, and flux calibration, either relative and absolute. Generally, single pointings are aligned throughout reference stars and later combined to produce the final MUSE mosaics.

P19 and Iodice et al. (2019) present 2D maps of the stellar kinematics for the edge-on lenticular galaxy FCC 170, showing also maps of the stellar populations including metallicities and ages. The main differences between these studies are the required target SNR of the Voronoi bin [P19 impose a target |$\rm SNR=40$| for FCC 170 and a minimum |$\rm SNR=1$| per spaxel, whereas Iodice et al. (2019) require a minimum |$\rm SNR=3$| per spaxel] and the methodology for deriving stellar populations based on two different methods: spectral fitting in P19 and measurement of line-strength indices in Iodice et al. (2019). P19 obtained kinematics 2D maps of FCC 170 by using the Penalized Pixel-Fitting algorithm (Cappellari & Emsellem 2004; Cappellari 2017), fitting the stellar spectra of the galaxy combined with a set of stellar population templates. The galaxy is initially binned in a Voronoi tessellation grid (Cappellari & Copin 2003) accounting for a required minimum SNR in each spaxel. After binning the flux image of the galaxy, P19 apply a set of single-stellar population templates from the MILES stellar library based on BaSTI isochrones (described in Vazdekis et al. 2015) to fit each Voronoi-binned spectrum and recover the kinematics of the stars, using a standard Gauss–Hermite expansion to parametrize the stellar LOS velocity distribution (LOSVD; Gerhard 1993; van der Marel & Franx 1993).

3 BUILDING N-BODY MODELS TAILORED TO FCC 170 OBSERVATIONS

3.1 Method

We construct initial conditions for our simulations resembling the actual FCC 170 galaxy, using the agama stellar-dynamical framework (Vasiliev 2019). It provides a wide range of tools for various tasks in stellar dynamics, in particular, several methods for constructing multicomponent equilibrium galaxy models and computing their observational properties. We use the iterative self-consistent modelling approach described in section 6.2 of that paper. Below, we briefly summarize its features and our fitting strategy.

and finally recompute the total potential Φ from the Poisson equation, |$\nabla^2\Phi = 4\pi G\sum_c \rho_c,$| where G is the gravitational constant. The whole process is then repeated several times, until the changes in the potential are negligible (≲1 per cent; Binney 2014).

We use different classes of DFs for disc and spheroidal galactic components, which are fully described in section 4 of Vasiliev (2019). The three discs (thin, thick, and NSD) are represented by QuasiIsothermal models with the following free parameters: scale length, scale height, central value of the radial velocity dispersion, and scale radius of its exponential decay. The radial and vertical density profiles generated by this DF are close to exponential and sech² (van der Kruit & Freeman 2011), respectively, as it is typical for stellar discs. The bulge DF belongs to the DoublePowerLaw family (Posti et al. 2015), which has as many as nine free parameters: power-law indices of asymptotic behaviour at small and large radii, steepness of the transition between the two regimes, corresponding characteristic spatial scale, amount of rotation, and four dimensionless coefficients describing the radial and vertical velocity anisotropy (the latter also implicitly determines the flattening of the density profile) at small and large radii. Finally, for the dark matter (DM) halo we do not have any kinematic constraints, so we choose a simple QuasiSpherical DF determined by its density profile, which is taken to follow a Navarro, Frenk & White (1996) model with an adjustable scale length and a truncation radius 10 times larger (it is essentially unconstrained by stellar kinematics). The masses of all five components are also free parameters, bringing their total number to 26.

We also add a central SMBH with a mass of 3 × 10⁷ M_⊙. It is slightly larger than the value adopted by Poci et al. (2021), but consistent with the SMBH mass–spheroid mass relation as given by Kormendy & Ho (2013). This value is fixed throughout the fitting process, since the expected kinematic signature of such an SMBH is confined to a very small spatial region (the sphere of influence), which is not resolved by the non-adaptive-optics integral-field unit observations.

3.2 Fitting models to observations

To evaluate the likelihood of a model with the given set of parameters, we need to compute the LOSVDs in each of the ∼8000 Voronoi bins, convert them to Gauss–Hermite moments, and then compare to the observational kinematic maps and to the surface brightness map. The parameters of the Gauss–Hermite expansion are the centre value v₀ and width σ of the base Gaussian (which are close to but not identical to the mean velocity and its dispersion, respectively), and two higher order moments h₃ and h₄ (which are related to skewness and kurtosis and quantify the shape of the LOSVD). Although agama can compute these quantities very accurately by integrating the DF over the LOS and the two sky-plane velocity components, this is very computationally expensive. Instead, we use the following strategy: the DFs of a fiducial model |$f_c^\mathrm{(fid)}$| are sampled into equal-mass N-body snapshots, and the positions/velocities |$\lbrace \boldsymbol{x}_k, \boldsymbol{v}_k\rbrace$| of these particles are used to compute the kinematic maps in the Voronoi bins (essentially replacing the deterministic multidimensional integration by Monte Carlo sampling). The rejection sampling procedure itself requires a few times larger number of DF evaluations than the number of output samples, and inevitably introduces some discreteness noise due to a finite number of particles. To mitigate both factors, in the subsequent evaluation of model likelihoods for different choices of parameters, we use the same set of sample points in the six-dimensional phase space, but reweigh their contribution to the kinematic maps by the ratio of the DF values in the current model to those in the fiducial one, |$f_c^\mathrm{(curr)}\big (\boldsymbol{J}(\boldsymbol{x}_k, \boldsymbol{v}_k;\,\,\Phi ^\mathrm{(curr)}) \big)/f_c^\mathrm{(fid)}\big (\boldsymbol{J}(\boldsymbol{x}_k, \boldsymbol{v}_k;\,\,\Phi ^\mathrm{(fid)}) \big)$|⁠, where the conversion from the phase space to action space also depends on the corresponding model potential. Finally, in the early stages of the parameter space exploration, we use a more approximate but ∼5× faster interpolation scheme for action computation, and only switch to the more accurate default (non-interpolated) approach once near the maximum-likelihood point. In the end, each model construction and likelihood evaluation on a 32-core workstation takes only |$\mathcal {O}(10)$| s with interpolated actions, or about 1 min in the default approach, making it possible to explore thousands of points in the parameter space.

We use a standard simplex (Nelder–Mead) method for minimizing the objective function, but find it to be struggling with the high dimensionality of the parameter space and inevitable noisiness of the likelihood function stemming from various numerical effects. To reduce the chance of being trapped in a local minimum, we restarted the optimization procedure multiple times from various initial points; the fits usually ended up in the same region, but the rate of convergence was fairly slow. Evidently, a more efficient optimization scheme is needed to make this fitting method practical.

Fig. 1 compares the kinematic maps of the ‘best-fitting’ (in a loosely defined sense) model with observations. The surface brightness map is well reproduced, even though the stellar density profile is not fitted independently, but comes out as an integral of the DF over velocity. The kinematic maps also qualitatively match the data, but not without apparent deviations in all four Gauss–Hermite moments. The most noticeable discrepancy occurs in the bulge region outside the NSD, up to a few arcseconds (corresponding to a few hundred pc): the model fails to reproduce a nearly flat plateau in |$v_0\approx \pm 60\, \mathrm{km\, s}^{-1}$| between 1 and 5 arcsec and a subsequent rather sharp rise of v₀ to more than twice this value in the disc region. Likewise, the central region of high velocity dispersion extends much farther in the model than in the observations. This suggests that the parametric DF family used in our fits does not have enough flexibility to fully represent the kinematic structure of the bulge and its transition to the disc. There are also noticeable deviations in the higher Gauss–Hermite moments in the outer parts of the disc, probably due to inadequacy of the exponential decline of the velocity dispersion with radius prescribed by the disc DF. Nevertheless, the DF-based models are suitable for the main purpose of this paper – the creation of equilibrium N-body models qualitatively resembling the observed galaxy and having a physically motivated multicomponent structure. A list of the main parameters we use in agama to set up our FCC 170 N-body model is tabulated in Table B1.

Figure 1.

Comparison of the observational data (left-hand column) with the DF-based model (central column) and (observed minus model) residuals normalized by the observational uncertainties (right-hand column) along with χ² of each residual map. From top to bottom, the rows show the maps of the stellar surface luminosity density and the four parameters describing the LOSVD in terms of Gauss–Hermite moments. The solid lines correspond to the radial profiles extracted along the galaxy major axis in the region bracketed by the horizontal dashed lines.

Open in new tabDownload slide

We also explored an alternative approach for the construction of dynamical models tailored to observations, namely the Schwarzschild (1979) orbit-superposition code forstand (also included in agama). A comparison between the orbit-based and DF-based models is provided in Appendix  A.

3.3 Comparison with the photometric decomposition

The NSD in FCC 170 is a thin nuclear structure at scales of a few tens of pc, as shown by Ledo et al. (2010), who catalogued a variety of NSDs with scale lengths of about ∼100 pc, including the one at the centre of FCC 170. The image of the NSD, along with its magnitude, ellipticity, and A4 parameters are shown in the right-hand panels of fig. A3 in Ledo et al. (2010). Morelli et al. (in preparation) produce an in-depth modelling of the surface brightness distribution of this NSD. They have analysed a Wide Advanced Camera Survey image of FCC 170 obtained from the Hubble Legacy Archive with the Scorza and Bender surface brightness decomposition (Scorza & van den Bosch 1998; Morelli et al. 2004) as implemented by Corsini et al. (2016). This method is based on the assumption that the isophotal discyness is the result of the superposition of a spheroidal component (which is either an elliptical galaxy or a bulge) and an inclined infinitesimally thin exponential disc. The two components are assumed to have both perfectly elliptical isophotes with constant but different ellipticities.

In Fig. 2, we present the contribution D/T of the NSD to the total stellar light measured in the central region along the major axis of FCC 170 (red solid line; Morelli et al., in preparation) and the total particle mass of its equivalent N-body model (blue dashed line). We also include the 1σ errors of the photometric decomposition with the shaded grey region. The selection of the bin size for the N-body model is done by adopting the characteristic scale height of the NSD (Table B1) as the constant bin height (0.21 arcsec), while the bin width is given by a logarithmic radius of constant step Δln (R/arcsec) = 0.14 to account for quick variations of the D/T at very short scales in contrast with a smoother slope at larger scales. The D/T radial profiles from the photometric decomposition and N-body model match each other rather well and show how good are our N-body simulations in reproducing such a very peculiar galaxy as FCC 170. The small (<5 per cent) discrepancies in the central region of the galaxy (r < 2 arcsec) are a consequence of the fact that the N-body model of the NSD is mainly based on kinematics, whereas the observational properties of the NSD result from photometry only.

Figure 2.

Radial profile of the disc-to-total fraction of the NSD in FCC 170 along the major axis of the galaxy from the photometric decomposition by Morelli et al. (in preparation) (black solid line) and our N-body model (red dash–dotted line) along with the 1σ errors for the photometric decomposition.

Open in new tabDownload slide

4 CHECKING THE STABILITY OF THIN DISCS IN ISOLATION

We check the stability and passive evolution of our model in isolation in order to quantify how thin-disc structures are affected by numerical heating or local dynamical instabilities. Particularly, disc galaxies could form bars (e.g. Abbott et al. 2017; Zana et al. 2018; Patsis & Athanassoula 2019) and spiral arms (e.g. Díaz-García et al. 2019; Sellwood & Carlberg 2019; Martinez-Medina, Pérez-Villegas & Peimbert 2022) with the formation of the latter being often induced by galactic bars (Garma-Oehmichen et al. 2021) or pseudo-bulges (Yu et al. 2022). Moreover, the stability test in isolation provides a good reference to assess how much the thin stellar discs are perturbed by a merging event, when we will consider the case of two interacting galaxies.

4.1 Creating and running the N-body simulation

After the best-fitting parameters have been found, we construct an N-body realization of the model by sampling particle positions and velocities from the DF. We use particles of the same mass for all stellar components, but retain the information about the component they belong to, and use a larger particle mass for the DM halo. The SMBH is represented by a single massive particle with a very small softening length.

In order to evolve our N-body model in isolation, we use the code gizmo (Hopkins 2015), which has an architecture that accounts for many different physical processes, including magneto-hydrodynamics and the black hole and supernova feedback. gizmo offers two different gravitational solvers – hybrid tree or the tree–particle mesh scheme – both offering automatic adaptivity of the gravitational resolution in structures that are collapsing or expanding. gizmo also permits to fix the softening lengths of each particle type independently, which facilitates to account for 2D scattering at different spatial scales.

Softening lengths have been chosen to be small enough to properly resolve the vertical structure of each component. This can normally be addressed by selecting a length at most half of the characteristic scale radius or scale height of a given component. We have chosen the softening lengths of each particle type to be: ε_DM = 50 pc, ε_bulge = ε_thin-disc = ε_thick-disc = 10 pc, ε_NSD = 5 pc, and ε_SMBH = 1 pc. We run the simulation by using 480 processor cores spread between 15 nodes, taking about 200 h of wall-clock time to reach a total integration time of 10 Gyr and producing 200 snapshots equally spaced in time by 50 Myr. Our main run hosts N_⋆ = 3.71 × 10⁶ stellar particles with M_⋆ = 7.32 × 10³ M_⊙ per particle and N_DM = 1.37 × 10⁷ DM particles with M_DM = 2.48 × 10⁵ M_⊙ per particle. To further test numerical heating on thin discs (see Section 4.2), we also run three additional N-body runs of FCC 170 by decreasing the number of DM particles by 10, 5, and 2 with respect to our main model.

In Fig. 3, we show the evolution of the model by plotting the surface mass density of stellar particles at different times. The bottom panels show the face-on view of the whole stellar structure and the top panels show the edge-on view. Investigating the stellar particles, we observe that the model remains in equilibrium throughout the whole simulation without forming a bar, spiral arms, or any dynamical instability. At 10 Gyr, the disc becomes slightly thicker, increasing the vertical size and partially losing its initial flattening, as a consequence of numerical heating. We explore in more detail the thickening of both the thin galactic disc and NSD in Section 4.2.

Figure 3.

Surface mass density maps of the FCC 170 N-body model at different times. Each column shows the edge-on (top panel) and face-on (bottom panel) views of the stellar component (i.e. NSD, bulge, and thin and thick discs) at a given time. Each panel is 40 kpc wide.

Open in new tabDownload slide

Last, when evolving the model in isolation, we find that the kinematics at t = 10 Gyr remains similar to that shown in Fig. 1, confirming that the galaxy remains stable for the whole run.

4.2 Numerical heating of thin-disc components

L19, L21, and W22 demonstrated that different mass components in N-body simulations lead to a numerical expansion of the less massive structural component and, consequently, to a flow of kinetic energy from the more massive components to the less massive ones. This is a numerical effect in which less massive components reduce their numerical two-body relaxation time (already significantly shorter than the physical one due to the limited number of particles used) by an additional μ = m₁/m₂ factor (with m₁ ≥ m₂), where m₁ and m₂ are the masses of particle type 1 and 2, respectively. As a consequence, we expect the natural size of every disc structure to be increased in both radial scale and height, as reported by L19, L21, and W22, as the mass of the stellar particles is way smaller than the mass of the DM particles. Increasing the number of particles in the other components with respect to the NSD with the aim to significantly reduce the μ fraction will alleviate the issue. W22 show that, when the number of particles of the DM halo rises up to ∼10⁷ for a prefixed individual stellar particle mass, numerical heating on thin-disc structures gets significantly reduced, becoming negligible on an ∼10 Gyr time-scale.

In Appendix  B, we show that numerical heating on the NSD is a direct consequence of the presence of stellar and DM particles with different masses. We replace the DM halo and the other stellar components except the NSD itself by a static analytical Hernquist (1990) potential with parameters chosen to approximate their combined gravitational force. The NSD and the central SMBH are then evolved as an N-body system in this external potential, and in this case the NSD does not experience any significant expansion or distortion.

To investigate the effects that the other particle types have on the NSD, we rebuild the FCC 170 model using the same set-up as in Table B1 but decreasing the number of DM particles by a factor of 10, 5, and 2. By running all these models in isolation, we observed how segregation depends on the differences in mass between the stellar and DM components. In Fig. 4, we show the evolution of thin-disc particles through time for four runs of the FCC 170 N-body model with 10⁶, 2 × 10⁶, 5 × 10⁶, and 10⁷ (with the latter the main model) DM particles. Fig. 5 shows the NSD particles for the same snapshots as in Fig. 4. The effect of numerical heating correlates with the mass of the DM particles getting significantly reduced when N_DM ≥ 10⁷, in agreement with the results found by W22 (see their figs 1, 2, and 3). Indeed, the variations on the NSD are rather larger than those for the thin-disc component, as the characteristic scale length and mass of the NSD are about 10 times smaller than those of the thin disc. Consequently, as the relaxation time is proportional to both the size and number of particles, we should expect the relaxation time of the NSD to be shorter than that of the thin disc, and numerical heating would reduce this relaxation time even more.

Figure 4.

Surface mass density maps of the thin-disc particles of the FCC 170 N-body model seen edge-on at different times. Each column corresponds to a run with a different particle number and mass of DM particles, as given in the first and second rows, respectively. Each panel is 20 kpc wide.

Open in new tabDownload slide

Figure 5.

Same as Fig. 4, but for the NSD particles.

Open in new tabDownload slide

If these results show that numerical heating can significantly affect smaller thin-disc structures, such artificial changes may not necessarily impede our planned investigation of the impact of intermediate/minor mergers on NSDs, as long as a more precise evaluation of numerical disc heating and of its evolution in time is at hand.

4.3 Quantifying numerical heating by mock imaging

To quantify the time changes in shape and surface brightness of our N-body model of FCC 170, we measure the profiles of ellipticity and surface mass density of the thin disc and the NSD, using the package photutils (Bradley et al. 2020). We fit a set of ellipses to edge-on projections of the thin disc and the NSD separately, excluding the other components.

In Fig. 6, we show the photometric analysis of our N-body model at different times plotting the surface mass density and ellipticity of all the stellar particles (blue lines), thin-disc particles (red lines), and NSD particles (black lines) along the major axis in edge-on projection. We also indicate the limiting values of ellipticity at large radii at t = 0 (dashed horizontal lines), in order to illustrate how this value decreases with time for each component. The dot–dashed lines show the values of surface mass density and ellipticity for the N-body model with 10⁶ DM particles. We observe that the surface mass density of the galaxy evolves very little with time for the model with 10⁷ DM particles, contrary to the lower resolution model. All stellar components gradually expand in radius and thickness and decrease their surface mass density, as a consequence of numerical heating. To quantify the impact of heating on the disc thickness, we analyse the ellipticity of the stellar distribution. In the model with 10⁷ DM particles, the ellipticity of the entire stellar distribution decreases by just 2 per cent over the course of 10 Gyr, with the structure of the kpc-scale thin disc changing even less, by 1 per cent at most. The impact on the NSD particles is more pronounced, with the NSD eventually becoming 22 per cent thicker after 10 Gyr. On shorter time-scales, however, the impact of numerical heating is more contained, with an artificial thickening of ∼5, 6, and 10 per cent after 1, 2, and 5 Gyr, respectively. Fig. 6 also further demonstrates how increasing the number of DM particles is critical in controlling the impact of numerical heating on smaller central structures, as adopting only 10⁶ DM particles leads not only to a much thicker but also considerably more extended NSD structure.

Figure 6.

Photometric analysis of the FCC 170 N-body model at different times. Left-hand column: edge-on view of the model along with some of the elliptical apertures adopted for the photometric analysis. Central column: surface mass density in arbitrary units along the major axis of the model for all stellar particles (blue lines), thin disc (red lines), and NSD (black lines). Right-hand column: ellipticity radial profiles along the major axis of the model for all stellar particles (blue lines), thin disc (red lines), and NSD (black lines). Dashed horizontal lines mark the value at which the ellipticity converges for each component at t = 0. In the central and right-hand columns, the solid and dot–dashed lines refer to the N-body models with 10⁷ and 10⁶ DM particles, respectively.

Open in new tabDownload slide

In Paper II, we will analyse how mergers affect the thickening of thin discs by comparing the time changes in ellipticity of the merging model with the N-body FCC 170 model evolved in isolation. For instance, if at the end of a 1:4 merger event lasting ∼2 Gyr (counting since the first pericentre passage of the secondary galaxy), we observe a decrease in the NSD ellipticity by 25 per cent, then we may conclude that the merger event affected the NSD beyond what could be artificially produced by numerical heating.

5 CONCLUSIONS

We presented a pure N-body model tailored to represent the edge-on lenticular galaxy FCC 170 based on integral-field spectroscopic data. It is constructed using a new approach for building self-consistent equilibrium models specified by DFs. The DF parameters are optimized to match the observed kinematics, including not only v and σ but also the high-order moments of the LOSVD. We track the evolution of the thin disc and NSD in the subsequent N-body simulation. We demonstrate that the model remains in equilibrium when evolved in isolation, and its properties do not change significantly over 10 Gyr. The only exception is a moderate increase in thickness of the NSD, caused by numerical heating from the much heavier DM halo particles. We kept this artificial thickening under control by using a sufficiently high number of halo particles.

In Paper II, we will explore the fragility of thin-disc structures in galactic encounters between our N-body model and a secondary spheroidal galaxy in a scenario that is consistent with the mergers observed in cosmological simulations. In this context, our N-body model offers an advantage compared with cosmological simulations (along with a higher resolution than various zoom-in state-of-the-art simulations), as we can control the conditions of the merger such as the initial distance between galaxies and the type of orbit.

ACKNOWLEDGEMENTS

We thank the reviewer, Curtis Struck, for the constructive feedback. This work was performed on the OzSTAR national facility at Swinburne University of Technology. The OzSTAR programme receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government. We are grateful for use of the computing resources from the Northern Ireland High Performance Computing (NI-HPC) service funded by the Engineering and Physical Sciences Research Council (EPSRC) (EP/T022175). Lorenzo Posti acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement no. 834148). Enrico Maria Corsini acknowledges support by Padua University grants Dotazione Ordinaria Ricerca (DOR) 2019–2022 and by Italian Ministry for Education University and Research (MIUR) grant Progetto di Ricerca di Interesse Nazionale (PRIN) 2017 20173ML3WW-001.

DATA AVAILABILITY

The data underlying this article can be made available upon request. The model results can be reproduced using publicly available codes.

REFERENCES

APPENDIX A: SCHWARZSCHILD MODEL OF FCC170

In addition to the DF-based models described in Section 3, we also constructed Schwarzschild (1979) orbit-superposition models of the same galaxy, using the forstand code (Vasiliev & Valluri 2020). These models are constrained by the observed kinematic and photometric properties of the galaxy in the form of maps of the Gauss–Hermite moments of the stellar LOSVD and multi-Gaussian parametrization of the surface-brightness distribution. They are similar to the ones presented by Poci et al. (2021), who relied on another implementation of this method from Bosch et al. (2008). Fig. A1 illustrates that these models fit the data largely down to the level of measurement uncertainty, by virtue of having an enormously larger number of (hidden) free parameters – orbit weights. We use 40 000 orbits to fit roughly the same number of observational constraints. The disadvantage of this method is that it does not discriminate a priori between different stellar components (i.e. the bulge, NSD, and thin and thick discs), which is critical for our application. Although it is possible to assign orbits to these components according to their properties (e.g. circularity), as shown in fig. 5 of Poci et al. (2021), this is beyond the scope of this study. Fig. A2 compares the internal kinematic properties of the DF-based and orbit-superposition models, demonstrating a satisfactory agreement. As our goal in this paper is not to find the best match to the observed galaxy, but rather to construct a suitable approximation with well-defined physical components, we opted to use the multicomponent DF-based models in the rest of our analysis.

Figure A1.

Same as Fig. 1, but for the Schwarzschild (1979) orbit-superposition model of FCC 170.

Open in new tabDownload slide

Figure A2.

Comparison of the internal kinematic properties of the DF-based (left-hand panel) and Schwarzschild (1979) orbit-superposition (right-hand panel) models of FCC 170. The radial profiles of the streaming velocity |$\overline{v_\phi }$| (green line), radial velocity dispersion σ_R (red line), and vertical velocity dispersion σ_z are shown together with the curve of circular velocity |$v_{\rm circ}\equiv \sqrt{R\, \mathrm{d}\Phi /\mathrm{d}R}$| (black line) of the total potential. Both methods produce qualitatively similar profiles, although the orbit-superposition model better matches the bump in the streaming velocity observed at ∼20 arcsec and its decline at large radii.

Open in new tabDownload slide

APPENDIX B: NUMERICAL HEATING OF THE NSD IN AN EXTERNAL ANALYTICAL POTENTIAL

To illustrate the effect of other galaxy components (primarily the DM halo) on the evolution of the NSD, we conducted the following experiment. We replaced all other components except the NSD and the SMBH by a static Hernquist (1990) potential, whose mass and scale radius are chosen to approximate the circular-velocity curve (equivalently, the cumulative mass profile) of these components combined (primarily the bulge, which is the dominant contribution in the spatial region occupied by the NSD). By doing so, we eliminate the numerical relaxation caused by other components. We then ran an N-body simulation of the NSD and the SMBH in gizmo with such external potential.

Fig. B1 shows the snapshots of the self-gravitating NSD component embedded in the external Hernquist potential with a black hole at three different times. As we can observe, the NSD remains unperturbed for the whole run. This suggests that the puffing up of the small razor-thin disc is a consequence of the much larger particle masses in the DM halo, as shown in L19 and L21, and it only occurs when representing the other components as a live N-body system. In the right-hand panel, we also plot the circular velocity of the NSD, of the analytical Hernquist profile, and of the bulge, as this last component dominates the central potential up to 1 kpc. Despite the discrepancies at short kpc scales between the analytical Hernquist and the bulge component, the NSD remains in perfect equilibrium for the entire run and any other choice of parameters for the Hernquist deviates the NSD out from the equilibrium stage.

Figure B1.

Left-hand panels: surface mass density maps of the NSD of the FCC 170 N-body model at different times, where all the other structural components are replaced by an analytical Hernquist (1990) mass model. Each column shows the edge-on (top panel) and face-on (bottom panel) views of the NSD particles at a given time. Each panel is 2 kpc wide. Right-hand panel: circular velocity of the bulge (blue dash–dotted line), NSD (orange dash–dotted line), and analytical Hernquist profile (green solid line) at different radii.

Open in new tabDownload slide