# The imprint of bursty star formation on alpha-element abundance patterns in Milky Way-like galaxies

Hanna Parul, Jeremy Bailin, Andrew Wetzel, Alexander B. Gurvich, Claude-André Faucher-Giguère, Zachary Hafen, Jonathan Stern, Owain Snaith

► To cite this version:

Hanna Parul, Jeremy Bailin, Andrew Wetzel, Alexander B. Gurvich, Claude-André Faucher-Giguère, et al.. The imprint of bursty star formation on alpha-element abundance patterns in Milky Way-like galaxies. Monthly Notices of the Royal Astronomical Society, 2023, 520, pp.1672-1686. 10.1093/mn-ras/stad206. insu-04848744

HAL Id: insu-04848744

https://insu.hal.science/insu-04848744v1

Submitted on 19 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The imprint of bursty star formation on alpha-element abundance patterns in Milky Way-like galaxies

Hanna Parul,<fnref n="1" />* Jeremy Bailin,<fnref n="1" /> Andrew Wetzel,<fnref n="2" /> Alexander B. Gurvich,<fnref n="3" />

Claude-André Faucher-Giguère,<fnref n="1" />,<fnref n="3" /> Zachary Hafen,<fnref n="1" />,<fnref n="4" /> Jonathan Stern<fnref n="5" /> and Owain Snaith<fnref n="6" />

Accepted 2023 January 17. Received 2022 December 12; in original form 2022 July 7

## ABSTRACT

Milky Way-mass galaxies in the FIRE-2 simulations demonstrate two main modes of star formation. At high redshifts star formation occurs in a series of short and intense bursts, while at low redshifts star formation proceeds at a steady rate with a transition from one mode to another at times ranging from 3 to 7 Gyr ago for different galaxies. We analyse how the mode of star formation affects iron and alpha-element abundance. We find that the early bursty regime imprints a measurable pattern in stellar elemental abundances in the form of a ‘sideways chevron’ shape on the [Fe/H] – [O/Fe] plane and the scatter in [O/Fe] at a given stellar age is higher than when a galaxy is in the steady regime. That suggests that the evolution of [O/Fe] scatter with age provides an estimate of the end of the bursty phase. We investigate the feasibility of observing of this effect by adding mock observational errors to a simulated stellar survey and find that the transition between the bursty and steady phase should be detectable in the Milky Way, although larger observational uncertainties make the transition shallower. We apply our method to observations of the Milky Way from the Second APOKASC Catalogue and estimate that the transition to steady star formation in the Milky Way happened 7 – 8 Gyrs ago, earlier than transition times measured in the simulations.

**Key words:** galaxies: abundances – galaxies: star formation – galaxies: evolution – methods: numerical.

## 1 INTRODUCTION

Star formation in galaxies sometimes has a bursty character, with large fluctuations in the star formation rate (SFR). In contrast to large nearby spiral galaxies with star formation confined to thin rotationally supported gas discs, high-redshift galaxies were gas-rich with irregular morphology dominated by star-forming clumps (Elmegreen et al. 2007). Scatter in the SFR in individual galaxies in the bursty phase is predicted to be ~0.3 dex (Flores Velázquez et al. 2021; Gurvich et al. 2022) comparable to the observed scatter in the population of galaxies on the ‘star formation main sequence’ (SFMS) (Orr et al. 2017; Sparre et al. 2017; Whitaker et al. 2012; Speagle et al. 2014; Leja et al. 2019). One way to test whether star formation is bursty is to study the ratio of SFR indicators sensitive to different time scales. A common choice is to use the Hα-FUV ratio: Hα emission traces star formation over the most recent ~ 5 Myr, while FUV is an indicator of star formation over ~10 – 100 Myr; therefore after a short burst of star formation, the ratio would exceed unity at first and rapidly decrease afterwards. This method was used, for example, in Guo et al. (2016) (with Hβ-FUV) to reveal bursty star formation in low-mass galaxies at 0.4 < z < 1 and in nearby

dwarfs (Weisz et al. 2012; Emami et al. 2019). Unlike dwarf galaxies which continue to be bursty throughout their life, simulations predict that Milky Way-mass galaxies change their star formation mode from highly variable at early times to calm and steady at later times. Sparre et al. (2017) investigated the scatter in the SFR-$M_*$ relation at z = 2 for the FIRE-simulated galaxies and found that due to bursty star formation, the scatter is higher for Hα-derived SFR than for FUV-based SFR and the difference between the two metrics is consistent with observations. They also showed that the scatter increases with decreasing mass. The response of the SFR indicators to realistic star formation histories in FIRE galaxies was also studied in Flores Velázquez et al. (2021). They were aiming to identify the width of the boxcar over which the true SFR should be averaged in order to match the SFR indicated from observations and found that for Hα this best-fitting boxcar-averaging time scale is short (~5 Myr) and relatively stable both in bursty and steady regimes, while for FUV it fluctuates strongly from 10 to 100 Myr during the bursty epoch and stays at the level of 10 – 20 Myr when SFR is steady.

Variable bursty star formation is considered to be one of the main predictions of simulations with explicitly resolved stellar feedback (Faucher-Giguère 2018a). Stinson et al. (2007) used simulations of isolated dwarf galaxies to show that supernova feedback is able to launch periodic episodes of star formation. Domínguez et al. (2015) studied the effect of bursty star formation in dwarf galaxies

<footnote n="1">¹Department of Physics and Astronomy, University of Alabama, Box 870324, Tuscaloosa, AL, 35487, USA</footnote>

<footnote n="2">²Department of Physics & Astronomy, University of California, Davis, 1 Shields Ave, CA 95616, USA</footnote>

<footnote n="3">³Department of Physics & Astronomy and CIERA, Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA</footnote>

<footnote n="4">⁴Department of Physics & Astronomy 4129 Frederick Reines Hall Irvine, University of California, Irvine, CA 92697, USA</footnote>

<footnote n="5">⁵School of Physics & Astronomy, Tel Aviv University, Tel Aviv 69978, Israel</footnote>

<footnote n="6">⁶GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, 92190, Meudon, France</footnote>

<footnote>* E-mail: hparul@crimson.ua.edu</footnote>
from the EAGLE simulations on the ionizing photon production and SFR derived from continuum measurement. Multiple studies (Pontzen & Governato 2012; Governato et al. 2012; Teysier et al. 2013) of simulated dwarf galaxies have focused on cusp-flattening processes of dark matter haloes accompanied by the presence of outflows and galactic fountains induced by bursty star formation. Short-scale variability of star formation is also observed in a diverse sample of galaxies from the FIRE simulations (Hopkins et al. 2018; Muratov et al. 2015; Anglés-Alcázar et al. 2017; Sparre et al. 2017; Ma et al. 2017; Faucher-Giguère 2018b). Of particular interest are the processes that accompany the transition from bursty to steady regime in Milky Way-mass galaxies. Muratov et al. (2015) demonstrated that feedback from bursts at early times is able to create ‘gusty’ superwinds and eject a significant mass of dense ISM which leads to temporary suppression of star formation. Yu et al. (2021) found that the end of the bursty phase correlates with the time when formation of stars with thin-disc kinematics begins to prevail over the formation of stars with thick-disc kinematics; they also found that a more extended steady phase results in a higher thin-disc fraction at z = 0. Gurvich et al. (2022) shows that formation of the thin disc occurs only after the end of the bursty phase. The transition from bursty to steady star formation also correlates with a change in CGM properties from cool (T ~ 10^4 K) and dynamic to hot (T ~ 10^6 K) and nearly hydrostatic (Stern et al. 2021). The connection between thin disc formation and a hot halo may be that the hot halo deposits gas onto the galaxy gently and already-aligned (Hafen et al. 2022). Bellardini et al. (2021) reported that a transition to steady star formation is complemented by smoothing of azimuthal abundance variations and establishment of large-scale radial gradients. Therefore, the transition from bursty to steady mode of star formation marks an important milestone in galactic life, and estimating the transition time based on analysis of observational data can help to test theories of galaxy formation and evolution.

For the Milky Way, the usual method of studying Hα and FUV-derived SFR is not suitable for tracing the evolution of burstiness because it provides only the current state of star formation in the galaxy. We also cannot simply count the number of stars as a function of age to resolve individual bursts since uncertainties in age measurements are comparable to or larger than the typical length of a burst.

One of the widely used ways to look into the past star formation history of the Galaxy is to use stellar elemental abundances. While high mass stars quickly deposit metals into the interstellar medium in form of stellar winds, mass-loss and supernovae explosions, low mass stars show negligible evolution of their surface abundances and preserve the chemical composition of their birthplace environment in their atmospheres. Different elements are produced on different time scales through various formation channels; therefore, the observed patterns of elemental abundance ratios contains the imprint of the star formation history (Tinsley 1979).

In this paper we use Milky Way-mass galaxies from FIRE-2 to show that the scatter of alpha-elements abundances at a given age is correlated with the burstiness of star formation and can be used to recover the mode of historical star formation and estimate the time of transition from bursty to steady regime. We apply our method to the stellar data from the Second APOKASC Catalogue to estimate the end of the bursty period in the Milky Way.

The paper is organized as follows. In Section 2 we describe the set of simulations, in Section 3 we define burstiness of star formation and investigate how it is manifested through age-abundance relations. We show how to estimate the transition time between bursty and steady

modes in Section 3.5 and apply our method to Milky Way data. We discuss our results in Section 4 and conclude in Section 5.

## 2 FIRE-2 SIMULATIONS

We use simulations of Milky Way-like galaxies from the FIRE-2 cosmological baryonic zoom-in simulations described in Hopkins et al. (2018). The simulations were run using the hydrodynamical code GIZMO (Hopkins 2015) in the mesh-free finite-mass Lagrangian Godunov mode. FIRE-2 simulates a number of heating and cooling processes including free-free emission, photoionization/recombination, Compton scattering, photoelectric, metal-line, molecular, fine structure, dust-collisional along with cosmic-ray heating and metallicity-dependent radiative cooling for gas at a range of temperatures T ~ 10 ~ 10^10 K. The simulations model the cosmic UV background based on Faucher-Giguère et al. (2009). Star formation occurs as described in Hopkins et al. (2018) in dense (n > 1000 cm⁻³), self-gravitating, Jeans-unstable molecular gas with a local 100 per cent efficiency per free-fall time. Every star particle is considered as a simple stellar population with Kroupa IMF which evolves along stellar population models from STARBURST99 (Leitherer et al. 1999). Essential for this paper is the detailed simulation of stellar feedback from supernovae of Type Ia (SNIa), core-collapse supernovae (CCSN), and mass-loss dominated by winds from O, B, and AGB stars. The rate of SNe Ia per unit stellar mass per Myr is adopted from Mannucci, Della Valle & Panagia (2006) as follows:

$$ 
\begin{aligned}
dN_{\text{Ia}}/dt &= 0 \text{ for } t_{\text{Myr}} < 37.53 \\
dN_{\text{Ia}}/dt &= 5.3 \times 10^{-8} + 1.6 \times 10^{-5} \exp \left\{ -\left[ \frac{(t_{\text{Myr}} - 50)}{10} \right]^2 \right\} \\
&\quad \text{for } t_{\text{Myr}} \ge 37.53
\end{aligned}
\quad (1) $$

These are normalized such that the expectation number of SNe Ia per stellar mass is 0.001 $M_\odot^{-1}$. For SNe II the rate is approximated as a piecewise-constant function:

$$ 
\begin{aligned}
dN_{\text{II}}/dt &= 0 \text{ for } t_{\text{Myr}} < 3.401 \text{ or } t_{\text{Myr}} > 37.53 \\
dN_{\text{II}}/dt &= 5.408 \times 10^{-4} \text{ for } 3.4 < t_{\text{Myr}} < 10.37 \\
dN_{\text{II}}/dt &= 2.516 \times 10^{-4} \text{ for } 10.37 < t_{\text{Myr}} < 37.53
\end{aligned}
\quad (2) $$

The yields for SN Ia are from Iwamoto et al. (1999) and for SN II from tables in Nomoto et al. (2006) Yields for continuous stellar winds are from the models described in van den Hoek & Groenewegen (1997), Marigo (2001), Izzard et al. (2004), and synthesized in Wiersma et al. (2009); the mass-loss rate is presented in the IMF-integrated form and has a mild dependence on the metallicity. The simulations trace 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), which allows us to constrain the star formation history from the abundance ratio of α-elements and iron. To account for enrichment of gas particles from neighbouring particles through the unresolved metal diffusion (e.g. Rayleigh–Taylor and Kelvin–Helmholtz instabilities and turbulent eddies) the simulations include a model examined in Shen, Wadsley & Stinson (2010):

$$ 
\begin{aligned}
\frac{\partial M_i}{\partial t} + \nabla \cdot (D \nabla M_i) &= 0, \\
D &= C_0 \| \mathbf{S} \|_f \mathbf{h}^2,
\end{aligned}
\quad (3) $$

where **h** is the resolution scale, **S** is the shear tensor, and $C_0$ is a constant calibrated from numerical simulations (Su et al. 2017).

In this paper we used 11 galaxies – 5 galaxies from the Latte suite (Wetzel et al. 2016; Garrison-Kimmel et al. 2017; Hopkins et al. 2018; Garrison-Kimmel et al. 2019a) which are broadly consistent with properties of M31 or the Milky Way and 3 Local Group-like
<table><caption>**Table 1.** Properties of galaxies from FIRE-2 simulations used in this paper.</caption><thead><tr><th>Simulation Name</th><th><i>M</i><sub>*</sub> [M<sub>☉</sub>]</th><th><i>M</i><sub>halo</sub> [M<sub>☉</sub>]</th><th><i>m</i><sub>i</sub> [M<sub>☉</sub>]</th><th><i>t</i><sub>B</sub> [O/Fe] [Gyr]</th><th><i>t</i><sub>B</sub> burst [Gyr]</th><th>Reference</th></tr></thead><tbody><tr><td>m12f</td><td>8.8 × 10<sup>10</sup></td><td>1.4 × 10<sup>12</sup></td><td>7070</td><td>4.71</td><td>4.98</td><td>A</td></tr><tr><td>m12i</td><td>7.0 × 10<sup>10</sup></td><td>9.7 × 10<sup>11</sup></td><td>7070</td><td>3.63</td><td>3.21</td><td>B</td></tr><tr><td>m12m</td><td>1.3 × 10<sup>11</sup></td><td>1.9 × 10<sup>12</sup></td><td>7070</td><td>4.35</td><td>3.79</td><td>C</td></tr><tr><td>m12c</td><td>6.4 × 10<sup>10</sup></td><td>1.9 × 10<sup>12</sup></td><td>7070</td><td>3.26</td><td>3.70</td><td>D</td></tr><tr><td>TheIma</td><td>7.9 × 10<sup>10</sup></td><td>3.4 × 10<sup>12</sup></td><td>4000</td><td>3.26</td><td>2.57</td><td>D</td></tr><tr><td>Louise</td><td>2.9 × 10<sup>10</sup></td><td>1.4 × 10<sup>12</sup></td><td>4000</td><td>5.08</td><td>5.53</td><td>D</td></tr><tr><td>Romeo</td><td>7.4 × 10<sup>10</sup></td><td>1.1 × 10<sup>12</sup></td><td>3500</td><td>6.53</td><td>6.62</td><td>D</td></tr><tr><td>Juliet</td><td>4.2 × 10<sup>10</sup></td><td>9.1 × 10<sup>11</sup></td><td>3500</td><td>4.36</td><td>4.48</td><td>D</td></tr><tr><td>m12b</td><td>9.4 × 10<sup>10</sup></td><td>1.2 × 10<sup>12</sup></td><td>7070</td><td>7.26</td><td>6.32</td><td>D</td></tr><tr><td>Romulus</td><td>1.0 × 10<sup>11</sup></td><td>1.7 × 10<sup>12</sup></td><td>4000</td><td>5.44</td><td>4.88</td><td>E</td></tr><tr><td>Remus</td><td>5.1 × 10<sup>10</sup></td><td>1.1 × 10<sup>12</sup></td><td>4000</td><td>5.44</td><td>5.83</td><td>E</td></tr></tbody></table>

We provide the stellar mass ($M_*$) within the central 20 kpc comoving at $z = 0$, the halo virial mass ($M_{\text{halo}}$) defined by $\rho_{200m}$, the mass resolution for baryonic particles ($m_i$), and transition times derived from analysis of scatter of [O/Fe] as a function of age ($t_{\text{B, [O/Fe]}}$) and from analysis of burstiness ($t_{\text{B, burst}}$) following the definition in Yu et al. (2021). *Notes.* The references are: A: Garrison-Kimmel et al. (2017), B: Wetzel et al. (2016), C: Hopkins et al. (2018), D: Garrison-Kimmel et al. (2019a), E: Garrison-Kimmel et al. (2019b).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Normalized SFR (purple line) overplotted with burstiness, $\sigma_{10}/\text{SFR}$ (blue line). The vertical line corresponds to the time when $\sigma_{10}/\text{SFR} > 500$ drops below 0.2 which separates bursty and steady modes.</figcaption></figure>

pairs from ‘ELVIS on FIRE’ (Garrison-Kimmel et al. 2019a, b). We limited our analysis to *in situ* stars, that is, those that formed within 20 kpc comoving of the host centre. Properties of the simulated galaxies are listed in Table 1.

# 3 ANALYSIS

## 3.1 Star formation histories

We begin with characterizing the mode of star formation as a function of time using the example of galaxy m12i. We present results for all 11 simulations in Appendix A. The purple line on Fig. 1 shows the instantaneous SFR (measured in 1 Myr bins) normalized by the SFR averaged over a 500 Myr rolling window. There are two distinct regimes in which star formation occurs. At early times ($t_{\text{lookback}} \gtrsim 3-4$ Gyr) star formation was ‘bursty’

and characterized by high variability with 10–40 Myr time scale. Recently, over the past 3–4 Gyr, star formation proceeded in a ‘steady’ mode with relatively constant SFR. Fig. 2 gives a closer look at these two modes – in the ‘bursty’ regime, star formation (blue line) on the top panel occurs in bursts with width $\approx 40-50$ Myr and intensity up to 10 times higher than the average value, while the bottom panel shows ‘steady’ regime variations which are much shorter and do not exceed the mean level of SFR by much.

To quantify the variability of star formation and define the time of transition between modes we follow the definition of burstiness in Yu et al. (2021): burstiness is computed as the ratio of the standard deviation of star formation measured in 10 Myr bins to the star formation averaged over 500 Myr, $\sigma_{10}/\text{SFR} > 500$ and is shown with a blue line on Fig. 1. The red vertical line on Fig. 1 marks the end of the bursty phase, $t_{\text{B, burst}}$, which is defined as time when
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution of [O/Fe] for a short interval during the bursty (top) and steady (bottom) modes of star formation. The blue line shows normalized SFR < SFR > $\Sigma_{500}$ for galaxy m12i. Peaks of star formation corresponds to a wider range of [O/Fe].</figcaption></figure>

burstiness falls below 0.2:

$$ \frac{\sigma_{10}(t_{\text{B, burst}})}{\langle SFR \rangle_{500}(t_{\text{B, burst}})} = 0.2. \qquad (4) $$

The burstiness for all galaxies, $\sigma_{10}/\langle SFR \rangle_{500}$ (red line), with vertical lines indicating the end of the bursty regime are shown on panels of Fig. A1 and estimated transition times are reported in Table 1.

## 3.2 Age-abundance relation

The distribution of stellar ages and abundance ratios such as [O/Fe] provides a way to infer historical star formation since stars of different masses inject different elements into the ISM on a different time

scales. The amount of oxygen in the galaxy increases due to CCSN explosions that occur on a relatively short time scales while iron formation is dominated by SNIae that need longer to evolve and explode (in the FIRE-2 implementation SNIae begin when the age of a star particle $t > 37.53$ Myr). Therefore episodes of active star formation are characterized by high [O/Fe] values and the transition to moderate star formation is accompanied by a decrease of [O/Fe] and an increase in iron content.

We now consider the elemental abundance ratio and metallicity for stars of different ages in the galaxy m12i as an example (see Appendix A for details on the other galaxies). Fig. 3 shows stellar metallicity [Fe/H] as a function of lookback time on the left panel and relative oxygen abundance [O/Fe] also as a function of lookback time on the right panel. The iron content in the galaxy increases with time
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Simulation m12i: evolution of [Fe/H] and [O/Fe], our tracer of [α/Fe], with lookback time. The vertical line marks the end of the bursty phase. The distribution of [O/Fe] during the bursty regime (to the right of the vertical line) has a noticeably wider scatter.</figcaption></figure>

with younger stars being more metal-rich while relative oxygen ratio stays nearly constant with minor decrease at low redshift. It should be noted here, that the absolute values of [O/Fe] are larger than those typically observed in the Milky Way and similar galaxies. This difference is likely due to yields and delay time distribution (DTD) of SN Ia implemented in FIRE-2. Gandhi et al. (2022) showed that the values of [Fe/H] increase for the runs with an updated DTD model from Maoz & Graur (2017) due to boost in number of SN Ia events but stay in agreement with observational constraints (see Fig. 8 in Gandhi et al. (2022)). In combination with possible changes in star formation history this could potentially bring down [O/Fe] to near-solar values without significantly affecting [Fe/H].

The most glaring feature of the time evolution of the relative oxygen abundance is the tightening of the distribution at times smaller than ≈ 3 Gyr, to the left of the vertical line which marks the end of the bursty phase (4). The width of the [Fe/H] distribution increases in the steady phase which is related to the development of radial metallicity gradients in the disc (Bellardini et al. 2021).

In Fig. 2 we zoom in on a short interval during the bursty (top panel) and steady (bottom panel) periods and plot both the normalized SFR and [O/Fe]-age distributions. This direct comparison reveals a perfect match between regions with wider distribution on the age-[O/Fe] plane and individual bursts of star formation at earlier lookback times.

The distribution of stars on the [Fe/H]-[O/Fe] diagram also behaves differently in steady and bursty regimes. In general, older populations born during the bursty regime, have a wide range of [O/Fe], while younger stars are concentrated along a narrow sequence. An interesting feature of the bursty mode is that the stars are arranged in a ‘sideways chevron’ shape with two tails spreading into areas of lower and higher [O/Fe]. An example of such a distribution is shown on Fig. 4 for stars born during the burst of star formation which happened 7.6 – 7.64 Gyr ago. Purple lines outlining the ‘chevron’ show the yield tracks for CCSN and AGB stars (top arm) and SNIa stars (lower arm). To compute yield tracks we started with $[Fe/H]_{\text{init}} = -0.55$ and $[O/Fe]_{\text{init}} = 0.48$, corresponding to the iron and relative oxygen abundance of the majority of the stars born during the burst, then we used the rate of SNe Ia, CCSN and mass-loss rate of OB/AGB stars along with ejecta yield masses to obtain the change in oxygen and iron content as a function of time (fits to the stellar evolution models and yields are described in Appendix (A) Hopkins et al. (2018)).

A consistent explanation for the presence of a ‘chevron’ shape in the [Fe/H]-[O/Fe] plane is as follows: an episode of intense star formation leads to rapid enrichment of the ISM with a significant amount of alpha-elements. This is responsible for the upper arm of the chevron, as confirmed by the fact that this arm is aligned with the yield track for CCSN and AGB stars (upper purple line on Fig. 4). However, in the bursty mode, new stars are born mostly in a small number of dense clumps (Oklopčić et al. 2017; Sparre et al. 2017; Faucher-Giguère 2018b) therefore strong stellar feedback following the burst of star formation leads to the generation of outflows, quick depletion of gas, and subsequent cessation of star formation. Gas from the outflows mixes with pristine gas in the halo and falls back down onto the galactic disc providing fuel for a new episode of star formation (Anglés-Alcázar et al. 2017) but with lower metallicity than the metallicity at the end of chevron arm. The lower arm of the chevron is inhabited by stars formed from the gas enriched by SN Ia explosions which are the main source of iron in the ISM; therefore the lower arm on the [Fe/H]-[O/Fe] diagram is aligned with yield vector from SN Ia (Fig. 4). The delay of SN Ia explosions relative to the start of star formation, along with the turbulent and ‘patchy’ (Hayward & Hopkins 2017) character of the ISM during the bursty epoch prevents effective mixing of chemical elements in the ISM during each burst and keeps the chevron arms separated. This is also in agreement with high azimuthal variation of elemental abundances in simulated gas discs at high redshifts explored by Bellardini et al. (2021). In the regime with lower instantaneous SFR, metals produced by each CCSN or SNIa event have enough time to mix with the surrounding gas before new star particles are formed, unlike in the bursty regime when star formation occurs on time scales shorter than the local mixing time. A similar effect was observed in dwarf galaxies in Patel et al. (2021), for the case of m11b galaxy they show diagonal stripes on the [Fe/H]-[Mg/Fe] plane resulting from individual bursts of star formation. Unlike in dwarf galaxies from Patel et al. (2021) we do not see strong age gradient along arms of the chevron because complex temporal and spatial structure of more massive galaxies makes it more difficult to pick out the age trend.

In summary: during the bursty regime, the age-[O/Fe] distribution is wider than during the steady regime, and stars on the [Fe/H]-[O/Fe] plane are organized in a chevron shape. These results suggest that the width of the [O/Fe] distribution measured at a particular time might be a good proxy for the regime in which star formation occurs.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** [Fe/H]–[O/Fe] diagram for the simulation m12i plotted for stars born during a short burst between 7.6 and 7.64 Gyr ago (shown on the subplot). Colour codes the density of data points in each pixel. The 2D distribution has a clear sideways chevron shape: the arms of the chevron follow yield tracks (purple line) calculated for enrichment solely from SNIa (lower arm) and SNII and OB/AGB stars (upper arm). The vertical line indicates the metallicity corresponding to the 75th percentile of stars born during this time interval.</figcaption></figure>

## 3.3 Measuring the scatter of the age-[O/Fe] distribution

To quantify the degree of scatter, referred to as $w_{84}$, we computed the difference between the 8th and 92nd percentiles of [O/Fe] for stars born within a certain time bin. The width of the bin, $\delta t$ depends on the time scale of variability that we want to capture. We tested bin widths of 70, 230, 400, and 900 Myr and found that the overall shape of the curve remains the same for different $\delta t$, although shorter bin width leads to a noisier measurement, especially at early times when each bin has relatively few stars with highly varied abundances. Therefore for the rest of the paper we will use $\delta t = 230$ Myr as it provides sufficient time resolution without information loss.

To increase the contrast between $w_{84}$ values measured in bursty and steady modes we filtered stars based on their metallicity [Fe/H]. Fig. 4 shows the [Fe/H]–[O/Fe] distribution of stars selected during a bursty phase within a narrow range of ages from 7.6 to 7.64 Gyr, with the vertical line indicating the 75th percentile on [Fe/H] for galaxy m12i. The distribution of the metal-rich stars (those to the right of the vertical line) born at this time interval are dominated by the open arms of the chevron, rather than by the overdensity at and to the left of its vertex, and therefore have a wider [O/Fe] distribution.

Fig. 5 shows the effect of the applied metallicity cut on the scatter, $w_{84}$, [O/Fe]. Blue line shows the burstiness, $\sigma_{10}/\text{SFR}_{500}$, yellow line – $w_{84}$, [O/Fe] measured for all stars, and the red line shows the evolution of $w_{84}$, [O/Fe] when only metal-rich stars at given age were considered.

The first thing to notice is that all three curves show the same trend and transition to a plateau at the similar time, around 3 Gyr ago, together with the end of bursty phase of star formation marked by a vertical line. However the measurement for metal-rich stars demonstrates a larger contrast between values measured in bursty and steady modes, providing an easily detectable signal for the end of the bursty regime.

## 3.4 Observational uncertainties

With realistic observational uncertainties, features in age-abundance space can get washed out (e.g. Snaith et al. 2016; Patel et al. 2021) which could potentially lead to a flatter shape of $w_{84}$ and complicate the estimation of the transition time. To estimate the effect of uncertainties we used following simple model.

For age uncertainties, we drew values from a normal distribution with standard deviation $\sigma_{\tau}$, where $\sigma_{\tau}/\tau = 0.1, 0.2, 0.3$, and $\tau$ is the age of the star. The chosen values were motivated by output uncertainties for ages provided in Sanders & Das (2018), where ages, masses, and distances were determined through the Bayesian approach based on the astrometric, photometric, and spectroscopic data. The resulting $w_{84}$ measurements are shown on Fig. 6 (left panel). Larger age uncertainties (blue line) make the transition shallower because of the redistribution of stars formed during the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Burstiness, $\sigma_{SFR}$ (blue curve), and the width of the [O/Fe] distribution, $w_{84}$, as a function of time, measured for stars above the 75th percentile of [Fe/H] at a given age (red curve), versus for all stars (yellow curve). The scatter is in excellent agreement with the burstiness, however $w_{84}$ measured for high-metallicity stars shows a larger difference between values in bursty and steady mode.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of $w_{84}$ for data with different implemented error models for the case of galaxy m121: left: age uncertainties of $\sigma_{\tau}/\tau = \{0.1, 0.2, 0.3\}$, and abundance uncertainty 0.01 dex, right: abundance uncertainty $\sigma_{[O/Fe]} = \{0.01, 0.03, 0.05\}$ dex, and age uncertainties of $\sigma_{\tau}/\tau = 0.1$.</figcaption></figure>

bursts of star formation into adjacent age bins. However even with 50 per cent age errors it is still possible to recover the transition time to the steady regime as a drop in $w_{84}$.

For abundances, we explored uncertainties of 0.01, 0.03, and 0.05 dex. The range of values for uncertainties was picked based on the median uncertainties for stars from the APOKASC-2 Catalogue (Pinsonneault et al. 2018): $\sigma_{[Fe/H], med} = 0.03$ and $\sigma_{[\alpha/Fe], med} = 0.01$. Fig. 6 (right panel) shows the effect of increased abundance error. The presence of larger observational errors increases the level of $w_{84}$, which is especially noticeable in the steady regime (lookback time < 3.5 Gyr) and suppresses the amplitude of variation; however, a rise beyond ≈ 3 – 4 Gyr is still apparent in all cases.

Another possible source of discrepancy between observations and simulation is the size of the sample. Small numbers of stars in a bin can introduce artificial increases to $w_{84}$ due to poor statistics. We tested samples containing 10 per cent and 1 per cent of stars randomly chosen from the initial sample, which resulted in $1.1 \times 10^6$ and $1.1 \times 10^5$ star particles, respectively. We found that reduced sample size does not change the overall shape of the curve, but amplifies the short-scale variations and increases the error of the width measurement.

To explore the effect of the selection function on our measurement, we used mock surveys of synthetic stars generated with the Ananke framework (Sanderson et al. 2020) where each simulation’s star
<figure><img src="image_8.png" /><figcaption>**Figure 7.** *w*<sub>84</sub> as a function of time measured for stars from the simulation (dashed line) and three synthetic surveys (solid lines) generated with Ananke framework (Sanderson et al. 2020). Shaded region shows the variation in measured *w*<sub>84</sub> due to the different survey depth. Both simulations and surveys include error model with σ<sub>τ</sub>/τ = 0.1 and σ<sub>[O/Fe]</sub> = σ<sub>[Fe/H]</sub> = 0.01 dex. The dotted line shows *w*<sub>84</sub> for a subset of evolved stars selected such that their temperature and surface gravity match the properties of stars in the APOKASC-2 Catalogue: 3800 K < *T*<sub>eff</sub> < 5400 K, and 1.0 < log *g* < 3.5. The similarity of all functions means that selection effects in observations should not limit the ability to recover the transition time from bursty to steady star formation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Transition times between the bursty and steady mode of star formation for all 11 galaxies. The transition time measured from burstiness analysis is on the x-axis, and the y-axis contains the transition time from analysis of [O/Fe] scatter with following error model applied: στ/τ = 0.1 and σ[O/Fe] = σ[Fe/H] = 0.01 dex. The dashed line shows the 1:1 relation; all galaxies lie close to the 1:1 line with deviations less than 1 Gyr (shown with shaded region).</figcaption></figure>

particle represented a single stellar population. The survey was
designed to resemble the Gaia DR2 catalogue in data structure,
magnitude limits and observational errors with a self-consistent
dust extinction map, computed directly from the gas metallicity
distribution in the simulation. We used surveys for galaxies m12i,

m12f, and m12m presented in Sanderson et al. (2020) for 3 solar viewpoints resulting in 9 surveys in total. For ages and abundances of stars in these surveys we applied the same simple error model as for star particles from the simulations.

Fig. 7 compares w<sub>84</sub> measured in synthetic surveys (solid lines) to the w<sub>84</sub> in the simulation (dashed line). To the star particles both in simulation and surveys we applied an error model with following uncertainties: σ<sub>τ</sub>/τ = 0.1 and σ<sub>[O/Fe]</sub> = σ<sub>[Fe/H]</sub> = 0.01 dex. Each solid line corresponds to the measurement made for stars from the 3 kpc region centred at three different solar viewpoints which have the same galactocentric distance (R<sub>☉</sub> = 8.2 kpc) but different azimuthal angle. The shaded region shows the variation in w<sub>84</sub> measured for different survey with depths ranging from 3 kpc to 15 kpc. In order to make a proper comparison with the observational data (e.g. APOKASC-2 Catalogue (Pinsonneault et al. 2018)) we restricted our analysis to subsample of synthetic stars with the same temperature and surface gravitiy as the evolved stars in the APOKASC-2: 3800 K < T<sub>eff</sub> < 5400 K and 1.0 < log g < 3.5 (dotted line on the Fig. 7) The widths calculated from the synthetic surveys are in excellent agreement with those measured from the simulation despite the fact that stellar age distribution ‘observed’ in the survey is dominated by young stars. This means that the correlation between the width of the [O/Fe] distribution and burstiness could be in principle observed in the Milky Way regardless of the presence of selection effects.

3.5 Transition times

As we showed in Fig. 5, the transition to the steady regime manifests
itself as a drop in burstiness, and is associated with a similar drop
in w<sub>84</sub>. One interesting feature of the left panel on Fig. 6, which
investigates the effect of age uncertainties on w<sub>84</sub>, is that all three lines
converge at the time close to the end of the bursty phase marked by the
vertical line. That suggests that we can calculate the transition time
based on w<sub>84</sub> by convolving the data with different age uncertainties
and finding the latest time at which they agree.

Fig. 8 summarizes the estimates of the transition times for all
simulations in the sample. Here the x-axis is the transition time
measured from burstiness analysis, while the y-axis contains the
transition time estimated from the width measurement. All galaxies
lie close to the 1:1 line with a deviation of 1 Gyr at most between
t<sub>w84</sub> and t<sub>SFR</sub>, meaning that analysis of the scatter of [O/Fe], w<sub>84</sub>, for
simulated galaxies provides a good estimate of the end of the bursty
epoch.

3.6 Milky Way analysis

In previous sections we showed that our method of inferring the
mode of star formation is robust in the presence of realistic age and
abundance uncertainties, sample size, and probed volume. We now
apply it to observational data to estimate the end of the bursty epoch
for the Milky Way.

The observational data sets we used are the Second APOKASC
Catalogue (Pinsonneault et al. 2018) containing stellar properties
for 6676 evolved stars with spectroscopic data that were obtained
as a part of the APOGEE project (Majewski et al. 2017) and
asteroseismic data from the *Kepler* mission (Borucki et al. 2010).
The catalogue includes evolutionary state, surface gravity, mean
density, mass, radius, along with the age, total metallicity, [M/H],
and alpha-elements abundance ratio, [α/Fe], that is required for our
analysis.

Left panel of Fig. 9 shows that the distribution of APOKASC-2
stars on the age-[O/Fe] plane has a prominent bimodality at earlier
<figure><img src="image_10.png" /><figcaption>**Figure 9.** [α/Fe] as a function of lookback time for stars from Second APOKASC Catalogue (left) and the high-metallicity subsample (right). The high-[α/Fe] sequence, which is present in the full sample, disappears for stars filtered by their metallicity, and therefore, does not dominate the measurement of $w_{84}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** $w_{84}$, [α/Fe] as a function of time for stars from the APOKASC-2 Catalogue. The decline around 8 Gyrs ago marks the end of the bursty epoch of star formation in the Milky Way. Beyond 10 Gyr, the measurement becomes compromised by poor statistics. $w_{84}$, measured for the FIRE-2 galaxies with observational uncertainties close to ones in APOKASC-2 ($\sigma_{\tau/\tau} = 0.1$ and $\sigma_{[O/Fe]} = \sigma_{[Fe/H]} = 0.01$ dex), are overplotted in grey.</figcaption></figure>

times (for stars older than 6 Gyr). Part of this bimodality is explained by the metal-poor thin disc stars that were the subject of discussion in several works (Haywood 2008; Haywood et al. 2013); it is suggested that these stars formed early in the outer parts of the thin disc and migrated to the solar neighbourhood. While they are not related to the bursty star formation, one might be concerned that this bimodality would increase the scatter of [O/Fe] at early times and mimic the signal that we associate with burstiness. However the right panel of Fig. 9 demonstrates that the higher sequence disappears when only metal-rich stars at given age are considered; therefore our measurement of $w_{84}$ measurement would not be dominated by the bimodality.

We measured the width of the [α/Fe] distribution as a function of time for metal-rich (above 75th percentile by [M/H] at a given age bin) stars. Fig. 10 shows the evolution of $w_{84}$, [α/Fe] measured for stars from the APOKASC-2 Catalogue (red line) along with $w_{84}$ for simulated galaxies (grey lines). The error model applied to the simulated galaxies includes age and abundance uncertainties that resemble uncertainties in observational data. For the APOKASC-2 we show the measurement only up to 11 Gyr due to the lack of old

stars however that does not affect our estimate of the transition time. The overall behaviour of $w_{84}$ for APOKASC-2 is similar to $w_{84}$ in the simulations from FIRE-2 – high value of $w_{84}$ at early times followed by a decline and transition to a plateau. The visual inspection of the observational $w_{84}$ suggests that the bursty phase of star formation in the Milky Way ended 7 – 8 Gyr ago.

## 4 DISCUSSION

For the case of simulated Milky Way-like galaxies, we have demonstrated that the scatter of the [O/Fe] distribution measured for metal-rich stars at a given age bin is correlated with the burstiness of historical star formation, which allows an estimation of the transition time from bursty to steady regime. We have also tested whether this result could be achieved with the measurements of the scatter for the [Fe/H] and [O/H] distributions and found that $w_{84}$, [Fe/H] and $w_{84}$, [O/H] do not provide as good proxy for burstiness as $w_{84}$, [O/Fe], especially for the galaxies from Local Group-like pairs. Our findings are in agreement with Carrillo et al. (2022) where they found that for the simulation m12i an increased scatter around the age-[X/Fe] trend at earlier times reflects the bursty mode of star formation.

For the Milky Way, the visual inspection of the evolution of $w_{84}$, [α/Fe] suggests that the bursty phase in our galaxy ended 7 – 8 Gyrs ago. A few recent works explored the connection between the morphology of the FIRE-2 galaxies and the mode of their star formation. Yu et al. (2021) showed that during the bursty mode stars were born with thick-disc kinematics and the formation of the thin disc correlates with the transition to the steady regime; they also estimated the end of the bursty phase for the Milky Way to be ~ 6.5 Gyr ago based on a comparison between the Milky Way and one of the FIRE-2 simulated galaxies, Romeo, which is later than our estimate. Gurvich et al. (2022) analysed the physics of gas disc settling and demonstrated that the rapid development of the rotationally supported gas disc coincides with the end of the bursty phase. Therefore we can expect that the transition time in the Milky Way might match the onset of thin disc formation. For example, Snaith et al. (2014) derived the Milky Way’s star formation history from fitting the age-[Si/Fe] relation with a chemical evolution model and suggested the formation epoch of the thick disc to be 13 – 9 Gyr ago, followed by a ~1 Gyr dip in star formation and then the formation of the thin disc for the last ~8 Gyr, which is in agreement with our estimates. Xiang & Rix (2022) used the data
from Gaia EDR3 and LAMOST to study the stellar age-metallicity distribution $p(\tau$, [Fe/H]) and also found that the distribution splits at ~8 Gyr with the younger part corresponding to the formation of the thin disc; Conroy et al. (2022) explored the chemistry, ages, and kinematics of high-$\alpha$ in situ stars and proposed that the oldest metal-poor in situ stars were born during the 'simmering phase' characterized by low star formation efficiency, which was followed by 'boiling phase' when the efficiency of star formation has increased dramatically. Their estimate of the end of high-$\alpha$ disc formation ($z \approx 1$) also agrees with our estimates for the end of bursty phase. Early (7 – 8 Gyrs) transition to steady star formation and, therefore, early end of the thick disc formation is also in qualitative agreement with the results of several chemical evolution models that predict early and rapid formation of halo and thick disc and more prolonged thin disc formation (Chiappini, Matteucci & Romano 2001; Micali, Matteucci & Romano 2013). A similar picture of the chaotic and isotropic early phase and subsequent rapid disc formation emerges from the analysis of kinematics and abundance data of in situ stars in Belokurov & Kravtsov (2022), who also show the larger scatter in elemental abundances for the low-metallicity (i.e. older) stars and emphasize that the scatter increases when Al, Si, N, and O ratios to either Fe or Mg are considered. To summarize, our estimate of the transition time to the steady phase is in line with other estimates of the end of the thick disc formation that are based on completely different methods.

One of the still remaining questions is whether the chevron shape of the [O/Fe]–[Fe/H] distribution and associated elevated level of scatter of alpha-elements abundances during the bursty phase is related to numerical effects in the simulation. The effect of implemented stellar feedback can be qualitatively understood in a following way: stronger feedback results in stronger bursts and faster depletion of star-forming gas reservoirs (see e.g. fig. 11 in Sparre et al. (2017)). This should lead to a more pronounced upper arm of the chevron and, therefore, increased scatter. For the case of stronger bursts we also expect them to be more isolated and scattered temporally, similar to the bursts in dwarf galaxies (Patel et al. 2021).

Insufficient metal mixing in simulations can artificially increase [O/Fe] scatter. In the FIRE-2 simulations we analysed in this paper, metal mixing is explicitly modelled to account for unresolved (or subgrid) turbulent eddies and metal diffusion between boundaries of neighbouring resolution elements. Hopkins et al. (2018) showed that this numerical diffusion term has negligible effect on star formation histories; however, it is important for proper modelling of chemical evolution and abundance ratios. This was investigated in Escala et al. (2018) for dwarf galaxies in FIRE-2. They found that including turbulent metal diffusion tends to reduce the intrinsic scatter in [α/Fe] versus [Fe/H] and allows agreement between the width of the metal distribution function and dispersion of α-element abundance ratios of the simulated and observed galaxies. Stronger mixing would lead to a further reduction of the chevron and tightening of the [O/Fe] – age relation. A different implementation of metal mixing in the simulations might therefore affect the magnitude of the effect we have identified, but is unlikely to change its qualitative existence.

Escala et al. (2018) compared intrinsic scatter in [α/Fe] versus [Fe/H] and its evolution between simulated and real dwarf galaxies. Due to the correlation of [Fe/H] with age they conclude that intrinsic scatter of [α/Fe] at fixed [Fe/H] might serve as a characteristic of ISM homogeneity. For the galaxies in their sample they found that evolution of [Si/Fe] scatter is characterized by presence of individual peaks resulting from starbursts which is in agreement with our finding for more massive galaxies. Outside of the peaks, their reported scatter stays near zero through all the time which is interpreted as a sign of

well-mixed ISM at almost all times. In contrast to Escala et al. (2018) our measurement of scatter is nowhere near zero during the bursty phase. Part of the explanation is that, unlike in dwarfs, in Milky Way-like galaxies bursts could affect the abundance distribution only in parts of the galaxy. Also, bursts in massive galaxies are more frequent than in dwarfs, therefore, in our analysis each considered time bin (230 Myr) covers a few burst episodes. In other words, instead of focusing on separate starbursts we study the evolution of the star formation regime on a longer time scale. The assumption of a well-mixed ISM during bursty periods needs some clarification: galactic discs during bursty epochs indeed lack the radial gradients that exist in present-day discs (Bellardini et al. 2021); however, formation of two distinct arms on the [Fe/H]–[O/Fe] plane requires the ISM to have clumps which might have different elemental abundance ratios than that of the smooth ISM around them, and therefore, abundances in the disc experience variations on smaller scales. That is confirmed by the work of Bellardini et al. (2021), who showed that azimuthal scatter of abundances in gas discs increase at higher redshift.

In our analysis we included only stars formed within the central 20 kpc of the host, therefore we mostly ignored the effect of mergers. Accreted stars could be filtered out based on the analysis of their chemodynamics so they would not affect the computation of $w_{84}$. However, mergers could induce strong bursts of star formation (Sparre & Springel 2016) and lead to a potential mismatch between burstiness and scatter of [O/Fe] as in case of m12b or elevated scatter in the steady regime as in the case of the recent merger in m12f. Therefore, a merger event occuring near the end of the bursty phase might affect the estimation of the transition time. We tested the effect of including *ex situ* stars by computing transition times for all stars within the central 100 kpc of the host at $z=0$ and found the effect negligible for almost all galaxies except for m12c where the difference was around 260 Myr, which is less than the typical uncertainty of our method.

The scatter of the alpha abundances in the Milky Way shows the same behaviour as in the FIRE-2 simulated galaxies, for which the scatter and star formation burstiness are correlated. We have suggested that we can therefore interpret the abundances of Milky Way stars to infer the history of its star formation character. However, it would be interesting to explore evolution of alpha abundance distribution in galaxies with wider range of star formation histories, including those simulations that do not show early bursty period for example, Park et al. (2021).

## 5 CONCLUSION

For 11 Milky Way-mass galaxies from the FIRE-2 simulations we analysed the distribution of [O/Fe] as a function of time for high-metallicity stars formed in a host galaxy. We found that the width of the alpha element abundance ratio distribution, $w_{84}$, is correlated with burstiness of star formation, $\sigma_{10}/<SFR_{500}>$ and this correlation is robust to observationally motivated age and abundance uncertainties and selection function. The wider distribution of alpha-elements during bursty periods is related to the behaviour of the stellar distribution on the [Fe/H]–[O/Fe] plane – an intense burst of star formation followed by depletion of the star-forming gas reservoir leads to the large amount of CCSN enriching the ISM mainly with alpha-elements, while iron mostly comes from the SNIae. For sufficiently intense bursts, this happens faster than the mixing time scale, resulting in the formation of a ‘sideways chevron’ shape on the [Fe/H]–[O/Fe] plane. Measuring the scatter of [O/Fe] for stars above 75th percentile in [Fe/H] at a given time captures the presence of chevron arms. In the steady regime of star formation, mixing
between different enrichment sites occurs more effectively and leads to a more uniform distribution of elements in the disc at a given radius. We estimate the end of the bursty phase as the time when $w_{84}$ as measured from data convolved with different age uncertainties converge. We find that this transition time is in excellent agreement with that derived directly from the star formation history. We applied this to the Milky Way and computed $w_{84}$ for data from the Second APOKASC Catalogue to estimate the end of the bursty star formation mode. We found that $w_{84}$ computed from the observational data set shows similar behaviour, and implies that the Milky Way transitioned from its bursty phase to its steady phase ≈ 7 – 8 Gyr ago, earlier than simulations.

## ACKNOWLEDGEMENTS

We thank the reviewer for providing us with constructive comments. We are also grateful for useful discussions with Jamie Tayar, Sanjib Sharma, and Eliot Quataert.

Simulations used in this paper were run on the Caltech computer cluster ‘Wheeler,’ allocations AST21010 and AST20016 supported by the NSF and TACC. JB was supported by HST grant number AR-12837. AW received support from: NSF via CAREER award AST-2045928 and grant number AST-2107772; NASA ATP grant number 80NSSC20K0513; HST grants AR-15809, GO-15902, and GO-16273 from STScI. CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant number 17-ATP17-0067; and by a Cottrell Scholar Award from the Research Corporation for Science Advancement.

## DATA AVAILABILITY

The FIRE-2 simulations are publicly available (Wetzel et al. 2022) at http://flathub flatironinstitute.org/fire. Additional FIRE simulation data are available at https://fire.northwestern.edu/data. A public version of the GIZMO code is available at http://www.tapir.calt.ech.edu/~phopkins/Site/GIZMO.html.

## REFERENCES

Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2017, MNRAS, 470, 4698

Bellardini M. A., Wetzel A., Loebman S. R., Faucher-Giguère C.-A., Ma X., Feldmann R., 2021, MNRAS, 505, 4586

Belokurov V., Kravtsov A., 2022, MNRAS, 514, 689

Borucki W. J. et al., 2010, Science, 327, 977

Carrillo A., Ness M. K., Hawkins K., Sanderson R., Wang K., Wetzel A., Bellardini M. A., 2023, ApJ, 942, 35

Chiappini C., Matteucci F., Romano D., 2001, ApJ, 554, L1044

Conroy C. et al., 2022, preprint (arXiv:2204.02989)

Domínguez A., Siana B., Brooks A. M., Christensen C. R., Bruzual G., Stark D. P., Alavi A., 2015, MNRAS, 451, 839

Elmegreen D. M., Elmegreen B. G., Ravindranath S., Coe D. A., 2007, ApJ, 658, L763

Emami N., Siana B., Weisz D. R., Johnson B. D., Ma X., El-Badry K., 2019, ApJ, 881, L71

Escala I. et al., 2018, MNRAS, 474, 2194

Faucher-Giguère C.-A., 2018a, Nat. Astron., 2, 368

Faucher-Giguère C.-A., 2018b, MNRAS, 473, 3717

Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, L1416

Flores Velázquez J. A. et al., 2021, MNRAS, 501, 4812

Gandhi P. J., Wetzel A., Hopkins P. F., Shappee B. J., Wheeler C., Faucher-Giguère C.-A., 2022, MNRAS, 516, 1941

Garrison-Kimmel S. et al., 2017, MNRAS, 471, 1709

Garrison-Kimmel S. et al., 2019a, MNRAS, 487, 1380

Garrison-Kimmel S. et al., 2019b, MNRAS, 489, 4574

Governato F. et al., 2012, MNRAS, 422, 1231

Guo Y. et al., 2015, ApJ, 800, L39

Guo Y. et al., 2016, ApJ, 833, L37

Gurvich A. B. et al., 2023, MNRAS, 519, 2598

Hafen Z. et al., 2022, MNRAS, 514, 5056

Hayward C. C., Hopkins P. F., 2017, MNRAS, 465, 1682

Haywood M., 2008, MNRAS, 388, 1175

Haywood M., Di Matteo P., Lehnert M. D., Katz D., Gómez A., 2013, A&A, 560, 109

Hopkins P. F., 2015, MNRAS, 450, 53

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439

Izzard R. G., Tout C. A., Karakas A. I., Pols O. R., 2004, MNRAS, 350, 407

Leitherer C. et al., 1999, ApJS, 123, 3

Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, L3

Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, MNRAS, 467, 2430

Majewski S. R. et al., 2017, AJ, 154, 94

Mannucci F., Della Valle M., Panagia N., 2006, MNRAS, 370, 773

Maoz D., Graur O., 2017, ApJ, 848, L25

Marigo P., 2001, A&A, 370, 194

Micali A., Matteucci F., Romano D., 2013, MNRAS, 436, 1648

Moreno J. et al., 2021, MNRAS, 503, 3113

Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691

Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nucl. Phys. A, 777, 424

Oklopćić A., Hopkins P. F., Feldmann R., Kereš D., Faucher-Giguère C.-A., Murray N., 2017, MNRAS, 465, 952

Orr M. E. et al., 2017, ApJ, 849, L2

Park M. J. et al., 2021, ApJS, 254, 2

Patel P. B., Loebman S. R., Wetzel A., Faucher-Giguère C.-A., El-Badry K., Bailin J., 2022, MNRAS, 512, 5671

Pinsonneault M. H. et al., 2018, ApJS, 239, 32

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Sanders J. L., Das P., 2018, MNRAS, 481, 4093

Sanderson R. E. et al., 2020, ApJS, 246, 6

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Snaith O. N., Haywood M., Di Matteo P., Lehnert M. D., Combes F., Katz D., Gómez A., 2014, ApJ, 781, L31

Snaith O. N., Bailin J., Gibson B. K., Bell E. F., Stinson G., Valluri M., Wadsley J., Couchman H., 2016, MNRAS, 456, 3119

Sparre M., Springel V., 2016, MNRAS, 462, 2418

Sparre M., Hayward C. C., Feldmann R., Faucher-Giguère C.-A., Muratov A. L., Kereš D., Hopkins P. F., 2017, MNRAS, 466, 88

Speagle J. S., Steinhardt C. L., Capak P. L., Silverman J. D., 2014, ApJS, 214, 15

Stern J. et al., 2021, ApJ, 911, L88

Stinson G. S., Dalcanton J. J., Quinn T., Kaufmann T., Wadsley J., 2007, ApJ, 667, L170

Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, MNRAS, 471, 144

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Tinsley B. M., 1979, ApJ, 229, L1046

van den Hoek L. B., Groenewegen M. A. T., 1997, A&AS, 123, 305

Weisz D. R. et al., 2012, ApJ, 744, L44

Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23

Wetzel A. et al., 2022, preprint (arXiv:2202.06969)

Whitaker K. E., van Dokkum P. G., Brammer G., Franx M., 2012, ApJ, 754, L29
Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574

Xiang M., Rix H.-W., 2022, Nature, 603, 599

Yu S. et al., 2021, MNRAS, 505, 889

## APPENDIX A: BURSTINESS AND ABUNDANCE PATTERNS FOR ALL SIMULATIONS

Each panel of Fig. A1 shows the evolution of the instantaneous SFR, measured in 1 Myr bins and normalized by the SFR averaged over 500 Myr rolling window (purple line). At early times ($t_{\text{lookback}} \gtrsim 5-6$ Gyr) star formation was ‘bursty’ and characterized by high variability with 10–40 Myr time scale. Recently, over the past 4–6 Gyr, star formation proceeded in a ‘steady’ mode with relatively constant SFR with the exception of a few isolated bursts like in m12f which are related to minor mergers. The red lines show burstiness, $\sigma_{10}/< \text{SFR} >_{500}$ and vertical lines mark transition times, defined as the time when burstiness falls below 0.2.

Fig. A2 shows iron abundance [Fe/H] and Fig. A3 shows relative oxygen abundance [O/Fe] as a function of lookback time for

all galaxies in the sample with vertical lines marking transition times derived from analysis of star formation history. Each galaxy demonstrates two distinct modes of enrichment with oxygen and iron. The scatter of the [O/Fe] distribution rapidly decrease after the end of the bursty phase.

Fig. A4 shows burstiness, $\sigma_{10}/< \text{SFR} >_{500}$ (blue curve), and the width of the [O/Fe] distribution, $w_{84}$, as a function of time, measured for stars above the 75th percentile of [Fe/H] at a given age (red curve), versus for all stars (yellow curve). For all galaxies the scatter is well correlated with burstiness, however, $w_{84}$ measured for metal-rich stars (red line) shows higher contrast between values measured in bursty and steady modes.

Fig. A5 shows the [O/Fe] scatter, $w_{84}$ versus time, computed using three different age ($\tau$) error models ($\sigma_{\tau}/\tau = 0.1, 0.2, 0.3$), and abundance uncertainty of 0.01 dex. The time of convergence of three $w_{84}$ functions (vertical-black line) is a close match to the transition time to the steady regime marked by the red-vertical line. Burstiness is overplotted in grey to emphasize the agreement between the scatter and mode of star formation.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** Star formation histories measured in 1 Myr bins using the ages of stars at $z=0$ (purple line) of all galaxies in our sample. All galaxies at high redshifts show a bursty mode of star formation, characterized by rapid variations with high amplitude, which is followed by steady star formation at more recent times (past ~5 Gyr). The red line shows burstiness, $\sigma_{10}/< \text{SFR} >_{500}$ and vertical lines mark transition times, defined as the time when burstiness falls below 0.2.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure A2.** Evolution of [Fe/H] with lookback time. The vertical line marks the end of the bursty phase as on the corresponding panel of Fig. A1.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure A3.** Age-[O/Fe] distribution for all 11 galaxies plotted for stars with metallicity higher than the 75th percentile value at a given time bin.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A4.** Burstiness, $\sigma_{10}/<\text{SFR}>_{500}$ (blue curve), and the width of the [O/Fe] distribution, $w_{84}$, as a function of time, measured for stars above the 75th percentile of [Fe/H] at a given age (red curve), versus for all stars (yellow curve). The scatter is well correlated with burstiness, however, $w_{84}$ measured for metal-rich stars (red line) shows higher contrast between values measured in bursty and steady modes.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A5.** Each panel shows the [O/Fe] scatter, $w_{84}$ versus time for a different simulated galaxy, computed using three different age error models and abundance uncertainty of 0.01 dex. The time of convergence of three $w_{84}$ functions (vertical-black line) is a close match to the transition time to the steady regime marked by the red-vertical line. Burstiness is overplotted in grey to emphasize the agreement between the scatter and mode of star formation.</figcaption></figure>
## APPENDIX B: CASE OF m12b

Despite a general good agreement between the burstiness and the scatter of [O/Fe], few galaxies show some mismatches between the $\sigma_{10}/< SFR >_{500}$. One such striking example is the simulation m12b, which have the peaks of burstiness and $\sigma_{[O/Fe]}$ separated by ~2 Gyr at high redshifts. To investigate the source of this discrepancy we compared the stellar distribution at snapshots 200 and 240 corresponding to lookback times 9.7 Gyr and 8.7 Gyr (Fig. B1). Both snapshots capture the galaxy in a bursty stage, so $\sigma_{SFR}$ and

<figure><img src="image_17.png" /><figcaption>**Figure B1.** Burstiness and [O/Fe] scatter for galaxy m12b with vertical lines indicating the times of snapshots 200 and 240. The decrease of the burstiness has a more monotonical behaviour, while $w_{84}$ shows a very high-amplitude peak at 9 Gyr which does not seem to match any particularly strong individual burst episode.</figcaption></figure>

$w_{84}$ are both relatively high; however, while burstiness declines on the interval bounded by these snapshots, the scatter of [O/Fe] rises dramatically. Both snapshots 200 and 240 have a chevron shape, clearly seen on the [Fe/H]–[O/Fe] diagram colour-coded by the age of the stars (top left panel on Fig. B2); however, the distribution for the later snapshot has a striking feature – a very long upper arm corresponding to a short and intense burst of star formation which occured in a very compact (less than 0.5 kpc) region near the centre.

This example supports the idea that the chevron shape and associated increased scatter in [O/Fe] as a function of age are related to episodes when intense star formation is confined to a number of dense clumps. The clumpy character of star formation in high-redshift galaxies has been confirmed in observations (e.g. Guo et al. 2015) and in FIRE simulations (Oklopčić et al. 2017). In m12b, the stellar distribution at snapshot 200 contains multiple concentrated clumps with very recent star formation. Snapshot 240, which is 1 Gyr later, has an overall smoother stellar distribution, punctuated by a small but intense centre of star formation, which rapidly produces new stars and deposits a large amount of alpha-elements into the ISM on time scales shorter than delay of SNIae. This is the source of the arm on the [Fe/H]–[O/Fe] diagram. m12b had a merger with a satellite in between snapshots 200 and 240; this merger might be a cause of the central burst of star formation given that galaxy interactions lead to fuel-driven enhancement of star formation in the central kiloparsec, as discussed in Moreno et al. (2021).

To sum up, the case of m12b with the decoupling between $w_{84}$ and $\sigma_{SFR}$ demonstrates that dispersion of [O/Fe], $w_{84}$, is very sensitive to spatially concentrated intense bursts of star formation even in cases where those bursts are not dominating the overall star formation mode of the galaxy.
<figure><img src="image_18.png" /><figcaption>**Figure B2.** Top row: [Fe/H]–[O/Fe] diagrams for galaxy m12b at snapshots 200 (**left**) and 240 (**right**) coloured by the age of the stars at the time of the snapshot, middle row: same as the top row but coloured by cylindrical radius, and bottom row: distribution of stars on the XY plane. The earlier snapshot shows the presence of multiple clumps of ongoing star formation, while the galaxy at the later snapshot has all its star formation confined in the compact central region which gives a rise to the extended upper 'arm' on the corresponding [Fe/H]–[O/Fe] diagram.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.