Cold Gas Subgrid Model (CGSM): a two-fluid framework for modelling unresolved cold gas in galaxy simulations

ABSTRACT

The cold (⁠|$\sim 10^{4}\, {\rm K}$|⁠) component of the circumgalactic medium (CGM) accounts for a significant fraction of all galactic baryons. However, using current galaxy-scale simulations to determine the origin and evolution of cold CGM gas poses a significant challenge, since it is computationally infeasible to directly simulate a galactic halo alongside the sub-pc scales that are crucial for understanding the interactions between cold CGM gas and the surrounding ‘hot’ medium. In this work, we introduce a new approach: the Cold Gas Subgrid Model (CGSM), which models unresolved cold gas as a second fluid in addition to the standard ‘normal’ gas fluid. The CGSM tracks the total mass density and bulk momentum of unresolved cold gas, deriving the properties of its unresolved cloudlets from the resolved gas phase. The interactions between the subgrid cold fluid and the resolved fluid are modelled by prescriptions from high-resolution simulations of ‘cloud crushing’ and thermal instability. Through a series of idealized tests, we demonstrate the CGSM’s ability to overcome the resolution limitations of traditional hydrodynamics simulations, successfully capturing the correct cold gas mass, its spatial distribution, and the time-scales for cloud destruction and growth. We discuss the implications of using this model in cosmological simulations to more accurately represent the microphysics that govern the galactic baryon cycle.

1 INTRODUCTION

Galactic evolution is driven by the exchange of baryons between the stellar disc and its surrounding halo. The interstellar medium (ISM) requires a consistent inflow of gas to form stars, which then expel chemically enriched gas back into the circumgalactic medium (CGM) at the end of their lifespans. In this way, the composition of the CGM provides a detailed record of a galaxy’s past and the physical processes that govern its evolution (Tumlinson, Peeples & Werk 2017; Faucher-Giguère & Oh 2023). In particular, cold (⁠|$\sim 10^4\, {\rm K}$|⁠) CGM gas is a crucial component fuelling ongoing star formation.

Cold gas is ubiquitous in the haloes around dwarf, |$L_*$|⁠, and massive galaxies across cosmic time (e.g. Thom et al. 2012; Werk et al. 2013; Cantalupo et al. 2014; Borisova et al. 2016; Mishra & Muzahid 2022; Qu et al. 2022). This cold gas plays an integral role in the galactic baryon cycle and is frequently observed both, in galactic outflows (e.g. Heckman et al. 2000; Nielsen et al. 2015; Veilleux et al. 2020; Burchett et al. 2021) and in accretion streams into the galaxy (Sancisi et al. 2008; Putman, Peek & Joung 2012; Richtler et al. 2018; Zhu et al. 2021; Kamphuis et al. 2022). Recent estimates indicate that cold CGM gas constitutes at least 25 per cent of the baryonic content in |$L_{\ast }$| galaxies (Werk et al. 2014), and in dwarf galaxies, it can account for as much as 90 per cent of galactic baryons (Zhang et al. 2020). Consequently, to fully understand the galactic baryon cycle and galaxy evolution, theoretical models must not only reproduce the observed properties of cold CGM gas but also provide robust predictions about the physical processes influencing cold gas formation, destruction, and its interactions with the hot phase.

Recent theoretical work and high-resolution idealized simulations have significantly advanced our understanding of the formation and evolution of cold gas across various contexts. For instance, it is now established that when radiative cooling is sufficiently efficient, cold gas can condense from thermally unstable hot gas (Field 1965; McCourt et al. 2012; Sharma et al. 2012; Girichidis et al. 2021) and subsequently precipitate onto the host galaxy (Voit & Donahue 2015). This thermal instability has been explored in diverse physical contexts including thermal conduction (Sharma, Parrish & Quataert 2010; Wagh, Sharma & McCourt 2014), magnetic fields (Ji, Oh & McCourt 2018), turbulence (Voit 2018), and cosmic rays (Butsky et al. 2020; Tsung, Oh & Bustard 2023). Additionally, high-resolution ‘cloud-crushing’ simulations have elucidated the conditions that lead to either the destruction or mass accretion of cold gas clouds (Armillotta et al. 2017; Gronke & Oh 2018, 2020; Li et al. 2020; Sparre, Pfrommer & Ehlert 2020; Kanjilal, Dutta & Sharma 2021; Farber & Gronke 2022; Fielding & Bryan 2022; Abruzzo, Fielding & Bryan 2023), as influenced by the thin, radiatively cooling boundary layer interface between the hot and cold gas phases (Fielding et al. 2020; Tan, Oh & Gronke 2021; Chen, Fielding & Bryan 2023a). Recent theoretical advances have begun to frame our understanding of what determines the size of cold gas clouds and under what circumstances these clouds form a mist of cloudlets or coalesce into larger clouds (Gronke & Oh 2020, 2023). However, without a cosmological context (e.g. Saeedzadeh et al. 2023), determining the relative impact of each of these processes on the observed structures of cold CGM gas is exceedingly difficult.

Studying cold CGM gas with galaxy-scale simulations presents significant challenges due to the vast range of relevant physical scales involved. Theoretically, the characteristic scale of cold CGM clouds is estimated to be around |$0.1-10$| pc (Gronke & Oh 2018; McCourt et al. 2018; Liang & Remming 2020). While it is difficult to observationally determine the minimum size of these clouds, there is evidence suggesting the existence of cold CGM clouds smaller than |$\sim 10-100$| pc (Hennawi et al. 2015; Stern et al. 2016; Rubin et al. 2018; Bish et al. 2019; Rudie et al. 2019; Werk et al. 2019; Zahedy et al. 2019; Chen et al. 2023b). The exact resolution required to accurately simulate cold CGM gas in a galactic context is still under debate. However, even simulations that employ innovative algorithms to maximize CGM resolution (reaching down to about 100 pc in the inner CGM) have not achieved convergence across all relevant cold gas properties (Hummels et al. 2019; Peeples et al. 2019; Suresh et al. 2019; van de Voort et al. 2019; Nelson et al. 2020; Ramesh & Nelson 2024).

If it is necessary to resolve scales significantly smaller than 100 pc, we may be decades away from the computational resources required to explicitly simulate cold CGM gas physics. For example, resolving the halo of an |$L_{\ast }$| galaxy out to |$200\, {\rm kpc}$| with parsec-scale spatial resolution would require roughly 10 million times more voxels than current state-of-the-art cosmological simulations. However, the current alternative of studying CGM gas using simulations that severely under-resolve its structure substantially restricts the conclusions we can draw regarding the origins of cold gas and its impact on galaxy evolution.

Instead of waiting for technological advances that enable direct resolution of the CGM at sub-pc scales, we will show that resolution issues can be effectively addressed using a physically motivated subgrid model for cold CGM gas. The use of subgrid models has a strong precedent in galaxy simulations, particularly for processes that cannot be directly resolved, such as star formation and supernova feedback (e.g. Cen & Ostriker 1992; Abadi et al. 2003; Oppenheimer & Davé 2006; Stinson et al. 2006; Hopkins, Quataert & Murray 2012; Keller et al. 2014), ISM phase structure (e.g. Springel & Hernquist 2003), and metal diffusion (e.g. Shen, Wadsley & Stinson 2010; Escala et al. 2017). More recently, several studies have developed subgrid prescriptions to model unresolved cold gas in multiphase supernova-driven winds (Huang et al. 2020, 2022; Smith et al. 2024).

In this work, we present a novel framework for explicitly modelling unresolved cold CGM gas in hydrodynamics simulations, conceptually similar to how dust is modelled in protoplanetary discs (e.g. Laibe & Price 2012a) or the more general multifluid cosmological hydrodynamics framework (Weinberger & Hernquist 2023). Our Cold Gas Subgrid Model (CGSM), treats this unresolved cold gas as a second fluid, tracking its total mass density and bulk momentum as collections of unresolved cold cloudlets. The properties of these cloudlets are determined by the properties of the resolved (‘hot’) gas phase, which we assume to be in thermal pressure equilibrium with the unresolved cold gas. The interaction terms between the standard and subgrid fluids are informed by the latest models for cold gas formation and destruction, as indicated by the idealized simulations mentioned above. As we will demonstrate, this subgrid approach effectively replicates the qualitative behaviour of cold gas in situations where traditional hydrodynamics simulations struggle due to inadequate resolution.

The remainder of this paper is organized as follows: in Section 2, we detail the physical basis of the CGSM and introduce the set of modified conservation equations governing the evolution of the two fluids. In Section 3, we use a series of idealized simulations to show that, in contrast to traditional hydrodynamics, our two-fluid model successfully captures the correct qualitative behaviour of cold CGM gas, even in cases of limit of very low resolution. The implications of the CGSM for cosmological simulations are discussed in Section 4, and we summarize our results in Section 5. Finally, in Appendix  A, we provide tests of the CGSM implementation in the enzo astrophysical simulation code (Bryan et al. 2014; Brummel-Smith et al. 2019).

2 COLD GAS SUBGRID MODEL TWO-FLUID DESCRIPTION

In this section, we detail the physical rationale and numerical implementation of our CGSM, as illustrated in Fig. 1. Briefly stated, the resolved (hot) gas phase is treated as a conventional hydrodynamics fluid, adhering to the equations outlined in Section 2.1 and modified in Section 2.3. Meanwhile, the unresolved (cold) gas phase is represented as a second fluid, cospatial with the hot fluid, following the set of modified governing equations described in Section 2.3. Given that the structure of the cold phase is not directly resolvable, we focus on evolving the total mass density and bulk momentum of this phase. We make informed assumptions about the physical characteristics of the unresolved cold gas based on the properties of its surrounding medium. This approach is qualitatively similar to modelling dust in protoplanetary discs as described in Laibe & Price (2012a), as well as the compressible multifluid hydrodynamics described in Weinberger & Hernquist (2023).

2.1 Traditional ‘resolved’ fluid equations

The evolution of the resolved gas fluid is described by the conventional hydrodynamics conservation equations for mass, momentum, and energy:

Here, |$\Phi$| represents the gravitational potential, while |$\rho _g$| and |${\boldsymbol v}_g$| denote the gas density and velocity, respectively. I is the identity matrix and the internal gas pressure, |$P_g$|⁠, is related to the internal thermal energy density |$\varepsilon _g = (\gamma _g-1)P_g$|⁠, where |$\gamma _g = 5/3$|⁠. |$S_g$| symbolizes thermal energy sources and sinks, such as energy inputs from supernova feedback events or energy losses due to radiative cooling. Throughout, the subscript ‘g’ indicates the properties of the resolved gas. The standard equations above are provided for context, leading to the modified conservation equations introduced in equations (10)–(17).

2.2 Unresolved cold gas fluid

In addition to evolving a ‘traditional’ resolved fluid, the CGSM explicitly evolves a second fluid of cold gas, comprised of numerous unresolved cloudlets, under the following assumptions.

First, we assume that the size of the unresolved cold gas cloudlets is significantly smaller than the resolution element size, |$r_{\rm cl} \ll \Delta x_{\rm cell}$|⁠, an assumption easily met in the simulated CGM with cell widths |$\Delta x \gtrsim 1$| kpc. With this disparity in scale, we make the simplifying assumption that the cold gas clouds are in thermal pressure equilibrium with the resolved hot phase and ignore any pressure gradients within the cold cloud interiors. This allows us to treat the unresolved gas as a pressure-free fluid. Given that the majority of thermal pressure originates from the hot gas phase, we can deduce the system’s thermal pressure from the properties of the resolved fluid.

Second, we assume that the unresolved cold fluid is at a constant temperature |$T_{\rm cl}$|⁠. Unless otherwise noted, we default to |$T_{\rm cl} = 10^4\, {\rm K}$|⁠. Assuming the two fluids are in thermal pressure equilibrium, the physical density of the unresolved cloudlets is given by:

where |$T_g$| is the temperature of the resolved gas phase, |$\mu$| is the mean molecular weight, |$m_p$| is the proton mass, and |$k_\mathrm{ B}$| is the Boltzmann constant. The subscript ‘|${\rm cl}$|’ is used throughout to refer to the properties of the cold gas fluid.

Third, we assume that the characteristic size of cold cloudlets is |$r_{\rm cl} \approx \mathrm{min}(c_s t_{\rm cool})$|⁠, motivated by the mist model of cold CGM gas (McCourt et al. 2018). The gas sound speed and cooling time are described by

and

for some radiative gas cooling rate |$\mathcal {L}$|⁠. For a wide range of CGM pressures, the minimum product of the sound speed and cooling time occurs at a temperature of |$T_{\rm cool, peak} \approx 10^{4.3}\, {\rm K}$| (Liang & Remming 2020). Therefore, we solve for the cold cloud radius at |$T = T_{\rm cool, peak}$|⁠,

It is important to note that both the gas sound speed and cooling time in the above equation are calculated at |$T = T_{\rm cool, peak}$|⁠. Assuming thermal pressure equilibrium, we solve for the gas density at the peak-cooling temperature, |$\rho _{\rm cool, peak}/\rho _g = T_g/T_{\rm cool, peak}$|⁠. In the simulations below, we assume |$T_{\rm cool, peak} = 10^{4.3}\, {\rm K}$|⁠.

Combining equations (4) and (7), we estimate the mass of a single cold gas cloudlet as

It then follows that the number of cold cloudlets in a given cell is given by,

where |$M_{\rm cl}$| is the total mass of unresolved cold gas in a cell.

2.3 Conservation equations for the CGSM

We explicitly evolve the total mass density (⁠|$\bar{\rho }_{\rm cl}$|⁠) and bulk velocity (⁠|$\boldsymbol {\bar{v}}_{\rm cl}$|⁠) of the unresolved cold fluid in each cell. Since the temperature of the cold gas cloudlets is assumed to be constant, we do not explicitly evolve the internal energy of the cold gas fluid. In the following sections, we introduce the new conservation equations for the unresolved fluid and outline the necessary modifications for the resolved fluid.

2.3.1 Conservation of mass

Equations (10) and (11) function similarly to the traditional mass conservation equation (equation 1), but they include an additional source term, |$\dot{\rho }_{\rm mix}$| to account for the mass transfer between the two fluids. This transfer occurs during the formation of cold clouds, their mass accretion, and their mass loss due to hydrodynamic instabilities and conduction. We provide specific prescriptions for these processes in Sections 2.4 and 2.5. In the equations above, the variable |$\bar{\rho }_{\rm cl} = {M}_{\rm cl, cell}/{V}_{\rm cell}$| represents the total mass of cold clouds in a single cell divided by the cell’s volume. It is important to note that this does not represent the physical density of the unresolved cold clouds, |$\rho _{\rm cl}$| (equation 4). Similarly, |$\boldsymbol {\bar{v}}_{\rm cl}$| denotes the aggregate velocity of the unresolved cold fluid in a given cell rather than the physical velocity of any individual cloudlet.

2.3.2 Conservation of momentum

As with the mass conservation equations, the revised momentum conservation equations include new source terms to account for momentum transfer between the cold and hot fluids. These transfers occur due to the formation, growth, and destruction of cold clouds (⁠|$\boldsymbol {\dot{p}}_{\rm mix}$|⁠), as well as from the drag force resulting from the relative motion of the two fluids (⁠|$\boldsymbol {\dot{p}}_{\rm drag}$|⁠).

The momentum transfer due to mixing of the two fluids is dependent on the mass transfer, |$\dot{\rho }_{\rm mix}$|⁠, and is expressed as follows:

We note that |$\dot{\boldsymbol p}_{\rm mix}$| only accounts for the net transfer of mass, either to the cold fluid of from the cold fluid. In reality, in the ensemble of unresolved cold clouds within a single cell, some clouds would be losing mass while others would be accreting mass. This could lead an underestimate of the true momentum transfer, for example, in the limit where |$\dot{\rho }_{\rm mix} = 0$|⁠. We leave exploring these effects to future work.

We model the drag force as a function of the relative velocity between the two fluids, |$\boldsymbol {v}_{\rm rel} = \boldsymbol {v}_g- \boldsymbol {\bar{v}}_{\rm cl}$|⁠,

where |$K_{\rm drag}$| is the drag coefficient given by:

as described in equation (41) of Laibe & Price (2012b).

2.3.3 Conservation of energy

Since the cold gas fluid has a fixed temperature, we do not explicitly track its energy field. However, the energy conservation equation for the traditional fluid is modified to include two additional source terms, |$\dot{\varepsilon }_{\rm mix}$| and |$\dot{\varepsilon }_{\rm drag}$|⁠. The energy lost from the resolved fluid due to mixing with the cold phase is calculated as

In addition, we consider energy loss due to frictional heating caused by drag between the two fluids, given by

We provide tests of the implementation of the modified conservation equations in the enzo simulation code in Appendix  A.

2.4 Cold cloud formation due to thermal instability

A primary mechanism for the spontaneous formation of cold gas from hot gas is thermal instability (Field 1965). While the cooling times of hot (⁠|$T\sim 10^6\, {\rm K}$|⁠) CGM gas are generally long, fluctuations in gas density and temperature can trigger a runaway cooling effect. This process results in cold gas condensing from the hot phase (e.g. McCourt et al. 2012; Sharma et al. 2012; Voit & Donahue 2015). In the CGSM, we model this phenomenon by first identifying cells likely to experience thermal instability. Then, we convert the energy that would have been lost due to radiative cooling in the resolved fluid to generate mass in the unresolved cold fluid.

At each time-step, we evaluate the cells of the resolved fluid for thermal instability using the following criteria:.

• The gas temperature falls within the thermally unstable range, |$10^4\, {\rm K} \le T \le 10^6\, {\rm K}$|⁠.

• The size of the cell exceeds the size necessary to resolve a single cold cloud, |$\Delta x_{\rm cell} \gt r_{\rm cl}$|

• The net energy change due to radiative cooling in the cell is negative, meaning it would cool in the absence of the two-fluid model.

If the hot gas cell is not thermally unstable, it is allowed to undergo radiative cooling (or heating) as usual. However if the gas in a cell is determined to be thermally unstable, we transfer mass from the hot gas to the cold subgrid fluid according to the following equation:

where |$\Delta t$| represents the simulation time-step. Assuming conservation of energy, the quantity of mass transferred depends on the amount of energy the hot gas cell is expected to lose due to radiative cooling, |$\Delta E_{\rm cool}$|⁠:

In the equations above, |$e_{\rm hot}\,\mathrm{ and}\, e_{\rm cold}$| are the specific energies (energy per unit mass) of the hot and cold gas. Note that |$\Delta \rho _{\rm mix, TI}$| is always positive because, by definition, |$e_g \gt e_{\rm cl}$|⁠. Therefore, in the case of thermal instability, mass is only ever transferred from the hot fluid to the cold.

We update each cell, i, at time t, as follows:

2.5 Cloud crushing

A cold cloudlet moving relative to a hot background medium may either lose or accrete mass, depending on the physical conditions of the two gas phases and their relative velocities (e.g. Klein, McKee & Colella 1994; Marinacci et al. 2010; Armillotta, Fraternali & Marinacci 2016; Gronke & Oh 2018; Li et al. 2020; Sparre et al. 2020; Banda-Barragán et al. 2021; Kanjilal et al. 2021; Abruzzo, Bryan & Fielding 2022).

In our model, we incorporate the mass transfer prescription as derived in equation (36) of Fielding & Bryan 2022, which effectively captures various regimes of cold gas growth and destruction observed in high-resolution simulations,

In the equation above, |$\chi = \rho _{\rm cl}/\rho _g$| represents the density contrast between the cold cloudlets and the hot phase, and |$\xi = r_{\rm cl}/({\boldsymbol v}_{\rm turb} t_{\rm cool, peak})$| determines whether the cold gas is accreting or losing mass, with |$\alpha = 1/4$| if |$\xi \ge 1$| and |$\alpha = 1/2$| otherwise. We assume that the turbulent velocity, |${\boldsymbol v}_{\rm turb} = 0.1 {\boldsymbol v}_{\rm rel}$|⁠, to be a constant fraction of the relative velocity between the two fluids.

In the limit of no cooling, equation (24) implies a cloud-crushing time,

2.6 Time stepping

Following the approach in Laibe & Price (2012a), we incorporate a time-step criterion that is a function of the drag coefficient, |$K_{\rm drag}$|⁠. The time-step requirement is given by

Next, we consider the time constraints imposed by the mass and momentum transfer due to thermal instability and cloud-crushing processes in the CGSM

The final time-step constraint for the CGSM is the minimum of the above constraints,

In the simulations below, we choose the constant pre-factor |$\epsilon _c = 0.2$|⁠.

3 PHYSICAL APPLICATIONS OF THE CGSM

Next, we use a series of idealized simulations to underscore the fundamental qualitative differences between traditional hydrodynamics and CGSM simulations. Specifically, we focus on examples of the formation, spatial distribution, destruction, and accretion of cold gas. As we will detail below, a recurring theme is that in the low-resolution limit typical of the CGM in galaxy simulations, traditional hydrodynamics significantly diverges from the ‘true’ behaviour of cold gas. In such instances, the CGSM is able to reproduce the essential qualitative behaviour of cold gas more accurately.

3.1 Cold gas formation in a single cell

For the initial application of our subgrid model, we choose a simple but important scenario: the formation of cold gas in a single cell. The initial conditions consist of a uniform gas with a temperature |$T_{g, 0} = 10^6\, {\rm K}$| and density |$\rho _{g, 0} = 10^{-27}{\rm g\, cm}^{-3}$|⁠. The gas undergoes radiative cooling following an analytic approximation of the cooling curve for a gas with a metallicity of |$0.3 Z_{\odot }$|⁠. While this metallicity is within the expected range for hot CGM gas (Bregman et al. 2018), we stress that the exact value has no bearing on the qualitative behaviour of the simulations discussed below. The gas is initially motionless, and experiences no external forces. Therefore, when using the CGSM, we focus solely on cloud formation due to thermal instability, omitting terms for cold cloud growth, destruction, and drag.

Fig. 2 compares the growth of cold gas mass in both the traditional hydrodynamics (orange) and CGSM (blue) simulations. In the traditional hydrodynamic approach, the gas cools uniformly until it the entire cell reaches the temperature threshold (⁠|$T \lt 5.05\times 10^4\, {\rm K}$|⁠) to be considered as ‘cold’. As a result, the evolution of cold gas mass follows a step function, centred near the cooling time determined by the initial gas properties. In contrast, the CGSM allows for a gradual increase in cold gas mass as it is progressively transferred from the resolved hot phase to the cold fluid. After one cooling cycle, approximately two-thirds of the total gas mass in the CGSM simulation is cold, compared to the traditional model where cold gas constitutes the entirety of the gas mass. After three cooling times, the proportion of cold gas in the CGSM simulation saturates at about 83 per cent of the total gas mass.

To provide context, the dashed line shows the cold mass growth rate from an idealized simulation (described in the following section) of thermal instability in which the |$\sim 20$| pc cold gas cloudlets are resolved. While there are some quantitative differences between the single-cell and fully resolved simulations, the latter underscores two key aspects of reality: (1) the growth of cold mass is a gradual process, and (2) the total cold gas mass eventually plateaus, falling short of encompassing 100 per cent of the total gas mass. Unlike the unresolved hydrodynamic simulation, the CGSM simulation successfully replicates this more realistic behaviour.

3.2 Spatial distribution of cold gas

Next, we consider the spatial distribution of cold gas in high-resolution hydrodynamics simulations, low-resolution hydrodynamics simulations, and low-resolution simulations using the CGSM. We simulate thermal instability in idealized 2D patches of CGM-like gas. The simulation setup is described in detail in Butsky et al. (2020); however, we summarize the key aspects here. We initialize a |$16\times 16$| kpc box with a static, uniform gas, setting the initial density at |$\rho _{g, 0} = 10^{-27}\, {\rm g\, cm}^{-3}$| and the initial temperature at |$T_{g, 0} = 10^6\, {\rm K}$|⁠. To seed the thermal instability, we introduce an isobaric perturbation with a white noise spectrum and 2 per cent amplitude into the uniform medium. In the CGSM, the subgrid cold fluid mass is initialized to zero. We ignore the effects of gravity, focusing solely on the in-situ formation of cold gas.

We model radiative gas cooling with a truncated power law, |$\Lambda (T) = \Lambda _0 T^{-2/3}$|⁠. This formulation ensures that the temperature of the cold gas remains at |$T_{\rm cl} = 5\times 10^4\, {\rm K}$|⁠.¹ We choose the constant, |$\Lambda _0 = 1.1\times 10^{-19}$| erg cm|$^3$| s|$^{-1}$| K|$^{2/3}$|⁠, so that the characteristic radius of cold gas clouds is 10 pc. A crucial aspect of simulating thermal instability involves incorporating a heating mechanism to maintain global equilibrium. We model this heating with a mass-weighted redistribution of the total energy lost to cooling. We refer the interested reader to Butsky et al. (2020) for additional details.

Fig. 3 compares the spatial distribution of cold gas at |$t = 2.5\, t_{\rm cool}$| across four different simulations. Progressing from left to right, the first panel shows the fractional cold gas mass (⁠|$M_{\rm cl}/M_{\rm cl, total}$|⁠) in a traditional hydrodynamics simulation resolved with |$4096^2$| cells, corresponding to a spatial resolution of |$\Delta x_{\rm cell} = 4\, {\rm pc}$|⁠. In this high-resolution simulation, the cold gas appears as a uniformly distributed mist of cloudlets throughout the CGM patch. The second panel depicts the cold gas mass fraction in a traditional hydrodynamics simulation with a much coarser |$16^2$| cell resolution, corresponding to a spatial resolution of |$\Delta x_{\rm cell} = 1\, {\rm kpc}$|⁠, which is typical for the CGM in galaxy-scale simulations. While this lower resolution simulation generates a total cold gas mass comparable to the high-resolution simulation (with less than a 10 per cent difference in cold gas masses), there is a marked contrast in the spatial distribution of cold gas clouds. Due to the limited resolution, the minimum size of cold gas clouds is constrained, resulting in fewer but larger clouds compared to the high-resolution simulation.

The distinct contrast in the spatial distribution of cold gas mass between resolved and under-resolved traditional hydrodynamics simulations is further demonstrated in the third panel of Fig. 3. Here, the results of the high-resolution simulation are mapped onto a |$16\times 16$| grid, mirroring the resolution of the under-resolved simulation. This pixelated version of the resolved simulation underscores that, on kiloparsec scales, the distribution of cold gas mass is almost homogeneous, a characteristic that the low-resolution traditional hydrodynamics simulation fails to qualitatively replicate. Conversely, as shown in the right panel, the CGSM is capable of accurately capturing the qualitative distribution of cold gas mass, even at low resolution.

3.3 Cloud crushing and growth

Next, we turn our attention to the regime in which cold clouds are destroyed within a hot wind, emphasizing how under-resolving the CGM in simulations using traditional hydrodynamics can substantially affect the lifespans of cold clouds.

To illustrate this, we first compare simple cloud-crushing simulations in the unresolved regime with theoretical predictions of cloud-destruction times. The simulation setup is a |$64\times 16\times 16$| kpc box with kpc resolution and periodic boundary conditions. The background is a uniform, hot medium, with |$T_g = 10^6\, {\rm K},\, \rho _g = 10^{-28}\, {\rm g\, cm^{-3}}$|⁠, and |$\boldsymbol {v}_g = [100, 0, 0]\, {\rm km\, s^{-1}}$|⁠. We place a single cold cloud near the source of the wind. The cold cloud has |$T_{\rm cl} = 10^4\, {\rm K},\, \rho _{\rm cl} = 10^{-26}\, {\rm g\, cm}^{-3}$|⁠, and is initially at rest. In the context of equation (24), these initial conditions correspond to |$\chi = 100$| and |$\xi = 0$|⁠, in the limit of no radiative cooling.

In the traditional hydrodynamics simulation, we replicate a severely under-resolved scenario where the total gas mass is concentrated in a single cell, resulting in a cloud of length |$\ell _{\rm cl} = 2\, r_{\rm cl} = 1$| kpc. In the CGSM simulation, we distribute the cold gas mass is evenly throughout the simulation domain. Given the temperature of the cold phase and the properties of the hot gas, we expect cold cloudlets to have |$r_{\rm cl} = 1\, {\rm pc}$|⁠. Both simulations are evolved for two cloud-crushing times.

Fig. 4 shows the gradual reduction of cold gas mass over time in both the traditional hydrodynamics simulation (top panel) and the CGSM simulation (bottom panel). At first glance, both models exhibit similar qualitative behaviour: the cold gas is progressively destroyed by the hot wind, following the theoretical cloud-crushing time-scale (equation 25). The low-resolution hydrodynamic simulation abruptly deviates from the analytic solution when the entire cloud is above the threshold temperature (⁠|$T_{\rm thresh} = 3\times 10^4\, {\rm K}$|⁠) to be considered cold. However, the most notable difference in behaviour lies in the respective time-scales. In the low-resolution traditional hydrodynamics simulation, the destruction time is nearly a thousand times longer due to the artificially enlarged size of the cold cloud. In contrast, the CGSM simulation accurately captures the shorter cloud-destruction time-scales, even at the same resolution.

Finally, we consider the low-resolution effects on cloud-growth time-scales. To do this, we start with the same initial conditions as in Fig. 4 and turn on radiative cooling as described in Section 3.2. In these cloud-crushing simulations, we do not artificially truncate the radiative cooling curve, and allow gas to cool down to |$10^4$| K. Given the expected cooling time of the mixed-temperature gas stripped off of the cold cloud (⁠|$T_{\rm mix} = \sqrt{T_gT_{\rm cl}}$|⁠), these initial conditions correspond to |$\xi \approx 30$|⁠, placing them firmly in the expected cloud-growth regime (Fielding & Bryan 2022). Yet, Fig. 5 shows that in the low-resolution limit, the cold cloud in the traditional hydrodynamics simulation loses mass, albeit, more slowly than it did without radiative cooling. While this simulation behaviour is expected at low resolutions and small domains (Gronke & Oh 2018), it once again underscores that with |$\sim$|kpc resolution, traditional hydrodynamics simulations struggle to capture the correct qualitative behaviour of cold CGM gas. Meanwhile, the CGSM allows simulations of the same resolution to capture the subgrid physics of cold mass accretion.

4 DISCUSSION

4.1 Requirements and limitations of simulating cold clouds as a pressureless fluid

A key assumption of the CGSM approach is that the unresolved cold clouds are significantly smaller than the size of the resolution element. Therefore, this approach is primarily suited for studying cold gas in the unresolved CGM of galaxy-scale and cosmological simulations. In this context, the implications of modelling the unresolved cold clouds as a pressureless fluid are that we ignore any internal pressure forces in the cold cloud interiors. Also, by imposing pressure equilibrium, we are ignoring any physical processes that introduce pressure anisotropies. For these reasons, the CGSM model is best suited for applications in which the internal pressure anisotropies and dynamics of individual cold clouds is dynamically unimportant.

By modelling cold gas a separate fluid, the CGSM model removes the prohibitively high-resolution requirements imposed by resolving the sharp density gradients between the cold and hot gas phases. However, the CGSM still requires that the velocity gradients within each of the fluids are well resolved.

4.2 Comparison to existing subgrid models

In the last few years, there have been several new subgrid models to address the resolution challenges facing galaxy-scale simulations. Models such as the physically evolved winds (PhEW; Huang et al. 2020, 2022) and Arkenstone (Smith et al. 2024) specifically address the difficulties of accurately simulating multiphase galactic winds driven by supernova feedback. These approaches are motivated by the results of high-resolution, idealized simulations of ISM patches, which have shown that while multiphase outflows predominantly carry their energy in the hot gas phase, the bulk of their mass is in the cold gas phase (e.g. Li & Bryan 2020; Kim et al. 2020a, b). Directly resolving this energy-mass partitioning is currently beyond the capabilities for galaxy simulations in a cosmological context.

To overcome this limitation, the PhEW and Arkenstone models introduce a novel type of wind particle during supernova feedback events. These wind particles, conceptually similar to the CGSM presented here, represent unresolved cold gas mass and engage in mass, momentum, and energy exchanges with the surrounding medium before they are ultimately recoupled with the ‘regular’ particles in the CGM.

A key difference between the PhEW and Arkenstone models and the CGSM lies in the spatial extent and continuity of the cold gas subgrid prescription. In the PhEW and Arkenstone frameworks, cold wind particles are generated exclusively during supernova feedback events. Once these particles recouple with the regular gas, the capacity to model subgrid cold gas ceases. In contrast, the CGSM allows every cell in the simulation to host arbitrarily small amounts of cold gas. This approach facilitates evolving a continuous distribution of subgrid cold gas throughout the entire simulation domain. Furthermore, the grid-based approach to the CGSM is particularly well suited for solving the conservation equations in low-density regions.

4.3 Missing physics and future applications

In this work, the CGSM intentionally omits certain physical processes, such as magnetic fields, conduction, turbulence, and cosmic rays, to focus on presenting a ‘proof of concept’ of the advantages of a subgrid approach in CGM simulations. While this missing physics would likely alter the quantitative aspects of the cold–hot gas interactions – such as the rate of thermal instability growth or the time-scales for cloud mass loss and accretion – the modular design of the CGSM makes it straightforward to incorporate changes to the expressions for |$\dot{\rho },\, \dot{\boldsymbol p}$|⁠, and |$\dot{\varepsilon }$|⁠. Future modifications will incorporate new physics or reflect updated methodologies based on the latest high-resolution simulations, such as using a power law to model unresolved cold cloud sizes or accounting for the turbulent velocity dispersion of unresolved clouds within a cell.

Importantly, changes to the quantitative details of the interaction terms do not compromise the qualitative advantages offered by the CGSM over standard hydrodynamics simulations in the low-resolution limit. The CGSM’s ability to (1) capture the gradual accretion of cold gas mass within a single cell, (2) generate smooth spatial distributions of cold gas, and (3) accommodate short cloud-destruction and growth time-scales remain fundamentally robust, irrespective of these potential updates.

The assumption in the CGSM that subgrid cold gas clouds are significantly smaller than the typical CGM resolution (⁠|$\Delta x \sim 1$| kpc) may not hold scenarios where extreme non-thermal pressure leads to the formation of significantly larger cold ‘clouds’ with characteristic sizes larger than |$\sim 1-100$| kpc, as seen in some CGM simulations that include cosmic ray physics (e.g. Salem, Bryan & Corlies 2016; Butsky & Quinn 2018; Buck et al. 2020; Ji et al. 2020; Butsky et al. 2022). However, in such cases, the non-thermal pressure can be factored into the approximation of the expected cold cloud sizes, as demonstrated in Butsky et al. (2020). When cold clouds are large enough to be resolved at a given resolution, the CGSM does not need to be applied in that region, and the mass of the cold subgrid fluid would simply be zero.

In future work, we plan to incorporate the missing physics described above and calibrate the CGSM for use in cosmological zoom-in simulations. The combination of a cosmological context and physically motivated treatment of subgrid cold-gas physics will enable us to better determine the origin and impact of cold gas in a variety of contexts, including the CGM, galactic winds and mass accretion, as well as high-velocity clouds in our own Milky Way.

5 SUMMARY

In response to the challenge of resolving cold CGM gas in galaxy simulations, we introduce a two-fluid framework for modelling the subgrid physics of unresolved cold gas. The CGSM is designed to explicitly evolve the total mass density and bulk momentum of unresolved cold gas cloudlets. It uses the properties of the resolved gas fluid to inform predictions about the physical state of cold gas. In this model, the unresolved cold fluid interacts with the resolved hot fluid, exchanging mass, momentum, and energy in accordance with the findings from high-resolution, idealized simulations (Fig. 1).

The CGSM offers several distinct benefits over traditional hydrodynamics methods in situations where the resolution is significantly lower than necessary to adequately resolve cold-gas structure. In contrast to traditional hydrodynamics simulations, which are limited to a single-phase, single-temperature gas within each cell, the CGSM allows for the presence of arbitrarily small amounts of cold gas throughout the simulation (Fig. 2). As a result, the CGSM is capable of producing more realistic spatial distributions of cold gas mass. This is in stark contrast to under-resolved traditional hydrodynamics simulations, which tend to accumulate cold gas mass in a limited number of large clouds, with sizes artificially inflated by the size of the low-resolution voxels (Fig. 3). Furthermore, where under-resolved hydrodynamics simulations predict artificially long cloud-destruction and accretion time-scales, the CGSM captures the expected behaviour of cold gas, even when operating at the same resolution (Figs 4 and 5).

These findings suggest that in the limit of low resolution – as is typical in the haloes of galaxy-scale simulations – traditional, single-fluid hydrodynamic simulations may be unreliable tools for determining the origin and evolution of cold CGM gas. Even if the simulations converge on certain cold-gas metrics, such as the total cold gas mass, by artificially inflating cold-gas sizes and evolution time-scales, we cannot rule out that such simulations are finding the ‘right answer’ for the wrong reasons.

Certainly, opting for a subgrid model comes with its own set of trade-offs. Fundamentally, this approach introduces new simulation parameters that require precise tuning, nuanced resolution requirements, and the assumption that the subgrid model accurately represents the ensemble of cold clouds. Using a subgrid model for cold CGM gas also means that the simulations cannot be used to study small-scale cold gas structure and evolution. However, for simulations with a prohibitively large dynamic range of physical and temporal scales, the subgrid approach is unavoidable. For example, there is a strong precedent for using subgrid models of star formation and stellar feedback in galaxy-scale simulations. While such simulations cannot be used to study stellar evolution or supernova remnants, the subgrid approach has been invaluable for understanding the effects of star formation and stellar feedback on galaxy evolution. We are now faced with a similar trade-off in studying galactic haloes.

For those seeking to study the flow of cool gas and its relationship to galaxy evolution on cosmological scales, the subgrid model approach is likely inevitable. The alternative – inferring cold gas properties from simulations where the resolution elements are orders of magnitude larger than the actual cold gas structures – is fundamentally flawed. Choosing the right approach hinges on a better understanding of which scales need to be resolved in order to accurately model cold-gas physics. Should it be determined that the required scales are small (⁠|$\lesssim$| pc) compared to the typical resolution in the simulated CGM (⁠|$\gtrsim$| 100s pc), then directly resolving cold CGM gas would be computationally infeasible with current technologies. In anticipation of such a scenario and rather than waiting potentially decades for the requisite computational advancements, the CGSM offers a means to effectively model cold-gas physics within the limitations of current resolution capabilities.

SOFTWARE

The CGSM was implemented in the enzo astrophysical simulation code (Bryan et al. 2014; Brummel-Smith et al. 2019). The analysis of the simulations relied heavily on the yt (Turk et al. 2011), matplotlib (Hunter 2007), and numpy (Harris et al. 2020) packages for the python (Perez & Granger 2007) programming language.

ACKNOWLEDGEMENTS

We thank the referee, Rainer Weinberger, for constructive comments and suggestions that improved this manuscript. The authors would also like to thank Greg Bryan, Drummond Fielding, Matthew Smith, and Peng Oh for insightful conversations that contributed to the development of ideas presented in this paper. ISB was supported by HST Legacy grant AR-15800, the DuBridge Postdoctoral Fellowship at Caltech, and by NASA through the Hubble Fellowship, grant HST-HF2-51525.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS 5–26555. CBH is supported by NSF grant AAG-1911233, and NASA grants 80NSSC23K1515, HST-AR-15800, HSTAR-16633, and HST-GO-16703. Support for PFH was provided by NSF Research grants 1911233 and 20009234, NSF CAREER grant 1455342, and NASA grants 80NSSC18K0562 and HST-AR-15800.001-A.

DATA AVAILABILITY STATEMENT

The data supporting this article and CGSM source code are available on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: TESTS OF MODEL BEHAVIOUR

In this section, we test the core functionality of the two-fluid model and its implementation in the enzo astrophysical simulation code (Bryan et al. 2014; Brummel-Smith et al. 2019).

In Fig. A1, we demonstrate the advection of the two-fluid model in the presence of a shock. For this test, we set up a modified Sod shocktube with both the regular fluid and the unresolved cold fluid initialized to the same values described below. We simulate the shocktube in a 1D domain with |$x \in [-0.5, 0.5]$|⁠, resolved by 200 cells. The initial conditions are given by |$\rho _g = 1, \bar{\rho }_{\rm cl} = 1,\,\mathrm{ and}\, P_g = 1$| for |$x \le 0$|⁠, and |$\rho _g = 0.125, \bar{\rho }_{\rm cl} = 0.125, \,\mathrm{ and}\,P_g = 0.1$| for |$x \ge 0$|⁠. The velocity of both the regular fluid and the cold subgrid fluid is initialized to zero everywhere. |$P_g = (\gamma -1)\varepsilon$| is the thermal gas pressure with |$\gamma = 5/3$|⁠. We evolve the shocktube for |$t = 0.2$| internal time units.

In the case of strong drag (⁠|$K_{\rm drag} = 1000$|⁠), the analytic solution is given by the black dashed line. In this case, the simulated shocktube follows the analytic solution well. There is no analytic solution for the case of weak drag. Instead, we repeat the numerical test in Laibe & Price (2012a) with |$K_{\rm drag} = 1$| and plot the results in the bottom panel of Fig. A1. The behaviour of the two fluids agrees well with the results in Paardekooper & Mellema (2006) and Laibe & Price (2012a).

In Fig. A2, we demonstrate the performance of the drag coupling with a modified version of the dustybox test in Laibe & Price (2012a). The physical setup of this problem is similar to that of the shocktube, only the gas properties are uniform throughout the entire space and the boundary conditions are periodic. The normal gas is initially at rest and the cold fluid is initialized with a velocity in the |$\hat{x}$|-direction. If the cold and regular fluids are coupled through a drag coefficient, then the velocity of the cold gas will decrease over time as it imparts momentum on the regular gas. For this test, we consider the simplest case, in which the drag coefficient, |$K_{\rm drag}$|⁠, is a constant.

We simulate this process for a variety of different drag coefficients, |$K_{\rm drag}$|⁠, and cold gas density ratios, |$\bar{\rho }_{\rm cl}/\rho _{\rm g}$|⁠. In all cases, the initial hot gas density and the velocities of the two fluids are: |$\rho _{\rm g} = 1, \boldsymbol {v}_{g} = 0, \,\mathrm{ and}\, \boldsymbol {\bar{v}}_{\rm cl} = 1$|⁠. When testing the effect of the drag coefficient, we keep |$\bar{\rho }_{\rm cl}/\rho _{\rm g} = 1$| constant and vary |$K_{\rm drag} \in [0.01, 0.1, 1, 10, 100]$|⁠. When testing the effect of the cold gas density ratio, we fix |$K_{\rm drag} = 1$| and vary |$\bar{\rho }_{\rm cl}/\rho _{\rm g} \in [0.01, 0.1, 1, 10, 100]$|⁠.

Fig. A2 shows the time evolution of the cold gas velocities for a variety of different drag coefficients (top) and ratios of cold gas density to regular gas density (bottom). The black dotted lines show the analytic solutions, which agree with the simulation results in all cases. When the drag coefficient is low or the cold gas density is high relative to the hot gas density, the velocity of the cold gas remains relatively unchanged over time. When the drag coefficient is high or the cold gas density is low, the cold gas velocity quickly reaches its equilibrium velocity.

Author notes