Stellar populations and origin of thick disks in AURIGA simulations

The origin of thick disks and their evolutionary connection with thin disks are still a matter of debate. We provide new insights into this topic by connecting the stellar populations of thick disks at redshift $z=0$ with their past formation and growth, in 24 Milky Way-mass galaxies from the AURIGA zoom-in cosmological simulations. We projected each galaxy edge on, and decomposed it morphologically into two disk components, in order to define geometrically the thin and the thick disks as usually done in observations. We produced age, metallicity and [Mg/Fe] edge-on maps. We quantified the impact of satellite mergers by mapping the distribution of ex-situ stars. Thick disks are on average $\sim 3$~Gyr older, $\sim 0.25$~dex more metal poor and $\sim 0.06$~dex more [Mg/Fe]-enhanced than thin disks. Their average ages range from $\sim 6$ to $\sim 9$~Gyr, metallicities from $\sim -0.15$ to $\sim 0.1$~dex, and [Mg/Fe] from $\sim 0.12$ to $\sim 0.16$~dex. These properties are the result of an early initial in-situ formation, followed by a later growth driven by the combination of direct accretion of stars, some in-situ star formation fueled by mergers, and dynamical heating of stars. The balance between these processes varies from galaxy to galaxy. Mergers play a key role in the mass assembly of thick disks, contributing an average accreted mass fraction of $\sim 22$\% in the analyzed thick-disk dominated regions. In two galaxies, about half of the geometric thick-disk mass was directly accreted. While primordial thick disks form at high redshift in all galaxies, young metal-rich thin disks, with much lower [Mg/Fe] abundances, start to form later but at different times (higher or lower redshift) depending on the galaxy. We conclude that thick disks result from the interplay of external processes with the internal evolution of the galaxy.


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 high detail the mass assembly of disk-dominated galaxies.While only bright thin disks were initially identified in galaxies, fainter thick disks were detected much later.They 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 ⋆ e-mail: francesca.pinna@iac.esdisk, 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, stars in the thick disk were found to have different properties from the ones in the thin disk.Thickdisk stars were associated with old ages, much lower metallicities and higher velocity dispersions than the ones 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 thick disks were enhanced in α elements (Fuhrmann 1998).Old ages, low metallicities and high [Mg/Fe] abundance 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 it was later showed (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 were 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 redshift, during the epoch of gas-rich frequent 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 redshift.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 the ones 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 would be mostly made up (more than 60% of their mass) by the debris of accreted satellites and would have an ex-situ origin (Abadi et al. 2003).Most of the thin disk would have formed during a later quiescent phase with no mergers.
After these scenarios initially came out from numerical studies, 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 the thick-disk dominated regions.Yoachim & Dalcanton (2005, 2008b) extracted the kinematics of the thick and the thin disks of a sample of nine edgeon late-type galaxies.They found a variety of kinematic properties, with high-mass galaxies (with a circular velocity above 120 km s −1 ) showing similar rotation curves in the thick and the thin disks, and low-mass galaxies showing strong differences.The substantial lag and counterrotation 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 deepenough archival observations.For all these reasons, at the moment, there are only 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 VI-MOS at 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 the thin disk supported a fast thick-disk formation at high redshift.Later studies with MUSE, at a much higher resolution, revealed more subtle differences between the thin and the thick disks.Comerón et al. (2016) presented 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 discs at high redshift.Pinna et al. (2019a,b) added [Mg/Fe] maps, as well as starformation histories coupled to 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 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 of accreted satellites.Stellarpopulation 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 inside-out later 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 above-mentioned recent integralfield spectroscopy observations revealed that mergers played a key role in the mass assembly of the thick and the thin disks.However, the relative importance of ex-situ and in-situ components, and therefore the balance between direct accretion of exsitu stars and in-situ star formation, which may 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, solving 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 investigated in cosmological simulations 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 project (Wang et al. 2015), reached a similar conclusion.This idea was also supported by some observational studies (e.g.Scott et al. 2021).VINTER-GATAN, a cosmological zoom simulation of a Milky Way-mass disk galaxy, also supported a similar picture (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 be made up of two components.One in-situ component (with similar properties to the "starburst sequence" shown by An et al. 2022) would have formed during a starburst triggered by the gas-rich merger related to the Gaia-Enceladus Sausage (GES, see also Ciucȃ et al. 2023).A second component would have formed with stars from a pre-existing thinner disk, which 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 disc region of the low-α sequence (Mackereth et al. 2019).An alternative scenario, not involving galaxy mergers, predicts that intense star formation in clumps naturally forms the high-α thick disk, while the low-α thin disk was the result of a more 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 to the abovementioned works.Martig et al. (2014a,b) analyzed the mono-age populations in seven galaxies from Martig et al. (2012), and found 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 thin component.While most of the very old stars were born already dynamically hot, later generations were heated with time, partially also 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 one of them with two disk components.Part of the sample, with lower galaxy mass and thick-disk mass and a more quiescent merger history, shows a continuous and smooth age vertical gradient (mono-age populations with increasing scale height with age).The rest of the sample, with higher galaxy masses, showed a thicker morphology, higher thick-disk mass and mass fractions, and flaring thick disks.Park et al. (2021) also decomposed morphologically the thick from the thin disks in 19 simulated galaxies, and found 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 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 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 oriented to the comparison with (IFS) observations of edge-on external galaxies, allowing to test the pictures drawn by observations and give a better interpretation of observed stellar-population properties, are currently lacking.In this work, based on a geometrical definition of the thick and the thin disks, we use zoom-in magnetohydrodynamical cosmological simulations from the AURIGA project (Grand et al. 2017).Projecting simulations to an edgeon view allows us to perform our analysis in a way similar to 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 show our results in Sect.4, discuss them in Sect. 5 and summarize our conclusions in Sect.6.

Our sample of 24 simulated galaxies
We have 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.

AURIGA zoom-in cosmological simulations
The AURIGA project1 (Grand et al. 2017) includes zoomin magneto-hydrodynamical cosmological simulations of 30 Milky Way-mass late-type galaxies, with stellar masses between 10 10 M ⊙ and slightly more than 10 11 M ⊙ .These were obtained by re-simulating 30 haloes from the parent largest-volume Dark Matter Only EAGLE simulation (Schaye et al. 2015), at intermediate resolution (L100N1504).These haloes were selected to have a virial mass 1 < M 200 /10 12 M ⊙ < 2 at z = 0, defined as the mass contained inside the virial radius defined as R 200 , i.e. the radius enclosing a mean mass volume density of 200 times the critical density for closure.The 30 haloes were also selected to be relatively isolated at redshift z = 0.In particular, they were randomly chosen among the ones in the lowest quartile in τ iso,max .The latter is the maximum value of the tidal isolation parameter τ iso,i = (M 200,i )/(M 200 ) × (R 200 /R i ) 3 , where M 200,i and R i are respectively the virial mass of, and distance to the ith halo in the simulation.
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 of 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 the thick disk (see Sect. 3.2 for details).The region between the two vertical magenta dashed 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.
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 magnetohydrodynamics 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 haloes in the parent simulation.Then, each dark-matter particle was substituted by 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 we use in this work 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 5 M ⊙ and ∼ 5 × 10 4 M ⊙ , respectively (Grand et al. 2017).The comoving gravitational softening length for stellar and high-resolution dark-matter particles was set to 500 h −1 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 density (and thus larger volumes) have larger softening lengths than the ones with higher density.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, kinematic and chemical properties, see e.g.Grand et al. 2016Grand et al. , 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 modelled explicitely, but only as part of a subgrid model (Springel & Hernquist 2003).Gas becomes thermally unstable and forms stars above a density threshold of 0.13 cm −3 .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 a supernova type II (SNII) feedback.In the first case, the star particle represents a single stellar population (SSP) of specific mass, age 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 supernova type Ia (SNIa) and asymptotic giant branch (AGB) stars are also modeled, in terms of metals and mass loss.BHs with a seed mass of 10 5 M ⊙ h −1 were introduced in the position of the densest gas cell in mostmassive haloes (> 5 × 10 10 M ⊙ h −1 ), 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, non-local 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.

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. 2019b,a;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 towards 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 stellarpopulation 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 (Sect.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.2kpc.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 stellarpopulation parameters.The photometric properties of the star particles were obtained by Grand et al. (2017) using stellarpopulation 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 have 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 10 M ⊙ and ∼ 10 11 M ⊙ .

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 the latter (examples in Comerón et al. 2016Comerón et al. , 2018Comerón et al. , 2019;;Pinna et al. 2019b,a;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. 2019b,a;Martig et al. 2021;Sattler et al. 2023).We selected a region with a radius of 0.8R opt and a vertical extension of 2h 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 −2 , while the parameter h scale was calculated as the standard deviation of the vertical positions of stellar particles at 0.5R opt , as defined in García de la Cruz et al. (2021).Thus, the analyzed region has a total size of 1.6R opt × 4h 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.6R opt × 4h 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.

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. 2015Comerón et al. , 2016;;Pinna et al. 2019b,a;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.8R 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.6R opt ×4h scale area of each galaxy into four quadrants with a radial range between R 0,disk and 0.8R 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): 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 ) 2 , 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 −2 and the points obtained averaging less 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 δz −2 tT,i , where δz 2 tT,i 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.
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) = Σ 0 sech 2 (|z|/2h z ).This suggests that we do not have two clear morphologically distinct disk components within the analyzed region in this galaxy.

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. 2019b,a;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 have 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 have 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.

Stellar-population maps
While maps of the kinematic parameters are included in Appendix B.1, we show in Fig. 3 to 5 the mass-weighted stellarpopulation maps for our sample of 24 galaxies.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 Fig. 3 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, 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 andBenítez-Llambay et al. 2018 for other simulations).This flaring extends to the outer regions of the geometrically defined thick disk, which presents a nega-A&A proofs: manuscript no.auriga_thickdisks_aanda tive 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).
In Fig. 4, 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 (Fig. 5).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.

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.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 Fig. 7.In some cases, such as Au12 or Au14, we reach up to 80% of accreted stars in some Voronoi bins in the thickdisk 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 Fig. 6 and 7).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 Fig. 3 and 4 as a younger and more metal-rich structure than its surroundings.It corresponds to slightly lower values of [Mg/Fe] in Fig. 5, although differences with the surroundings are not obvious.
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.6R opt × 4h scale , Sect.3.1).These mass fractions were also indicated on top of each panel of Fig. 7, 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 8 to 2×10 10 M ⊙ in the galaxy (4 × 10 8 to 10 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. 3. 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.Black horizontal 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).Black vertical dashed 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.
Article number, page 9 of 50 Fig. 4. 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.Black horizontal 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).Black vertical dashed 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.
Article number, page 10 of 50 Fig. 5. 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.Black horizontal 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).Black vertical dashed 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.
Article number, page 11 of 50 Fig.6.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.Black horizontal 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).Black vertical dashed 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.
Article number, page 12 of 50 Fig.7. 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.Black horizontal 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).Black vertical dashed 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.6R opt × 4h scale , 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.
Article number, page 13 of 50

Comparison of thick and thin disks
We show in Fig. 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 Fig. 8 and 9).The average [M/H] of thick discs range from subsolar to slightly supersolar values, while thin-disk average [M/H] is always supersolar, higher than the most metal-rich thick disk in the sample.
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 2 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 2 ∼ 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.18).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 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 2 ∼ 0.3).Galaxies with less significant accretion reach lower [Mg/Fe] values in their thin disks.We have 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.

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.2, 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 Fig. 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.2.While we show here only surface brightness, [M/H] and [Mg/Fe], velocity maps (not shown) helped the identification of (rotating-)disk formation.Au7, with a stellar mass of 4.9 × 10 10 M ⊙ at z = 0, had a slower evolution than Au18, which is slightly more massive (with a stellar mass of 8.0 × 10 10 M ⊙ ).
Au7 had formed only 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 black vertical dashed 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.4R 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.6R opt × 4h scale ).It had a stellar mass of about 5.4 × 10 9 M ⊙ at z = 1.7, about 10% of the final galaxy mass at z = 0.No solarmetallicity 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 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 towards 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 Table 2. Columns from left to right: galaxy name, ex-situ mass (mass of accreted stars) in the analyzed region of the galaxy (M accr,gal ), mass fraction of accreted stars in the analyzed region of the galaxy ( f accr,gal ), ex-situ mass in the thin-disk region (M accr,thin ), ex-situ mass fraction in the thin-disk region ( f accr,thin ), ex-situ mass in the thick-disk region (M accr,thick ), ex-situ mass fraction in the thick-disk region ( f accr,thick ).
(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

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 (Fig. B.15 and B.18).However, the primordial thicker components had very different sizes, compared to the final thick disks at z = 0, in different galaxies (see horizontal dashed lines in the figures).Some had already extended thick disks (vertically and radially, e.g.Au2, Au17, Au18), while others were smaller and relatively thinner (Au7) or thicker (e.g.Au3, Au10, Au12), compared to 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 line of sight and that are not possible to disentangle morphologically.These components include stellar haloes, 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 threedimensional 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 exsitu component of thick disks.Within the analyzed thick-disk region, we have disentangled the fractions of stars with thick-disk kinematic properties, with orbital circularities ϵ > 0.5 (following Yi et al. 2023), or associated to 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. (2016Grand et al. ( , 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 starformation histories of thick disks and will discuss further the insitu 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, Au18), started to form a thin metalrich 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.ies 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 re-gions 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 frac- tion 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 green and blue dashed vertical lines).
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 Fig.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 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 to note 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.

Spectroscopic observations of edge-on galaxies
We have used throughout the paper 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 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 metalpoor 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 edgeon 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. 2019b,a), while others show thick disks of similar ages of our youngest thick disks in AURIGA (Kasparova et al. 2016(Kasparova et al. , 2020)).
The fact that young stars strongly flare in AURIGA disks (Fig. 3), 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. 3, 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, either 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 discussion in Pinna et al. 2019b,a).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 supersolar 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 very well in agreement with values in the right panels in Fig. 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 affect the recovery of stellar-population parameters, due e.g. to projection, systematics related to the SSP models and their age, [M/H] or [Mg/Fe] resolution, 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 towards 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 more quantitatively comparable results 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. (2019b,a) 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. 7).On the other hand, the mass density of accreted stars is larger in a central region of few kpc both in all our simulated galaxies (Fig. 6), and in the observed ones.Since no kinematic cut was applied in observation or simulations, stars to the thickdisk dominated region may contain in both cases some contamination from stars with halo kinematics (see Sect. 5.1).

Conclusions
The origin of thick disks and the connection with thin disks are still a matter of debate in spite of the recent progress.Detailed observational studies (in particular from IFS) are so far limited to individual galaxies or very small samples, while simulation studies do not usually provide directly comparable results 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 a Voronoi binning.We fitted vertical surface brightness profiles in the V band with two disk components, allowing us for a geometrical definition of the thick and the 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 [Mg/Fe] enhanced than thin disks, indicating that most of their mass was formed at earlier times from gas 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 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 its growth.Two galaxies had half of their thick disks made up of ex-situ stars.In only four galaxies, satellite debris contributes 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 show a higher mass density).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 may be contaminated, significantly in some cases, from stars with halo kinematic properties, that cannot be disentangled in observations of edge-on external galaxies.The halo contamination is in general larger 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 insitu component, especially when they were relatively recent.-Mergers also cause additional thickening in preexisting disks, scattering cooler stars to higher 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 redshift only in some galaxies.In others, with a recent more active merger history, it 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 will discuss further 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, which will shed some 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  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.Black horizontal 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).Black vertical dashed 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.
Article number, page 24 of 50 Kurtosis (h 4 ) 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.Black horizontal 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).Black vertical dashed 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.
Article number, page 27 of 50

Notes. ( 1 )
Hubble types from Walo-Martín et al. (2021); (2) Stellar mass M * and optical radius R opt (defined as the radius where the B-band surface brightness drops below 25 mag arcsec −2 ) are from Grand et al. (2017); (3) R 0,disk , the minimum radius at which the thin disk dominates over the central component.When the galaxy is barred, this value is the length of the bar taken from Blázquez-Calero et al. (2020), following Walo-Martín et al. (2021).When the galaxy is unbarred, R 0,disk is the approximate radius where the disk starts to dominate over the bulge, according to Fig.4inGrand et al. (2017); (4) The h scale parameter was calculated here as the standard deviation of the vertical position of particles at 0.5R opt , following García de laCruz et al. (2021).

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

Fig. 8 .
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 (Table1).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 .
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.A red "+" symbol indicates solar metallicity in the middle panel.

Fig. 10 .
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.

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

Fig
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.Black horizontal 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).Black vertical dashed lines enclose the central region dominated by a bar or a classical bulge.Hubble types from Table1are 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
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.Black horizontal 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).Black vertical dashed lines enclose the central region dominated by a bar or a classical bulge.Hubble types from Table1are 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.

Table 1 .
Main properties of galaxies in our sample.
TableA.1.Parameters from the morphological decomposition into a thin and a thick disk.Columns from left to right: galaxy name, number of disk components (N disks ) identified in the analyzed 1.6R opt ×4h scale region, distance from the midplane beyond which the thick-disk light dominates over the thin-disk light (z tT ), its uncertainty (δz tT ), thin-disk scale height (h z,thin ), its uncertainty (δh z,thin ), thick-disk scale height (h z,thick ) and its uncertainty (δh z,thick ).Galaxy name N disks z tT (kpc) δz tT (kpc) h z,thin (kpc) δh z,thin (kpc) h z,thick (kpc) δh z,thick (kpc) Au1, with no clear double-disk structure, was fitted with only one disk component and h z,thick and δh z,thick refer to its scale height and uncertainty.
Article number, page 23 of 50