Cosmological evolution of metallicity correlation functions from the Auriga simulations

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\sim 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.


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 ★ E-mail: zefeng.li@anu.edu.auet 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 twopoint 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 et al. 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. (2021Li et al. ( , 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.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 summarise our work in Section 5.

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 ΛCDM cosmological setting, with Planck Collaboration et al. ( 2014) cosmological parameters Ω m = 0.307, Ω b = 0.048, Ω Λ = 0.693, and a Hubble constant of  0 = 100ℎ km s −1 Mpc −1 , where ℎ = 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 9 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 = 0.13 atoms cm −3 ; cells exceeding this threshhold form stars with a star formation timescale 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 halos (Au6, Au16, Au21, Au23, Au24, and Au27) have been simulated with a baryonic mass resolution of ∼ 6 × 10 3 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 4 M ⊙ , corresponding to a maximum softening length of 369 physical pc (level 4); note that the 6 halos simulated at level 3 were also simulated at level 4. We perform a comparison of the results obtained from the six halos 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 halos 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).

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 1 , but the choice of size in the directions parallel to the galactic plane requires some care.For the purposes of our mock observations of  = 0 galaxies (Section 3) we use a 40 × 40 kpc 2 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 2 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 3 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( O / H ), where  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 a star formation rate map at the same resolution by projecting the total masses of star-forming gas parcels and dividing by the star formation timescale  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. (2021Li et al. ( , 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 ( inj ) and correlation length ( 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  inj and  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.

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  = 0 derived from observations Li et al. (2021Li et al. ( , 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, PSF, spatial sampling) drawn from the AMUS-ING++ 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 ′ field of view corresponds to 40 × 40 kpc 2 , 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  [Nii]6584 (lower panel) lines over all pixels extracted from four sample galaxies in AMUSING++; see main text for details.The contours represent the areas that enclose 39% (1 in 2D Gaussian), 68% (2), and 86% (3) of the data; outside the outermost contour, we show individual pixels as black points.The black solid lines show our best flat-plus-linear fit, while the gray dashed lines denote a fixed signal-to-noise ratio of 3.
these correlation lengths are in reasonable agreement with the values expected for  = 0 galaxies (Section 3.4).

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 Figure 1, using the simulated galaxy Au6 as an example.
(i) In order to replicate observational errors on the simulated metallicity map, we first create mock emission line maps for the lines H and [Nii]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 star formation rate map using the calibration suggested by Kennicutt & Evans (2012), SFR/(M ⊙ yr −1 ) = 5 × 10 −42  H /(erg s −1 ).Next, we use the metallicity map and the PPN2 diagnostic to predict the [Nii]6584/H ratio in each pixel by solving the equation 12 + log(O/H) = 9.37 + 2.03 + 1.26 2 + 0.32 3 , where  =[Nii]6584/H in each pixel.Multiplying the resulting value of  by the H map produces an [Nii]6584 line map.We show the resulting H and [Nii] line maps in the left column, top two rows of Figure 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.
(ii) 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 (FWHM) of 1. ′′ 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 Figure 1.Blank pixels in the metallicity map correspond to locations where, as a result of noise, the flux of either the H or [Nii] line is negative and thus it is not possible to derive a metallicity.
(iii) 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-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 Figure 1.
(iv) Again following the procedure described in Li et al. ( 2023) we apply the adaptive binning algorithm adabin to the [Nii]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 [Nii] 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 Figure 1.
(v) Our final step is to mask the adaptively binned maps to avoid computing metallicities in locations that go beyond the true boundaries of ionised 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 Figure 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 MCMC approach that includes two additional parameters to describe observational nuisance effects: beam size  beam and the noise factor  .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  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 Table 1.Global properties and correlation lengths (in both the original simulations and mock observations for interior comparison) of the six Auriga simulations at level 3 (high spatial resolution).Columns are as follows: (1) Auriga ID; (2) stellar mass; (3) star formation rate from stellar cells in the recent 30 Myr; (4) star formation rate from star-forming gas cells, SFR g =  c  g / sf , assuming a typical mass fraction of cold star formation clouds of  c = 0.9 (see details in the bottom panel of Springel & Hernquist's Fig. 1) and a star formation timescale  sf = 2.2 Gyr (Grand et al. 2017); (5) half-stellar-mass radius (effective radius); (6) correlation length in the original simulation; (7) correlation length in the mock observation.For correlation lengths the central value is the 50th percentile of the posterior PDF, and the error bars show the 16th to 84th percentile range.

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 Mpc to 131 Mpc, where the spatial resolution of MUSE (0. ′′ 2) matches the 125 pc resolution of the Auriga simulations; the chosen galaxies are SDSSJ102131.91+082419.8 (labeled as ASASSN14ba in the AMUSING++ catalogue), MCG-03-07-040 (labeled as SN2005lu), ESO584-7 (labeled as SN2007ai), and NGC539 (labeled 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 [Nii]6584 emission lines.In total this yields 4 × 320 × 320 ≈ 4 × 10 5 distinct fluxes and uncertainties, which we use to estimate the typical noise properties given their exposure time (∼ 1 hour).
Figure 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  0 and  0 are parameters to be fit.The motivation for this functional form, which is consistent with the distribution shown in Figure 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-fit values where the first number in parentheses is for the H line and the second for the [Nii] line.We show these fits by the solid lines in Figure 2.
To generate our synthetic noise maps, we use these fits to predict the noise level in every pixel of both the H and [Nii] 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.

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 Figure 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 Figure 3 we show the two-point correlation functions of the both the true and mock-observed metallicity maps for six example halos.Figure 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  factor in the MCMC fit captures this effect, which is why it is critical to include it.
In Table 1 the column  corr,true reports the true correlation lengths measured directly from the simulations; these are reported without errors, because they come from the  2 fitting and we do not have the uncertainties of the original metallicity maps.The column  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  corr derived from the mock observations are within ±40% of the true ones.We can there- fore have confidence that the correlation lengths returned from mock observations are reasonably close to reality.

Comparison between mock observations and real observations
To examine how well our sample Auriga galaxies compare to observations, we show the relationship between  corr (from mock observations) and  * , SFR, and   for the six example halos in Figure 4.
The results demonstrate that the Auriga simulations are consistent with observational results, in that the Auriga simulations at  = 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 halos with the  The blue, cyan, green, yellow, orange, and red (in some panels hidden by the other two symbols) stars show the correlation lengths from the mock observations of the six halos, Au6, Au16, Au21, Au23, Au24, and Au27, respectively.The background heat map shows the distributions derived from 100 galaxies in the CALIFA survey (Li et al. 2021, using the same PPN2 metallicity diagnostic) and the circles show those of 219 galaxies in the AMUSING++ survey (Li et al. 2023).The upper limit on the bottom right shows a scale below which the physics is modified by numerical smoothing at level 3 in the Auriga simulations.largest stellar mass, SFR, and effective radius also tend to have larger correlation lengths.
In Figure 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 3 cells closest to a supernova event, so the characteristic size of a metal injection region is 4 g , where  g is the length scale of a gas cell.Cell sizes in Auriga are adaptive, so  g is not fixed, but we can estimate it by considering a gas cell with a density at the star formation threshold  0 = 0.13 H / cm −3 (mentioned in Section 2.1) and with the typical baryonic mass resolution at level 3 of the simulations,  b = 6 × 10 3 M ⊙ (Grand et al. 2017).Since  0  3 g =  b , a typical cell satisfying these conditions has  g ∼ 125 pc, and thus we estimate the metal injection scale to be ∼ 500 pc.

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  = 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 Figure 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  corr .We therefore obtain the history of  corr as a function of redshift for each Auriga galaxy.
Figure 5 shows the star formation histories (SFHs) and  corr evolution over cosmic time for 28 Auriga halos 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  sat > 10 10 M ⊙2 .The width of each grey band is 160 Myr, close to the interval between two snapshots.
As shown in Figure 5, correlation length fluctuations correlate roughly with star formation rate 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  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) and Li et al. (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  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  corr are perfect indicators of merger events likely due to the complex nature of merger (e.g.mass ratio, gas mass ratio, orientation).
To explore the co-variation of correlation length and SFR over cosmic time further, in the left panel of Figure 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- (green).The quantity shown is the fraction of lookback time that the galaxy spends at a given combination of SFR and  corr ; thus for example the outer purple 86% contour in the figure indicates that, if one chooses a random Auriga galaxy at a random lookback time corresponding to the redshift range  ∈ [2.6, 1.5] and plots its coordinates (SFR,  corr ), there is an 86% 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- corr does not shift significantly over cosmic time, and in general resembles the sequence traced out in this plane by observed  = 0 galaxies (c.f. the middle panel of Figure 4).Thus the (very broad) SFR- 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- corr relationship with redshift might seem surprising, but we remind the reader that what we are plotting here is the histories of particular halos that have been selected specifically to be Milky Way-analogues at  = 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 Figure 6 shows the tracks of two example halos, Au16 and Au24, in the SFR- 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- corr relationship, they do circulate with a clear pattern, one that is replicated in many other Auriga halos as well, though we show only two in this figure for clarity.The pattern is that galaxies do not move in both  corr and SFR simultaneously.Instead,  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 - corr decreases first (downwards movement) and SFR decreases next (leftwards movement).Consequently, the track forms a roughly clock-wise cycle in the  corr versus SFR plane.
To confirm this visual impression, we examine the crosscorrelation of SFR and  corr , defined in the usual way whereby the cross-correlation of two time sequences  () and  () is given by where ⟨•⟩ denotes an average over .A positive time lag means that  reacts first and  follows.Figure 7 shows the cross-correlations between SFR (as ) and  corr (as  ) for all the halos.In most cases, we see that SFR and  corr are positively correlated at zero time lag, consistent with the overall positively-sloped "main sequence" visible in the SFR versus  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  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  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 halos 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  corr is a galactic orbital period, while the natural response time for the SFR is the SFR timescale of ∼ 2 Gyr, which is much longer than an orbital period.This explains both why  corr fluctuates more, and why its fluctuations tend to lead SFR fluctuations.We note that, while the SFR timescale of ∼ 2 Gyr is hardwired into the Auriga star formation prescription, and the ∼ 200 Myr orbital timescale 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 timescales 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 vs. Σ g / dyn , where Σ SFR , Σ g , and  dyn denote star formation rate surface density, gas surface density, and local dynamical timescale, respectively).The relation also holds for high-redshift galaxies (e.g.Daddi et al. 2010).

CONCLUSIONS
In this study, we utilise 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 ( corr ) introduced in the Krumholz & Ting (2018) 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  = 0 Auriga snapshots, demonstrating that when we add realistic noise and beam smearing effects, the Auriga simulations produce galaxies with correlation lengths  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 star formation rate, over cosmological timescales.The mean correlation lengths of a population of galaxies at high- shows only marginal differences from the mean of those same galaxies at low-, and the galaxy-to-galaxy scatter of  corr also remains unchanged over cosmic time.We therefore conclude that the relationship between correlation length and star formation rate 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 star formation rate Sharda et al. (2021).
Furthermore, our analysis reveals an intriguing trend that fluctuations in correlation lengths precede fluctuations in star formation rate.This finding suggests that the scatter of the  corr distribution might be a result of galaxies being at different evolutionary stages, whereby galactic correlation lengths react more rapidly to external perturbations than star formation rates.The former evolve on timescales 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- corr depends in part on whether its star formation rate 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  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. of percent even at  = 1.By contrast, there is a very large difference in correlation length in the  = 2 snapshot, and examination of the metallicity maps at star formation history makes it clear why: at  = 2 the galaxy is on the cusp of a merger, leading to a substantial increase in star formation rate.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 twopoint 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 (CGM).This phenomenon drives the differences in the distributions of correlation lengths visible in Figure 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.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Illustration of the mock observations pipeline applied to the galaxy Au6.Each column shows one step in the procedure described in Section 3.4, respectively; these steps are (i) producing true emission line maps, (ii) convolution of the maps with a synthetic beam and addition of noise, (iii) pruning of low signal-to-noise pixels, (iv) reconstruction of the low signal-to-noise areas using the adabin algorithm; (v) masking the remaining low signal-to-noise regions of the adaptively-binned maps.From top to bottom the rows show H maps, [Nii]6584 maps, and metallicity maps.The first two rows correspond to the colour bar on the upper right, showing line flux emitted per unit area, and the bottom row corresponds to the lower right colour bar, showing metallicity.

Figure 2 .
Figure 2. Distribution of signal versus noise for the H (upper panel) and[Nii]6584 (lower panel) lines over all pixels extracted from four sample galaxies in AMUSING++; see main text for details.The contours represent the areas that enclose 39% (1 in 2D Gaussian), 68% (2), and 86% (3) of the data; outside the outermost contour, we show individual pixels as black points.The black solid lines show our best flat-plus-linear fit, while the gray dashed lines denote a fixed signal-to-noise ratio of 3.

Figure 3 .
Figure3.The two-point correlation functions of the six example simulated galaxies, and the comparison between the two-point correlation functions of the pure simulations (blue dots) and those of mock observations (red dots).The blue solid lines and red solid lines represent the best estimate of the parameters in the KT18 model.Specifically the effects from the mock observations are shown as the discontinuity between the first two red dots.

Figure 4 .
Figure 4. Correlation length versus stellar mass (left), SFR (middle), and   (right).The blue, cyan, green, yellow, orange, and red (in some panels hidden by the other two symbols) stars show the correlation lengths from the mock observations of the six halos, Au6, Au16, Au21, Au23, Au24, and Au27, respectively.The background heat map shows the distributions derived from 100 galaxies in the CALIFA survey(Li et al. 2021, using the same PPN2 metallicity diagnostic) and the circles show those of 219 galaxies in the AMUSING++ survey(Li et al. 2023).The upper limit on the bottom right shows a scale below which the physics is modified by numerical smoothing at level 3 in the Auriga simulations.

Figure 5 .
Figure 5. Cosmological evolution of star formation rates (blue solid lines) and correlation lengths (red dashed lines) for all the halos.The grey vertical bands indicate periods when major mergers take place.The horizontal solid and dashed lines represent the averaged values of correlation length during merger and non-merger periods, respectively, and are plotted only for galaxies that experience at least one major merger.

Figure 6 .Figure 7 .
Figure 6.Left: contours showing the distribution of all the halos in the SFR- corr plane, divided up in bins of redshift.To construct this figure, we take the  corr and SFR values from all snapshots in three temporal bins - ∈ [2.6, 1.5] (cosmic noon),  ∈ [1.5, 0.5] (disc settling), and  ∈ [0.5, 0.0] (disc settled) -and count the frequency in a 2D histogram in this plane with bins that are 0.14 dex wide in SFR and 0.12 dex wide in  corr .We then construct the contours shown from the 2D histograms.Each contour represents the boundary of the area within which 86% (darker and outer) and 39% (lighter and filled) data are included -these contour levels correspond to 2 and 1 for a 2D Gaussian.The contours are coloured according to the centres of the corresponding time intervals on the colour bar to the right.Right: the trajectory of two example halos, Au16 (solid) and Au24 (dashed), in the SFR- corr plane.Colours on the tracks indicate lookback time as shown on the colour bar to the right.The starting and ending positions are marked in circles.

Figure A1 .Figure B2 .
Figure A1.Cumulative distribution function of correlation lengths for the six Auriga halos that are simulated at both level 3 (blue) and level 4 (red) resolution.