# Supernova-driven turbulent metal mixing in high redshift galactic disks:
metallicity fluctuations in the interstellar medium and its imprints on metal
poor stars in the Milky Way

Anne Noer Kolborg Niels Bohr Institute, University of Copenhagen,Blegdamsvej
17, DK-2100 Copenhagen, Denmark  Department of Astronomy & Astrophysics,
University of California, Santa Cruz, CA 95064, USA  Davide Martizzi Niels
Bohr Institute, University of Copenhagen,Blegdamsvej 17, DK-2100 Copenhagen,
Denmark  Department of Astronomy & Astrophysics, University of California,
Santa Cruz, CA 95064, USA  Enrico Ramirez-Ruiz Department of Astronomy &
Astrophysics, University of California, Santa Cruz, CA 95064, USA  Niels Bohr
Institute, University of Copenhagen,Blegdamsvej 17, DK-2100 Copenhagen,
Denmark  Hugo Pfister Department of Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong, China  Niels Bohr Institute, University of
Copenhagen,Blegdamsvej 17, DK-2100 Copenhagen, Denmark  Charli Sakari
Department of Physics & Astronomy, San Francisco State University, 1600
Holloway Avenue, San Francisco, CA 94132 USA  Risa H. Wechsler Kavli Institute
for Particle Astrophysics and Cosmology and Department of Physics, Stanford
University, 382 Via Pueblo Mall, Stanford, CA 94305, USA  Kavli Institute for
Particle Astrophysics and Cosmology, SLAC National Accelerator Laboratory,
2575 Sand Hill Road, Menlo Park, CA 94025, USA

###### Abstract

The extent to which turbulence mixes gas in the face of recurrent infusions of
fresh metals by supernovae (SNe) could help provide important constraints on
the local star formation conditions. This includes predictions of the
metallicity dispersion amongst metal poor stars, which suggests that the
interstellar medium (ISM) was not very well mixed at these early times. The
purpose of this Letter is to help isolate via a series of numerical
experiments some of the key processes that regulate turbulent mixing of SN
elements in galactic disks. We study the gas interactions in small simulated
patches of a galaxy disk with the goal of resolving the small-scale mixing
effects of metals at parsec scales. By investigating the statistics of
variations of $\alpha$ elements in these simulations we are able to derive
meaningful constraints on the star formation history (SFH) of the Milky Way
(MW). We argue that the observed dispersion in stellar abundances of metal
poor halo stars are naturally explained by star formation conditions expected
in dwarf satellites or in an early low starforming MW progenitor. Metal
variations in phase-mixed stars give further credence to the idea that the MW
halo could have been assembled primarily by disrupted dwarf satellites.

††journal: ApJL

tablenum SIXtablenum

##  1 introduction

The distribution of chemical abundances in stars serves as a “fossil record”
of the galaxy’s evolutionary history and provides clues about the types of
nucleosynthetic processes that occurred early in the history of the Milky Way
(MW). Of particular interest for corecollapse enrichment is the $\alpha$
element composition of Galactic halo stars. In the metallicity regime [Fe/H]
$\approx-4.0$ to $-2.0$ (e..g Audouze & Silk, 1995; Macias & Ramirez-Ruiz,
2018), $\alpha$ elements have been found to be pure core-collapse products
with a fairly characteristic [$\alpha$/Fe] abundance ratio but with a star to
star bulk scatter in their [$\alpha$/H] concentrations (Cayrel et al., 2004;
Ryan et al., 1996; Frebel & Norris, 2015). The presence of these elements
demonstrates that their synthesis operated in a fairly robust manner in the
early MW while their abundance dispersions should yield clues about the
formation rate of the massive stellar progenitor systems that lived and
dwindled before them. The reason is that metal injection rates can help
establish how chemically unmixed and inhomogeneous the interstellar gas was
when these stars formed (Audouze & Silk, 1995; Krumholz & Ting, 2018). At
later times, these localized inhomogeneities are expected to be smoothed out
as more events take place and corecollapse supernova (SN) products are given
more time to migrate and mix throughout the progenitor galaxy (de Avillez &
Mac Low, 2002; Pan et al., 2013; Krumholz & Ting, 2018). To this end, a
comprehensive understanding of the inhomogeneous enrichment of the
interstellar medium (ISM) holds great promise for our attempts to decipher the
MW’s stellar assembly history (e.g. Wang et al., 2021) and, in particular, the
origin of halo stars. The relative contribution between stars formed in situ
(e.g. Eggen et al., 1962) and those formed ex-situ (e.g. Searle & Zinn, 1978)
to the halo of the Milky Way has been long debated. While recent studies of
the respective importance of these different star-formation channels to the
halo have primarily supported ex-situ assembly (e.g. Robertson et al., 2005;
Johnston et al., 2008), a comprehensive picture remains elusive and a series
of issues are still pressing. In particular related to the specific merger and
SFHs of the individual satellite galaxies.

Comprehensive numerical modeling aimed at addressing the mixing of metals in
the MW requires a suite of codes to treat the diverse range of physics and
timescales at play. Simulations performed in a cosmological context (e.g. Few
et al., 2014; Shen et al., 2017; Naiman et al., 2018) are required to capture
the large-scale mixing mechanisms including satellite mergers, gas inflows,
galactic winds and fountains, and shearing disk instabilities, all of which
are unquestionably important to comprehend the enrichment and dissemination of
metals. On the other hand, since cosmological simulations today are limited by
resolutions of a few tens to hundreds of parsecs, they inevitably entail “sub-
grid” models for turbulent mixing, star formation, and stellar feedback. As
such, ample uncertainties still exist at the hundred parsec scale. For
turbulent mixing the trade-offs are even more pronounced and require a large
range of scales to be resolved simultaneously (e.g. Colbrook et al., 2017; de
Avillez & Mac Low, 2002; Yang & Krumholz, 2012). As a consequence, a single
simulation of the full problem incorporating all of the aforementioned effects
would not only be prohibitively expensive, but also difficult to interpret
because of the complexity of the interplay between the various physical
mechanisms at different scales. Instead, in this Letter we study the
inhomogeneous enrichment of the ISM via a series of numerical experiments that
isolate the key processes that regulate the mixing of SN elements.

In addition to being computationally feasible, this approach enables a
thorough understanding of the relevant processes. To this end we present
models that study the gas interactions in small patches of a galaxy disk in
order to study the small-scale mixing effects of metals at parsec scales from
SN into the ISM. Investigating the distribution of $\alpha$ elements in these
simulations and comparing them extensively with the observational data can, in
turn, provide important constraints on the starforming conditions (Audouze &
Silk, 1995; Armillotta et al., 2018) within the early MW and/or the infalling
satellites that builtup the stellar halo. Given the relatively nascent stage
of modeling in the field, the results we present in this Letter amount to a
substantive improvement in our understanding of the physical mechanism leading
to the inhomogeneous enrichment of the ISM and the mixing of $\alpha$ elements
in the early MW.

Galaxy type | $\Sigma_{\text{gas}}$ | $\dot{\Sigma}_{\text{SFR}}$ | $\Gamma$ | $f_{\rm gas}$ | $\rho_{0}$ | $L$ | $z_{\rm eff}$ | $\langle\sigma_{v}\rangle$ | $t_{r}$ | $\kappa$
---|---|---|---|---|---|---|---|---|---|---
| [$\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{p}\mathrm{c}^{2}$] | [$\text{\,}\mathrm{M}_{\odot}\mathrm{/}\mathrm{k}\mathrm{p}\mathrm{c}^{2}\mathrm{/}\mathrm{M}\mathrm{y}\mathrm{r}$] | [$\text{\,}\mathrm{/}\mathrm{k}\mathrm{p}\mathrm{c}^{2}\mathrm{/}\mathrm{M}\mathrm{y}\mathrm{r}$] |  | [$\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{3}$] | [pc] | [pc] | [km/s] | [Myr] | [pc km/s]
MW progenitor, high SFR | 50 | $3\text{\times}{10}^{4}$ | $300$ | 0.100 | $3.47\text{\times}{10}^{-23}$ | $1000$ | 40 | 11 | 19 | 147
MW progenitor, low SFR | 5 | $1\text{\times}{10}^{3}$ | $10$ | 0.088 | $2.08\text{\times}{10}^{-24}$ | $1000$ | 80 | 8 | 39 | 213
Satellite galaxy | 5 | $8\text{\times}{10}^{2}$ | $8$ | 0.496 | $4.67\text{\times}{10}^{-25}$ | $4000$ | 365 | 12 | 120 | 1460
MW progenitor, low $\epsilon_{*}$ | 5 | $2\text{\times}{10}^{2}$ | $2$ | 0.088 | $2.08\text{\times}{10}^{-24}$ | $1000$ | 80 | 4 | 73 | 107
High redshift, high $\epsilon_{*}$ | 5 | $5.5\text{\times}{10}^{3}$ | $55$ | 0.088 | $2.08\text{\times}{10}^{-24}$ | $1000$ | 80 | 14 | 23 | 373
Table 1: Parameters for all the galaxy models examined in this project,
columns are: simulation name, gas surface density, star formation rate surface
density, surface density of rate of SNe
($\Gamma=\dot{\Sigma}_{\text{SFR}}/100M_{\odot}$), fraction of gas mass to
stellar mass, initial density of gas in the mid plane of the disk, simulation
box size, effective scale height of the disk, mean velocity dispersion of the
gas, characteristic relaxation time and turbulent diffusion coefficient. The
results of the supporting simulations (in grey) are presented in Appendix B.

##  2 Methods and Metal Mixing Experiments

We use the hydrodynamical, adaptive mesh refinement (AMR) code RAMSES
(Teyssier, 2002) to model the gas interactions in small patches of a galaxy
disk to study the small-scale mixing effects of metals from supernovae (SNe)
into the ISM. The simulation set-up is described in detail in Martizzi et al.
(2016) and here we summarize the salient features and provide details for
those aspects of the set-up that are different to their investigations.

The set-up studies small patches of galaxy disks at relatively high
resolution. In particular, we are able to capture the momentum injection and
thermal energy of individual SNe as the cooling phase of individual remnants
is well resolved. By changing the parameters of the gravitational potential,
gas fraction, gas density, and star formation efficiency, we are also able to
replicate the characteristics of different local galactic environments that
might be representative of the early MW or its accreted satellites.

For the purposes of this Letter we have chosen to focus on two different
possible MW progenitors as well as one likely satellite candidate. The MW
progenitors are chosen to model two different high redshift cases: one gas-
rich, highly starforming model (van Dokkum et al., 2013; Shen et al., 2015;
Naiman et al., 2018). The other one has a lower star formation rate (SFR) and
is motivated in particular by the recent findings of Wang et al. (2021) that
the presence of the Large Magellanic Cloud (LMC) satellite hints at a MW
progenitor with a less vigorous SFR. Finally, we include a dwarf-like
satellite galaxy in order to examine the viability of an ex situ birth
environment for the halo stars. In order to study how mixing is altered by
changes in the gas-to-light ratio, we have also run two additional versions of
the low SFR MW progenitor model with the same gas density but with two
different star formation efficiencies ($\epsilon_{\star}$). The results of
these additional simulations are included in Appendix B. Table 1 summarizes
the parameters of all the models used in this work.

The size of the computational domain changes between different models. It
increases for the satellite galaxy simulation to accommodate the larger scale
height of the disk which is caused by the reduced gravity of the galaxy. We
use a regular Cartesian grid (all cells within the box have the same physical
size, but not the same mass resolution) and fix the resolution at $2^{8}$
cells. This is sufficient to ensure that the cooling radius is resolved by at
least a few resolution elements for the majority of the supernova remnants
(SNR) in all models. The cooling radius is calculated using (Martizzi et al.,
2016):

$R_{\text{c}}\approx$14\text{\,}$\left(\frac{n_{\text{H}}}{$1\text{\,}\mathrm{c}\mathrm{m}^{3}$}\right)^{-3/7}\left(\frac{E_{\text{SNe}}}{10^{51}{\rm
erg}}\right)^{2/7}\left(\frac{Z}{Z_{\odot}}\right)^{-1/7}{\rm pc}$ (1)

In each model the percentage of SNe that are not resolved by at least 2 (5)
cells at their cooling size are: MW progenitor, high SFR 0.6% (6.1%), MW
progenitor low, SFR 0% (1.6%), and satellite galaxy 0% (3.8%).

The SNe are seeded following the ”fixed” distribution scheme used in Martizzi
et al. (2016); in this framework a scale height for explosions is defined and
the SNe have equal probability of happening anywhere in the volume contained
within this height and zero probability outside. In all the models used here,
this height is twice the gassous scale height of the disk
($z_{\text{SNe}}=2z_{\text{eff}}$), $z_{\text{eff}}$ is defined in Section
3.1. Each SN injects $M_{\text{ej}}=6.8M_{\odot}$ and energy $E_{\rm
SNe}=10^{51}{\rm erg}$, and this energy is partitioned into thermal energy and
momentum injection using the sub-grid model of Martizzi et al. (2016), which
also distributes the mass among the neighboring cells of the explosion. The
rate of SNe per unit area, $\Gamma$, is given by the SFR as
$\Gamma=\frac{\dot{\Sigma}_{\text{SFR}}}{$100\text{\,}\mathrm{M}_{\odot}$}$.

The evolutionary timescale ($t$) of each simulation is dictated by the SFR and
differs between simulations in order to account for the gentler evolution of
galaxies with lower surface densities of gas (efficient cooling is impaired)
and lower rate of SNe (feedback is reduced). The time between consecutive
snapshots ($\delta t$) is also adjusted to ensure that there are not too many
SNe happening in each snapshot, i.e. models with higher rates of SNe have
lower $\delta t$. This allows us to follow the interaction between the remnant
and the ISM more closely. The specific time intervals are: MW progenitor, high
SFR: $\delta t=$
$5\text{\times}{10}^{4}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$, MW
progenitor, low SFR: $\delta
t=$$5\text{\times}{10}^{4}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$,
and satellite: $\delta
t=$$1\text{\times}{10}^{5}\text{\,}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}$.

Figure 1: Evolution of the velocity dispersion, $\sigma_{v}$ (blue) and the
mean [Fe/H] abundance (purple line) as well as 1 and 3 $\sigma$ spread (shaded
regions) as a function of time. The relaxation time, $t_{r}$, for each model
is marked by the vertical gray line in each panel. The panel titles indicate
the simulation. Models with lower SFR have a more slowly evolving mean [Fe/H]
but tend to have smaller spread at same mean abundance as models with higher
SFR.

We monitor individual elements as passive scalars. For the purposes of this
Letter we focus on Mg (as our representative $\alpha$ element) and Fe. Metals
are introduced into the gas by SNe, which are assumed to be chemically
identical.111The reader is refer to Section 4.2 for a discussion on the
implications of our findings when relaxing this assumption. This allows us to
isolate the effects that turbulent mixing has on the abundance spread of
metals within the interstellar medium. We assume that each SN yields a mass of
$M_{\text{Mg}}=0.10M_{\odot}$ and $M_{\text{Fe}}=0.08M_{\odot}$. These values
are chosen to achieve a yield ratio of
$\frac{M_{\text{Mg}}}{M_{\text{Fe}}}=1.27$, equivalent to the best fitting
value for the observational set we wish to compare with (see Appendix A).
Fixing the Fe to Mg mass ratio allows us to directly compare the spread of
[Mg/H] in the simulations as a function of the mean [Fe/H] abundance and
compare it with that derived from observations. This approach is similar to
assuming that the mean abundances are fixed by the mean yields, which is a
commonly used conjecture when studying metal mixing (e.g. Dalcanton, 2007).

##  3 Simulating Supernova Feedback and Mixing in galaxy disk patches.

###  3.1 Steady State

The early time of the simulations are characterized by a relaxation from the
initial conditions. Initially the disk is supported by thermal pressure while
at later times the disk is supported in part by turbulent pressure. As a
result, the true effective scale height differs slightly from what one would
calculate assuming hydrostatic equilibrium with the unchanging gravitational
potential. To address this Martizzi et al. (2016) defined the effective scale
height, $z_{\text{eff}}$ as:

$z_{\text{eff}}=\frac{\Sigma_{\text{gas}}}{2\rho_{0}}=\frac{1}{2\rho_{0}}\int_{-L/2}^{L/2}\rho_{\text{gas}}(z)\,dz$

where $\Sigma_{\text{gas}}$ is the surface density of gas and $\rho_{0}$ is
the gas density in the disk midplane during initial conditions. Both
$\rho_{0}$ and $z_{\text{eff}}$ are listed in Table 1 for all models.
We use $z_{\text{eff}}$ to define a characteristic relaxation time, $t_{r}$,
for our disk models:

$t_{r}=\frac{4z_{\text{eff}}}{\langle\sigma_{v}\rangle},$ (2)

$\langle\sigma_{v}\rangle$ is the time averaged, mass-weighted velocity
dispersion of the gas. We calculate the mass weighted velocity dispersion as:

$\sigma_{v}=\sqrt{\frac{\sum_{i}^{N}(\vec{v_{i}}-\langle\vec{v}\rangle)^{2}M_{i}}{\sum_{i}^{N}M_{i}}},$

where $N$ is the number of cells in the simulation volume, $M_{i}$ is the mass
of gas in the cell, $\vec{v}_{i}$ is the gas 3D velocity, and $\langle
v\rangle$ is the mean velocity of the gas. The parameter
$\langle\sigma_{v}\rangle$ is set by the properties of the turbulence in the
disk and depends weakly on the surface density of gas (Martizzi et al., 2016).
The evolution of the velocity dispersion is shown in the top panels of Figure
1. We calculate the mean of the velocity dispersion over the entire evolution
time of the model and report this value in Table 1.

The factor $4$ in equation 2 is a choice we made to ensure that information
has been communicated effectively to all regions and that turbulent motions in
the bulk of the disk reach a statistical steady state (Martizzi et al., 2016).
This choice is validated by comparing the evolution of the velocity dispersion
before and after $t_{r}$, in almost all cases the velocity dispersion has a
fairly gradual (or even flat) development after the relaxation time and a much
more volatile development before (Figure 1). This dynamical timescale results
from complex interplay between the gas density, the shape of the gravitational
potential, and the SFR. For instance, when the gas density increases, the disk
tends to relax more promptly, because as gas cools it comes into a quasi-
equilibrium with the gravitational potential and higher density gas will cool
more efficiently thus reaching steady state sooner. This process is nicely
exemplified in the top and middle panels of Figure 1. In Appendix B we compare
models with the same gas density and gravitational potential but varying SFR,
here we find that $t_{r}$ increases with increasing SNe feedback. In all cases
we find that steady state is reached at $t\approx t_{r}$. In the remainder of
this work we only show results for $t>t_{r}$, when turbulent motions in the
bulk of the disk have reached a statistical steady state.

In both MW progenitor cases the velocity dispersion flattens out at or soon
after the relaxation time indicating that these models are dynamically stable
over long periods of time. The dwarf model however shows a continually
increasing velocity dispersion with time indicating that this model is in fact
not dynamically stable. Over the course of the simulation this model looses
$\approx 35\%$ of its initial mass in SNe driven winds. In comparison the two
MW progenitor models lose $\lesssim 5\%$. This behaviour is consistent with
what is routinely observed when modelling dwarf galaxy systems (e.g. Fragile
et al., 2003; Mac Low & Ferrara, 1999; Fielding et al., 2017). We see similar
behaviour in the MW progenitor with low SFR and high $\epsilon_{*}$ (see
Appendix B), which looses $\approx 10\%$ of its mass to SNe driven winds
during its lifetime.

###  3.2 Metallicity evolution and metal maps

Figure 2: Slice maps through the center of the computational domain for
different galaxy models, indicated by the panel titles. All snapshots are
taken when $\langle[\text{Fe}/\text{H}]\rangle$ = -2.2. The number of SNe,
$N_{\text{SNe}}$, required to reach this mean abundance is 5600, 1300, and
14,000, respectively. The top row shows an edge on view while the bottom row
shows a face on view of the disk. In each panel the left-hand side (blue-green
color map) is number density of H ($\log_{10}n_{\text{H}}$) while the right-
hand side (purple-yellow) is the local gas abundance of [Mg/H]. Each panel is
$6z_{\text{eff}}$ along each axis while the scale bars indicate the relative
sizes between panels. The time since the beginning of the simulation is shown.
These times correspond to $1.89t_{r}$, $3.87t_{r}$, and $1.41t_{r}$,
respectively. In all models, individual SNR are clearly resolved and the
increased homogeneity of the lower SFR models at the same
$\langle[\text{Fe}/\text{H}]\rangle$ can be easily observed.

Over time the material within the disk is progressively enriched as subsequent
SNe increase the mean metallicity of the gas. At the same time SNe drive
turbulent mixing which causes the metals to spread from their localized
injections sites. In general, models with fewer SNe take longer to arrive at
the same mean metallicity and as such their products are given more time to
migrate throughout the disk. This implies that galactic disks with lower SFRs
show more homogeneity at a comparable average abundance. This can be seen in
the evolution of the mean [Fe/H] abundance,
$\langle[\text{Fe}/\text{H}]\rangle$, shown in Figure 1.

$\left[\text{Fe}/\text{H}\right]$ is often used, for convenience, as a rough
proxy for age when thinking about the assembly of metal poor stars (e.g.
Cayrel et al., 2004; Frebel, 2010). This relation is shown in Figure 1 and is,
as expected, non-linear, thus a given galaxy disk patch spends significantly
longer time at higher abundance values than at lower ones. This relationship
is useful to keep in mind as we study the scatter in [Mg/H] of the gas at a
given average [Fe/H] abundance and compare it with that derived from stellar
spectroscopic observations.

In all simulations, an early, chemically unmixed and inhomogeneous galactic
disk is observed. As time evolves, these localized inhomogeneities within the
gas are smoothed out as more events take place and SN products are granted
more extended periods to migrate and mix across the disk. Figure 2 shows slice
maps for each of the three galaxy models considered in this work. Individual
SNR are clearly visible as pockets of lower density (higher pressure) gas with
high metal abundances. The influence of the density and density gradients is
also clearly seen when comparing between different models. The lower density
gas of the MW progenitor with low SFR allows each remnant to reach larger
physical sizes than remnants in the higher density environment of the MW
prognitor with high SFR, leading to a higher degree of homogeneity in the
[Mg/H] abundance of the gas at given $\langle[\text{Fe}/\text{H}]\rangle$ in
the lower SFR models. In Appendix B we show that for identical surface density
of gas (and galaxy potential) a lower rate of SNe also leads to a higher
degree of homogenization of the gas at a given mean metallicity, which is
caused by the slower evolution with time, giving each SNR longer to evolve
before the required metallicity is reached.

The nearly uniform density along the disk plane compel the remnants to stay
nearly circular in this plane, however, the density gradient perpendicular to
the disk causes them to expand more easily along the steep vertical pressure
gradient, which allows metal enriched outflows to reach the edge of the disk
and be ultimately expelled (Martizzi et al., 2015). A review of the time
series of these metal maps shows enriched gas getting expelled from the disk,
rising within the gravitational potential before cooling and returning back,
which is reminiscent of a mini galactic fountain. We note here that some gas
is able to leave the simulation domain and the outflow boundary conditions at
the top and bottom of the computational box does not allow it to return. As a
result, enriched material is steadily lost from the galactic disk patch. As
discussed in Section 3.1 this is a small effect for the MW progenitor models
but a considerable one for the dwarf galaxy model which cannot be considered
dynamically stable over longer periods of time and is unable to effectively
retain the metals injected by SNe. In Section 3.1 and Appendix B we show that
the MW progenitor model with low SFR and high $\epsilon_{*}$ exhibits the same
behaviour, although to a lesser degree than the dwarf-satellite model.

Having discussed the evolution and the salient features of the small patch
disk simulations we now turn our attention to quantifying the degree of metal
mixing and its relevance to spectroscopic observations of stars in the MW.

##  4 Metal Mixing and Comparison to Observations of Metal Poor Stars

In this Section we examine the mixing of metals by SNe in the surrounding gas.
First, we consider the local mixing around the sites where a SN recently
occurred. Second, we study the spread in abundances achieved in each of our
models and compare our simulation results with the spreads inferred from metal
poor halo stars in the MW with the goal to constrain the SFR in the birth
environments of these stars.

###  4.1 Understanding local mixing

In this Section we leverage the high spatial and temporal resolution of the
patch disk simulations to study the metal mixing properties of individual SNe.
In each simulation we record the time and position of all SNe, thus enabling
us to study their evolution and the associated metal dispersion across the
galaxy disk.

####  4.1.1 Metal mixing at the cooling scale

Figure 3: Top panel: Median (solid lines), 25th, and 75th percentile (shaded
region) of the mass of Mg enclosed within $R_{\text{cool}}$ around all SN
injection sites (average over $15\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$
intervals) in the MW progenitor with low SFR (purple) and the one with high
SFR (yellow). The solid black line marks the mass of Mg per injected SN and
thus provides a measurement of the degree of overlap between adjacent mass
injection events. Both models are marginally consistent with a single
enrichment model at low [Fe/H]. Bottom panel: cumulative distribution
functions (CDF) of the mixing masses derived for all injected SNe in the two
MW progenitor models (colors are the same as top panel). Shaded region denotes
the IQR for both distributions. The patchy distribution of elements in the
higher SFR model translates to a wider distribution of mixing masses for the
SNe.

The origin of the dispersion of Mg masses around injection sites depends
sensitively on the degree of overlap between the metals injected from
neighboring SNR. Intuitively, a population of SN is well mixed at a particular
length scale when the dispersion of Mg masses measured across the disk at this
particular scale is modest. Here we select the cooling radius around all of
the injection sites as our characteristic scale. In order to integrate the
metals into new stars, the swept up mass must first cool (Nomoto et al., 2013;
Macias & Ramirez-Ruiz, 2018); therefore the cooling mass yields robust lower
limit on the minimum amount of enriched gas that can be incorporated into the
next generation of stars. The metal content within a cooling mass thus
provides a simple test of the single enrichment star formation assumption
(e.g., Karlsson et al., 2013). It is useful to recall that the cooling radius
(and hence the cooling mass) depends sensitively on the mean density of the
surrounding gas (equation 1) but only weakly on its mean metallicity (Karpov
et al., 2020). The upper panel of Figure 3 shows the enclosed Mg mass within
the cooling scales of the SNe in both MW progenitor models. Both models
display an initial quasi-linear trend in the Mg mass with metallicity followed
by a change of slope at $\langle[\text{Fe}/\text{H}]\rangle\approx-2$. The
early behavior follows the expected evolution for overlapping SNR expanding in
a medium in which the enrichment of the surrounding material is dominated by
recently injected metals and not by the mass of the metals residing in the
ISM. At $\langle[\text{Fe}/\text{H}]\rangle\approx-2$, the mass of Mg
contained within the swept-up ISM begins to be a meaningful contribution to
the total Mg mass within a cooling radius.
In this formalism, the number of overlapping SN within a cooling mass can be
easily inferred by comparing the total enclosed mass of Mg with that expected
from a single SNR, which is illustrated by the solid black line (upper panel
of Figure 3).

As expected, the dispersion of Mg within a cooling mass is smaller for the low
SFR model. When examining the morphology of the individual SNR in Figure 2 it
is clear that this happens because neighboring remnants in this model show
significant overlap. This effect is compounded by the larger cooling radii of
the remnants, caused by the lower mean density of the surrounding ISM. In
contrast, the high SFR model, despite having a larger number of injected SNe
at the same $\langle[\text{Fe}/\text{H}]\rangle$, shows a larger dispersion in
the number of SNR that overlap at the chosen scale. Across the disk we see
regions with no overlapping remnants while other regions show significant
overlap between neighboring remnants. This emanates by virtue of the fact that
individual remnants in this model remain compact enough to be restricted to
scales smaller than the disk scale height. As such, individual remnants sample
a wide range of density scales and thus a corresponding wide range of cooling
masses.

Macias & Ramirez-Ruiz (2018) argued that the abundances of [Mg/H] of stars
with [Fe/H] $\lesssim-2$ could be consistent with enrichment from a single SN
yielding $M_{\text{Mg}}=$0.1\text{\,}\mathrm{M}_{\odot}$$. Our results suggest
that while it is unlikely, at the SFRs considered here, to form stars from gas
that has been enriched by a single event, stars with enrichment by a few
events are plausible at [Fe/H] $\lesssim-2$. This is in broad agreement with
the results of Welsh et al. (2021) who modelled the stochastic enrichment of
metal poor stars by the first SNe and found the expected number of enriching
first SNe to the most metal poor stars to be $N_{\rm SNe}=5_{-3}^{+13}$.

Our simulations suggest that the ISM in the galactic disk patches increased
rapidly in metallicity after less than 100 Myr. With the assumption that Pop
III stars can only form in gas with metallicity below a critical value it is
thus worth quantifying the pristine gas fraction below a threshold
metallicity. Following Pan et al. (2013) we calculate the fraction of gas with
[Fe/H] $<-6$ at the onset of steady state in each of our models and find
$\lesssim 0.5$ % (high SFR MW progenitor at
$t_{r}\approx$20\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$), $\lesssim 1.5$ %
(low SFR MW progenitor at
$t_{r}\approx$40\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$) and 0% (dwarf-like
satellite at $t_{r}\approx$120\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$). The
ubiquity of heavy elements in the Lyman alpha forest indicates that there has
been widespread diffusion from the sites of these early SNe, which is
consistent with our findings. Mackey et al. (2003) estimate that most of the
star-forming material in the intergalactic medium (IGM) has reached [Fe/H]
$<-3.5$ in the range $z\approx 15-20$. Constraining how these Pop III stars
form at high-$z$ is thus of great relevance.

####  4.1.2 Characteristic mixing masses

Enrichment of the ISM depends on how efficiently metals are mixed with the
ambient material swept up by the individual blast waves. At the time of shell
formation (i.e. cooling radius), we find the mass fraction of Mg to be
enhanced in relation to the surrounding medium for most injected remnants,
although that fraction decreases as the simulation evolves. These localized Mg
enhancements are then smoothed out as the enriched gas is given more time to
expand and mix throughout the disk.

Motivated by this, we calculate the scale at which the enclosed [Mg/H] around
each SN site is less than or equal to the mean abundance,
$\langle[\text{Mg}/\text{H}]\rangle$, of the disk. We refer to this integrated
mass as the mixing mass, which broadly represents the characteristic scale at
which SNR become chemically indistinguishable from the background gas. The
bottom panel of Figure 3 shows the cumulative distributions of the mixing
masses in each of the two MW progenitor models. A modest fraction of all SNR
reach their mixing masses within the scale height of the disk, only those that
do are shown. The fraction of the remnants that are confined within $\lesssim
2z_{\text{eff}}$ of the disk is 25% for the MW progenitor model with high SFR
and 20% for the model with low SFR.

The MW progenitor with low SFR shows a narrower range of mixing masses, which
is indicative of a more homogeneous distribution of metals within the disk. By
contrast, the high SFR model exhibits a more patchy distribution of elements
at the same $\langle[\text{Fe}/\text{H}]\rangle$, leading to a larger spread
in mixing masses. The smallest mixing masses within the distribution can be
understood as the masses of those SNe that happened to go off in more pristine
environments at a time when the disk has already achieved a large mean
metallicity, so less dilution is required to reach the background value.
Conversely, the highest mixing masses are caused by SNe that occur in regions
where the local abundance is higher than the disk average background abundance
and are thus harder to dilute.

It is noteworthy that the range of mixing masses we derive from these models
are similar to the range obtained by Macias & Ramirez-Ruiz (2018) for metal
poor stars in the same [Fe/H] range, which have been calculated under the
assumption of mixing with a perfectly pristine background. This assumes that
the cold gas that turns into new stars preserves the variations of [Mg/H] at a
given $\langle[\text{Fe}/\text{H}]\rangle$. Having made the case that the
spread of stellar abundances at a given [Fe/H] can be associated with the star
formation properties of the progenitor disk and their end products, it is
worth thinking further about the implications.

###  4.2 Starforming gas and comparison to observations

Figure 4: Top row: horizontally stacked, mass weighted distributions of [Mg/H]
in the cold, dense gas as a function of
$\langle[\mathrm{Fe}/\mathrm{H}]\rangle$. The columns are the different galaxy
models; panel titles indicate the simulation name. From left to right the
surface density of gas decreases. The lines are the mean (solid) and $1\sigma$
(dashed) of the metal-poor halo stars from Roederer et al. (2014), binned in
0.5 dex wide bins of [Fe/H]. Bottom row: comparison of the spreads of [Mg/H]
in the simulated patches and the metal poor halo stars. Shaded regions are the
$1\sigma$ spread of the distributions in the top row, dashed lines are
$1\sigma$ spread of the halo stars in 0.5 dex wide [Fe/H] bins (same as in top
row).

Motivated by the results presented in the previous Section, we examine the
spread of [Mg/H] as a function of the mean [Fe/H] predicted by our models and
compare it to the spread of stellar abundances seen for metal-poor halo stars.
Our framework does not self-consistently include star formation so, in order
to make a meaningful comparison we consider the cold gas within the disk that
is the most likely to form stars. We impose a temperature threshold of $T\leq
10^{4}{K}$, associated with the cooling function temperature floor, and a
density threshold of $\rho\gtrsim(\rho_{0}/4)$. At each simulation output we
calculate the mass-weighted distributions of [Fe/H] and [Mg/H] in the gas that
fulfils both of these requirements. 222We note that this gives a slightly
different mean abundance than what is considered in the previous Sections
where all the gas, not only the cold component, is used for the calculation.

In Figure 4 we show these distributions of [Mg/H] as as a function of
$\langle[\mathrm{Fe}/\mathrm{H}]\rangle$. The color map indicates the gas mass
($M_{\odot}$) in each bin of [Mg/H] abundance. The differences in [Fe/H] range
covered in each panel are due to the different speeds at which the patches are
enriched by SNe in combination with the different time scales over which the
models come to steady state. All models are initialized at the same [Fe/H] but
as some models relax more quickly the first distribution included may be at
lower metallicity compared to a model which relaxes more slowly (see Section
3.1). The black lines show the mean and $1\sigma$ spread of the [Mg/H]
abundance of the halo stars calculated in 0.5 dex bins in [Fe/H]. The metal
poor halo stars used for this comparison are taken from Roederer et al. (2014)
and the specific sample selection criteria is discussed in Appendix A. We opt
to compare our simulations with this sample because of its relatively large
size ($\approx 300$ metal poor halo stars) and the care taken by the authors
to document the uncertainty of their measurements. This allows us to use the
star to star [Mg/H] scatter as a function of [Fe/H] derived by the authors as
a strict limit for the intrinsic elemental variations.

Some obvious points should be emphasized. The degree of metal mixing of the
cold gas within the disk is, as we argue in the preceding Sections, directly
impressed on the extent of Mg spreading. This metal spreading is thus an
observable diagnostic of the star formation processes occurring within the
disk. In almost a literal sense, the stage is set for us to interpret the
observational properties concerning metal poor stars in the context of
turbulent mixing driven by SN feedback. Additional spread is to be expected
from the intrinsic distribution of Mg and Fe yields in SNe (e.g. Kobayashi et
al., 2006); this would serve to increase the overall spreads observed and
hence the simulations presented here provide us with a lower limit on the
metal dispersion at a given SFR. Thus we are able to constrain the SFR of the
birth environments of the metal poor stars. In addition, we observe large
abundance fluctuations in a few places (Figure 4), suggesting the overall
distribution function for fluctuations has a long, non-Gaussian tail as
expected for interstellar turbulence (Elmegreen & Scalo, 2004; Scalo &
Elmegreen, 2004).

It is evident from Figure 4 that the evolution of the metal variations across
the disk depends fairly strongly on the properties of the system. Even in the
simplest case of a constant Mg yield, we can infer that the properties of
metal poor stars points to stellar birth sites in a MW progenitor with a low
SFR or a satellite galaxy. Our findings give independent evidence that the
dispersion of [Mg/H] observed in metal poor stars are consistent with them
being either deposited into the stellar halo of the MW by mergers of low mass
dwarf galaxies (e.g. Bullock et al., 2001; Robertson et al., 2005; Johnston et
al., 2008) or formed in situ if the MW had a modest SFR as has been recently
suggested by Wang et al. (2021).

##  5 discussion

Motions in a turbulent flow endure over a broad range of length and time
scales. The characteristic length scales correspond to the fluctuating eddy
motions that exist with the largest scales bounded by the geometric dimensions
of the flow. The diffusion coefficient associated with turbulence can be
estimated as (Karlsson et al., 2013; Krumholz & Ting, 2018)

$\kappa\approx{z_{\rm eff}\langle\sigma_{\rm v}\rangle\over 3},$ (3)

where $z_{\rm eff}$ is the outer scale of the turbulence (scale height of the
disk) and $\langle\sigma_{v}\rangle$ is the gas turbulent velocity dispersion
on this scale. Table 1 gives the diffusion coefficient for all the models
presented in this work. For comparison, in the solar neighborhood
$\kappa\approx$300\text{\,}\mathrm{p}\mathrm{c}\;\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$
(Krumholz & Ting, 2018).

An important quality of metal mixing is the connection that exists between the
variance of the spatial metal concentration to the diffusion coefficient.
Krumholz & Ting (2018) used a stochastically forced diffusion framework to
model the spatial metal fluctuations in galaxies and derive a simple estimate
for a characteristic mixing length scale

$x_{\rm m}\approx\left({\kappa\over\Gamma}\right)^{1/4},$ (4)

where $\Gamma$ is the rate of SN events per unit area. Using the values for
$\kappa$ and $\Gamma$ listed in Table 1 along with equation (4), we derive
$x_{\rm m}=$26\text{\,}\mathrm{p}\mathrm{c}$$ and
$67\text{\,}\mathrm{p}\mathrm{c}$ for the high SFR and the low SFR models,
respectively. This simple estimate is in close agreement with the
characteristic mixing mass scales calculated in Section 4.1.2 for the high SFR
and the low SFR models at the fixed
$\langle[\text{Fe}/\text{H}]\rangle\approx-2.3$: $x_{\rm
mix}=$31\text{\,}\mathrm{p}\mathrm{c}$$ and $x_{\rm
mix}=$78\text{\,}\mathrm{p}\mathrm{c}$$, respectively. As we argue in Section
4, the SFR plays an important role in determining the scale of the metal
abundance fluctuations and can be inferred by studying the metallicity
dispersion in stars.333We note here that the typical abundance scatter
predicted by Krumholz & Ting (2018) is smaller than the one deduced from our
simulations and observed in metal poor stars. This is to be expected as their
model is manufactured to predict abundance dispersions for the MW at solar
metallicity.

The mixing mass scale measures a characteristic enlargement of the effective
size of a remnant. At this scale the metallicity fluctuations start to be
smoothed out. An akin pressure length scale (Draine, 2011) can be estimated by
calculating the characteristic radius at which the SN blast wave expansion
velocity becomes comparable to $\langle\sigma_{v}\rangle$:

$x_{P}\approx 67\;\rho_{0,-24}^{-0.37}\left({E_{\rm SN}\over 10^{51}{\rm
erg}}\right)^{0.32}\left(\langle\sigma_{v}\rangle\over 10{\rm
km/s}\right)^{-0.4}{\rm pc},$ (5)

where $\rho_{0}=10^{-24}\rho_{0,-24}\;{\rm g\;cm^{-3}}$. Taking $\rho_{0}$ as
the mid plane density in the galactic disk yields
$x_{P}=$20\text{\,}\mathrm{p}\mathrm{c}$$ and
$x_{P}=$59\text{\,}\mathrm{p}\mathrm{c}$$ for the high SFR and the low SFR
models, respectively. We caution the reader that these length scales do not
take into account galactic winds (e.g. Christensen et al., 2018) and that
$x_{P}$ and $x_{\rm mix}$ are comparable to the scale height of the disk. In
Section 4.1.2 we show that a sizable number of remnants are advected by the
galactic wind before they are able to effectively mix with the ISM, this
implies that the winds simply remove portions of the gas that are not widely
distributed and are highly correlated with the metal field. Therefore, the
effect of metal removal on the dispersion of the metal field is limited.
Similarly, we expect galactic shear (Yang & Krumholz, 2012; Petit et al.,
2015) only to be prominent on relatively longer timescales when compared with
the the timescale for mixing of gas by turbulent feedback from SNe (e.g., Pan
& Scannapieco, 2010; Pan et al., 2013).

Observations of metallicity dispersion in metal poor stars suggest that the
ISM is not well mixed on average at these early times. The observed dispersion
values in stellar abundances of metal poor stars are naturally explained when
considering the early spotty pollution of new elements from SNe because the
numbers of contaminating events are small (Audouze & Silk, 1995). For a
hierarchically structured ISM with star formation operating on a local
dynamical timescale, the size of OB stellar associations is typically
comparable with the metal mixing scale (e.g., Elmegreen, 1998) and the new
elements should continue to be patchy and variable for stars formed in
different clusters.

The statistics of variations of stars in chemical space provide strong
constraints on the SFH of the MW. According to our study, the observed
abundance spread of the metal poor halo stars (Roederer et al., 2014) is most
consistent with the star formation properties of dwarf satellites or an early
low starforming MW like the one described by Wang et al. (2021).

Our findings are in agreement with the idea that the metal poor stellar
content of the MW halo could in principle be explained entirely by disrupted
dwarf satellites (e.g. Bullock et al., 2001; Robertson et al., 2005; Johnston
et al., 2008; Frebel, 2010; Roederer et al., 2018), which has recently gained
stellar kinematic support (Santistevan et al., 2021).

Metal variations in halo stars thus represent an opportunity for chemical
tagging studies to further elucidate their formation origins. In principle
large spectroscopic surveys of MW stars such as Gaia-ESO (Gilmore et al.,
2012) provide us with a unique opportunity for reconstructing the SFH of the
MW by measuring the metal abundance spreads of stars both phase-mixed in the
Galactic field as well as residing in individual stellar streams. While
disentangling the SFHs of the many satellites that contributed to stars in
phase-mixed galactic fields might be challenging, one should be able to infer
their spread in star formation histories before infall. More optimistically,
we hope that for stars residing in individual streams, metal variations may be
key to understand the SFH of the progenitor galaxy. For surviving galaxies, it
might be possible to deduce if there are clear differences between the
formation history of stars belonging to the surviving self-bound satellite and
the debris stars of disrupted accretion event. Our understanding of the
assembly of the MW has come a long way. Nevertheless, it continues to offer
major puzzles and challenges. Metal variations thus provide us with an
exciting opportunity to derive, in combination with kinematic studies and
stellar tagging, the contribution of satellites to the accreted stellar halo
mass, as well as, inferring the specific SFHs of individual galaxies.

We thank S. Faber, E. Kirby, A. Ji, K. Brauer, P. Macias, J. Primack and D.
Zaritsky, for useful discussions. A.N.K. and E.R.-R. acknowledge support by
Heising-Simons Foundation, the Danish National Research Foundation (DNRF132)
and NSF (AST-1911206 and AST-1852393). This work made use of an HPC facility
funded by a grant from VILLUM FONDEN (project number 16599). H.P. acknowledges
support from the Danish National Research Foundation (DNRF132) and the Hong
Kong government (GRF grant HKU27305119, HKU17304821).

## References

  * Abohalima & Frebel (2018) Abohalima, A., & Frebel, A. 2018, ApJS, 238, 36, doi: 10.3847/1538-4365/aadfe9
  * Aoki et al. (2007) Aoki, W., Beers, T. C., Christlieb, N., et al. 2007, ApJ, 655, 492, doi: 10.1086/509817
  * Armillotta et al. (2018) Armillotta, L., Krumholz, M. R., & Fujimoto, Y. 2018, MNRAS, 481, 5000, doi: 10.1093/mnras/sty2625
  * Audouze & Silk (1995) Audouze, J., & Silk, J. 1995, ApJ, 451, L49, doi: 10.1086/309687
  * Beers & Christlieb (2005) Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531, doi: 10.1146/annurev.astro.42.053102.134057
  * Bullock et al. (2001) Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2001, ApJ, 548, 33, doi: 10.1086/318681
  * Cayrel et al. (2004) Cayrel, R., Depagne, E., Spite, M., et al. 2004, A&A, 416, 1117, doi: 10.1051/0004-6361:20034074
  * Christensen et al. (2018) Christensen, C. R., Dave, R., Brooks, A., Quinn, T., & Shen, S. 2018, doi: 10.3847/1538-4357/aae374
  * Colbrook et al. (2017) Colbrook, M. J., Ma, X., Hopkins, P. F., & Squire, J. 2017, MNRAS, 467, 2421, doi: 10.1093/mnras/stx261
  * Dalcanton (2007) Dalcanton, J. J. 2007, ApJ, 658, 941, doi: 10.1086/508913
  * de Avillez & Mac Low (2002) de Avillez, M. A., & Mac Low, M.-M. 2002, ApJ, 581, 1047, doi: 10.1086/344256
  * Draine (2011) Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium
  * Eggen et al. (1962) Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748, doi: 10.1086/147433
  * Elmegreen (1998) Elmegreen, B. G. 1998, in Astronomical Society of the Pacific Conference Series, Vol. 147, Abundance Profiles: Diagnostic Tools for Galaxy History, ed. D. Friedli, M. Edmunds, C. Robert, & L. Drissen, 278.  https://arxiv.org/abs/astro-ph/9712354
  * Elmegreen & Scalo (2004) Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211, doi: 10.1146/annurev.astro.41.011802.094859
  * Few et al. (2014) Few, C. G., Courty, S., Gibson, B. K., Michel-Dansac, L., & Calura, F. 2014, MNRAS, 444, 3845, doi: 10.1093/mnras/stu1709
  * Fielding et al. (2017) Fielding, D., Quataert, E., Martizzi, D., & Faucher-Giguere, C.-A. 2017, 5, 1, doi: 10.1093/mnrasl/slx072
  * Fragile et al. (2003) Fragile, P. C., Murray, S. D., Anninos, P., & Lin, D. N. C. 2003, The Astrophysical Journal, 590, 778–790, doi: 10.1086/375183
  * Frebel (2010) Frebel, A. 2010, Astronomische Nachrichten, 331, 474, doi: 10.1002/asna.201011362
  * Frebel & Norris (2015) Frebel, A., & Norris, J. E. 2015, Annual Review of Astronomy and Astrophysics, 53, 631, doi: 10.1146/annurev-astro-082214-122423
  * Gilmore et al. (2012) Gilmore, G., Randich, S., Asplund, M., et al. 2012, The Messenger, 147, 25
  * Johnston et al. (2008) Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689, 936, doi: 10.1086/592228
  * Karlsson et al. (2013) Karlsson, T., Bromm, V., & Bland-Hawthorn, J. 2013, Reviews of Modern Physics, 85, 809, doi: 10.1103/RevModPhys.85.809
  * Karpov et al. (2020) Karpov, P. I., Martizzi, D., Macias, P., et al. 2020, ApJ, 896, 66, doi: 10.3847/1538-4357/ab8f23
  * Kobayashi et al. (2006) Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, The Astrophysical Journal, 653, 1145, doi: 10.1086/508914
  * Krumholz & Ting (2018) Krumholz, M. R., & Ting, Y.-s. 2018, Monthly Notices of the Royal Astronomical Society, 475, 2236, doi: 10.1093/mnras/stx3286
  * Mac Low & Ferrara (1999) Mac Low, M.-M., & Ferrara, A. 1999, ApJ, 513, 142, doi: 10.1086/306832
  * Macias & Ramirez-Ruiz (2018) Macias, P., & Ramirez-Ruiz, E. 2018, The Astrophysical Journal, 860, 89, doi: 10.3847/1538-4357/aac3e0
  * Mackey et al. (2003) Mackey, J., Bromm, V., & Hernquist, L. 2003, ApJ, 586, 1, doi: 10.1086/367613
  * Martizzi et al. (2015) Martizzi, D., Faucher-Giguère, C. A., & Quataert, E. 2015, Monthly Notices of the Royal Astronomical Society, 450, 504, doi: 10.1093/mnras/stv562
  * Martizzi et al. (2016) Martizzi, D., Fielding, D., Faucher-Giguère, C. A., & Quataert, E. 2016, Monthly Notices of the Royal Astronomical Society, 459, 2311, doi: 10.1093/mnras/stw745
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206, doi: 10.1093/mnras/sty618
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, Annual Review of Astronomy and Astrophysics, 51, 457, doi: 10.1146/annurev-astro-082812-140956
  * Pan & Scannapieco (2010) Pan, L., & Scannapieco, E. 2010, ApJ, 721, 1765, doi: 10.1088/0004-637X/721/2/1765
  * Pan et al. (2013) Pan, L., Scannapieco, E., & Scalo, J. 2013, The Astrophysical Journal, 775, 34, doi: 10.1088/0004-637X/775/2/111
  * Petit et al. (2015) Petit, A. C., Krumholz, M. R., Goldbaum, N. J., & Forbes, J. C. 2015, Monthly Notices of the Royal Astronomical Society, 449, 2588, doi: 10.1093/mnras/stv493
  * Robertson et al. (2005) Robertson, B., Bullock, J. S., Font, A. S., Johnston, K. V., & Hernquist, L. 2005, ApJ, 632, 872, doi: 10.1086/452619
  * Roederer et al. (2018) Roederer, I. U., Hattori, K., & Valluri, M. 2018, AJ, 156, 179, doi: 10.3847/1538-3881/aadd9c
  * Roederer et al. (2014) Roederer, I. U., Preston, G. W., Thompson, I. B., et al. 2014, Astronomical Journal, 147, doi: 10.1088/0004-6256/147/6/136
  * Ryan et al. (1996) Ryan, S. G., Norris, J. E., & Beers, T. C. 1996, ApJ, 471, 254, doi: 10.1086/177967
  * Santistevan et al. (2021) Santistevan, I. B., Wetzel, A., Sanderson, R. E., et al. 2021, MNRAS, 505, 921, doi: 10.1093/mnras/stab1345
  * Scalo & Elmegreen (2004) Scalo, J., & Elmegreen, B. G. 2004, ARA&A, 42, 275, doi: 10.1146/annurev.astro.42.120403.143327
  * Searle & Zinn (1978) Searle, L., & Zinn, R. 1978, ApJ, 225, 357, doi: 10.1086/156499
  * Shen et al. (2015) Shen, S., Cooke, R. J., Ramirez-Ruiz, E., et al. 2015, Astrophysical Journal, 807, 115, doi: 10.1088/0004-637X/807/2/115
  * Shen et al. (2017) Shen, S., Kulkarni, G., Madau, P., & Mayer, L. 2017, MNRAS, 469, 4012, doi: 10.1093/mnras/stx1094
  * Teyssier (2002) Teyssier, R. 2002, Astronomy & Astrophysics, 385, 337
  * van Dokkum et al. (2013) van Dokkum, P. G., Leja, J., Nelson, E. J., et al. 2013, ApJ, 771, L35, doi: 10.1088/2041-8205/771/2/L35
  * Wang et al. (2021) Wang, Y., Nadler, E. O., Mao, Y.-Y., et al. 2021, arXiv e-prints, arXiv:2102.11876.  https://arxiv.org/abs/2102.11876
  * Welsh et al. (2021) Welsh, L., Cooke, R., & Fumagalli, M. 2021, MNRAS, 500, 5214, doi: 10.1093/mnras/staa3342
  * Yang & Krumholz (2012) Yang, C. C., & Krumholz, M. 2012, Astrophysical Journal, 758, 1, doi: 10.1088/0004-637X/758/1/48

##  Appendix A Observational data set

We compare our simulated results with the observations of the 313 metal poor,
$-4.5\leq[\mathrm{Fe}/\mathrm{H}]\leq-1.1$, MW halo stars deduced by Roederer
et al. (2014). After removing stars with upper limit measurements for either
[Fe/H] or [Mg/H], as well as, carbon enhanced metal poor (CEMP) stars, defined
as stars with [C/Fe] $>0.7$ (Aoki et al., 2007), we are left with 224 stars in
the sample. CEMP stars exhibit anomalous abundances of C relative to Fe,
several theories exist to explain this overabundance, one likely model is that
these stars have experienced mass transfer in a binary with an asymptotic
giant branch star (Beers & Christlieb, 2005). If this is the case their
measured abundances will not be representative of the gas from which they
formed and we thus remove them from this sample.

To determine the SNe yield ratio introduced in Section 2 we calculate the mass
ratios of Mg to H and Fe to H and fit a linear relationship. The result of
this is plotted in the left panel of Figure A.1. The best fit has slope:
$\frac{M_{\mathrm{Mg}}}{M_{\mathrm{H}}}=1.27\frac{M_{\mathrm{Fe}}}{M_{\mathrm{H}}}$,
fixing the mass of Mg per event to $0.10M_{\odot}$ and requiring all SNe to
obey this relationship the mass of Fe per event needs to be
$M_{\mathrm{Fe}}=0.08M_{\odot}$.

In the middle panel of Figure A.1 we plot the 224 stars we use from Roederer
et al. (2014), with their associated error estimates. We use the total
uncertainty from their table 12 and the method they outline to calculate
uncertainties for other ratios to get $\sigma_{[\mathrm{Mg}/\mathrm{H}]}$.
Also plotted are the mean and $1\sigma$ spreads we calculate in 0.5 dex wide
bins in [Fe/H]. The shaded region in this Figure corresponds the the ranges we
consider in Section 4.2 (see also Figure 4).

Finally, in the left panel of figure A.1 we show abundances reported to
Jinabase (Abohalima & Frebel, 2018) for MW halo stars marked as non-CEMP,
fulfilling $-3.5\leq[\mathrm{Fe}/\mathrm{H}]\leq-1.0$. When abundances are
reported several times for the same star only one data point is shown for that
star, this selection was made at random. For a subset of these stars with
multiple measurements we calculate the intrinsic star-to-star scatter between
reported abundances from different groups, this is indicated by the colored
error bars on the panel. This intrinsic scatter is indicative of the
systematic errors between different analyses. The existence and size of this
scatter is the reason we decided to compare our results with observations from
a single study to ensure that, to the largest degree possible, the spreads
reported were intrinsic rather than being dominated by systematic
uncertainties between different measurements. The mean [Mg/H] as a function of
mean [Fe/H] of the full sample from Jinabase is quite similar to the subset
from Roederer et al. (2014), the spread is slightly larger given the
systematic errors associated to different analyses. Having said this, our
conclusions are unchanged when considering this sample instead.

Figure A.1: Observations of metal poor halo stars in the MW. The left panel
shows the mass ratio of Mg to H as a function of the same for Fe to H for the
224 stars from Roederer et al. (2014), which we use as our comparison sample.
The solid line is the best linear fit to the data and is used to determine the
mass of Fe per SNe. The middle panel shows the measured abundances of [Mg/H]
and [Fe/H] of the stars, error bars indicate the total uncertainty as reported
by Roederer et al. (2014). The solid line is the mean in 0.5 dex wide bins in
[Fe/H] and the dashed lines are the $1\sigma$ spread in the same bins. These
are the same lines used in Figure 4, the overlapping region between the two
Figures is shown in gray. The right panel shows all abundance measurements of
non-CEMP stars in the MW halo from Jinabase (Abohalima & Frebel, 2018). The
colored crosses show the intrinsic dispersion of measurements for a selection
of stars. Dashed lines are the same as in the middle panel and have been added
to aid comparison.

##  Appendix B Exploring changes in the star formation efficiency

Here we investigate the effects of the star formation (SF) efficiency,
$\epsilon_{*}$, in our simulations. In particular, we explore the mixing
properties of SNe occurring in three versions of the MW model with varying
$\epsilon_{*}$. We use the MW progenitor with low SFR as our base simulation.
We then increase and decrease $\epsilon_{*}$ by a factor of $5$, respectively.
The parameters of these models are summarized in Table 1 in Section 2. In
Figure B.1 we compare the time evolution of $\sigma_{v}$ for the three
different simulations. The influence of $\epsilon_{*}$ on the velocity
dispersion is clearly discernible. When $\epsilon_{*}$ is more or less
extreme, there is a marked difference in the timescales for mixing of gas by
the turbulence driven by SNe. $\epsilon_{*}$ influences the extent to which
turbulence mixes and homogenizes the gas in the face of repeated SN events.
Altering $\epsilon_{*}$ by a factor of 5 leads to a difference of $\approx 2$
in $\sigma_{v}$ and $t_{r}$ (Table 1).

Figure B.1: Comparison of the turbulent velocity dispersion, $\sigma_{v}$, in
the three models using the same galaxy potential with varying rates of SNe per
unit volume. Time is normalized by the relaxation time, $t_{r}$ in all cases.
Both the base and low $\epsilon_{*}$ models reach a steady state and the
velocity dispersion stays relatively constant . The high $\epsilon_{*}$ case,
however, does not reach dynamical stability, similarly to what was observed
for the satellite model in Section 3.1.

The value of $\sigma_{v}$ is set primarily by the intensity of the SN heating,
which depends on the SN rate per unit volume, and the cooling rate of the
shocked gas. The effect of cooling is negligible between models as they have
similar gas density structures. When the energy input from SNe is equal to the
binding energy of the cold gas, like in the model with highest $\epsilon_{*}$
and in the satellite model in Section 3.1, the disk becomes dynamically
unstable and mass blow out occurs. Over its lifetime the high $\epsilon_{*}$
model loses $\approx 10\%$ of its initial mass.

In Section 3.2 we discuss the complexity of element injection and transport
and how our simulations capture the key qualitative feature that metallicity
statistics result from a competition between stochastic injection of SNe,
which produce metal inhomogeneity, and mixing by turbulence, which smooths out
the abundance variations. Notably we show that simulations with higher SN rate
per unit volume produce higher abundance dispersion at the same mean
metallicity. This relationship is even more clear for simulations in which the
galaxy potential and gas surface density are fixed. For example at
$\langle[\text{Fe}/\text{H}]\rangle=-2.5$ the [Fe/H] spread (measured as the
$1\sigma$ of the distribution) is: 0.22, 0.27 and 0.55 for the low
$\epsilon_{*}$, base model and high $\epsilon_{*}$, respectively. In Figure
B.2 we show a face-on view of the galaxy disc in the three models at the same
$\langle[\mathrm{Fe}/\mathrm{H}]\rangle$. As expected, the model with the
highest $\epsilon_{*}$ has a much more patchy distribution of metals at the
same mean abundance. This is because at the same mean abundance (i.e., a
similar amount of SNe) the individual SNR do not have as much time to spread
and mix as they do when the rate is lower. This leads to a more patchy
distribution of metals for the high $\epsilon_{*}$ model.

Figure B.2: Same as bottom row of panels in Figure 2 but for the three MW
progenitor models with varying $\epsilon_{*}$. The slices are taken at the
same mean metallicity as those in Figure 2. Because the total gas mass is the
same in the three models the time to reach the same level of enrichment
depends solely on $\epsilon_{*}$. The effects of rapid enrichment for the high
$\epsilon_{*}$ model are clearly evident in the right panel.

To study the abundance variations with $\epsilon_{*}$ we consider the
distribution of [Mg/H] in the cold, dense gas, which is selected in the same
way as in Section 4. In Figure B.3 we compare the 1d distributions of [Mg/H]
at $\langle[\mathrm{Fe}/\mathrm{H}]\rangle=-2.25$ in the three simulations.
While the mean abundance of [Mg/H] is similar in all three models, the shapes
of the distributions are very different as well as their corresponding
spreads. From Figure B.3 we conclude that the high $\epsilon_{*}$ model is
much less consistent with observations, thus reinforcing our conclusions about
a preference for models with a modest SFR.

Figure B.3: Comparison of the 1d distribution of [Mg/H] in the three models at
mean [Fe/H] = -2.25. Also shown are the $1\sigma$ spread of observations of
[Mg/H] at the same corresponding mean [Fe/H] metallicity (dashed lines).