© The Authors 2024

1. Introduction

Galactic internal structures are important tracers of global galaxy formation and evolution. In particular, properties of massive stellar disks naturally allow us to reconstruct, with great detail, the mass assembly of disk-dominated galaxies. Only bright thin disks were initially identified in galaxies, with fainter thick disks detected much later. The latter were first discovered in edge-on galaxies, defined as a second (thicker) component that was necessary to fit luminosity vertical profiles of galactic stellar disks (Burstein 1979). A thick disk was initially defined, therefore, as a morphological component with a larger scale height than the thin disk. While the midplane region was dominated by the thin disk, the thick-disk dominated the light beyond a certain distance from the midplane of the galaxy. Only later on, when this second and thicker component was also identified in the solar neighborhood of the Milky Way, were stars in the thick disk found to have different properties from those in the thin disk. Thick-disk stars were associated with old ages, much lower metallicities, and higher velocity dispersions than those in the thin disk (Gilmore & Reid 1983). About one decade later, it was shown that the geometrically defined thick and thin disks, in the solar neighborhood, formed two different sequences in the [Mg/Fe]-metallicity plane, and that thick disks were enhanced in α elements (Fuhrmann 1998). Old ages, low metallicities, and high [Mg/Fe] abundances suggested that the Milky Way thick disk was formed at early times from gas that was still poorly enriched. Assuming that thick disks in external galaxies also showed similar properties, as was later shown (e.g., Yoachim & Dalcanton 2008b), understanding their origin means shedding light on the early stages of galaxy formation and evolution. On the other hand, young, metal-rich, and kinematically cool thin disks trace later and relatively recent phases of the evolution of a galaxy. Thus, tracing the full history of a disk galaxy requires understanding the processes leading to the formation of both thick and thin disks, as well as the evolutionary connection between them.

Different scenarios have been proposed to explain the early formation of the thick disk and the later transition to star formation in the thin disk. One possible explanation resides in turbulent gas at high redshifts, during the epoch of frequent gas-rich mergers (Brook et al. 2004). Thick disks were formed already thick from that dynamically hot gas. This gas settled down with time, allowing for the formation of stars in much thinner layers at later times and giving birth to the thin disks. A different scenario suggests that stars were instead formed with low velocity dispersions, in a relatively thin disk at high redshifts. These stars were later dynamically heated through processes such as mergers or encounters with clumps, either violently or slowly with time (Quinn et al. 1993; Villalobos & Helmi 2008; Bournaud et al. 2009; Di Matteo et al. 2011). Stars that are now older, such as those in the thick disk, had more time to be heated and were more exposed to violent processes (e.g., related to mergers and intense star formation) than young stars (in the thin disk). According to a third scenario, thick disks are mostly made up (more than 60% of their mass) by the debris of accreted satellites of ex situ origin (Abadi et al. 2003). Most of the thin disk would have formed during a later quiescent phase, during which no mergers occurred.

Since these scenarios, based on numerical studies, were initially proposed, observational works have been alternatively supporting one or another. Early thick-disk spectroscopic studies in edge-on galaxies used two long slits, placed on and off the galaxy midplane, respectively in the thin- and thick-disk-dominated regions. Yoachim & Dalcanton (2005, 2008b) extracted the kinematics of the thick and thin disks of a sample of nine edge-on late-type galaxies. They found a variety of kinematic properties, with high-mass galaxies (with a circular velocity above 120 km s⁻¹) showing similar rotation curves in the thick and thin disks, and low-mass galaxies showing strong differences. The substantial lag and counter-rotation identified in lower-mass galaxies strongly supported an accretion origin for thick disks. Yoachim & Dalcanton (2008a) observed (also in a sample of nine late-type galaxies) thick disks that were much older than the young thin disks, both with similar low metallicities. On the other hand, the long-slit study by Kasparova et al. (2016) showed three thick disks with diverse properties and proposed that different galaxies can form their thick disks via different scenarios.

The advent of integral-field spectroscopy (IFS) has brought important advances in this field, providing the two-dimensional kinematic and stellar-population structure of thin and thick disks and covering the transition region between them. IFS observations of thick disks are challenging due to their low surface brightness, especially in late-type galaxies, and require long integration times. Furthermore, galaxies need to be seen edge on in order to dissect their vertical structure into different components. However, there are not many edge-on galaxies with deep enough archival observations. For all these reasons, there are currently few IFS thick-disk studies in the literature, and most of them include only one or two galaxies. The first attempt to map thick-disk properties with IFS was done with the VIsible Multi-Object Spectrograph (VIMOS) at the Very Large Telescope (VLT) at poor spatial resolution (only one Voronoi bin in the thick disk; Comerón et al. 2015). The sharp age and metallicity differences between the thick and thin disk supported a fast thick-disk formation at high redshifts. Later studies with the Multi Unit Spectroscopic Explorer (MUSE), at a much higher resolution, revealed subtler differences between the thin and thick disks. Comerón et al. (2016) present an edge-on lenticular galaxy with clear metallicity differences between the thick and thin disks, although very similar old ages. This suggested a fast internal formation of both thin and thick disks at high redshifts.

Pinna et al. (2019a,b) added [Mg/Fe] maps, as well as star-formation histories coupled with chemical properties, to thick-disk studies. They analyzed MUSE observations of three lenticular galaxies in the Fornax cluster, located in regions of the cluster with different densities. In all cases, thick disks had similar properties, which suggested that they formed and grew via similar mechanisms. These two studies were the first to show a younger, more metal-poor component, interpreted as having been accreted, in addition to an old and metal-poor dominant population, which may have formed in situ. The authors proposed a complex thick-disk formation scenario with a dominant in situ formation combined with a significant contribution from accreted satellites. Stellar-population properties of thin disks varied from galaxy to galaxy (possibly related to different impacts of the cluster environment). Flaring in the outer disk was interpreted as a potential sign of dynamical heating through mergers. A similar complex thick-disk formation scenario was supported by the MUSE analysis from Martig et al. (2021) of one massive disk-dominated (Sb) galaxy. Here, an important merger left a significant contribution in the thick disk, and probably drove an enhancement of star formation in the thin disk. On the other hand, the lack of evidence for retrograde material in the MUSE kinematics of eight additional late-type galaxies (only five of them with a clear thick-disk component) suggested a predominant internal origin for thick-disk stars, still allowing for some fraction of accreted material (Comerón et al. 2019). Finally, a recent MUSE study of a young Sc galaxy showed a massive, thicker, metal-poor disk and an inner, small, very thin and metal-rich disk, which might be a future thin disk in its first stages of formation (Sattler et al. 2023). These results support a first buildup of a thicker disk (which occurred late in this case, since the galaxy is overall young), before the later inside-out formation of a thin disk.

Despite the recent progress on this matter, the origin of thick disks and their evolutionary connection with thin disks are still partially unclear. The abovementioned recent IFS observations revealed that mergers played a key role in the mass assembly of the thick and thin disks. However, the relative importance of ex situ and in situ components, and therefore the balance between direct accretion of ex situ stars and in situ star formation, which can also occur from accreted gas, is still uncharted waters. In this context, numerical simulations are paramount since they can offer large samples of galaxies, with a wide variety of different properties and an extended coverage of thin and thick disks, resolving some of the issues related to the observational work. Furthermore, the full evolution of thick and thin disks and their properties is tracked at different snapshots of the simulation. Several studies have used cosmological simulations to investigate the chemical bimodality in the [α/Fe]-metallicity plane of (thick- and thin-) disk stars. Grand et al. (2018a), based on the six highest-resolution simulated galaxies of the AURIGA suite (Grand et al. 2017), proposed an early and fast formation of the high-α sequence, during an intense star-formation phase induced by gas-rich mergers. The high-α component is usually associated with the thick disk in both the Milky Way (e.g., Hayden et al. 2015) and external galaxies (Scott et al. 2021). One of these gas-rich mergers would have provided the metal-poor gas necessary to start the formation of the low-α sequence, associated with the thin disk. This would have been followed by a more quiescent gas-enrichment phase driving the growth of the thin disk. Buck (2020), using four galaxies from the Ultra High-Definition (UHD) simulation suite (Buck et al. 2020) of the NIHAO (Numerical Investigation of a Hundred Astrophysical Objects) project (Wang et al. 2015), reached a similar conclusion. This idea was also supported by some observational studies (e.g., Scott et al. 2021). VINTERGATAN, a cosmological zoom simulation of a Milky Way-mass disk galaxy, supported a similar picture as well (Agertz et al. 2021).

Grand et al. (2020), also based on AURIGA simulations, proposed a “dual origin” for the Milky Way thick disk, which would comprise two components. One in situ component (with similar properties to the “starburst sequence” shown by An et al. 2023) would have formed during a starburst triggered by the gas-rich merger related to the Gaia-Enceladus Sausage (see also Ciucă et al. 2024). A second component would have formed from stars from a preexisting thinner disk that was dynamically heated by the same merger. This second component explains the positive velocity–metallicity relation that was found in the Milky Way thick disk (Belokurov et al. 2020). The analysis based on Milky Way-mass galaxies in the EAGLE suite of cosmological simulations shows that the chemical bimodality found in the Milky Way might not be very common in disk-dominated galaxies (Mackereth et al. 2018). This bimodality tends to appear only when accretion episodes of different intensities happen at different times. Early fast gas accretion and subsequent bursty star formation would lead to the formation of the high-α sequence. The geometrically defined thick disk would result from the combination between this inner high-α sequence and the flaring outer disk region of the low-α sequence (Mackereth et al. 2019). An alternative scenario, not involving galaxy mergers, hypothesizes that intense star formation in clumps naturally forms the high-α thick disk, while the low-α thin disk is the result of a more evenly distributed star formation (Clarke et al. 2019).

Additional studies of thick disks in galaxies, without a special focus on the Milky Way, have drawn similar conclusions as the abovementioned works. Martig et al. (2014a,b) analyzed the mono-age populations in seven galaxies from Martig et al. (2012) and find that galaxies with a quiescent merger history tend to show a continuous distribution of scale heights correlated with ages. However, mergers are capable of causing an age bimodality between a thick and a thin component. While most of the very old stars were born already dynamically hot, later generations were heated with time, partially due to mergers. García de la Cruz et al. (2021), expanding the sample to 27 simulated galaxies from Martig et al. (2012), fitted the vertical distribution of stars in each of them with two disk components. Part of the sample, with lower galaxy and thick-disk masses and a more quiescent merger history, shows a continuous and smooth age vertical gradient (mono-age populations whose scale height increases with age). The rest of the sample, with higher galaxy masses, showed a thicker morphology, higher thick-disk masses and mass fractions, and flaring thick disks (see also Minchev et al. 2015). Park et al. (2021) also morphologically decomposed the thick from the thin disks in 19 simulated galaxies and find that geometrically defined thick disks are predominantly formed in situ, although they host a large fraction of accreted stars. The ex situ contributions in AURIGA stellar disks were quantified by Gómez et al. (2017a), who showed that accreted stars are distributed in a thicker region, are older, and are more metal poor than the in situ disk component. Yu et al. (2021) used a kinematic definition for thick disks (orbital circularities lower than 0.8) and investigated their origin in FIRE-2 (Feedback In Realistic Environments-2) simulations. While the thick disk is mostly formed in an early bursty phase of star formation, later mergers also contributed to its growth with an additional component of kinematically heated younger stars.

In summary, thick-disk studies from cosmological simulations broadly support a predominant in situ formation during a bursty phase of star formation, with mergers playing a key role, potentially combined with dynamical heating and/or accretion of stars. However, numerical studies focused on a comparison with (IFS) observations of edge-on external galaxies, which would allow us to test the pictures drawn by observations and better interpret observed stellar-population properties, are currently lacking. In this work, based on a geometrical definition of the thick and thin disks, we use zoom-in magnetohydrodynamic cosmological simulations from the AURIGA project (Grand et al. 2017)¹. Projecting simulations to an edge-on view allows us to perform our analysis in a way similar to that of present and future observations of edge-on galaxies.

This paper is structured as follows. Section 2 describes the simulations, the data projection and spatial binning, and the selected sample. Section 3 explains the methods used for the analysis, including the morphological decomposition into a thick and a thin disk. We present our results in Sect. 4, discuss them in Sect. 5, and summarize our conclusions in Sect. 6.

2. Our sample of 24 simulated galaxies

We selected our galaxy sample from the original suite of AURIGA magneto-hydrodynamical zoom-in cosmological simulations of galaxy formation and evolution. We describe in this section the simulations and our selected sample.

2.1. AURIGA zoom-in cosmological simulations

The AURIGA project² (Grand et al. 2017) includes zoom-in magneto-hydrodynamical cosmological simulations of 30 Milky Way-mass late-type galaxies, with stellar masses between 10¹⁰ M_⊙ and slightly more than 10¹¹ M_⊙. These were obtained by re-simulating 30 halos from the parent largest-volume Dark Matter Only EAGLE simulation (Schaye et al. 2015), at intermediate resolution (L100N1504). These halos were selected to have a virial mass 1 < M₂₀₀/10¹² M_⊙ < 2 at z = 0, defined as the mass contained inside the virial radius defined as R₂₀₀, the radius enclosing a mean mass volume density of 200 times the critical density for closure. The 30 halos were also selected to be relatively isolated at redshift z = 0. In particular, they were randomly chosen from those in the lowest τ_iso, max quartile, where τ_iso, max is the maximum value of the tidal isolation parameter τ_iso, i = (M_200, i)/(M₂₀₀)×(R₂₀₀/R_i)³; here, M_200, i and R_i are respectively the virial mass of and distance to the ith halo in the simulation.

These simulations were obtained using an improved version of the moving-mesh code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020), which includes magneto-hydrodynamics and collisionless dynamics in a cosmological context. The initial conditions to be used in the re-simulations, and in particular the initial distribution of dark-matter particles, were taken from the halos in the parent simulation. Then, each dark-matter particle was substituted with a pair of a dark-matter particle and a gas cell. A fixed mass was assigned to each gas cell based on the value of the adopted cosmic baryon fraction (see Grand et al. 2017). AURIGA simulations offer different resolution levels, and in this work we use the fiducial resolution at “level 4”, the only suite including the full sample of 30 galaxies. This level corresponds to typical masses of high-resolution dark-matter and baryonic-mass particles of ∼3 × 10⁵ M_⊙ and ∼5 × 10⁴ M_⊙, respectively (Grand et al. 2017). The comoving gravitational softening length for stellar and high-resolution dark-matter particles was set to 500 h⁻¹ cpc, with a physical gravitational softening length of 369 pc at redshift z < 1. The softening length for gas cells was scaled by the cell size, from a minimum value of physical softening length of 369 pc to a maximum value of 1.85 kpc. This ensures that gas cells with lower densities (and thus larger volumes) have larger softening lengths than those with higher densities High-density gas cells can be smaller than their softening length since each cell has a given target mass, so high-density regions are resolved with more cells than low-density regions. Cell sizes range from 80 to 300 pc in the region within the galaxies that we analyze in this work (Sect. 3.1).

The galaxy formation model, including star formation, stellar and black-hole (BH) feedback, and magnetic fields, led to realistic galaxies matching a large variety of properties found in observed galaxies (e.g., morphologies and properties of structures such as spiral arms or bars, galaxy sizes and masses, and kinematic and chemical properties; see, e.g., Grand et al. 2016, 2017). The interstellar medium is characterized by two phases with cold and dense gas clouds embedded in a hot medium. These two phases are not modeled explicitly, but only as part of a sub-grid model (Springel & Hernquist 2003). Gas becomes thermally unstable and forms stars above a density threshold of 0.13 cm⁻³. According to a Chabrier (2003) initial mass function, each star-forming gas cell is either converted to a star particle or selected to be a site for type II supernova feedback. In the first case, the star particle represents a single stellar population (SSP) of specific masses, ages, and chemical abundances. In the second case, a single wind particle is launched in a random direction, and travels until it releases its energy and metals (40% of the total metals in the original gas cell) in the final gas cell. The rest of the metals remain local. Stellar feedback from type Ia supernova and asymptotic giant branch stars are also modeled, in terms of metals and mass loss. BHs with a seed mass of 10⁵ M_⊙h⁻¹ were introduced in the position of the densest gas cell in most-massive halos (> 5 × 10¹⁰ M_⊙ h⁻¹), and acquire mass from gas cells or merging with other BHs. BH feedback injects thermal energy to surrounding gas cells following two modes: in local isotropic quasar mode, or in a randomly isotropic, nonlocal radio mode. A homogeneous magnetic field was introduced at the beginning of the zoom-in re-simulations and let evolve. The simulations were saved in 128 different snapshots corresponding to different redshifts. For further details on the implemented physical models or general properties of these simulations, we refer the reader to Grand et al. (2017), who extensively described them.

2.2. Projection, spatial binning, and sample selection

We projected to an edge-on view the 30 AURIGA galaxies, to allow for an approach as similar as possible to what was done in previous observational works on edge-on disk galaxies (e.g., Pinna et al. 2019a,b; Martig et al. 2021). To perform the projection, we followed the approach explained in Walo-Martín et al. (2021). We first aligned galaxies with the axes of the reference system, by aligning with the vertical axis Z the average vector of the angular momentum of young stars (younger than 3 Gyr) located within a sphere of a radius of 60 kpc around the center of potential. Afterward, we projected the properties of all stellar particles to the YZ plane along the line of sight (LOS), which is parallel to the X axis and lying on the XY plane. The center of the coordinate system was placed on the center of the galaxy, calculated as the potential minimum of the main subhalo.

Since the number of particles decreases toward the outskirts of galaxies, we applied a Voronoi binning (Cappellari & Copin 2003) following Walo-Martín et al. (2021), to ensure a minimum number of particles and a roughly constant number of particles per bin. We used three different spatial-resolution setups, designed for different purposes. One used a projection into pixels of size 0.1 kpc  ×  0.1 kpc, and was designed to produce mock images and allow for smooth morphological fits of the thin and the thick disks (Sect. 3.2). We included all pixels with at least one star particle, and we performed a Voronoi binning to a target number of particles of 100 per bin, following requirements in previous work (Schulze et al. 2018; Walo-Martín et al. 2020). The second setup used a pixel size of 0.5 kpc  ×  0.5 kpc, including all pixels with more than two particles, and used a Voronoi binning to at least 900 particles per bin. This second setup was designed to allow the storage of the full distribution of the stellar-population parameters in each Voronoi bin, plus the information about the origin of star particles (ex situ or in situ, defined in Sect. 3.3), still in files of a reasonable size. With a pixel size well beyond the softening length (Sect. 2.1), and a number of particles per bin well beyond the usual minimum requirement of 100 particles, this setup is optimal for the analysis of the stellar properties of the thin and thick disks. It was used for stellar-population maps and the spatial distribution of accreted stars (Sects. 4.1 and 4.2). Finally, when mapping different simulation snapshots, we aimed at preserving spatial resolution at higher redshifts when galaxies are smaller and contain a lower number of particles. For this, we used a pixel size of 0.2 kpc  ×  0.2 kpc. We binned to 16 particles per bin and included all pixels with at least one star particle. However, the very low density of star particles at the highest redshift in our maps (z = 3.5) required a pixel size of 0.5 kpc  ×  0.5 kpc coupled with a Voronoi binning to only 9 particles per bin (Sect. 4.4).

The projection and binning process provided the average properties of particles in each Voronoi bin in the given projection, including surface brightness in different bands and stellar-population parameters. The photometric properties of the star particles were obtained by Grand et al. (2017) using stellar-population synthesis models from Bruzual & Charlot (2003), considering each star particle as a SSP. From the original sample of 30 Milky Way-size galaxies in the AURIGA suite of simulations, which were identified with the prefix “Au” plus numbers from 1 to 30 (Grand et al. 2017), we excluded from our analysis six of them (Au4, Au11, Au13, Au20, Au29, Au30). These galaxies showed, in their edge-on projection, strong distortions, lack of a disky shape or in one case an ongoing major merger. The main properties of the selected sample are indicated in Table 1. Most galaxies in the sample are barred, while only five of them are unbarred. Stellar masses vary between ∼3 × 10¹⁰ M_⊙ and ∼10¹¹ M_⊙.

3. Analysis methods

3.1. Analyzed region

Galaxies in our sample, when projected and binned as explained in Sect. 2.2, are much more extended than regions that are usually observed in edge-on galaxies, due to the low surface brightness in the outskirts of edge-on galaxies (examples in Comerón et al. 2016, 2018, 2019; Pinna et al. 2019a,b; Martig et al. 2021). We selected for the analysis in this work a region that was comparable to the region typically observed in the photometric study of edge-on galaxies (Comerón et al. 2018) and to the typical coverage reached in published studies using two deep MUSE pointings (at a distance between 20 and 30 Mpc, Pinna et al. 2019a,b; Martig et al. 2021; Sattler et al. 2023). We selected a region with a radius of 0.8 R_opt and a vertical extension of 2 h_scale (Table 1) above and below the midplane. R_opt is the optical radius, defined as the radius where the B-band surface brightness drops below 25 mag arcsec⁻², while the parameter h_scale was calculated as the standard deviation of the vertical positions of stellar particles at 0.5 R_opt, as defined in García de la Cruz et al. (2021). Thus, the analyzed region has a total size of 1.6 R_opt × 4 h_scale in the edge-on projection. This region was used for the morphological decomposition in Sect. 3.2 and all calculations throughout the paper, while for display purposes we use a central region of 50 kpc  ×  20 kpc, equal for all galaxies and still similar to the 1.6 R_opt × 4 h_scale region. We show in Fig. 1 mock images (surface brightness in the V band) of this 50 kpc  ×  20 kpc region for the 24 galaxies in our sample. These mock images were obtained by projecting to the edge-on view the luminosity density of the star particles in the V band (see Sect. 2.2), and converting the maps to the surface brightness. The setup with 0.1 kpc  ×  0.1 kpc pixel size and a Voronoi binning with a target number of 100 star particles per bin (see Sect. 2.2) was used for these images.

Fig. 1. Mock edge-on images, in terms of surface brightness in the V band (SB), of our selected sample of 24 Milky Way-like galaxies from AURIGA simulations. For display purposes, we show for all galaxies an area of the same size, of radius 25 kpc and height 10 kpc. Surface brightness was projected into pixels of 0.1 kpc  ×  0.1 kpc size and was Voronoi binned to obtain at least 100 star particles per bin. The region between the two horizontal lines is where the light of the thin-disk component dominates over that of the thick disk (see Sect. 3.2 for details). The region between the two vertical dashed magenta lines is dominated by a central component (the bar in barred galaxies and a classical bulge in non-barred galaxies). Names of the AURIGA halos are indicated on top of each panel.

3.2. Morphological decomposition into a thick and a thin disk

In this paper, we aim at analyzing thick and thin disks defined in a purely geometrical way, similarly to what is usually done in observations of edge-on galaxies (e.g., Comerón et al. 2015, 2016; Pinna et al. 2019a,b; Martig et al. 2021). This definition consists of defining a region dominated by the thin disk, close to the midplane, and a region dominated by the thick disk, at larger heights, based on a morphological decomposition into two disk components (see, e.g., Comerón et al. 2018).

We decomposed morphologically our edge-on projected galaxies into two disk components. For that purpose, we fitted vertical surface brightness profiles in a similar way to what is usually done in observations (e.g., Comerón et al. 2018). We included in the fit the region described in Sect. 3.1, of radial extension 0.8 R_opt following Comerón et al. (2018). We excluded the region within a radius of R_0, disk (Table 1), dominated by the bar in barred galaxies and by a classical bulge in non-barred galaxies. We divided the projected 1.6 R_opt × 4 h_scale area of each galaxy into four quadrants with a radial range between R_0, disk and 0.8 R_opt, two on each side with respect to the galactic center, one above and one below the midplane. We additionally divided each quadrant into two equal radial bins, obtaining eight bins in total. For each one of the eight radial bins, we extracted a median vertical profile of the luminosity density in the V band. We modeled the vertical luminosity density with a double hyperbolic secant square (following the study regarding one Galactica simulation by Park et al. 2021):

$$\begin{array}{c}\hfill \mathrm{\Sigma }(z)={\mathrm{\Sigma }}_{0,\mathrm{thin}}{\mathrm{sech}}^2\left(\frac{|z|}{2h_{z,\mathrm{thin}}}\right)+{\mathrm{\Sigma }}_{0,\mathrm{thick}}{\mathrm{sech}}^2\left(\frac{|z|}{2h_{z,\mathrm{thick}}}\right),\end{array}$$(1)

where Σ_0, thin and Σ_0, thick are the luminosity densities in the midplane of the thin and the thick disks, and h_z, thin and h_z, thick are their respective scale heights.

We used a Bayesian approach to fit these four parameters and estimate their uncertainties. We used the emcee Python implementation (Foreman-Mackey et al. 2013) of the affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler from Goodman & Weare (2010). We used 100 walkers and 20 000 chains. We maximized a likelihood p of the form ln(p) = − 0.5∑_n(SB_obs, n − SB_mod, n)², where SB_obs, n is the observed surface brightness on the point n, and SB_mod, n is the modeled surface brightness from Eq. (1) for the point n. This is a Gaussian likelihood assuming equal variance for all the points. We excluded from the fits the points fainter than 26 mag arcsec⁻² and the points obtained from averaging fewer than eight Voronoi bins. We used uniform priors to constrain the parameters to the values possible in the analyzed region, the thin disk to be brighter than the thick disk in the midplane, and to have a shorter scale height than the thick disk. We show in Fig. 2 one fit example of the vertical surface-brightness profile of one radial bin. For each profile, we calculated the height z_tT, i where the thick disk starts to dominate over the thin disk, as the height where the two hyperbolic secant squares cross each other (Fig. 2). We averaged these values for all eight vertical profiles, weighting by $\delta z_{tT,i}^{-2}$, where $\delta z_{tT,i}^2$ is the variance of the distribution of the z_tT, i value from the MCMC approach. With this approach, we obtained the average height z_tT at which the thick disk starts to dominate in each galaxy. z_tT values are indicated in Fig. 1 as orange, dashed and horizontal lines, and in Table A.1 together with the scale heights of the thin and the thick disks.

Fig. 2. Fit of the vertical surface-brightness (SB) profile of one radial bin of the galaxy Au5. Each small gray circle is the measurement for one Voronoi bin, while larger magenta circles correspond to the median profile of the gray points. The magenta shading indicates the 1σ dispersion of the gray points. Dashed blue and orange lines correspond respectively to the thin- and thick-disk fits, and the solid red line is a combination of the two.

For one galaxy, Au1, we did not obtain good fits (no good convergence of the free parameters for most radial bins), when modeling the vertical luminosity density profiles with Eq. (1). We obtained good fits when using only one hyperbolic secant square: Σ(z) = Σ₀sech²(|z|/2 h_z). This suggests that we do not have two clear morphologically distinct disk components within the analyzed region in this galaxy.

3.3. Maps of galaxy properties when seen edge on

We analyze in this work the stellar-population properties that are typically analyzed in recent observational work (Pinna et al. 2019a,b; Martig et al. 2021; Scott et al. 2021; Sattler et al. 2023). Maps of the mass-weighted projected star-particle age, [M/H] and [Mg/Fe] were obtained using the low-resolution setup (0.5 kpc  ×  0.5 kpc pixel size and Voronoi binning with a target number of 900 star particles per bin; see Sect. 2.2). In this work, we corrected the original values of [Mg/Fe] by applying a factor of 2.5 to magnesium yields, following van de Voort et al. (2020). They established this correction due to the underproduction of magnesium in the yields adopted in the simulations, if compared to observations. Since magnesium is not a cooling agent in the simulations, the need for this correction is not expected to affect the physics of the formation and evolution of the galaxy. Since the aim of this paper is to assess internal and external processes, we classified star particles in our simulations as in situ or ex situ, using the distinction made by Gómez et al. (2017a). In situ stars were formed in the main subhalo, either in the same region where they are observed now or in a different region, while ex situ stars were formed in a satellite galaxy before its final disruption. This distinction allows us to calculate the fraction of accreted stars in each Voronoi bin and to map their spatial distribution.

4. Results

4.1. Stellar-population maps

Maps of the kinematic parameters are included in Appendix B.1. We describe here the mass-weighted stellar-population maps, presented for our full sample of 24 galaxies in Appendix B.2, in Figs. B.5–B.7. We show in this section, in Figs. 3–5, only four representative galaxies as examples. For display purposes, we show for all galaxies a region extended 50 kpc  ×  20 kpc, centered on the center of the galaxy (Sect. 2.2). In general, clear differences in the stellar-population parameters are visible between the thin-disk-dominated region (between the two horizontal lines located at z_tT, excluding the region between the two vertical lines indicating R_0, disk) and the thick-disk region (above and below the region between the horizontal lines). Age maps in Figs. 3 and B.5 show in all galaxies a much younger thin disk than the thick disk. In some cases, a strong positive radial gradient is observed in the displayed region of the thin disk, with the youngest stars concentrated in an inner region (e.g., Au10, Au12, Au17, Au22). These are mostly the smallest galaxies, since we used a fixed-size box for all galaxies. In some other galaxies, with larger R_opt, young stars cover a region out to the outer regions of our maps (e.g., Au2, Au3, and Au24). In general, the flared appearance in the age maps reflects that younger stars flare, because they are formed in a flaring star-forming gas disk (Grand et al. 2016; see also Kawata et al. 2017 and Benítez-Llambay et al. 2018 for other simulations). This flaring extends to the outer regions of the geometrically defined thick disk, which presents a negative age radial gradient. While thin disks show ages as young as 1 − 2 Gyr in some Voronoi bins, thick-disk regions show Voronoi bins with ages as old as 11 Gyr in the inner regions in some galaxies (e.g., Au9; see Sect. 4.3 for average ages).

Fig. 3. Four representative examples of mass-weighted stellar age maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Ages were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

Fig. 4. Four representative examples of mass-weighted total stellar metallicity [M/H] maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal black dashed lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Metallicities were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

Fig. 5. Four representative examples of mass-weighted stellar [Mg/Fe]-abundance maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. [Mg/Fe] abundances were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

In Figs. 4 and B.6, thin-disk regions are clearly metal rich, with negative radial gradients and values as high as ∼0.4 dex in the inner region. Thick disks show more metal-poor values, mostly subsolar (often as low as −0.4 dex) or close to solar values (see also Sect. 4.3). Thick disks in AURIGA simulations are more enhanced in [Mg/Fe] than thin disks (Figs. 5 and B.7). While the latter show values closer to solar, thick disks show values as high as ∼0.2 dex or beyond in some cases (e.g., Au9). Our stellar-population maps (in particular metallicity maps) reveal signs of mergers in numerous galaxies: warps like in Au21 and Au27, streams as in Au14, and other distortions.

4.2. Spatial distribution of accreted stars

We show in Fig. 6 maps of the mass density of the accreted stars, in the central 50 kpc  ×  20 kpc region of four galaxies in our sample. Maps of the full sample are shown in Fig. B.8. These maps show that accreted stars are concentrated mostly in the inner region. In some galaxies, there are additional specific regions with a high density of accreted stars, tracing the remnants of recent minor mergers (e.g., Au17; see also below). Interestingly, while accreted stars are denser in the inner midplane region, their mass fraction (in each Voronoi bin) is usually higher in the thick-disk region, as shown in Figs. 7 and B.9. In some cases, such as Au12 or Au14, we reach up to 80% of accreted stars in some Voronoi bins in the thick-disk region. In general, ex situ stars constitute a high fraction of mass in the thick-disk-dominated region, while accreted fractions in the thin-disk region are low in comparison, usually below 10% in most Voronoi bins. In Au17 we can identify a region (below the thin-disk-dominated region) dominated by accreted stars (high density and high mass fraction of accreted stars; see Figs. B.8 and B.9). This is the remnant of an accreted satellite, still in the process of completely merging (brighter than its surroundings in Fig. 1). It can be identified in Figs. B.5 and B.6 as a younger and more metal-rich structure than its surroundings. It corresponds to slightly lower values of [Mg/Fe] in Fig. B.7, although differences with the surroundings are not obvious.

Fig. 6. Four representative examples of maps of the mass surface density of accreted stars, in logarithmic scale. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps. Maps of the full sample are shown in Appendix B.3.

We provide in Table 2 the total mass and mass fraction of accreted stars, computed in the full analyzed region of the galaxy (of size 1.6 R_opt × 4 h_scale, Sect. 3.1). These mass fractions were also indicated on top of each panel of Figs. 7 and B.9, together with the name of the AURIGA halo. Fractions of accreted stars are below 10% for half of the sample, and relatively low in almost all galaxies (below 20%). Exceptions are Au7 with almost 30% of accreted stars and Au1 with ∼46%. However, these fractions are much larger in thick-disk regions, ∼22% on average, with ∼50% of accreted stars in Au7 and ∼55% in Au12. These fractions correspond to total accreted masses of about 7 × 10⁸ to 2 × 10¹⁰ M_⊙ in the galaxy (4 × 10⁸ to 10¹⁰ M_⊙ in the thick disks). Mass densities and fractions of accreted stars are larger if a more extended region (than the one used in this work; see Sect. 3.1) is taken into account, since the ex situ fraction per Voronoi bin increases the farther we go from the galactic midplane. Gómez et al. (2017a) showed that, in AURIGA disks defined based on circularity cuts (circularities larger than values from 0.7 to 0.9), the fraction of accreted stars is lower in dynamically colder disks. They also showed that a positive radial gradient of the accreted fraction is present in most galaxies. Although we do not apply any filter based on kinematic properties to the star particles, we find similar results (see Sect. 5.1 for an extended discussion).

Fig. 7. Four representative examples of maps of the mass fraction of accreted stars. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and the total mass fraction of accreted stars (calculated in the analyzed region, of size 1.6 Ropt × 4 hscale, Sect. 3.1) are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps. Maps of the full sample are shown in Appendix B.3.

4.3. Comparison of thick and thin disks

We show in Figs. 8 and 9 a comparison of mass-weighted average thick-disk age, [M/H] and [Mg/Fe] with the respective properties in the thin disks. Each point corresponds to one of the 24 galaxies in our sample. Thick disks are older, more metal poor and [Mg/Fe] enhanced than thin disks in all 24 galaxies in our sample, as suggested by maps in Sect. 4.1. They are from ∼1.5 to ∼4.5 Gyr older (on average ∼3 Gyr older), from ∼0.13 to ∼0.35 dex (on average ∼0.25 dex) more metal poor, and from ∼0.04 to ∼0.09 dex (on average ∼0.06) more [Mg/Fe]-enhanced than thin disks. However, the ages of thick disks are distributed over a wide range and some thick disks are younger than other thin disks in the sample (left panels in Figs. 8 and 9). The average [M/H] of thick disks range from subsolar to slightly super-solar values, while thin-disk average [M/H] is always super-solar, higher than the most metal-rich thick disk in the sample.

Fig. 8. Thick-disk versus thin-disk stellar-population properties. From left to right: age, [M/H], and [Mg/Fe] abundance. Points are color coded according to the total mass of the galaxy (Table 1). The one-to-one line is indicated in black and dashed, in each panel. A red “+” symbol indicates solar metallicity in the middle panel.

Fig. 9. Thick-disk versus thin-disk stellar-population properties. From left to right: age, [M/H], and [Mg/Fe] abundance. Points are color coded according to the total mass fraction of accreted stars in the analyzed region of the galaxy (Sect. 3.1). The one-to-one line is indicated in black and dashed in each panel. In the middle panel, the red “+” symbol indicates solar metallicity.

In Fig. 8, we have color coded each point according to the total stellar mass of the galaxy (Table 1). There is no clear trend of the thick-disk properties with galaxy stellar mass (with a R² correlation coefficient lower than 0.11). However, galaxies with the lowest masses show more similar properties in the thin and the thick disks (points closer to the one-to-one line). In Fig. 9, we have color coded points by the total mass fraction of accreted stars (f_accr, gal in Table 2). While the total accretion fraction does not show a clear trend with thick-disk [M/H] and [Mg/Fe], it shows an anticorrelation with thick-disk age (R² ∼ 0.5), with the oldest thick disks located in galaxies with the lowest accreted fractions. We find a similar age trend, but with a shallower slope, also with the accretion fractions in the thick disk. This means that the youngest thick disks tend to be hosted by galaxies with the largest fractions of accreted stars. This is usually the outcome of more recent mergers, involving satellite galaxies that have had more time to grow in mass (e.g., in Au7, see also Fig. 10), and heating stars originally formed in preexisting thin disks (e.g., in Au19, see Fig. B.23). Star-formation histories of these thick disks (which will be published in a companion paper) support this scenario, and also show that mergers provide additional dynamically hot gas, which leads to some later star formation at heights corresponding to the geometric thick disks (see also the discussion in Sect. 5.1). On the other hand, thin-disk [Mg/Fe] abundances (mildly) correlate with the accretion fraction in the right panel of Fig. 9 (R² ∼ 0.3). Galaxies with less significant accretion reach lower [Mg/Fe] values in their thin disks. We also checked for trends of thin- and thick-disk stellar-population properties with the presence of bars or with galaxy morphology in general, but no trends were found.

Fig. 10. Time evolution of the galaxy Au7 seen edge on. Surface brightness is shown in the left column, total metallicity [M/H] in the middle column, and [Mg/Fe] abundance in the right column. From top to bottom, the stellar component is shown for seven different snapshots. On top of the middle panels, the redshift and the look-back time (LBT) corresponding to the snapshot shown in that row are indicated. For comparison, horizontal dashed lines separate the regions where the thin or the thick disks dominate at z = 0, while the vertical lines indicate the region where the central component (bar or classical bulge) dominates at z = 0.

4.4. Time evolution

We analyze in this section the evolution of thick and thin disks in different snapshots of the simulations across time. We mapped all galaxies in our sample at seven different snapshots from redshift z = 3.5 to z = 0. Galaxies were aligned and projected edge on in each snapshot, according to the criteria explained in Sect. 2.2. Since the alignment of the galaxies is based on the angular momentum of young stars, which at high redshift is not as clear as at z = 0, it did not give an optimal result when the galaxy was not yet clearly rotating or in the presence of a merging satellite with (differently) rotating young stars (e.g., z = 3.5 in Fig. 10; see also Appendix B.4, e.g., Au26 at z = 1.7). For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size in each simulation snapshot. Since at high redshift the number of stellar particles is much lower, it was needed to adapt the pixel size (to 0.2 and 0.5 kpc) and the Voronoi binning for different redshifts (to 9 and 16 particles per bin, see Sect. 2.2).

We show here two specific examples, Au7 and Au18, respectively in Figs. 10 and 11. These two galaxies are representative of the sample as they show two typical different types of evolution histories (although each galaxy has its own peculiarities): one with a significant ex situ contribution and the other with a predominantly in situ growth. Snapshots for the rest of the galaxies are shown in Appendix B.4. While we show here only surface brightness, [M/H] and [Mg/Fe], velocity maps (not shown) helped in the identification of (rotating-)disk formations. Au7, with a stellar mass of 4.9 × 10¹⁰ M_⊙ at z = 0, had a slower evolution than Au18, which is slightly more massive (with a stellar mass of 8.0 × 10¹⁰ M_⊙).

Fig. 11. Same as Fig. 10 but for the galaxy Au18.

Au7 had formed only a few, metal-poor and [Mg/Fe]-enhanced stars at z = 3.5, and had no disky shape. At z = 1.7, it showed a small and relatively thick disk. It had a similar vertical extension to the future thin disk at z = 0 (indicated by the two black horizontal lines), and thus a significantly lower vertical extension than the thick disk at z = 0. However, this primordial disk was radially much less extended than both the thin and thick disks at z = 0 (see the vertical dashed black lines as reference), leading to a similar thick shape to its shape at z = 0. This primordial thick disk covered a region of approximate size 0.4 R_opt × h_scale, with the same shape (intended as the ratio between the vertical and horizontal size in projection) as our analyzed region at z = 0 (1.6 R_opt × 4 h_scale). It had a stellar mass of about 5.4 × 10⁹ M_⊙ at z = 1.7, about 10% of the final galaxy mass at z = 0. No solar-metallicity stars were yet present in this second snapshot. Later on, the disky shape and rotation were distorted by mergers (z = 1 and 0.79). Around z = 0.6 the thin disk was developing from inside out embedded in the thick disk. A final important merger, tilting the disk as visible in the z = 0.36 snapshot, provided a large amount of gas to allow for the full growth of the thin disk to the outskirts at z = 0.

Au18 formed its disks and grew much faster than Au7. At z = 3.5, it was already more extended than Au7. Its fast formation timescale led to a global higher [Mg/Fe] abundance, as well as a higher [M/H] in the inner region (solar values were already reached). The galaxy already started to show some amount of rotation. At z = 1.7, its stellar mass was already about ∼4.7 × 10¹⁰ M_⊙ and most of its total mass at z = 0 had already formed. It was approximately one order of magnitude more massive than Au7 at the same time. A primordial thinner metal-rich and a thicker metal-poor disks were already in place. From that time on, the thicker component evolved only slowly during the rest of the life of the galaxy. While in Au7 the low density at early stages and the disturbed shapes due to the merger later on make it difficult to fit surface brightness profiles for different redshifts, we have done so for Au18. These fits confirm that thick-disk scale height varied very little (within their uncertainties) across the galaxy life. Thick-disk [Mg/Fe] and [M/H] remain almost unchanged after z ∼ 1. On the other hand, the buildup of the metal-rich thin disk had just started, it grew inside out to its full final extension at z = 0, became very metal rich and its [Mg/Fe] abundance decreased toward solar values.

At z = 0, Au7 is thicker, as indicated by its h_scale, than Au18 (Table 1). The thick disk is younger in Au7 than in Au18, as displayed in Fig. 3. Half of the mass in the thick-disk region of Au7 was accreted (Table 2, as also suggested by the merger in Fig. 10 at z = 1 and 0.36). On the other hand, only ∼7% of the thick-disk mass in Au18 was accreted, and it was already mostly (in situ) formed at z = 1.7. This is connected to the fact that, in general, the oldest more [Mg/Fe]-enhanced thick disks are found in more massive galaxies with a less intense merger history (at least after z ∼ 1, see Figs. B.10–B.31).

5. Discussion

5.1. Origin of thick disks and their evolutionary connection with thin disks

We have shown in previous sections that geometrically defined thick disks are relatively old, metal poor and enhanced in [Mg/Fe] in AURIGA simulations. This indicates that most of their masses were formed a long time ago from gas that had not been significantly enriched. In fact, Sect. 4.4 showed that a primordial thick rotating disk (which later kept growing) formed in most galaxies more than 10 Gyr ago (around redshift z = 2). This happened in a comparable timescale to what was found for the Milky Way (e.g., Gallart et al. 2019; Belokurov & Kravtsov 2022; Conroy et al. 2022). Only in Au15 and Au19 we observe distorted shapes (and no clear rotation) in the snapshots at redshift z = 1.7 (Figs. B.20 and B.23). However, the primordial thicker components had very different sizes, compared to the final thick disks at z = 0, in different galaxies (see the horizontal dashed lines in the figures). Some had already extended thick disks (vertically and radially, e.g., Au2, Au17, and Au18), while others were smaller and either thinner (Au7) or thicker (e.g., Au3, Au10, and Au12) than the thick disk at z = 0.

In Sect. 4.2, we have seen that all thick disks in the sample are made up of in situ and ex situ components, whose mass ratios can vary significantly for different galaxies. Most of the mass in the thick-disk regions comes from in situ stars, with a mass fraction of accreted stars lower than 50% in all cases except in Au12 and Au7. However, a significant accreted component indicates that geometric thick disks grew thanks to the contribution from merging satellites across time. While the balance between in situ and ex situ formation strongly depends on the individual galaxy, mergers played a key role in the mass assembly of thick disks. We would like to add a word of caution about the implications on the discussion about thick-disk formation of a pure geometric definition of them (chosen because our aim is to provide clues for the interpretation of observational results of edge-on galaxies). Geometric thick disks are defined by the region where their light dominates. However, in that region, there might be a contribution from stars belonging to other structural or kinematic components overlapped in the LOS and that are not possible to disentangle morphologically. These components include stellar halos, often defined in the Milky Way or in simulations using kinematic properties (e.g., Naidu et al. 2020; Gallart et al. 2019; Monachesi et al. 2019), as well as flaring thin disks (Sect. 4.1). These are very challenging to be disentangled in observed edge-on galaxies.

Here, we take advantage of the combination of spatial three-dimensional information with kinematic properties in simulations, to quantify the contamination from stars belonging to other (non disk-like) kinematic components in the geometric thick disks. In general, stellar-halo stars are thought to have been mostly accreted from satellites (about 80% of halo stars were accreted in the Milky Way, Naidu et al. 2020, and about 90% in AURIGA galaxies, Monachesi et al. 2019). We take the opportunity to discuss a potential contamination of halo stars in the ex situ component of thick disks. Within the analyzed thick-disk region, we disentangled the fractions of stars with thick-disk kinematic properties, with orbital circularities ϵ > 0.5 (following Yi et al. 2024), or associated with spheroidal components such as the stellar halo (ϵ < 0.5). Stars with ϵ < 0.5 make up on average 38% of the mass in our geometric thick disks. However, it can be as low as 13%, it is less than 35% in half of the thick disks, while in about one third of the galaxies this fraction is more than 50% (up to 63% in one case). Thus, the degree of contamination can be significant depending on the specific case. Regarding the accreted stars in the thick-disk region, the mass fraction of the ones with ϵ < 0.5 can be as low as 14% and as high as 79%, on average 57% in our sample. Hence, the contamination from ex situ halo stars in the accreted component in thick-disk region can be insignificant in some cases (since it is much lower than typical differences in the accreted fraction from galaxy to galaxy), or dominant. This contamination, while it is challenging to be corrected, should be discussed in observations of edge-on galaxies, where thick disks are also defined geometrically. A potential way forward would come from combining a geometrical definition with de-projected kinematic properties extracted with orbit-superposition methods (Zhu et al. 2018). Despite this word of caution, our results show that the impact of mergers through direct accretion of stars to the thick-disk-dominated region is important.

In addition to star accretion, we expect mergers also to impact on the in situ component of thick disks, through the dynamical heating of stars closer to the midplane and in situ star formation after gas accretion. The vertical dynamical heating of AURIGA stellar disks by mergers was previously discussed by Grand et al. (2016, 2020) and Monachesi et al. (2019). Gas accretion has also a strong impact in the baryonic cycle and chemical evolution of galaxies, as shown in AURIGA simulations by Grand et al. (2019). About half of this gas comes from accreted satellite galaxies. In a companion paper, we will show the star-formation histories of thick disks and will discuss further the in situ star formation and ex situ contribution of thick disks. Thus, the processes responsible for the formation of thick disks are both internal and external, with a combination of the three main scenarios presented in Sect. 1. After the formation of a primordial thick disk at high redshift, this kept growing, through accretion, in situ star formation, and dynamical heating.

Regarding the thin disks, in some galaxies they initially formed at high redshift. Galaxies with a faster evolution at their early stages, implying that they also formed extended thick disks very early (e.g., Au2, Au17, and Au18), started to form a thin metal-rich disk in the inner region about 10 Gyr ago (z = 1.7). Our results suggest that while thick disks systematically start forming at early times, the formation of thin disks can either follow the primordial thick-disk formation or can be delayed until later times.

5.2. Comparisons with observations

5.2.1. Morphology of thin and thick disks: Stellar disks are thicker in simulations

We show in Fig. 12 the distributions, in terms of frequency histograms, of the thick and thin-disk scale heights of the 24 galaxies in our AURIGA sample, obtained as described in Sect. 3.2 and gathered in Table A.1. Thick disks mostly have scale heights between 1.6 and 3.5 kpc (with three outliers, see also below), with an average value (weighted by the uncertainties) of 2.6 kpc. Thin disks have scale heights between 370 and 810 pc, on average 620 pc (see also Table A.1). We compare here the thickness of our thick and thin disks with the 124 galaxies morphologically fitted by Comerón et al. (2018) that have two clear distinct regions dominated by the thick or the thin disk. We find on average thicker thick and thin disks in AURIGA simulations than in this observed sample. Thick disks in the sample from Comerón et al. (2018) have a smaller average scale height (∼1 kpc, magenta dashed vertical line in Fig. 12), but a maximum scale height of ∼3.7 kpc (magenta dotted vertical line). All except three thick disks from AURIGA have scale heights below this upper limit. We have not found any correlation between h_z, thick and the fraction of stars with a halo-like kinematics in the thick-disk-dominated region (Sect. 5.1), reason why we exclude a potential contamination from halo stars as responsible for the higher scale heights of these three galaxies. Observed thin disks in Comerón et al. (2018) have very similar upper limit (green dotted vertical line in Fig. 12) as our thin disks, but an average scale height of ∼200 pc, much lower than ours (respectively vertical dashed green and blue lines).

Fig. 12. Frequency histograms of hz, thin and hz, thick in our AURIGA sample, with 0.1 and 0.5 kpc bins, respectively (solid blue and red lines). These distributions have a weighted average value of 0.6 and 2.6 kpc, respectively, indicated as vertical dashed blue and red lines. Vertical green and magenta lines indicate the average (dashed) and maximum (dotted) thin- and thick-disk values from observations in the Comerón et al. (2018) sample.

While the fact that stellar disks in AURIGA are thicker than observed disks was already pointed out (Grand et al. 2017; Gómez et al. 2017b), this may be partially due to some limitations of current cosmological simulations (e.g., Martig et al. 2014b). One potential issue is numerical heating. Due to the resolution limitations of the simulations (in particular the mass of dark-matter particles, often more massive than baryonic particles), and to equipartition of energy, an additional two-body scattering increases random motions of the star particles (Ludlow et al. 2021). As said by Ludlow et al. (2021), AURIGA simulations are unlikely to be strongly impacted by this effect, as this is a more common problem of cosmological volume simulations, with more massive dark-matter particles. Grand et al. (2017) showed in their Figs. 24 and 25 that changing resolution does not have a strong impact on the thickness of disks. However, the comparison between scale heights at solar radius in the AURIGA simulations of different resolution levels (for simulations used here in Grand et al. 2017 and for higher resolution in Grand et al. 2018b) shows that resolution has some impact on disk scale heights.

Another potential culprit is the interstellar medium model. Gas disks are already very often too thick in simulated galaxies (e.g., Trayford et al. 2017). This is related to the fact that, to limit computational costs (and to avoid the spurious fragmentation of gas structures), the interstellar medium does not explore the full temperature – pressure parameter space in most current cosmological simulations. An artificial “pressure floor” leads to a longer Jeans length at the star-formation threshold, so that stellar disks cannot be much thinner than that. This results in thicker disks than in observed galaxies (Bahé et al. 2016; Trayford et al. 2017). This has been widely discussed in the literature (e.g., Gabor et al. 2016; Crain et al. 2017; Sullivan et al. 2018; Weinberger & Hernquist 2023) and it is illustrated for AURIGA simulations in Kelly et al. (2022, in their Fig. 11, compared to APOSTLE simulations). Nevertheless, it is worth noting that Marinacci et al. (2017) showed that the thickness of HI disks in AURIGA might be in agreement with HI observations. Finally, stellar feedback likely plays an important role in setting the scale heights of cool gas (e.g., Roškar et al. 2014; Marinacci & Vogelsberger 2016) and remains an active area of study.

5.2.2. Spectroscopic observations of edge-on galaxies

Throughout the paper we have used an approach that is oriented to facilitate the comparison of these results with observations of edge-on galaxies. As mentioned in Sect. 1, detailed observational studies of thick disks through IFS are limited to a few galaxies, in particular for spiral galaxies (Yoachim & Dalcanton 2008a; Comerón et al. 2015; Martig et al. 2021). We should be cautious with comparing our results, from AURIGA Milky Way-mass spiral galaxies, with observations of galaxies of different morphology, mass, and environment. Furthermore, observational stellar-population parameters of thick disks in different studies were sometimes derived with different methods and we should also be cautious in comparing them to each other. Despite all this, most observations displayed older, more metal-poor and [Mg/Fe]-enhanced thick disks than their corresponding thin disks, independently of the morphological type and environment of their host galaxy (Yoachim & Dalcanton 2008a; Comerón et al. 2015; Pinna et al. 2019a; Martig et al. 2021). For these reasons, the following discussion is mostly based on a qualitative comparison.

Thick disks in our study, with ages between 5 and 9 Gyr, are on average not as old as most observed thick disks. In particular, we compare qualitatively our sample to the only two edge-on spiral galaxies with clearly distinct thick and thin disks, with published two-dimensional stellar-population analyses from IFS data. One is NGC 5746, an Sb galaxy of similar mass to the most massive AURIGA galaxies, with a ∼12 Gyr-old thick disk (Martig et al. 2021). Another is the Sc galaxy ESO 533-4, of similar mass to the least massive AURIGA galaxies, and with a ∼10 Gyr-old thick disk (Comerón et al. 2015). Both thick disks are older than all AURIGA thick disks, and the difference is larger than the estimated uncertainties in observations (∼1 Gyr in Martig et al. 2021). Some lenticular galaxies show old thick disks of similar ages as NGC 5746 and ESO 533-4 (Comerón et al. 2016; Pinna et al. 2019a,b), while others show thick disks of similar ages of our youngest thick disks in AURIGA (Kasparova et al. 2016, 2020).

The fact that young stars strongly flare in AURIGA disks (Fig. B.5), contaminating the thick-disk heights with younger stars, could be contributing to the fact that the average age of thick disks is never as old as the inner thick disk (see, in Fig. B.5, the age difference between the inner and the outer Voronoi bins in the thick disk of, e.g., Au9). In observations, such clear and strong flares have not been detected in age maps, because the outer region of the disks is not covered by the data (e.g., Martig et al. 2021), because of the lower spatial resolution in the faint outskirts, or just because they are not as strong in most real galaxies. In some cases they were instead observed in metallicity maps, suggesting that age differences may be present but sometimes difficult to be recovered (e.g., more than metallicity differences; see the discussion in Pinna et al. 2019a,b). Finally, another reason why thick disks are younger in AURIGA could be that star formation occurs in thicker gas disks even at low redshift, as discussed in Sect. 5.2.1. This could also explain the fact that thick-disk metallicities are also on average higher than observed, reaching sometimes solar or slightly super-solar values. However, the two thick disks observed in NGC 5746 and ESO 533-4 show [M/H] values in the range covered by our simulation sample. While Comerón et al. (2015) did not provide [Mg/Fe] values for ESO 533-4, NGC 5746 shows higher values than our thick-disk sample.

Regarding thin disks, they have shown a larger variety of different stellar-population properties than thick disks in observations. Furthermore, Pinna et al. (2019a,b) showed that their formation and evolution, happening in general later than thick disks, is more affected by the environment of their host galaxy than thick disks. Some observed thin disks have similar properties to the ones in our simulations. However, very old thin disks have been also observed in lenticular galaxies in galaxy clusters (Comerón et al. 2016; Pinna et al. 2019b). These are absent in our sample, made up of relatively isolated spiral galaxies. On the other hand, the few measured thin-disk [Mg/Fe] values in observations (Pinna et al. 2019a,b; Martig et al. 2021) are in good agreement with the values in the right panels of Figs. 8 and 9. We have already mentioned several potential causes of the identified differences between stellar-population properties in observed and simulated thin and thick disks. In addition to that, observational effects that can impact the recovery of stellar-population parameters due to, for example, projection, systematics related to the SSP models and their age, [M/H] or [Mg/Fe] resolution, and instrument spatial and spectral resolution, as well as the presence of dust (although not expected in large amounts in thick disks), were not corrected in observations and not taken into account in simulations. Furthermore, recent studies have suggested that full-spectrum fitting techniques tend to bias the results toward older ages (e.g., Pinna et al. 2019b). While producing mock spectra is beyond the scope of this paper, it might be the future way to go to produce results that are more quantitatively comparable to observations.

While accreted stars are easy to be mapped in simulations where we can track the particles and classify them into in situ and ex situ, this is challenging to be done in observations. Pinna et al. (2019a,b) and Martig et al. (2021) identified accreted stellar populations from their distribution in the age - metallicity plane, and mapped their spatial distribution. While this method is completely different from what is done here, they also found that accretion has a much more significant impact (measured in terms of mass fraction) in the thick-disk region, similarly to our results (Fig. B.9). On the other hand, the mass density of accreted stars is higher in a central region of a few kiloparsecs both in all our simulated galaxies (Fig. B.8) and in the observed ones. Since no kinematic cut was applied in the observations or simulations, stars in the thick-disk-dominated region may in both cases contain some contamination from stars with halo kinematics (see Sect. 5.1).

6. Conclusions

The origin of thick disks and the connection with thin disks is still a matter of debate despite recent progress. Detailed observational studies (in particular from IFS) are so far limited to individual galaxies or very small samples, and simulation studies do not usually provide results that are directly comparable to observations. We analyzed the properties of thick and thin disks, defined purely geometrically, in 24 galaxies from the AURIGA suite of zoom-in cosmological simulations. We projected the star particles of the 24 galaxies in an edge-on view, applying Voronoi binning. We fitted vertical surface brightness profiles in the V band with two disk components, enabling a geometrical definition of the thick and thin disks. We analyzed the age, total metallicity [M/H], and [Mg/Fe] abundance of the galaxies seen edge-on. After classifying the star particles into in situ and ex situ components, we also mapped the distribution of accreted stars. From these results, we drew the following conclusions:

• Thick disks are older, more metal poor, and more [Mg/Fe] enhanced than thin disks, indicating that most of their mass was formed at earlier times from gas that was not significantly enriched. They are on average ∼3 Gyr older, ∼0.25 dex more metal poor, and ∼0.06 more [Mg/Fe]-enhanced than thin disks.

• In most galaxies, primordial thick disks were already in place 10 Gyr ago (z = 1.7), with a shape and size sometimes very different from those at z = 0, depending on the specific galaxy.

• Most thick disks had most of their mass formed in situ, with a significant fraction (on average ∼22%) of accreted stars contributing to their growth. Two galaxies had half of their thick disks made up of ex situ stars. In only four galaxies did satellite debris contribute less than 10%. These results suggest that, in general, accretion plays a very important role in the growth of geometric thick disks, although internal formation is the dominant mechanism.

• Accreted stars are mostly distributed in the inner region of galaxies (where they have higher mass densities). However, their mass fraction is much more significant in the thick-disk-dominated regions. A very similar trend was found in observations.

• Geometric thick disks can be contaminated, significantly in some cases, by stars with halo kinematic properties that cannot be disentangled in observations of edge-on external galaxies. The halo contamination is in general greater in the accreted component of geometric thick disks.

• Thick disks in galaxies with higher accreted fractions tend to be younger. Mergers contributed younger stars than the in situ component, especially relatively recent ones.

• Mergers also cause additional thickening in preexisting disks, scattering cooler stars to greater distances from the midplane and thus contributing to the growth of the thick disk.

• Thin disks, which are young, metal-rich, and poorly enhanced in [Mg/Fe] in all simulations, start to form at high redshifts only in some galaxies. In others, with recent more active merger histories, they formed only at late times.

To sum up, thick disks result from a combination of external processes with the internal evolution of the galaxy. The balance between these processes, as well as the time of formation of thin metal-rich disks, depends on the specific galaxy and its merger history. In one companion paper, we will show the star-formation histories of thick disks in this sample and further discuss the in situ star formation and ex situ contribution in thick disks. In another companion paper, we will investigate bimodalities in the [Mg/Fe] – [M/H] plane, their correspondence with geometrically defined thick and thin disks, and how bimodality can be recovered in observations. While this work will be essential for the interpretation of future observations of Milky Way-mass galaxies, shedding more light on the properties and origin of thick disks and their connection with thin disks, additional simulations with a larger variety of galaxy masses, morphological types, and environments will be needed to achieve a complete view of how thick and thin disks form. Finally, comparing stellar properties from the direct star-particle analysis with results from mock spectra obtained from the simulations, and comparing the mock spectra with observations, will be the way forward toward understanding IFS observations.

Acknowledgments

R.G. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through the Spanish State Research Agency, under the Severo Ochoa Program 2020–2023 (CEX2019-000920-S), and support from an STFC Ernest Rutherford Fellowship (ST/W003643/1). F.F. acknowledges support from STFC grants ST/T000244/1 and ST/X001075/1. F.A.G. acknowledges support from ANID FONDECYT Regular 1211370, the Max Planck Society through a “Partner Group” grant and ANID Basal Project FB210003. We would like to thank Sylvia Plöckinger for the useful discussions.

References

Appendix A: Results from the morphological decomposition into two disk components

We provide in Table A.1 the average heights (z_tT) beyond which the light of the thick disk dominates over the thin-disk light, as well as the scale heights of the thin and the thick disks (respectively h_z, thin and h_z, thick), resulting from the fits described in Sect. 3.2. Uncertainties (δz_tT, δh_z, thin and δh_z, thick), also provided in Table A.1, were calculated propagating the 1σ uncertainties obtained from the Bayesian approach for single radial bins. For Au1, fitted with only one disk component (Sect. 3.2) as also indicated in the second column of Table A.1, we provide the scale height and its uncertainty of this full-disk component.

Appendix B: Additional maps

B.1. Stellar kinematic maps in an edge-on projection

Kinematic parameters shown here were calculated by fitting the line-of-sight velocity distribution (LOSVD) of star particles with a Gauss-Hermite expansion, which is widely used for observations (van der Marel & Franx 1993; Gerhard 1993; Cappellari & Emsellem 2004). We show here the maps of the first four moments of the LOSVD, often shown in observations of edge-on galaxies (Guérou et al. 2016; Pinna et al. 2019b,a; Martig et al. 2021; Sattler et al. 2023). As in the rest of the maps shown in this paper, we map a region extended 25 kpc radially and 10 kpc in height, with a total size of 50 kpc  ×  20 kpc. Velocity (V) maps, in Fig. B.1, show strong rotation in the thin-disk region and slower rotation in the thick disks, with larger (e.g., in Au9) or smaller (Au3) thick-disk rotation lags, depending on the galaxy. Velocity dispersion (σ) maps show in Fig. B.2 that stars in the thin-disk midplane region are in general kinematically cooler than stars in the region dominated by thick disks. However, while differences between the two components are clear in most galaxies (see, e.g., Au5 and Au28), they are not clear in some others with similar thick- and thin-disk properties (such as Au8 and Au25; see also Sect. 5.1). In some galaxies, a much higher velocity dispersion is seen in the central region (see, e.g., Au10 and Au28), mostly related to the presence of a (boxy or classical) bulge and showing a clear “X” shape in Au14.

Fig. B.1. Mean-velocity (V) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

Fig. B.2. Velocity-dispersion (σ) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

Skewness (h₃) is shown in Fig. B.3 and is characterized by a clear anticorrelation with V in regions dominated by disk-like rotation (as observed by, e.g., Krajnović et al. 2008; Pinna et al. 2019b,a). The highest absolute values of h₃ can be found either in the thin- or the thick-disk regions, depending on the galaxy. In the central region of barred galaxies (see in particular Au2, Au6, Au9, Au14, Au16, Au18, Au23, Au24, Au27, Au28), h₃ is correlated with V, as it is expected in bar-dominated regions characterized by cylindrical rotation (see also, e.g., Pinna et al. 2019b). Kurtosis (h₄), in Fig. B.4, shows different patterns in different galaxies. Transitions from negative values (indicating a LOSVD with a broader symmetric profile than a pure Gaussian) to positive values (indicating a LOSVD with a narrower symmetric profile than a pure Gaussian), or vice versa, are observed from the inner to the outer regions of thin and thick disks (see, e.g., Au3, Au5, Au9, and Au14) or across the vertical direction (as in Au2, Au5, Au19, Au24, and Au28).

Fig. B.3. Skewness (h3) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

While it is not clear in mock images that thin disks flare (Fig. 1), kinematic maps (Figs. B.1 to B.4) show that a disk component defined by fast-rotating, dynamically cool stars flares. While most of these patterns in Figs. B.1 to B.4 were also found in observations, some features may be more visible in results from simulations, where the signal-to-noise is not affecting the penalization of higher moments (here h₃ and h₄) as it happens in observations, biasing toward lower values (Pinna et al. 2019b,a, 2021). In particular, the wealth of structures displayed here in the h₄ maps has not been observed in observations of thick and thin disks in edge-on galaxies.

Fig. B.4. Kurtosis (h4) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

B.2. Stellar-population maps for the full sample of 24 galaxies

We provide here maps of age (Fig. B.5), metallicity (Fig. B.6) and [Mg/Fe] abundance (Fig. B.7) for the full sample of 24 galaxies.

Fig. B.5. Mass-weighted stellar age maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Ages were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

Fig. B.6. Mass-weighted total stellar metallicity [M/H] maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal black dashed lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Metallicities were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

Fig. B.7. Mass-weighted stellar [Mg/Fe]-abundance maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. [Mg/Fe] abundances were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

B.3. Maps of accreted stars for the full sample of 24 galaxies

We show here maps of the mass surface density (Fig. B.8) and mass fraction (Fig. B.9) of accreted stars for the full sample of 24 galaxies.

Fig. B.8. Maps of the mass surface density of accreted stars, in logarithmic scale, of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

Fig. B.9. Maps of the mass fraction of accreted stars of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and the total mass fraction of accreted stars (calculated in the analyzed region, of size 1.6Ropt × 4hscale, Sect. 3.1) are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

B.4. Time evolution

For completeness, we show here the time evolution of all galaxies in our sample that were not shown in Sect. 4.4. Au1 (Fig. B.10) never grew an extended thin metal-rich disk and does not have a clear double-disk structure at z = 0 (Sect. 3.2). Very recently (see the transition between redshift z = 0.36 and z = 0), its global growth was rapidly enhanced by an important merger, which also led to the scattering and distortions of the more metal-rich stars. This galaxy has a mass fraction of accreted stars of ∼46% at z = 0 (Table 2). At early times, Au2 and Au9 formed a relatively metal-poor (thick) disk (which was already in place at z = 1.7), with a metal-rich core that later expanded and formed the thin disk in an inside-out fashion (Fig. B.11 and B.16). The snapshot at z = 1.7 in Fig. B.16 is an example of a failed alignment of the galaxy (see Sect. 4.4). Au3, Au5, Au12 and Au14 (respectively Figs. B.12, B.13, B.18, and B.19) show originally (z = 1.7) a thick structure that grew inside out, while the metal-rich disk was formed only later on. In Au5 and Au12 this clearly happened after a merger between z = 0.62 and 0.79, while at z = 0.36 this thin disk was still disturbed in Au12. Au14 formed it at z ∼ 1, it was disturbed by a merger around z = 0.62, but it had recovered its thin and smooth shape at z = 0.36. Au3 showed a metal-rich core since z = 3.5 (12 Gyr ago), when it still had a spheroidal morphology. Au6 formed a small, thick, metal-poor disk (which had at z = 3.5 the thickness of the thin disk at z = 0 but a much smaller radial extension), which grew very fast in radius and height. It started to form a metal-rich thin disk, inside out, around z = 1 (Fig. B.14). Au7 had a similar but slightly slower initial evolution, with at least two mergers (at z ≲ 1 and z ≲ 0.36), which disturbed the initial disky morphology and prevented the formation of a massive thin disk until very recent times (z < 0.36, Fig. 10). In Au8, the first metal-poor structure was spheroidal and a disky shape appeared only around z ∼ 1 (Fig. B.15). Au10, with a similar appearance to Au3 and Au5 at early times, did not grow much in radius, while it formed a thin metal-rich disk embedded in the thick disk (Fig. B.17).

Fig. B.10. Same as Fig. 10 but for the galaxy Au1.

Fig. B.11. Same as Fig. 10 but for the galaxy Au2.

Fig. B.12. Same as Fig. 10 but for the galaxy Au3.

Fig. B.13. Same as Fig. 10 but for the galaxy Au5.

Fig. B.14. Same as Fig. 10 but for the galaxy Au6.

Fig. B.15. Same as Fig. 10 but for the galaxy Au8.

Fig. B.16. Same as Fig. 10 but for the galaxy Au9.

Fig. B.17. Same as Fig. 10 but for the galaxy Au10.

Fig. B.18. Same as Fig. 10 but for the galaxy Au12.

Fig. B.19. Same as Fig. 10 but for the galaxy Au14.

Fig. B.20. Same as Fig. 10 but for the galaxy Au15.

Fig. B.21. Same as Fig. 10 but for the galaxy Au16.

Fig. B.22. Same as Fig. 10 but for the galaxy Au17.

Fig. B.23. Same as Fig. 10 but for the galaxy Au19.

Fig. B.24. Same as Fig. 10 but for the galaxy Au21.

Fig. B.25. Same as Fig. 10 but for the galaxy Au22.

Fig. B.26. Same as Fig. 10 but for the galaxy Au23.

Fig. B.27. Same as Fig. 10 but for the galaxy Au24.

Fig. B.28. Same as Fig. 10 but for the galaxy Au25.

Fig. B.29. Same as Fig. 10 but for the galaxy Au26.

Fig. B.30. Same as Fig. 10 but for the galaxy Au27.

Fig. B.31. Same as Fig. 10 but for the galaxy Au28.

In summary, AURIGA galaxies generally first formed a metal-poor, small spheroidal or thick-disk structure that then grew inside out and bottom up. They formed a metal-rich structure a bit later on (still very early in some galaxies, and much later in others), which grew inside out as a thin disk. All galaxies experienced the expected chemical enrichment (e.g., Worthey et al. 1992) with time. They started forming the first stars with low metallicity and high [Mg/Fe], and the latter decreased only when and where metallicity increases. In different regions at different snapshots, all galaxies show anticorrelation between [M/H] and [Mg/Fe].

All Tables

All Figures

Fig. 1. Mock edge-on images, in terms of surface brightness in the V band (SB), of our selected sample of 24 Milky Way-like galaxies from AURIGA simulations. For display purposes, we show for all galaxies an area of the same size, of radius 25 kpc and height 10 kpc. Surface brightness was projected into pixels of 0.1 kpc  ×  0.1 kpc size and was Voronoi binned to obtain at least 100 star particles per bin. The region between the two horizontal lines is where the light of the thin-disk component dominates over that of the thick disk (see Sect. 3.2 for details). The region between the two vertical dashed magenta lines is dominated by a central component (the bar in barred galaxies and a classical bulge in non-barred galaxies). Names of the AURIGA halos are indicated on top of each panel.

In the text

Fig. 2. Fit of the vertical surface-brightness (SB) profile of one radial bin of the galaxy Au5. Each small gray circle is the measurement for one Voronoi bin, while larger magenta circles correspond to the median profile of the gray points. The magenta shading indicates the 1σ dispersion of the gray points. Dashed blue and orange lines correspond respectively to the thin- and thick-disk fits, and the solid red line is a combination of the two.

In the text

Fig. 3. Four representative examples of mass-weighted stellar age maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Ages were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

In the text

Fig. 4. Four representative examples of mass-weighted total stellar metallicity [M/H] maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal black dashed lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Metallicities were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

In the text

Fig. 5. Four representative examples of mass-weighted stellar [Mg/Fe]-abundance maps. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. [Mg/Fe] abundances were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin. Maps of the full sample are shown in Appendix B.2.

In the text

Fig. 6. Four representative examples of maps of the mass surface density of accreted stars, in logarithmic scale. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps. Maps of the full sample are shown in Appendix B.3.

In the text

Fig. 7. Four representative examples of maps of the mass fraction of accreted stars. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and the total mass fraction of accreted stars (calculated in the analyzed region, of size 1.6 Ropt × 4 hscale, Sect. 3.1) are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps. Maps of the full sample are shown in Appendix B.3.

In the text

	Fig. 8. Thick-disk versus thin-disk stellar-population properties. From left to right: age, [M/H], and [Mg/Fe] abundance. Points are color coded according to the total mass of the galaxy (Table 1). The one-to-one line is indicated in black and dashed, in each panel. A red “+” symbol indicates solar metallicity in the middle panel.
In the text

	Fig. 9. Thick-disk versus thin-disk stellar-population properties. From left to right: age, [M/H], and [Mg/Fe] abundance. Points are color coded according to the total mass fraction of accreted stars in the analyzed region of the galaxy (Sect. 3.1). The one-to-one line is indicated in black and dashed in each panel. In the middle panel, the red “+” symbol indicates solar metallicity.
In the text

Fig. 10. Time evolution of the galaxy Au7 seen edge on. Surface brightness is shown in the left column, total metallicity [M/H] in the middle column, and [Mg/Fe] abundance in the right column. From top to bottom, the stellar component is shown for seven different snapshots. On top of the middle panels, the redshift and the look-back time (LBT) corresponding to the snapshot shown in that row are indicated. For comparison, horizontal dashed lines separate the regions where the thin or the thick disks dominate at z = 0, while the vertical lines indicate the region where the central component (bar or classical bulge) dominates at z = 0.

In the text

	Fig. 11. Same as Fig. 10 but for the galaxy Au18.
In the text

Fig. 12. Frequency histograms of hz, thin and hz, thick in our AURIGA sample, with 0.1 and 0.5 kpc bins, respectively (solid blue and red lines). These distributions have a weighted average value of 0.6 and 2.6 kpc, respectively, indicated as vertical dashed blue and red lines. Vertical green and magenta lines indicate the average (dashed) and maximum (dotted) thin- and thick-disk values from observations in the Comerón et al. (2018) sample.

In the text

Fig. B.1. Mean-velocity (V) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

Fig. B.2. Velocity-dispersion (σ) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

Fig. B.3. Skewness (h3) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

Fig. B.4. Kurtosis (h4) maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show a central area of the same 50 kpc  ×  20 kpc size for all galaxies. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

Fig. B.5. Mass-weighted stellar age maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Ages were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

In the text

Fig. B.6. Mass-weighted total stellar metallicity [M/H] maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal black dashed lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Metallicities were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

In the text

Fig. B.7. Mass-weighted stellar [Mg/Fe]-abundance maps of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. [Mg/Fe] abundances were integrated and projected into pixels of 0.5 kpc  ×  0.5 kpc size and Voronoi binned to obtain at least 900 star particles per bin.

In the text

Fig. B.8. Maps of the mass surface density of accreted stars, in logarithmic scale, of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and Hubble types from Table 1 are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

Fig. B.9. Maps of the mass fraction of accreted stars of the 24 galaxies in our sample, seen edge on. For display purposes, we show for all galaxies a central area of the same 50 kpc  ×  20 kpc size. Horizontal dashed black lines indicate the regions where the thin disk (within the two lines) and the thick disk dominate (above and below the region between the two lines). Vertical dashed black lines enclose the central region dominated by a bar or a classical bulge. Names of the AURIGA halos and the total mass fraction of accreted stars (calculated in the analyzed region, of size 1.6Ropt × 4hscale, Sect. 3.1) are indicated on top of each panel. Pixels of 0.5 kpc  ×  0.5 kpc size, Voronoi binned to a target number of particles of 900 star particles per bin, were used to plot these maps.

In the text

	Fig. B.10. Same as Fig. 10 but for the galaxy Au1.
In the text

	Fig. B.11. Same as Fig. 10 but for the galaxy Au2.
In the text

	Fig. B.12. Same as Fig. 10 but for the galaxy Au3.
In the text

	Fig. B.13. Same as Fig. 10 but for the galaxy Au5.
In the text

	Fig. B.14. Same as Fig. 10 but for the galaxy Au6.
In the text

	Fig. B.15. Same as Fig. 10 but for the galaxy Au8.
In the text

	Fig. B.16. Same as Fig. 10 but for the galaxy Au9.
In the text

	Fig. B.17. Same as Fig. 10 but for the galaxy Au10.
In the text

	Fig. B.18. Same as Fig. 10 but for the galaxy Au12.
In the text

	Fig. B.19. Same as Fig. 10 but for the galaxy Au14.
In the text

	Fig. B.20. Same as Fig. 10 but for the galaxy Au15.
In the text

	Fig. B.21. Same as Fig. 10 but for the galaxy Au16.
In the text

	Fig. B.22. Same as Fig. 10 but for the galaxy Au17.
In the text

	Fig. B.23. Same as Fig. 10 but for the galaxy Au19.
In the text

	Fig. B.24. Same as Fig. 10 but for the galaxy Au21.
In the text

	Fig. B.25. Same as Fig. 10 but for the galaxy Au22.
In the text

	Fig. B.26. Same as Fig. 10 but for the galaxy Au23.
In the text

	Fig. B.27. Same as Fig. 10 but for the galaxy Au24.
In the text

	Fig. B.28. Same as Fig. 10 but for the galaxy Au25.
In the text

	Fig. B.29. Same as Fig. 10 but for the galaxy Au26.
In the text

	Fig. B.30. Same as Fig. 10 but for the galaxy Au27.
In the text

	Fig. B.31. Same as Fig. 10 but for the galaxy Au28.
In the text