apostle–auriga: effects of stellar feedback subgrid models on the evolution of angular momentum in disc galaxies

ABSTRACT

Utilizing the apostle–auriga simulations, which start from the same zoom-in initial conditions of Local Group-like systems, but run with different galaxy formation subgrid models and hydrodynamic solvers, we study the impact of stellar feedback models on the evolution of angular momentum in disc galaxies. At |$z = 0$|⁠, auriga disc galaxies tend to exhibit higher specific angular momenta compared to their cross-matched apostle counterparts. By tracing the evolution history of the Lagrangian mass tracers of the in-situ star particles in the |$z = 0$| galaxies, we find that the specific angular momentum distributions of the gas tracers from the two simulations at the halo accretion time are relatively similar. The present-day angular momentum difference is mainly driven by the physical processes occurring inside dark matter haloes, especially galactic fountains. Due to the different subgrid implementations of stellar feedback processes, auriga galaxies contain a high fraction of gas that has gone through recycled fountain (⁠|${\sim } 65$| per cent) which could acquire angular momentum through mixing with the high angular momentum circumgalactic medium (CGM). In apostle, however, the fraction of gas that has undergone the recycled fountain process is significantly lower (down to |${\sim } 20$| per cent for Milky Way-sized galaxies) and the angular momentum acquisition from the CGM is marginal. As a result, the present-day auriga galaxies overall have higher specific angular momenta.

1 INTRODUCTION

Angular momentum plays a crucial role in shaping many fundamental properties of galaxies, such as size, kinematics, and morphology (e.g. Fall & Efstathiou 1980; Mo, Mao & White 1998; Romanowsky & Fall 2012; Genel et al. 2015; Teklu et al. 2015; Romeo, Agertz & Renaud 2023; Yang et al. 2023). Understanding the evolution of angular momentum is key to comprehending galaxy formation and evolution, especially for the rotation-supported disc galaxies (e.g. Fall & Efstathiou 1980; Mo et al. 1998). In the standard Lambda cold dark matter model, galaxies form through cooling and condensation of gas within their host dark matter haloes (White & Rees 1978; White & Frenk 1991). The initial angular momentum of dark matter haloes is generally acquired through tidal torquing interactions with surrounding inhomogeneities (Peebles 1969; Doroshkevich 1970; White 1984). The evolution of angular momentum of dark matter haloes is primarily driven by gravity and thus can be accurately studied using N-body simulations (e.g. Bullock et al. 2001; Bett et al. 2007; Liao, Chen & Chu 2017b). However, the angular momentum of galaxies is influenced not only by gravity but also by various complex baryonic processes during their evolution (e.g. Genel et al. 2015; Teklu et al. 2015; DeFelippis et al. 2017; Zjupa & Springel 2017; Liao et al. 2017a; Du et al. 2022; Rodriguez-Gomez et al. 2022; Du et al. 2024; Zeng et al. 2024), necessitating the use of hydrodynamic simulations to understand the angular momentum evolution of galaxies.

Early hydrodynamical simulations suffered from the so-called angular momentum catastrophe, where gas experienced over-cooling and drastic angular momentum loss, leading to the formation of compact clumps instead of an extended disc (e.g. Navarro, Frenk & White 1995; Navarro & Steinmetz 2000; Maller & Dekel 2002). Baryonic feedback is considered a necessary heating mechanism to regulate the star formation and prevent over-cooling (see e.g. Somerville & Davé 2015; Naab & Ostriker 2017; Crain & van de Voort 2023, for reviews). However, the limited resolution and the wide range of spatial and temporal dynamic scales in cosmological hydrodynamical simulations make it challenging to include ab initio galaxy formation processes based on first principles. As an alternative method, most baryonic processes have to be modelled using an effective parametrized subgrid approach. The effective stellar feedback is a key ingredient of subgrid physical models for reproducing disc galaxies. In general, the stellar feedback-driven outflows can remove low angular momentum gas, preventing the formation of dense spherical compositions in the early stages of galaxy formation (e.g. Scannapieco et al. 2008, 2012; Übler et al. 2014; Agertz & Kravtsov 2016). With the increased resolution and the adoption of more sophisticated subgrid models, modern hydrodynamical simulations can reproduce more extended disc galaxies (e.g. Governato et al. 2010; Teklu et al. 2015). In addition, baryonic feedback processes are crucial for shaping the observed properties of both circumgalactic (CGM) and intergalactic media around galaxies (e.g. Nelson et al. 2015; Kelly et al. 2022; Faucher-Giguère & Oh 2023).

As the representatives of the current state-of-the-art cosmological hydrodynamical simulations, both the eagle family (Crain et al. 2015; Schaye et al. 2015) and the illustris family [including the illustris (Vogelsberger et al. 2013, 2014), auriga (Grand et al. 2017) and illustrisTNG (Nelson et al. 2018; Pillepich et al. 2018) simulations]¹ effectively reproduce a large body of observed galaxy properties. For instance, they successfully reproduce diverse morphological types (Snyder et al. 2015; Correa et al. 2017; Tacchella et al. 2019; Rodriguez-Gomez et al. 2019), the global stellar mass assemble history (Furlong et al. 2015; Sparre et al. 2015; Qu et al. 2017), the colour bimodality distribution (Vogelsberger et al. 2014; Trayford et al. 2015), the galactic size-to-mass relation (Furlong et al. 2017; Genel et al. 2018; Yang et al. 2023), and so forth.

However, due to the adoption of different subgrid physics models and hydrodynamics solvers, there are discrepancies in galaxy properties between the eagle and illustris families. In particular, Yang et al. (2023) recently revealed significant differences in the correlation between disc galaxy sizes and halo spins in the two simulation families. For Milky Way-sized disc galaxies, both the eagle and illustris families exhibit positive correlations, but the correlation in the former is weaker. Conversely, in the case of dwarf disc galaxies, no correlation is observed in the eagle family, while a positive correlation is identified in the illustris family. This finding raises questions about the validity of classical semi-analytical disc formation models (e.g. Mo et al. 1998; Guo et al. 2011), which predict the positive correlation between disc size and halo spin. It becomes crucial to investigate the reasons behind the distinct stellar disc evolution (especially the evolution of angular momentum) resulting from the two types of subgrid models in these simulations.

In this work, we further investigate the formation histories of disc galaxies simulated by the eagle and illustris families, with a primary focus on examining the impact of different stellar feedback implementations on the angular momentum evolution of these simulated disc galaxies. We utilize two suites of high-resolution zoom-in simulations, apostle (Sawala et al. 2016; Fattahi et al. 2016) and auriga (Grand et al. 2017), derived from the eagle and illustris families, respectively. Specifically, for the auriga suite, we adopt the simulations started from the same initial conditions (ICs) as the apostle AP-V1 and AP-S5 runs but run with the auriga model (Kelly et al. 2022). This approach enables a direct matching of haloes/galaxies between the two families, facilitating an apples-to-apples comparison.

By selecting corresponding disc galaxies within the matched halo pairs and employing the mass element tracer technique, we perform a detailed particle-based Lagrangian analysis to compare the angular momentum evolution between the two simulation families.

This paper is organized as follows. In Section 2, we provide the details of the two simulations and the galaxy sample selection. The method to trace the evolution of particles’ angular momenta is described in Section 3. We present our results in Section 4 and summarize in Section 5.

2 SIMULATIONS AND GALAXY SAMPLES

2.1 Simulation overview

This work makes use of two suites of zoom-in hydrodynamical simulations representing Local Group-like volumes. The first set comprises two volumes, namely AP-V1 and AP-S5, selected from the apostle project (Sawala et al. 2016; Fattahi et al. 2016). The apostle Local Group-like regions are drawn from the dark matter-only simulation dove according to their mass and kinematics at present. For detailed information on the selection conditions, we direct the interested reader to Fattahi et al. (2016). The zoom-in ICs are created by ic_gen (Jenkins 2010) which employs the second-order Lagrangian perturbation theory.

The apostle simulations are run using a highly modified version of the gadget-3 code, which is an improved version of gadget-2 (Springel 2005). The gas properties are calculated with the pressure–entropy formulation of smoothed particle hydrodynamics (SPH) method (Hopkins 2013; Schaller et al. 2015). All the numerical set-up and parameters are the same as the eagle reference model (Crain et al. 2015; Schaye et al. 2015).

The second suite of simulations use the same ICs as the apostle suite, but is run with the auriga galaxy formation model (Grand et al. 2017). The auriga simulations are performed with the moving-mesh magneto-hydrodynamics code arepo (Springel 2010). For convenience, we still call this suite of simulations as auriga, but note that it adopts different ICs compared to the original auriga project (Grand et al. 2017).

Both two simulation sets adopt the WMAP-7 (Wilkinson Microwave Anisotropy Probe) cosmological parameters (Komatsu et al. 2011) which are |$\Omega _{\Lambda }=0.728$|⁠, |$\Omega _{\rm m}=0.272$|⁠, |$\Omega _{\rm b}=0.0455$|⁠, |$h=0.704$|⁠, |$\sigma _{8}=0.81$|⁠, and |$n_{\rm s}=0.967$|⁠. The initial gas (dark matter) particle/cell mass is about |$1.2(5.9)\times 10^5\, {\rm M}_{\odot }$|⁠, and the maximum softening length is 307 pc in both simulations. The friends-of-friends (FOF, Davis et al. 1985) and subfind (Springel et al. 2001; Dolag et al. 2009) algorithms are applied to identify halo and subhalo structures, respectively. In this work, the halo virial radius, |$R_{200}$|⁠, is defined as the radius within which the mean density is 200 times the critical density of the Universe. Unless otherwise specified, we use the star particles within the 3D aperture with a radius of |$0.1 R_{200}$| to compute the stellar properties. Both simulations contain 128 particle snapshots and corresponding group catalogs, with a time interval being about 100 Myr. We match the dark matter particles within subhaloes by particle IDs to identify the main progenitor branch of the merger tree across different snapshots. Specifically, for a subhalo in snapshot |$S_{n}$|⁠, we consider its dark matter particles within twice the half-total-mass radius, and define the subhalo in the previous snapshot (i.e. snapshot |$S_{n-1}$|⁠) that contains the most matched particles as the main progenitor.

Both apostle and auriga simulations use the treepm method to compute gravity. However, they adopt different hydrodynamical solvers (i.e. SPH for the former versus moving-mesh for the latter), which can potentially introduce differences in the baryon cycle (e.g. Kereš et al. 2012; Nelson et al. 2013). Nevertheless, we expect that the effects of hydrodynamical solvers are secondary compared to those of the baryonic subgrid models, which are elaborated in the next subsection, especially for the evolution of gas within the halo virial radius (see Scannapieco et al. 2012; Hayward et al. 2014; Schaller et al. 2015; Hopkins et al. 2018). As we will show in Section 4.2, the distributions of gas angular momenta exhibit noticeable similarities in both simulations when the gas enters the halo virial radius, with larger differences emerging gradually thereafter.

2.2 Subgrid models

Both apostle and auriga simulations incorporate key baryonic processes relevant to galaxy formation, such as radiative cooling (Wiersma, Schaye & Smith 2009a), star formation (Schaye & Dalla Vecchia 2008), stellar evolution and metal enrichment (Wiersma et al. 2009b), stellar feedback (Dalla Vecchia & Schaye 2012; Vogelsberger et al. 2013), and black hole mass growth and active galactic nuclei (AGNs) feedback (Springel, Di Matteo & Hernquist 2005; Booth & Schaye 2009; Rosas-Guevara et al. 2015). The detailed descriptions of the subgrid models used in apostle and auriga can be found in Schaye et al. (2015) and Grand et al. (2017), respectively. In this subsection, we mainly summarize the differences in the implementation of stellar feedback, which plays a significant role in influencing the evolution of baryonic angular momentum in the dwarf and Milky Way-sized galaxies studied in this work.

In auriga, a gas cell is assumed to be star-forming and thermally unstable when its density is higher than the threshold density of |$n_{\rm th} = 0.13~\rm {cm^ {-3}}$|⁠. Once entering the thermally unstable star-forming regime, this gas cell is stochastically converted into a star particle according to a probability that scales exponentially with time in units of the star-forming time-scale |$t_{\rm sf} = 2.2~{\rm Gyr}$|⁠. In apostle, the star formation is also implemented stochastically, but following the pressure-dependent Kennicutt–Schmidt law (Schaye & Dalla Vecchia 2008) with a metallicity-dependent density threshold value |$n_{\rm th} = 0.1 (Z/0.002)^{-0.64} \rm {cm^{-3}}$|⁠, where Z is the gas metallicity (i.e. the mass fraction of gas in elements heavier than helium). Although the density threshold has some differences, the star formation processes are mainly self-regulated by feedback (e.g. Schaye et al. 2010).

In both apostle and auriga, a star particle represents a single stellar population (SSP), and the Chabrier (2003) initial mass function is adopted to specify the distribution of stellar masses contained in an SSP. Hence, the amount of supernovae (SNe) energy per unit stellar mass released into surrounding gas is similar in two simulations. However, the SNe feedback implementation in the two simulations are different.

In auriga, the kinetic SNe feedback is implemented by converting a gas cell into a wind particle (Vogelsberger et al. 2013). The wind particle is kicked with an isotropically random direction, and the wind speed |$v_{\rm {wind}}$| is proportional to the local 1D dark matter velocity dispersion |$\sigma _{\rm {dm}}$| which is calculated from 64 nearest dark matter particles, |$v_{\rm {wind}} = 3.46\sigma _{\rm {dm}}$|⁠. After launch, a wind particle is decoupled from the hydrodynamical computation until either reaching a low-density region (5 per cent of star-forming threshold density) or exceeding a maximum time (2.5 per cent of the Hubble time at launch). Once the wind particle recouples, it deposits its mass, metals, momentum, and energy into the intersected gas cells.

In contrast, the apostle simulations adopt the thermal implementation of stellar feedback (Dalla Vecchia & Schaye 2012). Specifically, the SNe energy of an SSP is stochastically ejected to a few neighbouring particles within the SPH kernel size, resulting in an increase in the internal energy of these gas particles. Note that in apostle, the stellar feedback is performed only through the thermal mechanism. To avoid that the cooling time is shorter than the sound cross time, each heated gas particle is subject to the same high enough temperature increase, |$\Delta T_{\rm SN} = 10^{7.5}$| K. As demonstrated by Dalla Vecchia & Schaye (2012), such a sufficiently high temperature boost is required for the injected thermal energy to be efficiently converted into kinetic energy, leading to massive, large-scale outflows.

2.3 Galaxy samples

In this work, we focus our analysis on central disc galaxies, which are classified by applying a criterion based on the kinematics-based morphology parameter |$\kappa \gt 0.5$| (Yang et al. 2023). Here, following Sales et al. (2012), the |$\kappa$| parameter is defined as the ratio between the rotational energy, |$K_{\rm rot}$|⁠, and the total kinetic energy, K, i.e.

where |$\hat{\boldsymbol{L}}$| is the unit vector of the total stellar angular momentum, |$\hat{\boldsymbol{r}}_i$| is the unit vector of the position vector from the galaxy centre for star particle i, |$\boldsymbol{v}_i$| is its velocity vector in the centre of mass frame, and |$m_i$| is its particle mass. For each galaxy at |$z=0$|⁠, |$\kappa$| is computed using all star particles within |$0.1 R_{200}$|⁠.

In total, we identify 12 isolated disc galaxies from the high-resolution regions of two auriga runs. To ensure reliable resolution, we select only the four most massive isolated disc galaxies from each auriga run (eight galaxies in total), each containing at least 50 000 star particles. The subhalo IDs, virial radii (⁠|$R_{200}$|⁠), and virial masses (⁠|$M_{200}$|⁠) of the eight selected auriga galaxies are provided in columns 2–4 of Table 1. The galaxies from the two volumes are denoted as S5-1 to S5-4 and V1-1 to V1-4, respectively.

To offer an apples-to-apples comparison, we further match the eight |$z=0$| counterparts from the apostle runs by comparing the virial radii and galaxy centres (i.e. the position of the particle with the minimum gravitational potential energy, |$\boldsymbol {r}_{\rm min}$|⁠). Specifically, for each auriga galaxy, we loop over all apostle galaxies and compute the relative differences,

and

The matched apostle galaxy is defined as the one with the minimum value of |$\Delta =\sqrt{\Delta _R^2+\Delta _r^2}$|⁠. We have further verified the cross-matching results by comparing the dark matter particle IDs from the matched galaxy pairs. The FOF group IDs, the subfind group IDs, the virial radii, and the virial masses of the eight apostle galaxies can be found in columns 5–8 of Table 1.

To distinguish the auriga and apostle galaxies, we introduce the prefixes ‘Au-’ and ‘Ap-’ to the galaxy names, respectively. For example, ‘Au-S5-1’ and ‘Ap-S5-1’ represent the most massive galaxy pairs within the S5 volume in the two respective simulations. The baryonic properties of the galaxy pairs are summarized in Table 2.

The auriga galaxies tend to have slightly higher |$M_{200}$| (⁠|$R_{200}$|⁠) compared to the apostle counterparts, with the relative differences being |$\lesssim 20~{{\ \rm per\ cent}}$| (⁠|$\lesssim 10~{{\ \rm per\ cent}}$|⁠), as a result of the higher baryonic fractions in auriga (see Table 2).

In this study, we analyse a total of eight pairs of disc galaxies spanning from dwarf to Milky Way-sized systems, i.e. with their present-day virial masses |$M_{200}$| ranging from |$6.1 \times 10^{10}$| to |$1.6 \times 10^{12} {\rm M}_{\odot }$| and virial radii raning from |${\sim } 80$| to |${\sim } 250$| kpc. The primary motivation of this study is to examine the influence of different stellar feedback implementations on the evolution of angular momentum. For Milky Way-sized galaxies, the comparison can be partially contaminated by the presence of AGN feedback. By incorporating a sample of dwarf galaxies, which do not harbour supermassive black holes, we can disentangle the impact of AGN feedback from our overall conclusions.

To visually inspect the stellar components of the selected galaxies, we plot the face-on and edge-on projected maps of the matched galaxy pair from the volumes S5 and V1 in Figs 1 and 2, respectively. The face-on direction is determined by the direction of the stellar specific angular momentum (sAM) computed using star particles within |$2 r_{*, 1/2}$|⁠.² As illustrated in Figs 1 and 2, all samples exhibit distinct disc structures. In auriga, galaxies display a higher central stellar surface density and a more complex dynamic structure, including bars and spiral arms. This phenomenon may originate from the non-axisymmetric instability as a result of the higher baryonic mass to dark matter mass ratio (Ostriker & Peebles 1973; Efstathiou, Lake & Negroponte 1982) in auriga compared to apostle.

3 ANGULAR MOMENTUM TRACING METHOD

The Lagrangian point of view has been extensively utilized in previous simulation studies to examine the evolution of galaxy angular momentum (e.g. White 1984; Sugerman, Summers & Kamionkowski 2000; Porciani, Dekel & Hoffman 2002; DeFelippis et al. 2017; Liao et al. 2017a; Sheng et al. 2023). This approach involves tracing the mass elements comprising a |$z=0$| structure back to earlier snapshots and computing their angular momentum to construct the evolutionary history. It has proven to be insightful in understanding the origin and evolution of galaxy angular momentum, especially in modern galaxy formation simulations (e.g. DeFelippis et al. 2017). In this study, we also employ this approach, and the detailed methodology is described as follows.

3.1 Particle tracing

To understand how the angular momentum of the stellar disc evolves into the final state, we trace back the history of each star particles within |$0.1 R_{200}$| in the galaxy at |$z=0$|⁠. In apostle, the Lagrangian SPH solver is used, and each star particle is converted from a corresponding progenitor gas particle. Therefore, it is natural to directly track a mass element from its very beginning as a gas particle to its final state as a star particle by using its unique particle ID.

However, in auriga, the gas ‘particles’ correspond to mesh cells rather than actual Lagrangian particles in SPH. Although arepo uses a quasi-Lagrangian moving mesh scheme that minimizes mass fluxes between cells, there is still non-zero mass exchange between cells. Furthermore, the number of cells is not conserved in order to adjust the mass resolution when the refinement (de-refinement) operation is used to split (merge) gas cells. Therefore, each gas cell might not represent the same material elements during the simulation. The additional tracer particle must be adopted to track the mass flow. The auriga simulations include Monte Carlo tracer particles (Genel et al. 2013) to track the mass flow throughout the simulation in a Lagrangian way. Each gas cell contains one tracer particle at the beginning of the simulation. According to mass exchange or mesh refinement (de-refinement), tracer particles are then exchanged between neighbouring cells in a probabilistic way. Tracer particles can also be retained when their host gas cells have been transferred to star particles. In this paper, we select all tracer particles resided in the selected stellar composition of the galaxy at |$z=0$| in auriga, and then record the information of their corresponding host particles at each snapshot. For convenience, we collectively refer to both SPH gas particles in apostle and Monte Carlo tracers in auriga as ‘gas tracers’ in the following paragraphs.

3.2 Characterizing the evolution of angular momentum

Once we have the evolutionary history of each |$z=0$| star particle within |$0.1 R_{200}$|⁠, we know whether it has an in-situ (formed within the main progenitor branch of the merger tree) or ex-situ (formed outside) origin. Since our focus is on how the baryonic feedback from the target galaxy affects its angular momentum evolution, we analyse only the in-situ star particles and disregard those originating from ex-situ accretion. Table 2 presents the ex-situ stellar mass fractions, |$f_{\rm {ex}}$|⁠, in the rightmost column. Notably, for our galaxy samples, the majority of star particles have an in-situ origin (i.e. |$f_{\rm ex} \lesssim 15~{{\ \rm per\ cent}}$|⁠). Note that similar considerations are also adopted in DeFelippis et al. (2017) when studying the impact of galactic winds on stellar angular momentum in the illustris simulations and their ex-situ mass fractions are similar, i.e. |${\sim } 10~{{\ \rm per\ cent}}$|⁠.

For each tracing particle i, we compute its sAM at each snapshot,

where |$\boldsymbol {r}_i$| and |$\boldsymbol {v}_i$| are the position and velocity of particle i, |$\boldsymbol {r}_{\rm min}$| is the galaxy centre, and |$\boldsymbol {v}_{\rm c}$| is the centre of mass velocity computed from all particles in the galaxy (including dark matter, gas, stars, and black holes if available). The z-axis component of sAM, |$j_{i,z}$|⁠, is computed by projecting |$\boldsymbol {j}_{i}$| onto the face-on direction defined at |$z = 0$| (see Section 2.3 for details).

The gas tracers we analysed in this paper eventually fall into the gas disc and form stars. However, after entering the disc and crossing the star formation threshold for the first time, not all gas remains in this state and subsequently converts into star particles. Some gas is ejected from the high-density ISM region, temporarily leaves the star-forming state, and later returns. This phenomenon is known as galactic fountains (see Fraternali 2017, for a review). This motivates us to compare the sAM of particles from the apostle and auriga simulations at the following characteristic time (DeFelippis et al. 2017), as we quantify the evolution of galaxy angular momentum:

• Halo accretion time (⁠|$z_{\rm {acc}}$|⁠): the tracer belongs to the main progenitor of the halo for the first time.

• First star-forming time (⁠|$z_{\rm {fsf}}$|⁠): the gas element enters the central galaxy disc (the distance to the galaxy centre is less than |$0.1R_{200}$|⁠) and crosses the star-forming density threshold for the first time.

• Last star-forming time (⁠|$z_{\rm {lsf}}$|⁠): the gas element crosses the star-forming density threshold for the last time. Note that for gas tracers that have never exited the star-forming state before being converted into stars, their |$z_{\rm {fsf}}$| and |$z_{\rm {lsf}}$| are identical.

• Star formation time (⁠|$z_{\rm {sf}}$|⁠): the tracer is converted from the gas phase into the stellar phase.

• The present time (⁠|$z = 0$|⁠): the final snapshot of simulations.

4 RESULTS

4.1 Angular momentum properties at |$z=0$|

To compare the stellar angular momentum properties at |$z = 0$| and check if both simulations agree with observations, in Fig. 3, we plot the |$z = 0$| stellar sAM–mass (⁠|$j_{*}$|–|$M_{*}$|⁠) relation of the simulated galaxies together with the observations from Romanowsky & Fall (2012). Overall, both auriga and apostle simulations reasonably reproduce the observed relation. This is consistent with the recent results of Hardwick et al. (2023), who found that the stellar |$j_{*}$|–|$M_{*}$| relations from both eagle and illustrisTNG simulations agree well with the xGASS (extended GALEX Arecibo SDSS Survey) observations (Hardwick et al. 2022). It also echos previous literature that studies the |$j_{*}$|–|$M_{*}$| relation in either eagle (e.g. Lagos et al. 2017) or illustrisTNG (e.g. Du et al. 2022; Rodriguez-Gomez et al. 2022; Fall & Rodriguez-Gomez 2023) simulations.

We also notice that auriga galaxies tend to exhibit marginally higher stellar mass and sAM compared to their apostle counterparts. This trend aligns with the slightly higher rotation-to-kinetic energy ratios observed in auriga galaxies (see the measurement results of |$\kappa$| parameters in Table 2). What causes the auriga galaxies to show relatively higher stellar angular momentum at |$z=0$|? Does this difference already appear in the initial gas accretion, or does it gradually develop due to later baryonic processes occurring inside haloes? To answer these questions, in the next subsection, we trace the Lagrangian tracers back to |$z = z_{\rm acc}$| and compare the sAM distributions in the two simulations.

4.2 Distribution of trace-back angular momentum at |$z_{\rm acc}$|

We track the in-situ star particles within |$0.1 R_{200}$| of the |$z=0$| galaxies back to |$z = z_{\rm acc}$|⁠, the time when the gas tracers first enter the main halo progenitor, and plot the probability distribution functions (PDFs) of the magnitude of their sAM, |$\log _{10} j$|⁠, as filled histograms in Fig. 4. The median sAM of each distribution is marked using a vertical solid line. Remarkably, the distributions from the the two simulations are relatively similar at |$z = z_{\rm {acc}}$|⁠. This outcome arises from the similar large-scale tidal torquing, as these simulations start from the same ICs, which determines the evolution of angular momentum in the early stage (e.g. Peebles 1969; Doroshkevich 1970; White 1984; Danovich et al. 2015).

As a comparison, we plot the j distributions at the endpoint of the gas evolution, |$z = z_{\rm sf}$| (i.e. the moment when the gas tracers are converted into star particles) in Fig. 4 (unfilled histograms). Once again, the median sAM of each distribution is denoted by a vertical dashed line. In contrast to the distributions at |$z = z_{\rm acc}$|⁠, notable differences emerge between the apostle and auriga simulations, with auriga tracers exhibiting higher j values. From |$z_{\rm acc}$| to |$z_{\rm sf}$|⁠, the average sAM loss amounts to approximately 0.70 (0.97) dex in the auriga (apostle) simulations.

This comparison hints that the difference in the final stellar sAM between the two simulations mainly originates from the processes within haloes, rather than from initial accretion. After being accreted into the main halo progenitor and before being converted into stars, gas can experience substantial changes in angular momentum due to the energy and momentum received from feedback processes. The observed disparity in sAM distributions at |$z = z_{\rm sf}$| suggests that different subgrid models for these physical processes may lead to distinct sAM evolution in these simulations.

4.3 Distributions of characteristic time

We employ some characteristic events which defined in Section 3.2 to delineate the entire evolution history. To quantify the typical values for the characteristic time, we plot their cumulative distribution functions (CDFs) in Fig. 5. The auriga and apostle runs are distinguished by solid and dashed lines, respectively. The median time for each CDF is marked by a vertical line segment.

For example, in the S5-1 galaxy pair, the distributions of halo accretion time (black) are relatively similar in the two simulations, with the median halo accretion time being |$t_{\rm {acc}}^{\rm Au} = 3.82$| Gyr and |$t_{\rm {acc}}^{\rm Ap} = 4.25$| Gyr (where t denotes the cosmic age) for auriga and apostle. For the first star-forming time (red), the median values for auriga and apostle being |$t_{\rm fsf}^{\rm Au} = 5.38$| Gyr and |$t_{\rm fsf}^{\rm Ap} = 5.57$| Gyr, respectively. However, the distributions for the last star-forming time (blue) and the star formation time (green) show significant differences between the two simulations. In auriga, the median last star-forming time is |$t_{\rm lsf}^{\rm Au} = 7.74$| Gyr and the median star formation time is |$t_{\rm sf}^{\rm Au} = 8.56$| Gyr. In contrast, in apostle, these median values are |$t_{\rm lsf}^{\rm Ap} = 5.95$| Gyr and |$t_{\rm sf}^{\rm Ap} = 6.89$| Gyr.

Similar trends are observed for other galaxy pairs. The substantial differences in the CDFs of |$z_{\rm lsf}$| and |$z_{\rm sf}$| between the two simulations indicate that the evolution of gas tracers after |$z = z_{\rm fsf}$| differs markedly in auriga and apostle, with gas taking more time to form stars after crossing the star-forming threshold in auriga. As we will detail in Section 4.5, this difference relates to the distinct behaviours of recycled fountain flows in these two simulations.

4.4 Evolution of angular momentum from |$z=z_{\rm acc}$| to |$0$|

To study the entire evolution history of angular momentum, we compute the sAM of each tracer at different characteristic events defined in Section 3.2. As we have shown in Section 4.2, the distributions of sAM at |$z = z_{\rm acc}$| from the two simulations are relatively similar. Therefore, in this subsection, we focus on the sAM change relative to the sAM at |$z = z_{\rm acc}$|⁠.

Specifically, for each tracer in a galaxy, we compute the sAM loss with respect to its |$j_i$| at |$z = z_{\rm acc}$|⁠, i.e. |$\Delta \log _{10} j_i (z) = \log _{10} j_i (z)-\log _{10} j_i (z_{\rm acc})$|⁠. The median sAM loss of this galaxy is then defined as the median |$\Delta \log _{10} j (z)$| of all tracers. In Fig. 6, the median sAM loss computed over eight galaxies is plotted using solid connected lines in the top panel, with the error bar representing the 16th–84th percentiles. In addition, to examine the evolution of the z-axis component of sAM (i.e. the component along the direction of the stellar sAM at the present day |$z = 0$|⁠), we plot the median relative loss of |$j_z$| normalized to |$|j_z|$| at the halo accretion event (solid connected lines) in the bottom panel.

We can observe the following evolution process of angular momentum from Fig. 6:

From |$z=z_{\rm acc}$| to |$z_{\rm fsf}$|. After being accreted into dark matter haloes, each gas tracers in auriga (apostle) on average lose approximately 0.71 (0.82) dex of their sAM magnitudes when entering the central galactic disc and crossing the star-forming density threshold for the first time (i.e. at |$z = z_{\rm fsf}$|⁠). During this evolution, tracers lose |${\sim }16 ({\sim }64)$| per cent of |$j_z$| in auriga (apostle). As pointed out by DeFelippis et al. (2017), the decrease in these tracers’ sAM results from the transport of angular momentum to other components (e.g. dark matter and baryonic) and the cancellation with other baryons which have been accreted with different sAM directions. The larger sAM loss in the apostle simulations implies that the angular momentum transfer and cancellation are more effective in the eagle model.

From |$z = z_{\rm fsf}$| to |$z_{\rm lsf}$|. The average sAM magnitude of gas tracers in auriga increases by |${\sim } 0.12$| dex, and the average increase in |$j_z$| is |${\sim } 60$| per cent. In contrast, the average sAM magnitude in apostle remains almost constant. As we will detail in Section 4.5, this process is related to the so-called galactic fountain, where gas flows are ejected by stellar feedback and recycled back, acquiring angular momentum through mixing with CGM and newly accreted gas (e.g. Brook et al. 2012; Grand et al. 2019). The almost constant sAM in this phase in apostle reflects that its recycled fountain is significantly weaker compared to auriga.

From |$z = z_{\rm lsf}$| to |$z_{\rm sf}$|. The average sAM magnitudes of gas tracers in both simulations decrease by |${\sim } 0.1$| dex, while the average |$j_z$| is almost unchanged. This indicates that during this stage, as the gas tracers continue to cool, condense into the galactic disc, and finally turn into stars, the sAM components perpendicular to the stellar disc direction are gradually dissipated into the interstellar medium (ISM), leaving only the sAM component along the disc direction.

From |$z = z_{\rm sf}$| to |$0$|. The stellar sAM are almost conserved in both simulations.

To summarize, by tracing the evolution history of the star particles in a |$z = 0$| disc galaxy, we find that the sAM difference at |$z = 0$| between the apostle and auriga simulations mainly develops between the halo accretion and last star-forming phases, which is closely related to the galactic fountain. In the following subsection, we delve into more details of the effects of fountain flows.

4.5 Effect of recycled fountain flows

Previous studies have shown that after being ejected from the ISM, the fountain flows mix with the low-density, high-temperature CGM gas, which usually has high angular momentum. As a result, galactic fountains effectively redistribute angular momentum, with the recycling gas acquiring angular momentum from the CGM (e.g. Brook et al. 2012; Valentini et al. 2017; Grand et al. 2019). This mechanism can explain the increase of sAM for gas tracers between |$z = z_{\rm fsf}$| and |$z_{\rm lsf}$| in simulations.

To examine the above scenario, we quantify the fraction of recycled flows in the gas tracers. Specifically, we record the number of times each gas tracer leaves star-forming state after |$z = z_{\rm fsf}$|⁠. Based on whether they leave the star-forming state at least once, gas tracers are divided into non-recycled and recycled categories.³ The number fractions of recycled and non-recycled tracers for each galaxy pair are summarized in Fig. 7. Note that in this figure, the galaxy pairs are plotted from left to right in descending order by virial mass.

In apostle, the majority of tracers are categorized as non-recycled gas, with the fractions of non-recycled gas being higher in more massive galaxies. Specifically, the fraction of non-recycled gas is |${\sim } 80$| per cent for Milky Way-sized galaxies, while it is |${\sim } 50$| per cent for dwarf galaxies. In contrast, the majority of tracers in auriga are recycled gas, with the fraction of non-recycled (recycled) gas being |${\sim } 35$||$({\sim } 65)$| per cent for all galaxies. The fractions of non-recycled or recycled gas do not show a clear mass dependence in auriga.

This comparison highlights the differences in galactic fountains between the two simulations, which are a direct outcome of their different stellar feedback subgrid models. Previous studies (e.g. Kelly et al. 2022; Wright et al. 2024) have shown that in galaxies below |${\sim } 10^{12}~{\rm M}_{\odot }$|⁠, the eagle stellar feedback model tends to eject gas particles from the ISM to several virial radii, reducing recycling flows within the virial radius and impeding cosmic gas accretion. In contrast, in the auriga simulations, the ejected gas flows driven by stellar feedback typically stall before reaching the galaxy’s virial radius, resulting in a higher degree of recycled fountain flows. As we focus on gas that eventually forms stars, in apostle, this gas usually does not experience strong ejections from stellar feedback. In auriga, however, even gas ejected by stellar feedback can recycle back to the galactic disc and form stars.

To illustrate the recycled flow, we randomly select 10 gas tracers with the same first star-forming time, |$z_{\rm fsf}$|⁠, from Galaxy S5-1 in both simulations. The time evolution of their distances to the galaxy centre (rescaled by the virial radius, i.e. |$r/R_{200}$|⁠) and sAM before being converted into stars are plotted in Fig. 8. In Au-S5-1, the majority of the selected gas tracers (6 over 10) are ejected to the CGM beyond |$0.1 R_{200}$| and then recycle back to form stars. These outflow-recycled gas tracers, which have low angular momenta at |$z = z_{\rm fsf}$| (agreeing with Brook et al. 2011), acquire angular momenta through mixing with the high angular momentum CGM gas, eventually settling in the relatively outer disc. This agrees with the findings of previous studies (e.g. Brook et al. 2012; Grand et al. 2019). In Ap-S5-1, the majority of the selected gas tracers form stars quickly after |$z = z_{\rm fsf}$|⁠, with only 2 over 10 being weakly ejected outside the galactic disc and forming stars shortly after recycling back. In addition, we notice that after |$z_{\rm fsf}$|⁠, the recycled gas tracers in auriga take more time to evolve into stars compared to those in apostle, providing an explanation for the observed differences in the CDFs of characteristic time in Section 4.3.

Finally, to investigate directly the effect of recycled fountain flows on angular momentum evolution, we consider only the recycled gas tracers and plot their angular momentum evolution in Fig. 6 using dashed lines. We can clearly see that between |$z = z_{\rm fsf}$| and |$z_{\rm lsf}$|⁠, the increase in sAM in auriga is larger compared to the ‘all tracers’ case. Intriguingly, we also observe an increase in sAM in apostle, reflecting that the recycled flows can indeed boost the angular momentum regardless of simulations. The almost constant sAM in the ‘all tracers’ case of apostle is a result of the dominant fraction of non-recycled tracers.

5 SUMMARY

In this paper, we have studied the impact of stellar feedback subgrid models on the evolution of disc galaxies’ angular momenta by using two suits of zoom-in simulations, apostle and auriga. By tracing the evolution history of the |$z = 0$| star particles from a Lagrangian perspective, we characterize the evolution of the particle tracers’ angular momenta at five different events: halo accretion, first star-forming, last star-forming, star formation, and the present state. We summarize the following:

• At |$z = 0$|⁠, both the apostle and auriga disc galaxies closely follow the observed stellar angular momentum–mass relation. However, the auriga galaxies exhibit slightly higher sAM and higher stellar masses compared to their apostle counterparts (Fig. 3).

• By tracing the Lagrangian tracers of the in-situ star particles in |$z = 0$| galaxies back to the halo accretion time, we find that the sAM distributions of the two simulations are relatively similar. This suggests that the main driver of the present-day difference in sAM is not initial accretion but rather the later processes occurring inside haloes (Fig. 4).

• After being accreted into dark matter haloes (⁠|$z = z_{\rm acc}$|⁠), gas tracers in auriga (apostle) on average lose |${\sim } 0.71({\sim } 0.82)$| dex of their sAM when entering the central galactic disc and reaching the star-forming condition for the first time (⁠|$z = z_{\rm fsf}$|⁠). From |$z = z_{\rm fsf}$| to |$z_{\rm lsf}$|⁠, the average sAM of gas tracers in auriga increases by |${\sim }0.12$| dex, while it remains almost constant in apostle. In the remaining evolution stages (i.e. last star-forming–star formation–present-day), the average angular momentum evolution of the Lagrangian tracers is similar in the two simulations (Fig. 6).

• Different implementations of stellar feedback in the two simulations lead to distinct recycled fountain flows, which further affect the angular momentum evolution history. auriga galaxies contain a higher fraction (on average |${\sim } 65$| per cent) of recycled gas tracers, and this gas acquires angular momentum through mixing with the high angular momentum CGM. In contrast, the Lagrangian tracers of the |$z=0$|apostle stellar discs are dominated by non-recycled gas (i.e. up to |${\sim } 80$| per cent for Milky Way-sized galaxies) and the acquisition of angular momentum from the CGM is negligible (Figs 7 and 8). As a result, auriga disc galaxies exhibit higher angular momentum at |$z = 0$| compared to apostle galaxies.

Our results reveal the different angular momentum evolution histories of disc galaxies in the apostle and auriga galaxy formation subgrid models. Future investigations into the impact of baryonic physics subgrid models (such as stellar feedback) on the masses and sizes of stellar discs will be valuable, as they will further help to understand the origin of differences in the disc size–halo spin relation observed in various hydrodynamical simulations (Yang et al. 2023).

ACKNOWLEDGEMENTS

We would like to thank Simon D. M. White for useful discussions and advice. We also thank the anonymous referee for their helpful comments. We acknowledge the supports from the National Natural Science Foundation of China (grant no. 11988101) and the K. C. Wong Education Foundation. SS acknowledges support from NSFC grant (no. 12273053), CAS Project for Young Scientists in Basic Research Grant (no. YSBR-062), and the science research grants from the China Manned Space Project with no. CMS-CSST-2021-B03. AF is supported by a UKRI Future Leaders Fellowship (grant no. MR/T042362/1). RJW acknowledges the supports by the Forrest Research Foundation and the European Research Council via ERC Consolidator Grant KETJU (no. 818930).

This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, Durham University, and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

We thank the developers of the open-source python packages that were used in the data analysis of this work, including matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and h5py (Collette et al. 2020).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES