Cosmological gas accretion history onto the stellar discs of Milky Way-like galaxies in the Auriga simulations -- (II) The inside-out growth of discs

We investigate the growth of stellar discs in Milky Way-mass galaxies using the magnetohydrodynamical simulations of the Auriga Project in a full cosmological context. We focus on the gas accretion process along the discs, calculating the net, infall and outflow rates as a function of galactocentric distance, and investigate the relation between them and the star formation activity. The stellar distributions of around 70% of the simulated galaxies exhibit an ``inside-out'' pattern, with older (younger) stellar populations preferentially located in the inner (outer) disc regions. In all cases, we find a very tight correlation between the infall, outflow and net accretion rates, as well as between these three quantities and the star formation rate. This is because the amount of gas which is ultimately available for star formation in each radial ring depends not only on the infall rates, but also on the amount of gas leaving the disc in outflows, which directly relates to the local star formation level. Therefore, any of these rates can be used to identify galaxies with inside-out growth. For these galaxies, the correlation between the dominant times of accretion/star formation and disc radius is well fitted by a linear function. We also find that, when averaged over galaxies with formation histories similar to the Milky Way, the simulated accretion rates show a similar evolution (both temporally- and radially-integrated) to the usual accretion prescriptions used in chemical evolution models, although some major differences arise at early times and in the inner disc regions.


INTRODUCTION
The physical processes that lead to the formation and govern the evolution of galaxies have long been one of the main quests in modern astrophysics.In particular, significant attention has been given to understand the formation of the Milky Way (MW) and similar spirals, as it is for the MW that we have access to the most detailed observations.The star formation history is one of the key processes that influence the final properties of galaxies, and results from the complex interplay of the various physical processes which occur during galaxy evolution.In this context, the variety of observed morphologies is a natural result of galaxy formation with different merger and accretion histories, resulting into a combination of disc and bulge components with varying relative significance.While the formation of bulges is thought to occur early and in short time-scales, the discs ★ E-mail: fiza@iafe.uba.ar are expected to form over longer periods, fed by gas accretion at increasingly large distances.
In this context, the stellar discs of spiral galaxies are thought to grow from the "inside-out" (Larson 1976;Fall & Efstathiou 1980), as old (young) stars preferentially populate the inner (outer) regions of galaxy discs.This is suggested by a number of observations (see, e.g.Frankel et al. 2019;Prantzos et al. 2023), although inside-out growth might also depend on galaxy mass (Pan et al. 2015).Most observations, however, rely on the present-day positions of stars instead of their unknown birth sites, but stars are subject to orbit mixing processes, therefore losing memory on their initial location.Unlike dynamical memory, stars retain their intrinsic chemical patterns during their whole evolution, thus providing information on the properties of stars at their birth time.As a result, the measurement of detailed chemical abundances of stars in the MW and other galaxies has opened up the possibility to investigate in detail the formation of the stellar components of spiral galaxies, especially their discs, and compare with the observed chemical abundances.
In particular, chemical evolution models (CEMs) of the MW show that observed galaxy properties can be reproduced, but the predicted stellar metallicity gradients in the discs are directly affected by the assumptions made on gas accretion rates as a function of time and radius within the disc, as gaseous material will ultimately give rise to the formation of new stellar populations (e.g.Matteucci 2012).In order to reproduce observational results, a first, rapid episode of gas accretion is required in CEMs to explain the formation of the spheroidal components of the Galaxy.In contrast, accretion onto the thin-disc component is required to occur at decreasing rates during longer timescales, which is usually modelled assuming cosmic timedependent exponential or Gaussian functions (see e.g.Iza et al. 2022, and references therein).Moreover, to account for the fact that less enriched (i.e.older) stars are formed towards the inner regions of the Galaxy, these models also assume that gas accretion timescales correlate with radius, thus delaying the onset of star formation in the disc outskirts (e.g Chiappini et al. 2001).While this inside-out formation assumption for the disc is key in these models, further ingredients such as radial gas flows and/or variable efficiency of star formation might also play a role (see, e.g.Palla et al. 2020).The combination of various assumptions is then needed in order to successfully reproduce the observed distribution of present-day chemical abundances in the Galactic disc and the photometric properties of galaxies similar to the MW (Boissier & Prantzos 1999).However, the validity of these assumptions has not always been tested using more realistic simulations in a cosmological context that model gas flows into and out of galaxies self-consistently over their whole evolution.
Together with advances in observational studies and chemical evolution models, numerical simulations have made significant progress during the last decades, being able to reproduce the formation of galaxy discs from cosmological initial conditions (e.g.Okamoto et al. 2005, Scannapieco et al. 2008, Guedes et al. 2011, Stinson et al. 2013, Aumer et al. 2013, Marinacci et al. 2014, Wang et al. 2015, Hopkins et al. 2018).In combination with detailed modelling of the chemical evolution of galaxies, these are a powerful tool to investigate the relation between star formation, gas accretion and morphology across cosmic time, and to study in more detail the formation of stellar discs.In cosmological simulations, relations between gas accretion and morphology have already been found, for instance, by Aumer et al. (2014).As the gas component is the fuel for star formation, the star formation rate (SFR) of a galaxy will be naturally linked to its accretion history, even though processes such as feedback can make this relation somewhat complex.In any case, if discs form from the inside-out, one might expect that gas accretion has an inside-out behaviour as well.Recent simulation studies suggest that discs, at least at the MW scale, form from the inside-out, even though not all of them necessarily follow this pattern (e.g.Sommer-Larsen et al. 2003;Scannapieco et al. 2009;Aumer et al. 2014;Gómez et al. 2017;Nuza et al. 2019).
In Iza et al. (2022), the first paper of this series (Paper I), we investigated the temporal evolution of gas accretion onto the discs of simulated MW-mass galaxies.For this, we used the original simulations of the Auriga Project (Grand et al. 2017), which consist on a sample of 30 zoom-in, MW-mass galaxy haloes simulated with the magnetohydrodynamic code arepo (Springel 2010) that were selected to be relatively isolated at the present day.The sample was separated in two different galaxy populations characterised by the behaviour of their late-time net gas accretion, with the vast majority (group G1) displaying a gentle, decaying exponential-like behaviour consistent with the expectation of CEMs, and the remaining (group G2) an increasing net accretion related to a more active merger history at late times which lead to the partial/total destruction of their discs, even though most of them ended up with new well-formed disc-like components at  = 0. We also found a strong correlation between the gas inflow onto the disc region and the SFR, especially at early times, confirming the close relationship between both quantities.Furthermore, outflows also showed a strong correlation with SFR, highlighting the link between SFR and the generation of outflows following star formation, despite the fact that these processes are related in a non-trivial manner, affecting the circulation of gas in the disc-halo interface.
In this work, we focus on the evolution and growth of stellar discs as a function of radius and evaluate whether the simulations are consistent with an inside-out formation scenario by analysing both the distribution of stellar ages along the plane of the disc as well as the radial dependency of the gas accretion rates.We quantify the insideout behaviour of each galaxy through an inside-out parameter which provides information on the dominant times of accretion for each radial ring and can be computed for the infall, outflow, net accretion and star formation rates.Since accretion laws are a fundamental ingredient in CEMs that aim to describe the distribution of metals in galaxies, we also compare the gas accretion results obtained from simulations with the usual assumptions adopted in CEMs.
The organisation of this paper is as follows: in Section 2 we briefly revisit the Auriga simulations, and discuss the methods used in the analysis of the "inside-out" behaviour; in Section 3 we focus on our main results: the distribution of stellar ages along the disc, the radial dependency of the inflow, outflow and net gas accretion rates and the star formation rate (SFR), and present the average behaviour over sub-samples of the Auriga galaxies, and in Section 4 we show a simple comparison of the obtained accretion rates, integrated over time and radius, with the expected accretion law from CEMs adopting the parameters of our simulated MW-like galaxy sample.Finally, in Section 5, we summarise and discuss our main results.

The Auriga Simulations
The Auriga simulations (Grand et al. 2017) consist on a suite of 30 Milky Way-mass disc galaxies simulated in a full cosmological setting at high resolution using the zoom-in technique.The simulations were run with the quasi-Lagrangian, moving-mesh magnetohydrodynamic (MHD) code arepo (Springel 2010).Gravitational forces in arepo are computed with a standard TreePM method, employing a fast Fourier Transform (FFT) technique and an oct-tree algorithm for long-and short-range forces, respectively.The MHD equations for the collisional component are discretized on a dynamic unstructured Voronoi mesh.Further details about the implementation and usage of the arepo code can be found in Springel (2010).
The physical model adopted in the Auriga simulations includes a variety of physical processes such as metal-line cooling, a uniform ultraviolet background field for reionization, star formation based on the Springel & Hernquist (2003) model, magnetic fields (Pakmor et al. 2014(Pakmor et al. , 2017(Pakmor et al. , 2018)), energetic and chemical feedback from Type II supernovae, metal return from Type Ia supernovae and asymptotic giant branch stars and AGN feedback (Vogelsberger et al. 2013;Marinacci et al. 2014;Grand et al. 2017).
For this study, we use the level 4 resolution runs of Grand et al. (2017) with a mass resolution of ∼ 3 × 10 5 M ⊙ and ∼ 5 × 10 4 M ⊙ for dark matter and baryons, respectively.The adopted cosmological pa-Table 1. Galactic properties at  = 0.The columns are: (1) galaxy name, (2) virial radius  200 , (3) virial mass  200 , (4) subhalo stellar mass  ★ , (5) subhalo gas mass  gas , (6) disc-to-total mass fraction, (7) disc radius  d and (8) disc height ℎ d .The symbol † in the first column identifies galaxies that include a treatment for stochastic tracer particles., where ℎ = 0.6777.The 30 Auriga galaxies have been selected for re-simulation at  = 0 from a dark matter-only simulation (Schaye et al. 2015), with the conditions of having virial masses in the range 1-2 × 10 12 M ⊙ , and being relatively isolated at the present time, with no other halo with a mass higher than 3% of its own mass located within 9 times its virial radius.From all haloes satisfying these conditions, the final Auriga sample was randomly selected from the most isolated quartile.Further details on the Auriga simulations can be found in Grand et al. (2017).

Selection and properties of the Auriga galaxies
In this work, we extend the analysis of gas accretion rates on the discs presented in Paper I, where we focused on the temporal dependencies of the infall, outflow and net rates.Here, we investigate the corresponding radial dependencies, with focus on the spatial build-up of gaseous and stellar discs and the possibility that some or most of the Auriga discs formed from the inside-out in agreement with usual assumptions for the MW.
In Paper I, we separated the Auriga galaxies in two groups according to the behaviour of the net accretion rates: galaxies with declining (constant/increasing) late-time accretion rates were classified as part of group G1 (G2).Furthermore, G1 galaxies were identified as MW-analogues, having quiet merger histories, while the majority of galaxies in G2 experience significant merger events after the initial formation period.In both groups, however, galaxies ended up as disc-dominated galaxies at the present.As part of our classification, we also identified 6 galaxies with strongly perturbed stellar discs that could affect our analysis methods and were excluded: Au1, Au5, Au19, Au28, Au29 and Au30.For consistency with Paper I, we exclude these galaxies in the present work as well.As a result, we are left with a total of 24 (19 in G1 and 5 in G2) galaxies, with virial masses in the range ∼ 1-1.7 × 10 12 M ⊙ and stellar masses in the range ∼ 4-12 × 10 10 M ⊙ .
Table 1 lists the main properties of the galaxies studied in this work.We note that all of them have virial masses consistent with the commonly accepted mass of ∼ 10 12 M ⊙ for the MW and can be thought of as "Milky Way-like" galaxies in terms of their mass but they have been chosen to be relatively isolated at  = 0, unlike in the real environment of our Galaxy.The table also provides information on other  = 0 properties, including an estimation of the radius and height of the galactic discs (see Section 2.3.1 and Paper I for details on these calculations).
An important aspect of these simulations is that, owing to the quasi-Lagrangian nature of arepo, it is not possible to trace the trajectories of gas elements through cosmic time, and thus relevant information for our analysis of gas accretion is lost.However, 7 of the selected galaxies have been re-simulated including a treatment for stochastic tracer particles (Genel et al. 2013;DeFelippis et al. 2017), making the study of the history of gas elements possible.These galaxies are identified with the symbol † in Table 1.As in a standard Lagrangian scheme, tracer particles can be tracked back through cosmic time allowing us to assess the evolution of gas accretion.Further details on the simulations with tracer particles can be found in Grand et al. (2019) and Fragkoudi et al. (2021).Projected density maps for the Auriga galaxies at  = 0, for the stellar (left) and gaseous (right) components.For each galaxy we show the face-on view (top) and edge-on view (bottom).The colour maps span five orders of magnitude in projected density using a logarithmic scale.The face-on view shows a cubic region of 100 ckpc on a side and the edge-on view a region of dimensions 100 × 100 × 20 ckpc 3 .

Defining the stellar disc
Computing gas accretion rates onto the discs of the Auriga galaxies as a function of cosmic time requires a consistent identification of the disc region in the simulations and a proper evaluation of its morphological properties, not only in the present, but throughout its whole evolution.In this section, we outline the procedure used to define the disc region where inflow and outflow patterns are measured.
We first calculate, for each galaxy at a given cosmic time, the inertia tensor of all stars in the inner 10 ckpc and take the principal axis of inertia to coincide with the -axis, in a way that the latter also corresponds to the direction of the angular momentum vector of stars.Thus, after properly rotating the reference frame of the galaxy, the galactic disc is contained in the  plane.
In Paper I we introduced a simple definition to calculate disc sizes that includes two parameters: a disc radius,  d , and a disc height, ℎ d .The first is defined, at each time, as the radius that encloses 90% of the total host stellar mass.In a similar fashion, the disc height for the positive ( > 0) and negative ( < 0) regions in the -direction above and below the mid plane is computed as the distance that encloses 90% of the stellar mass and we consider the disc height as the mean of both values.In this way, ℎ d refers to the height above and below the disc plane resulting in a disc width 2ℎ d .For more details on these calculations we refer the reader to our Paper I. This method works well for all galaxies at all times, except for very early epochs when discs have not formed yet and strong asymmetries are present in the stellar mass distribution.We include, therefore, some corrections for cosmic time less than 4 Gyr together with the changes shown in disc sizes and accretion rates when considering other stellar mass fractions to define the disc (see Paper I for details).Figures 1 and 2 show the face-on and edge-on stellar and gaseous density maps for the simulated galaxies.

Calculation of infall, outflow and net accretion rates
In this section, we describe the method used to quantify the accretion rates onto the galactic discs.To properly consider differences between galactic disc sizes and their time evolution, we divide each disc region in ten radial bins of width 0.1 d and height ℎ d .Both of these parameters were calculated in Paper I for all simulations to which we refer the reader for further details.
In the simulations that do contain tracer particles, we compute, for each radial ring, the number of tracer particles that were previously outside the stellar disc that are later found inside in the next snapshot.Since each tracer is considered to carry a fixed amount of mass, the total infalling mass between snapshots corresponds to the number of infalling tracers multiplied by the tracer mass.Outflowing mass, on the other hand, is calculated considering the amount of tracers that were previously inside a given ring that are later found outside the stellar disc.To compute the corresponding gas rates, the inflowing/outflowing gas masses are divided by the time span of the two considered snapshots.
Note that the inflow and outflow rates do not include radial mass flows.Including radial flows has, by definition, the effect of increasing the inflow/outflow rates, but does not produce any important variation to the results obtained in this work, particularly with our findings related to the inside-out growth of the stellar discs.An analysis of radial flows of gas inside the discs of the Auriga galaxies can be found in Okalidis et al. (2021).
The net accretion rate, on the other hand, is related to the change of gas mass inside each ring, which is, in turn, affected by both infalling/outflowing gas particles and star formation, being the latter responsible of gas depletion within the disc.To quantify  net we use the formula detailed in Paper I, which uses information of the cells and can, therefore, be applied to all the galaxies in the sample.Then, for each radial ring, the net gas accretion rate is given by where  gas is the gas mass in the stellar disc,  and  − 1 are indexes labelling consecutive snapshots,  ★ is the amount of stellar mass formed between them and  is the cosmic time.This formula accounts for the loss of gas mass in the stellar disc owing to feedback processes removing material from the disc and star formation.Furthermore, in Paper I we showed that this way of calculating the net accretion rates when tracer particles are not present yields similar results when compared to the "exact", particle-by-particle, method.

Quantifying inside-out growth with the mass-weighted time
After calculating the infall, outflow and net accretion rates, we would like to study the validity of the inside-out formation scenario within the context of cosmological simulations.In this picture, usually accepted in the formation of MW-like galaxies, the stellar disc of galaxies forms their inner regions first and the outer regions at later times.This is usually quantified by studying the radial profiles of stellar age with cosmic time, but as we showed in the previous paper of this series, there is a strong correlation between the amount of inflowing material and star formation rate since the latter is induced by the presence of new material that has fallen onto the disc.
We introduce a parameter that allows us to quantify the characteristic time at which a mass rate (the inflow rate, for example) is concentrated onto a specific disc region.Then, we define the massweighted time  as where  is a given mass rate (an amount of mass per unit time, used as a weight),  is the cosmic time,   is a reference time from which we start the calculation and  0 is the age of the universe.As shown in Paper I, stellar discs in the Auriga simulations are well defined only after the first ∼ 4 Gyr of evolution.Hence, in what follows, we adopt   = 4 Gyr for the starting point in the integrals above.
Note that this parameter is essentially the mean value of the cosmic time weighted by the corresponding mass rate.Therefore, for a given mass rate,  indicates the mean time at which the mass rate is concentrated and, as we will show in Section 3, it can be used to quantify the inside-out behaviour of each galaxy.

RESULTS
In this Section we discuss whether the simulated discs in the Auriga galaxies formed from the inside-out.In the context of galaxy formation, an inside-out scenario usually refers to an anti-correlation between the age of stars and their location onto the disc: in discs formed from the inside-out increasingly younger stars populate the outer regions of discs.If stellar migration is not significant, insideout will imply that stars in the inner (outer) regions of discs formed early (late).As star formation is mainly determined by the availability of dense gas, a relation would then be expected between the amount of accreted gas as a function of radii and the star formation rate along the disc, which would be valid across cosmic history.The following sections explore the inside-out formation scenario in the auriga simulations, in terms of the stellar populations and the accretion history.

Inside-out formation of the stellar discs
Fig. 3 shows the mean stellar age of stars as a function of projected radius at  = 0 (normalised by the corresponding disc radii), for all stars (black lines) as well as separating stars into disc and spheroidal components (blue and red lines, respectively) 1 .In approximately 70% of the galaxies, we observe a clear inside-out pattern for the disc components: the mean stellar age anti-correlates with radius.This means that the oldest stars tend to be located in the inner regions, while they get increasingly younger towards larger radii.In many other cases, we observe a similar behaviour up to about half the disc radius and the trend inverts in the outermost disc regions.As expected, the spheroidal components in most cases show no correlation between mean stellar age and radius, and their populations are significantly older compared to the discs (Monachesi et al. 2019).Nevertheless, some galaxies (e.g.Au9, Au11, Au17, and Au22) show that stellar populations in the inner regions are -on average -younger compared to the outer regions.This results from ongoing star formation in the inner parts of these galaxies 2 .
In the next sections, we investigate whether the behaviour observed for galaxies is consistent or not with an inside-out behaviour, and if this dichotomy could be explained by differences in their accretion patterns.

Inflow and outflow rates
The radial dependency of gas inflow and outflows rates onto the stellar discs of the 7 simulations including tracer particles are analysed in this Section.As shown in Paper I, the inflow and outflow rates follow very similar patterns, with the former being always dominant over the latter.For this reason, we only show the results obtained for the case of the inflow rates.Fig. 4 shows the evolution of the inflow rate, separated into various curves which represent the accretion onto different radial bins normalised with the corresponding disc radius ( d ) of each galaxy.Each radial ring is referenced by its maximum radius, and colour-coded using the colour map indicated in the legend.
In all cases except for Au13 and Au17, we observe that, for late times, the inflow rates are higher for larger radii, while the opposite occurs at early times (note that we focus on times  ≳ 2 Gyr, i.e. when galaxies start to form their discs).This inversion typically takes place in the range  ∼ 4 − 8 Gyr.At times ≳ 10 Gyr, the infall rates for the different radial bins differ in more than an order of magnitude, while at earlier times the differences are less important.Au13 and Au17, on the other hand, show a more complex behaviour.In the first case, infalling rates tend to be higher in the outer regions except for times between ∼ 11 Gyr and the current epoch.In Au17 a similar inversion seems to take place around ∼ 6 Gyr and up to ∼ 11 Gyr.We come back to possible reasons for the behaviour of these galaxies in the next sections.
The mass-weighted time  defined in Section 2.3.3 can be used to estimate the typical times at which inflows occur for different radial bins and identify galaxies that are consistent with an inside-out formation scenario.In Fig. 5, we show  as a function of ring radius 1 Disc stars were identified as those that are rotating in disc-like orbits.In particular, stars were considered as part of the disc if they have a circularity parameter  > 0.5 and part of the spheroid if |  | ≤ 0.5.As explained in Paper I (see also Scannapieco et al. 2009), the circularity parameter is defined as the ratio between the angular momentum of each star in the -direction (i.e.perpendicular to the disc plane) and the angular momentum expected for a circular orbit at the star's radius.For this reason, disc stars will have  ≈ 1 while a spheroid would be visible in the circularity distribution around  ≈ 0 if it is non-rotating. 2We discarded contamination of disc particles into the spheroidal component as a possible reason for this behaviour, which could in principle result from our simple separation of the disc and spheroidal components.

Age [Gyr]
All for the 7 simulations previously discussed.In this case, we include the results for the inflow rates of Fig. 4, as well as for the outflows.All galaxies except for Au13 and Au17 exhibit a positive correlation between  and radius for both inflows and outflows, confirming the results of the previous figure which suggested an inside-out pattern.This means that, as we go from the inner to the outer regions of the discs, the dominant times at which inflows/outflows occur move from earlier (∼ 5 − 8 Gyr) to later times (≳ 9 Gyr).Au13 and Au17, on the other hand, show nearly constant and linearly decaying profiles, respectively.For Au13, the accretion onto the disc seems to concentrate at ∼ 9 Gyr regardless of the radius; in the case of Au17, accretion onto the inner (outer) region is concentrated at ∼ 8-9 Gyr (∼ 7-8 Gyr) (see Fig. 7).
Another useful quantity to measure the level of the inside-out Temporal evolution of the gas inflow rate for different radial bins.Each line in this figure represents a different radial ring in the disc plane (the width is 0.1 d for each ring); the first line (labelled as 0.1 d in the legend), for example, represents the ring that extends from    = 0 to    = 0.1 d .The ten lines shown in each panel cover the complete radial extension of the disc for each galaxy.In this figure, the inside-out behaviour is evidenced by the relative heights of the curves for each radial bin; for Au26, for example, inflow rates are roughly similar for all bins before ∼ 7 Gyr but higher rates are present for the outer regions near the present.pattern is the slope of the -radius relation (the radius is normalised to the disc radius of each galaxy at any given time), which we refer to as the inside-out parameter .This parameter is calculated using a linear fit of the mass-weighted time as a function of radius and is given by where  is the radius on the disc plane (normalised by the disc radius, as we consider in the following sections),  is the mass-weighted time as defined in Eq. ( 2),  is the amount of data points, and ⟨⟩ denotes the average.The values of  are included in each panel of Fig. 5 along with their uncertainties.The typical  values are between ∼ 2 and 5 Gyr, and the  values calculated from infalling and outflowing material are consistent with each other3 .In the case of Au13 and Au17, whose behaviour differ from that of the rest of the galaxies, we find an  ≲ 0 value.In particular,  ≈ 0 for Au13, both in terms of the inflow and outflow rates, while Au17 exhibits negative  values.
It is worth noting that 5 out of 7 galaxies are consistent with an inside-out formation scenario which, as discussed above, provides a simple picture for the disc assembly of MW-type galaxies which is usually assumed in models of the formation of the MW.Furthermore, galaxies identified as inside-out in terms of the stellar disc (see Fig. 3) are also identified as inside-out from the infall rates, and the same is valid for galaxies that do not seem to have formed in an inside-out manner.
Our simulations, which involve MW-mass galaxies with different formation histories, show that galaxy-to-galaxy variations can translate into different values of the inside-out parameter  and that galaxies inconsistent with an inside-out behaviour can, in fact, also form (see Nuza et al. 2019).Estimating the relative fraction of galaxies with an inside-out pattern would require higher statistics.In the next section we will assess the previous point using the net rates of the Auriga sample, for which we have a higher number of simulated galaxies.

Net rates
As explained in Paper I and Section 2.3.2, the full Auriga sample can be used to estimate net inflow rates as a function of time4 .In this section, we investigate the net inflow rates as a function of the disc radius.Fig. 6 shows the net accretion rates as a function of time for the simulated galaxies, separated in various radial bins covering the disc extension at each time.Similarly to our results for the infall rates, we find that most of the galaxies are consistent with an insideout pattern of the net accretion rates: accretion at early (late) times occurs preferentially at the innermost (outermost) radial bins.This behaviour can be quantified using the mass-weighted time  versus radius relation, which we show in Fig. 7, and the corresponding values

Au26
Inflows Outflows Figure 5. Mass-weighted time  as a function of normalised radius for the re-simulated haloes.Each panel shows  calculated adopting the inflow rate (red) and the outflow rate (blue) as a weight, as well as a linear fit for each case (dashed lines, corresponding colours).On the bottom right of each panel we indicate the slopes of the fits, along with its uncertainty (corresponding colours).Increasing (decreasing) profiles indicate that outer regions concentrate accretion at later (earlier) times than inner regions.Therefore, inside-out galaxies are evidenced by positive slopes.
for the slope of this relation, i.e., the already mentioned inside-out parameter .We find that the net accretion rates are in general compatible with an inside-out formation scenario, with  values between ≈ 1-3 Gyr.Some galaxies have, however, inside-out parameters consistent with zero or negative values.In Table 2 we show , its standard deviation  and the correlation parameter between the variables for all galaxies.Using the  parameter calculated with the net accretion rates, we classify galaxies as "inside-out" if  ≥  th , assuming  th = 0.8 Gyr as a threshold value.With this choice, 16 out of the 24 galaxies fall in this category, wile the remaining 8 galaxies exhibit  <  th , 3 of which have  < −0.5 Gyr5 .
As discussed in Sec.2.2, it is worth noting that, in Paper I, we separated our galaxy sample into two groups: G1 formed by galaxies identified as MW analogues (systems with well-defined, stable discs at  = 0 and no strong perturbations during the last 8 Gyr) and G2 composed of galaxies with well-defined discs at  = 0, but with episodes of partial/total destruction of the discs after the first 4 Gyr of evolution.Given the small number statistics in group G2, we cannot confirm any systematic difference in the behaviour of galaxies between both groups (see Table 2) in terms of their net accretion rates as a function of radius.However, it is worth noting that 4 out of the 5 galaxies in G2 have very clear inside-out patterns, with  > 2 Gyr, and are among the galaxies with highest  values.However, we have checked that if we restrict the sample to times after the formation of the disc that survives up to the present time, the obtained  values decrease and become similar to the values obtained for galaxies in G1.
Among the galaxies that are incompatible with an inside-out formation scenario (Au4, Au10, Au11, Au13, Au14, Au17, Au22 and Au23), the most extreme cases in terms of the  parameter are Au4 and Au17.The latter has the highest slope (in absolute value), with  = −1.0Gyr, and has several distinctive features compared to the rest of the sample: it hosts a massive black hole, it has the highest accretion rate at early times, its stellar component formed very early on, is one of the most compact galaxies, and its environment is the densest one at very early times.Au4, on the other hand, has an  = −0.7 Gyr, but in this case it suffered a merger approximately 3 Gyr ago, shows a perturbed stellar distribution until the present time and is a very compact galaxy.Au10, Au11, Au13 and Au14 have all suffered merger or interaction events and also show perturbed stellar distributions at late times, particularly after 10 Gyr.In the case of Au11, a merger is occurring at  = 0, as can be seen from Fig. 1.It is worth noting that mergers might affect the orientation of the galaxy and the  calculation.Finally, Au22 and Au23 have  values close to the assumed threshold of 0.8, show quiescent formation histories and are similar to our sample of inside-out galaxies.
While most of the galaxies inconsistent with an inside-out behaviour have experienced recent mergers or have nearby satellites at the present day, no other systematics have been found in terms of virial masses, stellar masses, disc-to-total mass ratio, black hole mass, or even the presence/absence of a bar, as we show in Appendix A. However, it is worth noting that the vast majority of them have high mean stellar formation times, that are systematically higher than those of the rest of the sample (Fig. A).This might suggest that discs have not yet fully formed or had enough time to grow and generate an inside-out pattern.Furthermore, in Appendix B we show that, although these galaxies are not formally considered as "inside-out" globally (i.e., for our standard interval between 4 Gyr and the present time), there are certain periods in which they show an inside-out be-Table 2. Inside-out parameter  for each galaxy.We indicate the value of , along with its uncertainty, for the inflow-dominated times of the net accretion rate (an its correlation coefficient), for the SFR, and for the inflows and outflows for the simulations that include tracer particles.

Galaxy
haviour, and that these periods often coincide with epochs of smooth disc growth, as opposed to intervals with an active merger activity.

Star formation rates and relation to accretion rates
In this section, we investigate in more detail the relation between the star formation and the accretion rates along the discs.Note that while we have already seen in Fig. 3 which galaxies are compatible with an inside-out formation of the stellar disc, for consistency we also quantify this behaviour using the same methods as for the accretion rates.Fig. 8 shows the radial profile of the mass-weighted time calculated adopting the star formation rate as a weight of the simulated galaxies, along with the corresponding  values.For reference, we overplot the results corresponding to the -radius relation of the net accretion rates.We find that the star formation and net accretion rates are correlated, with very similar patterns and consistent  values.This result shows that gas inflows contribute with material which directly participates in the formation of new stars in the discs.It is worth noting that, while in Paper I we also found a strong correlation between the evolution of the net accretion and star formation rates integrated over the disc extent, we find that these two quantities are also correlated in terms of their radial distributions along the discs.
In fact, we find a correlation not only among the  parameters of the net accretion and star formation rates,  Net and  SFR , respectively, but also with those of the infall and outflow rates,  Inflows and  Outflows .This is expected, as the combination of infall and outflows induced by feedback linked to star formation determines the net amount of gas which is deposited in the discs and can serve as a fuel for new stellar generations.The correlations between the  values obtained for the various rates is shown in Fig. 9, where we also indicate the identity relation.Note that the infall and outflow rates can only be calculated for the simulations using tracer particles, and therefore the number of data points is smaller than for the full sample.From this figure we can observe that all correlations are close to the identity relation, most notably the one determined by the net accretion rate and the SFR.The most important deviations from the identity relation is found for galaxies with the lowest  values.We find no systematic differences in the behaviour of discs formed in our groups G1 and G2.
The results of this section suggest that the inside-out behaviour of gas accretion, at least for galaxies at the MW-mass scale, is responsible for the inside-out formation of the stellar discs.
In the remainder of this section and the rest of this work, we focus on the galaxies that are compatible with an inside-out behaviour, which represent roughly 70% of our galaxy sample, and study their average behaviour in terms of the infall, outflow, net accretion and star formation rates along the discs.Fig. 10 shows the average massweighted time  as a function of radius for the inflow, outflow, and net accretion rates and for the SFR.The values obtained for the average inside-out parameters , in all three cases, are also shown, together with the standard deviation.The gas inflow rate presents the highest  value, followed by the outflow and the net rates.Note that the number of galaxies for which inflow and outflow rates can be calculated is significantly lower compared to the full sample.
Similarly to our previous analysis, the average mass-weighted time  for all rates can be modelled as a function of the normalised disc radius by a simple linear relation of the form where  *   is the radius on the disc plane normalised by the size of the disc,  is the slope (the inside-out parameter), and  0 the intercept.As long as  is positive (or higher than a given threshold, as we previously commented) the behaviour can be thought of as "inside-out".Table 3 shows the values of the parameters that result from a least squares linear fit of the mass-weighted time as a function of normalised radius, averaged over the G1 galaxies.Note that the values presented in the table are not timescales related to an exponential decay but rather mean times around which the accretion is concentrated at any given normalised disc radius.

Au27
Figure 7. Mass-weighted time  for all galaxies calculated adopting the net accretion rate as a weight, and a linear fit of the data (indicated as a dashed line).The slope of the fit, along with its uncertainty, is shown in each panel and can also be consulted in Table 2.

COMPARISON WITH CHEMICAL EVOLUTION MODELS
Given the average behaviour6 of gas accretion rates in our sample of simulated MW-like galaxies in G1 we would like to compare the  presented in Fig. 10 does not show strong deviations from that of each individual galaxy; therefore, in this section, we assume average functions for a simpler comparison.
results obtained from simulations with those of chemical evolution models (CEMs) which are semi-analytical schemes developed with the aim of providing a description for the distribution of metals in galaxies and, particularly, in the MW.These models employ a series of analytical hypotheses that simulate certain aspects of the process of galaxy formation and evolution such as star formation and, of particular interest in this work, the accretion of gas onto the disc.One of the scenarios that have been proposed in the literature is the so-called "Double Infall Model" (hereafter DIM), built on top of the

Au27
Figure 8. Mass-weighted time  computed adopting the net accretion rate (green, repeated here for comparison) and the star formation rate (purple) as a weight.The purple dashed line shows a linear fit for the SFR curve and we indicate the slope and its uncertainty in each panel.
general assumption that the gas accretion process in disc-like galaxies takes place in two distinct phases (e.g.Chiappini et al. 2001).
In this context, the accretion of gas is given by the sum of two components linked to the formation of the stellar halo, bulge, and thick disc, on one hand, and to the thin disc, on the other.The rate of mass accretion per unit area as a function of cosmic time  and disc radius  is usually expressed as (see e.g.Matteucci 2012): where the first term represents the accretion linked to the spheroidal/thick disc components and the second represents the accretion onto the thin disc.The radial functions () and () are chosen to reproduce the total mass surface density of each stellar component at the present time (see below),  H and  D are typical timescales for the formation of the spheroidal and thin disc components which may depend on , and  max is the time at which the maximum accretion value occurs for the disc.
These models take into account some common assumptions that Table 3. Parameters of the linear fit the mass-weighted time as a function of normalised radius for the average behaviour of G1 (MW-like) galaxies.We indicate the value of the slope  and intercept  0 , along with its uncertainties, for each of the mass rates we analysed.
Mass rate  (Slope)  0 (Intercept) [Gyr] [Gyr] Inflow Rate 3.7 ± 0.2 6.1 ± 0.1 Outflow Rate 2.9 ± 0.1 6.9 ± 0.1 Net Rate 2.3 ± 0.1 7.6 ± 0.1 Star Formation Rate 2.6 ± 0.1 7.6 ± 0.1 are motivated by observational results.In particular, the formation timescale of the spheroidal component's half mass,  H , is assumed to be independent of radius and equal to 0.8 Gyr since the formation of the bulge is believed to occur in a violent, rapid process at early times.Similarly, the time of maximum accretion,  max , is commonly chosen to be either 1 or 2 Gyr and corresponds to the end of the spheroidal/thick disc phase.
The formation timescale of the thin disc,  D , is chosen in consistency with the inside-out formation scenario.This is done by taking the characteristic rate of gas accretion as an increasing function of the disc radius in order to properly reproduce the metallicity distribution of stars in the solar neighbourhood and the bulge.In Chiappini et al. (2001) this dependency is obtained by assuming the following linear function: As already mentioned, the normalisation profiles () and () are chosen to reproduce, respectively, the total mass surface density of the spheroidal and thin disc components at  = 0 or present-day time  G .It can be shown that the first normalisation function can be .Left: Mean mass-weighted time (red), calculated adopting the inflow rates as a weight, for all reruns excluding Au13 and Au17.The shaded region indicates the ±1 region (this is the standard deviation of the data, not of the mean).The grey curves in the background correspond to each individual galaxy considered in the mean.Middle left: Same as the panel on the left but for the outflow rates.Middle right: Same as the panel on the left but for the net accretion rates calculated with the cells.In this case, we consider the fifteen galaxies that show an inside-out formation scheme: Au2, Au3, Au6, Au7, Au9, Au12, Au15, Au16, Au18, Au20, Au21, Au24, Au25, Au26, Au27.Right: Same as the panel on the left but for the star formation rates.We include the same galaxies shown in the net accretion rate panel.dependence with time and radius when radial-and time-integrated values are considered.When focusing on the radially-integrated values, the differences observed in the evolution of the rates are less than one order of magnitude for times greater than ∼ 4 Gyr, after which the discs of the Auriga galaxies can be considered as "well-developed" in terms of the disc-to-total mass ratio.Furthermore, there is considerable agreement on the slope of the temporal decay of rates, which are similar at late times in both cases.On the time-integrated panels, there is a considerably better agreement between the three curves shown.This is especially evident after the first ∼ 4 Gyr of evolution, where all curves are well within the 3 regions shown in each panel.However, some major differences arise at early times and in the inner regions.In the radially-integrated panels, there is a considerable increase in the rates calculated with the fiducial disc timescale of the DIM (purple line), reaching values of order 100 M ⊙ yr −1 .A similar behaviour can be observed in the time-integrated panels, where the accretion rates reach values of order 10 M ⊙ yr −1 near ∼ 2 kpc in the rate case.
Finally, we stress that the sharp discontinuity seen in the timeintegrated panels of Fig. 11 at 2 kpc and, owing to the size of the disc decreasing to almost zero at the beginning of the simulation, as a sharp decrease in the accretion of the model at times smaller than ∼ 2 Gyr.We note, however, that at times less than 4 Gyr the galactic discs are not necessarily formed (see Iza et al. 2022) and, although we considered the inner regions of the galaxy as part of the disc for the calculation, they tend to be dominated by bulges instead of rotationally supported structures (Grand et al. 2017).
This comparison, finally, shows similar values between the accretion rates calculated from simulation and the results obtained from the DIM.Time-and disc-integrated values show similar trends when analysing the radial and temporal dependency, respectively.Most differences between the model and the simulations arise at times or regions when or where the disc may not be completely defined.

CONCLUSIONS
In this work we have investigated the inflow, outflow and net accretion rates onto the discs of MW-mass galaxies, simulated in a full cosmological setting in the context of the ΛCDM model.The simulations are part of the Auriga Project (Grand et al. 2017), and were run with the magnetohydrodynamical, moving-mesh code arepo (Springel 2010).This is the second paper of a series where we analyse the gas flows at the disc-halo interface.In Paper I, we investigated the temporal dependencies of the net accretion rates, integrated onto the radial extent of the discs, as well as for a subset of 9 simulations which included a treatment of tracer particles that allowed us to separately calculate the inflow and outflow rates onto the discs.Here, we focus on the radial dependencies of the inflow, outflow, net accretion and star formation rates onto the disc region, with the aim of verifying whether the usual assumption of inside-out formation of discs for MW-like galaxies is valid in our set of simulations.For this, we excluded from our analysis galaxies whose discs are strongly perturbed during evolution, leaving a total of 24 simulated galaxies, 7 of which have also been run with tracer particles.
More than 70% of the simulated galaxies were found to be consistent with an inside-out formation scenario of the stellar discs, showing an anti-correlation between the age of stars and their radial location along the disc.In order to investigate the origin of this pattern, we calculated the gas flow and star formation rates as a function of radius, and characterised their radial dependencies using the  parameter which measures, for various radial bins, the dominant times of each considered rate.The slope of this relation, denoted as the "inside-out parameter" , allows us to identify systems compatible with an inside-out behaviour and it was computed for the infall and outflow rates ( Inflows and  Outflows , respectively), the net accretion rate ( Net ) and the star formation rate ( SFR ).
A very tight correlation was found between  Net and  SFR , as well as between these two and  Inflows and  Outflows for the simulations with tracer particles.This is expected, as the connections between the star formation and gas flow rates are a natural consequence of the interplay between star formation, the corresponding feedback, and accretion to the discs which determine, at each time and radius, the amount of gas available for star formation.Assuming a conservative choice for the identification of galaxies consistent with an insideout behaviour ( ≥ 0.8 Gyr), we found that the discs of 16 out of the 24 galaxies of the sample follow this trend.Our results show that the inside-out formation of discs is a natural consequence of an analogous accretion process in MW-mass galaxies formed in the context of the ΛCDM model.
For galaxies whose discs do not show an inside-out growth, we found that the relation between accretion/star formation rates and radius is in most cases flat, and, in a few cases, these have negative slopes pointing to an "outside-in" formation law.Most of these galaxies show an active merger activity at late times suggesting that the associated discs did not have enough time to grow in an inside-out manner.The number of galaxies in either of these cases is, however, too small in our galaxy sample to draw any strong conclusion concerning this at this point, but this is certainly worth exploring in future works.
The accretion laws obtained for the simulated galaxies are in relative good agreement with the assumptions usually made in CEMs, in particular in relation to the inside-out growth.Comparisons between the gas accretion rates obtained in our simulations with those expected from classical CEMs for a set of galaxy parameters consistent with our MW-like galaxy sample (G1) show similar slopes for the accretion rate (radially-integrated along the discs) versus cosmic time but with simulated results higher by roughly a factor of 2. At very early times, however, when the galactic discs are not yet formed, differences between CEMs and simulations are larger.The accretion rates (temporally-integrated across galactic history) versus disc radius show more compatible results except in the very inner regions.The relative agreement between the CEM model and the simulations is owing to the fact that the CEM is fed with the mean galactic parameters describing our MW-like galaxy sample, thus validating the general picture of the DIM.However, when comparing simulations to the CEM expectation using the  D () relation adopted in Chiappini et al. (2001) the differences between the model and simulations are larger.This indicates that, although there is always a positive correlation between the accretion timescale  and disc radius, and that the rate of gas flows in the Solar neighbourhood is similar for insideout galaxies in CEMs and simulations, the trends are quantitatively different.Our analysis suggests that galaxy formation simulations of MW-like galaxies in a cosmological context predict a less steep rise in the accretion timescales as one moves farther from the galactic centre.  .Evolution of the disc-to-total mass ratio for the galaxies that do not show a clear inside-out behaviour throughout their history: Au4, Au10, Au11, Au13, Au14, Au17, and Au22.The green (red) bars at the bottom of each panel indicate the time intervals where we obtain an inside-out (not inside-out) behaviour determined by computing the parameter .For some of these galaxies, the inside-out periods coincide with times of disc growth, either after an initial chaotic formation stage (e.g., Au22) or after a (partial) destruction of the disc owing to a merger event (e.g., Au14).

Figure 1 .
Figure1.Projected density maps for the Auriga galaxies at  = 0, for the stellar (left) and gaseous (right) components.For each galaxy we show the face-on view (top) and edge-on view (bottom).The colour maps span five orders of magnitude in projected density using a logarithmic scale.The face-on view shows a cubic region of 100 ckpc on a side and the edge-on view a region of dimensions 100 × 100 × 20 ckpc 3 .

Figure 2 .
Figure 2. As Fig. 1 for the rest of the Auriga galaxies.

Figure 3 .
Figure3.Average stellar age in Gyr as a function of radius on the disc plane (normalised to the disc radius) at  = 0.Each panel shows the profile of a different galaxy and colours indicate the average age over all stars in black, stars with circularity below 0.5 in red (representative of the spheroid) and stars with circularities larger than 0.5 (representative of the disc) in blue, respectively.The vertical dashed-dotted line highlights the radius that encloses 80% of the stellar mass, followingIza et al. (2022).Decreasing profiles indicate an "inside-out" distribution of stars, with older stars in the inner regions and younger stars in the outer regions.

Figure 6 .
Figure 6.Temporal evolution of the gas (inflow-dominated) net rate (using cells) for different radial bins.Each line in this figure represents a different radial ring in the disc plane (the width is 0.1 d for each ring); the first line (labelled as 0.1 d in the legend), for example, represents the ring that extends from    = 0 to    = 0.1 d .The ten lines shown in each panel cover the complete radial extension of the disc for each galaxy.

Figure 9 .
Figure 9. Correlations between the inside-out parameter  calculated with different mass rates (net accretion rate, inflow rate, outflow rate and SFR).Each panel presents a different correlation, indicating the 1:1 line in black.Each point represents a given galaxy, colour-coded according to the group, as shown in the legend.The correlations shown here are close to the identity relation, with most deviations observed for low  values.
Figure10.Left: Mean mass-weighted time (red), calculated adopting the inflow rates as a weight, for all reruns excluding Au13 and Au17.The shaded region indicates the ±1 region (this is the standard deviation of the data, not of the mean).The grey curves in the background correspond to each individual galaxy considered in the mean.Middle left: Same as the panel on the left but for the outflow rates.Middle right: Same as the panel on the left but for the net accretion rates calculated with the cells.In this case, we consider the fifteen galaxies that show an inside-out formation scheme: Au2, Au3, Au6, Au7, Au9, Au12, Au15, Au16, Au18, Au20, Au21, Au24, Au25, Au26, Au27.Right: Same as the panel on the left but for the star formation rates.We include the same galaxies shown in the net accretion rate panel.

Figure 11 .
Figure 11.Comparison between the accretion calculated from the simulations and the Double Infall Model, commonly used in CEMs to model the accretion onto galactic discs.The curves are the results of the simulation (for the inflow-dominated times of the net accretion rate in green and for the inflow rate in blue) and the results of the model (using disc time-scale calculated from simulation data in red and one commonly used in CEMs in purple).We show the ±3 regions for each line in the corresponding colour.The top left panel shows the disc-integrated rate as a function of time, the top right panel shows the disc-integrated flux as a function of time, the bottom left panel shows the time-integrated rate as a function of disc radius and the bottom right panel shows the time-integrated flux as a function of disc radius.The calculation of each of these values is detailed in the text.

Figure A1 .
Figure A1.Correlation between the inside-out parameter  calculated for the net accretion rates and the following galactic properties, from the top left to the bottom right: bulge mass  bulge (calculated as the mass of stars inside 3 ckpc), disc stellar mass  disc , disc-to-total mass fraction D/T, stellar formation time  ★ , virial mass  200 , subhalo stellar mass  ★ , bar mass  bar , and black-hole mass  bh .All the values were calculated at  = 0.The mass of the bar was calculated as  bar =  ★ − 2 ★ (  ≥ 1) − 2 ★ (  ≤ 0), where  is the circularity parameter defined in Paper I. To calculate  ★ we use the time at which half the total stellar mass at  = 0 is accumulated in the galaxy.