<figure><img src="image_1.png" /></figure>

HAL
open science

# On the Significance of the Thick Disks of Disk Galaxies

Sukyoung K. Yi, J. K. Jang, Julien Devriendt, Yohan Dubois, San Han, Taysun Kimm, Katarina Kraljic, Minjung Park, Sebastien Peirani, Christophe Pichon, et al.

► To cite this version:

Sukyoung K. Yi, J. K. Jang, Julien Devriendt, Yohan Dubois, San Han, et al. On the Significance of the Thick Disks of Disk Galaxies. The Astrophysical Journal Supplement Series, 2024, 271, 10.3847/1538-4365/ad0e71. insu-04490419

HAL Id: insu-04490419

https://insu.hal.science/insu-04490419

Submitted on 10 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the Significance of the Thick Disks of Disk Galaxies

Sukyoung K. Yi<fnref n="1" />, J. K. Jang<fnref n="1" />, Julien Devriendt<fnref n="2" />, Yohan Dubois<fnref n="3" />, San Han<fnref n="1" />, Taysun Kimm<fnref n="1" />, Katarina Kraljic<fnref n="4" />, Minjung Park<fnref n="5" />, Sebastien Peirani<fnref n="3" />,<fnref n="6" />, Christophe Pichon<fnref n="3" />, and Jinsu Rhee<fnref n="1" />,<fnref n="7" />

## Abstract

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_* > 10^{10} M_\odot$ 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 $\text{sech}^2$ components, respectively. The NewHorizon2 data show that the selection of thick-disk stars based on a single [$\alpha$/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.

*Unified Astronomy Thesaurus concepts: Galaxy disks (589); Galaxy stellar disks (1594); Galaxy structure (622); Galaxy chemical evolution (580); Hydrodynamical simulations (767)*

## 1. 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 are their thin disks. 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 the secular diffusion of the 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 MW (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 MW 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 a $\text{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). Dynamical 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 have noted that stars in the solar neighborhood are bimodal in [$\alpha$/Fe] abundance (e.g., Bensby & Feltzing 2006; Reddy et al. 2006; Lee et al. 2011). A naive projection of the [$\alpha$/Fe] bimodality onto kinematic properties leads to a suggestion that the thick-disk fraction is roughly one-third within 1 kpc from the galactic midplane (Lee et al. 2011). This makes the thick disk more significant than previously believed. Does the [$\alpha$/Fe] bimodality indicate a thin-to-thick 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

<figure><img src="image_1.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea; starbrown816@yonsei.ac.kr</footnote>

<footnote n="2">²Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="5">⁵Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, Université Côte d'Azur, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote n="7">⁷Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote>Received 2023 May 31; revised 2023 October 27; accepted 2023 November 16; published 2024 February 12</footnote>
# On the Significance of the Thick Disks of Disk Galaxies

Sukyoung K. Yi<fnref n="1" />, J. K. Jang<fnref n="1" />, Julien Devriendt<fnref n="2" />, Yohan Dubois<fnref n="3" />, San Han<fnref n="1" />, Taysun Kimm<fnref n="1" />, Katarina Kraljic<fnref n="4" />, Minjung Park<fnref n="5" />, Sebastien Peirani<fnref n="3" />,<fnref n="6" />, Christophe Pichon<fnref n="3" />, and Jinsu Rhee<fnref n="1" />,<fnref n="7" />

## Abstract

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_* > 10^{10} M_\odot$ 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 $\text{sech}^2$ components, respectively. The NewHorizon2 data show that the selection of thick-disk stars based on a single [$\alpha$/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.

*Unified Astronomy Thesaurus concepts: Galaxy disks (589); Galaxy stellar disks (1594); Galaxy structure (622); Galaxy chemical evolution (580); Hydrodynamical simulations (767)*

## 1. 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 are their thin disks. 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 the secular diffusion of the 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 MW (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 MW 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 a $\text{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). Dynamical 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 have noted that stars in the solar neighborhood are bimodal in [$\alpha$/Fe] abundance (e.g., Bensby & Feltzing 2006; Reddy et al. 2006; Lee et al. 2011). A naive projection of the [$\alpha$/Fe] bimodality onto kinematic properties leads to a suggestion that the thick-disk fraction is roughly one-third within 1 kpc from the galactic midplane (Lee et al. 2011). This makes the thick disk more significant than previously believed. Does the [$\alpha$/Fe] bimodality indicate a thin-to-thick 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

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, Republic of Korea; starbrown816@yonsei.ac.kr</footnote>

<footnote n="2">²Department of Physics, University of Oxford, Keble Road, Oxford, OX1 3RH, UK</footnote>

<footnote n="3">³Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="5">⁵Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, Université Côte d'Azur, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote n="7">⁷Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote>Received 2023 May 31; revised 2023 October 27; accepted 2023 November 16; published 2024 February 12</footnote>
between 0% and 100%, and the MW 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; VeraCiro 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). The accretion of gas clouds may induce star formation away from the galactic disk midplane and thus occupy 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 MW 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 (NH2) simulation (S. K. Yi et al. 2024, in preparation), which is a lower-resolution twin that started from the same initial conditions, but has a more extensive chemical follow-up (see Section 2). The combination of the high resolution of NH and detailed chemical information of NH2 makes this investigation viable.

## 2. Data

The main data used were from the NH simulation (Dubois et al. 2021). NH is a medium-sized (20 Mpc diameter), cosmological hydrodynamic simulation of galaxy formation that 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 an MW-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_\odot$. Therefore, typical MW-mass galaxies are realized by millions of star particles. The temporal resolution is also exceptionally high: 867 snapshots were stored with a mean time step 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_\odot$ 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, the surface brightness-stellar mass relation, the star formation rate density evolution, the Kennicutt-Schmidt relation, the galaxy mass-size relation, and the 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 is run with a factor of 2 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 of thick disks 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 Type Ia and II supernovae (SNe). 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_\odot$ 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_\odot$) are taken from Woosley & Weaver (1995). The chemical yields of SNe 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 ($\propto 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_\odot$ is expected to have at most 13 SN Ia explosions during its lifetime.

Our galaxy mass threshold ($M \ge 10^{10} M_\odot$) allows at least one million star particles in each galaxy. Of the 48 NH galaxies, we remove three galaxies that are severely disturbed morphologically at $z=0.17$. We finally have 27 disk galaxies of $M_* \ge 10^{10} M_\odot$ with $v/\sigma \ge 1$ for stars in the last snapshot of $z=0.17$, where $v/\sigma$ 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 2 coarser spatial resolution reduced $v/\sigma$. Besides, its enhanced stellar feedback caused the final stellar mass to be smaller in NH2. We present images of the sample galaxies from NH and NH2 in Figure 1 and in Figures A1 and A2 in the Appendix.

## 3. Results

### 3.1. Profile Sampling

With the 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_\odot$ for star particles is also essential for dynamically
between 0% and 100%, and the MW 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; VeraCiro 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). The accretion of gas clouds may induce star formation away from the galactic disk midplane and thus occupy 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 MW 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 (NH2) simulation (S. K. Yi et al. 2024, in preparation), which is a lower-resolution twin that started from the same initial conditions, but has a more extensive chemical follow-up (see Section 2). The combination of the high resolution of NH and detailed chemical information of NH2 makes this investigation viable.

## 2. Data

The main data used were from the NH simulation (Dubois et al. 2021). NH is a medium-sized (20 Mpc diameter), cosmological hydrodynamic simulation of galaxy formation that 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 an MW-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_\odot$. Therefore, typical MW-mass galaxies are realized by millions of star particles. The temporal resolution is also exceptionally high: 867 snapshots were stored with a mean time step 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_\odot$ 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, the surface brightness-stellar mass relation, the star formation rate density evolution, the Kennicutt-Schmidt relation, the galaxy mass-size relation, and the 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 is run with a factor of 2 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 of thick disks 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 Type Ia and II supernovae (SNe). 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_\odot$ 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_\odot$) are taken from Woosley & Weaver (1995). The chemical yields of SNe 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 ($\propto 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_\odot$ is expected to have at most 13 SN Ia explosions during its lifetime.

Our galaxy mass threshold ($M \ge 10^{10} M_\odot$) allows at least one million star particles in each galaxy. Of the 48 NH galaxies, we remove three galaxies that are severely disturbed morphologically at $z=0.17$. We finally have 27 disk galaxies of $M_* \ge 10^{10} M_\odot$ with $v/\sigma \ge 1$ for stars in the last snapshot of $z=0.17$, where $v/\sigma$ 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 2 coarser spatial resolution reduced $v/\sigma$. Besides, its enhanced stellar feedback caused the final stellar mass to be smaller in NH2. We present images of the sample galaxies from NH and NH2 in Figure 1 and in Figures A1 and A2 in the Appendix.

## 3. Results

### 3.1. Profile Sampling

With the 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_\odot$ for star particles is also essential for dynamically
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The face-on and edge-on images of sample galaxies from NH (panels (a) and (c); also in the Appendix in Figure A1 as NH ID: 8) and from NH2 (panels (b) and (d); also in the Appendix in Figure A2 as NH2 ID: 2). The red, green, and blue colors correspond to the SDSS i-, r-, and g-band fluxes, respectively.</figcaption></figure>

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 \pm 0.1 R_d$ and the vertical distance from the midplane $|z| \le 3$ kpc using the sech² 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 three components:

$$ \rho(Z) = \sum_{i=1}^{3} \rho_i \text{sech}^2\left(\frac{|Z|}{2z_i}\right), \qquad (1) $$

where $i = 1, 2, 3$ are for the “thin,” “thick,” and “thicker” disk components and $\rho_i$ and $z_i$ are the density and vertical scale
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The face-on and edge-on images of sample galaxies from NH (panels (a) and (c); also in the Appendix in Figure A1 as NH ID: 8) and from NH2 (panels (b) and (d); also in the Appendix in Figure A2 as NH2 ID: 2). The red, green, and blue colors correspond to the SDSS i-, r-, and g-band fluxes, respectively.</figcaption></figure>

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<sub>d</sub>, using all the star particles inside R<sub>90</sub>, where R<sub>90</sub> stands for the radius within which 90% of the stellar mass resides. We found a linear relation between R<sub>d</sub> and the galaxy’s stellar mass, M*. We use R<sub>d</sub> based on this fit, and we measure the thin- and thick-disk contributions at the radial distance from the galaxy center

2.0 ± 0.1R<sub>d</sub> and the vertical distance from the midplane |z| ≤ 3 kpc using the sech<sup>2</sup> 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 three components:

$$ \rho(Z) = \sum_{i=1}^{3} \rho_i \text{sech}^2\left(\frac{|Z|}{2z_i}\right), \quad (1) $$

where i = 1, 2, 3 are for the “thin,” “thick,” and “thicker” disk components and ρ<sub>i</sub> and z<sub>i</sub> are the density and vertical scale
<figure><img src="image_3.png" /><figcaption>**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 Sérsic index 4. (b) The vertical luminosity profile at $(2.0 \pm 0.1)R_d$ and $|z| \le 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.</figcaption></figure>

height of each component, respectively. We use the Bayesian information criterion (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 the Appendix (Figures A3 and 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 Sérsic profiles (Sérsic 1963). A simple combination of a disk with a Sérsic 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 $\text{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 that 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 three-component 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 compared 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_{\text{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_{\text{circ}}$ has been measured from the stellar particles inside $R_{90}$ in the case of the NH simulation.

The thick-disk luminosity fraction was measured by integrating over the range of $0.5-2.0 R_d$. The mean value of the integrated luminosity fractions of the thick disks of the NH galaxies is $34\% \pm 10\%$, where the error is 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\% \pm 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
<figure><img src="image_5.png" /><figcaption>**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 Sérsic index 4. (b) The vertical luminosity profile at $(2.0 \pm 0.1)R_d$ and $|z| \le 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.</figcaption></figure>

height of each component, respectively. We use the Bayesian information criterion (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 the Appendix (Figures A3 and 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 Sérsic profiles (Sérsic 1963). A simple combination of a disk with a Sérsic 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 $\text{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 that 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 three-component 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 compared 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_{\text{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_{\text{circ}}$ has been measured from the stellar particles inside $R_{90}$ in the case of the NH simulation.

The thick-disk luminosity fraction was measured by integrating over the range of $0.5–2.0 R_d$. The mean value of the integrated luminosity fractions of the thick disks of the NH galaxies is $34\% \pm 10\%$, where the error is 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\% \pm 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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The measured properties of the thin and thick disks of the NH galaxies in the *r* band compared with the observational data. The 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 matching colors show the median and first–third quantiles of the three data sets.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**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 center in units of the disk scale length in the final snapshot. The ordinate shows the mean position of the stars in vertical distance from the galaxy disk’s midplane 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 median positions of the thin ($|z|/z_{\text{thick}} < 0.5$; blue), thick ($1.5 < |z|/z_{\text{thick}} < 2.5$; green), and thicker ($5 < |z|/z_{\text{thick}} < 6$; red) components are shown with the first and third quantiles, where $z_{\text{thick}}$ is roughly 0.8 kpc.</figcaption></figure>

disks, the properties of the thick disks in the NH galaxies 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.

## 3.2. 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 midplane. We fit the vertical *r*-band brightness profiles of the NH galaxies using two sech² components in the cylindrical region at $2 \pm 0.1 R_d$. We consider the stars in the midplane at vertical distance $|z|/z_{\text{thick}} \le 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 midplane “spatially thick” disk stars, based on the same arguments. We select two locations to sample thick-disk stars: $1.5 < |z|/z_{\text{thick}} < 2.5$ and $5 < |z|/z_{\text{thick}} < 6$, where $z_{\text{thick}} \approx 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 and thick-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 midplane within one to two vertical scale heights of the thin disk in the final snapshot ($|z| < z_{\text{thin}}^{z=0.17}$), where the mean value of $z_{\text{thin}}^{z=0.17}$ 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.

## 3.3. Chemical Sampling

As we mentioned in Section 1, the mass fraction of the MW thick disk is estimated to be 33% when $[\alpha/\text{Fe}]$ is used as the sole
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The measured properties of the thin and thick disks of the NH galaxies in the *r* band compared with the observational data. The 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 matching colors show the median and first–third quantiles of the three data sets.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**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 center in units of the disk scale length in the final snapshot. The ordinate shows the mean position of the stars in vertical distance from the galaxy disk’s midplane 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 median positions of the thin ($|z|/z_{\text{thick}} < 0.5$; blue), thick ($1.5 < |z|/z_{\text{thick}} < 2.5$; green), and thicker ($5 < |z|/z_{\text{thick}} < 6$; red) components are shown with the first and third quantiles, where $z_{\text{thick}}$ is roughly 0.8 kpc.</figcaption></figure>

disks, the properties of the thick disks in the NH galaxies 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.

## 3.2. 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 midplane. We fit the vertical *r*-band brightness profiles of the NH galaxies using two sech² components in the cylindrical region at $2 \pm 0.1 R_d$. We consider the stars in the midplane at vertical distance $|z|/z_{\text{thick}} \le 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 midplane “spatially thick” disk stars, based on the same arguments. We select two locations to sample thick-disk stars: $1.5 < |z|/z_{\text{thick}} < 2.5$ and $5 < |z|/z_{\text{thick}} < 6$, where $z_{\text{thick}} \approx 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 and thick-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 midplane within one to two vertical scale heights of the thin disk in the final snapshot ($|z| < z_{\text{thin}}^{z=0.17}$), where the mean value of $z_{\text{thin}}^{z=0.17}$ 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.

## 3.3. Chemical Sampling

As we mentioned in Section 1, the mass fraction of the MW thick disk is estimated to be 33% when $[\alpha/\text{Fe}]$ is used as the sole
<figure><img src="image_7.png" /></figure>

**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 MW stellar data adopted from the APOGEE (Jönsson et al. 2020) and GALAH (Buder et al. 2021) projects. The dashed-dotted lines show the zero-points and the white box shows the region that panel (b) zooms in on. (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 MW 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-disk (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.

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 α elements (e.g., O and Mg). We

<figure><img src="image_8.png" /></figure>

**Figure 6.** The color images of the thin-disk (top row) and thick-disk (bottom row) stars classified by [α/Fe] as illustrated in Figure 5(b) (NH2 ID: 2). Seemingly, the [α/Fe] = 0.06 cut is effective as a thin–thick disk classifier.

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]
versus [Fe/H] plane. Overall, there is a reverse correlation
between these properties. This is naturally expected for a largely
<figure><img src="image_9.png" /><figcaption>**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 MW stellar data adopted from the APOGEE (Jönsson et al. 2020) and GALAH (Buder et al. 2021) projects. The dashed-dotted lines show the zero-points and the white box shows the region that panel (b) zooms in on. (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 MW 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-disk (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.</figcaption></figure>

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 α elements (e.g., O and Mg). We

<figure><img src="image_10.png" /><figcaption>**Figure 6.** The color images of the thin-disk (top row) and thick-disk (bottom row) stars classified by [α/Fe] as illustrated in Figure 5(b) (NH2 ID: 2). Seemingly, the [α/Fe] = 0.06 cut is effective as a thin–thick disk classifier.</figcaption></figure>

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]
versus [Fe/H] plane. Overall, there is a reverse correlation
between these properties. This is naturally expected for a largely
<figure><img src="image_9.png" /><figcaption>**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/|e_{max}|$, $J_z/J_{circ}(e)$, and $J_p/J_{circ}(e)$). The dashed ellipses show the thin-disk (Component 1, marked in blue), thick-disk (Components 2 and 4, marked in green), and spheroid (Components 3 and 5, marked in red) components. The projected histograms associated with the three groups (blue, green, and red curves) along with the total (black) are also presented.</figcaption></figure>

“closed-box” system, considering the chemical yields of SNe 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 MW (e.g., Bensby et al. 2014).

Figure 5(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.06 effectively 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| \le 1$ kpc for easy comparison with the observational estimate. In the case of the MW, bimodality appears to be separable by a cut of [α/Fe] ≈ 0.2, and the number ratio between the high-[$\alpha$/Fe]-peak and the low-[$\alpha$/Fe]-peak stars is 1:2 (Lee et al. 2011), leading to a thick-disk mass fraction of 33%.

The apparent bimodality of the MW 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 MW, the bimodality of the alpha abundance found in the MW disk could be interpreted as a result of a galaxy merger that happened in the past. This

<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

explanation has already been 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
<figure><img src="image_11.png" /><figcaption>**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/|e_{max}|$, $J_z/J_{circ}(e)$, and $J_p/J_{circ}(e)$). The dashed ellipses show the thin-disk (Component 1, marked in blue), thick-disk (Components 2 and 4, marked in green), and spheroid (Components 3 and 5, marked in red) components. The projected histograms associated with the three groups (blue, green, and red curves) along with the total (black) are also presented.</figcaption></figure>

“closed-box” system, considering the chemical yields of SNe 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 MW (e.g., Bensby et al. 2014).

Figure 5(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.06 effectively 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| \le 1$ kpc for easy comparison with the observational estimate. In the case of the MW, bimodality appears to be separable by a cut of [α/Fe] ≈ 0.2, and the number ratio between the high-[$\alpha$/Fe]-peak and the low-[$\alpha$/Fe]-peak stars is 1:2 (Lee et al. 2011), leading to a thick-disk mass fraction of 33%.

The apparent bimodality of the MW 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 MW, the bimodality of the alpha abundance found in the MW disk could be interpreted as a result of a galaxy merger that happened in the past. This

<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

explanation has already been 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
<figure><img src="image_11.png" /><figcaption>**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.75 divides kinematically thin and thick disk components. The percentage of stars in each section is also given. (b) The same information, but for a different galaxy. This galaxy has more than three main [α/Fe] components in the disk, instead of two.</figcaption></figure>

α 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 MW, but the exact value that
divides the bimodality is different: the MW's cut appears to
occur at [α/Fe] = 0.2, while that of the NH2 galaxy occurs 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 (only
panels (b) and (c)) for the NH2 galaxies in the Appendix
(Figures A5 and A6). One can see that the 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 those with [α/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. A 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 Sloan Digital Sky Survey (SDSS)
observations. The proximity between the two values was most
likely achieved by chance.

3.4. 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, 2020;
Obreja et al. 2019; Jang et al. 2023) 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: the total energy *e*, the angular momentum
along the rotational axis *J*<sub>*z*</sub>, and the remaining angular
momentum *J*<sub>*p*</sub>, where *J*<sub>*p*</sub> = *J*<sub>tot</sub> − *J*<sub>*z*</sub>.

Figure 7 shows a sample exercise of the GMM on the galaxy
in Figure 1. The diagram shows the energy *e* normalized by the
maximum energy, and *J*<sub>*z*</sub> normalized by the maximum angular
momentum that a particle can have for a specific binding
energy *e*, *J*<sub>circ</sub>(*e*). The former roughly indicates the mean
distance from the galaxy center; the smaller, the closer. The
latter is also known as circularity (ϵ ≡ *J*<sub>*z*</sub>/*J*<sub>circ</sub>(*e*)), where 1
indicates a perfectly circular rotation with respect to the
rotational axis, −1 denotes a perfectly counterrotating motion,
and 0 represents a dispersion-dominant motion. We also note
that a 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 three to 15 most significant
components by mass contribution, i.e., N<sub>comp</sub> = 3–15. We then
choose the smallest value of N<sub>comp</sub> that satisfies
BIC < BIC(N<sub>comp</sub> = 15) + δ, where δ is tolerance. We use
δ = 0.15, following the approach of Du et al. (2019), and like
some other steps in the GMM, this choice is somewhat
arbitrary. Following these steps, different values of N<sub>comp</sub> are
selected for different galaxies. There is no strong guideline as
to what values should be taken to separate thin and thick disks,
so we adopt 0.75 as a center value between 0.5 and 1.0. A
similar method was exercised by Du et al. (2019). The BIC
analysis found N<sub>comp</sub> = 5 as the optimal choice for the galaxy
in Figure 7. The numerals in the diagram show the mass rank.

Components 1, 2, and 4 have $\bar{\epsilon} > 0.5$, where $\bar{\epsilon}$ denotes the mean value of $\epsilon$ of all the members of each component detected. Therefore, they follow the coherent rotational motion. The other components (3 and 5) are centered around $\bar{\epsilon} \approx 0$, meaning that they do not follow the coherent rotation of the galaxy. By using a cut in the mean circularity of each cluster, one may assign the components to the thin disk ($\bar{\epsilon} > 0.75$) or the thick disk ($0.5 < \bar{\epsilon} \le 0.75$). In this segmentation, the thin disk is therefore defined to be mostly made of quasi-circular in-plane 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
<figure><img src="image_13.png" /><figcaption>**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.75 divides kinematically thin and thick disk components. The percentage of stars in each section is also given. (b) The same information, but for a different galaxy. This galaxy has more than three main [α/Fe] components in the disk, instead of two.</figcaption></figure>

α 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 MW, but the exact value that
divides the bimodality is different: the MW's cut appears to
occur at [α/Fe] = 0.2, while that of the NH2 galaxy occurs 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 (only
panels (b) and (c)) for the NH2 galaxies in the Appendix
(Figures A5 and A6). One can see that the 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 those with [α/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. A 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 Sloan Digital Sky Survey (SDSS)
observations. The proximity between the two values was most
likely achieved by chance.

3.4. 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, 2020;
Obreja et al. 2019; Jang et al. 2023) 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: the total energy *e*, the angular momentum
along the rotational axis *J*<sub>*z*</sub>, and the remaining angular
momentum *J*<sub>*p*</sub>, where *J*<sub>*p*</sub> = *J*<sub>tot</sub> − *J*<sub>*z*</sub>.

Figure 7 shows a sample exercise of the GMM on the galaxy
in Figure 1. The diagram shows the energy *e* normalized by the
maximum energy, and *J*<sub>*z*</sub> normalized by the maximum angular
momentum that a particle can have for a specific binding
energy *e*, *J*<sub>circ</sub>(*e*). The former roughly indicates the mean
distance from the galaxy center; the smaller, the closer. The
latter is also known as circularity (ϵ ≡ *J*<sub>*z*</sub>/*J*<sub>circ</sub>(*e*)), where 1
indicates a perfectly circular rotation with respect to the
rotational axis, −1 denotes a perfectly counterrotating motion,
and 0 represents a dispersion-dominant motion. We also note
that a 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 three to 15 most significant
components by mass contribution, i.e., N<sub>comp</sub> = 3–15. We then
choose the smallest value of N<sub>comp</sub> that satisfies
BIC < BIC(N<sub>comp</sub> = 15) + δ, where δ is tolerance. We use
δ = 0.15, following the approach of Du et al. (2019), and like
some other steps in the GMM, this choice is somewhat
arbitrary. Following these steps, different values of N<sub>comp</sub> are
selected for different galaxies. There is no strong guideline as
to what values should be taken to separate thin and thick disks,
so we adopt 0.75 as a center value between 0.5 and 1.0. A
similar method was exercised by Du et al. (2019). The BIC
analysis found N<sub>comp</sub> = 5 as the optimal choice for the galaxy
in Figure 7. The numerals in the diagram show the mass rank.

Components 1, 2, and 4 have $\bar{\epsilon} > 0.5$, where $\bar{\epsilon}$ denotes the mean value of $\epsilon$ of all the members of each component detected. Therefore, they follow the coherent rotational motion. The other components (3 and 5) are centered around $\bar{\epsilon} \approx 0$, meaning that they do not follow the coherent rotation of the galaxy. By using a cut in the mean circularity of each cluster, one may assign the components to the thin disk ($\bar{\epsilon} > 0.75$) or the thick disk ($0.5 < \bar{\epsilon} \le 0.75$). In this segmentation, the thin disk is therefore defined to be mostly made of quasi-circular in-plane 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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The correspondence between different sampling results based on the NH2 galaxy shown in Figure 9(a) (NH2 ID: 2). (a) The circularity distribution of the three components by GMM. This galaxy has one thin-disk component and one thick-disk component detected by GMM. (b) The [α/Fe] distribution of the thin- and thick-disk stars sampled 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) The same as panel (b), but for [Fe/H]. (d) The 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.</figcaption></figure>

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 the “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 the GMM or in any other method. For reference, we provide the GMM detection maps of additional NH galaxies in the Appendix (Figures A7 and 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 and 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., $\epsilon > 0.75$). 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 < \epsilon \le 0.75$) at all. Panel (a) shows the stars virtually in the whole galaxy at 1 kpc – 2.5R<sub>d</sub> and |z| ≤ 15 kpc. Kinematically thick disk stars (i.e., B = 15.3%) take up only 29.3% 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 kinematically thick disk stars (B+E) is 22.9%. If we measure these values near the solar circle (e.g., at $2 \pm 0.5R_d$ and $|z| \le 1$ kpc), they become 35.0% (chemically thick) and 9.3% (kinematically thick). It should be noted that we used a single circularity cut of 0.75 for thin-/thick-disk classification in Figure 7 for illustrative purposes only. In the case of the galaxy in the upper panel, the GMM detection found that 37.2% of the disk stars were classified as thick-disk stars, compared to 38.2% [(B+E)/(B+C+E+F)] in this figure; thus, there is practically little difference. In this particular galaxy, using the analysis described above, 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 induced by minor mergers 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 that then diffuse away in orbital space. Stellar chemical information is the time capsule of galactic archeology. We provide [α/Fe] versus circularity diagrams for additional NH2 galaxies in the Appendix (Figure A9).
<figure><img src="image_14.png" /><figcaption>**Figure 10.** The correspondence between different sampling results based on the NH2 galaxy shown in Figure 9(a) (NH2 ID: 2). (a) The circularity distribution of the three components by GMM. This galaxy has one thin-disk component and one thick-disk component detected by GMM. (b) The [α/Fe] distribution of the thin- and thick-disk stars sampled 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) The same as panel (b), but for [Fe/H]. (d) The 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.</figcaption></figure>

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 the “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 the GMM or in any other method. For reference, we provide the GMM detection maps of additional NH galaxies in the Appendix (Figures A7 and 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 and 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., $\epsilon > 0.75$). 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 < \epsilon \le 0.75$) at all. Panel (a) shows the stars virtually in the whole galaxy at 1 kpc – 2.5R<sub>d</sub> and |z| ≤ 15 kpc. Kinematically thick disk stars (i.e., B = 15.3%) take up only 29.3% 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 kinematically thick disk stars (B+E) is 22.9%. If we measure these values near the solar circle (e.g., at $2 \pm 0.5R_d$ and $|z| \le 1$ kpc), they become 35.0% (chemically thick) and 9.3% (kinematically thick). It should be noted that we used a single circularity cut of 0.75 for thin-/thick-disk classification in Figure 7 for illustrative purposes only. In the case of the galaxy in the upper panel, the GMM detection found that 37.2% of the disk stars were classified as thick-disk stars, compared to 38.2% [(B+E)/(B+C+E+F)] in this figure; thus, there is practically little difference. In this particular galaxy, using the analysis described above, 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 induced by minor mergers 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 that then diffuse away in orbital space. Stellar chemical information is the time capsule of galactic archeology. We provide [α/Fe] versus circularity diagrams for additional NH2 galaxies in the Appendix (Figure A9).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The key properties of thick disks of NH galaxies compared between the photometric vertical profile fit and the kinematic decomposition at $2 R_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.</figcaption></figure>

We check the correspondence between the results from different sampling methods. Figure 10 shows a comparison for the solar circle region ($2 \pm 0.5 R_d$ and $|z| < 15$ kpc) of a sample NH2 galaxy. Panel (a) shows the circularity distribution of all the stars of the galaxy. The thin-disk, thick-disk, and spheroidal components were classified by the GMM and BIC analysis. The $[\alpha/Fe]$ distribution of the disk components based on this criterion is shown in panel (b). As shown in Figure 9, the low-$[\alpha/Fe]$ peak traces kinematically thin disk stars well, whereas the high-$[\alpha/Fe]$ peak is heavily contaminated by kinematically thin disk stars. Kinematically thick disk stars have various values of $[\alpha/Fe]$ because they are an ensemble of stars born at different times and are thus heated for different lengths of time.

<figure><img src="image_14.png" /><figcaption>**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\sigma$, $1.0\sigma$, and $1.5\sigma$ 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.</figcaption></figure>

Panel (c) shows their metallicity distributions. Panel (d) shows the age distributions. The kinematically thin disk is predominantly young, whereas the kinematically thick disk stars are clearly more biased to older ages. Their vertical distribution is shown in panel (e). In this galaxy, the kinematically thin disk is dominant all the way out to 5 kpc from the midplane; thus, one has to go beyond 5 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 Sections 3.1 and 3.2.

Figure 11 shows the comparison of the key properties of thick disks between the vertical luminosity profile fits and the kinematic detection. Panel (a) shows the vertical scale height ratios measured at $2 R_d$. The profile fits give roughly 50% higher values. Panel (b) shows the thick-disk luminosity fraction. The thick-disk luminosity fraction is $13\% \pm 10\%$ (or $25 \pm 15\%$ in mass). When we integrate over $0.5-2 R_d$, the thick-disk fraction becomes $20\% \pm 11\%$ in r-band luminosity or $35\% \pm 15\%$ 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\sigma$ contour shows that 68% of the stars are born with $\epsilon \gtrsim 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
<figure><img src="image_15.png" /><figcaption>**Figure 11.** The key properties of thick disks of NH galaxies compared between the photometric vertical profile fit and the kinematic decomposition at $2 R_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.</figcaption></figure>

We check the correspondence between the results from different sampling methods. Figure 10 shows a comparison for the solar circle region ($2 \pm 0.5 R_d$ and $|z| < 15$ kpc) of a sample NH2 galaxy. Panel (a) shows the circularity distribution of all the stars of the galaxy. The thin-disk, thick-disk, and spheroidal components were classified by the GMM and BIC analysis. The $[\alpha/Fe]$ distribution of the disk components based on this criterion is shown in panel (b). As shown in Figure 9, the low-$[\alpha/Fe]$ peak traces kinematically thin disk stars well, whereas the high-$[\alpha/Fe]$ peak is heavily contaminated by kinematically thin disk stars. Kinematically thick disk stars have various values of $[\alpha/Fe]$ because they are an ensemble of stars born at different times and are thus heated for different lengths of time.

<figure><img src="image_16.png" /><figcaption>**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\sigma$, $1.0\sigma$, and $1.5\sigma$ 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.</figcaption></figure>

Panel (c) shows their metallicity distributions. Panel (d) shows the age distributions. The kinematically thin disk is predominantly young, whereas the kinematically thick disk stars are clearly more biased to older ages. Their vertical distribution is shown in panel (e). In this galaxy, the kinematically thin disk is dominant all the way out to 5 kpc from the midplane; thus, one has to go beyond 5 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 Sections 3.1 and 3.2.

Figure 11 shows the comparison of the key properties of thick disks between the vertical luminosity profile fits and the kinematic detection. Panel (a) shows the vertical scale height ratios measured at $2 R_d$. The profile fits give roughly 50% higher values. Panel (b) shows the thick-disk luminosity fraction. The thick-disk luminosity fraction is $13\% \pm 10\%$ (or $25\pm 15\%$ in mass). When we integrate over $0.5-2 R_d$, the thick-disk fraction becomes $20\% \pm 11\%$ in r-band luminosity or $35\% \pm 15\%$ 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\sigma$ contour shows that 68% of the stars are born with $\epsilon \gtrsim 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
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The mean change of circularity of the star particles between their birth and now ($\Delta\epsilon = \epsilon_{\text{now}} - \epsilon_{\text{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.</figcaption></figure>

more dispersed at present (i.e., $z = 0.17$). The $1\sigma$ contour extends down to $\epsilon \approx 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, focusing 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 the 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 ($\epsilon < 0.75$). 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 ($\epsilon \ge 0.75$) are all negative, meaning that most of them become kinematically thicker with time. The magnitude of the change is so large that the first four points are below the red reference line, meaning that a large fraction of the thin-born stars becomes thick in the last snapshot. Many of the kinematically hot stars (e.g., $\epsilon \lesssim 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 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.

We summarize 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. The mean value of the thick-disk fraction is 20% ± 11% in r-band luminosity or 35% ± 15% in mass, which is consistent with the observations. Orbital properties are not well traced by the alpha abundance.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** The mean change of circularity of the star particles between their birth and now ($\Delta\epsilon = \epsilon_{\text{now}} - \epsilon_{\text{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.</figcaption></figure>

more dispersed at present (i.e., $z = 0.17$). The $1\sigma$ contour extends down to $\epsilon \approx 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, focusing 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 the 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 ($\epsilon < 0.75$). 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 ($\epsilon \ge 0.75$) are all negative, meaning that most of them become kinematically thicker with time. The magnitude of the change is so large that the first four points are below the red reference line, meaning that a large fraction of the thin-born stars becomes thick in the last snapshot. Many of the kinematically hot stars (e.g., $\epsilon \lesssim 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 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.

We summarize 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. The mean value of the thick-disk fraction is 20% ± 11% in r-band luminosity or 35% ± 15% 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 the 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.

## 4. Conclusions

We give a brief summary first. We used the NH simulation in this investigation. With its high resolutions in space and mass, we resolved the thin disks of massive spiral galaxies, making 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 NH2 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 most massive ($M_* \ge 10^{10} M_\odot$). In addition, we used a stellar $v/\sigma > 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*² 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 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 midplane and spatially thick disk stars away from the midplane. 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 midplane. 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 MW. 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 the [α/Fe] distribution, while there usually was a good correlation between age and [α/Fe], as in Buck (2020). 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 \approx 2$. In such an effectively “closed box” of a chemical evolution testbed, [α/Fe] by and large monotonically decreases with time. Thus, the [α/Fe] bimodality, more specifically the high-[α/Fe] peak, hints at a separate event of star formation from the main star formation history of the MW disk. We found a galaxy in our sample that resembles the MW 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 MW 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 also were significantly contaminated by the stars of other kinematic properties (such as thin-disk stars or dispersion-dominant stars).

We used an unsupervised machine-learning technique, the GMM, 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, this 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 there was a 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 20% ± 11% in *r*-band luminosity or 35% ± 15% in the 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 MW.

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
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 the 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.

## 4. Conclusions

We give a brief summary first. We used the NH simulation in this investigation. With its high resolutions in space and mass, we resolved the thin disks of massive spiral galaxies, making 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 NH2 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 most massive ($M_* \ge 10^{10} M_\odot$). In addition, we used a stellar $v/\sigma > 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*² 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 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 midplane and spatially thick disk stars away from the midplane. 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 midplane. 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 MW. 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 the [α/Fe] distribution, while there usually was a good correlation between age and [α/Fe], as in Buck (2020). 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 \approx 2$. In such an effectively “closed box” of a chemical evolution testbed, [α/Fe] by and large monotonically decreases with time. Thus, the [α/Fe] bimodality, more specifically the high-[α/Fe] peak, hints at a separate event of star formation from the main star formation history of the MW disk. We found a galaxy in our sample that resembles the MW 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 MW 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 also were significantly contaminated by the stars of other kinematic properties (such as thin-disk stars or dispersion-dominant stars).

We used an unsupervised machine-learning technique, the GMM, 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, this 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 there was a 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 20% ± 11% in *r*-band luminosity or 35% ± 15% in the 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 MW.

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 it 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 signatures 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 MW disk. However, this does not necessarily mean that the merger was the dominant contributor to the MW 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 + 1D orbital and chemical space, in the spirit of the bottom panel of Figure 13, so as to predict the mean effect of such a mixture of cosmic accretion events. This could then be compared to inhomogenous kinetic theory (Fouvry et al. 2017), which predicts such a diffusion rate. The GMM would prove extremely useful in this context.

## Acknowledgments

Y.D., S.K.Y., and S.H. ran the NewHorizon and New-Horizon2 simulations. J.R. and T.K. implemented the chemical evolution in the NewHorizon2 simulation. J.K.J. performed most of the hands-on analysis, and S.K.Y. wrote the draft. C.P., J.D., S.P., K.K., and M.P. contributed to the completion of the manuscript. We thank Young Sun Lee for making the SDSS data available to us and for useful discussion. We also thank Seyoung Jeon for his careful proofreading.

This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci; KSC-2017-G2-0003, KSC-2020-CRE-0055, and KSC-2020-CRE-0280 by KISTI; and as a "Grand Challenge" project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. The large data transfer was supported by KREONET, which

is managed and operated by KISTI. S.K.Y. acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769). T.K. was supported in part by the National Research Foundation of Korea (NRF-2020R1C1C1007079). J.R. was supported by the KASI-Yonsei Postdoctoral Fellowship and was supported by the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2023-1-830-00), supervised by the Ministry of Science and ICT. This study was funded in part by the NRF-2022R1A6A1A03053472 grant and the BK21Plus program. Part of the drafting of the paper was carried out during S.K.Y.'s participation in the KITP workshop on Cosmic Web (2023 February), which was supported in part by the National Science Foundation under grant No. NSF PHY-1748958.

## Appendix

We provide supplementary information on the galaxies from the NH and NH2 simulations in this section. They are not critical to the discussion presented in the main text, but are useful as a reference. Table A1 provides the basic properties of the NH galaxies used or presented in this study. Figures A1 and A2 show the face-on and edge-on images of the sampled NH and NH2 galaxies, respectively. Figure A3 shows the vertical profile fits to the NH galaxies measured at the radial distance $1 \pm 0.1 R_d$ from the galactic center, where $R_d$ stands for the scale length of the disk, and the vertical distance $|z| \le 3$ kpc from the galactic midplane. Figure A4 shows the same, but for the measurement at $2 \pm 0.1 R_d$, for comparison. Figure A5 shows the [α/Fe] versus [Fe/H] maps of nine massive NH2 galaxies that appeared in Figure A2. Figure A6 shows the corresponding [α/Fe] versus age maps to the galaxies presented in Figure A5. Figures A7 and A8 show the GMM phase-space diagrams in the same format as Figure 7 in the main text, but for additional NH galaxies. Figure A9 shows [α/Fe] versus circularity at present for the NH2 galaxies, except those already shown in Figure 9 in the main text.
fast. The stars are kinematically heated by force fluctuations along their orbit. Perturbation causes an increase of entropy and thus it 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 signatures 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 MW disk. However, this does not necessarily mean that the merger was the dominant contributor to the MW 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 + 1D orbital and chemical space, in the spirit of the bottom panel of Figure 13, so as to predict the mean effect of such a mixture of cosmic accretion events. This could then be compared to inhomogenous kinetic theory (Fouvry et al. 2017), which predicts such a diffusion rate. The GMM would prove extremely useful in this context.

## Acknowledgments

Y.D., S.K.Y., and S.H. ran the NewHorizon and New-Horizon2 simulations. J.R. and T.K. implemented the chemical evolution in the NewHorizon2 simulation. J.K.J. performed most of the hands-on analysis, and S.K.Y. wrote the draft. C.P., J.D., S.P., K.K., and M.P. contributed to the completion of the manuscript. We thank Young Sun Lee for making the SDSS data available to us and for useful discussion. We also thank Seyoung Jeon for his careful proofreading.

This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci; KSC-2017-G2-0003, KSC-2020-CRE-0055, and KSC-2020-CRE-0280 by KISTI; and as a "Grand Challenge" project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. The large data transfer was supported by KREONET, which

is managed and operated by KISTI. S.K.Y. acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769). T.K. was supported in part by the National Research Foundation of Korea (NRF-2020R1C1C1007079). J.R. was supported by the KASI-Yonsei Postdoctoral Fellowship and was supported by the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2023-1-830-00), supervised by the Ministry of Science and ICT. This study was funded in part by the NRF-2022R1A6A1A03053472 grant and the BK21Plus program. Part of the drafting of the paper was carried out during S.K.Y.'s participation in the KITP workshop on Cosmic Web (2023 February), which was supported in part by the National Science Foundation under grant No. NSF PHY-1748958.

## Appendix

We provide supplementary information on the galaxies from the NH and NH2 simulations in this section. They are not critical to the discussion presented in the main text, but are useful as a reference. Table A1 provides the basic properties of the NH galaxies used or presented in this study. Figures A1 and A2 show the face-on and edge-on images of the sampled NH and NH2 galaxies, respectively. Figure A3 shows the vertical profile fits to the NH galaxies measured at the radial distance $1 \pm 0.1 R_d$ from the galactic center, where $R_d$ stands for the scale length of the disk, and the vertical distance $|z| \le 3$ kpc from the galactic midplane. Figure A4 shows the same, but for the measurement at $2 \pm 0.1 R_d$, for comparison. Figure A5 shows the [α/Fe] versus [Fe/H] maps of nine massive NH2 galaxies that appeared in Figure A2. Figure A6 shows the corresponding [α/Fe] versus age maps to the galaxies presented in Figure A5. Figures A7 and A8 show the GMM phase-space diagrams in the same format as Figure 7 in the main text, but for additional NH galaxies. Figure A9 shows [α/Fe] versus circularity at present for the NH2 galaxies, except those already shown in Figure 9 in the main text.
<table><caption>Table A1
General Properties of the NH Galaxies</caption><thead><tr><th rowspan="2">ID<sup>a</sup></th><th rowspan="2">M<sub>*,3R<sub>e</sub></sub><sup>b</sup></th><th rowspan="2">V<sub>circ</sub><sup>c</sup></th><th rowspan="2">R<sub>e</sub><sup>d</sup></th><th rowspan="2">R<sub>d</sub><sup>e</sup></th><th rowspan="2">R<sub>90</sub><sup>f</sup></th><th colspan="2">z<sub>thin</sub><sup>g</sup></th><th colspan="2">z<sub>thick</sub><sup>h</sup></th><th colspan="2">(L<sub>thick</sub>/L<sub>tot</sub>)<sup>i</sup></th><th>(M<sub>thick</sub>/M<sub>tot</sub>)<sup>j</sup></th></tr><tr><th>1R<sub>d</sub></th><th>2R<sub>d</sub></th><th>1R<sub>d</sub></th><th>2R<sub>d</sub></th><th>phot</th><th>kin</th><th>kin</th></tr></thead><tbody><tr><td>1</td><td>11.21</td><td>304.59</td><td>4.31</td><td>6.83</td><td>14.30</td><td>0.26</td><td>1.41</td><td>0.57</td><td>2.73</td><td>0.18</td><td>0.11</td><td>0.31</td></tr><tr><td>2</td><td>11.16</td><td>343.62</td><td>3.24</td><td>6.77</td><td>12.35</td><td>0.18</td><td>0.81</td><td>0.24</td><td>1.07</td><td>0.28</td><td>0.01</td><td>0.02</td></tr><tr><td>3</td><td>11.08</td><td>283.95</td><td>3.28</td><td>6.54</td><td>13.50</td><td>0.24</td><td>1.16</td><td>0.21</td><td>0.84</td><td>0.25</td><td>0.05</td><td>0.13</td></tr><tr><td>4</td><td>10.89</td><td>258.21</td><td>5.31</td><td>5.98</td><td>14.29</td><td>0.31</td><td>1.19</td><td>0.57</td><td>1.92</td><td>0.25</td><td>0.08</td><td>0.16</td></tr><tr><td>5</td><td>10.79</td><td>214.50</td><td>5.24</td><td>5.70</td><td>15.55</td><td>0.38</td><td>1.07</td><td>0.64</td><td>2.41</td><td>0.26</td><td>0.21</td><td>0.39</td></tr><tr><td>6</td><td>10.78</td><td>235.98</td><td>5.57</td><td>5.68</td><td>15.47</td><td>0.19</td><td>0.99</td><td>0.36</td><td>1.23</td><td>0.34</td><td>0.06</td><td>0.09</td></tr><tr><td>7</td><td>10.77</td><td>226.11</td><td>5.21</td><td>5.65</td><td>15.53</td><td>0.23</td><td>0.89</td><td>0.26</td><td>0.87</td><td>0.30</td><td>0.02</td><td>0.03</td></tr><tr><td>8</td><td>10.57</td><td>194.36</td><td>2.53</td><td>5.06</td><td>13.06</td><td>0.22</td><td>0.90</td><td>0.26</td><td>0.92</td><td>0.20</td><td>0.11</td><td>0.18</td></tr><tr><td>9</td><td>10.50</td><td>169.03</td><td>3.58</td><td>3.43</td><td>10.99</td><td>0.15</td><td>0.72</td><td>0.25</td><td>0.99</td><td>0.52</td><td>0.38</td><td>0.43</td></tr><tr><td>10</td><td>10.49</td><td>181.28</td><td>4.46</td><td>4.83</td><td>11.39</td><td>0.30</td><td>0.84</td><td>0.28</td><td>1.21</td><td>0.07</td><td>0.21</td><td>0.37</td></tr><tr><td>11</td><td>10.42</td><td>160.28</td><td>4.03</td><td>4.19</td><td>10.89</td><td>0.22</td><td>0.79</td><td>0.54</td><td>1.54</td><td>0.26</td><td>0.22</td><td>0.29</td></tr><tr><td>12</td><td>10.40</td><td>159.31</td><td>5.12</td><td>4.40</td><td>12.12</td><td>0.21</td><td>0.75</td><td>0.27</td><td>0.92</td><td>0.38</td><td>0.27</td><td>0.40</td></tr><tr><td>13</td><td>10.37</td><td>169.11</td><td>3.82</td><td>3.36</td><td>12.05</td><td>0.16</td><td>0.88</td><td>0.29</td><td>1.23</td><td>0.41</td><td>0.30</td><td>0.56</td></tr><tr><td>14</td><td>10.36</td><td>155.43</td><td>3.25</td><td>2.75</td><td>9.21</td><td>0.11</td><td>0.71</td><td>0.20</td><td>1.00</td><td>0.48</td><td>0.29</td><td>0.61</td></tr><tr><td>15</td><td>10.34</td><td>168.51</td><td>3.03</td><td>3.00</td><td>13.17</td><td>0.15</td><td>0.81</td><td>0.22</td><td>0.89</td><td>0.22</td><td>0.15</td><td>0.33</td></tr><tr><td>16</td><td>10.28</td><td>148.47</td><td>4.60</td><td>3.66</td><td>11.90</td><td>0.17</td><td>0.80</td><td>0.28</td><td>1.07</td><td>0.48</td><td>0.19</td><td>0.25</td></tr><tr><td>17</td><td>10.24</td><td>125.72</td><td>3.95</td><td>2.31</td><td>9.94</td><td>0.36</td><td>0.86</td><td>0.67</td><td>1.88</td><td>0.21</td><td>0.39</td><td>0.31</td></tr><tr><td>18</td><td>10.19</td><td>125.94</td><td>5.90</td><td>3.99</td><td>12.83</td><td>0.32</td><td>1.06</td><td>0.49</td><td>1.63</td><td>0.32</td><td>0.30</td><td>0.38</td></tr><tr><td>19</td><td>10.19</td><td>127.86</td><td>3.64</td><td>2.51</td><td>9.76</td><td>0.21</td><td>0.77</td><td>0.45</td><td>1.37</td><td>0.30</td><td>1.00</td><td>1.00</td></tr><tr><td>20</td><td>10.16</td><td>152.94</td><td>3.42</td><td>3.35</td><td>8.01</td><td>0.20</td><td>0.89</td><td>0.25</td><td>1.11</td><td>0.30</td><td>0.19</td><td>0.28</td></tr><tr><td>21</td><td>10.14</td><td>129.34</td><td>6.57</td><td>3.85</td><td>12.63</td><td>0.19</td><td>0.73</td><td>0.43</td><td>0.99</td><td>0.44</td><td>0.10</td><td>0.15</td></tr><tr><td>22</td><td>10.08</td><td>132.25</td><td>4.70</td><td>3.65</td><td>10.92</td><td>0.19</td><td>0.68</td><td>0.28</td><td>0.84</td><td>0.38</td><td>0.26</td><td>0.33</td></tr><tr><td>23</td><td>10.07</td><td>107.29</td><td>7.59</td><td>3.63</td><td>14.90</td><td>0.17</td><td>0.79</td><td>0.29</td><td>0.92</td><td>0.46</td><td>0.22</td><td>0.30</td></tr><tr><td>24</td><td>10.06</td><td>110.89</td><td>5.36</td><td>3.60</td><td>12.15</td><td>0.26</td><td>0.83</td><td>0.35</td><td>1.12</td><td>0.43</td><td>1.00</td><td>1.00</td></tr><tr><td>25</td><td>10.04</td><td>127.39</td><td>4.20</td><td>3.52</td><td>10.16</td><td>0.23</td><td>0.72</td><td>0.55</td><td>1.23</td><td>0.51</td><td>0.19</td><td>0.33</td></tr><tr><td>26</td><td>10.01</td><td>131.04</td><td>3.96</td><td>3.47</td><td>10.81</td><td>0.25</td><td>0.85</td><td>0.30</td><td>0.85</td><td>0.28</td><td>0.43</td><td>0.36</td></tr><tr><td>27</td><td>10.00</td><td>122.00</td><td>4.80</td><td>3.44</td><td>13.36</td><td>0.24</td><td>0.99</td><td>0.30</td><td>1.08</td><td>0.25</td><td>0.26</td><td>0.45</td></tr></tbody></table>

Notes.

a The identification number of the NH galaxy.

b The stellar mass inside three effective radii (log-scale) in units of solar mass.

c The circular velocity of the galaxy in kilometers per second.

d The effective radius in kiloparsecs.

e The disk scale radius in kiloparsecs.

f The radius that contains 90% of the total stellar mass in kiloparsecs.

g The photometric vertical scale height of the thin disk, measured at $1R_d$ or $2R_d$ in kiloparsecs.

h The photometric vertical scale height of the thick disk, measured at $1R_d$ or $2R_d$ in kiloparsecs.

i The r-band luminosity fraction of the thick disk from the photometric or the kinematic decomposition, measured at 0.5–2R<sub>d</sub>.

j The mass fraction of the thick disk from the kinematic decomposition, measured at 0.5–2R<sub>d</sub>.
<table><caption>Table A1
General Properties of the NH Galaxies</caption><thead><tr><th rowspan="2">ID<sup>a</sup></th><th rowspan="2">M<sub>*,3R<sub>e</sub></sub><sup>b</sup></th><th rowspan="2">V<sub>circ</sub><sup>c</sup></th><th rowspan="2">R<sub>e</sub><sup>d</sup></th><th rowspan="2">R<sub>d</sub><sup>e</sup></th><th rowspan="2">R<sub>90</sub><sup>f</sup></th><th colspan="2">z<sub>thin</sub><sup>g</sup></th><th colspan="2">z<sub>thick</sub><sup>h</sup></th><th colspan="2">(L<sub>thick</sub>/L<sub>tot</sub>)<sup>i</sup></th><th>(M<sub>thick</sub>/M<sub>tot</sub>)<sup>j</sup></th></tr><tr><th>1R<sub>d</sub></th><th>2R<sub>d</sub></th><th>1R<sub>d</sub></th><th>2R<sub>d</sub></th><th>phot</th><th>kin</th><th>kin</th></tr></thead><tbody><tr><td>1</td><td>11.21</td><td>304.59</td><td>4.31</td><td>6.83</td><td>14.30</td><td>0.26</td><td>1.41</td><td>0.57</td><td>2.73</td><td>0.18</td><td>0.11</td><td>0.31</td></tr><tr><td>2</td><td>11.16</td><td>343.62</td><td>3.24</td><td>6.77</td><td>12.35</td><td>0.18</td><td>0.81</td><td>0.24</td><td>1.07</td><td>0.28</td><td>0.01</td><td>0.02</td></tr><tr><td>3</td><td>11.08</td><td>283.95</td><td>3.28</td><td>6.54</td><td>13.50</td><td>0.24</td><td>1.16</td><td>0.21</td><td>0.84</td><td>0.25</td><td>0.05</td><td>0.13</td></tr><tr><td>4</td><td>10.89</td><td>258.21</td><td>5.31</td><td>5.98</td><td>14.29</td><td>0.31</td><td>1.19</td><td>0.57</td><td>1.92</td><td>0.25</td><td>0.08</td><td>0.16</td></tr><tr><td>5</td><td>10.79</td><td>214.50</td><td>5.24</td><td>5.70</td><td>15.55</td><td>0.38</td><td>1.07</td><td>0.64</td><td>2.41</td><td>0.26</td><td>0.21</td><td>0.39</td></tr><tr><td>6</td><td>10.78</td><td>235.98</td><td>5.57</td><td>5.68</td><td>15.47</td><td>0.19</td><td>0.99</td><td>0.36</td><td>1.23</td><td>0.34</td><td>0.06</td><td>0.09</td></tr><tr><td>7</td><td>10.77</td><td>226.11</td><td>5.21</td><td>5.65</td><td>15.53</td><td>0.23</td><td>0.89</td><td>0.26</td><td>0.87</td><td>0.30</td><td>0.02</td><td>0.03</td></tr><tr><td>8</td><td>10.57</td><td>194.36</td><td>2.53</td><td>5.06</td><td>13.06</td><td>0.22</td><td>0.90</td><td>0.26</td><td>0.92</td><td>0.20</td><td>0.11</td><td>0.18</td></tr><tr><td>9</td><td>10.50</td><td>169.03</td><td>3.58</td><td>3.43</td><td>10.99</td><td>0.15</td><td>0.72</td><td>0.25</td><td>0.99</td><td>0.52</td><td>0.38</td><td>0.43</td></tr><tr><td>10</td><td>10.49</td><td>181.28</td><td>4.46</td><td>4.83</td><td>11.39</td><td>0.30</td><td>0.84</td><td>0.28</td><td>1.21</td><td>0.07</td><td>0.21</td><td>0.37</td></tr><tr><td>11</td><td>10.42</td><td>160.28</td><td>4.03</td><td>4.19</td><td>10.89</td><td>0.22</td><td>0.79</td><td>0.54</td><td>1.54</td><td>0.26</td><td>0.22</td><td>0.29</td></tr><tr><td>12</td><td>10.40</td><td>159.31</td><td>5.12</td><td>4.40</td><td>12.12</td><td>0.21</td><td>0.75</td><td>0.27</td><td>0.92</td><td>0.38</td><td>0.27</td><td>0.40</td></tr><tr><td>13</td><td>10.37</td><td>169.11</td><td>3.82</td><td>3.36</td><td>12.05</td><td>0.16</td><td>0.88</td><td>0.29</td><td>1.23</td><td>0.41</td><td>0.30</td><td>0.56</td></tr><tr><td>14</td><td>10.36</td><td>155.43</td><td>3.25</td><td>2.75</td><td>9.21</td><td>0.11</td><td>0.71</td><td>0.20</td><td>1.00</td><td>0.48</td><td>0.29</td><td>0.61</td></tr><tr><td>15</td><td>10.34</td><td>168.51</td><td>3.03</td><td>3.00</td><td>13.17</td><td>0.15</td><td>0.81</td><td>0.22</td><td>0.89</td><td>0.22</td><td>0.15</td><td>0.33</td></tr><tr><td>16</td><td>10.28</td><td>148.47</td><td>4.60</td><td>3.66</td><td>11.90</td><td>0.17</td><td>0.80</td><td>0.28</td><td>1.07</td><td>0.48</td><td>0.19</td><td>0.25</td></tr><tr><td>17</td><td>10.24</td><td>125.72</td><td>3.95</td><td>2.31</td><td>9.94</td><td>0.36</td><td>0.86</td><td>0.67</td><td>1.88</td><td>0.21</td><td>0.39</td><td>0.31</td></tr><tr><td>18</td><td>10.19</td><td>125.94</td><td>5.90</td><td>3.99</td><td>12.83</td><td>0.32</td><td>1.06</td><td>0.49</td><td>1.63</td><td>0.32</td><td>0.30</td><td>0.38</td></tr><tr><td>19</td><td>10.19</td><td>127.86</td><td>3.64</td><td>2.51</td><td>9.76</td><td>0.21</td><td>0.77</td><td>0.45</td><td>1.37</td><td>0.30</td><td>1.00</td><td>1.00</td></tr><tr><td>20</td><td>10.16</td><td>152.94</td><td>3.42</td><td>3.35</td><td>8.01</td><td>0.20</td><td>0.89</td><td>0.25</td><td>1.11</td><td>0.30</td><td>0.19</td><td>0.28</td></tr><tr><td>21</td><td>10.14</td><td>129.34</td><td>6.57</td><td>3.85</td><td>12.63</td><td>0.19</td><td>0.73</td><td>0.43</td><td>0.99</td><td>0.44</td><td>0.10</td><td>0.15</td></tr><tr><td>22</td><td>10.08</td><td>132.25</td><td>4.70</td><td>3.65</td><td>10.92</td><td>0.19</td><td>0.68</td><td>0.28</td><td>0.84</td><td>0.38</td><td>0.26</td><td>0.33</td></tr><tr><td>23</td><td>10.07</td><td>107.29</td><td>7.59</td><td>3.63</td><td>14.90</td><td>0.17</td><td>0.79</td><td>0.29</td><td>0.92</td><td>0.46</td><td>0.22</td><td>0.30</td></tr><tr><td>24</td><td>10.06</td><td>110.89</td><td>5.36</td><td>3.60</td><td>12.15</td><td>0.26</td><td>0.83</td><td>0.35</td><td>1.12</td><td>0.43</td><td>1.00</td><td>1.00</td></tr><tr><td>25</td><td>10.04</td><td>127.39</td><td>4.20</td><td>3.52</td><td>10.16</td><td>0.23</td><td>0.72</td><td>0.55</td><td>1.23</td><td>0.51</td><td>0.19</td><td>0.33</td></tr><tr><td>26</td><td>10.01</td><td>131.04</td><td>3.96</td><td>3.47</td><td>10.81</td><td>0.25</td><td>0.85</td><td>0.30</td><td>0.85</td><td>0.28</td><td>0.43</td><td>0.36</td></tr><tr><td>27</td><td>10.00</td><td>122.00</td><td>4.80</td><td>3.44</td><td>13.36</td><td>0.24</td><td>0.99</td><td>0.30</td><td>1.08</td><td>0.25</td><td>0.26</td><td>0.45</td></tr></tbody></table>

Notes.

a The identification number of the NH galaxy.

b The stellar mass inside three effective radii (log-scale) in units of solar mass.

c The circular velocity of the galaxy in kilometers per second.

d The effective radius in kiloparsecs.

e The disk scale radius in kiloparsecs.

f The radius that contains 90% of the total stellar mass in kiloparsecs.

g The photometric vertical scale height of the thin disk, measured at $1R_d$ or $2R_d$ in kiloparsecs.

h The photometric vertical scale height of the thick disk, measured at $1R_d$ or $2R_d$ in kiloparsecs.

i The r-band luminosity fraction of the thick disk from the photometric or the kinematic decomposition, measured at 0.5–2R<sub>d</sub>.

j The mass fraction of the thick disk from the kinematic decomposition, measured at 0.5–2R<sub>d</sub>.
<figure><img src="image_16.png" /><figcaption>**Figure A1.** The face-on and edge-on images of the 16 most massive disk galaxies in the NH simulation. The blue, green, and red components in the mock images correspond to the SDSS *g*, *r*, and *i* bands.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A1.** The face-on and edge-on images of the 16 most massive disk galaxies in the NH simulation. The blue, green, and red components in the mock images correspond to the SDSS *g*, *r*, and *i* bands.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure A2.** The same as Figure A1, but for the sample of 10 NH2 galaxies.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure A2.** The same as Figure A1, but for the sample of 10 NH2 galaxies.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A3.** The BIC-selected vertical fitting result, measured at $1 \pm 0.1 R_d$ for the same NH galaxies in Figure A1. The data points (gray plus symbols), thin disks (blue lines), thick disks (green lines), extra components (red lines), and total fits (magenta lines) are presented.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure A3.** The BIC-selected vertical fitting result, measured at $1 \pm 0.1 R_d$ for the same NH galaxies in Figure A1. The data points (gray plus symbols), thin disks (blue lines), thick disks (green lines), extra components (red lines), and total fits (magenta lines) are presented.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure A4.** The same as Figure A3, but for the measurement performed at $2 \pm 0.1 R_d$ for reference.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure A4.** The same as Figure A3, but for the measurement performed at $2 \pm 0.1 R_d$ for reference.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure A5. Similar to Figure 5(b) in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure A5. Similar to Figure 5(b) in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure A6. Similar to Figure 5(c) in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure A6. Similar to Figure 5(c) in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure A7.** Same as Figure 7 in the main text, but for additional NH galaxies for reference.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure A7.** Same as Figure 7 in the main text, but for additional NH galaxies for reference.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure A8. The same as Figure A7, but for different galaxies.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure A8. The same as Figure A7, but for different galaxies.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure A9. The same as Figure 9 in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure A9. The same as Figure 9 in the main text, but for the rest of the NH2 sample.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure A10. The same as Figure 12 in the main text, but for the 16 most massive NH galaxies.</figcaption></figure>

## ORCID iDs

Sukyoung K. Yi https://orcid.org/0000-0002-4556-2619

J. K. Jang https://orcid.org/0000-0002-0858-5264

Julien Devriendt https://orcid.org/0000-0002-8140-0422

Yohan Dubois https://orcid.org/0000-0003-0225-6387

San Han https://orcid.org/0000-0001-9939-713X

Taysun Kimm https://orcid.org/0000-0002-3950-3997

Katarina Kraljic https://orcid.org/0000-0001-6180-0245

Minjung Park https://orcid.org/0000-0002-8435-9402

Christophe Pichon https://orcid.org/0000-0003-0695-6735

Jinsu Rhee https://orcid.org/0000-0002-0184-9589

## References

Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21

Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Bensby, T., & Feltzing, S. 2006, MNRAS, 367, 1181

Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71

Binney, J., & Lacey, C. 1988, MNRAS, 230, 597

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bovy, J., & Rix, H. W. 2013, ApJ, 779, 115

Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894

Buck, T. 2020, MNRAS, 491, 5435

Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150

Burstein, D. 1979, ApJ, 234, 829

Chabrier, G. 2005, in The Initial Mass Function 50 Years Later, ed. E. Corbelli, F. Palla, & H. Zinnecker (Dordrecht: Springer), 41

Comerón, S., Elmegreen, B. G., Knapen, J. H., et al. 2011, ApJ, 741, 28

Comerón, S., Salo, H., & Knapen, J. H. 2018, A&A, 610, A5

Du, M., Ho, L. C., Debattista, V. P., et al. 2020, ApJ, 895, 139

Du, M., Ho, L. C., Zhao, D., et al. 2019, ApJ, 884, 129

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Fouvry, J.-B., Pichon, C., Chavanis, P.-H., & Monk, L. 2017, MNRAS, 471, 2642

Gilmore, G., & Reid, N. 1983, MNRAS, 202, 1025

Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 616, A86

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439

Jang, J. K., Yi, S. K., Dubois, Y., et al. 2023, ApJ, 950, 4

Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120

Juric, M., Ivezic, Z., Brooks, A., et al. 2008, ApJ, 673, 864

Kazantzidis, S., Bullock, J. S. R., R. A., Z., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254

Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145

Lee, Y. S., Beers, T. C., An, D., et al. 2011, ApJ, 738, 187
<figure><img src="image_27.png" /><figcaption>Figure A10. The same as Figure 12 in the main text, but for the 16 most massive NH galaxies.</figcaption></figure>

## ORCID iDs

Sukyoung K. Yi https://orcid.org/0000-0002-4556-2619

J. K. Jang https://orcid.org/0000-0002-0858-5264

Julien Devriendt https://orcid.org/0000-0002-8140-0422

Yohan Dubois https://orcid.org/0000-0003-0225-6387

San Han https://orcid.org/0000-0001-9939-713X

Taysun Kimm https://orcid.org/0000-0002-3950-3997

Katarina Kraljic https://orcid.org/0000-0001-6180-0245

Minjung Park https://orcid.org/0000-0002-8435-9402

Christophe Pichon https://orcid.org/0000-0003-0695-6735

Jinsu Rhee https://orcid.org/0000-0002-0184-9589

## References

Abadi, M. G., Navarro, J. F., Steinmetz, M., & Eke, V. R. 2003, ApJ, 597, 21

Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Bensby, T., & Feltzing, S. 2006, MNRAS, 367, 1181

Bensby, T., Feltzing, S., & Oey, M. S. 2014, A&A, 562, A71

Binney, J., & Lacey, C. 1988, MNRAS, 230, 597

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bovy, J., & Rix, H. W. 2013, ApJ, 779, 115

Brook, C. B., Kawata, D., Gibson, B. K., & Freeman, K. C. 2004, ApJ, 612, 894

Buck, T. 2020, MNRAS, 491, 5435

Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150

Burstein, D. 1979, ApJ, 234, 829

Chabrier, G. 2005, in The Initial Mass Function 50 Years Later, ed. E. Corbelli, F. Palla, & H. Zinnecker (Dordrecht: Springer), 41

Comerón, S., Elmegreen, B. G., Knapen, J. H., et al. 2011, ApJ, 741, 28

Comerón, S., Salo, H., & Knapen, J. H. 2018, A&A, 610, A5

Du, M., Ho, L. C., Debattista, V. P., et al. 2020, ApJ, 895, 139

Du, M., Ho, L. C., Zhao, D., et al. 2019, ApJ, 884, 129

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Fouvry, J.-B., Pichon, C., Chavanis, P.-H., & Monk, L. 2017, MNRAS, 471, 2642

Gilmore, G., & Reid, N. 1983, MNRAS, 202, 1025

Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2018, A&A, 616, A86

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439

Jang, J. K., Yi, S. K., Dubois, Y., et al. 2023, ApJ, 950, 4

Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120

Juric, M., Ivezic, Z., Brooks, A., et al. 2008, ApJ, 673, 864

Kazantzidis, S., Bullock, J. S. R., R. A., Z., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254

Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145

Lee, Y. S., Beers, T. C., An, D., et al. 2011, ApJ, 738, 187
Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Maeder, A., & Meynet, G. 2000, A&A, 361, 159

Maoz, D., Mannucci, F., & Brandt, T. D. 2012, MNRAS, 426, 3282

Martínez-Lombilla, C., & Knapen, J. H. 2019, A&A, 629, A12

Norris, J. E., & Ryan, S. G. 1991, ApJ, 380, 403

Obreja, A., Dutton, A. A., Macciò, A. V., et al. 2019, MNRAS, 487, 4424

Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2

Peebles, P. J. E. 1969, ApJ, 155, 393

Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74

Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1

Reddy, B. E., Lambert, D. L., & Allende Prieto, C. 2006, MNRAS, 367, 1329

Roškar, R., Debattista, V. P., & Loebman, S. R. 2013, MNRAS, 433, 976

Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269

Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61

Sellwood, J. A., Nelson, R. W., & Tremaine, S. 1998, ApJ, 506, 590

Sérsic, J. L. 1963, BAAA, 6, 41

Spitzer, L., & Schwarzschild, M. 1951, ApJ, 114, 385

Teyssier, R. 2002, A&A, 385, 337

Tinsley, B. M. 1980, FCPh, 5, 287

van der Kruit, P. C. 1981, A&A, 99, 298

VeraCiro, C., D'Onghia, E., Navarro, J., & Abadi, M. 2014, ApJ, 794, 173

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Yoachim, P., & Dalcanton, J. J. 2006, AJ, 131, 226

Yoachim, P., & Dalcanton, J. J. 2008, ApJ, 682, 1004
Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Maeder, A., & Meynet, G. 2000, A&A, 361, 159

Maoz, D., Mannucci, F., & Brandt, T. D. 2012, MNRAS, 426, 3282

Martínez-Lombilla, C., & Knapen, J. H. 2019, A&A, 629, A12

Norris, J. E., & Ryan, S. G. 1991, ApJ, 380, 403

Obreja, A., Dutton, A. A., Macciò, A. V., et al. 2019, MNRAS, 487, 4424

Park, M. J., Yi, S. K., Peirani, S., et al. 2021, ApJS, 254, 2

Peebles, P. J. E. 1969, ApJ, 155, 393

Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493

Quinn, P. J., Hernquist, L., & Fullagar, D. P. 1993, ApJ, 403, 74

Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1

Reddy, B. E., Lambert, D. L., & Allende Prieto, C. 2006, MNRAS, 367, 1329

Roškar, R., Debattista, V. P., & Loebman, S. R. 2013, MNRAS, 433, 976

Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269

Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61

Sellwood, J. A., Nelson, R. W., & Tremaine, S. 1998, ApJ, 506, 590

Sérsic, J. L. 1963, BAAA, 6, 41

Spitzer, L., & Schwarzschild, M. 1951, ApJ, 114, 385

Teyssier, R. 2002, A&A, 385, 337

Tinsley, B. M. 1980, FCPh, 5, 287

van der Kruit, P. C. 1981, A&A, 99, 298

VeraCiro, C., D'Onghia, E., Navarro, J., & Abadi, M. 2014, ApJ, 794, 173

Woosley, S. E., & Weaver, T. A. 1995, ApJS, 101, 181

Yoachim, P., & Dalcanton, J. J. 2006, AJ, 131, 226

Yoachim, P., & Dalcanton, J. J. 2008, ApJ, 682, 1004