On the Significance of the Thick Disks of Disk Galaxies

Thick disks are a prevalent feature observed in numerous disk galaxies, including our own Milky Way. Their significance has been reported to vary widely, ranging from a few percent to 100% of the disk mass, depending on the galaxy and the measurement method. We use the NewHorizon simulation, which has high spatial and stellar mass resolutions, to investigate the issue of the thick-disk mass fraction. We also use the NewHorizon2 simulation, which was run on the same initial conditions, but additionally traced nine chemical elements. Based on a sample of 27 massive disk galaxies with M * > 1010 M ⊙ in NewHorizon, the contribution of the thick disk was found to be 20% ± 11% in r-band luminosity or 35% ± 15% in mass to the overall galactic disk, which seems in agreement with observational data. The vertical profiles of 0, 22, and 5 galaxies are best fitted by 1, 2, or 3 sech2 components, respectively. The NewHorizon2 data show that the selection of thick-disk stars based on a single [α/Fe] cut is contaminated by stars of different kinematic properties, while missing the bulk of kinematically thick disk stars. Vertical luminosity profile fits recover the key properties of thick disks reasonably well. The majority of stars are born near the galactic midplane with high circularity and get heated with time via fluctuations in the force field. Depending on the star formation and merger histories, galaxies may naturally develop thick disks with significantly different properties.


INTRODUCTION
Well over half of the massive galaxies are said to be "spiral" in morphology, but what is perhaps more striking than the beautiful spiral patterns is their thin disk.Most of the massive spiral galaxies are thin disks in essence.The vertical scale height of the Milky Way (MW) disk is only about 300 pc (Juric et al. 2008), while the radial scale length of the thin disk is an order of magnitude larger (Bland-Hawthorn & Gerhard 2016).
The thin nature of disk galaxies is often attributed to the assembly history of angular momentum from the neighborhood during the disk development: tidal torque theory (Peebles 1969).In addition, a coherent assembly of angular momentum and secular diffusion of the * Email address: starbrown816@yonsei.ac.kr galaxy's orbital structure over time results in a disk galaxy (Binney & Lacey 1988;Pichon et al. 2011;Halle et al. 2018).In this view, the morphology of a galaxy is heavily influenced by the "cosmic web".
A thicker component of the galactic disk was found in the Milky Way (Gilmore & Reid 1983) and other galaxies (e.g., Burstein 1979;Yoachim & Dalcanton 2008;Comerón et al. 2018).The initial star count studies conducted on the Milky Way disk found that the vertical distribution of disk stars showed a break and thus was better fitted by two exponential components rather than one (e.g., Gilmore & Reid 1983).The mass fraction of the "thick disk" was initially derived to be only a couple of percent, and its vertical scale height was measured to be several times that of the thin disk.The MW thick disk was later fitted with sech 2 function (van der Kruit 1981), and the fraction of the thick disk stars was estimated to be about 4% (Juric et al. 2008).Dynami-cal modeling applied to the star count data considering the orbital motions of stars embedded in a dark matter halo suggested 14% for the thick disk mass contribution (Bovy & Rix 2013), which makes the thick disk no longer such a negligible component as was initially hinted.An explanation is required for such prominent structures.
A spectroscopic effort for finding thick disk stars was made in parallel.A flurry of observations has noted that stars in the Solar Neighborhood are bimodal in [α/Fe] abundance (e.g., Bensby & Feltzing 2006;Reddy et al. 2006;Lee et al. 2011).A naive projection of the [α/Fe] bimodality onto kinematic properties leads to a suggestion that the thick disk fraction is roughly a third within 1 kpc from the galactic mid-plane (Lee et al. 2011).This makes the thick disk more significant than previously believed.Does the [α/Fe] bimodality indicate a thin-tothick disk mass ratio?Why does it suggest a different mass fraction for the thick disk from vertical profile fit studies?
Thick disks are found in other galaxies, too.Yoachim & Dalcanton (2006) inspected 34 edge-on disk galaxies with a wide range of mass.They measured the vertical scale heights and luminosity fractions of their thick disks.The vertical scale height ratios between the thick and thin disks were found to be 2-3.The luminosity fraction of the thick disk ranges widely between 0 and 100 %, and the Milky Way values are compatible with these data.
The origin of thick disks has been a subject of debate.The stars in spiral galaxies were probably born predominantly in the molecular clouds within thin disks.Such a pre-existing disk of stars would be kinematically heated through various sources of perturbation such as minor mergers (Quinn et al. 1993;Kazantzidis et al. 2008), spiral arms and bars (Sellwood & Carlberg 1984;Roškar et al. 2013;Vera-Ciro et al. 2014), or giant molecular clouds (Spitzer & Schwarzschild 1951).The accreted stars from satellite galaxies are likely to be kinematically hotter occupying a thicker region of the disk (Abadi et al. 2003)."Resonant thickening" of the disk by small satellites has also been suggested (Binney & Lacey 1988;Sellwood et al. 1998).Accretion of gas clouds may induce star formation away from the galactic disk midplane and thus occupying a thicker region of the disk (Norris & Ryan 1991).Alternatively, gas-rich galaxy mergers could result in active star formation away from the thin disk.The stars formed this way would have a higher probability of occupying a thicker region of the disk (Brook et al. 2004;Agertz et al. 2021).
It is clear that the issue of thick disks is outstanding.How significant is it?Why is its significance in the Milky Way estimated to differ according to different methods?What are the main mechanisms behind it?Can we understand the range of thick disk contributions in other galaxies?These are the questions we want to answer through this investigation.
To answer these questions, we use the NewHorizon (NH) simulation (Dubois et al. 2021).We also use the NewHorizon2 simulation (NH2,Yi et al. in prep) which is a lower-resolution twin that started from the same initial conditions but has a more extensive chemical followup (see Section 2).The combination of the high resolution of NH and detailed chemical information of NH2 makes this investigation viable.

DATA
The main data used was the NH simulation (Dubois et al. 2021).NH is a medium-sized (20 Mpc diameter), cosmological hydrodynamic simulation of galaxy formation which was run with RAMSES (Teyssier 2002).It resolves dense regions spatially down to 34 pc (at z = 0) which allows for sampling the vertical structure of a Milky Way-size disk in 10 bins or more, which is ideal for tackling the issues of the galactic disk structure.Its stellar mass resolution is 10 4 M ⊙ .Therefore, typical Milky Way-mass galaxies are realized by millions of star particles.The temporal resolution is also exceptionally high: 867 snapshots were stored with a mean timestep of 15 Myr.This allows a precise inspection of the motions of the star particles and gas.The NH volume contains 48 galaxies of stellar mass above 10 10 M ⊙ and so provides good statistics for drawing a general picture on thick disks.Because of its relatively small volume, only some limited comparison against observations is possible.It has been shown that NH reproduces some of the key properties of galaxies: the stellar mass function of galaxies, surface brightness-stellar mass relation, star formation rate density evolution, Kennicutt-Schmidt relation, galaxy mass-size relation, and mass-metallicity relation to name a few.Readers are referred to Dubois et al. (2021) for a detailed description of NH.
We also use the NH2 simulation.It starts from the same initial conditions as NH but has run with a factor of two coarser spatial resolutions (the best resolution of 68 pc) allowing for quicker execution than NH.There are other minor differences between NH and NH2.For example, NH2 has a 50% higher stellar feedback efficiency in the hope of matching the empirical stellar to halo mass relation better.The most significant advantage of NH2 is that it traces nine chemical elements (H, D, C, N, O, Mg, Si, S, and Fe) and provides useful chemical information for star formation history studies.This allows us to enhance the understanding we earlier had on thick discs based on the NH simulation alone (Park et al. 2021).In this investigation, we use the NH data for most analyses.However, when the chemical information is inspected, we use the NH2 galaxies.We do not mix NH and NH2 for the same galaxy because, although NH and NH2 start from the same initial conditions, a few processes adopt a random choice, e.g., star formation (Rasera & Teyssier 2006;Kimm et al. 2017), and thus cause a difference in the final properties of the galaxy.
For the chemical evolution of galaxies, we take into account the contributions from stellar winds, as well as SN Ia and II.To calculate the yields from a simple stellar population (SSP) for each stellar particle, we use the Starburst99 model (Leitherer et al. 1999), assuming a Chabrier initial mass function (Chabrier 2005).The Geneva stellar wind model (Schaller et al. 1992;Maeder & Meynet 2000) is used to calculate the tabulated chemical yields from the stellar wind as a function of the age and metallicity of the SSP (see Leitherer et al. 2014 for more details).Therefore, each stellar particle is involved in chemical enrichment based on its evolving age and metallicity.SN II explosions from 8-50 M ⊙ stars are responsible for the chemical elements returned by each SSP.The SN II yields are based on Kobayashi et al. (2006), where the yields from intermediate-mass stars (8-13 M ⊙ ) are taken from Woosley & Weaver (1995).The chemical yields of SN Ia are obtained from Iwamoto et al. (1999).To estimate the SN Ia frequency for each stellar particle, we adopt the delay time distribution model described in Maoz et al. (2012).This model assumes a power-law decline of the SN Ia frequency over time (∝ t −1 ), which allows us to calculate the SN Ia event rate for a given age of the stellar particle.Based on the normalization values provided in Maoz et al. (2012), each stellar mass particle of 10 4 M ⊙ is expected to have at most 13 SN Ia explosions during its lifetime.
Our galaxy mass threshold (M ≥ 10 10 M ⊙ ) allows at least one million star particles in each galaxy.Of the 48 NH galaxies, we remove 3 galaxies that are severely disturbed morphologically at z = 0.17.We finally have 27 disk galaxies of M * ≥ 10 10 M ⊙ with v/σ ≥ 1 for stars in the last snapshot of z = 0.17, where v/σ stands for the ratio between the rotation speed and velocity dispersion.Similarly, we have selected 10 massive disk galaxies from the NH2 simulation.We have fewer galaxies in the final NH2 sample than in NH primarily for two reasons.First, the factor of two coarser spatial resolution reduced v/σ.Besides, its enhanced stellar feedback caused the final stellar mass to be smaller in NH2.We present the images of the sample galaxies from NH and NH2 in Figure 1 and Figures A1, A2 in Appendix A.

Profile Sampling
With its 34 pc "best" spatial resolution of NH in the adaptive mesh refinement scheme, we effectively resolve the disks of massive galaxies.In addition, its mass resolution of 10 4 M ⊙ for star particles is also essential for dynamically realizing the galactic disk in numerical simulations.Only when we could reproduce the thin disk in the first place would our attempt to understand the thick disk be justified.
To measure the structural properties of a disk, we need to choose a location in the disk that is representative of the disk.We first measure the scale length of the disk, R d , using all the star particles inside R 90 , where R 90 stands for the radius within which 90% of the stellar mass resides.We found a linear relation between R d and the galaxy's stellar mass, M * .We use R d based on this fit, and we measure the thin and thick disk contributions at the radial distance from the galaxy center 2.0±0.1Rd and the vertical distance from the mid-plane |z| ≤ 3 kpc using the sech 2 function (e.g., Yoachim & Dalcanton 2008;Comerón et al. 2011;Park et al. 2021).
We attempt to fit the vertical mass and luminosity profiles of NH galaxies by allowing up to 3 components: where i = 1, 2, 3 are for the "thin", "thick", and "thicker" disk components, and where ρ i and z i are the density and vertical scale height of each component.We use the Bayesian information criteria (BIC) to choose the best-fit model, where the BIC is a measurement that uses both the log-likelihood of the fit (the goodness of the fit) and the number of parameters.As we use more parameters in statistical analysis, the BIC value will be penalized more.Based on the BIC analysis, we found that 0, 22, and 5 galaxies are best fitted by 1, 2, or 3 vertical components, respectively.We will discuss this later in Section 4. We provide the BIC-selected fits for a few massive NH galaxies in Appendix (Figure A3 & Figure A4).
Figure 2 shows the sample two-component fits to the profiles of the sample galaxy presented in Figure 1.Panel-(a) shows the radial surface brightness profile fit to the face-on images using the Sersic profiles (Sérsic 1963).A simple combination of a disk with a Sersic index n = 1 and a spheroid with n = 4 reproduces the radial profile reasonably.The two lower panels show the vertical profile fits in mass and luminosity using the sech 2 function.While they appear somewhat different from each other, the mass and luminosity profiles are equally well-fitted by a model with two components that could be called the "thin" and "thick" disks.The disk properties of this model galaxy seem compatible with those of the observed galaxies (see below).
We now restrict ourselves to a two-component scenario to make our results directly comparable to previous studies which were predominantly based on the assumption of two disk components.In other words, we present the two-component fits to the vertical profiles of the NH galaxies.In the case of galaxies for which 3component fits were selected as the best by the BIC, the mass fractions of the third components are very small.
Figure 3 shows the vertical scale height ratios between the thin and thick disk components and the luminosity fractions of the thick disks of the NH galaxies com-pared with the observed data (Yoachim & Dalcanton 2006;Comerón et al. 2018;Martínez-Lombilla & Knapen 2019).This figure is an updated version of Figure 4 of Park et al. (2021).It is based on the galaxies at z = 0.17 and on the measurements at a new position on the disk mentioned above.The thick disk properties of the NH galaxies are in reasonable agreement with the observed data, particularly those of Comerón et al. (2018).There is a hint of a negative mass (V circ ) dependence in the luminosity fraction, in Panel-(b), when the observed data of Comerón et al. (2018) and the NH galaxies are considered, where the circular velocity V circ has been measured from the stellar particles inside R 90 in case of the NH simulation.
The thick disk luminosity fraction was measured by integrating over the range of 0.5-2.0Rd .The mean value of the integrated luminosity fractions of the thick disks of the NH galaxies is 34 ± 10%, where the errors are the standard deviation among galaxies.When we fit the mass profiles as in the bottom panel of Figure 2, the contribution of the thick disks becomes 54±10% instead.This indicates that the thick disks play a significant role in the overall mass distribution of the galaxies.
In summary, the NH simulation resolves the thin disks of spiral galaxies.The vertical distribution of stars indicates the presence of thicker disk components in all of our sample galaxies.When we consider a two-component assumption consisting of thin and thick disks, the properties of the thick disks in the NH galax-ies show good agreement with the observed data.Thick disks contribute roughly one third of the disk in r-band luminosity or half in mass.Thick disks are not minor.

Spatial Sampling
We would like to study the properties of thin and thick disk stars and their origins.The Solar Neighborhood is composed of both thin and thick disk stars, and it is not trivial to know which star belongs to which component.One way of evading this difficulty, at least in part, is to sample the stars spatially.Based on the definition of thin and thick disks, it is natural to assume that the fraction of thick disk stars decreases with the vertical distance from the disk's mid-plane.We fit the vertical r-band brightness profiles of the NH galaxies using two sech 2 components in the cylindrical region at 2 ± 0.1R d .We consider the stars in the mid-plane at vertical distance |z|/z thick ≤ 0.5 "spatially-thin" disk stars because they are dominated by the thin-disk component according to the two-component fit.We consider the stars away from the mid-plane "spatially-thick" disk stars based on the same arguments.We select two locations to sample thick disk stars: 1.5 < |z|/z thick < 2.5 and 5 < |z|/z thick < 6, where z thick ≈ 1 kpc.An advantage of this sampling method is that it can be easily applied to observed edge-on galaxies.Park et al. (2021) showed from their sample of 18 galaxies that there is little difference in terms of the birth position between the spatially-thin and -thick disk stars.We present in Figure 4 the birth positions of the spatially-thin andthick disk stars of the 27 NH galaxies and confirm their earlier report.Both axes are normalized to the scale values measured in each radial bin in the final snapshot.Both spatially-thin and -thick disk stars are born near the galactic mid-plane within 1-2 vertical scale heights of the thin disk in the final snapshot (|z| < z z=0.17 thin ), where the mean value of z z=0.17 thin is roughly 0.2 kpc.Spatial sampling is useful for learning the properties of thin and thick disk stars, provided that the contamination is not severe.However, it does not provide us with the mass contribution of thick disk stars.

Chemical Sampling
As we mentioned in Section 1, the mass fraction of the Milky Way thick disk is estimated to be 33% when [α/Fe] is used as the sole classifier, while star counts suggest much lower values.We would like to check what the simulation models suggest.In the NH2 simulation, we monitor the chemical evolution of galaxies in multiple elements including alpha elements (e.g., O, and Mg).We derive [α/Fe] based on the formula: [α/Fe] = [(O+Mg+Si)/3Fe], where the solar abundance, [O/H] ⊙ , [Mg/H] ⊙ , [Si/H] ⊙ , and [Fe/H] ⊙ are adopted from Asplund et al. (2009).We present the density distribution of the chemical properties of the star particles of one sample galaxy in Figure 5. Panel-(a) shows the [α/Fe] vs. [Fe/H] plane.Overall, there is a reverse correlation between these properties.This is naturally expected for a largely "closed-box" system, considering the chemical yields of SN Ia and II (Tinsley 1980).The majority of stars have high values of [Fe/H], and as can be seen later, they are the stars that belong to the disk.The highest density peak has [Fe/H] ≈ 0.1 and [α/Fe] ≈ 0.03 which roughly agree with the values of the disk stars of the Milky Way (e.g., Bensby et al. 2014).
Figure 5   The apparent bimodality of the Milky Way disk in alpha abundance is often attributed to two separate populations of disk stars and even to two distinct origins.Figure 5-(c) shows the age distribution of the same stars shown in Panel-(b).It is evident that the alpha abundance distribution is directly linked to the age distribution and hence the star formation history.In the case of this particular galaxy, there was a minor merger with a mass ratio of 5:1 around 6.2 Gyr ago (in coalescence time), which caused a starburst with rough properties of [Fe/H] ≈ −0.2, and [α/Fe] ≈ 0.09.If we consider this model galaxy as an analog of the Milky Way, the bimodality of the alpha abundance found in the Milky Way disk could be interpreted as a result of a galaxy merger that happened in the past.This explanation was already suggested by Agertz et al. (2021) based on the Vintergatan simulation.However, the reality is more complex.We will discuss in the next section why a high alpha peak cannot be interpreted as the result of a merger event alone.
The case of the NH2 galaxy shown in Figure 5 has a bimodal distribution similar to that of the Milky Way, but the exact value that divides the bimodality is different: Milky Way's cut appears to occur at [α/Fe]=0.2while that of the NH2 galaxy does at 0.06.The detailed shape of the bimodality (or multimodality) is specific to the unique star formation history of the galaxy in question.We provide the [α/Fe]-[Fe/H]-age maps (just Panels b and c) for the NH2 galaxies in Appendix (Fig- ures A5, A6).One can see that galaxies show a wide variety of the [α/Fe] distribution.
Figure 6 shows the stellar density distribution of the two components based on a cut of [α/Fe] = 0.06.For consistency in terms of nomenclature, we call the stars with [α/Fe] ≤ 0.06 "chemically-thin" and [α/Fe] > 0.06 "chemically-thick".Panel-(a) shows the chemically-thin disk stars and Panel-(b) shows the chemically-thick disk stars.Chemically-thick disk stars appear to occupy a thicker region of the disk compared to chemically-thin disk stars.Cursory inspection of this figure suggests that [α/Fe] is effective for finding thin and thick-disk star candidates.In the next section, we further inspect this finding in connection with kinematic information.The fraction of chemically-thick disk stars near the Solar Circle is 34% in mass (or 8% in luminosity) in this galaxy, compared to 33% in the MW according to the SDSS observations.The proximity between the two values was most likely achieved by chance.

Kinematic Sampling
An important advantage of using simulation data is that we have detailed kinematic information for the galaxies.Thanks to the high mass, spatial, and temporal resolutions of NH, each galaxy is resolved by millions of star particles with detailed orbital information.Recent studies (Du et al. 2019(Du et al. , 2020;;Jang et al. 2022) have utilized the Gaussian Mixture Model (GMM), which is an unsupervised machine learning technique that finds a given number of groups or kinematic components in the kinematic phase space.It uses three quantities: total energy e, angular momentum along the rotational axis J z , and the remaining angular momentum J p , where J p = J tot − J z .
Figure 7 shows a sample exercise of GMM on the galaxy in Figure 1.The diagram shows the energy e normalized by the maximum energy, and J z normalized by the maximum angular momentum that a particle can have for a specific binding energy e, J circ (e).The former roughly indicates the mean distance from the galaxy center; the smaller, the closer.The latter is also known as circularity (ϵ ≡ J z /J circ (e)), where 1 indicates a perfectly-circular rotation with respect to the rotational axis, −1 denotes a perfectly-counter-rotating motion, and 0 represents dispersion dominant motion.We also note that circularity departure from 1 does not always indicate vertical heating because the radial migration of a star confined to a razor-thin disk may change its circularity without increasing its vertical displacement.However, we found that such cases are rare in our simulations.
In this exercise, we identify the five most significant components by mass contribution, and the numerals in the diagram show the mass rank.Components 1, 2, and 4 have ε > 0.5, where ε denotes the mean value of ϵ of all the members of each component detected.Therefore, they follow the coherent rotational motion.The other components (3 and 5) are centered around ε ≈ 0, meaning that they do not follow the coherent rotation of the galaxy.By using a cut in circularity, one may assign the components to the thin disk (ε > 0.85) or the thick disk 0.5 < ε ≤ 0.85.In this segmentation, the thin disc is therefore defined to be mostly made of quasi-circular inplane orbits, while the thick component corresponds to orbits with higher vertical or radial velocity dispersion.Figure 8 shows their color maps which clearly have distinct spatial distributions.Components 1, 2, and 4 have disky distributions, whereas the other components have dispersed distributions.Kinematically-identified components appear to recover the morphology of separate components closely.While we detect five components in this exercise, it is possible to explore an arbitrary number of components, in which case even minor features such as the stellar halo and tidal debris could be detected.We consider the information derived this way "ground truth" in this investigation, but it should be noted that the classification of thin or thick disks based on a single cut in circularity is somewhat arbitrary regardless of whether it is through GMM or in any other method.For reference, we provide the GMM detection maps of additional NH galaxies in Appendix (Figure A7 and Figure A8).
Figure 9 shows the star particles in the plane of alpha abundance and circularity.The top panel shows the same galaxy shown in Figures 5 & 6.As shown in Figure 5, a bimodal distribution is visible in [α/Fe].A criterion of [α/Fe]= 0.06 roughly separates the bimodality in this galaxy.The low [α/Fe] clump is dominantly composed of stars with high circularity.Therefore, such a simple cut based on [α/Fe] would effectively select "kinematically-thin" disk stars (e.g., ϵ > 0.85).However, the situation for the thick disk selection is not so straightforward.Chemically-thick disk stars ([α/Fe]> 0.06) are not dominated by "kinematically-thick" disk stars (e.g., 0.5 < ϵ ≤ 0.85) at all.Panel-(a) shows the stars virtually in the whole galaxy at 2±0.5R d and |z| ≤ 15 kpc.Kinematically-thick disk stars (i.e., B=24.9%) take up only 47.6% of the chemically-thick disk stars.The mass fraction of the chemically-thick disk stars (A+B+C) is 52.4%, while that of the kinematicallythick disk stars (B+E) is 41.0%.If we measure these values near the Solar Circle (e.g., at 2 ± 0.5R d and |z| ≤ 1 kpc), they become 35.0%(chemically thick) and 25.3% (kinematically thick).The use of [α/Fe] as a thick disk indicator substantially overestimates the thick disk fraction.
Figure 9-(b) shows the case of another galaxy.While the galaxy in Panel-(a) shows two peaks in [α/Fe], the one in Panel-(b) shows three or even more peaks.While the lowest-α peak corresponds to the young thin disk, the other peaks are the results of older stellar populations.The two peaks at [α/Fe] ≈ 0.04 and 0.08 are associated with unusually-high star formation events that occurred around 4 and 6 Gyr ago, respectively.The thick disk of this galaxy has significantly different properties from the galaxy in Panel-(a), and likewise, a variety of thick disk properties are possible in different galaxies.This illustrates nicely how chemistry allows us to trace back bursts of star formation which then diffuse away in orbital space.Stellar chemical information is the time capsule of Galactic archaeology.We provide [α/Fe] vs. circularity diagrams for additional NH2 galaxies in Appendix (Figure A9).
We check the correspondence between the results from different sampling methods.Figure 10 shows a comparison for the Solar Circle region (2 ± 0.5R d and |z| < 15 kpc) of a sample NH2 galaxy.We here adopt a simple criterion based on circularity instead of using the more elaborate GMM technique: 0.5 < ϵ ≤ 0.85 for thick disk stars as illustrated in Panel-(a).The [α/Fe] distribution based on this criterion is shown in Panel-(b).As shown in Figure 6, the low-[α/Fe] peak traces kinematically-thin disk stars well, whereas the high-[α/Fe] peak is heavily contaminated by kinematically-  (d) shows the age distributions.The kinematically-thin disk is predominantly young, whereas the thick disk is composed of stars with a variety of ages.Their vertical distribution is shown in Panel-(e).In this galaxy, the kinematically-thin disk is dominant all the way out to 4 kpc from the mid-plane; thus one has to go beyond 4 kpc to sample thick-disk stars.This is why we presented such distant stars from the galactic plane (red symbols) in Figure 4. Panel-(f) shows a two-component vertical profile fit for comparison.The thick disk fraction is substantially different from that in Panel-(e) in this galaxy, which explains the large scatter in the luminosity fraction of the thick disk stars (Figure 11).This highlights the issue regarding the validity of the profile fit and in return the spatial sampling discussed in Section 3.1 and Section 3.2.
Figure 11 shows the comparison of the key properties of thick disks between the vertical luminosity profile fits and the kinematic detection based on the 0.5 < ϵ ≤ 0.85  shows the thick disk luminosity fraction.The agreement is somewhat better; the thick disk luminosity fraction is 28 ± 17% (or 43 ± 14% in mass).When we integrate over 0.5-2 R d , the thick disk fraction becomes 34 ± 15% in r-band luminosity and 48 ± 13% in mass.The scatter is large, partly due to the fact that the galaxies are in different dynamical stages.Besides, the classification of the thick disk depends on the selection criterion, while there is no clear ground that divides thin and thick disks.The origin of the offset between the two estimates is not obvious to guess at the moment except that we are tempted to say that the profile fits do not recover the true kinematic properties fully.
Figure 12 shows the circularity of all stars in a galaxy when they were born and at present.Stars are typically born in highly circular orbits.The 1σ contour shows that 68% of the stars are born with ϵ ≳ 0.75.This is because star formation in typical massive disk galaxies occurs predominantly within thin gas disks.The circularity distribution of the same stars is much more dispersed at present (i.e., z = 0.17).The 1σ contour extends down to ϵ ≈ 0.4.This summarizes the general pattern of kinematical heating which corresponds to a slow diffusion of the orbital structure driven by force fluctuation along the unperturbed trajectories.In orbital space, this diffusion can in principle operate along three directions, radial migration or "churning" (change of guiding center), "blurring" (increase in radial eccentricity), and "heating" (increase in the vertical oscillation of stars).In Figure 12 we look for simplicity at one projection of these, focussing on a decrease in circularity.This has also been demonstrated in Figure 4 based on the spatial/geometric sampling of thin and thick-disk stars.We provide the circularity at present and at birth for additional NH galaxies in Appendix (Figure A10).
Figure 13 shows the mean change of the circularity of the star particles binned by their birth circularity.The top three panels show the values for three different galaxies, and the last panel shows the mean values for the whole sample of 27 NH galaxies.The two vertical lines show the rough demarcation for the thin and thick disk components.The slanted line shows the reference: a circularity change below this line will make the star particle kinematically thick (ϵ < 0.85).The first three panels highlight that the changes in the kinematic properties vary substantially from galaxy to galaxy, and it is important to have a large sample to obtain a general picture.We then discuss the mean properties shown in Panel-(d).It is useful to remember that most of the stars of the NH disk galaxies are born as kinematically thin, as illustrated in Figure 12.The circularity changes of the stars born as kinematically-thin (four points of ϵ ≥ 0.85) are all negative, meaning that most of them become kinematically-thicker with time.The magnitude of the change is so large that all points are below the red reference line, meaning that a large fraction of the thinborn stars becomes thick in the last snapshot.Many of kinematically-hot stars (e.g., ϵ ≲ 0.35) become kinematically colder, probably due to the resonance with the galactic rotation.However, it should be noted that the mass fraction of stars born with such low circularity values is very small (see Figure 12).
The (ensemble average) stack in Panel-(d) is of clear dynamical interest as it corresponds to a (finite time) proxy of a drift coefficient for circularity.Such a coefficient can be predicted from the first principle within the context of kinetic theory (Binney & Lacey 1988).Kinetic theory captures the fact that environmentally or internally driven force fluctuations operate on stars The comparison between the circularity at birth and the circularity today for all the stars of one NH galaxy.A one-to-one correlation is shown by the diagonal line.The black curves show the 0.5, 1.0, and 1.5-σ contours.The distribution is heavily biased.Most stars are born with a large value of circularity on a thin rotating disk and are gradually heated to show a much wider distribution in the end.
moving to zeroth order along their unperturbed orbit, and induce a slow resonant distortion of their trajectory: the mean orbital structure diffuses on secular timescales.Eventually one would need to compute all drift (vector) and (tensor) diffusion coefficients, which is beyond the scope of this paper.
Figure 14 shows the birthplaces of the kinematicallythin and -thick disk stars sampled by the ϵ = 0.85 cut (between thin and thick).Panel-(a) shows the case of the NH2 galaxy with [α/Fe] bimodality, which was presented in Figures 5 and 6.We recall that the Vintergatan galaxy suggested that its thick disk stars formed closer to the galactic center than its thin disk counterpart (Agertz et al. 2021).In the case of NH2 galaxy, the kinematically-thick disk stars were born systematically closer to the galactic center and farther from the midplane than the kinematically-thin disk stars; but the difference is extremely small.Although it is not shown here, when we inspect the stars that were born between 5 and 7 Gyr of lookback time during the minor merger as we mentioned in Section 3.3, they were born near the galactic mid-plane, too.We suspect that the differences in the details of the merger geometry and galaxy properties caused the differences between the Vitergatan and the NH2 simulations.Panel-(b) shows the mean properties of 10 NH2 disk galaxies.The birth positions of thin and thick disk stars are almost indistinguishable.Whether we use GMM (not shown here) or single-ϵ cut to select thin and thick disk stars does not make a notable difference in this diagram.
We summarise this section as follows.We detected kinematic groups/components based on the kinematic properties of the star particles in the simulation.The use of reasonable criteria for energy and angular momentum successfully detected thin and thick disk particles; however, the separation remained somewhat arbitrary.If we use ϵ = 0.85 as one of the critical conditions, the mean value of the thick disk fraction in r-band luminosity is 34 ± 15% or 48 ± 13% in mass, which is consistent with the observations.Orbital properties are not well traced by the alpha abundance.Kinematic information, such as energy and angular momentum, is useful for studying the physical processes that generated the orbital properties of stars using theoretical/simulation models.However, it is not trivial to measure them observationally.Fortunately, the vertical luminosity profile fits provide key kinematic properties of thick disks reasonably closely.Most stars are born near the galactic plane and get dynamically heated with time.Thick disks are naturally generated as an ensemble of stars that are heated for different lengths of time after their birth.

CONCLUSIONS
We give a brief summary first.We used the NewHorizon simulation in this investigation.With its high resolutions in space and mass, we resolved the thin disks of massive spiral galaxies, which makes our investigation of thick disks possible.NH also has a high temporal resolution which allows accurate measurements of the kinematic properties of star particles.In addition, we used the NewHorizon2 simulation which further provides the evolution of nine chemical elements.This allows a comparison between the results from different approaches that are used to estimate the significance of thick disks.To ensure that the models reasonably resolve the thin disks, we used only the massive (M * ≥ 10 10 M ⊙ ).In addition, we used a stellar v/σ > 1 cut to select disk galaxies.We have 27 galaxies from NH and 10 galaxies from NH2 in our sample.
We first performed vertical profile fitting in r-band using up to three sech 2 components with different scale heights.Based on the BIC statistical assessment, 0, 22, and 5 galaxies were found to prefer one, two, or three components, respectively.This may have a profound implication for the study of thick disks.The majority of the stars in spiral galaxies are probably born within a thin cold disk and get gradually heated with time through many processes.If there were n rather distinctly more active periods of star formation in the galactic disk either via secular regulation, cyclic gas inflow, or episodic merger events, we may find n disk components with different scale heights, not necessarily because the n components were born with a different scale height but more likely because they were kinematically heated for different lengths of time.Therefore, it is not surprising to find more than two disk components in some galaxies.This has indeed been addressed in the observational study of Comerón et al. (2018).The vertical scale heights of the thin and thick disks and the luminosity contribution of the thick disks of the NH model galaxies were found to be in reasonable agreement with the observed data of external edge-on galaxies.
Assuming that the vertical profile fit gives reliable information regarding the distribution of thin and thick disk stars, we first attempted to sample thick disk stars based on position.We sampled spatially-thin disk stars near the galactic mid-plane and spatially-thick disk stars away from the mid-plane.A strong advantage of this sampling technique is that it could be compared with the observed data of external galaxies easily.A disadvantage on the other hand is that both spatially-thin and -thick disk stars sampled this way are bound to be contaminated by each other.We checked where these stars were born initially.We confirm the earlier result of Park et al. (2021) that both thin and thick disk stars were born near the mid-plane.Their birthplaces were almost indistinguishable.
We then sampled thin and thick disk stars based on alpha abundance, [α/Fe].Some observational studies have suggested a connection between the kinematic properties and [α/Fe] properties of disk stars in the case of the Milky Way.The bimodality in [α/Fe] is often considered as evidence for the distinct origin of the thick disk from the thin disk.Besides, a naive projection of the [α/Fe] bimodality would indicate a substantially higher prominence of the thick disk than that suggested by profile fit studies.The NH2 galaxies showed a variety in the pattern of [α/Fe] distribution, while there usually was a good correlation between age and [α/Fe].The NH simulation (and thus NH2 as well) is for an average field environment.The majority of its galaxies have not undergone major mergers since they developed a galactic scale disk or since z ≈ 2. In such an effectively "closed box" of chemical evolution testbed, [α/Fe] by and large monotonically decreases with time.Thus, the [α/Fe] bimodality, more specifically the high [α/Fe] peak hints for a separate event of star formation from the main star formation history of the Milky Way disk.We found a galaxy in our sample that resembles the Milky Way in terms of [α/Fe] bimodality.The high [α/Fe] peak stars were indeed born during the high star formation phase caused by a gas-rich galaxy merger accretion 6 Gyr ago.The high [α/Fe] peak found in the Milky Way disk may have had a similar origin.
In the NH2 galaxy with an [α/Fe] bimodality, the stars in the low and high-[α/Fe] peaks roughly appeared to be thin and thick in spatial distribution at first glance.But this correspondence was found to be weak by the following kinematic analysis.The low-[α/Fe]-peak stars were indeed dominantly kinematically-thin disk stars following the coherent rotation of the disk.However, the high-[α/Fe]-peak stars were not just kinematically-thick disk stars but were significantly contaminated by the stars of other kinematic properties (such as thin disk stars or dispersion-dominant stars).
We used the unsupervised machine learning technique, Gaussian Mixture Model, to detect kinematical groups of stars in each galaxy using energy and angular momentum.It was reasonably straightforward to identify the most significant components in the simulated galaxy.Whether GMM-detected components are "distinct" is debatable; yet, it is probably more reliable than applying a single value cut of any kind.When we classified the particles using reasonable conditions of kinematic properties, we found good correspondence between the kinematic decomposition and profile-fit decomposition of thin and thick disks.Kinematically-selected thin-disk stars were usually low in [α/Fe] and young in age, hence good correspondence.However, the thick-disk counterpart did not represent any single dominant group of stars in age or chemical properties.This is because, while today's kinematically-thin disk stars are predominantly young stars bearing similar chemical properties, the thick disk counterpart is a cumulative population of stars that were born at different epochs with different chemical backgrounds.
Applying an arbitrary yet reasonable criterion to the kinematic data of the 27 NH galaxies indicated that 34± 15% in r-band luminosity or 48±13% in mass of the disk was found to be kinematically hot.These measurements are in reasonable agreement with the observational data of edge-on disk galaxies and the Milky Way.
The presence of thick disks is probably not surprising after all.Star formation in a disk galaxy that is not involved in a significant galaxy interaction probably occurs predominantly within the thin cold gas disk.So the youngest stars likely share similar kinematic properties to the galactic gas disk, rotating fast.The stars are kinematically heated by force fluctuations along their orbit.Perturbation causes an increase of entropy and thus must build up with events and with time.If the sources of perturbation were numerous and effectively smooth, disk thickening would likely be smooth.If the dominant source of perturbation was n dramatic events such as galaxy mergers, their signature will be n discrete thick disk populations.If only one important merger happened in the detectable past, it may show just one thick disk component, which could be the case of the Milky Way disk.However, this does not necessarily mean that the merger was the dominant contributor to the Milky Way thick disk.It is probably more likely that its added contribution made the thick disk more conspicuous.Depending on the detailed nature of the perturbation events, a galaxy may have one dominant thin disk, two disks, or n disk components.
Eventually one should aim to ensemble average orbital diffusion in 3+1 D orbital and chemical space, in the spirit of the bottom panel of Figure 13, so as to predict the mean effect of such mixture of cosmic accretion event.This could then be compared to inhomogenous kinetic theory (Fouvry et al. 2017) which predicts such diffusion rate.The GMM would prove extremely useful in this context.

Figure 1 .
Figure1.The face-on and edge-on images of sample galaxies from NH (Panels a, c, also in Appendix FigureA1as NH ID: 8) and from NH2 (Panels b, d, also in Appendix FigureA2as NH2 ID: 2).Red, green, and blue colors correspond to SDSS i-, r-, and g-band fluxes, respectively.

Figure 2 .
Figure 2. The two-component fit to the profiles of the galaxy shown in Figure 1 (NH ID: 8).The plus symbols show the actual profiles and the continuous lines show the fits.(a) The radial profile of the galaxy (face on) is well reproduced by a combination of two components: an exponential disk and a spheroid with a Sersic index 4.(b) The vertical luminosity profile at (2.0 ± 0.1)R d and |z| ≤ 3 kpc is well reproduced by a two-component model.(c) The vertical mass profile appears slightly different from the mass profile but is still closely reproduced by a two-component model.

Figure 3 .
Figure 3.The measured properties of the thin and thick disks of the NH galaxies in r-band compared with the observational data.Circles denote the NH galaxies.(a) The ratios in the vertical scale height between thick and thin disks.(b) The luminosity fraction of the thick disks.The dashed lines and shades with a matching colors show the median and first-third quantiles of the three data sets.

Figure 4 .
Figure 4.The positions of the spatially-thin and -thick disk stars of all the 27 NH galaxies when they are observed (final snapshot) and at their birth.The abscissa shows the mean radial position of the stars from the galactic centre in the units of the disk scale length in the final snapshot.The ordinate shows the mean position of the stars in the vertical distance from the galaxy disk's mid-plane in units of the thin disk scale height in the final snapshot.The three spatially separated groups of stars in the final snapshot (circles) were barely distinguishable in terms of birth position (star symbols).The thin (|z|/z thick < 0.5; blue), thick (1.5 < |z|/z thick < 2.5; green), and thicker (5 < |z|/z thick < 6; red) component's median positions are shown with the first and third quantiles, where z thick is roughly 0.8 kpc.
-(b) shows a zoom-in map of a region of Panel-(a).The smoothed map clearly shows a bimodal distribution.A simple cut of [α/Fe]=0.06effectively separates the bimodal distribution, and the mass fraction of the "chemically-thick" (0.06 ≤ [α/Fe] < 0.14) component becomes 34% when we measure it at 2R d and |z| ≤ 1 kpc for easy comparison with the observational estimate.In the case of Milky Way, bimodality appears to be separable by a cut of [α/Fe] ≈ 0.2, and the number ratio between the high-[α/Fe] peak and the low-[α/Fe] peak

Figure 5 .
Figure 5.The stellar population properties of the star particles of a sample NH2 galaxy (NH2 ID: 2).(a) The density map of all the star particles in terms of chemical properties.The cyan contours show the Milky Way stellar data adopted from the APOGEE (Jönsson et al. 2020) and GALAH (Buder et al. 2021) projects.The dash-dotted lines show the zero points, and the white box shows the region that is zoomed in Panel-(b).(b) The zoom-in map of the box in Panel-(a).This smoothed distribution appears to be bimodal which is similar to what is found in the Milky Way disk.One may divide the bimodal distribution at the valley marked by the horizontal dashed line at [α/Fe]=0.06,which in turn classifies the stars into thin-(below the line) and thick-disk (above the line) stars.(c) [α/Fe] vs. age for the star particles in Panel-(b).The bimodality in [α/Fe] is clearly associated with the age distribution and hence the star formation history.The coalescence time for the gas-rich minor merger at 6.2 Gyr of lookback time, discussed in the main text, is marked.

Figure 7 .
Figure 7.The kinematic properties of all the star particles of a sample NH galaxy.This GMM exercise identified the five components in order of mass fraction based on three properties (e/|emax|, Jz/Jcirc(e), and Jp/Jcirc(e)).The dashed ellipses show the thin disk (Component 1 marked in blue), thick disk (Components 2 and 4 in green), and spheroid (Components 3 and 5 in red) components.The projected histograms associated with the three groups (blue, green, and red curves) along with the total (black) are also presented.

Figure 8 .
Figure 8.The face-on and edge-on images of the components in Figure 7.The mass fraction of each component is shown in each panel.

Figure 9 .
Figure 9.The chemical and kinematic properties of the star particles of two NH2 galaxies.(a) The same galaxy as shown in Figures 5 and 6.The [α/Fe] bimodality shown in Figure 5 is also visible here.The horizontal dashed line shows the arbitrary cut that separates the low-and high-[α/Fe] stars from the bimodality.The vertical line at ϵ = 0.5 divides rotation-dominant (ϵ ≥ 0.5) and dispersion-dominant (ϵ < 0.5) stars.The other vertical line at ϵ = 0.85 divides kinematically-thin and -thick disk components.(b) The same information but for a different galaxy.This galaxy has more than three main [α/Fe] components in the disk instead of two.

Figure 10 .
Figure 10.The correspondence between different sampling results based on the NH2 galaxy shown in Figure 9-(a) (NH2 ID: 2).(a) A simple circularity-based kinematic sampling scheme.(b) The [α/Fe] distribution of the thin and thick-disk stars sampled by circularity as shown in (a).The two vertical lines show the mean values of the two disk components.As shown in Figure 9, contamination is severe in the case of the thick disk.(c) Same as Panel-(b) but for [Fe/H].(d) Same as Panel-(b) but for age.(e) The vertical distribution of the thin and thick disk stars sampled in (a).(f) The vertical r-band luminosity profile fitted by a two-component model as described in Section 3.2.

Figure 11 .
Figure 11.The key properties of thick disks of NH galaxies compared between the photometric vertical profile fit and the kinematic decomposition at 2R d .(a) The ratio of the vertical scale heights between the thin and thick disks.(b) The r-band luminosity contribution of the thick disks to the total disk luminosity.
Figure 12.The comparison between the circularity at birth and the circularity today for all the stars of one NH galaxy.A one-to-one correlation is shown by the diagonal line.The black curves show the 0.5, 1.0, and 1.5-σ contours.The distribution is heavily biased.Most stars are born with a large value of circularity on a thin rotating disk and are gradually heated to show a much wider distribution in the end.

Figure 13 .
Figure 13.The mean change of circularity of the star particles between their birth and now (∆ϵ = ϵnow − ϵ birth ).Panels (a), (b), and (c) are for three sample galaxies, and Panel (d) shows the mean values for the whole sample of the 27 NH galaxies.The demarcation lines for the thin and thick disk components are marked by vertical dashed lines.The slanted line shows the criterion above which the star in the bin becomes kinematically thin.Considering that most of the stars are born with a high value of circularity, most of the disk stars become kinematically hotter compared to their birth.

Figure A1 .
FigureA1.The face-on and edge-on images of the sixteen most massive disk galaxies in the NH simulation.The blue, green, and red components in the mock images correspond to the Sloan Digital Sky Survey (SDSS) g, r, and i bands.

Figure A3 .
Figure A3.The BIC-selected vertical fitting result, measured at 1 ± 0.1 R d for the same NH galaxies in Figure A1.The data points (gray plus), thin disk (blue line), thick disk (green line), extra component (red line), and total fit (magenta line) are presented.

Figure A7 .
Figure A7.Same as Figure 7 in the main text but for additional NH galaxies for reference.

Figure A9 .
Figure A9.Same as Figure 9 in the main text but for the rest of the NH2 sample.

Figure A10 .
Figure A10.Same as Figure 12 in the main text but for the 16 most massive NH galaxies.