© The Authors 2025

1. Introduction

Observations show that Milky Way-like galaxies have formed stars at an almost constant star formation rate (SFR) over the past few gigayears (Binney et al. 2000; Wyse 2009; Madau & Dickinson 2014; Haywood 2014; Gondoin 2023). Given the fact that the gas depletion time in these galaxies is ≈1 − 2 Gyr (Saintonge et al. 2013; Tacconi et al. 2018), this would not be possible without a substantial accretion of gas onto the star-forming disc over time. Also, the observed metallicities of stars in the solar neighborhood require a consistent low metallicity source of accretion onto the galactic disc, known as the G-dwarf problem (e.g. van den Bergh 1962; Alibés et al. 2001). Thus, the accretion of gas from the surrounding environment is crucial for understanding galaxy formation and evolution.

Numerical simulations are a fundamental instrument in advancing our understanding of galaxy formation and evolution (see Vogelsberger et al. 2020, and references therein), offering key insights into the mechanisms driving gas inflows and star formation. Simulation results suggest that gas accretion occurs primarily in two main modes: cold and hot accretion. In the cold accretion mode (Kereš et al. 2005; Dekel et al. 2009; van de Voort et al. 2011; Nelson et al. 2016), which is typically dominant at high redshift (z ≳ 2), cold filaments from the intergalactic medium can survive their journey through the galaxy halo and deliver cold gas directly into the star-forming disc. On the other hand, in the hot accretion mode (Dekel & Birnboim 2006; Kereš et al. 2009; Nelson et al. 2013; Stern et al. 2020), gas that is falling into the dark matter potential well is shock-heated to the virial temperature, forming a halo of hot gas that subsequently must cool before being accreted onto the galaxy. This mode is generally prevailing at lower redshifts (z ≲ 2) for galaxies like the Milky Way and might be fundamental to the build-up of discs (Sales et al. 2012; Übler et al. 2014; Hafen et al. 2022; Yu et al. 2023).

This picture becomes more complex over several gigayears of evolution, as the environment around the galactic disc develops a complex, multi-phase structure, generating the so-called circumgalactic medium (CGM). The CGM is confirmed to be a multiphase medium by observations (e.g. Tumlinson et al. 2011; Putman et al. 2012; Anderson et al. 2016; Prochaska et al. 2017), containing cold (T ≲ 10⁴ K), warm (10⁴ ≲ T ≲ 10⁶ K) and hot (T ≳ 10⁶ K) gas. Furthermore, this component of the galactic ecosystem contains a significant fraction of the baryons associated with Milky Way-like galaxies (up to 50%; see Werk et al. 2014; Bregman et al. 2018). This makes the CGM the primary reservoir from which star-forming galaxies can get the fresh gas they need to sustain star formation. Hot gas in the CGM has been observed mainly through X-ray emission both in external galaxies (Anderson & Bregman 2011; Walker et al. 2015; Zhang et al. 2025; Locatelli et al. 2024) and in the Milky Way (Yoshino et al. 2009; Bluem et al. 2022; Bhattacharyya et al. 2023) as well as in the Large Magellanic Cloud (Gulick et al. 2021).

Despite these observations, one problem that is still not fully understood is how the hot gas in the CGM is cooling to accrete into the galactic disc. Although gas can cool through the spontaneous growth of thermal instabilities, this scenario has been ruled out by linear perturbation analysis (see e.g. Binney et al. 2009; Nipoti 2010); however, it could be still relevant in case of a substantial rotation of the gas in the inner regions of the corona (Sormani & Sobacchi 2019). Another interesting proposal is that gas may cool thanks to the mixing with the galactic fountains, gas clouds of mainly cold gas that are ejected from star-forming regions (i.e. the disc) by stellar feedback processes. These gas flows travelling through the CGM could enable a fraction of the hot gas condensate into the cold fountain clouds, increasing the amount of gas that then falls back onto the disc; therefore, it feeds the galaxy with low metallicity gas that can be used to form new stars. This scenario has gained support from hydrodynamical simulations both at parsec scales (Marinacci et al. 2010, 2011; Armillotta et al. 2016) and galactic scales (Hobbs & Feldmann 2020; Barbani et al. 2023).

Stellar feedback processes, such as supernova explosions, play a crucial role in shaping galaxy evolution and are responsible for launching galactic outflows that expel metal-rich gas from the galactic disc. These include galactic winds (Veilleux et al. 2005), which are strong outflows that are able to carry gas out of the galaxy entirely, and galactic fountains (Shapiro & Field 1976), which recirculate gas closer to the disc. The latter, in particular, are the primary focus of this work. The processes that involve the ejection of gas from the galactic disc and the subsequent accretion are generating the so-called baryon cycle. In particular, discovering the spatial distribution of the gas accretion and of the galactic fountains is crucial for the evolution of Milky Way-like galaxies. Studying gas accretion and outflows in external galaxies is notoriously difficult, but in recent years, the distribution of gas accretion has been analysed in observations (Zheng et al. 2017; Bish et al. 2019; Li et al. 2023; Zhang et al. 2025; Chastenet et al. 2024) and in numerical simulations (Übler et al. 2014; Nuza et al. 2019; Iza et al. 2021, 2024; Trapp et al. 2022; Roy et al. 2024).

The accretion of gas is expected to happen predominantly in the outer regions of the disc due to angular momentum conservation and to a larger cross-section of the external disc regions; therefore, a mechanism that move the accreted gas to the centre of the galaxy is needed to keep forming new stars. Lacey & Fall (1985) proposed three mechanisms that could generate these radial flows: if the specific angular momentum of the infalling gas is lower than that of the gas in the disc, we would expect a radial flow towards the centre. Other mechanisms proposed by Lacey & Fall (1985) include: the presence of viscosity in the gas layer and gravitational interaction between the gas in the disc and spiral arms and bars. The H I 21 cm emission line is the best way to investigate these radial motions in the outskirts of spiral galaxies, but it is still unclear if the radial inflows are sufficient to bring enough gas to the centre to sustain the star formation (Wong et al. 2004; Schmidt et al. 2016; Di Teodoro & Peek 2021). Finally, these radial flows are necessary also to reproduce the metallicity gradients (e.g. Mott et al. 2013) and the exponential density profiles (e.g. Wang et al. 2014) observed in spiral galaxies.

Understanding how CGM is accreted into the galactic disc is still puzzling. The goal of this work is to understand how the CGM can become available for the formation of new stars in the centre of the disc, therefore we investigate the distribution of vertical and radial outflow and inflows over the disc and how they are connected to the structure of the disc itself. To achieve this goal, we performed high-resolution hydrodynamical N-body simulations with the Stars and MUltiphase Gas in GaLaxiEs (SMUGGLE) model, an explicit ISM and stellar feedback model for the moving-mesh code AREPO. We simulated a set of nine star-forming Milky Way-like galaxies, including a hot gaseous halo surrounding the galactic discs. We varied the gas disc mass and size to quantify the key properties of vertical and radial gas inflows and outflows and investigate their connection to the evolution of the whole galactic disc. The paper is organised as follows. In Sect. 2, we describe the numerical methods and how we built the initial conditions for the simulations. In Sect. 3, we present our main findings, comparing them with previous works in Sect. 4. Finally, Sect. 5 provides a brief summary of our results.

2. Numerical methods

The simulations in this work were performed with the moving mesh code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020), a hydrodynamical N-body code that solves the equations of hydrodynamics on an unstructured Voronoi mesh that is able to move with the flow of the gas. On this moving mesh, AREPO integrates Euler equations with a finite-volume Godunov method that employs an exact Riemann solver. Gravitational accelerations were computed with an oct-tree method (Barnes & Hut 1986) and gravitational dynamics is evolved with a leapfrog algorithm. The code also uses a mesh refinement scheme, ensuring that each gas cell keeps its mass within a factor of two from a predetermined target value; together with the moving nature of the mesh this guarantees an automatic adjustment of the resolution in regions of high density. More details on AREPO can be found in Springel (2010) and Weinberger et al. (2020).

2.1. The SMUGGLE model

Coupled with the AREPO code, we used the SMUGGLE model, a sub-grid model, which incorporates baryonic processes that are fundamental for the evolution of star-forming galaxies. The major processes implemented by the model are the cooling and heating of the gas, star formation and evolution, and stellar feedback. We briefly describe these key constitutive elements of SMUGGLE below and we refer to Marinacci et al. (2019) for a complete description of the model.

SMUGGLE includes primordial and metal cooling (Vogelsberger et al. 2013). The cooling rates are tabulated as a function of density, temperature, and redshift using CLOUDY (Ferland et al. 1998), while also considering the presence of a UV background (Faucher-Giguère et al. 2009). An important feature of the model, allowing the gas to reach high densities and form new stars, is the generation of a low-temperature ISM phase (T ≈ 10 K) due the inclusion of low-temperature metal lines, fine structure, and molecular cooling (see Hopkins et al. 2018a). This low-temperature cooling is facilitated by the implementation of gas self-shielding from the UV background, which follows the parameterisation described in (Rahmati et al. 2013). In addition to stellar feedback (see below), gas heating can occur through cosmic rays (Guo & Oh 2008) and photo-electric heating (Wolfire et al. 2003), contributing to the thermal balance of T ≈ 50 K and T ≈ 8000 K gas in the ISM.

New stars can be formed in high-density (ρ > ρ_th = 100 cm⁻³) and gravitationally bound gas cells. If a gas cell is eligible for star formation, it forms stars at a rate Ṁ_⋆ = ϵ_SF M_g/t_ff, where M_g is the gas cell mass, t_ff is the free-fall time, and ϵ_SF(=0.01) is the star formation efficiency per free-fall time. Instead of forming stars in a continuous way, star formation works stochastically, transforming gas cells in star particles with a probability that is compatible with Ṁ_⋆. These star particles represent a single stellar population that follows a Chabrier (2003) initial mass function (IMF). As they evolve, star particles inject mass and metals – interpolated from a set of stellar evolution tables (see Vogelsberger et al. 2013, for details) – into the ISM, changing its chemical composition. Additionally, they deliver energy and momentum to the interstellar gas through stellar feedback processes, which are implemented by considering the following channels:

• (i)

  Supernova (SN) feedback is crucial in regulating star formation and the ISM structure through injection of momentum and energy (Martizzi et al. 2015; Kim & Ostriker 2018), while also being capable of ejecting multi-phase gas into the CGM (Fielding & Bryan 2022) and thereby distributing gas and metals into the galactic halo. To implement type II SN feedback, at each time step, the number of SNe is computed from the adopted Chabrier (2003) IMF (considering stars in the mass range 8 < M/M_⊙ < 100). Instead, the temporal distribution of type Ia SN events is parameterised with a delay time distribution (Maoz et al. 2012), from which the number of events is derived. From these rates, the injection of energy, momentum, mass, and metals into the ISM can be computed by assuming an energy of 10⁵¹ erg per SN. One problem in modeling SN feedback that arises due to insufficient resolution of galaxy scale simulations is how to properly describe the Sedov-Taylor phase of the SN remnant, in which most of the momentum imparted to the ISM is generated. To accurately capture the momentum generation, the model explicitly accounts for the momentum boost occurring during the Sedov-Taylor phase with a formulation derived from individual, high-resolution SN simulations (Cioffi et al. 1988, see also Hopkins et al. 2018b).

• (ii)

  Radiative feedback is produced by massive young stars (stellar particles with an age ≲5 Myr) and it is important for the thermal and dynamical state of the ISM. Indeed, several works have suggested that it can disrupt giant molecular clouds, thereby making SNe more effective (e.g. Walch et al. 2012; Smith et al. 2021). In SMUGGLE, this feedback channel includes photoionisation and radiation pressure. Photionisation is treated probabilistically and each gas cell surrounding a young star particle has a probability of being photoionised that depends on the ionising photon rate and the cell recombination rate. If a gas cell is photoionised, its temperature is set to 1.7 × 10⁴ K and its cooling is stopped for a duration equal to the time step of the (photo-)ionising star particle. The pressure generated from radiation produced by young stars is a source of momentum that is directly injected in the cells surrounding the star particle.

• (iii)

  Stellar winds are produced by OB and AGB stars, which contribute to feedback at different times due to the different evolutionary scales of these types of stars (e.g. Krumholz & Matzner 2009). To account for the effect of stellar winds, SMUGGLE computes the mass loss for each star particle at any given time step. This quantity is subsequently used to derive the energy and the momentum injected into the surrounding gas particles.

For all feedback channels, mass, linear momentum, and energy are coupled to the gas in a fully conservative manner, ensuring their conservation.

By including the aforementioned key physical processes, SMUGGLE is able to capture the complex multi-phase structure of the ISM, while simultaneously regulating star formation to match the observed levels in star-forming galaxies. Moreover, the explicit modeling of stellar feedback enables the self-consistent generation of galactic-scale gas flows to and from the disc. Owing to these capabilities, SMUGGLE has been successfully and extensively used to investigate various topics related to the evolution of galaxies and their ISM/CGM. These include, among others, the properties of giant molecular clouds (Li et al. 2020), the formation and evolution of young stellar clusters (Li et al. 2022), the development of cored dark matter profiles in dwarf galaxies (Burger et al. 2022a,b; Jahn et al. 2023), kinematics of stellar bars in spiral galaxies (Beane et al. 2023), impact of feedback on the distribution of O VI absorbers (Zhang et al. 2024a), formation of polycyclic aromatic hydrocarbons (Narayanan et al. 2023), evolution and distribution of super-bubbles in simulated disc galaxies (Li et al. 2024), and the effects of AGN feedback in galaxies with a multi-phase ISM (Sivasankaran et al. 2025).

2.2. Initial conditions

The generation of initial conditions is based on the method presented in Springel et al. (2005), which enables the creation of a multi-component galaxy in approximate equilibrium (see also Barbani et al. 2023, for a description of this method tailored to the type of simulations examined in this work). We considered a galaxy including the following components: a dark matter halo and stellar bulge, a thick stellar disc, a thick gaseous disc, and a galactic hot corona (i.e. in the hot phase of the galaxy CGM).

The dark matter halo and the bulge were modelled with a Hernquist profile:

$$\begin{array}{c}\hfill {\rho }_{\ast }(r)=\frac{M_{\ast }}{2\pi }\frac{a_{\ast }}{r{(r+a_{\ast })}^3},\end{array}$$(1)

where = dm for the dark matter halo and = b for the bulge. Also, a_* is the scale radius and M_* is the total mass. The dark matter was modelled as a static gravitational potential to reduce computational time.

The radial distribution of both stars and gas follow an exponential profile, as in observations (Freeman 1970). The surface density of this profile is given by

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\ast }(R)=\frac{M_{\ast }}{2\pi h_{\ast }^2}\exp (-R/h_{\ast }),\end{array}$$(2)

where =⋆ for stars and = g for gas. h_* is the scale length and Σ_* is the surface density of the stellar/gaseous discs. Vertically, the stellar disc has a sech² profile:

$$\begin{array}{c}\hfill {\rho }_{\star }(R,z)=\frac{{\mathrm{\Sigma }}_{\star }(R)}{2z_{\star }}{\mathrm{sech}}^2\left(\frac{z}{z_{\star }}\right),\end{array}$$(3)

where z_⋆ is the scale height and Σ_⋆ is the surface density of the stellar disc. The gas vertical profile is determined by iteratively solving the equation for vertical hydrostatic equilibrium $\frac{1}{{\rho }_{\mathrm{g}}}\partial P/\partial z=-\partial \mathrm{\Phi }/\partial z$, where Φ is the total gravitational potential and P is the gas pressure. In this way, we can obtain the gas density ρ_g, which is then used to evaluate

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{g}}(R)={\displaystyle {\int }_{-\infty }^{+\infty }}{\rho }_{\mathrm{g}}(R,z)\mathrm{d}z,\end{array}$$(4)

and the process is repeated until the desired gas surface density profile Σ_g(R) is obtained.

The metallicity in the gaseous disc is modeled as an exponential profile based on observation of the metallicity gradient in the Milky Way (Lemasle et al. 2018), with a slope of −0.0447 dex kpc⁻¹ and intercept of 0.3522 dex; thus, the metallicity at the solar radius of R ∼ 8 kpc is 1 Z_⊙.

To investigate the interaction between the disc and the surrounding medium, we added the hot component of the CGM (i.e. the galactic corona) around the galaxy. Here, we briefly explain how this component is included in the initial conditions. The coronal density profile is described by a β model, taking the form:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{cor}}(r)={\rho }_0{[1+{\left(\frac{r}{r_{\mathrm{c}}}\right)}^2]}^{-\frac{3}{2}\beta },\end{array}$$(5)

where β = 2/3 (Jones & Forman 1984; Moster et al. 2011), ρ₀ is the central density, and r_c is the core radius. Figure 1 shows the number density profile for the galactic corona. Points with errorbars are observational data for the density of the Milky Way corona (Putman et al. 2021), which are obtained from dwarf galaxies observations by evaluating the coronal density necessary to strip the gas of the dwarf at the perigalacticon. Hence, they should be considered as lower limits for the density of the corona. The β model profile of the corona has formally infinite mass, but we truncated the resulting density distribution at 10 R₂₀₀. To obtain the central density ρ₀ and the core radius r_c we applied a Markov chain Monte Carlo algorithm using an affine-invariant ensemble sampler (EMCEE package, Foreman-Mackey et al. 2013). From the Markov chain Monte Carlo (MCMC) fit we found a core radius of r_c = 83 kpc and a central density of ρ₀ = 3.3 × 10³ M_⊙ kpc⁻³. The coronal gas particles are initialised with an azimuthal rotational velocity v_ϕ = αv_c, where v_c = (R∂Φ/∂R)^1/2 is the circular velocity due to the total gravitational potential Φ. Also, α is set to 0.4 to ensure a rotational velocity v_ϕ ≈ 90 km s⁻¹ close to the galactic disc (Marinacci et al. 2011). Then, the temperature profile is obtained solving an equation analogous to the hydrostatic equilibrium equation using the effective potential, Φ_eff,

Fig. 1. Number density of the coronal gas as a function of radius. The black dashed line is the median model obtained by fitting a beta model to Milky Way coronal gas observations (upward pointing triangles with errorbars) using a Markov chain Monte Carlo algorithm. The coloured shaded area represents the 3σ region. Observations, which should be considered as lower limits for the coronal density, are taken from Putman et al. (2021).

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{eff}}=\mathrm{\Phi }(R,z)-{\displaystyle {\int }_{\infty }^R}\frac{v_{\varphi }^2(R\prime )}{R\prime }\mathrm{d}R\prime .\end{array}$$(6)

The effective equilibrium of each system was verified by running adiabatic simulations, that is, without dissipative processes such as radiative cooling. The metallicity of the coronal gas is set to Z = 0.1 Z_⊙, in line with observational estimates (e.g. Bogdán et al. 2017). Finally, we note that the profiles adopted to initialise all components of the systems in our simulations are axisymmetric. Therefore, the initial conditions have no azimuthal variations.

Table 1 lists the structural parameters of our fiducial Milky Way-like galaxy model, the other eight galaxies investigated in this work are variations of this fiducial model. In Table 2, we summarise the different parameter variations that we consider in our runs. Specifically, we scaled the mass of the gaseous disc M_g, by a factor f_M, and the gaseous disc scale length, h_g, by a factor f_G, with respect to the fiducial values 9 × 10⁹ M_⊙ and 7.6 kpc; in practice, our simulation set explores variations of a factor of two with respect to these reference values. The gravitational softening length is set at a minimum of ϵ_gas = 10 pc for the gas and at ϵ_⋆ = 21.4 pc for the stars for all the nine runs. In the fiducial run there are: N_disc = 3.2 × 10⁶ particles in the stellar disc, N_bulge = 8 × 10⁵ in the bulge, and N_gas ≃ 1.1 × 10⁷ gaseous particles, with N_gas, d = 3.2 × 10⁵ gaseous particles in the disc and N_gas, cor ≃ 1.1 × 10⁷ in the corona. This results in a mass resolution of m_gas = 1.1 × 10⁴ M_⊙, m_disc = 1.5 × 10⁴ M_⊙, and m_bulge = 2 × 10⁴ M_⊙ for the stellar disc and bulge particles, respectively. The number of gas particles in the disc in the simulations where the mass of this component is varied with respect to the reference value is changed accordingly to ensure the same mass resolution. All nine galaxies achieved an approximate state of equilibrium within ∼0.3 Gyr and were evolved for a total time span of 2 Gyr.

3. Results

We go on to discuss the main results of our simulations. In Sect. 3.1 we present the structure of the galaxies and of the gas in our simulations and we analyse the connection to the SFR. In Sect. 3.2, we focus on the vertical ejection/accretion of gas from and into the galactic disc, studying also the net accretion rate and how it changes between the nine simulations. In Sect. 3.3, we analyse the radial gas motions in the disc plane and their importance for the baryon cycle. Finally, in Sect. 3.4 we look at the temporal evolution of the disc size and its link to the vertical outflow and inflow rates.

3.1. Galaxy and outflow and inflow structure and distribution

We begin our discussion by presenting the general properties and the appearance of the gaseous discs in our runs. Figure 2 shows the gas projected density distribution (edge-on and face-on), after 1 Gyr of evolution, of all nine simulations considered in this work. Each row differs in the mass of the gaseous disc, whereas in columns only the initial scale length of the gaseous disc is varied. The simulations are arranged, from left to right, as follows (see Table 2 for their main properties): gal_m05r05, gal_m05 and gal_m05r2 (first row), gal_r05, gal and gal_r2 (second row), gal_m2r05, gal_m2 and gal_m2r2 (third row). In the second row, the central panel shows the fiducial galaxy; in the left panel the disc scale length is reduced by a factor f_G = 1/2 relative to the reference case, whereas in the right panel the scale length is scaled by f_G = 2 compared to the fiducial simulation. The first row displays simulations with a gaseous disc mass halved with respect to the second row (f_M = 1/2), while the third row shows simulations with gaseous disc mass doubled relative to the reference value (f_M = 2). Brighter colours represent higher surface density as indicated in the colorbar, which is kept the same in all panels to facilitate the comparison. The panels have the same physical size, 70 × 70 kpc for the face-on projections and 70 × 15 kpc for the edge-on projections. For subsequent analysis (see Sects. 3.2 and 3.3) it is useful to introduce a characteristic radius R_d representing the size of the gaseous disc. We define R_d as the radius at which the gas surface density drops below 1 M_⊙ pc⁻² (see also Di Teodoro & Peek 2021). It is important to note that this definition generally results in R_d differing from the disc scale length h_g set in the initial conditions. This characteristic radius is shown as a white dashed circle in Fig. 2 for the gal_r05, gal and gal_r2 simulations (second row). The adopted R_d definition captures well the disc extent in all the simulations analysed here.

Fig. 2. Gas column density face-on (top panels of each row) and edge-on (bottom panels of each row) projections computed using the SMUGGLE model at t = 1 Gyr for all the nine simulations examined in this work, as labelled in each panel. The simulations differ in the initial value of the scale-length and of the mass of the gaseous disc (see Table 2). Each panel is 70 kpc across and in projection depth with a total number of 1024 × 1024 pixels that give a resolution of ∼70 pc. Brighter colours correspond to higher column densities, as indicated in the colorbar. In the gal_r05, gal and gal_r2 projection panels, the dashed white circle shows the disc size, Rd, defined as the radius at which the gas surface density drops below 1 M⊙ pc−2 (not shown in the other projections for clarity).

All nine simulations produce stable discs, with no evidence of bar formation. The resulting disc morphology is similar in all runs despite the different scale lengths, which however influence the overall disc extent. The gaseous disc scale heights range from 0.3 to 0.8 kpc across all simulations and remains approximately constant over 2 Gyr. This scale height increases with radius due to the decreasing gravitational acceleration in the outer regions of the disc, leading to disc flaring. Looking at the face-on projections (top panels) we can appreciate the formation of a multi-phase ISM distributed in low density cavities and high density filaments, structures observed in recent JWST observations (Watkins et al. 2023; Thilker et al. 2023). Low-density cavities contain hot gas (T ≈ 10⁶ K) and are generated by SNe explosions in the galactic discs. These explosions heat and push the gas towards the edge of the cavities increasing, there, the cold gas surface density. Dense cold clumps are the environment in which new star formation is triggered. The newly formed stars can again inject momentum and energy in the ISM through feedback processes, disrupting cold gas clouds, eventually forming new cavities. These self-regulating mechanisms, implemented in SMUGGLE, can maintain the SFR stable over several Gyr.

Clustered SN explosions can generate large cavities (that can reach sizes up to 3 kpc, Egorov et al. 2017) in the ISM, also known as superbubbles (e.g. Reyes-Iturbide et al. 2014). These superbubbles may reach the vertical edge of the cold gaseous disc, eventually breaking through and pushing gas from the star-forming disc into the CGM. Stellar feedback processes in SMUGGLE are able to eject gas from the galactic disc, as clearly shown in the edge-on projections of Fig. 2 (bottom panels). The majority of these outflows appears to be confined in the first 1 − 2 kpc above the disc, but some of the ejected gas can reach distances up to 10 − 15 kpc. This ejected gas is a manifestation of the so-called galactic fountain (Shapiro & Field 1976; Fraternali & Binney 2008), as typically it has a velocity that is insufficient to escape the galactic potential. As a result, such gas is bound to the galaxy and will eventually fall back onto the disc. Fountain flows are present in all nine simulations as their generation is a natural consequence of the formation of new stars and the associated feedback. The distribution and the intensity of these galactic fountains is analysed in the following sections.

To demonstrate that SMUGGLE is able to regulate the SFR in our simulated galaxies, Fig. 3 shows the evolution of this quantity as a function of time for the entire simulation set. The early evolution of the SFR should be treated with caution. While we assume the corona is in equilibrium, cooling processes disrupt this balance, causing pressure loss and rapid accretion onto the disc. After about 0.3 Gyr, the SFR becomes stabilised due to the self-regulation between star formation and feedback in the SMUGGLE model. The corona also becomes more stable, though it continues to accrete onto the disc through radiative cooling. All the runs display, after this transient period in the first ∼300 Myr (not shown in the figure), an approximately constant SFR spanning the range 2 − 10 M_⊙ yr⁻¹. The average level of SFR depends on the initial properties of the disc and SFR variations for each simulation are confined within a factor of ∼2 from this average value. The gal simulation starts with a higher SFR that quickly stabilises at around 4 − 5 M_⊙ yr⁻¹. More specifically, the gal_r05 simulation starts with a higher SFR because of its higher gas surface density (it has the same gas mass contained in a smaller disc) but ends up with a lower SFR compared to the fiducial gal simulation. gal_r2 has the opposite behaviour, the SFR grows up over time surpassing the gal_r05 simulation at t ≈ 1 Gyr. Simulations with a larger disc mass have in general a higher SFR. Looking at the simulations with smaller disc mass: gal_m05r05 simulation has a SFR lower than the f_M = 1 counterpart gal_r05, this does not happen for gal_m05 and gal_m05r2. In this case, the SFR is lower in the first 0.8 Gyr and becomes slightly larger for the rest of the simulation.

Fig. 3. Star formation rate as a function of time for the nine simulations analysed in this work. The SFR has a wide range of values (≈2 − 10 M⊙ yr−1) across the different simulations, but it remains approximately constant over time for each simulated system. The average SFR depends on the disc initial properties, with variations in each simulation staying within a factor of ∼2 from the average value.

To better highlight the presence of a galactic fountain cycle in our simulations, Fig. 4 shows edge-on density-weighted projections of the metallicity (left panel), radial velocity (middle panel) and temperature (right panel) of the gas for the fiducial run gal after 0.5 Gyr of evolution. Each panel has a physical size of 70 × 70 kpc. The colour mapping encodes metallicities between 0.2 Z_⊙ and 2 Z_⊙, radial velocities between −150 km s⁻¹ and 150 km s⁻¹ and temperatures between 10² K and 10⁷ K (see the colourbar in each panel). As expected, the central region of the disc has a high metallicity (1 − 2 Z_⊙), which is consistent with the fact that the disc region is where the majority of the star formation is concentrated and stellar evolution can spread the bulk of metals in the ISM. The gas ejected by the galactic fountain has a metallicity consistent with the star-forming disc (≈1 Z_⊙), as it is expelled from that region by the action of stellar feedback. The CGM gas has an initial metallicity of Z = 0.1 Z_⊙, but we clearly notice the presence of metal-enriched gas (≈0.5 Z_⊙) at ≈30 kpc from the disc. This is the result of the mixing of disc gas, which is brought to these distances by the galactic fountain cycle, with the pristine coronal material. The middle panel shows the spherical radial velocity, red-coloured gas is moving radially away from the centre of the galaxy while blue-coloured gas is moving towards it. The galactic fountain cycle is clearly visible: gas that was pushed outside the galaxy (positive radial velocity) is then later re-accreted onto the disc (negative radial velocity). Looking at the metallicity projection we can notice that the gas moving towards the galaxy generally features a lower metallicity than the outflowing gas, which is consistent with the mixing scenario described above. Another important difference between the gas from the disc and the CGM is its temperature (right panel). The ISM is for the majority cold (T ≲ 10⁴ K) and this also the typical temperature of outflowing gas. As for the metallicity, when the fountain gas travels in the CGM it mixes and accretes material from it, forming an intermediate gas phase that subsequently cools and accretes onto the star-forming disc. As suggested in previous studies (see, e.g. Marinacci et al. 2010; Hobbs & Feldmann 2020; Barbani et al. 2023), this fountain-driven gas accretion may represent an important source of fresh gas for star-forming galaxies, directly impacting their evolution.

Fig. 4. Density-weighted maps of the gas metallicity, radial velocity and temperature (from left to right) in edge-on projections of the fiducial simulation at t = 0.54 Gyr. Each panel is 70 kpc across and in projection depth with a total number of 1024 × 1024 pixels that give a resolution of ∼70 pc. These projections show the stark difference between the disc and the CGM properties.

The presence of the galactic fountain cycle in each simulation is highlighted in Fig. 5 which shows density-weighted gas spherical radial velocity projections for the nine simulations (as indicated in each panel). The colour mapping is the same as central panel in Fig. 4. Each successive row (from top to bottom) increases the mass of the gaseous disc by a factor of 2, while each column from left to right represents a doubling of the disc scale length. The figure emphasises that simulations with a higher SFR (f_M = 2 simulations) have stronger ejections of gas from the disc, this results in clouds of cold gas that are pushed at larger distances from the disc and falling back into it after a larger time.

Fig. 5. Density-weighted maps of the gas radial velocity in edge-on projections of the nine simulations at t = 0.5 Gyr. Each panel is 70 kpc × 40 kpc and with a projection depth of 70 kpc with a resolution of ∼70 pc. The velocities go from −150 km s−1 (inflows, blue) to 150 km s−1 (outflows, red) as in central panel of Fig. 4. These projections clearly show the presence of galactic fountains in the simulations.

Figure 6 presents a density-weighted projection of metallicity for the gal simulation at t = 0.8 Gyr, considering only high-metallicity gas (Z > 0.9 Z_⊙), with the same colour mapping of the central panel of Fig. 4. The inflowing gas (blue-coloured) is enriched in metals, highlighting that some fraction of it originates from the disc. Therefore, a significant fraction of the inflowing gas stems from the galactic fountain process and is on its journey back onto the disc. This reinforces the idea that metal-enriched gas is recycled within the system, as there are no external sources of metallicity in the simulation.

Fig. 6. Density-weighted maps of the high-metallicity (Z > 0.9 Z⊙) gas radial velocity in edge-on projection of the gal simulation at t = 0.8 Gyr. The velocities range from −150 km s−1 (inflows, blue) to 150 km s−1 (outflows, red) as in central panel of Fig. 4. The projection clearly highlights how the accreted gas is not pristine material coming from the CGM, but is rather a mixture of the latter with gas that was previously ejected from the disc.

3.2. Vertical galactic fountain distribution

Galactic fountains play a fundamental role in galaxy evolution, by interacting with the CGM and bringing new gas to the disc (e.g. Marasco et al. 2022; Li et al. 2023). In this section we analyse how gaseous outflows and inflows are distributed over the galactic disc. The outflow and inflow rates are computed inside a slab centred at a height, z₀, and with a thickness, Δz, above and below the plane of the disc (see also Barbani et al. 2023). If a gas cell inside the slab meets all the required conditions (see below), it is counted as inflowing and outflowing gas and its contribution to the total mass rate is summed up as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{in}/\mathrm{out}}={\displaystyle \sum _i}\frac{m_i{v}_{z,i}}{\mathrm{\Delta }z},\end{array}$$(7)

where m_i is the mass of the gas cell and v_z, i is the z-component of its velocity. We set a velocity threshold of v_thres = 30 km s⁻¹; this value was chosen as an optimal compromise, ensuring that turbulent motions (tens of km s⁻¹) were not misclassified as inflows or outflows, while still accounting for most outflows with velocities in the range 50 km s⁻¹ < v_z < 100 km s⁻¹. If v_z, iz_i > 0 and |v_z, i|> v_thres, where z_i is the vertical distance from the disc plane of the gas cell, the gas in the cell was considered as a part of an outflow, whereas if v_z, iz_i < 0 and |v_z, i|> v_thres, it was counted as inflowing gas.

Figure 7 shows the outflow rate per unit area as a function of cylindrical radius normalised by R_d for the gal_r05 (left panel), gal (middle panel) and gal_r2 (right panel) simulations, computed at |z₀| = 0.5 kpc. This height corresponds to a ∼1 − 2 gaseous disc scale height and ensures that we focus on the region near the mid-plane, where inflows and outflows interact most directly with the star-forming disc. The disc size, R_d, is defined as the radius at which the gas surface density drops below 1 M_⊙ pc⁻² (see Fig. 2), resulting in R_d values of ≈20 kpc for the fiducial galaxy, ≈15 kpc for gal_r05 and ≈27 kpc for gal. These disc sizes are used to normalise the x-axis, which extends in the range 0 ≤ R/R_d ≤ 2. The gas has been divided in cold (T < 10⁴ K, dark blue shaded area), warm (10⁴ < T < 5 × 10⁵ K, blue shaded area) and hot (T > 5 × 10⁵ K, light blue shaded area). The outflow rate profiles have been computed at each time in the simulation with a temporal resolution Δt = 10 Myr (temporal difference between each snapshot of the simulation), the radial bins width is ΔR = 0.02 in R/R_d units. For each panel we computed the 25^th and 75^th percentiles of the resulting radial profiles and displayed their extent in the plot as the shaded coloured area, whereas the dashed black line indicates the median value. We notice that the outflows are dominated by the cold phase, which is on average at least one order of magnitude larger than the warm phase and two orders of magnitude larger than the hot one. Thus, the cold phase dominates the mass budget in the gas swept by SNe. The outflow rates are decreasing with the radius, following the distribution of gas and stars. At R ≲ 0.4 R_d, where more star formation is concentrated, the outflows are stronger. As expected, the distribution of the outflows is connected to the extension of the disc. All the simulations show a similar shape distribution with a break point at the edge of the disc R_d, beyond this radius the outflow rates drop significantly. In the fiducial simulation, the outflow rate is higher than 10⁻¹ M_⊙ yr⁻¹ kpc⁻² until R_d ≈ 20 kpc at all times. In gal_r05, that has a halved disc scale length, the outflows arrive to ≈15 kpc and in gal_r2 they extend up to ≈30 kpc. At t = 1 Gyr these three runs have a very similar SFR (4 − 5 M_⊙ yr⁻¹; see Fig. 3) and the intensity of the outflow rate in the first kpc is similar, but in gal_r2 the outflow are more extended. Therefore, the global outflow rate is higher in this system, with direct consequences on its evolution. For instance, more outflows can interact with the hot gas in the CGM, consequently increasing also the accretion rate. During the 2 Gyr evolution there is not a large change in the outflow rate profile and therefore the shaded areas are very close to the median value.

Fig. 7. Outflow rate per unit area as a function of the cylindrical radius, R (in units of the disc size, Rd) for gal_r05 (left panel), gal (central panel), and gal_r2 (right panel) simulations. Each shaded area indicates the outflow rate profiles between the 25th and the 75th percentile, whereas the black dashed lines represent the median values. The rates are divided in cold (T < 104 K, dark blue shaded area), warm (104 < T < 5 × 105 K, blue shaded area) and hot (T > 5 × 105 K, light blue shaded area) gas and computed at a height |z0| = 0.5 kpc inside a slab over the plane of the disc with a thickness Δz = 0.1 kpc and with |vz|> 30 km s−1.

Figure 8 shows inflow rate per unit area as a function of cylindrical radius normalised by R_d for the same simulations and in the same temperature range of Fig. 7. As this quantity is being computed at a relatively close distance to the disc mid-plane, namely, |z₀| = 0.5 kpc, the accreted gas has time to cool down before entering the gaseous disc, therefore the inflowing gas is mostly cold. Similarly to the outflow rates, the extension of the inflow rate distribution is connected to the size of the disc, more extended discs have more extended inflows, showing a break in the distribution at ≈R_d, after which the inflow rates decline rapidly. This happens because a large fraction of the inflows is gas that was beforehand ejected from the disc and subsequently falls back into it. The analogy between outflow and inflow rate profiles is again suggesting an essential role of the galactic fountains in setting the spatial distribution of the gas onto the disc. In the gal simulation the inflow rate arrives to a lower radius with respect to gal_r2, indicating that the majority of the accreted gas comes from the recycling of the ejected gas in a galactic fountain cycle. The other 6 simulations (not shown) behave in a similar way.

Fig. 8. Same as Fig. 7, but for inflow rates.

Our goal is to examine the net vertical mass inflow of gas onto the disc, considering that a significant portion of this inflowing gas was previously ejected from the disc as part of the galactic fountain process. Figure 9 shows the average net inflow rate, computed using Eq. (7) as Ṁ_net = Ṁ_in−Ṁ_out, divided in cold (red line), warm (orange line) and hot (yellow line) gas for the nine simulations analysed in this work (as indicated in each panel). The rates are computed at a height of |z₀| = 0.5 kpc from the disc plane with a Δz = 0.1 kpc width. Together with the net inflow rate we show the SFR radial profile (dashed grey line). As in the previous figure, the disc is divided into radial bins of width 0.02 in R/R_d units.

Fig. 9. Net vertical inflow rate, computed as Ṁnet = Ṁin−Ṁout, as a function of cylindrical radius R (in units of Rd) for the nine simulations divided in cold (T < 104 K, red line), warm (104 < T < 5 × 105 K, orange line) and hot (T > 5 × 105 K, yellow line) gas. The dashed grey line is the SFR profile. The solid lines represent the profile averaged over the entire simulation and the shaded areas are the standard deviations. The rates are evaluated inside a slab over the plane of the disc with a width Δz = 0.1 kpc at a height |z0| = 0.5 kpc. Each point represents the gas net inflow rate within a circular annulus with a radial width of 0.02 in R/Rd units and considering only gas cells with |vz|> vth = 30 km s−1. In each panel we have indicated the average net vertical rate at R < 0.4 Rd, this is roughly the region of the disc where 90% of the SFR is concentrated in all nine simulations. The net inflows can extend beyond Rd supplying the outer disc with fresh gas for star formation.

Figure 9 emphasises the connection between accretion from the CGM and galactic disc structure. In all the simulations, the net inflow rates exhibit a bell-shaped distribution. These accretion rates extend up to ≈1.2 R_d in all the simulations, which corresponds to ≈24 kpc in the f_G = 1 simulations, ≈18 kpc in the f_G = 0.5 simulations and ≈32.5 kpc in the f_G = 2 simulations, following the respective discs sizes, R_d. It is worth to notice that the peak of the total net inflow rate is not located at the centre of the galaxy, but is instead skewed towards larger radii. The peak distance from the centre is increasing with disc size and is typically around 0.75 R_d. The grey dashed line shows the SFR as a function of radius, which has a central peak and also a bell shape but is more concentrated in the inner region of the disc, where the gas surface density is higher. If the gas were to fall back to roughly the same radial position from which it was ejected, the inflow distribution would closely resemble the SFR distribution. Instead, the majority of the gas is accreted in the outer regions of the disc. This happens because the orbital paths of galactic fountain clouds increase with radius, due to the progressively weaker gravitational attraction moving towards the external regions of the galaxy (e.g. Li et al. 2023). Additionally, the region of the CGM affected by the fountain clouds passage – and where the warm gas phase is generated – expands with their orbital path, facilitating an easier accretion into the disc. The total net accretion rate Ṁ_tot (R < R_d) is indicated in each panel. This total accreted gas would be enough to sustain the SFR (Fig. 3), a fraction of this accretion in the external regions is used to grow the disc inside-out (see Sect. 3.4). While most of the gas is accreted in the outskirts of the disc, a fraction of it falls in the central region (R < 0.4 R_d), in each panel of Fig. 9, we indicated the average value of the net rate in this region. Overeall, 90% of the SFR is located at R < 0.4 R_d, thus gas accreted in these central regions can be directly available for the formation of new stars, whereas the gas accreted at R > 0.4 R_d cannot be immediately used for star formation as the average gas density is not high enough. This average net rate, while not negligible, is still not sufficient to entirely sustain the SFR of the galaxy (see Figs. 3 and 11). Therefore, there needs to be an alternative mechanism that feeds the star formation in the centre. The accretion budget and the comparison with the SFR are discussed in Sect. 3.3.

When examining the net inflow rate separated in hot, warm and cold components, we notice several interesting features. We have seen that the outflow rate of warm and hot gas is almost negligible (Fig. 7), indicating that the net warm and hot inflow primarily originates from the CGM rather than just being gas previously ejected from the disc. In fact, as shown in Barbani et al. (2023), the warm gas phase is generated by the passage of the cold fountain clouds through the hot CGM. In all simulations, the hot gas has an almost flat distribution that slightly decreases with increasing radius, with values ranging from 0 to ∼0.15 M_⊙ yr⁻¹. This likely represents gas directly accreted from the hot CGM via a cooling flow, showing no clear correlation with the cold and warm gas distributions. The warm and cold gas, however, are more closely connected, with their net inflow peaks occurring at similar radii. The net cold inflow rate shows more discontinuous peaks and it is not consistently positive¹. As shown in Fig. 7, nearly all of the ejected gas is in the cold phase, leading to regions of the disc where the outflow rate exceeds the inflow rate (resulting in a negative net rate) and vice versa. This does not happen for the net warm inflow rate, which remains consistently positive and generally higher than the net cold inflow rate. The spatial distribution of cold and warm inflows appears interconnected, suggesting that cold outflows may enhance the formation of inflowing warm gas. As discussed extensively in Barbani et al. (2023), mixing between the cold fountain clouds and the hot CGM can reduce the cooling time of the gas mixture, allowing the gas to cool more rapidly and accrete onto the disc, thereby supplying fresh gas to sustain star formation.

The different disc sizes and gaseous masses also result in variations in the net inflow distribution. The simulations with the smaller discs (gal_r05, gal_m05r05, gal_m2r05) show the lowest net inflow rates, this is a combination of both a lower SFR (see average SFR indicated in each panel of Fig. 11) and a smaller extent of the disc. The hot gas accretion rate increases with disc size, reflecting the slow accretion of the galactic corona; a larger disc is more capable of capturing and accreting more hot gas. The cold gas rate, that derives from stellar feedback, is also lower in smaller discs due to their lower SFR, which in turn reduces the generation of warm gas, formed from the mixing and condensation of CGM into the galactic fountain cloud. Consequently, the rates increase for larger discs, which also present on average a higher SFR. Therefore, if the mass of the gaseous disc decreases, gas accretion diminishes for the same reasons.

3.3. Radial gas motions

Figure 9 illustrates that the majority of the gas is accreted at a significant distance from the galactic centre (≈0.75 R_d), near the edge of the gaseous disc. Star forming galaxies need to acquire gas from the external environment to keep forming stars with an almost constant SFR (e.g. Cignoni et al. 2006; Isern 2019; Mor et al. 2019), as the gas within the disc is not sufficient. The majority of stars form in the central regions of the disc, where gas density is the highest. Therefore, the accreted gas must travel from the disc outskirts to the centre in order to be accessible for the formation of new stars (Schmidt et al. 2016; Di Teodoro & Peek 2021).

A key aspect of understanding radial gas motion within the disc is the distribution of angular momentum in both the gaseous disc and the CGM. Figure 10 presents the specific angular momentum profiles for the nine different galaxy models, computed as j = J/M, where J is the total angular momentum and M is the total mass of each respective component. The angular momentum is evaluated near the disc (|z|< 1 kpc), with CGM gas selected by metallicity Z < Z_⊙ and disc gas by Z ≥ Z_⊙. The figure reveals a consistent trend: in the inner regions the specific angular momentum of the gaseous disc is systematically higher than that of the CGM throughout the entire time span of the simulations. Within R_d = 1, the disc angular momentum increases linearly closely following the trend expected given the disc rotation curve, indicating full rotational support. Beyond this point, the disc component fades, and its angular momentum declines. A similar trend is observed in the inner CGM, though at lower values, suggesting some degree of co-rotation with the gaseous disc and a spin up of the corona in the regions close to the disc, likely driven by galactic fountains, as also indicated by cosmological simulations (e.g. Grand et al. 2019). At larger radii (R/R_d ≳ 1), the disc angular momentum diminishes, while the CGM angular momentum initially decreases before slowly increasing outwards, marking the transition to a regime where the coronal gas ceases to co-rotate with the disc. The impact of galaxy properties is evident in the figure. In the f_g = 2 simulations, increasing the total gas mass of the disc leads to a noticeable increase in the disc specific angular momentum at all radii. In the other cases, changes in gas mass do not significantly affect the absolute angular momentum values. Variations in disc scale length also influence the angular momentum distribution, with larger scale lengths generally corresponding to higher disc angular momentum at fixed gas mass. Despite these differences, the higher angular momentum of the disc relative to the CGM in the inner regions remains a robust feature across all models.

Fig. 10. Specific angular momentum, computed as j = J/M, where J is the total angular momentum and M is the total mass of each respective component, for the gaseous disc and the CGM computed as a function of cylindrical radius R (in units of Rd) for the nine simulations. The angular momentum is measured near the disc (|z|< 1 kpc), with CGM gas selected by Z < Z⊙ and disc gas by Z ≥ Z⊙. In each panel, the grey lines represent the specific angular momentum profile at a given time averaged in a time span of 100 Myr. The black line is the average specific angular momentum and the blue (red) shaded area represents the standard deviation of the disc (CGM) specific angular momentum distribution. Each panel represents a different combination of disc mass and scale length, showing how these parameters influence the angular momentum distribution. Across all models, the disc maintains a systematically higher specific angular momentum than the CGM in the inner regions following the rotation curve, while at R/Rd ≳ 1, the disparity between the two components becomes more variable, marking the transition where the disc fades.

To investigate the radial migration of gas we have computed the radial mass rate in each radial bin as ${\dot{M}}_{R,\mathrm{in}/\mathrm{out}}={\displaystyle \sum _i}{v}_{R,i}{m}_i/\mathrm{\Delta }R$, where v_R, i is the cylindrical radial velocity of the i-th gas particle and ΔR is the radial bin size. The particles are selected within a slab centred on the disc plane (z₀ = 0 kpc) with a thickness of Δz = 0.5 kpc to ensure that only disc particles are included. The radial outflow (inflow) rate is determined by considering particles with v_R > 0 (v_R < 0). The net mass rate is computed as Ṁ_R,net =Ṁ_R,in−Ṁ_R,out. The radial velocity is computed by averaging the cylindrical radial component of the velocity weighted by the mass of each gaseous cell. For both quantities, the disc is divided in annuli with a spatial extent ΔR = 0.02 in R/R_d units and the profiles are computed in the radial range 0 ≤ R/R_d ≤ 2.

Figure 11 shows the net radial mass rate of the gas for the nine simulations, as indicated in each panel. The grey lines represent the net inflow rate profile at individual times, while the black line shows the average profile over a 2 Gyr time span, with the shaded area indicating the standard deviation. The average value of the gas mass rate is positive throughout the entire disc, indicating a significant radial mass flow towards the inner regions. When the accreted material has a specific angular momentum that is different than the one in the disc, radial flows will inevitably arise owing to angular momentum conservation (Mayor & Vigroux 1981). Since the gas accreted from the CGM has a rotational velocity that is lower than that of the disc, it will naturally tend to migrate towards the inner regions of the disc, forming a radial flow that feeds the star formation. This is apparent in Fig. 10, which shows that the specific angular momentum of the CGM is systematically lower than that of the disc in the inner regions. This discrepancy drives angular momentum redistribution and facilitates radial transport of gas towards the centre.

Fig. 11. Gas net radial mass rate for the simulations analysed in this work. In each panel, the grey lines represent the mass rate profile at a given time averaged in a time span of 100 Myr. The black line is the average mass rate and the blue shaded area represents the standard deviation. On the bottom of each panel the net radial mass rate averaged within the radial range 0 ≤ R/Rd ≤ 1 and the average SFR, computed in the temporal range 0.3 < t < 2 Gyr, are indicated.

Generally, the gas net mass rate is near zero in the disc outskirts and it moves towards more positive values (indicating inflows) in the inner regions of the disc across all simulations, eventually approaching zero again at the centre. This mass inflow becomes less prominent in the inner regions for two main reasons: (i) the higher gas density at the centre and (ii) more importantly, the fact that a portion of the gas is consumed to form stars in these central regions. A consistent pattern emerges in all the simulations, a peak in the mass rate that is generally found roughly at the edge of the disc (R ≈ R_d), except for simulations gal_m2r05 and gal_m05r05. The cumulative star formation rate of the galaxy follows a similar pattern, exhibiting a decline for radii R/R_d ≲ 0.4, where the large majority of the SFR is located. This similarity in trends emphasises that star formation primarily takes place in the inner disc regions, with the radial mass flux diminishing towards the centre as it feeds star formation. In all simulations, the SFR remains higher than the radial net inflow rate. This discrepancy is compensated by vertical inflows from the CGM, which, together with radial mass transport, sustain the SFR in the disc (see below). When comparing simulations with the same gas disc scales, we observe values of Ṁ_R,net within a factor of 2, except for the f_G = 1 simulations where the differences in the average mass rate value are almost negligible. The primary difference arises when the size of the gaseous disc is changed: the gas mass rate increases with larger disc sizes. A larger disc extends into regions where the CGM has lower angular momentum, facilitating easier inward motions towards the centre of the disc.

Overall, while the simulations show a non-negligible net inflow rate, this rate is generally lower than the average SFR. For instance, the net radial inflow rate is 1.30 M_⊙ yr⁻¹ for the fiducial run, then 0.71 M_⊙ yr⁻¹ and 2.33 M_⊙ yr⁻¹ for gal_r05 and gal_r2, respectively, compared to the average SFR of 4.5 M_⊙ yr⁻¹, 4.29 M_⊙ yr⁻¹, and 4.5 M_⊙ yr⁻¹ for these cases. Thus, radial mass motion towards the centre alone may not fully sustain the SFR over time. However, it is important to note that not all gas is accreted at the edge of the disc; a significant fraction, approximately 25% of the total accreted gas, is directly deposited in the inner regions (R < 0.4 R_d; see Fig. 9), where 90% of star formation takes place. Combined with radially transported gas, this could be sufficient to sustain the global SFR of the galaxy. Therefore, for the SFR to be fueled by these two phenomena in conjunction we would need to have Ṁ_SFR =Ṁ_z,net(R < 0.4 R_d) + Ṁ_R,net ≈ SFR. In fact, this is true for most of the nine simulations analysed in this work (see Table 3). For instance, in gal_m05r05 we find Ṁ_SFR ≈1.92 M_⊙ yr⁻¹, which has to be compared with the average SFR ≈ 2.43 M_⊙ yr⁻¹. Even though strictly speaking the rate of replenishment of gas is below the SFR, the two values are comparable. Moreover, we have to consider that Ṁ_SFR and SFR are averaged over a 0.3−2 Gyr time span, meaning that our analysis is smoothing what occurs over the evolution of the galaxy. Therefore, these values should be interpreted with some caution, but nonetheless they provide insight into whether the formation of new stars can be sustained by gas accretion from the CGM. Other systems show a similar behaviour. Overall, as seen in Table 3, the majority of the simulations have Ṁ_SFR ≈ SFR, however there are some exceptions: in gal_r05 we find Ṁ_SFR ≈ 1.42 M_⊙ yr⁻¹ which accounts for less than half of the average SFR (≈4.29 M_⊙ yr⁻¹). Similarly, in gal_m2r05 and gal_m2, the SFR deviates from Ṁ_SFR, as it is partially sustained by the larger initial gas reservoir in the disc at the start of the simulations.

Figure 12 shows the average radial velocity, with values ranging between −4 and −6 km s⁻¹. All the simulations exhibit a similar radial velocity profile: the velocity remains almost constant within the gaseous disc (R < R_d), rapidly decreases outside the disc and then start to slightly increase again (except in gal_m2r05 and gal_m2). This profile shape appears to be independent of the disc structural characteristics. However, the magnitude of the velocity slightly changes between discs of different sizes. Specifically, the radial velocity increases with disc size, similarly to the mass rate. The average velocity inside the disc increases between f_M = 1 and f_M = 0.5 simulations, when the mass of the gaseous disc is smaller, probably because the lower average ISM density allows gas to flow more easily toward the inner regions. The velocity is higher in the outskirts, beyond the galactic disc, where gas is accreted radially on the edge of the galaxy, allowing the disc to grow in size over time (see Fig. 14). Moving inwards, the radial velocity decreases and approaches zero at the centre of the disc.

Fig. 12. Gas cylindrical radial velocity for the simulations analysed in this work. In each panel, the grey lines represent the (mass weighted) radial velocity profile at a given time averaged in a time span of 100 Myr. The black line is the average cylindrical radial velocity and the green shaded area represents the standard deviation. On the bottom left corner the radial velocity averaged within the radial range 0 ≤ R/Rd ≤ 1 is indicated.

Figure 13 presents the radial net rate of the gas as a function of time for the three simulations with the standard gaseous disc mass (gal_r05, gal and gal_r2), the other 6 simulations are not shown as they all have a very similar trend. For each simulation, we computed the total radial net rate by averaging it over the whole galactic disc (i.e. at R < R_d and |z|< 0.5 kpc) at each time (with a temporal resolution ≈10 Myr). The trend is consistent across all simulations, alternating between periods of negligible or negative net rate and periods of positive net rate. This fluctuation is likely connected with the stellar feedback cycle; when SNe explode, the net rate is dominated by the outflowing gas. As negative stellar feedback reduces star formation, thereby decreasing the injection of energy from new SNe, fresh gas can flow towards the centre once again, fueling star formation and repeating the cycle. Overall, the gas predominantly flows towards the centre (positive net rate), suggesting that the radial gas flow phenomenon is independent of the gaseous disc size, given that the rotational velocity of the coronal gas has the same values at the beginning of all simulations.

Fig. 13. Radial gas net rate as a function of time for the gal_r05 (dotted line), gal (solid line) and gal_r2 (dashed line) simulations. At each time the radial net rate averaged within the disc, at R < Rd and |z|< 0.5 kpc, is displayed. Overall the radial net rate is mostly positive, indicating a flow of gas towards the centre of the disc.

To summarise, the simulations reveal a general inflowing trend that, while insufficient to sustain the star formation in the centre alone, is nonetheless crucial for feeding the inner regions of the galaxy, working in conjunction with vertical gas accretion in the central regions (R < 0.4 R_d).

3.4. Inside-out evolution

The vertical accretion in the outer regions of the disc, coupled with the subsequent radial motions, drives the gas from the periphery towards the disc centre. Inevitably, this accretion increases the gas density in progressively more outer regions, consequently enhancing star formation activity in these external areas. Over time, this mechanism enables the disc to grow from the inside to the outside, creating what is known as inside-out evolution (e.g. Firmani & Avila-Reese 2009; Brooks et al. 2011; Pezzulli & Fraternali 2016; Li et al. 2023). To explore the relationship between gas accretion and disc growth, we analysed the evolution of different scale lengths over time and their connection to the total outflow and inflow rates. The top panels of Fig. 14 show the evolution of the scale lengths of stars (orange line), gas (light brown line) and SFR (brown line) for the nine simulations as a function of time. These scale lengths are determined by computing the radius which encloses 90% of the stars/gas mass and of the SFR, with each value smoothed using a moving average over a 100 Myr interval. From left to right each panel shows the f_G = 0.5, f_G = 1 and f_G = 2 simulations respectively, with solid lines indicating the f_M = 1 runs, dashed lines indicating the f_M = 0.5 runs and dotted lines indicating the f_M = 2 runs. All the simulations show a periodic trend with oscillations that can be linked to the typical time scale of a fountain cycle in Milky Way-like galaxies (t_fount ≈ 100 Myr, e.g. Fraternali 2014). Over a span of 2 Gyr all three scale lengths have increased.

Fig. 14. Top panels: Star (orange line), gas (light brown line), and star formation rate (brown line) scale length for the gal_r05 (solid), gal_m05r05 (dashed), gal_m2r05 (dotted, left panel), gal (solid), gal_m05 (dashed), gal_m2 (dotted, central panel), and gal_r2 (solid), gal_m05r2 (dashed), gal_m2r2 (dotted, right panel) simulations as a function of time. These scale lenghts are computed as the radius that encloses 90% of stars/gas mass or SFR. Bottom panel: total outflow (blue) and inflow (light blue) rates of the same simulations of the top panels as a function of time. The simulations with fM = 0.5 are indicated with a dashed line and the simulations with fM = 2 with a dotted line. In the central panels the peak in the inflow rate of the gal simulation are highlighted with dashed grey vertical lines. Each of the scale is connected and they have an up-and-down trend linked to the gas outflow-inflow cycle.

We start by focusing on the fiducial run gal (solid lines, central panel): the gas scale length has expanded from ≈12 kpc to ≈14 kpc, the star scale length from ≈8 kpc to ≈9 kpc and the SFR scale length from ≈6 kpc to ≈9 kpc. The gas scale length remains consistently larger than the stellar scale length by about a factor of ≈2. Notably, both the gas and stellar discs exhibit a concurrent growth, expanding steadily over a 2 Gyr period. The stellar scale shows higher peak frequencies, with approximately two peaks for every peak in the gas scale length. Interestingly, the peaks in the SFR do not coincide with those of the gas and star scales; rather, they tend to align with depressions in the gas scale length. The SFR exhibits the greatest amplitude variation, whereas the stellar scale remains more stable, with peaks that nearly disappear in the other simulations. An increase in the SFR spatial extent implies that more stars are being formed in the outer regions of the disc compared to previous times, leading to more outflows ejected from these areas. Consequently, these episodes of gas ejection may reduce the gas distribution scale length, with a fraction of these outflows eventually falling back onto the disc, thus increasing the gas scale and repeating the cycle again. The gal_m05 (central panel, dashed lines) and gal_m2 (central panel, dotted lines) simulations, with a gaseous disc mass halved and doubled, respectively, have similar behaviours. The values of the three scales are similar to gal, with the gal_m05 scales slightly lower and the gal_m2 scales slightly larger than the fiducial run scales. Also in these cases the peaks in the gas scale and in the SFR scale are misaligned. Similar features can be seen in the other simulations. The simulations with f_G = 0.5 (left panel) show lower values in all the three scales as expected and they are all increasing over the 2 Gyr of the simulations. In the f_G = 2 runs (right panel) the scales are higher than in the f_G = 1 runs, gal_r2 gas scale is passing from ≈18 kpc to ≈20 kpc, whereas in gal_m05r2 the gas scale stays almost constant at 16 kpc. Overall, the increasing trend in these scale lengths in all the simulations over time indicates the growth in the size of the star-forming disc, which is consistent with an inside-out evolution of the disc.

In the bottom panels, we display the total outflow (dark blue line) and inflow (light blue line) rates for the nine simulations, computed at z₀ = 0.5 kpc, Δz = 0.1 kpc and R < R_d. In general, the inflow rate is larger than the outflow rate. In the fiducial run gal (central panel, solid line) the net inflow rate is Ṁ_net =Ṁ_in−Ṁ_out ≈ 10 M_⊙ yr⁻¹. This net inflow is larger than the average SFR observed in the fiducial simulation (see Fig. 3), confirming that the formation of new stars is largely sustained by gas accretion from the CGM, the same applies to the other simulations. An important aspect to notice is that the outflow and inflow rates exhibit similar trends but with an offset in time, with depressions in the outflows corresponding to peaks in the inflows, and vice versa. This pattern arises due to the galactic fountain cycle: a certain amount of gas is ejected from the disc as a galactic fountain due to clustered type II SN explosions, creating a peak in the outflow rate. After a time interval of approximately t_fount, the fountain gas falls back onto the disc, leading to a peak in the inflow rate. Periods of higher inflow rates result in an increased and more extended SFR, thereby enlarging the SFR scale length. This occurs because gas accretion is not concentrated in the disc inner regions but is rather shifted with respect to the peak of the SFR (see Fig. 9). Vice versa, when the outflows are dominant, the SFR scale length diminishes.

In all the simulations we find that the peaks in the inflow rate roughly coincide with the peaks in gas and stellar scale lengths, we highlight these peaks for gal with dashed vertical lines in the central panel. In contrast, the peaks in the outflow rate correspond to the SFR scale peaks (see also Zhang et al. 2024b, for a similar effect on dwarf-mass galaxies). This shows us the intricate connection between outflows/inflows and SFR, which is mediated by stellar feedback and the consequent gas circulation. As discussed earlier, vertical gas inflows are distributed throughout the disc plane, with the majority of the gas accreting at ≈0.75 R_d (Fig. 9). This accretion pattern increases gas density in these external regions, leading to an expansion in the size of the gaseous disc and consequently also the stellar disc. The SFR peaks, along with the associated stellar feedback, generate outflows, which are fundamental to the cyclical process of disc growth and evolution.

4. Comparison with previous works

Despite the challenges, gas accretion has been studied both in the Milky Way and in external galaxies. Evidence of cold gas accreting in the Milky Way is observed through intermediate velocity clouds (IVCs, 40 < v_LSR < 90 km s⁻¹) and high velocity clouds (HVCs, v_LSR > 90 km s⁻¹) (e.g. Wakker & van Woerden 1997; Lehner & Howk 2010), gas clouds with anomalous kinematic with respect to the rotation of the disc. HVCs usually have higher velocities, lower metallicities (0.1 < Z/Z_⊙ < 1) and are located farther from the disc (up to 10 kpc) with respect to IVCs, that have higher metallicities (Z ≈ Z_⊙) and are usually close to the disc (within 2 kpc). Historically, these differences have been interpreted as a dichotomy in the origin of these two populations of clouds: the IVCs are produced by the cycle of gas inside galactic fountains and the HVCs are coming from the accretion of gas from the CGM (Peek et al. 2008). However, in recent years this picture has been challenged and it has been proposed that HVCs and IVCs are both manifestation of the galactic fountain phenomenon at different velocities (Marasco et al. 2022; Lehner et al. 2022). The IVCs are outflows of gas ejected at lower velocity, therefore they are located closer to the disc and can maintain a metallicity similar to the ISM. HVCs are pushed farther away from the disc and therefore they can mix and decrease their metallicity. Zhang et al. (2024a) employed the SMUGGLE model in simulations of isolated Milky Way-like galaxies, without the inclusion of a CGM, to study O VI absorbers. Their findings suggest that key observational properties of low-velocity (v_LSR ≲ 100 km s⁻¹) O VI absorbers, such as column density and scale height, are well reproduced by the inclusion of galactic fountains driven by stellar feedback. However, the absence of the CGM in their simulations may lead to an underestimation of the column density of high-velocity (v_LSR ≳ 100 km s⁻¹) O VI absorbers.

Extrapolating the accretion rate values of observed high- and intermediate-velocity clouds to the whole disc, observations find a maximum accretion rate for HVCs of 0.4 M_⊙ yr⁻¹ (Putman et al. 2012) and of 1.3 − 4.3 M_⊙ yr⁻¹ for IVCs (Röhser et al. 2016), which would be enough to sustain the SFR in the Milky Way (SFR ≈ 1.5 − 3 M_⊙ yr⁻¹, Chomiuk & Povich 2011; Licquia & Newman 2015; Elia et al. 2022). HVCs and IVCs are clearly visible in the gas projections of our simulations (see Figs. 4 and 5). Assuming that 40% of the total inflow rate, which averages between 10−20 M_⊙ yr⁻¹ (see Fig. 14), is attributed to HVCs and IVCs (Putman et al. 2012; Trapp et al. 2022), we find values consistent within a factor of 2 with observed data. It is important to note that the nine simulated star-forming galaxies have significantly different SFR (higher than the Milky Way), which will be reflected in the resulting inflow rates.

Li et al. (2023) conducted a study on the extraplanar gas of the spiral galaxy NGC 2403, employing a dynamical galactic fountain model. Their findings suggest that fountain-driven accretion can be responsible for an inside-out evolution, effectively cooling the CGM at large radii. In the absence of galactic fountains, cooling would occur predominantly at the centre of the CGM, where gas density is the highest. However, their model applied to observations demonstrates that most gas is accreted at larger radii, a result that aligns closely with the outcomes of our simulations. They observed outflow/inflow rates per unit area similar to those in our work. In their case, the rates peak at the centre. In contrast, our simulations show a central ‘hole’ due to the height at which our rates are computed. Li et al. (2023) derived the rates at z = 0 kpc, a condition that we cannot replicate in our simulations. They derived a coronal gas inflow rate whose shape can be compared to our net vertical inflow rate (Fig. 9). We both find a bell-shaped distribution, with the peak skewed towards larger radii. As already mentioned, this is likely due to the increased orbital path of the galactic fountains with radius and to the fact that the SFRD profile is declining together with the radius.

Different works have focused on the analysis of radial mass flows and their relation to galaxy evolution. Observationally, there is not a clear evidence of the existence of radial mass fluxes strong enough to sustain the formation of new stars in the central regions of disc galaxies. Some references find radial inflows that are comparable with the SFR of the host galaxy, whereas in other works such inflows are not present. For instance, Schmidt et al. (2016) studied 10 local spiral galaxies using the H I Nearby Galaxy Survey. They derived average radial gas mass rates of ≈1 − 3 M_⊙ yr⁻¹, which will be sufficient to explain the level of SFR in such galaxies. On the other hand, Di Teodoro & Peek (2021) analysed a sample of 54 local spiral galaxies and found small inflow rates with velocities of a few km s⁻¹ in only half the galaxies of the sample. They computed an average inflow rate over the entire sample of ( − 0.3 ± 0.9) M_⊙ yr⁻¹ at radii larger than the optical disc (inside which most of the star formation occurs). They concluded that the radial inflows are not sufficient to sustain the star formation across the disc and therefore that a different mechanism in addition to secular accretion is necessary. We find that net radial mass rate value, averaged across our nine simulations, is ≈1.72 M_⊙ yr⁻¹ with a standard deviation of ≈0.86 M_⊙ yr⁻¹, which lies somewhat in between the results of Schmidt et al. (2016) and Di Teodoro & Peek (2021). The net radial mass rates estimated from our numerical experiment can sustain a significant fraction of the SFR observed in simulations. This gas radial motion combined with gas accretion in the central regions of the disc can maintain an active level of star formation in the simulated systems.

Finally, in recent years, some authors have been using cosmological simulations of Milky Way analogues to investigate vertical and radial gas inflows. For example, Trapp et al. (2022) studied gas accretion using FIRE-2 cosmological zoom-in simulations of 4 Milky Way-like galaxies. They found that gas is accreted preferentially in the outskirts of the galaxy, from where the gas moves towards the centre with a mass rate comparable to the galaxy SFR. In their simulations, gas vertical inflows present a bell-shaped distribution and gas particles join the disc at 0.5 − 1 R_d, similarly to our results. Radial mass rate and velocities present trends broadly consistent with our findings, but in our case the radial gas motion is not sufficient to entirely sustain the star formation activity of the galaxy. Moreover, Trapp et al. (2022) show that most accretion occurs co-rotating and parallel to the disc plane, with only a subdominant contribution from HVCs and IVCs, whereas our results indicate that the majority of the gas is accreted vertically onto the disc, while still having significant radial inflows that move the gas towards the disc centre. However, in Trapp et al. (2022) gas accretion is not mainly connected to a tight galactic fountain structure. The differences in the origin of these features are likely due to the different implementations and strengths of stellar feedback in the simulations. In our simulations, using the SMUGGLE model, stellar feedback primarily drives galactic fountains. These fountains result in gas being ejected to moderate heights above the disc and subsequently falling back, often in coherent, vertically aligned trajectories that concentrate near the disc edge. In contrast, the FIRE-2 simulations used in Trapp et al. (2022) employ a more energetic stellar feedback model that produces stronger outflows that reach much larger distances into the CGM. As a result, the inflowing gas in their simulations is more chaotic in terms of its origins and is not necessarily tied to a galactic fountain cycle.

5. Summary and conclusions

In this work, we analyse the baryon cycle in nine spiral galaxies with different masses and disc sizes, surrounded by a hot CGM, using state-of-the-art N-body hydrodynamical simulations. In particular, we focus on the vertical and radial gas mass flows and on their connection to the structure of the disc. We analyse the spatial distribution of outflows (forming the galactic fountain cycle) and inflows linking the vertical accretion and the formation of stars in the centre through radial flows in the plane of the disc. In doing so, we have employed the moving mesh code AREPO in conjunction with the explicit ISM and stellar feedback model SMUGGLE. We summarise the main results as follows.

• (i)

  We studied the inflows and outflows per unit area divided in cold (T < 10⁴ K), warm (10⁴ < T < 5 × 10⁵ K), and hot (T > 5 × 10⁵ K) gas. We found that the cold phase at |z₀| = 0.5 kpc is dominating both gas accretion and ejection followed by the warm and the hot phases. In each simulation, we find a similar distribution, with a break point roughly corresponding to the disc size R_d, after which the outflow and inflow rates decline rapidly, showing how the outflow-inflow distribution is linked to the structure of the disc.

• (ii)

  We found that gas is accreted vertically onto the disc in a bell-shaped distribution, with the peak skewed towards larger radii (≈0.75 R_d) compared to the SFR distribution. This offset is due to the orbital paths of galactic fountains and their interaction with the CGM, which increases the accreted material. The hot gas exhibits a flatter distribution and appears to be independent from the warm and cold gas, likely resulting from CGM accretion via a (subdominant in terms of mass) hot cooling flow. In contrast, the cold and warm gas distributions are closely linked. The cold gas is predominantly generated by the galactic fountain cycle, while the warm gas forms at the CGM-disc interface owing to the passage of the galactic fountain clouds. This result was also found in Barbani et al. (2023) and confirmed by the simulations in this work.

• (iii)

  We examined the radial gas inflows in the disc. We acknowledge the presence of such net inflows with average values ≈0.7 − 3.8 M_⊙ yr⁻¹ across the nine simulations. The gas is flowing inwardly and outwardly in a periodic trend with an overall inflowing flux in all the simulations over 2 Gyr. The shape of the radial net inflow distribution presents similarities in all the simulations. The mass rate has its peak at ≈R_d, decreasing as it goes inside, owing to the higher gas density and to the formation of stars in the central region of the disc. We found that these radial inflow rates have radial velocity of ≈ − (3 − 6) km s⁻¹ with an almost constant value within R_d. The radial inflows alone are not sufficient to sustain the star formation in the centre, but they can do do in conjunction with vertical inflows which fall at R < 0.4 R_d, where 90% of the SFR is concentrated. Thus, while the gas accreted from the CGM lands preferentially in the outer regions of the disc, the SFR can be fed by radial motions generated by the difference in angular momentum between the CGM and the gas in the disc and by the gas vertically accreted in the centre.

• (iv)

  Gas accretion concentrated in the external regions of the disc causes an inside-out evolution scenario, where the galaxy is growing from the inside to the outside. This is confirmed by the slow growth of both the star and gas disc sizes over time in all nine simulations. In particular, gas and stellar discs grow together with an up-and-down trend linked to stellar evolution. Furthermore, the accretion and ejection of gas are deeply linked to the periodic growth of the disc: peaks in the inflow rate are connected to periods where the disc is growing, while it is the opposite case for the outflow rate.

The analysis we have carried out in this work has enabled us to link together different crucial aspects of gas dynamics in spiral galaxies, collectively known as the baryon cycle. Fresh gas is predominantly accreted from the circumgalactic medium in the outer regions of the disc. This infalling gas has a lower specific angular momentum compared to the gas already present in the disc, leading to an inward flow of gas due to the conservation of angular momentum. At the galaxy centre, where gas density is the highest, a portion of this gas is converted into stars, sustaining the SFR. These newly formed stars impart feedback to the surrounding ISM, eventually driving cold gas clouds out of the disc and forming galactic fountains. As these fountains travel through and interact with the CGM, they create a warm gas phase that eventually cools down and accretes onto the disc at larger radii, thus continuing the baryon cycle. Looking ahead, the inclusion of additional physics, such as magnetic fields, could offer new insights on accretion dynamics. Additionally, studying baryon cycle processes within a full cosmological framework will be crucial for advancing our understanding of spiral galaxies evolution. We defer these aspects to a future work.

Acknowledgments

We acknowledge the CINECA award HP10C9836Y (project IsCb9_MWOUTF) under the ISCRA initiative, for the availability of high performance computing resources and support. We also acknowledge the use of computational resources from the parallel computing cluster of the Open Physics Hub (https://site.unibo.it/openphysicshub/en) at the Department of Physics and Astronomy, University of Bologna. RP acknowledges support from PRIN INAF 1.05.01.85.01. FM acknowledges funding from the European Union – NextGenerationEU under the HPC project “National Centre for HPC, Big Data and Quantum Computing” (PNRR – M4C2 – I1.4 – CN00000013 – CUP J33C22001170001). HL is supported by the National Key R&D Program of China No. 2023YFB3002502, the National Natural Science Foundation of China under No. 12373006, and the China Manned Space Program through its Space Application System. LVS is grateful for partial financial support from NSF-CAREER-1945310 and NSF-AST-2107993 grants.

References

All Tables

All Figures

Fig. 1. Number density of the coronal gas as a function of radius. The black dashed line is the median model obtained by fitting a beta model to Milky Way coronal gas observations (upward pointing triangles with errorbars) using a Markov chain Monte Carlo algorithm. The coloured shaded area represents the 3σ region. Observations, which should be considered as lower limits for the coronal density, are taken from Putman et al. (2021).

In the text

Fig. 2. Gas column density face-on (top panels of each row) and edge-on (bottom panels of each row) projections computed using the SMUGGLE model at t = 1 Gyr for all the nine simulations examined in this work, as labelled in each panel. The simulations differ in the initial value of the scale-length and of the mass of the gaseous disc (see Table 2). Each panel is 70 kpc across and in projection depth with a total number of 1024 × 1024 pixels that give a resolution of ∼70 pc. Brighter colours correspond to higher column densities, as indicated in the colorbar. In the gal_r05, gal and gal_r2 projection panels, the dashed white circle shows the disc size, Rd, defined as the radius at which the gas surface density drops below 1 M⊙ pc−2 (not shown in the other projections for clarity).

In the text

Fig. 3. Star formation rate as a function of time for the nine simulations analysed in this work. The SFR has a wide range of values (≈2 − 10 M⊙ yr−1) across the different simulations, but it remains approximately constant over time for each simulated system. The average SFR depends on the disc initial properties, with variations in each simulation staying within a factor of ∼2 from the average value.

In the text

	Fig. 4. Density-weighted maps of the gas metallicity, radial velocity and temperature (from left to right) in edge-on projections of the fiducial simulation at t = 0.54 Gyr. Each panel is 70 kpc across and in projection depth with a total number of 1024 × 1024 pixels that give a resolution of ∼70 pc. These projections show the stark difference between the disc and the CGM properties.
In the text

Fig. 5. Density-weighted maps of the gas radial velocity in edge-on projections of the nine simulations at t = 0.5 Gyr. Each panel is 70 kpc × 40 kpc and with a projection depth of 70 kpc with a resolution of ∼70 pc. The velocities go from −150 km s−1 (inflows, blue) to 150 km s−1 (outflows, red) as in central panel of Fig. 4. These projections clearly show the presence of galactic fountains in the simulations.

In the text

Fig. 6. Density-weighted maps of the high-metallicity (Z > 0.9 Z⊙) gas radial velocity in edge-on projection of the gal simulation at t = 0.8 Gyr. The velocities range from −150 km s−1 (inflows, blue) to 150 km s−1 (outflows, red) as in central panel of Fig. 4. The projection clearly highlights how the accreted gas is not pristine material coming from the CGM, but is rather a mixture of the latter with gas that was previously ejected from the disc.

In the text

Fig. 7. Outflow rate per unit area as a function of the cylindrical radius, R (in units of the disc size, Rd) for gal_r05 (left panel), gal (central panel), and gal_r2 (right panel) simulations. Each shaded area indicates the outflow rate profiles between the 25th and the 75th percentile, whereas the black dashed lines represent the median values. The rates are divided in cold (T < 104 K, dark blue shaded area), warm (104 < T < 5 × 105 K, blue shaded area) and hot (T > 5 × 105 K, light blue shaded area) gas and computed at a height |z0| = 0.5 kpc inside a slab over the plane of the disc with a thickness Δz = 0.1 kpc and with |vz|> 30 km s−1.

In the text

	Fig. 8. Same as Fig. 7, but for inflow rates.
In the text

Fig. 9. Net vertical inflow rate, computed as Ṁnet = Ṁin−Ṁout, as a function of cylindrical radius R (in units of Rd) for the nine simulations divided in cold (T < 104 K, red line), warm (104 < T < 5 × 105 K, orange line) and hot (T > 5 × 105 K, yellow line) gas. The dashed grey line is the SFR profile. The solid lines represent the profile averaged over the entire simulation and the shaded areas are the standard deviations. The rates are evaluated inside a slab over the plane of the disc with a width Δz = 0.1 kpc at a height |z0| = 0.5 kpc. Each point represents the gas net inflow rate within a circular annulus with a radial width of 0.02 in R/Rd units and considering only gas cells with |vz|> vth = 30 km s−1. In each panel we have indicated the average net vertical rate at R < 0.4 Rd, this is roughly the region of the disc where 90% of the SFR is concentrated in all nine simulations. The net inflows can extend beyond Rd supplying the outer disc with fresh gas for star formation.

In the text

Fig. 10. Specific angular momentum, computed as j = J/M, where J is the total angular momentum and M is the total mass of each respective component, for the gaseous disc and the CGM computed as a function of cylindrical radius R (in units of Rd) for the nine simulations. The angular momentum is measured near the disc (|z|< 1 kpc), with CGM gas selected by Z < Z⊙ and disc gas by Z ≥ Z⊙. In each panel, the grey lines represent the specific angular momentum profile at a given time averaged in a time span of 100 Myr. The black line is the average specific angular momentum and the blue (red) shaded area represents the standard deviation of the disc (CGM) specific angular momentum distribution. Each panel represents a different combination of disc mass and scale length, showing how these parameters influence the angular momentum distribution. Across all models, the disc maintains a systematically higher specific angular momentum than the CGM in the inner regions following the rotation curve, while at R/Rd ≳ 1, the disparity between the two components becomes more variable, marking the transition where the disc fades.

In the text

Fig. 11. Gas net radial mass rate for the simulations analysed in this work. In each panel, the grey lines represent the mass rate profile at a given time averaged in a time span of 100 Myr. The black line is the average mass rate and the blue shaded area represents the standard deviation. On the bottom of each panel the net radial mass rate averaged within the radial range 0 ≤ R/Rd ≤ 1 and the average SFR, computed in the temporal range 0.3 < t < 2 Gyr, are indicated.

In the text

Fig. 12. Gas cylindrical radial velocity for the simulations analysed in this work. In each panel, the grey lines represent the (mass weighted) radial velocity profile at a given time averaged in a time span of 100 Myr. The black line is the average cylindrical radial velocity and the green shaded area represents the standard deviation. On the bottom left corner the radial velocity averaged within the radial range 0 ≤ R/Rd ≤ 1 is indicated.

In the text

	Fig. 13. Radial gas net rate as a function of time for the gal_r05 (dotted line), gal (solid line) and gal_r2 (dashed line) simulations. At each time the radial net rate averaged within the disc, at R < Rd and |z|< 0.5 kpc, is displayed. Overall the radial net rate is mostly positive, indicating a flow of gas towards the centre of the disc.
In the text

Fig. 14. Top panels: Star (orange line), gas (light brown line), and star formation rate (brown line) scale length for the gal_r05 (solid), gal_m05r05 (dashed), gal_m2r05 (dotted, left panel), gal (solid), gal_m05 (dashed), gal_m2 (dotted, central panel), and gal_r2 (solid), gal_m05r2 (dashed), gal_m2r2 (dotted, right panel) simulations as a function of time. These scale lenghts are computed as the radius that encloses 90% of stars/gas mass or SFR. Bottom panel: total outflow (blue) and inflow (light blue) rates of the same simulations of the top panels as a function of time. The simulations with fM = 0.5 are indicated with a dashed line and the simulations with fM = 2 with a dotted line. In the central panels the peak in the inflow rate of the gal simulation are highlighted with dashed grey vertical lines. Each of the scale is connected and they have an up-and-down trend linked to the gas outflow-inflow cycle.

In the text