A Unified Model for the Coevolution of Galaxies and Their Circumgalactic Medium: The Relative Roles of Turbulence and Atomic Cooling Physics

The cosmic gas mass accretion rate into halos is

Here, is the gross inflow rate of DM and baryonic mass. ⁹ f_b = 0.158 is the universal baryon fraction from Planck Collaboration et al. (2016). f_UV suppresses cosmic accretion below the universal value preferentially in lower-mass halos due to photoionization from the cosmic UV background. To estimate f_UV, we first use Appendix B of Kravtsov et al. (2004) to compute the redshift-dependent “filtering halo mass” at which 50% of baryons are prevented from accreting and then Equation (1) of Okamoto et al. (2008) to compute f_UV depending on the ratio of the halo’s current mass to the filtering mass. The filtering mass steadily increases from ∼10⁸ M_⊙ at z ∼ 4 to ∼8 × 10⁹ M_⊙ at z ∼ 0 so that f_UV → 0 for progressively lower-mass halos at later times. The f_prev parameter accounts for any additional suppression of baryon accretion due to pre-heating, SN feedback, etc. In general, as we will show later, we adopt f_prev ≈ 1 for MW-mass halos so that they experience no suppression but ≈0.3 for dwarfs such that they accrete only 30% of the universal baryon fraction times the total mass accretion rate.

We assume that halo gas accretion brings in a total energy

associated with gas free-falling into the halo. Here, μ = 0.59 is the mean molecular weight and T_vir is the halo virial temperature:

with

being the virial velocity of the DM halo.

We introduce a free parameter that partitions this inflowing energy into thermal versus kinetic forms:

In this subsection we describe the details of how we model the thermal and turbulent structural components of the CGM as well as the gas cooling and accretion rate into the ISM. The relevant CGM structural parameters are illustrated in the inset panel of Figure 1.

2.2.1. Thermal Component

We assume flexible power-law profiles for the CGM number density and temperature:

where n₀ and T₀ are the density and temperature, respectively, at R_vir. The power-law slopes α_n and α_T are free parameters—we typically fix α_n = − 3/2 and α_T = 0 corresponding to a steady-state cooling flow solution in an isothermal potential (Stern et al. 2019), and also very close to what is assumed in our companion paper (Carr et al. 2023). The normalizations are then derived by requiring that the integrals of the profiles match the state variables M_CGM and at a given time (the lower limit of 0.1R_vir is somewhat arbitrary but consistent with our assumed CGM–galaxy boundary for computing fluxes in the FIRE-2 simulations in Section 4):

These can be solved analytically to get

Note that, contrary to the assumptions of traditional SAMs, T₀ is generally not equal to the halo virial temperature T_vir. In principle, we could also similarly set up profiles for the CGM metallicity and turbulent velocity. However, for simplicity we assume that the CGM metallicity is constant with radius. We also adopt a single global characteristic turbulent velocity since introducing a radial dependence for v_turb would require additional complications in order to relate it to our assumed model for the circular velocity profile of the underlying DM halo.

Then the radiative cooling rate for the thermal component is

where Λ(r) is the cooling function. We adopt the Wiersma et al. (2009) cooling tables, which assume collisional ionization equilibrium (CIE) plus photoionization from the Haardt & Madau (2001) UV background, thus causing Λ to depend on density, temperature, metallicity, and redshift. The cooling tables only go up to z = 8.989—above this redshift (or for any other extrapolation) we set Λ = 0. These cooling tables include photoionization heating of gas that is already in the CGM whereas our prescription for f_UV above inhibits the accretion of gas into the halo in the first place (see also Benson et al. 2002). Note that our model does not account for the possibility that cold clouds condense out of the warm/hot phase thereby lowering its density and cooling rate, nor for the scenario that these cold clouds can contribute energy and mass back to the warm/hot phase as they freefall toward the galaxy (Maller & Bullock 2004; Murray & Lin 2004; Forbes & Lin 2019; Fielding et al. 2020a).

2.2.2. Turbulent Kinetic Component

The turbulent kinetic energy of the CGM is assumed to dissipate on a timescale (Stone et al. 1998; Mac Low 1999)

so that

t_turb is known as the characteristic eddy turnover time, R_turb is the characteristic size of the largest turbulent eddies, and v_turb is their corresponding turbulent velocity (this can also be thought of as the specific turbulent kinetic energy of the CGM). We take

We do not have a priori knowledge of what R_turb should be and how it may vary with halo mass, redshift and/or CGM properties. Thus, we make the simple assumption that

In other words, we assume that the sizes of the largest turbulent eddies in the CGM are proportional to the halo virial radius and that this size sets the timescale on which the turbulent cascade proceeds. The proportionality factor may or may not be a function of time depending on the driving sources of the turbulence. In Section 4.3.5 we will present our fiducial function that relates R_turb and R_vir as a function of time.

2.2.3. CGM Mass Cooling and Inflow Rate

So far we have dealt with CGM energetics, but now we must consider the rate at which the CGM loses mass via accretion into the ISM. We compute the ISM accretion rate as

Here, t_cool,eff is the “effective” timescale on which the CGM would radiate away its thermal energy given the net cooling rate (Equation (21)) minus the turbulence dissipation rate:

There is another limiting timescale for ISM accretion that we call the “effective” freefall time:

After CGM gas has had sufficient time to cool (or while it is cooling), it will take this additional amount of time to freefall into the galaxy from some characteristic CGM radius, which we take to be . ¹⁰ This is the radius where the halo circular velocity profile reaches its maximum with R_s = R_vir/c being the scale radius of a Navarro–Frenk–White (NFW) halo with concentration c. Thus, in the absence of turbulent pressure support, the cooling gas would inflow on a timescale . However, Equation (28) multiplies this by a smooth function that increases the inflow timescale when the ratio of the turbulent pressure support () to the local gravitational potential energy density () is large, and reverts to the baseline freefall time otherwise. The exact functional form is arbitrary but achieves our desired effect in a smooth manner. To estimate , we use an analytic approximation for the circular velocity profile of an NFW halo with a given concentration and V_vir (Equation (11.26) from Mo et al. 2010, with x ≡ r/R_vir):

Another way to compute would be to radially integrate over 4π r² ρ(r)/(t_cool,eff(r) + t_ff,eff(r)). Alternatively, we could compute a “cooling radius” R_cool within which t_cool,eff < t_ff,eff and then set following the approach originally laid out by White & Frenk (1991) and later adopted by most SAMs. While both of these methods would naturally take into account the fact that cooling and freefall times are longer at larger radii and negate the need for arbitrarily choosing a single characteristic freefall radius, it would introduce additional complications for parameterizing a v_turb(r) profile and relating it to the underlying V_circ(r) profile of the DM halo to compute t_ff,eff(r). Instead here we adopt the simpler approach of treating the CGM cooling and freefall times with single global parameters.

We model the small-scale physics of star formation and SN-driven winds in a very simple, flexible way. The star formation rate is modeled using a single depletion time parameter for the entire ISM mass:

The SN-driven mass outflow rate is simply

where η_M is the wind mass loading factor. The corresponding energy outflow rate of the SN-driven wind is

where

is the Bernoulli velocity that quantifies the specific energy of the wind. Here, v is the wind velocity, and c_s is the wind sound speed corresponding to its temperature. Note that

where erg/(100M_⊙) is the SN energy produced per 100M_⊙ of stars formed, consistent with the Kroupa (2001) initial mass function (IMF). Thus, we will use v_B and η_E/η_M interchangeably to refer to the wind specific energy. Finally we introduce a free parameter that partitions the outflowing wind energy into thermal and kinetic forms (analogous to how we partitioned the cosmic accretion energy):

It follows from Equation (33) that corresponds to a Mach number of one. Various functional forms can be adopted for the free parameters t_dep, η_M, v_B (or equivalently η_E), and . We will present our fiducial functional forms in Section 4.

The energy pumped into the CGM from SN winds and cosmic accretion may over-pressurize it compared to its own binding energy ¹¹

This overpressurization may happen if radiative cooling cannot keep up with heating and/or if the turbulent energy is not dissipating fast enough. In this case we envision that some fraction of the CGM will become unbound from the halo and drive outflows at R_vir to naturally decrease the level of overpressurization (see also Carr et al. 2023)

The total excess energy outflow rate at R_vir is

where t_dyn,halo = R_vir/V_vir as usual. ¹² For the thermal versus kinetic energy partitioning of the halo outflows, , we assume the current thermal energy support fraction of the CGM (note that this is not a free parameter):

so that

Finally the mass outflow rate at R_vir is just the total energy outflow rate divided by the specific energy of the halo outflows, e_out,halo, which by default we take as equal to the specific energy of the CGM:

Because we are forcing the halo energy outflows to have the same thermal versus kinetic split as the CGM, halo outflows do not change the CGM temperature. However, this can be modified by introducing another multiplicative free parameter to rescale the denominator to much higher or lower specific energies relative to the CGM.

We use the instantaneous recycling approximation to track the production and flow of all metals combined (i.e., we do not track different species or production channels). We follow Section 10.4.2 of Mo et al. (2010), which is exactly equivalent to the ODEs implemented in some SAMs (Cole et al. 2000; De Lucia et al. 2004; Lagos et al. 2018). Following Tinsley (1980), we assume that new stars form with the same metallicity as the ISM (Z_ISM) so that

Note that Equations (7) and (8) already include the (1 − f_rec) factor to account for the instantaneous return fraction of ISM metals via stellar winds and SNe. The addition of new metals to the ISM from stellar nucleosynthesis is given by

The parameter y_Z is defined as the ratio of new metal mass ejected from stars divided by the total mass locked up in low-mass stars and stellar remnants. Thus, by definition, y_Z includes a factor of (1 − f_rec) in its denominator (see Equation (10.120) of Mo et al. 2010). Note that some SAMs define y_Z differently by excluding the factor of (1 − f_rec) in its denominator and then also removing it from Equation (44) (Somerville et al. 2015; Lagos et al. 2018), but this is identical to our approach. We assume (1 − f_rec)y_Z ≈ 0.02, consistent with a Kroupa (2001) IMF.

The metal accretion rate into the halo is assumed to be

where Z_in,halo is a free parameter for the inflow metallicity. For pristine accretion, this would be Z_in,halo = 0 although we expect some pre-enrichment from recycling of galactic winds and/or outflows from nearby halos. Note that some SAMs model a separate “ejected” component into which they deposit winds and metals ejected from halos, and then compute the recycling rate of this ejected gas back into the halo (Henriques et al. 2013; White et al. 2015). However, this introduces additional complications such as wind escape fractions, mass-dependent recycling times, and contribution of additional metals beyond those produced by the central galaxy that we do not include in this work.

The metal accretion rate into the ISM from CGM cooling is taken to be proportional to the CGM metallicity:

This assumes perfect mixing of inflows and outflows in the CGM and ignores the possibility that metal-enriched material is expected to cool more efficiently (Hobbs et al. 2015).

The metal mass carried by SN-driven winds is assumed to be proportional to the ISM metallicity:

In principle the wind metallicity could be different from the ISM metallicity if, for example, there was minimal entrainment and the wind contained mostly SN ejecta. This would require introducing a new wind enrichment factor or metal loading factor (e.g., Kim et al. 2020; Sharda et al. 2021; Carr et al. 2023), but we defer this complication to future work. Finally, we assume that the metal outflow rate from the halo is directly proportional to the CGM metallicity:

Table 1 lists all of our model parameters and whether they are fixed or allowed to vary. In total our model has 16 parameters governing cosmic accretion, CGM structure, star formation, and stellar feedback, and chemical evolution. However, most of these parameters except one (R_turb) are either fixed based on simple physical arguments or on our calibration to the FIRE-2 simulations (see Section 4). In many cases, a parameter is not simply a constant but rather is given by a functional form with a number of additional arguments that quantify the redshift and/or halo mass dependence of that parameter. In these cases, we refer the reader to the relevant sections quoted in the table for the exact parameterizations.

Table 1. List of Model Parameters

Parameter	Fixed/Free	Value	Meaning
Cosmic accretion			Sections 2.1, 4.3.3, 4.3.4
fUV	Fixed	Okamoto et al. (2008)	Suppression factor for cosmic accretion due to UV background
fprev	Fixed*	Logistic	Suppression factor for cosmic accretion due to preventative feedback
	Fixed*	Power law	Fraction of accretion energy that is thermal rather than kinetic
CGM structure			Sections 2.2 and 4.3.5
αn	Fixed	−1.5	Slope of CGM density power law
αT	Fixed	0.0	Slope of CGM temperature power law
Rturb	Free	Logistic	Sizes of largest CGM eddies (controls turbulence dissipation rate)
rff	Fixed		Radius at which to compute the effective freefall time of the CGM
Star formation and stellar feedback			Sections 2.3, 2.4, and 4.3.1
tdep	Fixed*	Power law	ISM depletion time
ηM	Fixed*	Power law	Mass loading of SN-driven winds
vB	Fixed*	Power law	Specific energy of SN-driven winds
	Fixed*	Logistic	Fraction of SN wind energy that is thermal rather than kinetic
	Fixed		Fraction of halo outflow energy that is thermal rather than kinetic
eout,halo	Fixed	ECGM/MCGM	Specific energy of halo outflows
Chemical evolution			Sections 2.5 and 4.3.2
frec	Fixed	0.4	Instantaneous stellar recycling fraction
(1 − frec)yZ	Fixed	0.02	Nucleosynthetic yield from star formation and stellar evolution
Zin,halo	Fixed*	Power law	Metallicity of cosmic accretion

Notes. Fixed* means the parameter was fixed based on FIRE-2. If the value is not a constant and is instead given as a reference or functional form, then there are a number of additional arguments that control the redshift and/or halo mass dependence of the parameter (see the relevant section for the exact parameterization). In total we have 16 parameters, but effectively all of these except one (R_turb) are fixed based on the FIRE-2 simulations or simple physical arguments.

Download table as:  ASCIITypeset image

We solve the system of ODEs defined in Equation (1) using the Python scipy.integrate.solve_ivp ODE solver. Specifically, we use the adaptive-time-step implicit method “BDF” (backward differentiation formula), which gives identical results compared to the adaptive-time-step explicit methods “RK23” and “RK45” (Runge–Kutta) but with fewer overall iterations required. Our results are insensitive to the choice of initial conditions for the state variables as long as they are reasonably small and close to, but not exactly, zero (this is because the very early evolution is driven by the rapid cosmic assembly of halos).

Many of the individual terms in our ODEs require external time-dependent inputs related to the halo assembly history. These include the halo mass, halo radius, halo concentration, and gross DM accretion rate. Since we use measurements of these halo properties from numerical simulations (see Section 4) and since those measurements can be noisy (the typical spacing between simulation outputs is ∼2–10 Myr), we need to smooth and interpolate those input time series with a sufficiently high-order function to ensure optimally adaptive-time-stepping and stability for our ODE solver. For this, we make use of scipy.interpolate.UnivariateSpline to fit a smoothing spline of degree 5 with a smoothing factor of 2. Furthermore, when comparing our predicted time series for any given quantity to the measurements of that property from simulations, we smooth the simulation time series with a Gaussian whose standard deviation is 10 (in units of number of adjacent data points for simplicity) using scipy.ndimage.Gaussian_filter1d. This is done because our model is based on ODEs, which are inherently smooth and do not account for time delays and stochasticity.

T₀ and v_turb, respectively, quantify the specific thermal energy and specific kinetic energy of the CGM in our model (see Equations (20) and (24), respectively). These two parameters hence reflect the amount of cooling, heating, and overpressurization in the CGM, which in turn will affect star formation and supernova feedback. Since these two parameters are also novel compared to previous SAMs (which fix T₀ = T_vir and neglect CGM turbulence), it is important to understand their equilibrium values (where and ) since our model will tend to evolve toward those values.

As an illustrative example, we adopt model parameters that are reasonable for an MW-mass halo at z = 0:

1. 1.  
   M_vir = 10¹² M_⊙, R_vir = 275 kpc, c_NFW = 15. These together give V_vir = 125 km s⁻¹ and T_vir = 5.6 × 10⁵ K.
2. 2.  
   yr⁻¹
3. 3.  
   M_CGM = 5 × 10¹⁰ M_⊙ and Z_CGM = 0.3Z_⊙
4. 4.  
   α_n = − 3/2 and α_T = 0.0
5. 5.  
   R_turb = 0.5R_vir
6. 6.  
   and
7. 7.  
   η_M = 0.5 and v_B = 300 km s ⁻¹ which together imply
8. 8.  
   e_out,halo = e_CGM.

Then we set up a 2D grid spanning a range of v_turb and T₀ values. At each point in this grid, we recompute and , which in turn will change the instantaneous , , and . Since we also want to see the effect of changing v_turb and T₀ on the SFR and winds, we recompute . This is an approximation because we are ignoring the ODEs linking to and hence not self-consistently predicting the SFR. However, setting the SFR equal to the ISM accretion rate should be appropriate for our case of an MW-mass halo at z = 0 (more generally, in lower-mass halos and high-redshift MW progenitors, we expect the SFR to only be a small fraction of the ISM accretion rate). Finally, we are then in a position to compute and across our grid.

Figure 2 illustrates how and vary across our v_turb − T₀ grid. The narrow strips where and , respectively, pinpoint the range of possible equilibrium values of T₀ and v_turb for our choice of z = 0 MW-like parameters. The equilibria are generally stable: for small perturbations away from the strips, the model will be pushed back toward equilibrium. The exception is low-temperature solutions (T₀ ≲ 5 × 10⁵ K) around which , indicating that those equilibria are unstable and that small perturbations would drive the model to a hotter temperature. The intersection of the two equilibrium regions for both and gives the combination of equilibrium T₀ and v_turb that the full model will tend to evolve toward (and once reaching these values, the model will stay there except for forcing terms from cosmological accretion). For our choice of z = 0 model parameters, this intersection happens at T₀ ≈ 8 × 10⁵ K (slightly hotter than T_vir) and v_turb ≈ 130 km s⁻¹ (comparable to V_vir).

Zoom In Zoom Out Reset image size

Having illustrated the existence of equilibrium solutions in our model, a natural follow-up question is to ask how these equilibria depend on our choice of model parameters. Figure 3 is similar to Figure 2, but now we vary a few model parameters alongside T₀ while fixing v_turb = V_vir to see the effect on , and then repeat the variation alongside v_turb while fixing T₀ = T_vir to see the effect on . The two model parameters that we explore are the wind specific energy and the largest eddy turnover scale. We expect our model to be quite sensitive to these two parameters (in addition to redshift, halo mass, CGM metallicity, etc., but we leave those for Section 5.6, where we will show how model equilibria evolve naturally along individual halo assembly histories). The following subsections examine each of these parameter variations in turn.

Zoom In Zoom Out Reset image size

3.2.1. Wind Specific Energy

3.2.2. Largest Eddy Turnover Scale

We use the “core” suite of the FIRE-2 cosmological hydrodynamical “zoom-in” simulations (Hopkins et al. 2018). Our sample includes three ultrafaint dwarfs (m10q, m10y, and m10z, which have M_vir ∼ 10¹⁰ M_⊙ at z = 0), six intermediate-mass dwarfs (m11a, m11b, m11c, m11q, m11v, and m11f, which have M_vir ∼ 10¹¹ M_⊙ at z = 0), and three MW-mass halos (m12i, m12f, and m12m, which have M_vir ∼ 10¹² M_⊙ at z = 0). This is the same sample that we used in Pandya et al. (2020, 2021) except we exclude the anomalously late-forming low-mass dwarf m10v, which has effectively no star formation until z ≲ 0.5 (see Section 2 of those papers for more details and references about the simulations). We use the same Rockstar halo catalogs and consistent-trees merger trees (Behroozi et al. 2013a, 2013b) that we generated and described in Pandya et al. (2020). We adopt the Bryan & Norman (1998) definition of virial overdensity.

As in Pandya et al. (2020), we split each FIRE-2 halo into two zones: gas and star particles within 0.1R_vir are used to compute the stellar and ISM mass and metallicity of the central galaxy in each snapshot whereas the CGM mass, thermal energy and kinetic energy are computed using all gas particles between 0.1 and 1.0 R_vir. We do not attempt to exclude gas associated with satellites orbiting in the CGM: both their cold gas and hot outflows can affect the global thermodynamics of the host CGM. The global CGM thermal energy is computed by summing over the thermal energy of the individual gas particles:

where the subscript i runs over all CGM gas particles, and m_i and T_i are, respectively, the mass and temperature of the ith particle. Similarly the global CGM kinetic energy is computed as

where v_i ² is the norm of the halo-centric particle velocity vector. Note that we do not attempt to decompose the CGM velocity field into a turbulent component so our measurement of includes contributions from bulk and rotating flows. Hence we are measuring an upper limit to the actual turbulent kinetic energy that we care about for comparison to our model. Finally, the stellar, ISM and CGM metal masses are computed very similarly but summing over only the metal mass fractions of the relevant particles. The metallicities are then computed as the ratio of the total zone metal mass divided by the total zone mass further normalized to solar metallicity Z_⊙ = 0.02.

We also compute gas inflow and outflow rates at the galaxy and halo scale following Pandya et al. (2020, 2021). Complementing those earlier papers, we have implemented single-adjacent-snapshot particle tracking, which we will extend in the future to operate over multiple snapshots for the purpose of characterizing halo outflow recycling fractions and timescales (adapting the methodology of Anglés-Alcázar et al. 2017a; Hafen et al. 2019). After classifying particles into different zones as described above for any two adjacent snapshots, we compute the intersection of the particle ID arrays to identify which particles crossed zones. For example, inflowing particles at the halo scale are those that crossed from >R_vir in the first snapshot to <R_vir in the later snapshot. We can then estimate the mass inflow rate into the CGM by adding up the masses of all of the crossing particles and dividing by the time step between these two adjacent snapshots (see also Wright et al. 2020). We repeat this procedure to compute the mass outflow rate from the CGM at R_vir and the mass inflow and outflow rates at 0.1R_vir. The gross DM accretion rate at R_vir is also computed in this way.

We additionally compute the metal mass flow rates by summing only over the metal mass fractions of the crossing particles at both interfaces. We also compute the energy transported by galactic winds at 0.1R_vir and at R_vir following Pandya et al. (2021) by multiplying the mass of each crossing particle by its Bernoulli velocity squared, which is a measure of the energy available to drive outflows. We also compute the energy inflow rate at the galaxy and halo scale using the same v_B definition but for particles that are crossing inward. Finally, since is the sum of the specific kinetic energy flux and the specific enthalpy flux, we can compute the fraction of the energy flux that is already thermalized by taking the ratio of the specific enthalpy flux to the Bernoulli velocity squared ( following Equation (33)). This constrains and for our model.

Our results are very similar to the mass fluxes reported in Pandya et al. (2020, 2021) with the caveat that here we do not impose any cut on v_B for outflowing particles, unlike in Pandya et al. (2021). Thus our outflow rates correspond to simply using a v_rad > 0 km s, ⁻¹ threshold. This is done to ensure that when we assume the FIRE-2 wind loading factors in our model, we will conserve mass during the gas flow cycle. Otherwise, a significant cut on v_B for winds as in Pandya et al. (2021) may miss up to half of the outflowing gas going into the inner CGM, and this would make it impossible for our model to reproduce the time evolution of the CGM and ISM mass measured in the simulations. As a consistency check, we have verified that integrating our measured mass fluxes gives a nearly identical CGM and ISM mass as a function of time as measured from the particle data. This exercise increases our confidence in our fluxes and bulk measurements in the sense that if our model can accurately reproduce the fluxes as a function of time, then it should naturally also reproduce the time evolution of the integrated CGM and ISM properties.

4.3.1. Star Formation and SN-driven Winds

In our model, the small-scale physics of SF and SN-driven winds is represented by three parameters: the depletion time of the entire ISM mass (t_dep), the wind mass loading factor (η_M), and the wind specific energy as quantified by the Bernoulli velocity (v_B). Figure 4 shows our fit to these parameters using the FIRE-2 data. We parameterize t_dep as

We parameterize η_M as

and the Bernoulli velocity as

Although we could equivalently parameterize η_E instead of v_B, we find that v_B can more easily be fit with a simple power law whereas η_E requires a more complicated fit (see also Figure 9 of Pandya et al. 2021).

Zoom In Zoom Out Reset image size

4.3.2. Chemical Evolution

Recall that our model has two free parameters governing chemical evolution: (1) the nucleosynthetic yield, which we take to be (1 − f_rec)y_Z = 0.02 as appropriate for a Kroupa (2001) IMF, and (2) the metallicity of gas flowing into the halo. In principle we could also introduce another free parameter for the wind metal enrichment factor, but we avoid that complication and assume winds have the same metallicity as the ISM. This may be particularly justified for FIRE-2 for which we found that the winds were so heavily mass-loaded that their average metallicity probably tracks that of the entrained ISM (Pandya et al. 2021). Figure 4 shows how Z_in,halo scales with V_vir and redshift for the FIRE-2 halos. We find that

The halo inflow metallicities in FIRE-2 are nonzero and depend strongly on both V_vir and redshift. The nonzero metallicities likely reflect both recycling of galactic winds on large scales as well as chemical pre-enrichment by neighboring halos. The redshift dependence at fixed V_vir can be understood as there simply being more metals produced and hence available to be recycled at later times. The V_vir dependence likely arises for multiple reasons: (1) wind metallicities are higher in more-massive halos as shown in the bottom panel, (2) more-massive halos have higher SFRs and hence higher wind mass outflow rates even though their mass loading factors are lower, and (3) more-massive halos have more satellites that can pollute their local environment before accretion. Still, it is striking that in FIRE-2, the MW-mass halos at z ∼ 0 have inflowing metallicities at R_vir of the order of ∼0.1Z_⊙.

4.3.3. Preventative Feedback for Halo Gas Accretion

The leftmost panel of Figure 5 shows our preventative feedback parameter f_prev from Equation (9). The data points are measurements of the ratio at R_vir for individual FIRE-2 halos at z = 0 and z = 3. We use a simple generalized logistic function to approximately capture this trend:

In Carr et al. (2023), we self-consistently predict f_prev using the ratio of E_out,halo to E_in,halo (see also Section 5.2 of Pandya et al. 2020). In order for that model to capture the complexity of FIRE-2, we would need to introduce a free parameter that quantifies how much of the outflowing energy couples to and heats the accreting gas. This parameter may have a complicated dependence on halo mass, redshift, etc. with scatter. Thus we opt for a simpler parameterization in this work.

Zoom In Zoom Out Reset image size

4.3.4. Thermal versus Kinetic Flux Partitioning

The middle two panels of Figure 5 parameterize the fraction of cosmic accretion energy and SN wind energy that is thermal for the FIRE-2 halos. We assume

where captures the redshift dependence of the slope. At low V_vir, this power law can lead to values larger than one; in this case, we simply set to prevent negative . The decrease in with V_vir and redshift reflects the importance of cold filamentary accretion in the FIRE-2 halos.

For we use a generalized logistic function:

It is interesting that ultrafaints have thermally supported winds whereas classical dwarfs and MW-mass galaxies all cluster around , indicating efficient mixing and/or cooling. In the Appendix, we will show how our predictions change if we set and so that we are in the purely thermal limit of our model (as in our companion paper, Carr et al. 2023).

4.3.5. Characteristic Size of Largest Turbulent Eddies

The characteristic size of the largest turbulent eddies in the CGM is one of the main uncertainties of our model since we do not have existing constraints on this observationally or from cosmological simulations. The rightmost panel of Figure 5 shows our assumed generalized logistic function for R_turb(t) normalized by R_vir(t):

where again is the radius of as used in Equation (28). The slope (2) and pivot (7 Gyr) of our logistic function were arbitrarily chosen and control when and how quickly R_turb drops from . Note that if we replace with R_vir, then we would simply have R_turb(t) = R_vir(t). We imagine that at early times, large-scale structure formation continuously drives turbulence such that , which is achieved as R_turb → ∞ . On the other hand, at later times, we assume that SN winds are the primary drivers of turbulence and that this occurs mainly in the inner halo such that . Our model is quite sensitive to the exact choice of these asymptotic values—smaller R_turb at early times leads to faster turbulence dissipation and thus less turbulent pressure support, higher accretion rates into the ISM, and more star formation.

Figure 6 shows the time evolution of DM, CGM, ISM, and stellar mass for the 12 individual core FIRE-2 halos as measured from the particle data and as predicted by our model. The model trajectories generally follow the time series of mass measurements from the FIRE-2 particle data to within a factor of a few. However, there are exceptions: the CGM mass tends to be systematically underestimated by a factor of 2 in the model relative to FIRE-2, and the stellar and ISM masses of some dwarfs are larger than FIRE-2 by a factor of up to ∼10. Some of these discrepancies may be due to our simple average parameter fits that do not capture scatter between individual halos (Figures 4 and 5). More importantly, for parameters that could not easily be measured from the simulations such as the turbulence dissipation timescale or an effective freefall radius for computing the ISM accretion rate, our choices may not reflect what is happening in the simulations. By comparing the flow rates of mass, energy, and metals between the model and FIRE-2 in the next few subsections, we can learn more about the behavior of the new model and possible causes for the discrepancies in the mass assembly histories.

Zoom In Zoom Out Reset image size

Since we roughly reproduce the time evolution of the masses of the different components of the FIRE-2 halos, it is natural to ask what our model implies for the baryon fractions of halos. Figure 7 shows the average halo baryon fraction since z = 2 as measured from the particle data and as predicted by our model. We have verified that the baryon fractions are roughly constant since z = 2, so taking the average is a good summary statistic. We reproduce the trend in FIRE-2, which is that dwarfs have ≲10% of the universal baryon fraction f_b = 0.158. In our model, this happens because of two things: (1) the preventative feedback parameter f_prev suppresses cosmic accretion, and (2) the CGM can become overpressurized, which ejects previously accreted baryons. Since f_prev provides a rough upper limit on the halo baryon fraction, having f_b < f_prev indicates that CGM overpressurization is important, and this is indeed the case for the lowest-mass halos in our model.

Zoom In Zoom Out Reset image size

Matching the evolution of bulk masses in itself is not sufficient to claim that our model reproduces FIRE -2 since there are many different ways to get to the same mass (Pandya et al. 2020). Figure 8 shows that we also roughly reproduce the time evolution of the individual mass flow rates underlying the and time derivatives. Since we parameterized f_prev and η_M from FIRE-2, it is perhaps not surprising that we reproduce the cosmic halo accretion rate and ISM wind mass-loss rate. However, the cooling rate and halo outflow rate are genuine predictions of our model and match FIRE-2 quite well. The main exception is that our cooling rates tend to be on the higher side for the dwarfs, but these are sensitive to our predictions for the turbulence dissipation rate and choice of an effective freefall radius for computing . It is possible to achieve better agreement with FIRE-2 by increasing R_turb or r_ff since the former would provide more turbulent pressure support (because the turbulence would decay more slowly) and the latter would increase the effective freefall time of CGM gas, in turn lowering , SFR, , and probably also . But even with our fiducial choices, the results of the simple model show promise for reproducing mass flows in FIRE-2.

Zoom In Zoom Out Reset image size

Having shown that our model is capable of roughly reproducing the mass budgets and mass flow rates of FIRE-2, we now turn to the metal budgets and metal flow rates. Metals are an important additional dimension to predict because they will unlock more observables that we can use to eventually test our model.

Figure 9 compares the time evolution of metal flow rates from our model to FIRE-2. We see very good agreement, as with the overall mass flow rates. Figure 10 plots the bulk metallicities of our three components (CGM, ISM, and stars) as a function of halo mass in comparison to FIRE-2. We show the comparison for z = 0 and a representative high redshift z = 3. The model tracks the general trends of CGM, ISM, and stellar metallicity with halo mass at both low and high redshift, which is noteworthy given the simplicity of our chemical evolution prescriptions. One discrepancy is that the CGM metallicity of the model MW halos is ∼4× higher than FIRE-2 whereas the ISM and stellar metallicities agree to better than a factor of 2. The high CGM metallicity of the MW halos at z = 0 partially reflects the CGM mass being a factor of ∼2 lower in the model compared to FIRE-2 (see Figure 6) with the remaining excess possibly attributable to our assumption that the wind metallicity equals the ISM metallicity (it may be lower if there was entrainment of inner CGM gas or if the ISM is inhomogeneous and winds are launched from less-enriched regions). The other main discrepancy is that the lowest-mass halos at z = 3 have higher metallicities in the model compared to FIRE-2. This may reflect a breakdown of the instantaneous recycling approximation at high redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 11 shows the time evolution of the energy inflow rate from cosmic accretion and SN-driven winds as well as the halo-scale energy outflow rate due to CGM overpressurization. Here again the inflow rate of energy from cosmic accretion agrees to better than a factor of 2 between the model and FIRE-2, whereas the wind and halo energy outflow rates show some discrepancies. First, while the wind energy outflow rates for most halos agree with FIRE-2 to within a factor of a few, the lowest-mass dwarf shows an excess in the model. This directly follows from its excess CGM cooling rate and excess SFR, which could be resolved by increasing the turbulence dissipation time and/or effective freefall time to limit the ISM accretion rate as discussed in Section 5.3. Second, the halo energy outflow rates also agree with the FIRE-2 measurements to within a factor of a few except at very early times where the model tends to have no halo outflows. This lack of early halo outflows in the model may be due to additional contributions from other sources that are not included (e.g., more turbulence driven by satellite motions). Despite these differences, it is encouraging that the model follows the general trends of the simulations.

Zoom In Zoom Out Reset image size

Analogously, Figure 12 shows the time evolution of CGM thermal energy, kinetic energy, and binding energy both as predicted by our model and as measured in the FIRE-2 particle data for three representative halos. Our model roughly reproduces the trends measured from the FIRE-2 data. All three halos show a transition from an early, cool, kinetic-dominated phase to a thermally supported warm/hot halo at later times. In the low-mass dwarf, this transition happens very early since cooling can only balance heating from UVB photoionization, SN winds, turbulence dissipation, and cosmic accretion at a slightly supervirial temperature. The intermediate-mass dwarf transitions more gradually and at lower redshift since it is at the peak of the cooling curve and hence can offset more of the SN/turbulent heating. The MW-mass halo transitions abruptly at an intermediate redshift in our model but more gradually in FIRE-2. As shown by Stern et al. (2021), the transition in FIRE-2 is radially dependent with the outer CGM undergoing the transition first. This complexity is not reflected in our model, which treats the entire CGM as a single zone (thus leading to a sharper transition) and is also not captured by our measurements of a single global for the FIRE-2 halos in Figure 12 (we will discuss this further in Sections 6.1 and 6.3). Nevertheless, our model does capture the essential idea that the CGM may undergo a “thermalization” process.

Zoom In Zoom Out Reset image size

Figure 13 compares the average global temperature and turbulent velocity of the CGM to the halo virial velocity for our halos at both low and high redshift. At late times, the CGM temperature roughly tracks the virial temperature of the dark matter halo as expected, but at early times the two are decoupled with T₀ ≪ T_vir. The turbulent velocity is comparable to the halo virial velocity, or even exceeds it, in the dwarfs at both low and high redshift. In the case of the MW-mass halos, this is only true for their progenitors: by late times the turbulent specific energy decays to become only a fraction of the halo specific energy.

Zoom In Zoom Out Reset image size

We saw in the previous subsection that the CGM in our model transitions from an early, cool, turbulent phase to a warm/hot thermally supported phase at late times. Here we analyze the time evolution of our model equilibria to provide insight on why this CGM phase transition happens. This builds on Section 3 where we showed how to analyze the equilibria of our model using idealized state variables and parameters for a z = 0 MW-mass halo. The difference is that here the state variables such as CGM mass and metallicity are self-consistently evolved in time, and we have also determined our model parameters from FIRE-2 or simple physical arguments.

The top panel of Figure 14 shows the time evolution of v_turb and c_s (or equivalently T₀) predicted by our model when we run it on the merger tree of a single MW-mass FIRE-2 halo (m12m). Also shown are the dependence of and on varying T₀ and v_turb, respectively, without any other parameter variations. We urge the reader to download our supplementary movie showing the time evolution of this figure, which otherwise shows only a single annotated frame at z = 0.

We see that at early times, the CGM temperature in our model is quite cool (T₀ ∼ 10⁴–10⁵ K) but turbulent (v_turb ≳ V_vir). By z = 0, the CGM temperature rises to ≳T_vir and the turbulent velocity has decayed to <V_vir. Thus this is another way to visualize the CGM phase transition that we described in the previous section, with the phase transition in this case happening abruptly at t ≈ 7 Gyr (z ≈ 0.75). At z = 0 as shown in this frame, there is only one equilibrium for T₀ and v_turb and the model lives very close to it. The T₀ ≈ 10⁶ K equilibrium is set by a balance between the thermal energy sink terms (primarily though also plays a role) and the heating source terms (dissipation and accretion dominate over the wind thermal energy at z = 0). Similarly, the v_turb ≈ 90 km s⁻¹ equilibrium is set primarily by a balance between turbulence dissipation (sink term) and turbulence driving by SN winds and cosmic accretion.

In the movie version of Figure 14 (see caption for download link), we can compare the time evolution of the model to its equilibria. In general, the equilibria themselves evolve with time, and the model generally follows this secularly evolving equilibrium. For example, at early times there is only a single equilibrium temperature, and it steadily increases due to the forcing terms from cosmological accretion and the increasing star formation rate. However, just before and after the phase transition, we see bifurcations in the solution to our coupled ODE system: multiple equilibria can appear and disappear depending on the shape of the cooling function and where it intersects with the sum of the various heating terms. Just before the phase transition, the original cool equilibrium vanishes since the heating terms exceed the available cooling, and the model is then forced to quickly evolve toward the next remaining hotter equilibrium. These bifurcations are largely driven by the fact that the drops steadily with time up until the phase transition. Since , this steady drop in cooling reflects the gradual decrease in the mean cosmic density as shown in Figure 15. We have verified that both Λ and Z_CGM remain roughly constant with time and that it is indeed the decline in the CGM density that drives the reduction in .

Zoom In Zoom Out Reset image size

This general time-dependent picture also applies to our dwarf halos, but the results are different (see our supplementary movies for dwarfs analogous to Figure 14). In the intermediate-mass dwarfs, the equilibrium temperature remains close to the peak of the cooling curve until very late times, so there is either no phase transition or it happens very late. In these classical dwarfs, the heating terms are not able to overcome cooling so we expect roughly ∼10⁵ K temperatures. In contrast, low-mass dwarfs become thermal-pressure-dominated at very early times in our model with their CGM temperature exceeding the halo virial temperature by a factor of ∼2. These ultrafaint-scale halos have low enough densities that their CGM is susceptible to photoionization heating from the UVB, which leads to a reduced net cooling rate that cannot balance heating from SN winds, turbulence dissipation, and cosmic accretion except at a slightly supervirial temperature. In addition, we found in FIRE-2 and thus prescribed in our model that rather than 0.5 for these M_vir ∼ 10¹⁰ M_⊙ halos (see Figure 5), and this contributes to them forming thermally supported halos quite early.

We have presented a new model that is capable of predicting the time evolution of the global thermodynamic state of the CGM. This was accomplished by self-consistently linking both thermal and turbulent kinetic energy flows in the CGM to cosmic accretion and SN-driven galactic winds. In Carr et al. (2023), we used the purely thermal limit of this model to understand the shape of the SMHM relation, which has been empirically constrained by Behroozi et al. (2019; see also the recent review by Wechsler & Tinker 2018). There, we found that the SMHM relation should be quite insensitive to variations in the mass loading factor of galactic winds alone because the CGM in our new framework is self-regulated. Increasing the mass loading factor while keeping the energy loading factor fixed leads to an increase in the CGM density and thus a decrease in the CGM specific energy, which in turn increases the mass accretion rate back into the ISM without achieving the intended decrease in SFR and hence SMHM ratio. In contrast, increasing the specific energy of the winds can lead to an increase in the CGM specific energy and therefore a suppression in the cooling rate and SFR. Although we do not perform a parameter space exploration and forward modeling of the SMHM relation in this paper with the additional turbulent CGM component, we expect very similar conclusions regarding the self-regulating nature of the CGM and implications for the SMHM relation.

What is unique about our model is that it is single-phase but that phase changes naturally with time (we again urge the reader to watch our supplementary movies corresponding to Figure 14 to better understand how the individual terms contribute to the evolution of model equilibria). The phase dichotomy is such that at early times, turbulence dominates the CGM energy budget because densities and thus radiative cooling rates are very high while dissipation rates are relatively low. Hence it is only natural that should remain relatively low at early times while remains elevated for longer. This has the striking implication that the CGM–galaxy connection at early times should primarily be set by the turbulent rather than thermal properties of the CGM since the turbulence is the main source of pressure support and hence self-regulation for both the CGM and star formation. Of course, strong individual bursts of star formation can also cause jumps in the CGM thermal energy due to large-scale heating but the subsequent radiative cooling should be rather efficient, at least in the high-density conditions of the early Universe (whereas the turbulence may be more long-lived). ¹³

We find that turbulence plays an important role at early times in the CGM of MW-mass halo progenitors as well as down to low redshifts for intermediate-mass (“classical”) dwarfs. However, ultrafaint dwarfs seem to have a much shorter-lived early turbulent phase and instead transition to a warm thermal-pressure-dominated CGM at quite high redshift. At first, this may appear natural because these low-mass dwarfs are assumed to have much smaller eddy sizes, but they also have proportionately smaller turbulent velocities so that the eddy turnover timescale (i.e., the ratio R_turb/v_turb) itself remains similar to that of more-massive halos. Instead, it is simply that these halos lie below the peak of the cooling curve, and their low CGM densities make them susceptible to photoionization heating from the UV background, which together causes heating to exceed the available radiative cooling at an earlier time compared to the later phase transitions seen for more-massive halos. The inefficient star formation histories and low SMHM ratios of ultrafaint progenitors thus appear to be intimately connected to the inability of their CGM to cool efficiently at early times. However, turbulence may still be important in limiting ISM accretion and excess star formation at very early times (z ≳ 3) in ultrafaint progenitors when their CGM would otherwise have very short cooling times and be susceptible to thermal instabilities.

It is curious how closely connected the CGM phase transition is to the continuous decrease of the mean cosmic density ρ_m for all halo masses we consider. We note that the CGM masses of our halos are relatively constant with time except for an initial steep rise as halos first build up. And yet the CGM number density n(r) steadily declines with time, reflecting the importance of the growth of halo radii as the Universe expands. Since , this means that the CGM cooling time gets progressively longer with time until it becomes the dominant term in limiting the ISM accretion rate and hence star formation and turbulence driving. This effect is even more pronounced in the ultrafaint progenitors for which the ratio of specific wind energy to specific CGM binding energy is preferentially larger. Thus they have reduced baryon fractions due to both preventative feedback suppressing halo gas accretion as well as their more easily overpressurized CGM ejecting a substantial fraction of previously accreted baryons. In turn, the lower overall CGM densities of dwarfs contributes to their earlier CGM phase transition. With that said, density is not the only contributor to the phase transition: our assumption for R_turb directly sets the turbulence decay timescale so variations in this uncertain parameter also play an important role.

While CGM turbulence plays a crucial role in regulating galaxy formation at early times, it is the physics of atomic cooling that governs the phase transition and CGM thermodynamics at late times. In particular, the temperature, density, and metallicity will set the ionization state of different metals and in turn dictate the maximum achievable radiative cooling rate. When this cooling rate can no longer balance the heating from SN-driven winds, turbulence dissipation, and cosmic accretion, we get a bifurcation in the solution to our ODEs in the sense that cool equilibria vanish and the system quickly evolves to the next hotter equilibrium state as imposed by the cooling curve. Even in the intermediate-mass dwarfs that retain significant CGM turbulent pressure support down to low redshift, it is only possible because they just happen to converge to a CGM (and virial) temperature that is close to the peak of the cooling curve (∼10⁵ K). This results in relatively low t_cool,eff/t_ff,eff ∼ 1 and places classical dwarfs squarely in the thermal instability regime, allowing star formation and turbulence driving to continue unimpeded.

Another way to appreciate our model is in the context of CGM virialization, though we argue “thermalization” is a more apt description since even at early times when T₀ ≪ T_vir, we still have v_turb ∼ V_vir, which satisfies the virial theorem (see Figure 14). Birnboim & Dekel (2003) argued that gas accreting into halos more massive than ∼10¹¹ M_⊙ gets shock-heated to the virial temperature of the halo, but that, in the absence of feedback, such shocks do not form in lower-mass halos or at z > 2 owing to efficient radiative cooling of the infalling gas. Subsequent studies using cosmological simulations demonstrated that the nonshocked gas accretes along dense, cold filaments and that this “cold mode accretion” directly feeds central galaxies and powers their star formation (Kereš et al. 2005, 2009; Dekel et al. 2009; Faucher-Giguère et al. 2011; van de Voort et al. 2011). This picture is complicated by the contribution of galactic winds to the hot CGM (e.g., van de Voort & Schaye 2012; Fielding et al. 2017), the angular momentum and baryon fraction of halos (Stern et al. 2020), and potentially unresolved hydrodynamical instabilities that may cause cold gas filaments to disintegrate before they can reach the central galaxy (e.g., Mandelker et al. 2020). More recently, Stern et al. (2021) extended this standard CGM virialization (thermalization) picture to take into account radial dependence and directionality. They found that the CGM of MW-mass FIRE-2 halos virializes (thermalizes) gradually from the “outside-in” over a period of several gigayears, completing by z ≈ 0.5 (this is consistent with our own measurements of from the particle data in Figure 12). This “inner CGM virialization” (thermalization) is correlated with the formation of a thin disk, the transition from bursty to steady star formation, and damping of the strength of SN-driven winds (see also Anglés-Alcázar et al. 2017b; Pandya et al. 2021; Gurvich et al. 2023; Hafen et al. 2022). Thus the conditions of the inner CGM are important in setting the properties of the central galaxy. In contrast, intermediate-mass FIRE-2 dwarfs retain a thermally unstable inner CGM down to z ∼ 0 (even though their outer CGM is thermalized) and do not show signs of a prominent stable gaseous disk (compare Figures 4 and 5 in Stern et al. 2021).

Although our model treats the entire CGM as a single zone and is therefore not designed to address the radial dependence of CGM thermalization (a radially continuous 1D model would be needed for that) and consequently predicts a more abrupt phase transition than measured in FIRE-2, we are consistent with the overall picture of CGM thermalization. For our fiducial choice of model parameters calibrated from the FIRE-2 simulations, MW-mass halos transition from an early, turbulent, thermally unstable CGM to a warm/hot thermally supported CGM with temperature only slightly above T_vir by z ≈ 0.5 as in FIRE-2 (although we note that these results would change for different parameter choices). The CGM of intermediate-mass dwarfs does not experience this phase transition until very late times, if at all. However, we predict that the CGM of ultrafaint dwarfs should experience the phase transition at much higher redshift, contrary to the simple statement that any halo with M_vir ≲ 10¹¹ M_⊙ should always experience cold mode accretion. The low baryon fractions, low metallicities, and reduced dense ISM gas fractions of these dwarfs invites additional complications such as significant CGM photoionization heating and easier SN wind breakout, which our model accounts for phenomenologically with our parameterized loading factors and UV background-dependent cooling function.

In the future, we will use our new framework to explore whether the observable properties of galaxies themselves are correlated with the observable properties of their CGM. This will require forward modeling predictions for large populations and assessing the scatter in galaxy properties at fixed CGM properties and vice versa. We expect the CGM–galaxy connection to manifest in several observable scaling relations such as the mass–metallicity relations for both stars and ISM gas (and possibly also the CGM, although that may be harder to untangle observationally, but see Faerman et al. 2022), the stellar-to-halo-mass relation and ISM gas fractions (see our initial exploration of these in Carr et al. 2023), and the star-forming main sequence and the cause of its scatter (e.g., Rodríguez-Puebla et al. 2016). It will be very informative to see how sensitive our predictions are for these and other galaxy–CGM scaling relations as we vary our free parameters, especially for model realizations far from our fiducial FIRE-calibrated one. We also plan to investigate whether adding CGM-related observational constraints to the usual set of galaxy population metrics at z = 0 breaks parameter degeneracies and shrinks the allowed parameter space of our model.

As reviewed by Carr et al. (2023), our model builds on but also goes significantly beyond previous “bathtub” and semianalytic modeling approaches. Bathtub models generally only deal with the ODEs for the evolution of ISM mass, stellar mass, and the metallicities of these two components, and a special subset of bathtub approaches further enforce an equilibrium condition on the ISM mass (i.e., set assuming inflows balance star formation and outflows). These models generally neglect the physics of the CGM and instead have free parameters that roughly control how much cosmic accretion ends up as cold gas in the galaxy and becomes available for star formation. While this approach has achieved success in reproducing some key observed scaling relations of galaxies at a range of redshifts (see the recent review by Tacconi et al. 2020) and even summarizing the gas flow cycle in hydrodynamical simulations (e.g., Neistein et al. 2012; Mitchell & Schaye 2022), it leaves much to be desired in terms of clearly elucidating the detailed underlying physics. In contrast, we are explicitly accounting for multiple physical processes that regulate the CGM not only via an ODE for its mass evolution but also its thermal and turbulent kinetic energy evolution. The closest model to ours is the purely thermal limit version that we have put forward in Carr et al. (2023), which we would converge to as and (see also Cousin et al. 2015, which similarly tracks the thermal energy of the hot CGM phase). One difference is in how we set our model parameters: Carr et al. (2023) used a fixed power law for the mass loading factor, adopted an observationally inferred relation for the ISM depletion time, and fit for the parameters of a power law for the energy loading factor with halo mass using the Behroozi et al. (2019) SMHM relation as a constraint. In contrast, we constrain these and other free parameters directly from the FIRE-2 simulations as advocated in Pandya et al. (2020, 2021). In the future we plan to do an exhaustive parameter space exploration and comparison to observations in the spirit of Carr et al. (2023) but that will require sophisticated inference machinery that is beyond the scope of this paper.

SAMs include many more physical processes than both bathtub models and our new approach such as satellite orbital dynamics and gas stripping, merger-induced starbursts, recycling of halo outflows back into the CGM, growth of and feedback from supermassive black holes, multiphase ISM partitioning, and galaxy structural evolution (see Somerville & Davé 2015, for a recent review). It is possible that including some of these additional processes could affect the thermodynamic evolution of the CGM within our new framework. For example, the extra energy injected into the CGM from winds and turbulence stirred by satellites and black holes could affect the nature of our CGM phase transition, especially in more-massive halos. And while the standard SAM treatment of many of these extra physical processes may itself be subject to uncertainties, we expect our new CGM model to integrate naturally within the foundation of existing SAMs. Furthermore, Pandya et al. (2020) showed that existing CGM prescriptions in some SAMs can predict dramatically different CGM properties compared to simulations like FIRE-2 (e.g., orders-of-magnitude lower CGM masses for dwarfs in their Figure 6). Even though our new model shows its own discrepancies relative to FIRE-2, it is an improvement over previous approaches, and it follows the overall trends of the simulations. We therefore intend for the new model presented in this work to become the backbone of a next-generation SAM that we ourselves are building to which many of the aforementioned uncertain physical processes will be added in piece by piece.

We point out that there have been many previous efforts to revamp the way that SAMs traditionally model the CGM. As nicely reviewed by Lu et al. (2011), most existing implementations of the CGM in SAMs trace back to White & Frenk (1991) who assumed that the thermodynamics of the CGM traces that of the underlying dark matter in the sense that the CGM temperature everywhere must be equal to the halo virial temperature. A singular isothermal density profile (n ∝ r⁻²) for the CGM is typically assumed, and the metallicity is also assumed to be the same everywhere, usually with solar abundance ratios. For simplicity and historical reasons, the Sutherland & Dopita (1993) cooling function is commonly assumed to compute the radiative cooling rate, but this assumes CIE even though photoionization heating by the UV background can be important in halo outskirts and especially around dwarfs (see also Benson et al. 2002). The turbulent kinetic component is altogether generally neglected. These simplifying assumptions allow one to compute a so-called “cooling radius” within which the cooling time of the gas is shorter than some long timescale such as the Hubble time. That cooled gas is then assumed to freefall into the ISM on the halo dynamical timescale. Whenever R_cool > R_vir, it is assumed that the CGM is undergoing “cold mode” accretion and that filaments are directly free-falling into the ISM from large scales (e.g., Guo et al. 2011). There are further extensions built on top of this approach such as tracking CGM angular momentum evolution (e.g., Stevens et al. 2017; Hou et al. 2018; Lagos et al. 2018) and allowing for the simultaneous cosmic accretion of cold and hot gas (Benson & Bower 2011; Lu et al. 2011; Cousin et al. 2015).

What both Carr et al. (2023) and we have done is to show how the thermodynamic evolution of the CGM can be decoupled from that of the dark matter. By introducing the and ODEs as well as free parameters for the density, temperature, and turbulent velocity structure of the CGM, we can track energy flows in the CGM and predict its global thermodynamic state (see also Cousin et al. 2015, for another approach in the purely thermal limit). We roughly converge to the assumption that the temperature of the CGM is ≈T_vir at late times in MW-mass halos, but at early times when cooling rates are very high, we generally predict subvirial temperatures that are then self-consistently fed into the Wiersma et al. (2009) cooling function, which accounts for CGM photoionization. During these early times when CGM cooling times are extremely short, we also suggest that turbulent pressure can significantly limit ISM accretion and hence prevent excess early star formation. SN-driven winds are the dominant source of CGM turbulence except at very early times before SF kicks in, when cosmic accretion would be the sole driver (see bottom-right panel of the animation for Figure 14). This may have implications for the normalization and slope of the faint end of the stellar mass function and mass–metallicity relation for both stars and the ISM. In addition, whereas SAMs generally assume a phenomenological function to predict a mass outflow rate from the halo depending on the virial velocity of the halo, our formalism provides a prediction of the excess energy and mass that must be vented by an overpressurized CGM to remain in a quasi-HSE. As suggested by the models of Lu et al. (2015), Pandya et al. (2020), and Carr et al. (2023), this outflowing energy from the halo may act in concert with other mechanisms to pre-heat the gas outside of low-mass halos and prevent its accretion, further contributing to the reduced baryon fractions and SMHM ratios of dwarfs.

Finally, we briefly remark on the existing family of 1D CGM models that are very compelling in their ability to describe the properties of the CGM in both simulations and observations. These models envision three different physical scenarios for the CGM: steady-state cooling flows (see Fabian 1994; Stern et al. 2019, and references therein), HSE (e.g., Faerman et al. 2017, 2020; Qu & Bregman 2018), and precipitation (e.g., McCourt et al. 2012; Sharma et al. 2012; Voit et al. 2015). The advantage of these models is that they can predict CGM observables starting with very clear explanations for the underlying physical principles. Their disadvantage is that these models do not simultaneously model the galaxy formation process and are generally only applied at a single instant (though their parameters can be varied to describe CGM conditions in different mass halos or at a range of redshifts). Our approach is complementary in that we explicitly take into account time dependence via our system of coupled ODEs that self-consistently links energy flows in the CGM to cosmic accretion and SN-driven winds. However, one limitation of our approach is that we currently prescribe rather than predict the radial structure of the CGM: the slopes of the density and temperature profiles in our model are free parameters, whereas HSE, precipitation, and cooling flow models make clear predictions for these radial gradients. Another limitation is that our model is single-phase and does not account for scatter in the thermodynamic properties of the CGM as a function of radius whereas precipitation models at least attempt to capture the multiphase aspect of halo gas (see also Esmerian et al. 2021). Extending our model to include a multiphase CGM is nontrivial but will be the subject of future work.

In the future it would be interesting to forward model observables of the CGM such as column densities of various ions as a function of impact parameter, X-ray luminosities, and the Sunyaev–Zel’dovich effect. We could do this with our fiducial assumed density and temperature profiles but we could also feed the scatter in CGM properties predicted by our model at fixed halo mass into the HSE, precipitation, and cooling flow models to see what those frameworks would predict for the CGM structure and related observables. For an initial exploration of this approach, we refer the reader to Faerman et al. (2022), who took the CGM masses and metallicities for a large population of MW-mass halos from the Santa Cruz SAM (Somerville et al. 2015), generated several CGM observables by exploring the parameter space of the Faerman et al. (2020) HSE model, and placed constraints on those HSE parameters by comparing to observations of the MW CGM (see also Qu & Bregman 2018).

Although we believe our new approach to be significant step toward a more self-consistent and predictive SAM of galaxy formation, our model is still subject to many limitations and uncertainties. Here we briefly provide a nonexhaustive list of possible issues and hence avenues for future work (grouped into a few representative categories):

1. 1.  
   Turbulence: the main uncertainty is that we do not know how R_turb should vary with halo mass, redshift, and CGM/galaxy conditions. In addition, our calculation of the turbulence dissipation rate is based on the largest eddy turnover time (Equation (28)), which is appropriate for subsonic turbulence but not necessarily for supersonic turbulence in which shocks may allow the turbulence to dissipate even faster (but see Mac Low & Ferrara 1999). Relatedly, we do not distinguish between bulk flows and turbulence because we argue that any differences are effectively averaged over in our smooth model and because our model does predict bulk outflows of kinetic energy from the halo when v_turb > V_vir. However, it would be good to check this with a stochastic model for bulk flows and explicitly introduce free parameters to separate the two processes. Finally, it would be insightful to develop predictive models for and since these parameters largely control the amount of turbulence driving. We note that previous studies like Birnboim & Dekel (2003) predicted based on accretion shock calculations, but they neglected feedback, which can dramatically alter outer CGM properties especially around dwarfs.
2. 2.  
   Atomic cooling physics: we are assuming the cooling function of Wiersma et al. (2009), which is a significant improvement over the Sutherland & Dopita (1993) tables generally adopted by SAMs, because it allows us to take into account CGM photoionization heating. However, it would be good to check how our radiative cooling rates would be affected if we use more recent cooling tables (e.g., Ploeckinger & Schaye 2020, which uses the more recent UV background model by Faucher-Giguère 2020). In the same vein, there may be updated parameterizations for the fraction of gas around dwarfs that is photoionized by the UV background and hence prevented from accreting (compared to our assumed values from Okamoto et al. 2008, which uses the old Haardt & Madau 2001 UV background model). Finally, nonequilibrium ionization processes are not captured by our assumed cooling function but may have significant effects on our predictions (Tumlinson et al. 2017).
3. 3.  
   Multiphase gas: our CGM model is single-phase even though we expect gas to exist at a range of temperatures. It may be interesting to separately track the already-cooled gas and explore various models for its evolution (e.g., cold cloud scenarios; Maller & Bullock 2004; Faerman & Werk 2023). We also do not currently allow for the possibility that some fraction of cosmic accretion enters the CGM as filaments and directly deposits cold gas into the ISM without being subject to our normal energy flow cycle, which may be especially important at high redshift (e.g., Mandelker et al. 2020). The inclusion of multiphase gas may also affect the nature of the phase transition including if, when, and how abruptly it happens.
4. 4.  
   Radial structure: we only track the global thermodynamic state of the CGM, but it is likely that the inner and outer CGM have quite different properties. Indeed Stern et al. (2021) found that the thermalization of the CGM in FIRE-2 has a radial dependence and directionality. Including even a steady-state (non-time-varying) radially resolved 1D CGM model (e.g., Stern et al. 2019) within our framework could help capture a more gradual CGM thermalization process and also account for changes in the gravitational potential energy as material moves in the CGM. A radially resolved model would also negate the need for introducing a single arbitrary radius at which to define the effective freefall time and instead allow us to compute as a radial integral instead of Equation (26).
5. 5.  
   Chemical evolution: our instantaneous recycling assumption is successful in roughly reproducing the CGM, ISM, and stellar metallicities of the FIRE-2 halos as a function of time. We assumed that winds have the same metallicity as the ISM, but this is likely not the case in reality (or even in detailed simulations like FIRE). We also parameterized the halo inflow metallicity from FIRE-2, but this could be modeled more self-consistently with an outer halo wind recycling model. More generally, we should switch to time-dependent, multiyield functions to track the production and inflows/outflows of individual elements, which would allow us to self-consistently predict abundance ratios of the gas and stars and assess the impact of assuming solar abundance ratios on CGM radiative cooling rates.
6. 6.  
   Stochastic effects: our star formation model is continuous since we simply define a single ISM depletion time parameter, but star formation is a stochastic process, especially in dwarfs, in FIRE-2, and other high-resolution simulations (e.g., Muratov et al. 2015; Christensen et al. 2016; Sparre et al. 2017; Faucher-Giguère 2018; Iyer et al. 2020; Pandya et al. 2020; Gurvich et al. 2023). Thus, we should account for stochastic bursts of feedback and hence heating effects on the CGM. In particular, winds and halo outflows should occur on a variety of timescales. Related to this is the effect of recycling: some fraction of winds may recycle in the inner halo on a rapid timescale as fountain flows and some fraction of halo outflows may also recycle back into the CGM (Anglés-Alcázar et al. 2017a).
7. 7.  
   Other missing physics: connection to galaxy structure and specifically disk formation (Forbes et al. 2019), multiphase gas and turbulence in the ISM (Ginzburg et al. 2022; Forbes et al. 2023), a self-consistent star formation model, cosmic rays, and magnetic fields as additional sources of nonthermal pressure support beyond turbulence alone, additional energy input from supermassive black hole feedback, which we expect to be important for extending this model to group/cluster scales, wind heating, dynamical friction heating, turbulence driving, stochastic feedback from major mergers, and contribution of multiphase CGM gas by satellites.
8. 8.  
   Measurement uncertainties: in any analysis of simulations or observations, there will be uncertainties in derived quantities due to the definitions and techniques that are adopted. Our measurements of model parameters and galaxy/CGM properties in the FIRE-2 simulations rely, in some cases, on arbitrary choices for our chosen definition of where the CGM begins and ends (0.1–1.0R_vir), our neglect of satellite contributions to host CGM heating and cold gas content, our choice of spherical boundaries through which to track mass, energy, and metal fluxes into and out of galaxies/halos (0.1R_vir and R_vir), our decision to not impose any further velocity cuts except v_rad = 0 as the split between inflowing and outflowing particles, etc. Also, while our approach of adding up the thermal energy of all CGM particles should be robust, our measurement of the total CGM kinetic energy is an upper limit to the turbulent CGM energy component alone since bulk flows and rotation can be important. We also have not yet performed multisnapshot particle tracking to understand the longer-term evolution of particles such as whether inflows actually get to the galaxy and what fraction of halo outflows recycle back into the CGM versus become unbound from the halo forever (but see Hafen et al. 2020).
9. 9.  
   Statistical inference: we have not exhaustively explored our parameter space but instead have fixed most of our parameters to the values we measure in the FIRE-2 simulations. However, many of these parameters likely suffer from degeneracies and uncertainties, and probably vary across suites of simulations with different implementations of key physical processes. It would be good to measure the parameters of our model in many different simulations and assess the scatter as a measure of uncertainty on galaxy formation-related processes. Longer term, it would be informative to use sophisticated techniques such as implicit likelihood inference (Cranmer et al. 2020) to explore the highly multidimensional parameter space of our model using both simulations and observations as constraints. In particular, it would be good to know what extra value CGM observations bring for breaking parameter degeneracies in our model and SAMs more generally compared to just using galaxy scaling relations alone.