Cosmological evolution of metallicity correlation functions from the Auriga simulations Open Access

ABSTRACT

We study the cosmological evolution of the two-point correlation functions of galactic gas-phase metal distributions using the 28 simulated galaxies from the Auriga Project. Using mock observations of the z = 0 snapshots to mimic our past work, we show that the correlation functions of the simulated mock observations are well matched to the correlation functions measured from local galaxy surveys. This comparison suggests that the simulations capture the processes important for determining metal correlation lengths, the key parameter in metallicity correlation functions. We investigate the evolution of metallicity correlations over cosmic time using the true simulation data, showing that individual galaxies undergo no significant systematic evolution in their metal correlation functions from z ∼ 3 to today. In addition, the fluctuations in metal correlation length are correlated with but lag ahead fluctuations in star formation rate. This suggests that re-arrangement of metals within galaxies occurs at a higher cadence than star formation activity, and is more sensitive to the changes of environment, such as galaxy mergers, gas inflows/outflows, and fly-bys.

1 INTRODUCTION

Understanding the baryonic processes responsible for metal (elements heavier than hydrogen and helium) transportation is crucial in galaxy evolution. Metals can be found in stars, their birth places where they are synthesized through nucleosynthesis; they can also be found in interstellar medium (ISM), into which they are injected by stars through various mechanisms including supernova explosions and stellar winds. Once in the ISM, metals mix with the surrounding gas and diffuse, leading to changes in the overall metallicity distribution within galaxies, participating in next-generation star formation, and thus forming a crucial part of the baryon cycle (Tinsley 1980).

Both observers and theorists have studied this process. On the observational side, gas-phase metallicity is most commonly measured by the abundance of oxygen, the most abundant metal in the Universe. The invention of integrated field units (IFUs, e.g. Croom et al. 2012; Sánchez et al. 2012; Bundy et al. 2015; Erroz-Ferrer et al. 2019) has made it possible to measure the spatial distribution of oxygen abundance in nearby galaxies (e.g. Rosales-Ortega et al. 2011; Sánchez-Menguiano et al. 2016). IFU studies have revealed that metallicity gradients are ubiquitous and that their steepness is correlated with other galaxy properties such as stellar mass (e.g. Belfiore et al. 2017; Ho et al. 2018; Poetrodjojo et al. 2018; Sánchez-Menguiano et al. 2018; Kreckel et al. 2019). On the theoretical side, numerical models have begun to explore the origin and evolution of metallicity gradients (e.g. Di Matteo et al. 2009; Ma et al. 2017; Sharda et al. 2021; Tissera et al. 2022), but have also gone beyond simply 1D statistics such as gradients to examine chemical mixing in a broader view. The models contain various mechanisms, including bar driven mixing (Di Matteo et al. 2013), spiral-arm driven mixing (Grand et al. 2016; Orr et al. 2023), thermal instabilities (Yang & Krumholz 2012), gravitational instabilities (Petit et al. 2015), cosmological accretion (Ceverino et al. 2016), and supernova-driven turbulence (de Avillez & Mac Low 2002; Colbrook et al. 2017).

Exploring beyond simple metallicity gradients requires new statistical tools to characterise 2D metallicity maps. One of the simplest and most straightforward to measure from observations is the two-point correlation function. Two-point correlations are able to provide unique information of ISM turbulence and chemical mixing by decoding metal fields of galaxies. Krumholz & Ting (2018, hereafter KT18) propose a model to predict this quantity that considers metal injection and diffusion. This prediction has inspired a number of observational studies (Kreckel et al. 2020; Metha, Trenti & Chu 2021; Williams et al. 2022) that measure the two-point correlations (or closely related statistics) of PHANGS-MUSE galaxies. While much of this effort has focused on very nearby galaxies observed at extremely high resolution, more recently, Li et al. (2021, 2023) extended the method to larger but more distant galaxy samples in CALIFA (Sánchez et al. 2012) and AMUSING++ (López-Cobá et al. 2020), respectively, and demonstrated its robustness against beam smearing, choice of metallicity diagnostic, and binning used to remove the overall metallicity gradient. The samples in these studies are large enough to permit for the first time an examination of statistical distributions of metallicity correlation functions across the local galaxy population, and the correlations between metal distributions and other galactic properties. Most strikingly, these studies have uncovered a robust correlation between the galaxies’ metallicity correlation lengths – the characteristic size scales of their metallicity correlations after removing the large-scale gradient – and their other bulk properties such as stellar mass, size, and star formation rate (SFR). This discovery raises immediate questions: How and when in the course of galaxy evolution were these relationships established? Have they varied over cosmic time?

While there has been extensive work on the cosmological evolution of metallicity gradients using cosmological zoom-in simulations (e.g. Di Matteo et al. 2009; Torrey et al. 2012; Ma et al. 2017; Bellardini et al. 2021; Metha & Trenti 2023), there has yet to be a similar study of higher order metallicity statistics such as the metallicity correlation function, despite the growing body of observational literature. In this paper, we present a pioneering study aimed at filling this gap. The correlation function is of interest because it encodes the fundamental physics of chemical mixing in the ISM, the process responsible for distributing metals from their birth sites. This process, in turn, depends on the evolution of ISM turbulence across cosmic time, another process of great interest. Our specific aims are to (1) determine if the simulations are able to reproduce the distribution of two-point correlation functions found in observed local galaxies, (2) use the simulations to study the cosmological evolution of two-point correlations so that we can understand when and how they are established, and (3) guide future observational work aimed at measuring metallicity correlations beyond the local Universe.

This paper is organized as follows. In Section 2, we describe the Auriga simulations and the procedures by which we extract metallicity maps from them. In Section 3, we discuss the pipeline of mock observations we use to test whether the simulations are consistent with local galaxy observations. In Section 4, we describe our findings regarding the cosmological evolution and origin of metallicity correlations. Finally, we discuss and summarize our work in Section 5.

2 SIMULATIONS

2.1 Description

In this work, we examine 28 simulated galaxies from the Auriga Project (Grand et al. 2017, named as AuN, where N is the halo number). The Auriga simulations are high-resolution cosmological zoom-in simulations using the magnetohydrodynamic (MHD) code arepo (Springel 2010). arepo is a quasi-Lagrangian, moving-mesh code that tracks the evolution of MHD cells and collisionless particles in a Lambda cold dark matter cosmological setting, with Planck Collaboration XVI (2014) cosmological parameters Ω_m = 0.307, Ω_b = 0.048, Ω_Λ = 0.693, and a Hubble constant of H₀ = 100h km s⁻¹ Mpc⁻¹, where h = 0.6777. The simulations include primordial and metal-line cooling, a uniform ultraviolet (UV) background field for reionisation, star formation, magnetic fields, active galactic nuclei, energetic and chemical feedback from Type II supernovae, and mass-loss/metal return due to Type Ia supernovae and asymptotic giant branch (AGB) stars, accounting for nine elements (H, He, C, N, O, Ne, Mg, Si, Fe; Grand et al. 2017).

Star formation in the Auriga simulations occurs in cells that exceed a density threshold of ρ₀ = 0.13 atoms cm⁻³; cells exceeding this threshold form stars with a star formation time-scale of τ_SF = 2.2 Gyr (following Springel & Hernquist 2003). The Auriga simulations assume that every star cell represents a simple stellar population (SSP) with a specified age and metallicity that is directly inherited from the surrounding gas. The distribution of stellar masses present in each SSP follows a Chabrier (2003) initial mass function.

The Auriga simulation suite includes runs from four different particle resolutions, and in this work we discuss two of them. Six haloes (Au6, Au16, Au21, Au23, Au24, and Au27) have been simulated with a baryonic mass resolution of ∼6 × 10³M_⊙, corresponding to maximum gravitational softening length of 184 physical pc; these are resolution level 3 in the Auriga terminology. By contrast 28 have a baryonic mass resolution of ∼5 × 10⁴M_⊙, corresponding to a maximum softening length of 369 physical pc (level 4); note that the 6 haloes simulated at level 3 were also simulated at level 4. We perform a comparison of the results obtained from the six haloes that are simulated at both resolutions in Appendix  A. There, we show that the results derived at the two resolutions are statistically similar, and thus for the bulk of our analysis in this paper we will use the 28 haloes available at lower resolution due to the greater statistical power they offer. The exception is in Section 3, where we construct mock observations, and the higher spatial resolution is useful for studying observational effects (e.g. beam smearing).

2.2 Data extraction and analysis

For each simulation snapshot, we extract a box centred on the Galactic Centre for analysis; our boxes extend to ±10 kpc on either size of the Galactic plane¹, but the choice of size in the directions parallel to the Galactic plane requires some care. For the purposes of our mock observations of z = 0 galaxies (Section 3), we use a 40 × 40 kpc² box, which is well matched to the field of view (FoV) of MUSE for local galaxy observations, but for the purposes of studying the evolution of the correlation function over cosmological times (Section 4) we instead use a smaller 20 × 20 kpc² region to avoid contamination by mergers and fly-bys, which are much more common at higher redshift; we discuss the choice of box size in more detail in Appendix  B.

Our first step is to resample the star-forming gas cells within the box into 125 × 125 × 125 pc³ cells. We then convert the 3D box into a 2D map by integrating the element (hydrogen and oxygen) mass in all the cells along a line-of-sight normal to the Galactic plane. We obtain a face-on oxygen metallicity map in a commonly used form from the definition of metallicity, 12 + log (N_O/N_H), where N_O is the column density of oxygen nuclei. The end result is a projected metallicity map with a spatial resolution of 125 pc. To complement this, we also produce an SFR map at the same resolution by projecting the total masses of star-forming gas parcels and dividing by the star formation time-scale τ_SF used in the simulations (see Section 2.1).

To quantitatively compare multiple metallicity maps, we extract two-point correlations from metallicity maps, and estimate corresponding correlation lengths. We do this in several steps, following the procedure outlined by Li et al. (2021, 2023). First, we subtract the radially averaged metallicity in each annulus to obtain a metallicity residual map. Next we compute the two-point correlation function of the residual map, which we fit with the functional form for the two-point correlation function predicted by the KT18 model. This model predicts the correlation function in terms of two free parameters: injection width (w_inj) and correlation length (l_corr). Injection width describes the size of the initial bubble formed in explosion events (e.g. supernovae), and is usually too small (≲ 100 pc; Li et al. 2023) to measure in either observations or cosmological simulated galaxies. Correlation length describes a characteristic length of ISM chemical mixing and is a key parameter of the KT18 model. We fit the simulation correlation function to the KT18 functional form using a least-squares approach with w_inj and l_corr as free parameters. We refer to fits on the pure simulation data as the ‘true’ values to distinguish them from the values derived from mock observations as we describe next.

3 MOCK OBSERVATIONS IN THE LOCAL UNIVERSE

Before using the simulations to study the evolution of metallicity correlations over cosmic time, we first confirm that the simulations can reasonably reproduce the metallicity distributions of galaxies at z = 0 derived from observations Li et al. (2021, 2023). In this section, we create mock observations of the Auriga galaxies for direct comparison with previous results, using realistic observational effects [e.g. noise, point spread function (PSF), spatial sampling] drawn from the AMUSING++ survey (López-Cobá et al. 2020). This survey uses MUSE observations, and at the mean 129 Mpc distance of AMUSING++ targets, MUSE’s 1 × 1 arcmin² FoV corresponds to 40 × 40 kpc², which motivates our choice of region to extract from the the Auriga snapshots in Section 2.2. For the purpose of this comparison, we select the six Auriga galaxies at level 3 (high spatial resolution) since they provide the best intrinsic spatial resolution, and thus the most stringent test of the effects of beam-smearing. We list the properties of these six galaxies in Table 1.

We describe the pipeline we use to create the mock observations in Section 3.1, and provide details of our noise model in Section 3.2. We then demonstrate that the mock observations faithfully recover the real correlation lengths in the simulations (Section 3.3), and that these correlation lengths are in reasonable agreement with the values expected for z = 0 galaxies (Section 3.4).

3.1 Pipeline for mock observations

Given our extracted region, we generate mock observations using a pipeline that consists of the following five steps; the final three of these steps closely follow the analysis method described by Li et al. (2023). We also illustrate these steps in Fig. 1, using the simulated galaxy Au6 as an example.

• In order to replicate observational errors on the simulated metallicity map, we first create mock emission-line maps for the lines Hα and [N ii]λ6584, required by the commonly used Pettini & Pagel (2004) metallicity diagnostic (hereafter PPN2). We choose PPN2 because it is one of the simplest diagnostics and requires only two emission lines. To produce synthetic maps in the two required lines, we first derive an Hα map from the SFR map using the calibration suggested by Kennicutt & Evans (2012), SFR/(M_⊙yr⁻¹) = 5 × 10⁻⁴²L_Hα/(ergs⁻¹). Next, we use the metallicity map and the PPN2 diagnostic to predict the [N ii]λ6584/Hα ratio in each pixel by solving the equation 12 + log (O/H) = 9.37 + 2.03x + 1.26x² + 0.32x³, where x =[N ii]λ6584/Hα in each pixel. Multiplying the resulting value of x by the Hα map produces an [N ii]λ6584 line map. We show the resulting Hα and [N ii] line maps in the left column, top two rows of Fig. 1, along with the true metallicity map (left column, bottom row). Empty pixels correspond to locations where the simulations include no fluid elements dense enough to be star forming.

• Next, we convolve the simulated flux maps with a Gaussian PSF and add noise. We convolve the two emission-line maps using a Gaussian beam with a full width at half-maximum of 1|${_{.}^{\prime\prime}}$|0, which is the median PSF size for the AMUSING++ sample (Li et al. 2023). We then generate a noise map using the method described in Section 3.2, and add this to the beam-convolved line maps. At this point, our maps represent a reasonable facsimile of the data products to which we have access from the observations; we illustrate these line flux maps, together with the metallicities one would derive by directly applying the PPN2 diagnostic to them, in the second column in Fig. 1. Blank pixels in the metallicity map correspond to locations where, as a result of noise, the flux of either the Hα or [N ii] line is negative and thus it is not possible to derive a metallicity.

• Our third step is to mimic the data quality cuts that Li et al. (2023) apply to the AMUSING++ data for our simulated maps; they show that these quality cuts are required in order to extract a reliable correlation length from the metallicity field, since in their absence the correlation is corrupted by noise. Specifically, we divide our simulated line and noise maps from step (ii) to produce maps of signal-to-noise ratio (S/N) for both lines, and we mask pixels where the S/N is below 3. We show the masked maps in the third column of Fig. 1.

• Again following the procedure described in Li et al. (2023), we apply the adaptive binning algorithm adabin to the [N ii]λ6584 map, which is the weaker of the two lines in the diagnostic. This algorithm recursively groups pixels into larger and larger blocks in order to increase the S/N; each region is coarsened until it reaches the target S/N of 3. We refer readers to Li et al. (2023) for full details of the algorithm. Once we have re-binned the [N ii] map, we reconstruct the Hα map with the same binning in order to ensure that we only ever compute line ratios at matched spatial resolution. We show the binned maps and the metallicity map derived from them in the fourth column of Fig. 1.

• Our final step is to mask the adaptively binned maps to avoid computing metallicities in locations that go beyond the true boundaries of ionized gas emission in the target galaxy. To achieve this, we mask pixels where in the original, unbinned Hα map, the fluxes are detected at S/N <3. We show the final, masked maps in the final column of Fig. 1.

We estimate the true correlation length in the KT18 model on the metallicity maps at step (i) using a least-squares fit, as discussed in Section 2.2. For the metallicity map at step (v), we instead follow Li et al. (2023) by adopting an Markov chain Monte Carlo (MCMC) approach that includes two additional parameters to describe observational nuisance effects: beam size σ_beam and the noise factor f. Beam size has a Gaussian prior, centred at the known beam size and with dispersion varying in different observations; its purpose in the model is to account for the fact that beam-smearing introduces artificial correlations in the inferred metallicity at small scales. The f factor accounts for the decrease of two-point correlations at non-zero separations caused by noise. Its effect is to reduce the correlation at non-zero separations by a factor of f, which, as Li et al. (2021) show, is how noise affects a measured two-point correlation function. In an extreme case where the noise is significantly larger than signal, the two-point correlation function approaches a δ function, the two-point correlation function of pure noise.

3.2 Noise estimates

The mock pipeline described in the previous section requires estimates of the noise level in the line flux maps. Here, we describe the process by which we construct these maps. We first choose four galaxies from AMUSING++ (López-Cobá et al. 2020) spanning distances ranging from 127 to 131 Mpc, where the spatial resolution of MUSE (0|${_{.}^{\prime\prime}}$|2) matches the 125 pc resolution of the Auriga simulations; the chosen galaxies are SDSSJ102131.91+082419.8 (labelled as ASASSN14ba in the AMUSING++ catalogue), MCG-03-07-040 (labelled as SN2005lu), ESO584-7 (labelled as SN2007ai), and NGC539 (labelled as SN2008gg). For each of these galaxies, we extract one spectrum per spaxel, from which we obtain the flux intensity and its uncertainty for both the Hα and [N ii]λ6584 emission lines. In total, this yields 4 × 320 × 320 ≈ 4 × 10⁵ distinct fluxes and uncertainties, which we use to estimate the typical noise properties given their exposure time (∼1 h).

Fig. 2 shows the distribution of signal and noise for these spaxels. To model the relation between signal and noise at the given distance and spatial resolution, we fit the noise as a function of the signal with a function of the form

where N₀ and S₀ are parameters to be fit. The motivation for this functional form, which is consistent with the distribution shown in Fig. 2, is that a roughly constant background noise level dominates when the signal is weak, while Poisson noise should dominate when the signal is strong. Performing a simple least-squares fit of this model to the measured S/N data yields best-fitting values

where the first number in parentheses is for the Hα line and the second for the [N ii] line. We show these fits by the solid lines in Fig. 2.

To generate our synthetic noise maps, we use these fits to predict the noise level in every pixel of both the Hα and [N ii] maps, and then we draw a noise value for that pixel from a Gaussian distribution with zero mean and a dispersion equal to the noise level.

3.3 Comparison between mock observations and true correlation lengths

We compare the correlation lengths present in the original metallicity maps to those we recover from mock observed metallicity maps (the bottom left and the bottom right panels in Fig. 1); our goal is to establish that the correlation lengths recovered from observed metallicity maps are reasonably accurate estimates of the true correlations, despite the effects of finite telescope resolution and sensitivity. In Fig. 3, we show the two-point correlation functions of the both the true and mock-observed metallicity maps for six example haloes. Fig. 3 illustrates that making a mock observation of the simulation results in suppressing the two-point correlations at all separations larger than zero due to noise, which artificially de-correlates the true metallicity map. It is a common phenomenon that appears when comparing theoretical two-point correlations and real ones. Due to inevitable noise, the two-point correlation function of an observed metallicity map will deviate from that of the true metallicity map. The overall effect of noise is to decrease the correlation function by a constant factor at all separations larger than zero. The f factor in the MCMC fit captures this effect, which is why it is critical to include it.

In Table 1, the column l_{corr, true} reports the true correlation lengths measured directly from the simulations; these are reported without errors, because they come from the χ² fitting and we do not have the uncertainties of the original metallicity maps. The column l_{corr, mock} reports the result from the MCMC fitting, and the reported uncertainties correspond to the 16th to 84th percentile range. Table 1 demonstrates that in most cases the value of l_corr derived from the mock observations is within |$\pm 40~{{\ \rm per\ cent}}$| of the true ones. We can therefore have confidence that the correlation lengths returned from mock observations are reasonably close to reality.

3.4 Comparison between mock observations and real observations

To examine how well our sample Auriga galaxies compare to observations, we show the relationship between l_corr (from mock observations) and M_*, SFR, and R_e for the six example haloes in Fig. 4. The results demonstrate that the Auriga simulations are consistent with observational results, in that the Auriga simulations at z = 0 have correlation lengths well within the range observed for galaxies of similar properties. We warn that, because all the Auriga galaxies are chosen to be Milky Way analogues, this test covers only a limited dynamic range in galaxy properties. However, it is interesting to note that, despite this limited dynamic range, we do recover hints of the correlations seen in the real data, i.e. the Auriga haloes with the largest stellar mass, SFR, and effective radius also tend to have larger correlation lengths.

In Fig. 4, we also indicate by downward arrows our estimate of the metal injection scale in a level 3 Auriga simulation, which represents the scale below which the correlation length will be dominated by numerical resolution; we see that all the measured simulation correlation lengths lie well above this value. To estimate this scale, we first note that in the Auriga simulations, metals are injected into the 4³ cells closest to a supernova event, so the characteristic size of a metal injection region is 4l_g, where l_g is the length-scale of a gas cell. Cell sizes in Auriga are adaptive, so l_g is not fixed, but we can estimate it by considering a gas cell with a density at the star formation threshold ρ₀ = 0.13 H/cm⁻³ (mentioned in Section 2.1) and with the typical baryonic mass resolution at level 3 of the simulations, m_b = 6 × 10³ M_⊙ (Grand et al. 2017). Since |$\rho _0 l_{\rm g}^3 = m_{\rm b}$|⁠, a typical cell satisfying these conditions has l_g ∼ 125 pc, and thus we estimate the metal injection scale to be ∼500 pc.

4 ANALYSIS: THE EVOLUTION OF METALLICITY CORRELATIONS OVER COSMIC TIME

Having established that the metallicity distributions in the Auriga galaxies are a reasonable match to those observed in comparable galaxies at z = 0, we now use the Auriga sample to examine the evolution of the two-point correlation function of metallicities over cosmic time. Since we have also established that simulated observations of the Auriga simulations yield good approximations to the true correlation length, at least for observational parameters typical of current local Universe measurements, the subsequent discussion will only focus on the true two-point correlation functions, not on simulated observations – that is, we base all analysis from this point forward on true metallicity maps such as the one shown in the first column, bottom row of Fig. 1. We have 128 such maps per Auriga halo, at time intervals of ≈160 Myr. For each map of this form, we fit a noise-free KT18 model, following the procedure outlined in Section 3.1, to extract a correlation length l_corr. We therefore obtain the history of l_corr as a function of redshift for each Auriga galaxy.

Fig. 5 shows the star formation histories (SFHs) and l_corr evolution over cosmic time for 28 Auriga haloes at level 4 resolution. Correlation lengths are smoothed in the adjacent three snapshots (corresponding to intervals of ≈±160 Myr on either side of a given time). We also indicate as grey vertical bands in the figure periods of time when massive mergers are taking place, using the merger list identified by Gargiulo et al. (2019). The intervals shown indicate all merger events between the primary Auriga galaxy being simulated and satellites of mass M_sat > 10¹⁰M_⊙.² The width of each grey band is 160 Myr, close to the interval between two snapshots.

As shown in Fig. 5, correlation length fluctuations correlate roughly with SFR fluctuations, but in general are of larger amplitude. Nor is the correlation perfect: while fluctuations of the correlation length in some galaxies resemble amplified versions of the SFR fluctuation, for example, in Au16, there are other galaxies where trends in l_corr and SFR appear to be weakly correlated at best, for example in Au3 and Au19. During merger periods correlation lengths tend to be larger on average than during non-merger periods; we illustrate these averages as solid horizontal lines for merger periods and dashed horizontal lines for non-merger periods. This is consistent with the findings of Li et al. (2021, 2023), who note that correlation lengths tend to be larger in observed galaxies that show signs of being in the process of merging. This increase in correlation length during merger periods is particularly pronounced in some galaxies, for instance Au7, Au21, Au29, and Au30. SFR is also highly sensitive to mergers, which can trigger shock waves through tidal forces, and the common response to l_corr and SFR to merger events clearly accounts for at least some of the correlation between them; Au24 and Au27 provide clear examples. However, it is also clearly not always the case that mergers induce large changes in SFRs or correlation lengths, and sometimes one responds but not the other. For example, both Au18 and Au19 show enhanced SFRs during mergers, but no corresponding increases in correlation length. The examples above demonstrate that neither SFR nor l_corr are perfect indicators of merger events likely due to the complex nature of merger (e.g. mass ratio, gas mass ratio, and orientation).

To explore the co-variation of correlation length and SFR over cosmic time further, in the left panel of Fig. 6, we show a contour plot of the distribution of galaxies in the correlation length–SFR plane broken up into in three temporal bins: high- (purple), medium- (cyan), and low-z (green). The quantity shown is the fraction of lookback time that the galaxy spends at a given combination of SFR and l_corr; thus, for example, the outer purple 86 per cent contour in the figure indicates that, if one chooses a random Auriga galaxy at a random lookback time corresponding to the redshift range z ∈ [2.6, 1.5] and plots its coordinates (SFR, l_corr), there is an 86 per cent chance that these coordinates will lie within the contour. From high to low redshift, we see that the contours are largely consistent. The distribution of the Auriga galaxies in SFR–l_corr does not shift significantly over cosmic time, and in general resembles the sequence traced out in this plane by observed z = 0 galaxies (cf. the middle panel of Fig. 4). Thus, the (very broad) SFR–l_corr relationship appears to be essentially invariant over redshift, at least within the set of galaxies sampled by Auriga, and the evolution of individual galaxies appears to consist primarily of wandering about this relationship. At first, the fact that there is no overall secular trend in where galaxies live on the SFR–l_corr relationship with redshift might seem surprising, but we remind the reader that what we are plotting here is the histories of particular haloes that have been selected specifically to be Milky Way-analogues at z = 0, not the distribution of all galaxies.

However, we caution that the Auriga simulations are by construction all Milky Way analogues, and thus cover a limited range of stellar mass and morphology. We cannot determine with certainty if this conclusion will continue to hold over a wider range of galaxy properties. However, we conjecture that it will, since the observations do include galaxies with wide range of masses and morphologies at z  = 0, and these appear to be part of the same continuous distribution as Milky Way-like galaxies, with no strong dependence on morphology or other galaxy structural parameters (Li et al. 2021).

The right panel of Fig. 6 shows the tracks of two example haloes, Au16 and Au24, in the SFR–l_corr plane, with time flowing from purple dots to yellow dots as indicated by the colour bar. The figure demonstrates that, while on average galaxies circulate rather than migrating systematically toward one end or the other of the SFR–l_corr relationship, they do circulate with a clear pattern, one that is replicated in many other Auriga haloes as well, though we show only two in this figure for clarity. The pattern is that galaxies do not move in both l_corr and SFR simultaneously. Instead, l_corr increases first (corresponding to upwards movement in the figure) and SFR increases next (rightwards movement). The same trend is visible when both of them decrease – l_corr decreases first (downwards movement) and SFR decreases next (leftwards movement). Consequently, the track forms a roughly clockwise cycle in the l_corr versus SFR plane.

To confirm this visual impression, we examine the cross-correlation of SFR and l_corr, defined in the usual way whereby the cross-correlation of two-time sequences X(t) and Y(t) is given by

where 〈 · 〉 denotes an average over t. A positive time lag means that Y reacts first and X follows. Fig. 7 shows the cross-correlations between SFR (as X) and l_corr (as Y) for all the haloes. In most cases, we see that SFR and l_corr are positively correlated at zero time lag, consistent with the overall positively sloped ‘main sequence’ visible in the SFR versus l_corr plane. However, it is also clear that in most cases the cross-correlations have a positive slope, meaning that the correlation is stronger at positive time lag, corresponding to l_corr changing first and the SFR reacting slightly later.

This is an interesting phenomenon because we see for the first time the causality of the interplay between l_corr and SFR. This cannot solely be merger-driven, since we see the same general positive trend in the cross-correlation function for the two haloes that have no mergers (Au17 and Au25) as for all the others. The positive time lag is probably a result of the correlation length and SFR reacting to perturbations on different time-scales – the natural response time for l_corr is a galactic orbital period, while the natural response time for the SFR is the SFR time-scale of ∼2 Gyr, which is much longer than an orbital period. This explains both why l_corr fluctuates more, and why its fluctuations tend to lead SFR fluctuations.

We note that, while the SFR time-scale of ∼2 Gyr is hardwired into the Auriga star formation prescription, and the ∼200 Myr orbital time-scale is implicitly fixed by the choice to simulate Milky Way analogues, the general physical point we make here is true more generally: All observed galaxies have star formation time-scales roughly an order of magnitude longer than their orbital times. This statement is equivalent to the second form of the well-known Kennicutt relation (Σ_SFR versus Σ_g/t_dyn, where Σ_SFR, Σ_g, and t_dyn denote SFR surface density, gas surface density, and local dynamical time-scale, respectively). The relation also holds for high-redshift galaxies (e.g. Daddi et al. 2010).

5 CONCLUSIONS

In this study, we utilize the Auriga simulations to investigate the evolution of the spatial distributions of chemical elements in galaxies using the two-point correlation function as a tool. To quantitatively compare the two-point correlations of different galaxies, and of individual galaxies as they evolve over time, we study the correlation length (l_corr) introduced in the KT18 model for 2D galactic abundance distributions. As a parameter that reveals galaxy chemical mixing mechanisms, the correlation length provides a unique window into galaxy past evolution history.

We first confirm that the Auriga simulations have correlation lengths in galactic metal fields comparable to local observation. We first carry out simulated observations of the z = 0 Auriga snapshots, demonstrating that when we add realistic noise and beam smearing effects, the Auriga simulations produce galaxies with correlation lengths l_corr consistent with observed values for galaxies of similar properties, indicating that the Auriga simulations capture the dominant chemical mixing processes involved in galaxy evolution. We further show that correlation lengths recovered from these simulated observations are reasonably close to the true values present in the simulations.

We find that for an individual galaxy, there is no significant secular evolution in its correlation length, or in its correlation length versus SFR, over cosmological time-scales. The mean correlation lengths of a population of galaxies at high z show only marginal differences from the mean of those same galaxies at low-z, and the galaxy-to-galaxy scatter of l_corr also remains unchanged over cosmic time. We therefore conclude that the relationship between correlation length and SFR is not the result of a build-up over cosmic history. Instead, it appears to be an equilibrium relationship that is established in a time much less than a Hubble time. Galaxies can move along this relationship in response to perturbations in their environments, but the relationship itself is essentially invariant (xxx). This picture is very similar to the one proposed for the origin of galaxy metallicity gradients and their relationship with other galaxy properties such as mass and SFR Sharda et al. (2021).

Furthermore, our analysis reveals an intriguing trend that fluctuations in correlation lengths precede fluctuations in SFR. This finding suggests that the scatter of the l_corr distribution might be a result of galaxies being at different evolutionary stages, whereby galactic correlation lengths react more rapidly to external perturbations than SFRs. The former evolve on time-scales comparable to the galactic orbital period, while the latter are change only over many orbital periods, so at a given time where a galaxy lies in the SFR–l_corr depends in part on whether its SFR has had time to ‘catch up’ to the changes in correlation length induced by whatever perturbed it most recently.

Finally, in the Auriga simulations we have presented, although gas cells in the star-forming disc can be as small as a few tens of pc, metals are injected across 64 cells and therefore the injection width, the other key parameter of the KT18 model with typical values of ≲ 100 pc (Li et al. 2023) is not resolved everywhere. However, in principle, the value of w_inj should vary between elements with different nucleosynthetic origins (e.g. N versus O). One should be able to find signatures of this effect in both observations and higher resolution simulations. In future work, we intend to apply two-point correlations to higher resolution zoom-in simulations, to investigate if they are able to recover the imprints of elements’ differing nucleosynthetic origins on their present-day spatial distributions.

ACKNOWLEDGEMENTS

RJJG acknowledges support from an STFC Ernest Rutherford Fellowship (ST/W003643/1). EW and JTM acknowledge support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. MRK acknowledges support from the Australian Research Council through award FL220100020.

DATA AVAILABILITY

The scripts and plots for this paper will be shared on reasonable request to the corresponding author. The arepo code is publicly available in Weinberger, Springel & Pakmor (2020).

Footnotes

References

APPENDIX A: COMPARISON OF RESULTS DERIVED FROM LEVEL 3 AND LEVEL 4

The majority of the Auriga simulations use ‘level 4’ resolution, for which the maximum softening length is 369 physical pc, but a subset of six haloes (Au6, Au16, Au21, Au23, Au24, and Au27) were also simulated at a resolution with a maximum softening length of 184 physical pc (‘level 3’ resolution). In order to test the resolution dependence of our results, we compute correlation lengths from both sets of simulations. In Fig. A1, we show the cumulative distribution function (CDF) of correlation lengths of the six haloes that are simulated at both level 3 and level 4; the CDF that we show here is the distribution of correlation lengths over all all snapshots for all six haloes, and thus represents the distribution in time. The two distributions are clearly very similar, and a quantitative comparison using a two-sided KS test indicates that we cannot reject the null hypothesis that the level 3 and level 4 simulation results were drawn from the same underlying distribution (p = 0.083). This suggests that our analysis is robust against changing resolutions.

APPENDIX B: COMPARISON OF RESULTS DERIVED FROM DIFFERENT BOX SIZES

As discussed in the main text, there is some subtlety in the choice of FoV around each galaxy to examine when computing the two-point correlation function. A larger FoV of 40 × 40 kpc² at the median ∼129 Mpc distance of the AMUSING++ sample is well matched to the MUSE FoV, and thus is well suited to making mock observations, but we find that a smaller 20 × 20 kpc² FoV is preferable for studying cosmological evolution. To motivate that choice and examine its impact, in Fig. B1, we show the CDF of correlation lengths of all the 28 haloes at level 4 computed using both the larger and smaller FoV. We find that the larger FoV choice leads to a distribution that is shifted significantly to larger correlation lengths. The two-sided KS test indicates that we reject the null hypothesis that the larger and smaller FoV samples were drawn from different parent distributions (p = 6.6 × 10⁻¹⁶). This suggests that the choice of box size will affect our analysis.

Given that the choice of FoV matters, it is important to understand what causes the difference. For this purpose, we show a detailed example of Au24 in Fig. B2. In this figure, we show in the upper panels the metallicity residual maps at four redshifts (z = 2, 1, 0.5, and 0), and in the lower panel the SFH and the correlation length as a function of redshift. We see that at the two lower redshifts, z = 0 and 0.5, the correlation lengths for the two different FoV sizes are nearly identical, and they differ by only tens of per cent even at z = 1. By contrast, there is a very large difference in correlation length in the z = 2 snapshot, and examination of the metallicity maps at SFH makes it clear why: At z = 2, the galaxy is on the cusp of a merger, leading to a substantial increase in SFR. Because the merger has not yet occurred, however, we can clearly see in the larger FoV two distinct galaxies with different mean metallicities. When we then compute the two-point correlation of resulting map, we get a very large value that in effect corresponds to the projected separation between the two pre-merger galaxies, rather than describing the metal field within either of them; a similar phenomenon is seen in observations of interacting galaxies (Li et al. 2023). The smaller FoV avoids this effect because it is not large enough to include gas in satellite galaxies or metal-poor gas from the circumgalactic medium. This phenomenon drives the differences in the distributions of correlation lengths visible in Fig. B1, and it motivates us to choose the smaller box size for our cosmological analysis because we wish to focus on correlation lengths within galaxies, and therefore to the extent possible to avoid contamination by pre-merger interacting systems.