# From giant clumps to clouds - II. The emergence of thick disc kinematics
from the conditions of star formation in high redshift gas rich galaxies

Floor van Donkelaar et al.
1Lund Observatory, Department of Astronomy and Theoretical Physics, Lund
University, Box 43, SE-221 00 Lund, Sweden
Floor van Donkelaar1,2, Oscar Agertz1 and Florent Renaud1
1 Lund Observatory, Department of Astronomy and Theoretical Physics, Box 43,
SE-221 00 Lund, Sweden
2Center for Theoretical Astrophysics and Cosmology, Institute for
Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057
Zürich, Switzerland<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

High redshift disc galaxies are more gas rich, clumpier, and more turbulent
than local Universe galaxies. This early era of galaxy formation imprints the
distribution and kinematics of the stars that we observe today, but it is not
yet well established how. In this work, we use simulations of isolated Milky
Way-mass disc galaxies to study how kinematic properties of stars change when
varying the gas fraction. This allows us to quantify the roles played by
internal processes, e.g. gas turbulence and gravitational scattering off
massive gas clumps, in establishing the observed stellar velocity dispersions
and orbital eccentricities. We find that models with gas fractions $>20$ per
cent feature a turbulent and clumpy interstellar medium (ISM), leading to
zero-age stellar velocity dispersions $\sim 20-30~{}{\rm km\,s}^{-1}$ and high
mean orbital eccentricities. Low eccentricities cannot arise from these
physical conditions. For gas fractions below $20$ per cent, the ISM becomes
less turbulent, with stellar velocity dispersions $<10~{}{\rm km\,s}^{-1}$,
and nearly circular orbits for young stars. The turbulence present in gas-rich
high redshift galaxies hence acts as a ‘barrier’ against the formation of thin
discs. We compare our findings to the Milky Way’s age-velocity dispersion
relation and argue that velocity dispersions imprinted already at star
formation by the ISM contribute significantly at all times. Finally, we show
that observed orbital eccentricities in the Milky Way’s thick and thin discs
can be explained entirely as imprints by the star-forming ISM, rather than by
mergers or secular processes.

###### keywords:

galaxies: evolution – galaxies: structure – Galaxy: disc – methods: numerical

††pubyear: 2021††pagerange: From giant clumps to clouds - II. The emergence of
thick disc kinematics from the conditions of star formation in high redshift
gas rich galaxies–From giant clumps to clouds - II. The emergence of thick
disc kinematics from the conditions of star formation in high redshift gas
rich galaxies

##  1 Introduction

The Milky Way’s past evolution can only be studied through galactic
archaeology (Bland-Hawthorn & Gerhard, 2016). Using astrometric, spectroscopic
and IFU surveys, the detailed kinematics and evolution of the Milky Way, and
local disc galaxies in general, can now be studied in unprecedented detail
(e.g. Gaia, APOGEE, 4MOST, MaNGA, SAMI, Gaia Collaboration et al., 2016;
Majewski et al., 2017; de Jong et al., 2019; Bundy et al., 2015; Bryant et
al., 2015). At higher redshift, large galaxy surveys provide insight on the
kinematic state of gas and stars in young, intensively star-forming galaxies
(see review by Glazebrook, 2013). Yet, how high redshift galaxies evolve into
the well ordered spirals we observe in the local Universe is a central problem
in astrophysics (e.g. Renaud et al., 2021b).

The Milky Way has a thin and thick disc component, with the thick disc thought
to carry signatures of its high redshift formation. It is not only
geometrically thicker (Reddy et al., 2003) and more radially compact than the
thin disc, but also richer in $\alpha$-elements (Bensby et al., 2011; Bensby
et al., 2014; Hayden et al., 2015), and contains older stars ($\sim 8-12$ Gyr,
Bensby et al., 2014; Feuillet et al., 2019). Stellar velocity dispersions are
observed to monotonically increase with the age of stellar populations in the
Milky Way, with the thick disc observed to have velocity dispersions $\gtrsim
30{\,\rm{km\,s^{-1}}}$ (Mackereth et al., 2019), compared to the $\sim
10-20{\,\rm{km\,s^{-1}}}$ for the thin disc (Nordström et al., 2004; Hayden et
al., 2020). This age-velocity dispersion relation (AVR) can arise from
dynamical heating of stars, at least to ages of $\sim 8$ Gyr (Gustafsson et
al., 2016; Yu & Liu, 2018; Mackereth et al., 2019; Ting & Rix, 2019). However,
constant orbital heating is thought to be insufficient to explain the
relatively high velocity dispersions of the thick disc (e.g. Sellwood, 2014;
Mackereth et al., 2019; Aumer & Binney, 2017). As such, a number of mechanisms
for creating thick discs have been proposed in the literature, including
vertical disc heating by satellite encounters (Quinn et al., 1993; Kazantzidis
et al., 2009), accretion of satellite stars (Abadi et al., 2003; Read et al.,
2008), star formation in turbulent gas rich discs (Bournaud & Elmegreen, 2009;
Martig et al., 2014), bursty star formation at high redshift (Yu et al.,
2021), and secular formation by radial migration of kinematically hot stars
from the inner to the outer disc (Sellwood & Binney, 2002; Loebman et al.,
2011).

Sales et al. (2009) demonstrated, using numerical simulations, how several of
the above listed thick disc formation scenarios result in different
distributions of stellar orbital eccentricities. As such, orbital
eccentricities provide complementary information to the AVR, and can be used
to constrain the scenarios. Li et al. (2018) used Gaia and APOGEE data to show
that the thick disc is characterised by a higher fraction of stars on
eccentric orbits in comparison to the thin disc. Beraldo e Silva et al. (2021)
combined this data with stellar ages and found an overlap in terms of
eccentricities for old ( > 10 Gyr) stars in the thin and thick discs. This is
suggests a co-formation scenario of the discs, in contrast to purely
sequential formation.

The degree to which the observed thick disc properties is an imprint of the
physical conditions of the star-forming ISM they are born, and hence
originates from a purely in-situ process, is a debated topic. Observations of
disc galaxies out to $z\sim 3$ reveals that the velocity dispersion of the ISM
increases with redshift (e.g. Gnerucci et al., 2011; Kassin et al., 2012;
Wisnioski et al., 2015; Mieda et al., 2016), and thick disc stellar kinematics
can, in principle, be accounted for if the ISM kinematics were directly mapped
onto the stars.

In this work we aim to pinpoint the role played by the turbulent and clumpy
ISM in gas rich high redshift galaxies in establishing thick disc kinematics.
While, in principle, cosmological simulations of galaxy formation are adequate
tools to study this scenario (see e.g. Brook et al., 2004; Martig et al.,
2014; Agertz et al., 2021; Renaud et al., 2021b; Yu et al., 2021), such
simulations suffer from the lack of control, as multiple physical processes
are simultaneously active. We overcome this problem by carrying out several
non-cosmological simulations of Milky Way-mass disc galaxies. The key
parameter here is the gas fraction of the galaxies, which affects the nature
of gravitational instability and turbulent state of the disc (Romeo, 1992;
Romeo & Agertz, 2014; Renaud et al., 2021a). This approach allows us to self-
consistently model the development of a turbulent ISM from disc instabilities
and stellar feedback, while excluding effects such as mergers and gas
accretion.

This paper is organized as follows. In Section 2 we describe the numerical
method, including our choice of galaxy formation physics and the simulation
setup. In Section 3 we present the simulations results, with a focus on the
dependence of stellar kinematics on the gas fraction and how our results
connect to the observed properties of the Milky Way’s thin and thick discs.
Finally we discuss our results in Sections 4 and conclude in 5.

Table 1: Properties of the different simulations, all measured at 300 Myr unless otherwise indicated. Simulation names are indicated in column 1. Column 2: gas fraction at 300 Myr; Column 3: star formation rates with the error margins representing the spread in values between 250 and 300 Myr; Column 4: mean velocity dispersion, $\pm 1$ standard deviation, of stars younger than 10 Myr; Column 5: approximate redshifts in comparison to Milky Way-like galaxies found by van Dokkum et al. (2013); Column 6: approximate redshifts to allow for birth dispersion to explain the Milky Way’s age-velocity dispersion relation (see Section 3.2.1); Column 7: mean eccentricities of stars, $\pm 1$ standard deviation.  Name |  $f_{\rm gas}$ (at 300 Myr) | SFR | $\sigma_{\rm z,birth}$ | Redshift | Redshift | Eccentricity
---|---|---|---|---|---|---
| per cent | [M⊙ yr-1] | [km s-1] | (SFR matching) | (AVR matching) |
fg10 | 9 | 1 $\pm$ 0.25 | 5.2 $\pm$ 0.4 | 0.0 - 0.2 | 0.0 - 0.2 | 0.10 $\pm$ 0.06
fg15 | 16 | 5 $\pm$ 0.8 | 8.6 $\pm$ 1.8 | 0.5 - 1.7 | 0.1 - 0.4 | 0.33 $\pm$ 0.11
fg20 | 21 | 15 $\pm$ 2 | 17.9 $\pm$ 3.5 | 1.4 - 2.6 | 0.4 - 0.8 | 0.49 $\pm$ 0.16
fg30 | 29 | 25 $\pm$ 5 | 27.3 $\pm$ 3.8 | 1.8 - 2.6 | 0.8 - 2 | 0.53 $\pm$ 0.21
fg45 | 46 | 50 $\pm$ 10 | 25.8 $\pm$ 4.3 | 1.8 - 2.6 | 0.8 - 2 | 0.47 $\pm$ 0.22
Figure 1: Stellar (upper) and gas (under) density projections of the simulated
galaxies at 300 Myr. The initial gas fraction increases, from left to right,
from 10 to 70 per cent. The gas density projections show the presence of
flocculent spiral structure in all galaxies. In the higher gas fraction
regimes, there is a clear presence of massive clumps and more irregular ISM.

##  2 Method

We use the hydrodynamic+$N$-body code RAMSES (Teyssier, 2002) to model Milky
Way-mass galaxies in isolation. The fluid dynamics of the baryons is
calculated using a second-order unsplit Godunov method, with the equation of
state of an ideal mono-atomic gas with an adiabatic index $\gamma=5/3$. The
code accounts for metallicity dependent cooling by using the tabulated cooling
functions of Sutherland & Dopita (1993) for gas temperatures $>10^{4}$ K, and
rates from Rosen & Bregman (1995) for cooling down to lower temperatures. The
collisionless dynamics of stellar and dark matter particles is evolved using
the particle-mesh technique (Hockney & Eastwood, 1981), with gravitational
accelerations computed from the gravitational potential on the mesh using a
multi-grid method (Guillet & Teyssier, 2011).

We achieve high resolution in high density regions using adaptive mesh
refinement. A cell is split if its baryonic mass (gas and stars) exceeds 4014
${\,\rm M}_{\odot}$. In addition, a cell is allowed to refine if it contains
more than 8 dark matter particles. This allows the local force softening to
closely match the local mean inter-particle separation, which suppresses
discreteness effects (e.g. Romeo et al., 2008). The maximum refinement level
is set to allow for a mean physical resolution of $\sim 9{\,\rm pc}$ in dense
gas.

The adopted star formation and feedback physics is presented in Agertz et al.
(2013, see also ) Briefly, star formation is treated as a Poisson process,
sampled using star particles of $1000{\,\rm M}_{\odot}$, occurring on a cell-
by-cell basis according to the star formation law:

$\dot{\rho}_{\star}=\epsilon_{\rm ff}\frac{\rho_{\rm g}}{t_{\rm
ff}}\quad{\mbox{for}}\quad\rho>\rho_{\rm SF}.$ (1)

Here $\dot{\rho}_{\star}$ is the star formation rate (SFR) density, $\rho_{\rm
g}$ the gas density, $t_{\rm ff}=\sqrt{3\pi/32G\rho_{\rm g}}$ is the local
free-fall time of the gas, $\rho_{\rm SF}=100~{}{\rm cm}^{-3}$ is the star
formation threshold, and $\epsilon_{\rm ff}$ is the local star formation
efficiency per free-fall time of gas in the cell. Observationally,
$\epsilon_{\rm ff}$ averages at about $1$ per cent in Milky Way giant
molecular clouds (GMCs Krumholz & Tan 2007), albeit with a spread of several
dex (Murray, 2011; Lee et al., 2016). Grisdale et al. (2019) demonstrated how
high efficiencies ($\epsilon_{\rm ff}\sim 10$ per cent) on scales of parsecs
provides a close match to observed efficiencies on scales of individual GMCs,
as well as for the properties of the cold ISM (Grisdale et al., 2017, 2018).
Motivated by these findings, we adopt $\epsilon_{\rm ff}=10$ per cent. Each
formed star particle is treated as a single-age stellar population with a
Chabrier (2003) initial mass function. We account for injection of energy,
momentum, mass, and heavy elements over time from core-collapse supernovae
(SNe) and SNIa, stellar winds, and radiation pressure on the surrounding gas.
Each mechanism depends on the stellar age, the mass and stellar metallicity
(through the metallicity dependent age-mass relation of Raiteri et al. 1996,
calibrated on the stellar evolution code STARBURST99, Leitherer et al. 1999).

SNe explosions are modelled as discrete events, and we follow Kim & Ostriker
(2015, see also ) in injecting the momentum generated during the Sedov-Taylor
phase if the cooling radius111The adopted definition of the cooling radius,
for a supernova explosion with energy $E_{\rm SN}=10^{51}$ erg, is $r_{\rm
cool}=30n_{0}^{-0.43}(Z/Z_{\odot}+0.01)^{-0.18}$ pc (e.g. Cioffi et al., 1988;
Thornton et al., 1998; Kim & Ostriker, 2014) is not resolved by at least 6
grid cells. Otherwise, we inject $10^{51}$ ergs of thermal energy (see Agertz
et al., 2015, for details) and allow for the hydrodynamic solver to track the
buildup of momentum.

###  2.1 Galaxy simulations

The simulations have galactic parameters chosen to mimic a present-day Milky
Way-mass disc galaxy. The simulations are carried out in isolation, with the
simulated galaxy positioned at the center of box with sides 600 kpc. While
this means we are neglecting cosmological effects such as mergers and inflows,
it has the benefit of allowing for high numerical resolution, as well as being
a controlled experiment. The initial conditions are based on the isolated disc
galaxy in the AGORA project (Kim et al., 2014, 2016), following the methods
described in Hernquist (1993) and Springel (2000). A Navarro-Frenk-White
profile (Navarro et al., 1996) is adopted for the dark matter halo, with a
concentration parameter $c=10$, and a virial mass222These choices yield a
virial circular velocity v200 = 150 km/s and virial radius r200 = 205 kpc. of
M${}_{200}=1.1\times 10^{12}$ M⊙.

The initial stellar and gaseous components follow exponential surface density
profiles with scale lengths $r_{\rm d}=3.4$ kpc and scale heights $z_{\rm
d}=0.1r_{\rm d}$. The gas disc is initialised at a temperature of
$T=10^{4}~{}$K. The total disc mass is $4.5\times 10^{10}~{}{\rm M_{\odot}}$,
with the fraction of mass in the gas disc varying between simulations. This
allows us to isolate the effect of the gas fraction, which changes the nature
of gravitational instability and ISM turbulence on the resulting stellar
kinematics (see Renaud et al. 2021a, hereafter Paper I). This approach and the
numerical recipe is the same as in Paper I, but we use a different, larger set
of initial conditions. The gas fraction is defined as

$f_{\rm gas}=\frac{M_{\rm gas}}{M_{\star}+M_{\rm gas}},$ (2)

where $M_{\star}$ and $M_{\rm gas}$ is the total stellar and gas mass in the
disc and the bulge, respectively. The gas fraction here acts as a proxy for
cosmological epoch, with initial disc gas fractions of 10, 20, 30, 50 and 70
per cent. The bulge-to-disc mass ratio is 0.125 and the bulge mass profile is
that of Hernquist (1990) with scale-length 0.1$r_{\rm d}$. The halo and
stellar disc are represented by 106 particles each, and the bulge consists of
105 particles.

The galaxies were simulated for 300 Myr, with the initial $\sim 100$ Myr
evolved at lower spatial resolution and without stellar feedback. This allows
for a gentle transition from the initial conditions into a new equilibrium
(see also Ejdetjärn et al. in prep.). All the analysis is carried out at the
end of the simulation by only considering the stars formed after this
relaxation period. Given the choice of initial conditions, we emphasise that
our study is relevant for understanding the origin of stellar kinematics in
typical $L^{\star}$ disc galaxies, and useful for a comparison to Milky Way-
mass disc galaxies an their progenitors.

We assign approximate cosmological epochs to our simulated galaxies by
comparing their star formation rates333Computed as the average SFR between
250-300 Myr after simulation start. (SFRs) to observationally inferred star
formation histories of Milky Way-mass galaxy progenitors by van Dokkum et al.
(2013). Dokkum et al. used deep photometric catalogs from the 3D-HST and
CANDELS Treasury surveys up to $z\sim 2.3$ to determine SFRs for galaxy
populations that match the abundance of present-day Milky Way-mass galaxies.
Each simulation hence has a corresponding redshift range which we summarise,
together with other galaxy properties, in Table 1.

Figure 2: The mean birth velocity dispersion (upper panel) and eccentricity
(lower panel) as a function of the galaxy gas fraction. The vertical error
bars correspond to the standard deviation of the measurements. The horizontal
error bars display the range of gas fractions each simulation went through.
The colored bar indicates the transition in kinematic properties around
$f_{\rm gas}\approx 20$ per cent (recently highlighted by Renaud et al.,
2021a).

##  3 Results

We begin the analysis by studying stellar and gas density maps of the
simulations at 300 Myr, shown in Figure 1. At this point in time, star
formation and feedback have decreased the initial gas fractions (computed in a
slab covering $\pm 2.5{\,\rm kpc}$ above and below the mid-plane of the
disc444Almost identical values for the gas fractions are found when
considering the cold gas only ($T\leq 10^{4}$ K).) to $9,16,21,29$ and $46$
per cent. We henceforth use this evolved gas fractions to label our
simulations: fg10, fg15, fg20, fg30, and fg45 respectively (see Table 1).

With increasing gas fraction, a larger portion of the ISM reaches high, star
forming densities (see Paper I). In the highest gas fraction galaxy, the SFR
settles to $\sim 50~{}$ M⊙ yr-1, in contrast to $\sim 1~{}$ M⊙ yr-1 in the
most gas poor case (see Table 1). The resulting increase in stellar feedback
activity in the gas rich galaxies produces vigorous galactic outflows, as is
evident from the edge-on views of the gas density fields in Figure 1. Together
with gravitational instabilities, feedback drives ISM turbulence, with higher
velocity dispersions measured in the cold star forming ISM in the gas rich
galaxies ($\sim 30{\,\rm{km\,s^{-1}}}$ in fg30) compared to the gas poor ones
($<10{\,\rm{km\,s^{-1}}}$ in fg10, see Ejdetjärn et al. in prep. for a
detailed analysis).

For $f_{\rm gas}\lesssim 20$ per cent, a spiral structure characterizes the
galaxy morphology, while massive clumps dominate the morphology at higher gas
fractions, in line with observations of gas-rich high redshift galaxies (see
e.g. Elmegreen et al., 2009; Dessauges-Zavadsky et al., 2019). This transition
in terms of morphology around a gas fraction of $\sim$ 20 per cent, for a disc
galaxy of this size and mass, is expected from the evolution of the regimes of
gravitational instability (see Paper I), and leaves distinct imprints in
stellar kinematics, which we turn to next.

###  3.1 Imprints on stellar kinematics

The upper panel of Figure 2 shows the mean vertical stellar velocity
dispersion of all stars in the disc as a function of gas fraction, at the same
simulation time as in Figure 1 ($300$ Myr). Only stars younger than 10 Myr are
considered, which allows us to quantify the imprint on stellar kinematics at
star formation, rather than secular heating.

As seen in the figure, stars born in higher gas fraction environments have a
significantly higher vertical velocity dispersion at birth, with
$\sigma_{z,{\rm birth}}$ plateauing at $\sim 25-30{\,\rm{km\,s^{-1}}}$ for
$f_{\rm gas}\gtrsim 30$ per cent. As discussed above, a gas fraction of 20 per
cent marks a transition (indicated by the vertical red band), with lower gas
fractions leading to a less turbulent, more ordered state in the disc, with
$\sigma_{z,{\rm birth}}\lesssim 10{\,\rm{km\,s^{-1}}}$ for $f_{\rm
gas}\lesssim 20$ per cent.

Figure 3: The full eccentricity distributions for the different simulations,
all stars are included.

The vertical velocity dispersion is not the only quantity that carries a
distinct imprint of the ISM conditions. The bottom panel of Figure 2 shows the
mean eccentricity555For each star, the eccentricity is computed as
$e=\frac{r_{\rm a}-r_{\rm p}}{r_{\rm a}+r_{\rm p}}$ , with $r_{\rm a}$ and
$r_{\rm p}$ being the apocenter and pericenter of an individual orbit,
respectively. of stellar orbits as a function of the gas fraction. The
vertical error bars indicate the standard deviation of the full eccentricity
distributions shown in Figure 3. Low gas fractions entail low eccentricities,
with fg10 featuring, on average, close to circular orbits with mean stellar
eccentricities $\sim 0.1$.

As the gas fraction is increased, a similar transition as previously discussed
occurs; the eccentricity distribution shifts to higher values, plateauing at
an average of $\sim 0.5$ for all simulations with $f_{\rm gas}\gtrsim 20$ per
cent. In these gas rich galaxies, the eccentricity distributions are also
significantly broader and flatter compared to the low gas fraction
counterparts (see Figure 3), with an almost equal likelihhood to find stars
with eccentricities anywhere in between 0.2 and 0.8 for $f_{\rm gas}>20$ per
cent. The origin of this marked shift towards higher eccentricities is likely
the increase in number, and masses, of dense clumps that gravitationally
scatter newly born stars onto eccentric orbits. However, a more detailed
analysis, tracing individual stars back to their birth environments, is
necessary to more conclusively characterise this process, which we leave for
future work. Without a significant change to the galactic potential or violent
scattering events in the later life of the galaxy, such as galaxy mergers,
these orbital properties will remain to be observed today as imprints of the
gas-rich, high-redshift Universe.

###  3.2 Implications for the Milky Way

####  3.2.1 The age-velocity dispersion relation

As discussed in Section 1, recent astrometric and spectroscopic surveys such
as Gaia (Gaia Collaboration et al., 2016) and APOGEE (Majewski et al., 2017)
have allowed for accurate positions and 3D velocities of millions of stars in
the Milky Way. Mackereth et al. (2019) used this information, together with
estimated stellar ages (using a Bayesian neural network model trained on
asteroseismic ages) to construct an AVR for the thick and thin discs.

The symbols in Figure 4 show the AVR from Mackereth et al. for stars in the
Solar neighbourhood (i.e. galactocentric radii 7-9 kpc). Our simulations are
shown as shaded regions with the dashed line being a non-linear least squares
fit of the form $\sigma_{z,{\rm birth}}\propto\exp(\rm{age})$. The horizontal
ranges of the simulation data refer to the approximate redshifts obtained from
matching the simulated SFRs to the observed SFR-redshift relation for Milky
Way-mass progenitors (see Table 1). We note that this is a qualitative
comparison; we do not account for any galactocentric radial dependencies in
the AVR that exist in the observational or simulation data, nor do we attempt
to account for radial migration (Schönrich & Binney, 2009) that could shape
the AVR we observe today (but see Minchev et al., 2012; Agertz et al., 2021).

As seen in Figure 4, stellar velocity dispersions imprinted at birth are not
sufficient to account for the observed trend. While this discrepancy can be
mitigated by accounting for the large observational age uncertainties for old
stars666While uncertainties on stars younger than 5 Gyr are on the order of
$\sim 30$ per cent in the method of Mackereth et al. (2019), estimated ages of
stars with true ages $>10$ Gyr can be under-estimated by as much as 3.5 Gyr.
Such large uncertainties can mitigate the discrepancies, at least at $z>1$,
between our simulations and the observations., we explore two alternative ways
of explaining the difference between the observed and simulated AVR. While
both options can contribute to the observed velocity dispersions, we consider
them one by one independently for simplicity.

  1. 1.

The method for associating redshifts to the simulations is incorrect, and
higher gas fractions were present at recent times. The corresponding redshifts
for matching the observed AVR with birth dispersions alone are shown as semi-
transparent boxes in Figure 4 (see also Table 1).

  2. 2.

Additional secular heating plays a key role in shaping the AVR.

Figure 4: The Milky Way’s AVR for the thin and thick disc (gray symbols,
Mackereth et al., 2019), together with birth velocity dispersions from the
simulations, where the associated cosmic epochs are chosen to match the Milky
Way data (semi-transparent boxes). The non-transparent boxes show the
simulated data at epochs obtained by matching SFRs to observations of Milky
Way progenitors, see main text for details. The black dashed line is an
exponential fit to this data, and by adding secular heating to this relation,
an AVR is obtained (red line) that jointly matches the Milky Way’s thick and
thin discs.

The first option cannot be ruled out. In fact, by allowing for $f_{\rm
gas}\approx 20$ and $\approx 30$ per cent at $z\sim 0.5$ and $1$,
respectively, hence matching the AVR, our results are in closer agreement with
the cold gas reservoirs observed in high redshift galaxies (see a compilation
of molecular gas fraction over time by Saintonge et al., 2013). It is
therefore possible to explain the observed AVR of the Milky Way by solely
invoking birth velocity dispersions and an evolving galactic gas fraction.
However, we note that the evolution of the gas fraction in galaxies depends on
a large number of factors, including the stellar mass and environment (e.g.
Santini et al., 2014; Schreiber et al., 2015; Genzel et al., 2015; Tacconi et
al., 2018), as well as observational selection biases. To relate the specific
gas fraction evolution of the Milky Way to observed (average) trends is
therefore non-trivial, and beyond the scope of this work.

If we instead retain the originally assigned cosmic epochs to the simulations,
secular heating via gravitational scattering of stars off GMCs and massive
clumps, spiral arms, and bars occurring over many billions of years, must
contribute to the AVR (see review by Sellwood, 2014). Aumer & Binney (2017)
used $N$-body simulations of galactic discs to study secular heating, and
demonstrated that their simulated AVRs are well fitted by a relation that
combines secular and birth velocity dispersions, of the form

$\sigma(t)=\sigma_{0}\Big{(}\frac{{\rm age}}{10\text{
Gyr}}\Big{)}^{\beta}+\sigma_{z,{\rm birth}}(t)[{\rm km/s}].$ (3)

Here $\sigma_{0}$ scales the secular heating term, $t$ is the age of the
stars, the power law index $\beta$ describes the efficiency of vertical
heating, and $\sigma_{z,{\rm birth}}(t)$ is the birth velocity dispersion of
stars. The last term is often assumed to be a constant, e.g.
$6{\,\rm{km\,s^{-1}}}$ (Aumer et al., 2016b). However, our results show that
this underestimates the contribution of birth velocity dispersions in the
Milky Way’s more gas rich past.

The red line in Figure 4 shows the total heating using Equation 3 with
$\sigma_{0}=22.2{\,\rm{km\,s^{-1}}}$, $\beta=0.65$ (compatible with the
results from Aumer et al., 2016b) and $\sigma_{z,{\rm birth}}(t)$ obtained
from our simulations (dashed line). This AVR provides a good match to the
Milky Way data, and jointly reproduces the thick and thin disc data points in
the Solar neighbourhood. This joint match cannot be achieved without the
contribution from a varying $\sigma_{z,{\rm birth}}(t)$ in Eq. 3.

We note that outliers exists to the mean relation, e.g. 6 Gyr old thick disc
stars. Such populations could be a signature of co-eval stars originating from
different birth environments (e.g. inner vs. outer disc, with different
[$\alpha$/Fe]) in terms of gas fractions, or regions with different strengths
of secular heating (Aumer et al., 2016a; Agertz et al., 2021), that
subsequently mix radially into the Solar neighbourhood (e.g. Minchev et al.,
2013; Frankel et al., 2018; Minchev et al., 2018; Mikkola et al., 2020).

In summary, we find that birth velocity dispersions alone can explain the
Milky Way’s AVR. The AVR is then intimately linked to the Milky Way’s gas
fraction evolution, and the associated level of ISM turbulence and
fragmentation physics. On the other hand, if the Milky Way’s gas fraction was
lower at early times, as suggested from matching simulated SFRs to
observations of Milky Way mass progenitors, a combination of evolving birth
velocity dispersions and secular processes must be invoked to explain the
scaling and shape of the AVR. However, the birth velocity dispersion is even
in this case always a crucial component, being responsible for at least 30 per
cent of the total heating at all epochs. Furthermore, it accounts for at least
50 per cent of the vertical kinematics for stars belonging to the old
($\gtrsim 9$ Gyr) thick disc, as well as stars younger than $\sim 1$ Gyr as
secular heating has not yet had time to significantly affect their kinematics.

Figure 5: Normalised distributions of Orbital eccentricity distributions from
Li et al. (2018) in gray for the thick (top) and thin (bottom) discs. Coloured
distributions are from two of our simulations, fg30 and fg10 respectively. The
gas rich fg30 simulation can account for the entire high eccentricity tail of
the Milky Way’s thick disc distribution but shows a deficit of low
eccentricity populations.

####  3.2.2 Formation environments of the thin and thick discs constrained by
eccentricities

Li et al. (2018) used Gaia and APOGEE data to study eccentricities of stars in
the Milky Way’s thin and thick discs. The discs were defined by their
[$\alpha$/Fe] chemical abundances, analogous to Mackereth et al. (2019).
Figure 5 shows the eccentricity distributions from Li et al. (2018) for the
thick (top) and thin (bottom) disc. The thick disc is characterised by a
higher fraction of stars on eccentric orbits in comparison to the thin disc.

The sensitivity of orbital eccentricities to the gas fraction allows us to to
constrain the formation environments of the discs. However, we refrain from
directly fitting the distributions obtained from our simulation suite to the
observations; the sensitivity to adopted galactocentric radial cuts and
unaccounted processes, such as radial migration over a long period, makes such
an approach degenerate. Instead, we overlay the eccentricity distribution from
fg30 (‘gas rich disc’) and fg10 (‘gas poor disc’) to perform a more
qualitative comparison. Only stars present at galactocentric radii
$7-9~{}{\,\rm kpc}$ in the simulations are considered.

As seen in Figure 5, the tail of high eccentricities ($\gtrsim 0.4$) in the
thick disc distribution is well matched by stars that have formed in the gas
rich disc. This means that cosmological processes, such as mergers, are not
necessary to explain such high eccentricities. We re-emphasize that this is
not due to slow secular processes (e.g. radial migration, Sales et al., 2009),
but an imprint of the turbulent and clumpy ISM. However, to explain the
existence of eccentricities $\lesssim 0.1-0.2$, lower gas fractions are
required. This may indicate that the thick disc’s gas fraction decreased
rapidly in the early Universe, while retaining a high [$\alpha$/Fe]. The
redshifts obtained from our SFR matching are compatible with this notion, as
$f_{\rm gas}$ decreases from 50 to 20 per cent in less than a Gyr at $z\sim
1.5$ (see Table 1). Alternatively, as discussed in Section 3.2.1, different
gas fractions may have been present simultaneously in different regions of the
galaxy when the [$\alpha$/Fe]-enhanced thick disc formed.

The thin disc distribution is well matched by the gas poor disc for
eccentricities $<0.2$, but higher gas fractions are required to fit the entire
distribution. However, it is unlikely for the thin disc to have formed any
significant fraction of its stars when $f_{\rm gas}\gtrsim 20$ per cent: such
environments produce too many stars with higher than observed eccentricities.

The existence of a population of low eccentricity stars, originating from low
gas fraction conditions, in both the thin and thick disc, is compatible with
them having overlapping formation times. Indeed, Beraldo e Silva et al. (2021)
combined stellar age estimates, orbital eccentricities and [$\alpha$/Fe] from
APOGEE and Gaia, to argue for a co-eval formation scenario $9-10$ Gyr ago. At
this epoch, the Milky Way’s inner and outer disc must have had similar gas
fractions, but different levels of $\alpha-$enhancement. Such an inner-outer
disc dichotomy, in terms of [$\alpha$/Fe], is recovered in cosmological
simulations of galaxy formation around the time of the last major merger
(Agertz et al., 2021; Renaud et al., 2021b; Renaud et al., 2021c), with
subsequent radial mixing bringing stars from both regions into the Solar
neighbourhood.

##  4 Discussion

The physical origin of the diversity of orbital and kinematical properties
observed in the stellar populations in the solar neighbourhood has been a long
lasting debate in the literature. Models and simulations have emphasised the
role of dynamical heating by a bar and/or spiral arms (e.g. Grand et al.,
2016; Aumer et al., 2016a), or via radial migration (Schönrich & Binney,
2009). Cosmological simulations also point out that the transition from a
merger-dominated growth of the galaxy to a more quiescent, secular evolution
coincides with the transition from formation in the turbulent thick disc,
subject to repeated tidal excitations to that in the dynamically colder thin
disc (Renaud et al., 2021b; Bird et al., 2021). However, simulations run
without the cosmological context have already pointed out comparable
transitions driven by the intrinsic evolution of the regime of star formation
(e.g. Clarke et al., 2019; Khoperskov et al., 2021). Our results complement
this by showing that the shapes of the AVR and the distribution of orbital
eccentricities can be explained by the intrinsic, natural, evolution of the
gas content of the galaxy, without the need to invoke secular dynamical
processes like migration, nor external effects like mergers.

The decrease of the gas fraction expected along galaxy evolution changes the
regime of disc instabilities and, in turn, of the properties of the ISM
turbulence (Paper I, Ejdetjärn et al. in prep). A similar gas fraction
dependence is expected777Hayward & Hopkins (2017) showed that for $f_{\rm
gas}\gtrsim 30$ per cent, feedback is likely to be trapped in the dense ISM,
leading to high levels of turbulence. At lower gas fractions, feedback ejecta
is able to leak out of the galaxy, hence not injecting momentum into the ISM,
with lower levels of turbulence as a result. when stellar feedback is the
dominant turbulence driver (Hayward & Hopkins, 2017). As such, the existence
thin galactic discs is intimately connected to the evolution of galactic
baryon fractions (see e.g. Agertz & Kravtsov, 2016) and the way in which gas
is accreted and lost in feedback driven outflows. Yu et al. (2021) used
cosmological simulations of Milky Way-mass galaxies to show that thin galactic
discs cannot form unless the gas cooling time in the circum-galactic medium
(CGM) is longer than the dynamical time (see also Stern et al., 2021). The
authors argued that for short CGM cooling times, gas is accreted more rapidly,
leading to burstier star formation and thick disc formation. In light of our
work, it is not necessarily star formation bursts that lead to thick discs,
but rather the high levels of ISM turbulence and the presence of massive
clumps, that in turn imprint thick disc stellar kinematics. In this picture,
the state of the CGM is a crucial ingredient, as it dictates, together with
the strength of galactic outflows, the ability to form and maintain gas rich
turbulent discs.

Our result serve as a proof-of-concept that such an “historical”
reconstruction can suffice to reproduce the observables. However, key aspects
neglected here for simplicity would likely alter our conclusions. In
particular, our simulations do not represent different snapshots along the
evolution of the same galaxy, as the mass and size are unchanged, while a real
galaxy would grow. A major contribution to this growth is that of galaxy
mergers, which are know to stir the ISM and boost the turbulence precisely at
the epochs the most fertile in star formation. By not capturing this aspect,
our model underestimate the velocity dispersion of the stars born during the
merger-dominated phase (e.g. $z\gtrsim 1$ for the Milky Way).

Finally, the accretion of gas along cosmological filaments, also neglected in
our work, is suspected to influence the injection scale of turbulence (Agertz
et al, in prep.), alter the relation between gas fraction and velocity
dispersion, and thus skew the AVR. Therefore, our results should be
interpreted as an illustration of the importance of the secular evolution of
the gas content on the dynamics of the stellar populations, among other coeval
processes.

The sharp transition seen for gas fractions $\approx 20$ per cent constitutes
a strong indication of how fast this process contributes to the switch between
dynamically hot to cold populations. We note that, interestingly, in the Milky
Way such a transition is estimated to occur approximately at the same epoch as
the last major merger, and that the two changes of regimes likely conspire in
creating the observed dynamical bimodality of stellar populations.

##  5 Conclusions

Using a series of hydrodynamical simulations of isolated, Milky Way-mass disc
galaxies, we study the impact of the gas fraction on the stellar kinematics.
Our main conclusions are as follows.

  * •

Newly formed stars in disc galaxies with high gas fractions feature high
vertical velocity dispersions. This is a direct imprint of the turbulent ISM.
Simulations with $f_{\rm gas}\gtrsim 30$ per cent generate vertical stellar
velocity dispersions $\sigma_{z}\sim$ 30 km/s already at star formation, well
above the observed kinematics in the Milky Way’s thin disc, but compatible
with the thick disc (Mackereth et al., 2019). In contrast, simulations with
$f_{\rm gas}\lesssim 20$ per cent form stars with $\sigma_{z}\lesssim$ 10
km/s, in line with young stars in the Milky Way’s thin disc.

  * •

The $\sigma_{z}$ imprinted by the ISM is in principle sufficient for
explaining the Milky Way’s age-velocity dispersion relation (AVR). However, by
assigning a lookback time to the simulations, by matching their SFRs to those
of observed Milky Way progenitors, additional secular heating is likely
required to explain the AVR. However, the birth velocity dispersion is always
a crucial component for total level of stellar heating, especially in the
early Universe ($>9-10$ Gyr in the past).

  * •

Orbital eccentricities of young stars depend on the disc gas fraction. Stars
born in galaxies with $f_{\rm gas}\sim 10$ per cent feature stars on mostly
circular orbits, with mean orbital eccentricities $\sim 0.1$. As the gas
fraction is increased, stronger gravitational scattering allows for wider
eccentricity distributions with higher mean values. For $f_{\rm gas}>20$ per
cent, we find a ‘saturation effect’, with similar (wide and flat) eccentricity
distributions, and mean eccentricities of young stars $\sim$ 0.5.

  * •

We find that the high eccentricity tail of the Milky Way’s thick disc
distribution in the Solar neighbourhood can be fully explained with star
formation in a turbulent high gas fraction disc ($f_{\rm gas}>30$ per cent).
Both the thin and the thick disc require some contribution from a gas-poor
regime to fully characterise their eccentricity distributions. This notion is
compatible with a thin/thick disc co-formation scenario (e.g. Agertz et al.,
2021; Renaud et al., 2021b; Renaud et al., 2021c; Beraldo e Silva et al.,
2021) in which a range of ISM conditions existed simultaneously, with radial
migration later mixing stars across the galaxy.

We conclude that the decrease of gas fraction, as expected along galaxy
evolution, changes the regime of disc instabilities, properties of the ISM
turbulence, and in turn the imprinted kinematical properties of newly formed
stellar populations. A key result is that this is not a gradual change, but
rather a rapid transition in kinematical properties around a gas fraction
$\sim 20$ per cent.

Although our controlled experiments of Milky Way-sized galaxies highlight the
impact of a specific parameter on the emerging stellar dynamics, other aspects
like the disc growth, gas accretion, and galaxy mergers would also play
important roles. Our conclusions thus call for further studies where the
impact of such processes on a galaxy’s AVR and orbital eccentricities can be
quantified.

## Acknowledgements

We acknowledge support from the Knut and Alice Wallenberg Foundation and the
Royal Physiographic Society of Lund. OA acknowledges support from the Swedish
Research Council grant 2019-04659. Simulations were performed in part using
computational resources at LUNARC, the centre for scientific and technical
computing at Lund University, on the Swedish National Infrastructure for
Computing (SNIC) allocation 2018/3-649, as well as allocation LU 2018/2-28
thanks to financial support from the Royal Physiographic Society of Lund.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abadi et al. (2003) Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 597, 21
  * Agertz & Kravtsov (2016) Agertz O., Kravtsov A. V., 2016, ApJ, 824, 79
  * Agertz et al. (2013) Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25
  * Agertz et al. (2015) Agertz O., Romeo A. B., Grisdale K., 2015, MNRAS, 449, 2156
  * Agertz et al. (2020) Agertz O., et al., 2020, MNRAS, 491, 1656
  * Agertz et al. (2021) Agertz O., et al., 2021, MNRAS, 503, 5826
  * Aumer & Binney (2017) Aumer M., Binney J., 2017, MNRAS, 470, 2113
  * Aumer et al. (2016a) Aumer M., Binney J., Schönrich R., 2016a, MNRAS, 459, 3326
  * Aumer et al. (2016b) Aumer M., Binney J., Schönrich R., 2016b, MNRAS, 462, 1697
  * Bensby et al. (2011) Bensby T., Alves-Brito A., Oey M. S., Yong D., Meléndez J., 2011, ApJ, 735, L46
  * Bensby et al. (2014) Bensby T., Feltzing S., Oey M. S., 2014, A&A, 562, A71
  * Beraldo e Silva et al. (2021) Beraldo e Silva L., Debattista V. P., Nidever D., Amarante J. A. S., Garver B., 2021, MNRAS, 502, 260
  * Bird et al. (2021) Bird J. C., Loebman S. R., Weinberg D. H., Brooks A. M., Quinn T. R., Christensen C. R., 2021, MNRAS, 503, 1815
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529
  * Bournaud & Elmegreen (2009) Bournaud F., Elmegreen B. G., 2009, ApJ, 694, L158
  * Brook et al. (2004) Brook C. B., Kawata D., Gibson B. K., Freeman K. C., 2004, ApJ, 612, 894
  * Bryant et al. (2015) Bryant J. J., et al., 2015, MNRAS, 447, 2857
  * Bundy et al. (2015) Bundy K., et al., 2015, ApJ, 798, 7
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Cioffi et al. (1988) Cioffi D. F., McKee C. F., Bertschinger E., 1988, ApJ, 334, 252
  * Clarke et al. (2019) Clarke A. J., et al., 2019, MNRAS, 484, 3476
  * Dessauges-Zavadsky et al. (2019) Dessauges-Zavadsky M., et al., 2019, Nature Astronomy, 3, 1115
  * Elmegreen et al. (2009) Elmegreen D. M., Elmegreen B. G., Marcus M. T., Shahinyan K., Yau A., Petersen M., 2009, ApJ, 701, 306
  * Feuillet et al. (2019) Feuillet D. K., Frankel N., Lind K., Frinchaboy P. M., García-Hernández D. A., Lane R. R., Nitschelm C., Roman-Lopes A. r., 2019, MNRAS, 489, 1742
  * Frankel et al. (2018) Frankel N., Rix H.-W., Ting Y.-S., Ness M., Hogg D. W., 2018, ApJ, 865, 96
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Genzel et al. (2015) Genzel R., et al., 2015, ApJ, 800, 20
  * Glazebrook (2013) Glazebrook K., 2013, Publ. Astron. Soc. Australia, 30, e056
  * Gnerucci et al. (2011) Gnerucci A., et al., 2011, A&A, 528, A88
  * Grand et al. (2016) Grand R. J. J., Springel V., Gómez F. A., Marinacci F., Pakmor R., Campbell D. J. R., Jenkins A., 2016, MNRAS, 459, 199
  * Grisdale et al. (2017) Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093
  * Grisdale et al. (2018) Grisdale K., Agertz O., Renaud F., Romeo A. B., 2018, MNRAS,
  * Grisdale et al. (2019) Grisdale K., Agertz O., Renaud F., Romeo A. B., Devriendt J., Slyz A., 2019, arXiv e-prints,
  * Guillet & Teyssier (2011) Guillet T., Teyssier R., 2011, Journal of Computational Physics, 230, 4756
  * Gustafsson et al. (2016) Gustafsson B., Church R. P., Davies M. B., Rickman H., 2016, A&A, 593, A85
  * Hayden et al. (2015) Hayden M. R., et al., 2015, ApJ, 808, 132
  * Hayden et al. (2020) Hayden M. R., et al., 2020, MNRAS, 493, 2952
  * Hayward & Hopkins (2017) Hayward C. C., Hopkins P. F., 2017, MNRAS, 465, 1682
  * Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  * Hernquist (1993) Hernquist L., 1993, ApJS, 86, 389
  * Hockney & Eastwood (1981) Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles. New York: McGraw-Hill, 1981
  * Kassin et al. (2012) Kassin S. A., et al., 2012, ApJ, 758, 106
  * Kazantzidis et al. (2009) Kazantzidis S., Zentner A. R., Kravtsov A. V., Bullock J. S., Debattista V. P., 2009, ApJ, 700, 1896
  * Khoperskov et al. (2021) Khoperskov S., Haywood M., Snaith O., Di Matteo P., Lehnert M., Vasiliev E., Naroenkov S., Berczik P., 2021, MNRAS, 501, 5176
  * Kim & Ostriker (2014) Kim C.-G., Ostriker E. C., 2014, preprint,  (arXiv:1410.1537)
  * Kim & Ostriker (2015) Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99
  * Kim et al. (2014) Kim J.-h., et al., 2014, ApJS, 210, 14
  * Kim et al. (2016) Kim J.-h., et al., 2016, ApJ, 833, 202
  * Krumholz & Tan (2007) Krumholz M. R., Tan J. C., 2007, ApJ, 654, 304
  * Lee et al. (2016) Lee E. J., Miville-Deschênes M.-A., Murray N. W., 2016, ApJ, 833, 229
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Li et al. (2018) Li C., Zhao G., Zhai M., Jia Y., 2018, ApJ, 860, 53
  * Loebman et al. (2011) Loebman S. R., Roškar R., Debattista V. P., Ivezić Ž., Quinn T. R., Wadsley J., 2011, ApJ, 737, 8
  * Mackereth et al. (2019) Mackereth J. T., et al., 2019, MNRAS, 489, 176
  * Majewski et al. (2017) Majewski S. R., et al., 2017, AJ, 154, 94
  * Martig et al. (2014) Martig M., Minchev I., Flynn C., 2014, MNRAS, 443, 2452
  * Martizzi et al. (2015) Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504
  * Mieda et al. (2016) Mieda E., Wright S. A., Larkin J. E., Armus L., Juneau S., Salim S., Murray N., 2016, ApJ, 831, 78
  * Mikkola et al. (2020) Mikkola D., McMillan P. J., Hobbs D., 2020, MNRAS, 495, 3295
  * Minchev et al. (2012) Minchev I., Famaey B., Quillen A. C., Dehnen W., Martig M., Siebert A., 2012, A&A, 548, A127
  * Minchev et al. (2013) Minchev I., Chiappini C., Martig M., 2013, A&A, 558, A9
  * Minchev et al. (2018) Minchev I., et al., 2018, MNRAS, 481, 1645
  * Murray (2011) Murray N., 2011, ApJ, 729, 133
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Nordström et al. (2004) Nordström B., et al., 2004, A&A, 418, 989
  * Quinn et al. (1993) Quinn P. J., Hernquist L., Fullagar D. P., 1993, ApJ, 403, 74
  * Raiteri et al. (1996) Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105
  * Read et al. (2008) Read J. I., Lake G., Agertz O., Debattista V. P., 2008, MNRAS, 389, 1041
  * Reddy et al. (2003) Reddy B. E., Tomkin J., Lambert D. L., Allende Prieto C., 2003, MNRAS, 340, 304
  * Renaud et al. (2021a) Renaud F., Romeo A. B., Agertz O., 2021a, arXiv e-prints, p. arXiv:2106.00020
  * Renaud et al. (2021b) Renaud F., Agertz O., Read J. I., Ryde N., Andersson E. P., Bensby T., Rey M. P., Feuillet D. K., 2021b, MNRAS, 503, 5846
  * Renaud et al. (2021c) Renaud F., Agertz O., Andersson E. P., Read J. I., Ryde N., Bensby T., Rey M. P., Feuillet D. K., 2021c, MNRAS, 503, 5868
  * Romeo (1992) Romeo A. B., 1992, MNRAS, 256, 307
  * Romeo & Agertz (2014) Romeo A. B., Agertz O., 2014, MNRAS, 442, 1230
  * Romeo et al. (2008) Romeo A. B., Agertz O., Moore B., Stadel J., 2008, ApJ, 686, 1
  * Rosen & Bregman (1995) Rosen A., Bregman J. N., 1995, ApJ, 440, 634
  * Saintonge et al. (2013) Saintonge A., et al., 2013, ApJ, 778, 2
  * Sales et al. (2009) Sales L. V., et al., 2009, MNRAS, 400, L61
  * Santini et al. (2014) Santini P., et al., 2014, A&A, 562, A30
  * Schönrich & Binney (2009) Schönrich R., Binney J., 2009, MNRAS, 399, 1145
  * Schreiber et al. (2015) Schreiber C., et al., 2015, A&A, 575, A74
  * Sellwood (2014) Sellwood J. A., 2014, Reviews of Modern Physics, 86, 1
  * Sellwood & Binney (2002) Sellwood J. A., Binney J. J., 2002, MNRAS, 336, 785
  * Springel (2000) Springel V., 2000, MNRAS, 312, 859
  * Stern et al. (2021) Stern J., et al., 2021, ApJ, 911, 88
  * Sutherland & Dopita (1993) Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253
  * Tacconi et al. (2018) Tacconi L. J., et al., 2018, ApJ, 853, 179
  * Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
  * Thornton et al. (1998) Thornton K., Gaudlitz M., Janka H.-T., Steinmetz M., 1998, ApJ, 500, 95
  * Ting & Rix (2019) Ting Y.-S., Rix H.-W., 2019, ApJ, 878, 21
  * Wisnioski et al. (2015) Wisnioski E., et al., 2015, ApJ, 799, 209
  * Yu & Liu (2018) Yu J., Liu C., 2018, MNRAS, 475, 1093
  * Yu et al. (2021) Yu S., et al., 2021, MNRAS, 505, 889
  * de Jong et al. (2019) de Jong R. S., et al., 2019, The Messenger, 175, 3
  * van Dokkum et al. (2013) van Dokkum P. G., et al., 2013, ApJ, 771, L35