© The Authors 2025

1. Properties of the circumgalactic medium

The CGM encompasses an extensive gas-rich envelope surrounding most galactic discs, spanning from dwarfs to massive ellipticals and extending to the virial radius and beyond (Tumlinson et al. 2017). Observational studies suggest that the CGM can harbour a gas mass surpassing that of a galactic disc, underscoring its significance as a dynamic reservoir. The CGM plays a pivotal role in galactic evolution, serving both as a source of gas accretion for star formation and as a sink for energy, mass, and momentum, which are injected by feedback-driven galactic winds (Faucher-Giguère & Oh 2023).

The CGM exhibits a complex multiphase structure characterised by significant variations in density, temperature, and ionisation state (Chen et al. 2010; Prochaska et al. 2011; Sameer et al. 2024). It consists of a hot diffuse phase with temperatures exceeding 10⁶ K, as observed in X-ray studies (Anderson & Bregman 2010; Gupta et al. 2012); a warm phase at around 10⁵ K, traced by UV absorption lines (Wakker et al. 2012); and a cold phase near 10⁴ K (Wisotzki et al. 2018), which may account for up to 50% of the galactic halo's baryonic gas mass (Werk et al. 2014). Cold dense gas is present in various environments such as multiphase outflows (Veilleux et al. 2020), high-velocity clouds (Bregman 1980; Wakker & van Woerden 1997; Putman et al. 2003; Maller & Bullock 2004; Richter et al. 2017), and accretion streams that supply material to galactic discs (Rubin et al. 2010; Stern et al. 2024). Despite its significance, the origin of the cold phase within the CGM remains a fundamental question in the context of galactic evolution that is yet to be resolved.

A common idea regarding the origin of this phase is that the cold gas is transported from the interstellar medium (ISM) into the CGM via galactic winds and outflows triggered by supernovae (SNe) and active galactic nuclei (AGNe). In such turbulent settings, cold clouds can be disrupted and fragmented by the ambient hot wind (McCourt et al. 2018; Gronke & Oh 2018; Sparre et al. 2019, 2020; Li et al. 2020), reshaped or protected by magnetic draping layers (Dursi & Pfrommer 2008; McCourt et al. 2015; Jung et al. 2023; Ramesh et al. 2024), or dissociated by the intense ultraviolet (UV) radiation field emanating from the galactic centre (Decataldo et al. 2019), which are processes that significantly alter the survival characteristics of cold clouds. Another theory suggests that the gas condenses locally via TI (Field 1965; McCourt et al. 2012; Sharma et al. 2012), a process initiated by local density perturbations, and eventually this condensation results in the precipitation of cold clouds onto the galactic disc (Voit et al. 2015, 2017).

Observations indicate that the cold phase of the CGM is out of pressure equilibrium with the surrounding hot medium (Werk et al. 2014), which may suggest additional non-thermal pressure sources (McQuinn & Werk 2018)¹ to maintain the overall stability of the CGM (Tumlinson et al. 2017). Cosmic rays, as a pervasive and energetically significant component of galactic ecosystems, are increasingly recognised for their role in shaping the dynamics of the CGM (Ferrière 2001; Owen et al. 2023; Ruszkowski & Pfrommer 2023). In the galactic disc, these relativistic particles can provide non-thermal pressure support comparable to or exceeding that of thermal gas, magnetic fields, or turbulence (Boulares & Cox 1990; Naab & Ostriker 2017), thus influencing large-scale dynamics and the stability of gas phases. The CRs acquire their high energies through diffusive acceleration processes at supernova-induced shocks (Marcowith et al. 2016) or from AGN-driven feedback (Guo & Oh 2008; Jacob & Pfrommer 2017a, b) and are transported from their sources within the ISM into the CGM by galactic winds, as shown in global galaxy simulations (Uhlig et al. 2012; Salem & Bryan 2014; Ruszkowski et al. 2017; Buck et al. 2020; Dashyan & Dubois 2020; Hopkins et al. 2020; Thomas et al. 2023, 2025) or in local high-resolution ISM setups (Girichidis et al. 2016; Farber et al. 2018; Sike et al. 2024).

Although previous studies have highlighted the potential of CRs to modulate TI (Butsky et al. 2020; Beckmann et al. 2022) and to influence the multiphase structure of the CGM (Butsky & Quinn 2018; Ji et al. 2020; Thomas et al. 2025), the interplay between CR transport mechanisms and the resulting cold gas morphology remains poorly understood. Cosmic rays have the ability to both suppress and enhance TI depending on the underlying transport physics and the properties of the ambient medium (Wiener et al. 2017a; Hopkins et al. 2021). Furthermore, CR heating driven by interactions with Alfvén waves may alter cooling rates, thereby affecting the conditions under which cold gas forms and survives (Zweibel 2013; Ruszkowski et al. 2017).

Addressing the challenges related to CR modulated TI in the CGM necessitates detailed simulations that capture the full complexity of CR physics. By systematically exploring the effects of varying CR pressure, transport mechanisms, and transport speeds, we aim to develop a comprehensive understanding of how CRs shape the CGM. In this study, we focus on the impact of CRs on the morphology and observable properties of cold gas formed via TI in CGM environments. To achieve this, we performed a series of three-dimensional cosmic ray magnetohydrodynamics (CRMHD) simulations using the moving-mesh code Arepo. We examined the outcomes of simulations implementing different CR transport mechanisms across a range of environments, from those dominated by thermal pressure to those heavily influenced by CR pressure. The outline of this paper is as follows. In Sect. 2, we provide a brief explanation of the theoretical background of the process of TI and the nature of CR physics. In Sect. 3, we present our numerical and physical setup along with the initial conditions for the simulated CGM-like environments. We discuss the suppression of TI employing purely advective CR transport in Sect. 4. In Sect. 5, we investigate and analyse the revival of TI in the context of two-moment CR transport. In Sect. 6, we detail the importance of CR transport speed on the onset of TI and compare the purely diffusive CR transport to the two-moment approach in this context. Section 7 is dedicated to emphasising the critical role of numerical resolution in simulating CRs in the CGM. We close this paper with a discussion of our results and the limitations of the employed setup in Sect. 8, and we present our conclusion in Sect. 9. In Appendix A, we address the issue of numerical resolution in a simulation of TI, while we study the influence of the ratio of cooling to free-fall time on the emerging multiphase structure in our TI setup in Appendix B.

2. Cosmic ray modulated thermal instability

In this section, we provide a concise summary of the key physical processes that regulate CR-modulated TI. We also outline the metrics used for our analysis throughout this study.

2.1. Thermal instability

Thermal instability is a self-amplifying process driven by the dependence of the cooling function on temperature and density in astrophysical plasmas (Field 1965). In an initially stable gas distribution, small density enhancements accelerate cooling, leading to a reduction in thermal pressure support relative to the surrounding hot environment. This pressure imbalance causes the over-dense regions to condense further, increasing the local density and amplifying the cooling process. This positive feedback loop results in the growth of density perturbations, culminating in the formation of cold, dense structures within the medium.

The cooling time-to-free-fall time ratio, τ=t_cool/t_ff, has been identified as a critical parameter for initiating this runaway process (McCourt et al. 2012; Sharma et al. 2010, 2012; Voit et al. 2015). In a stratified atmosphere, as cooling gas parcels sink deeper into the gravitational potential, the ambient thermal pressure increases, inducing compressive heating. For τ≲1, cooling dominates over compressive heating, allowing TI to grow and form dense condensations. However, when τ≳10, the onset of TI is suppressed as a result of buoyancy.

Several additional factors influence the development of TI beyond τ. For instance, the amplitude of initial density perturbations can significantly impact its growth (Choudhury et al. 2019; Esmerian et al. 2021). Magnetic fields, which modify gas dynamics, can either suppress or enhance instability depending on their configuration (Ji et al. 2018; Fournier et al. 2024). Similarly, turbulence (Voit 2018; Mohapatra et al. 2023), metallicity (Das et al. 2021), thermal conduction (Koyama & Inutsuka 2004; McCourt et al. 2012; Wagh et al. 2014), and of course, CRs (Ruszkowski & Pfrommer 2023) all play essential roles in regulating the onset and progression of TI.

Figure 1 presents 2D slices of gas mass density (left panel) and temperature (right panel) from a simulation without CR physics. Magnified insets show regions where TI is actively developing. These insets highlight the size and structure of the emerging cold gas clouds. The simulation employs τ=0.3, corresponding to a regime where radiative cooling dominates over compressive heating. This configuration serves as the fiducial case for this study. We confirm findings of previous studies, detailed in Appendix B, which explore TI in the absence of CR physics. Building on this foundation, the focus of the present work shifts to the influence of CRs on the formation and evolution of TI in CGM-like environments.

Fig. 1. Simulation domain overview. The figure shows a slice through the x–z plane at y=0 of our tall-box setup. The box extends from −3H to +3H in the z-direction and from −0.5H to +0.5H in the x- and y-directions, where H=30.11 kpc. Periodic boundary conditions are applied in the x- and y-directions, while the z-direction features an outflow boundary condition. The simulation is conducted without CRs and assumes a target mass of mtarget=22 M⊙. The computational domain is initialised with 256×256×1536 equally spaced mesh-generating points. The left elongated panel illustrates the resulting mass density, and the right elongated panel displays the temperature. The five smaller panels provide magnified views of various structures of different sizes, shapes, densities, and temperatures formed due to TI.

2.2. Cosmic ray physics

We used the CRMHD framework developed by Thomas & Pfrommer (2019). While we summarise the key physical concepts of this two-moment method for CR transport here, we refer the reader to Thomas et al. (2023) for a more detailed and illustrative description of the underlying CR physics.

The transport of CRs, being charged particles, is influenced by their interaction with magnetic fields. Their gyroradii (∼0.25 AU for giga-electronvolt CRs in microgauss magnetic fields) are much smaller than the typical length scales of the CGM dynamics we are investigating here. Consequently, we can assume that they are tied to magnetic field lines. Because magnetic field lines also follow the gas motion, we can also readily assume that CRs are advected with the gas flow. On top of this, CRs are free to propagate anisotropically along magnetic field lines but interact with small perturbations of the magnetic field which determine how fast CRs are transported along magnetic field lines. Notably, CRs can excite Alfvén waves which are plasma waves that propagate at the Alfvén speed, v_a. When the mean CR transport speed exceeds the local Alfvén speed (v_cr>v_a), CRs will excite and amplify resonant Alfvén waves through the gyroresonant instability (Kulsrud & Pearce 1969; Shalaby et al. 2023). These Alfvén waves act as magnetic perturbations that affect the CR's gyromotion along the field lines leading to effective scattering of the CRs and a transfer of energy and momentum to the waves. Conversely, when the mean CR transport speed is lower than the Alfvén speed (v_cr<v_a), the excitation of the gyroresonant instability is suppressed. Instead, CRs can gain energy from Alfvén waves through a second-order Fermi process (Ko 1992) resulting in their acceleration.

The scattering of CRs off Alfvén waves redistributes their momentum and energy, effectively limiting their transport speed to a value close to v_a. In regions with high scattering frequencies, ν_cr, CRs remain strongly coupled to the waves, resulting in highly constrained motion along the magnetic field lines. This regime is characterised by CR streaming, where CR propagation is regulated by the local magnetic field properties and the efficiency of wave generation and damping. In contrast, if CR-wave scattering is inefficient, CRs are less confined by Alfvén waves and can diffuse more freely along magnetic field lines. The transport of CRs becomes primarily diffusive, which we characterise by the CR diffusion coefficient ${\kappa }_{\text{cr}}\sim {\nu }_{\text{cr}}^{-1}$. In this regime, CRs experience reduced wave-particle interactions resulting in an increased effective CR transport speed, v_cr, which redistributes them faster throughout the CGM.

Cosmic rays transfer part of their energy directly to the surrounding gas through hadronic and Coulomb interactions (Pfrommer et al. 2017). Additionally, there is an indirect channel for energy transfer mediated by Alfvén waves. These waves interact with the plasma and undergo non-linear Landau damping (Miller 1991) gradually transferring their energy to the gas. This damping reduces the wave amplitude, lowers the CR scattering rate, and also heats the surrounding gas. In this work, we focus exclusively on non-linear Landau damping, which is considered the dominant damping mechanism in CGM-like environments.

2.3. Cold gas metrics

We utilised four key metrics to assess the results of our simulations: density fluctuations, cold mass fraction, cold mass volume, and cold mass flux. Each parameter offers valuable insights into different aspects of the underlying physical processes.

The density fluctuations,

$$\frac{\delta \rho }{\rho }=\frac{\rho -\overline{\rho }}{\overline{\rho }},$$(1)

quantify local perturbations in gas mass density, ρ, relative to the global volume-weighted average background density, $\overline{\rho }={\sum }_i({\rho }_i{V}_i)/{\sum }_i{V}_i$, where V_i is the volume of the ith computational cell. This metric helped us identify regions where the local density significantly deviates from the background, which is crucial for understanding the dynamics of structure formation and collapse.

The cold mass fraction,

$$f_m=\frac{m_{\text{cold}}}{M},$$(2)

represents the proportion of cold gas mass, m_cold, relative to the total mass, M, within the analysed domain. This metric provides insights into the efficiency of gas cooling and conversion from the hot into the cold phase by TI.

The cold mass volume,

$$f_V=\frac{V_{\text{cold}}}{V},$$(3)

describes the spatial volume occupied by cold gas, V_cold, relative to the total volume, V, offering a measure of the extent and distribution of cooler and denser regions. This parameter provides information about the size and mean density of collapsing regions.

The cold mass flux,

$$\frac{{\dot{m}}_{\text{cold}}}{{\dot{m}}_{\text{ff}}}=\frac{{\rho }_{\text{cold}}}{{\rho }_0}\frac{{\mathit{v}}_{\text{cold}}}{{\mathit{v}}_{\text{ff}}},$$(4)

characterises the rate at which cold gas is accreting towards the midplane. Here, ρ_cold is the density of the cold gas, and v_cold denotes the velocity of infalling cold material. The reference values ρ₀ and v_ff=H/t_ff are normalisation factors where ρ₀ is the initial density in the centre of the simulation domain, v_ff is a velocity defined by the scale height H, and the free-fall time t_ff.

For the purpose of our analysis, a gas parcel is classified as cold if its temperature is at most 5×10⁴ K. To mitigate contamination from potential boundary effects, we assess these metrics within a specified subrange of the computational domain, bounded by 0.25H≤|z|≤2.75H. This range is chosen to avoid unphysical behaviours near the simulation boundaries that could artificially influence the analysis. Within this region, we calculate each metric for every snapshot of the simulation, enabling a consistent comparison across time and different simulation setups. We run each simulation for 8 t_cool, using a fiducial value of τ=t_cool/t_ff=0.3. Typically, we restrict our analysis up to t=7 t_cool≈2 t_ff as the system's dynamics are expected to be significantly altered beyond this period.

2.4. Cold cloud identification

To identify individual cold clouds in our simulations during post-processing, we utilised a watershed-like segmentation algorithm that leverages the Voronoi tessellation inherent to Arepo's computational mesh. The process begins by reconstructing the Delaunay triangulation (Virtanen et al. 2020) based on the mesh-generating points of the snapshot. In this triangulation, nodes correspond to the geometric centres of the original Voronoi cells, and edges represent connections between neighbouring cells.

For each node, we determined and stored the indices of all its neighbouring nodes. The nodes that satisfy specific temperature and density criteria, T_node≤T_thresh and ρ_node≥ρ_thresh, were flagged as cloud candidates. From these candidates, spatially connected groups of nodes were identified using a graph-based clustering approach (Csardi & Nepusz 2005). Any connected group containing at least five candidates was classified as a ‘cloud.’ By construction, these clouds are disjointed and non-overlapping. We estimated the effective radius of each cloud by first calculating the total volume of all associated cells. Assuming a spherical cloud geometry, the cloud radius was calculated by $r_{\text{cloud}}={\left(3V_{\text{cloud}}/4\pi \right)}^{1/3}$, where V_cloud is the total volume of the cloud's cells. The position of each cloud was determined by its centroid, which we calculated as the mass-weighted average of the positions of all computational cells associated with the cloud. An example of the cloud-identification algorithm is illustrated in Fig. 2.

Fig. 2. Projection in x–z of the hydrogen column density, NH. Red circles represent the radii of clouds detected by the cloud-identification algorithm and red dots indicate their respective centroids. We note that overlapping circles result from the projection of clouds located at different y positions and do not necessarily represent physical overlaps in three-dimensional space.

3. Simulation setup

3.1. Numerical framework

We performed 3D CRMHD simulations with the moving-mesh code Arepo (Springel 2010; Pakmor et al. 2016) using standard parameters for mesh regularisation (Vogelsberger et al. 2012; Weinberger et al. 2020), and the magnetic divergence cleaning method described in Powell et al. (1999) using the implementation of Pakmor & Springel (2013). The equations of the two-moment CRMHD formalism as described in Thomas & Pfrommer (2019) are implemented using a finite-volume method (Thomas et al. 2021). In this approach, the CRMHD equations apply the P1 Eddington approximation for closure, which yields stable and consistent outcomes when compared to closure schemes of higher orders (Thomas & Pfrommer 2022).

The Arepo code provides the ability to refine and de-refine computational cells based on arbitrary user-defined criteria. We used the standard Lagrangian refinement scheme of Arepo to keep the mass enclosed by a computational cell nearly constant. If the mass of a cell deviates from its target value² by more than a factor of two, the cell is refined or de-refined until the criterion is satisfied again. On top of that, in order to avoid significant volume deviations and thus large density gradients between neighbouring cells, computational cells are refined if one of their Voronoi neighbours has a volume that is smaller by a factor of 5. To manage computational resources efficiently, we set the minimum volume for a computational cell to ${\left(H/2N\right)}^3$, where N represents the number of cells³ used to resolve one scale height, H, in the initial mesh.

Cosmic ray transport is computed assuming a reduced speed of light, c_red=3000 km s⁻¹. We subcycle the CR transport along the magnetic field lines ten times for each call of the MHD solver. This choice ensures that the characteristic speeds of sound and Alfvén waves, which typically peak at ∼400 km s⁻¹, are well separated from the reduced speed of light, which is thus still the fastest signal speed in the simulation. We implement CR cooling following the approach by Pfrommer et al. (2017) such that the effective CR transport velocity, f_cr /(γ_crɛ_cr), remains constant during the CR cooling process. We note that in the low-density environment of our simulations (n∼10⁻³ cm⁻³), CR cooling has a minimal impact⁴ because the simulation time is short (∼1 Gyr) compared to the CR cooling time (∼10 Gyr) as shown by Enßlin et al. (2011).

3.2. Gravity

Although gravity fundamentally influences TI (Donahue & Voit 2022), the exact form of the gravitational potential is not critical for the subsequent evolution (Choudhury & Sharma 2016). We integrated a simple vertical gravity profile as proposed by McCourt et al. (2012):

$$g(z)=g_0\frac{z/a}{\sqrt{1+{\left(z/a\right)}^2}}{\hat{e}}_z=g(z){\hat{e}}_z,$$(5)

where z is the vertical distance from the midplane, and a=H/10 serves as a softening parameter for gravity that ensures that the profile transitions smoothly to zero near the centre of the simulation box, that is, for |z|<a. At large radii, the gravitational acceleration reaches a nearly constant value, determined by $g_0=c_{\text{iso}}^2/H$, where c_iso is the initial isothermal sound speed at the midplane. The corresponding free-fall time for a gas parcel at altitude z is

$$t_{\text{ff}}=\sqrt{\frac{2z}{g(z)}}.$$(6)

In our simulations, we do not account for self-gravity of the gas. This simplification is justified in CGM-like, low-density environments where the characteristic timescale for gravitational collapse of a self-gravitating gas cloud, t≃(Gρ)^−1/2, considerably exceeds the timescale for the collapse due to compression induced by TI. Self-gravity becomes important only if the mass of the initial density perturbation exceeds the Jeans mass (Li et al. 2020).

3.3. Cooling and heating

To establish a thermodynamic environment consistent with the conditions in the CGM, we use a realistic cooling function that is designed to model the relevant cooling processes. We account for three main cooling channels: Lyman-α cooling by H and He atoms, cooling from various metal lines, and bremsstrahlung cooling due to free-free emission. We capture these processes in the collision-ionisation equilibrium (CIE) approximation and compile a look-up table with the resulting cooling rates. This table is then used to calculate the cooling rate in the simulation. To model the Lyman-α cooling process, we solve the ionisation equilibrium of H and He following Cen (1992) and account for the ionisation due to an extragalactic ultraviolet background (UVB). We use the redshift z=0 UVB of Puchwein et al. (2019) and the ionisation cross-sections of Verner et al. (1996) to calculate the respective ionisation rates of H and He. The Lyα cooling rates are taken from Cen (1992). For metal line cooling, we calculate CIE cooling rates using Chianti (Dere et al. 1997) for the elements C, N, O, Si, Ne, Fe, and Mg assuming solar metallicity using the abundances of Asplund et al. (2009). These atoms are selected because they are the strongest metal coolants. Cooling by bremsstrahlung due to free-free emission from free electrons is calculated using the description of Ziegler (2018). The effective cooling rate in our simulation can be formulated as

$$\mathcal{C}=n_{\text{H}}^2\mathrm{\Lambda }(n_{\text{H}},T),$$(7)

where the volumetric cooling function, Λ(n_H, T), compiles all the information regarding radiative cooling and heating processes outlined above. We introduce a cut-off temperature for the cooling function at T_cut=10⁴ K, which is achieved by employing a tanh profile that causes the cooling function to smoothly transition to zero at around T_cut across a temperature range of ΔT=100 K. Gas with temperatures below T_cut remains in a cool state and undergoes no reheating due to any physically relevant processes.

Observational data suggest that the CGM is close to a state of global equilibrium (Werk et al. 2014; Miller & Bregman 2015). However, this balance is likely dynamic (Tumlinson et al. 2011) rather than static. These findings imply a global heating mechanism that prevents the CGM from a cooling catastrophe. Although, the exact nature of the heating process remains uncertain, Zhuravleva et al. (2014) validated its plausibility by demonstrating that turbulent heating can offset radiative cooling across a range of radii in both the Perseus and Virgo clusters. Their work provided observational evidence that turbulent dissipation may serve as a significant energy source, potentially stabilising the intracluster medium against runaway cooling. This mechanism is also relevant on galactic scales, where turbulence, driven by SN feedback, AGN outflows, and cosmic inflows, plays a similar role in balancing radiative losses (Voit 2018; Buie et al. 2020; Faucher-Giguère & Oh 2023). To address these observations, we introduce a heating mechanism that counterbalances cooling globally while permitting deviations from thermal equilibrium locally (McCourt et al. 2012). We subdivided the vertical extent of the simulation box into N/2 bins of equal height; collected the energy lost by radiative cooling in each of these bins at every time step,

$$\mathrm{\Delta }{E}_{\mathrm{th},\mathrm{bin}}=\sum m_i\mathrm{\Delta }{u}_i;$$(8)

and redistributed this energy inside the same bin by applying a mass-weighted scheme:

$$\mathrm{\Delta }{E}_{\text{th},i}=\frac{m_i}{M_{\text{bin}}}\mathrm{\Delta }{E}_{\mathrm{th},\mathrm{bin}}.$$(9)

Here, u_i denotes the internal specific energy, and m_i is the mass of cell i. The term M_bin is the cumulative mass of the bin containing the cell i.

3.4. Initial conditions

The simulation domain models a gas column symmetrically extending from the “galactic midplane” far into the CGM. The computational box has the dimensions 1H×1H×6H along the x, y, and z axes, where H=30.11 kpc in our fiducial run. We place the origin of the coordinate system in the centre of the box (cf. Fig. 1). The x and y boundaries are periodic while the boundaries in z-direction allow gas to exit the computational domain. Here, H denotes the characteristic scale height of the plasma, determined such that the ratio of cooling time to free-fall time, τ=t_cool/t_ff, is equal to the desired value at z=±H. The initial grid is uniform with a fiducial resolution of 128×128×768 gas cells. To assess the convergence of our simulations, we conduct a subset of runs at resolutions corresponding to a quarter, half, and twice the fiducial value (see Table 1).

Density profile. The initial gas distribution in each simulation follows an isocooling stratification (Butsky et al. 2020), implying the initial cooling time to be constant throughout the domain:

$$t_{\text{cool}}=\frac{{\varepsilon }_{\text{th}}}{n_{\text{H}}^2\mathrm{\Lambda }(n_{\text{H}},T)}=\mathrm{const}.$$(10)

Combining this condition with the assumption of hydrostatic equilibrium (HSE), dP_tot/dz=−ρ(z) g(z), we derived the differential temperature profile for an isocooling atmosphere:

$$\frac{\text{d}T}{\text{d}z}=\frac{\mu m_{\text{p}}}{\eta k_{\text{B}}}g(z){\left(1+\frac{1-{\mathrm{\Lambda }}_T}{1+{\mathrm{\Lambda }}_{n_{\text{H}}}}\right)}^{-1},$$(11)

where μ is the mean molecular weight, m_p denotes the proton mass, and k_B is Boltzmann's constant. The quantities Λ_T=∂lnΛ/∂lnT and ${\mathrm{\Lambda }}_{n_{\text{H}}}=\partial ln\mathrm{\Lambda }/\partial ln{n}_{\text{H}}$ describe the logarithmic derivatives of the cooling function with respect to temperature and density, respectively. η=1+X_{cr, 0}+X_{mag, 0}+X_{kin, 0} accounts for the contributions of thermal, CR, magnetic, and kinetic pressure to the total pressure, P_tot=ηP_th. We integrate Eq. (11) to obtain the vertical temperature profile and use the condition in Eq. (10) to calculate the corresponding pressure profile for an isocooling atmosphere in HSE. We present some of these initial density profiles for varying X_{cr, 0} in Fig. 3.

Fig. 3. Initial density profiles of the simulation box for various Xcr, 0. The solid lines show the implemented isocooling density profile, the dashed lines show the corresponding isothermal density profile with T0=106 K for comparison. For a higher initial Xcr, 0, more mass is contained in the atmosphere because the additional CR pressure supports the HSE without contributing to the gravitational force.

Velocity field. The onset of TI requires the presence of local density perturbations. Previous studies have typically introduced such fluctuations by embedding isobaric density variations within the HSE (McCourt et al. 2012; Butsky et al. 2020). In contrast, our approach is based on observational evidence of turbulence in the CGM, as indicated by the broadening of absorption lines in quasar spectra, which suggests gas cloud motions (Tumlinson et al. 2013), and by Mg ii absorption systems, which further support the presence of turbulence in the CGM (Huang et al. 2016). Turbulent motions generate converging flows that change the local density of the plasma on short timescales. To more accurately model the formation of thermally unstable regions, we initialise our simulations with a turbulent velocity field that follows a Kolmogorov spectrum (Kolmogorov 1941) within the range 4≤k H/2π≤256, which transitions to white noise on larger scales. We adjust this velocity field to ensure a kinetic-to-thermal pressure ratio of X_kin=P_kin/P_th=0.3 at every height. In our fiducial CR-dominated setup (described below), this method yields a velocity dispersion of approximately σ_v∼33 km s⁻¹ near the centre and σ_v∼30 km s⁻¹ at the edges of the simulation box.

Magnetic field. Turbulence amplifies the magnetic field predominantly via the small-scale dynamo mechanism (Brandenburg & Subramanian 2005; Pakmor et al. 2017; Pfrommer et al. 2022) and determines the orientation of the magnetic field lines. To encapsulate the dynamics between turbulent gas motion and magnetic fields in our simulations, we introduce a turbulent magnetic field that mirrors the scaling of the previously mentioned velocity field. Recent numerical investigations of the CGM in cosmological galaxies by Pakmor et al. (2020) have shown an approximately constant relative magnetic pressure, X_mag=P_mag/P_th, within the virial region. Thus, we vary the magnetic field strength in height following the trend of the thermal pressure profile to replicate this simulated trend. A Helmholtz decomposition is employed to ensure the magnetic field's divergence-free nature by eliminating its compressive component. This procedure yields a purely solenoidal magnetic field, which is subsequently normalised by a constant factor to ensure an average relative magnetic pressure of X_mag=0.01 throughout the simulation domain.

Cosmic rays. We include CR pressure, P_cr, in the initial conditions as a constant fraction of the gas pressure, P_th, in every computational cell, characterised by the parameter X_cr=P_cr/P_th. Introducing this additional non-thermal pressure modifies the isocooling density profile compared to a setup without CRs. Specifically, the non-thermal pressure gradient necessitates an adjustment in the gravitational force to maintain equilibrium, resulting in an atmosphere with a higher mass content. Figure 3 shows the resulting initial density profiles for various choices of X_{cr, 0}. We explore a range of initial values for X_{cr, 0} spanning 0, 0.03, 0.3, 3, and 30 to characterise different CGM environments, ranging from purely thermal to slightly CR-supported to strongly CR-dominated atmospheres. For runs modelling CR transport by pure diffusion, we vary the constant CR diffusion coefficient, κ₀, to simulate slow (3×10²⁷ cm² s⁻¹), moderate (3×10²⁸ cm² s⁻¹), and fast (3×10²⁹ cm² s⁻¹) diffusion rates. In simulations employing the two-moment approach for CR transport, CRs initially stream along the magnetic field at the Alfvén speed, v_a, characterised by the initial CR flux density, f_cr=v_a(ɛ_cr+P_cr). We initialise the energy density of Alfvén waves with ɛ_a=10⁻³ɛ_cr.

Numerical values. In the midplane of our simulation box, we initialised the values for gas density and temperature with ρ₀=2×10⁻²⁷ g cm⁻³ and T₀=10⁶ K, respectively. These parameters in combination with the isocooling condition yielded a global cooling time for the gas of t_cool≈108 Myr. Given our fiducial ratio of cooling time to free-fall time, τ=0.3, we derive a scale height of H=30.11 kpc. We summarise the parameters for the different simulations in Table 1.

4. Suppression of thermal instability by advective cosmic rays

We start our investigation of the impact of CRs on TI by examining purely advective CR transport, which represents the limiting case of inefficient anisotropic CR transport. In this scenario, CRs are exclusively advected with the thermal gas and have no motion relative to the medium. CR transport can be approximated by such advection models when the timescale for CR transport is significantly longer than the timescales governing other relevant physical processes. This situation is partially realised in environments where strong bulk gas flows are present, such as at the launching sites of galactic outflows or within supernova remnants (Ruszkowski & Pfrommer 2023; Armillotta et al. 2024). All simulations in this section are initialised with τ=0.3, X_{mag, 0}=0.01, X_{kin, 0}=0.3, and different values for X_{cr, 0}. We run each simulation until t=8 t_cool.

Figure 4 illustrates how the pressure support of advective CRs changes the morphology of cold gas in the CGM. We show projections in y-direction of the hydrogen column density, N_H, after the simulations run for 7 t_cool. From left to right, we systematically vary the initial CR pressure, X_{cr, 0}, ranging from purely thermal (X_{cr, 0}=0) to CR-pressure dominated (X_{cr, 0}=3) scenarios. As the gas cools, denser regions lose their thermal pressure support more rapidly than their surroundings. These denser gas parcels sink and undergo compression due to the rising ambient pressure of their vicinity. In the simulation without CRs, pressure equilibrium is approximately sustained due to the short cooling time (τ=0.3). This leads to ongoing compression until the gas temperature reaches a new equilibrium. However, the presence of CR pressure changes this scenario considerably. Even a minor CR pressure support of X_{cr, 0}=0.03 clearly alters the morphology of the collapsing regions. The inability of CRs to escape the contracting clouds in this model leads to an elevated non-thermal pressure support that counteracts the collapse and, with an increasing contribution of CRs, gradually limits the extent to which the gas can be compressed.

Fig. 4. Projections of the hydrogen column density, NH, employing purely advective CR transport. The panels show the x−z plane of a section extending from z=0 to z=2H and the depth of the projection is 1 H. The snapshot is at t=7 tcool≈2 tff. All simulations were initialised with different relative CR pressure, Xcr, 0. The density of the clouds decreases with increasing CR pressure, while the cloud size increases. In the CR pressure dominated run with Xcr, 0=3, cold gas condensation is entirely suppressed. This illustrates that advective CRs alter the morphology of the CGM significantly and suppress the formation of cold clouds.

The onset of TI is completely suppressed in the CR-pressure dominated run (X_{cr, 0}=3). While this suppression correlates with increased CR pressure, it is not the CR pressure alone that causes it. Instead, the suppression arises from the method used to generate the initial density perturbations required to trigger collapse. In all simulations, the initial ratio of kinetic to thermal pressure, X_{kin, 0}=P_{kin, 0} / P_{th, 0}, is fixed to 0.3, regardless of the CR pressure. As X_{cr, 0} increases, the thermal pressure makes up a smaller fraction of the total pressure, and thus the same X_{kin, 0} corresponds to a smaller absolute level of kinetic energy relative to the total pressure. This results in significantly weaker initial density perturbations compared to runs with lower X_{cr, 0}. These weaker perturbations are more easily smoothed out by turbulent heating introduced via the initial velocity field, further reducing fluctuations and ultimately suppressing the development of TI.

This behaviour is clearly illustrated in Fig. 5, which presents (from left to right) the density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux for simulations with purely advective CR transport. The left panel shows that the stabilising effect of CRs becomes apparent as the initial density fluctuations decrease with increasing X_{cr, 0}. Notably, the cold mass fraction converges to nearly identical values across all simulations permitting TI, as expected, since the cooling function depends solely on gas density and temperature and is independent of CR pressure. Once a region enters TI, the subsequent runaway cooling proceeds unaffected by CRs. However, the volume of cold gas increases by a factor of 10−20 compared to the non-CRs run, which is driven by the additional CR pressure counteracting the compression of thermally unstable gas. This influence is further underscored by the cold mass flux analysis: the increased volume and reduced density of cold gas result in stronger buoyancy forces, which reduce the cold mass flux in simulations with higher CR pressure support.

Fig. 5. Time evolution of cold gas metrics for simulations employing advective CR transport with varying Xcr, 0. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. We evaluated each computational cell in the region 0.25 H≤|z|≤2.75 H. The dashed vertical line marks the simulation time corresponding to the results shown in Fig. 4.

In previous work, Butsky et al. (2020) found that the onset of TI is entirely independent of CR pressure, which appears to contradict our results. However, we highlight that our approach to generating density perturbations differs significantly from theirs. While Butsky et al. (2020) initialise their simulations with a hydrostatic atmosphere containing small, pre-seeded density inhomogeneities, we evolve these perturbations self-consistently from an initial turbulent velocity field. Each procedure performs best for the distinct aspects of TI analysis: the method of Butsky et al. (2020) optimally captures the linear growth of TI whereas our more dynamical approach aims to mimic the more realistic evolution of astrophysical systems subject to turbulence.

However, the idealised scenario of purely advective CR transport only partially applies in realistic environments, where CRs usually can actively diffuse or stream along magnetic field lines. The transition from purely advective to active CR transport fundamentally alters the CR pressure distribution, reduces the stabilising effect of CRs, and reshapes the dynamics of TI. We address this topic in the following section.

5. Revival of thermal instability by self-confined cosmic rays

In realistic astrophysical environments, CRs interact with the surrounding medium in more complex ways than simply being advected with the gas. Active CR transport in the form of diffusion and streaming allows CRs to decouple from the gas and move along magnetic field lines at speeds that can exceed the local gas velocity. In the following, we discuss the implications of active CR transport on the formation of cold gas.

We begin our analysis with a general overview of our simulations. In Fig. 6, we present a gallery showcasing slices of various quantities from the simulation utilising two-moment CR transport with an initial relative CR pressure of X_{cr, 0}=3. The top row shows various gas quantities, namely density ρ, temperature T, cooling time τ_cool, and gas z-velocity v_z. The bottom row displays the CR pressure P_cr, the relative CR pressure X_cr=P_cr/P_th, the relative magnetic pressure X_mag=P_mag/P_th, and the effective CR diffusion coefficient κ_eff discussed below (see Sect. 6.1). This indicates that active CR transport induces substantial changes in the condensation process of cold gas compared to the advective CR scenario (cf. Fig. 4).

Fig. 6. Gallery showing 2D slices of various quantities from the simulations employing two-moment CR transport with an initial CR pressure fraction of Xcr, 0=3. All snapshots are at t=7 tcool. The top row illustrates quantities of the thermal gas, while the bottom row presents variables related to CR transport.

Starting from isocooling initial conditions, a two-phase medium emerges, characterised by a volume-filling, dilute hot phase (ρ≲10⁻²⁷ g cm⁻³, T≳10⁶ K) and a dense cold phase (ρ≳10⁻²⁵ g cm⁻³, T∼10⁴ K), the latter being composed of numerous small (∼100 pc) clouds with very short central cooling times (t_cool≲100 kyr). These overdense cold clouds fall towards the centre of gravity with high velocities so that the densest clumps reach infall speeds of ≳90 km s⁻¹. Interestingly, the upwards motions are caused by turbulence as well as the motions induced by our assumed mass-weighted heating that is globally offsetting radiative cooling. The CR pressure exhibits a relatively uniform distribution across the box (we point out the linear colour scale) with moderate peaks within and around the collapsing regions due to adiabatic compression of CRs. This compression effect also applies to the magnetic pressure, represented by X_mag, which reaches values of about one to eight in and around the cold clouds as a result of the flux-frozen characteristic of the magnetic field. The range in X_mag is much larger in comparison to that of X_cr, which is a direct consequence of fast CR transport out of the collapsing clumps. In fact, the effective CR diffusion coefficient spans an immense range of more than seven orders of magnitude with values between κ_eff∼10²⁷ cm² s⁻¹ and κ_eff>10³² cm² s⁻¹, with the largest values realised in the converging regions, which are also characterised by large values of X_mag.

Having identified substantial differences between outcomes under initially dominating CR pressurisation between the two-moment setup and the purely advective CR scenario, we extend our analysis to include a wider range of initial CR pressure values. To this end, we repeat the simulation suite from Sect. 4 and include a strongly CR-pressure-dominated atmosphere with X_{cr, 0}=30. All simulations are initialised with τ=0.3, X_{mag, 0}=0.01, X_{kin, 0}=0.3, and are evolved for eight cooling times.

Figure 7 demonstrates that under realistic conditions of active CR transport, variations in the initial CR pressure have minimal impact on the morphology of cold gas in the CGM. The figure shows gas mass density slices at t=7 t_cool, with panels illustrating a systematic variation in X_{cr, 0} from purely thermal conditions (X_{cr, 0}=0) to strongly CR pressure-dominated regimes (X_{cr, 0}=30). Remarkably, the number and size of condensing clouds remain nearly unchanged across all simulations, showing minimal dependence on the initial CR pressure. Even under the influence of large CR pressures (X_{cr, 0}=30), TI progresses without significant suppression.

Fig. 7. Slices showing the mass density in the x−z plane centred at y=0 of simulations employing two-moment CR transport and different relative CR pressures, Xcr, 0. All snapshots are at t=8 tcool. The initial CR pressure barely influences the morphology of the cold gas in terms of cloud size or density.

Figure 8 further supports these findings, showing the temporal evolution of cold gas metrics. As described in Sect. 4, the initial density fluctuations decrease with increasing X_{cr, 0} as a consequence of the reduced influence of the initial kinetic pressure, which is constant (X_{kin, 0}=0.3) across all runs. This effect is particularly relevant in setups where the CR pressure is significantly higher than the kinetic pressure, that is, in simulations with X_{cr, 0}≥3. Consequently, cold gas forms delayed (t∼4 t_cool) compared to thermally dominated runs (t∼2 t_cool). After the simulations evolved for t∼7 t_cool, nearly all metrics converge to similar values across the simulations. The only exception is the cold mass flux in the run with X_{cr, 0}=30, which remains offset by approximately a factor of five in comparison to the other simulations. We attribute this effect to the exceptionally high CR pressure, which counteracts gravity and consequently reduces the infall velocity of the cold clouds.

Fig. 8. Time evolution of cold gas metrics for simulations employing two-moment CR transport with varying Xcr, 0. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. We evaluate each computational cell in the region 0.25 H≤|z|≤2.75 H. All cold gas metrics show a very similar evolution and are almost independent of the initial CR pressure.

Our analysis demonstrates that the amount of CRs has only a small impact (within a factor of two) on the final amount of cold gas or its volume filling fraction, see Fig. 8. This finding is somewhat unexpected, because CR pressure is typically thought to play a significant role in shaping gas dynamics by contributing to additional non-thermal support and influencing cooling instabilities (Butsky et al. 2020). The lack of substantial morphological changes suggests that the interplay between CR transport, gas dynamics, and cooling is governed by other processes. We attribute this outcome to the high CR transport speed, as reflected by the high values of the effective CR diffusion coefficient (see bottom right panel of Fig. 6), which enables CRs to escape collapsing regions before significantly influencing TI. In the following section, we explore this behaviour further by examining the effect of varying CR transport speeds, focusing specifically on changes in the effective CR diffusion coefficient, κ_eff, and its impact on the development of TI.

6. Impact of cosmic ray transport speed on thermal instability

The results from the previous section suggest that CRs escape too quickly from collapsing regions to meaningfully affect the condensation of cold gas. Therefore, in this section, we explore the significance of the CR transport speed and begin by analysing how CRs are transported in the presence of TI.

Relative to astrophysical plasmas, the CR distribution primarily propagates via streaming and diffusion along magnetic field lines. In the streaming regime, CRs travel at the local Alfvén speed. As they propagate, they excite Alfvén waves that pitch-angle scatter the CRs. In turn, these plasma waves lose their energy through wave damping processes. On the other hand, (anisotropic) diffusion describes CR transport along magnetic field lines down their pressure gradient. This can be described as a random walk with a single step representing a CR pitch-angle scattering event of off small-scale magnetic fluctuations. In 1-moment CR transport models, both processes are assumed to be in steady-state: in the ISM, the CR scattering rate is typically constant and independent of local physical conditions, producing “canonical” values of κ_cr∼3×10²⁷⁻²⁸ cm² s⁻¹ (Ruszkowski & Pfrommer 2023), while perfect CR streaming assumes a tight coupling between CRs and gyroresonant Alfvén waves, which results in a streaming velocity equal to the Alfvén speed. However, these assumptions hold only in environments where Alfvén-wave generation and damping are perfectly balanced so that the energy in these waves is (i) constant and maintains a fixed scattering rate, and (ii) sufficient to continuously couple CR motion to the waves. In realistic environments, these two conditions are not permanently maintained. The generalisation to more realistic CR transport is provided by the two-moment CR transport approach.

6.1. Cosmic ray transport in the two-moment picture

The CR transport in the framework of CRMHD is also characterised by a combination of streaming and diffusion. The CRs stream along magnetic field lines at the CR streaming speed (Thomas & Pfrommer 2019)

$${\mathit{v}}_{\text{st}}={\mathit{v}}_{\text{a}}\frac{{\nu }_{+}-{\nu }_{-}}{{\nu }_{+}+{\nu }_{-}},$$(12)

where v_a denotes the local Alfvén speed while ν₊ and ν₋ are the scattering rates of CRs with respect to co- and counter-propagating Alfvén waves, respectively. If one of these wave types is damped away and its corresponding scattering frequency vanishes, this equation approaches the limiting case of 1-moment streaming where CRs stream down their pressure gradient, ${\mathit{v}}_{\text{st}}=-{\mathit{v}}_{\text{a}}\text{sign}\left(b·\nabla P_{\text{cr}}\right)$ (Zweibel 2013; Pfrommer et al. 2017; Thomas et al. 2023).

When the transport velocity of CRs exceeds v_a, the excess velocity is attributed to CR diffusion relative to the Alfvén waves. This diffusive component is quantified by the intrinsic CR diffusion coefficient (Thomas & Pfrommer 2019):

$${\kappa }_{\pm }=\frac{c^2}{3{\nu }_{\pm }}=\frac{8}{9\pi }\frac{\gamma m{c}^3}{\mathit{eB}}\frac{B^2}{{\varepsilon }_{\mathrm{a},\pm }},$$(13)

with its total value given by

$${\kappa }_{\text{cr}}={\left({\kappa }_{+}^{-1}+{\kappa }_{-}^{-1}\right)}^{-1},$$(14)

where the + and − signs of the subscripts correspond to quantities relative to co-propagating and counter-propagating Alfvén waves, respectively. Here, B is the magnetic field strength, e is the elementary charge, c is the speed of light, and γ=2 represents the typical Lorentz factor of giga-electronvolt CRs with rest mass m. The energy density content of Alfvén waves, ɛ_a, controls the level of CR diffusivity: in regions of high ɛ_a, CRs scatter frequently and are tightly coupled to the Alfvén waves, reducing their effective diffusive transport. Conversely, in regimes with low ɛ_a, CRs adopt higher diffusion coefficients, making their transport predominantly diffusive.

In Fig. 9, we show the CR distribution in the total intrinsic CR diffusion coefficient, κ_cr, and X_cr at t=6 t_cool. Simulations that are initialised with different values of X_{cr, 0} result in X_cr values that range from approximately 0.5 to 10 times their initial relative CR pressure. Two distinct trends emerge from the figure. First, the intrinsic diffusion coefficients span an extensive range, from the slow-diffusion regime with κ_cr∼2×10²⁷ cm² s⁻¹ to values exceeding 10⁴⁷ cm² s⁻¹. We observe that the fraction of CRs in the high-diffusivity regime increases with increasing X_{cr, 0} and even dominates the CR population in simulations with CR-pressure dominated atmospheres (X_{cr, 0}≥3). This is a result of the stabilising effect of the increasing CR pressure which smooths the CR pressure distribution and reduces the CR gradients. This suppresses the generation of Alfén waves and consequently reduces the CR scattering rate. According to Eq. (13), CRs with such high diffusion coefficients scatter, on average, less than once per gigayear. Since the duration of our simulations is of order this timescale (t_sim∼1 Gyr), these CRs are currently not scattered by Alfvén waves and the CR force parallel to the magnetic field line is suppressed. We note that this is a temporary state because these CRs can migrate to regions with lower diffusion coefficients, where they resume to scatter with Alfvén waves and exert the parallel forces again.

Fig. 9. Intrinsic CR diffusion coefficient, κcr, for simulations utilising two-moment CR transport with varying initial values of Xcr, 0. Scatter points represent the distribution of CRs as a function of κcr and Xcr at t=6 tcool, with darker colours indicating higher CR-energy-weighted probability density in a logarithmic colour scheme. The corresponding 1D probability density of κcr is displayed with a linear scaling on the right-hand side. In the regime of efficient CR scattering (small κcr), larger Xcr, 0 values increase the overall CR energy density in the CGM allowing for more CR energy to be converted into gyroresonant Alfvén waves and larger CR scattering rates, thereby reducing κcr (see Eq. (13)). In regions of a smooth CR distribution, the driving of Alfvén waves is inefficient in comparison to wave damping processes so that CRs are poorly coupled to the plasma, leading to a very large value of κcr.

Second, in regions of efficient CR scattering, the distribution of κ_cr shifts toward smaller values with increasing X_cr. This behaviour is in line with the self-confinement picture of CR transport: in regions with higher X_cr, the greater abundance of CRs amplifies the generation of gyroresonant Alfvén waves. These waves, in turn, increase CR scattering rates, confining CRs more effectively and thereby reducing their effective diffusivity. This interplay between CR abundance, wave generation, and scattering underscores the critical role of self-confinement in regulating CR transport properties across the CGM.

As stated above, κ_cr only captures the diffusive motion of CRs while ignoring their streaming motion. Therefore, it is not comparable to diffusion coefficients that arise from diffusion-only models. To mitigate this, we introduced the effective CR diffusion coefficient (Thomas et al. 2023),

$${\kappa }_{\text{eff}}=\frac{f_{\text{cr}}}{b·\nabla {\varepsilon }_{\text{cr}}},$$(15)

which characterises the entire CR transport, streaming or diffusing, using a single quantity. Here, the term b·∇ɛ_cr is the gradient of CR energy density along the magnetic field and f_cr is the CR energy flux that is evolved in our simulations as one of the CRMHD quantities. Defining such an effective CR diffusion coefficient is useful because in the limit of a purely diffusive CR transport, both the intrinsic and effective diffusion coefficient are the same (i.e., κ_cr=κ_eff). In general, the effective quantity κ_eff is the diffusion coefficient needed to mimic the actual CR transport in a hypothetical pure-diffusion model. Figure 10 illustrates the two-dimensional CR distribution in κ_eff and X_cr at t=6 t_cool. The distribution spans an immense range, from low values around κ_eff∼10²⁵ cm² s⁻¹ to highly diffusive CR populations with κ_eff∼10³⁵ cm² s⁻¹, covering ten orders of magnitude. There is a noticeable trend to higher values of κ_eff for increasing X_cr.

Fig. 10. Same as Fig. 9 but for the effective CR diffusion coefficient, κeff.

To represent this finding with a single metric, we examine the CR-energy-weighted effective diffusion coefficient, ${\overline{\kappa }}_{\text{eff}}$, which emphasises the dominant contribution of high-energy CRs. At each snapshot between 1 t_cool and 7 t_cool, we determine the CR-energy-weighted median of κ_eff and we compute the mean across these values. Even for simulations with moderate CR pressure support (i.e., X_{cr, 0}=0.03 and X_{cr, 0}=0.3), the value of ${\overline{\kappa }}_{\text{eff}}$ is increased by a factor of ∼10 in comparison to the “canonical” ISM CR diffusion coefficient in diffusion-only models. For the two cases, we find ${\overline{\kappa }}_{\text{eff}}=1.3\times {10}^{29}{\mathrm{cm}}^2\hspace{0.5em}{\mathrm{s}}^{-1}$ and ${\overline{\kappa }}_{\text{eff}}=1.6\times {10}^{29}{\mathrm{cm}}^2\hspace{0.5em}{\mathrm{s}}^{-1}$, respectively. The CR-dominated runs with X_{cr, 0}=3 and X_{cr, 0}=30 result in even higher values of ${\overline{\kappa }}_{\text{eff}}=4.3\times {10}^{29}{\mathrm{cm}}^2\hspace{0.5em}{\mathrm{s}}^{-1}$ and ${\overline{\kappa }}_{\text{eff}}=8\times {10}^{29}{\mathrm{cm}}^2\hspace{0.5em}{\mathrm{s}}^{-1}$, respectively. In the following, we utilise ${\overline{\kappa }}_{\text{eff}}$ to examine the impact of the effective CR transport speed on TI.

6.2. Non-Fickian cosmic ray diffusion

The simulations in the previous section employed the two-moment framework for CR transport and showed a highly variable effective CR diffusion coefficient and CR transport speed that depend on their environment in the TI-affected CGM. Quantifying the impact of the CR transport speed on the TI-induced cloud collapse becomes difficult due to the simultaneous presence of both transport modes. To simplify the analysis while maintaining the possibility for direct comparison, we conduct a series of simulations employing purely diffusive CR transport, with a constant CR diffusion coefficient, κ₀. A standard approach to describe such diffusion-only models is through Telegrapher's equation (Gombosi et al. 1993; Litvinenko & Schlickeiser 2013; Thomas & Pfrommer 2019):

$${\tau }_{\text{cr}}\frac{{\partial }^2{\varepsilon }_{\text{cr}}}{\partial t^2}+\frac{\partial {\varepsilon }_{\text{cr}}}{\partial t}=\nabla ·\left({\kappa }_{0}bb·\nabla {\varepsilon }_{\text{cr}}\right),$$(16)

where κ₀ is the constant CR diffusion coefficient, and τ_cr=3κ₀/c² defines the CR diffusion timescale. The Telegrapher's equation extends the usual Fickian diffusion model by introducing a second-order time derivative, effectively accounting for the finite propagation speed of CRs. Unlike the purely diffusive approach, which assumes instantaneous diffusion, the Telegrapher's equation describes CR transport as a advection-like propagation combined with diffusion. This non-Fickian description of CR transport (Malkov & Sagdeev 2015; Litvinenko & Noble 2016; Rodrigues et al. 2019) results in a hyperbolic equation, for which a numerical discretisation in the context of finite volume schemes is more straightforward in comparison to the parabolic anisotropic diffusion equation. Furthermore, we adapted our solver for the two-moment equations to simulate this constant diffusion coefficient model. This allows for a straightforward comparison between simulations using the non-Fickian approximation with a constant diffusion coefficient and the two-moment self-confinement model. The underlying physical approximation of the two models is different, and there is no satisfactory explanation for the constant diffusion model in the self-confinement picture of CR transport. The absence of explicit coupling to Alfvén waves means that scattering must be assumed to originate from an external process, not driven by the CRs themselves. Thus, in the diffusion picture, the motion of CRs along magnetic field lines is in the direction opposite to the CR pressure gradient (see the right-hand side of Eq. (16)).

To investigate the influence of CR transport speed on TI, we conducted a subset of simulations from Sect. 5, this time utilising the diffusion-only model for the transport of CR energy. As before, all simulations are initialised with τ=0.3, X_{mag, 0}=0.01, and X_{kin, 0}=0.3, and each run is evolved for eight cooling times. Additionally, we use three different values for the constant CR diffusion coefficient each representing a different regime: slow diffusion with κ₀=3×10²⁷ cm² s⁻¹, the “canonical” value (κ₀=3×10²⁸ cm² s⁻¹), and fast diffusion with κ₀=3×10²⁹ cm² s⁻¹, and compare the results to a run employing two-moment CR transport.

Figure 11 demonstrates that fast and active CR transport accelerates the onset of collapse. The figure shows projections of the hydrogen column density, N_H. All simulations are initialised with a relative CR pressure of X_{cr, 0}=3, but different CR transport mechanisms. The snapshots are taken at t=6 t_cool. There is a consistent trend of increased gas condensation with higher (effective) CR diffusion coefficient. In the purely advective case, no cold gas forms, as discussed in Sect. 4. However, cold gas starts to precipitate already for slowly diffusing CRs. This shows that as soon as there is an escape channel for CRs out of the collapsing TI perturbations, these perturbations lose the CR-pressure support that is present in the overly confining CR-advection models.

Fig. 11. Projections of the hydrogen column density, NH, for simulations utilising different CR transport models. All runs were initialised with Xcr, 0=3, and the snapshots are at t=7 tcool. The figure demonstrates that the onset of collapse is governed by the CR transport velocity rather than the CR pressure.

As the (effective) CR diffusion coefficient increases, the simulations show a gradual increase in dense gas formation. This condensation is fastest in the two-moment transport model. Notably, the diffusion model with κ₀=3×10²⁹ cm² s⁻¹ produces results comparable to those of the two-moment transport setup. We attribute this similarity to the close match of the constant diffusion coefficient in this simulation to the effective diffusion coefficient of the two-moment model, κ_eff=4.3×10²⁹ cm² s⁻¹. Consequently, both models exhibit analogous CR transport behaviour and cold gas condensation patterns.

Figure 12 confirms these findings in a more quantitative way. In each of the four metrics we observe a clear transition from the two-moment model, which represents the scenario of fastest CR transport, to the advective model, constituting fully ineffective CR transport. The two-moment model exhibits a nearly identical evolution to the fast-diffusion model across all metrics. By t=6 t_cool, both models have generated significantly more cold gas compared to the setups with slower diffusion. Specifically, they produce approximately ten times more cold gas than the model with κ₀=3×10²⁷ cm² s⁻¹ and five times more than the model with κ₀=3×10²⁸ cm² s⁻¹. This further supports our earlier finding that higher CR transport velocities play a critical role in facilitating cold gas condensation.

Fig. 12. Time evolution of cold gas metrics for simulations with different CR transport models and Xcr, 0=3. Metrics are evaluated for computational cells within 0.25 H≤|z|≤2.75 H. The dashed vertical line marks the time corresponding to the snapshot shown in Fig. 11. Lower CR transport velocities delay the onset of cloud collapse and reduce the amount of condensed cold gas.

6.3. Escape of cosmic rays from collapsing clouds

Cosmic rays can counteract TI by providing additional pressure support but only as long as they remain trapped within the collapsing regions so that they can build up a pressure gradient. Therefore, when the confinement time of CRs inside cold clouds is prolonged, their influence on TI is expected to be more pronounced. To quantify this effect, we defined the average time required for CRs to escape from a cold cloud as

$$t_{\text{cr}}=\frac{r_{\text{cloud}}^2}{\kappa },$$(17)

where r_cloud is the effective cloud radius (see Sect. 2.4 for the definition), and κ is either the effective CR diffusion coefficient in two-moment CRMHD (κ_eff) or the pure diffusion coefficient in our non-Fickian diffusion model (κ₀). For this definition, we continue to approximate the clouds as spheres and assume that CRs diffuse freely through the clouds with their respective diffusion coefficients.

In the absence of self-gravity, the collapse timescale of a cold cloud is typically governed by the cooling timescale, t_collapse∼t_cool, as the cloud's internal dynamics is primarily driven by radiative cooling and pressure balance rather than gravitational collapse. In this context, t_cool sets the pace at which the cold cloud loses energy and contracts due to thermal pressure imbalance. The cloud will collapse as long as it cools faster than it can regain pressure support. For typical conditions in the CGM, where cold clouds are embedded in a hot, diffuse medium, t_cool ranges from several kyr to hundreds of Myr, depending on the cloud's density, temperature, and metallicity.

To compute the cooling time of a cold cloud, we adopt its mean density, ρ_cloud=M_cloud/V_cloud, and mean internal specific energy, U_{th, cloud}=∑_i(m_i u_i)/M_cloud, as representative values. Here, M_cloud and V_cloud are the total mass and the total volume of the cloud, respectively, while m_i and u_i denote the mass and the internal specific energy of the computational cell i within the cloud. By averaging the cloud's thermodynamic properties, we account for both the faster cooling of the dense, cold core and the slower evolution of the outer layers, characterised by lower densities and higher temperatures in an average sense.

When t_cr<t_collapse, CRs escape from the cold clouds rapidly enough that their impact on the collapse process is minimal. Conversely, when CR transport occurs over timescales longer than the collapse time (t_cr≳t_collapse), CRs can set up a pressure gradient that points toward the cloud centre, which is able to significantly counteract TI to eventually halt cloud collapse.

Figure 13 illustrates the impact of varying CR transport velocities on CR-cloud interactions. The left panel shows the distribution of t_cr/t_collapse, while the right panel presents X_cr, both evaluated for various cloud radii. The analysis includes eight consecutive snapshots starting from the moment the cloud finder identifies the first cloud. Two consecutive snapshots are separated by Δt=10 Myr, which is significantly larger than the typical values of t_cool∼0.5 Myr and $t_{\text{cr}}\sim 0.02\left(3\times {10}^{29}{\mathrm{cm}}^2\hspace{0.5em}{\mathrm{s}}^{-1}/{\kappa }_{\text{eff}}\right)$ Myr. This indicates that the clouds in the analysed snapshots are statistically independent. Consequently, some clouds may appear multiple times in the dataset, each data point sampling a different episode in their evolution with a different radius or thermodynamic state. The grey dashed line in the left panel represents the threshold where the CR diffusion time equals the collapse time: clouds to the left of this line enable CRs to escape faster than the collapse time of the cloud, whereas clouds to the right can retain CRs, which then slow down (or even halt) cloud collapse.

Fig. 13. Ratios of the CR transport timescale to the cloud collapse timescale (left) and the relative CR pressure within cold clouds (right), depicted for various cloud radii. The 1D probability density of the corresponding variable is displayed above each axis with a linear scaling. We analyse cold gas clouds within the region −0.25H≤x, y≤0.25H and 0.25H≤z≤1.75H (We focus on this reduced region to avoid double counting clouds across periodic boundaries and to ensure consistent analysis across resolutions (see Sect. 7) because the memory-intensive cloud finder makes analysing the full domain computationally very expensive at high resolutions.), considering eight consecutive snapshots following the identification of the first cloud by the cloud finder. Different colours represent different CR transport models (as detailed in the figure legend). All simulations are initialised with Xcr, 0=3. The dashed vertical line in the left panel separates regions where clouds collapse is slower (tcr/tcollaspe<1) or faster (tcr/tcollaspe>1) than CRs escape these clouds. Fast transport enables CRs to escape collapsing clouds on timescales shorter than the cloud collapse time, reducing their relative pressure. In contrast, slow transport traps CRs within collapsing regions, increasing their pressure support.

The timescale ratio t_cr/t_collapse exhibits a broad range, spanning ≲10⁻³ to ≳20. This variability arises from three factors: first, differences in cloud collapse timescales, t_collapse, governed by internal thermodynamic properties for individual clouds; second, variations in cloud radii affecting $t_{\text{cr}}\propto r_{\text{cloud}}^2$; and third, in simulations with two-moment CR transport, a wide range of κ_eff values influences $t_{\text{cr}}\propto {\kappa }_{\text{eff}}^{-1}$.

Fast CR transport in the two-moment and diffusion model with κ₀=3×10²⁹ cm² s⁻¹, enables rapid CR escape and yields a median t_cr/t_collapse∼6×10⁻². For lower κ₀, slower CR escape increases the median ratio that reaches t_cr/t_collapse∼0.8 (∼7) for κ₀=3×10²⁸ cm² s⁻¹ (κ₀=3×10²⁷ cm² s⁻¹). This behaviour influences CR pressure support during collapse as shown in the right panel. Slow CR transport traps CRs and enhances the pressure difference between the cloud to its surroundings and delays collapse. In contrast, fast CR transport reduces the CR pressure gradient and lowers X_cr within clouds.

6.4. Cosmic ray motorways

While our previous analyses suggest that CRs can be transported rapidly enough to escape contracting clouds without exerting significant pressure to resist collapse, the exact nature of their escape pathways remains unclear. CR transport is inherently anisotropic as CRs primarily drift along magnetic field lines either by streaming or diffusion. The efficiency of CR escape therefore depends critically on the magnetic field topology in and around the collapsing region because open field lines facilitate rapid CR transport, whereas closed or convoluted lines impede it.

In Fig. 14, we visualise the magnetic field topology surrounding one of the cold clouds from our simulations (which represents a typical situation as we confirmed). The magnetic field strength is considerably amplified inside the cloud (∼1 μG) compared to the ambient gas (∼0.1 μG) because the magnetic field lines are flux-frozen into the collapsing gas: when a cloud contracts and its volume decreases then the field lines are compressed and the magnetic field strength increases in proportion to the decrease in cloud size. The magnetic field lines permeate the cloud, align with its morphology, and guide the CR motion. Within the collapsing cloud, these field lines remain open and provide a clear and direct pathway – essentially a “CR motorway” – from the cloud's interior to the ambient medium. This open field line structure facilitates the efficient escape of CRs and allows them to stream or diffuse outward along these magnetic field lines, bypassing the dense, collapsing gas. These direct paths enable CRs to rapidly escape, thus preventing significant pressure build-up inside the cloud, and reducing the CR-mediated influence on TI.

Fig. 14. Visual impression of the magnetic field topology in and around one of the clouds. This cloud has formed through TI and is currently falling down. The volume rendering highlights dense gas and the magnetic field lines are coloured according to the local magnetic field strength. The shown magnetic field lines connect the inside of the cloud to the ambient gas above the cloud. The magnetic field inside the cloud is stronger (achieving B∼1 μG) compared to the magnetic field in the ambient gas.

6.5. Impact of cosmic rays on gas phases

Figure 15 presents the mass-weighted n_H–T diagrams for the entire simulation suite. The columns represent different CR transport models and the rows correspond to varying initial values of X_{cr, 0}. Within each row, differences highlight the impact of CR transport mechanisms under identical initial conditions and reveal their influence on gas phases. In contrast, differences within a column show the effects of increasing CR pressure support alongside changes arising from different initial conditions. Above each column and row, we provide 1D probability density distributions for the respective variables. This enables a clear comparison of CR transport models at fixed X_{cr, 0} or across different X_{cr, 0} for the same transport model. The colour gradient represents the probability density within each histogram bin, with all snapshots analysed at t=6 t_cool. The results reveal several key trends:

Fig. 15. Mass-weighted nH–T diagrams from our simulation suite. Columns represent different CR transport methods, while rows correspond to varying initial CR pressure fractions, Xcr, 0. The 1D probability densities of the respective variables are shown above each axis: columns reflect the resulting nH distributions for a fixed CR transport method across different Xcr, 0 and rows show the resulting T distributions for a fixed Xcr, 0 across different CR transport methods. The phase diagrams are constructed using 70 bins of equal (logarithmic) width for both variables. All simulations are analysed at t=6 tcool.

Cosmic ray transport mechanisms. Provided the CR pressure inside a collapsed cloud starts to dominate, CRs in a purely advective CR transport model suppress gas condensation in comparison to active CR transport models. Active CR transport can enhance gas condensation by up to a factor of ∼20 relative to advective models. However, for marginal CR pressure support (X_{cr, 0}=0.03), both transport approaches produce similar temperature distribution profiles.

Impact of cosmic ray propagation speed. Faster CR transport promotes the formation of denser and colder gas. This effect is particularly pronounced in CR pressure-dominated atmospheres (X_{cr, 0}>1) with active CR transport, whereas differences in thermally dominated CGMs are minimal. Simulations using the two-moment model or a diffusion coefficient of κ₀=3×10²⁹ cm² s⁻¹ yield nearly identical results.

7. Critical role of resolution

In this section, we examine the effects of numerical resolution in TI simulations of the CGM and focus on its implications for cloud sizes and the resulting effects on CR-cloud interactions.

7.1. Theoretical cloud sizes

Theoretical models predict that cold gas clouds should have characteristic sizes determined by the cooling length, l_cloud∼c_s t_cool, where c_s is the local sound speed (Voit & Donahue 1990; McCourt et al. 2018). This relationship reflects the physical scale below which perturbations become thermally unstable and collapse into dense structures driven by cooling and compression by ambient pressure (instead of self-gravity) before significant thermal equilibration can occur. Under typical CGM conditions, this scale ranges from a few parsec to about a kiloparsec and depend on the gas density, temperature, and other environmental factors.

Resolving this scale in global simulations of galaxies or halos is particularly challenging because it often lies well below the achievable spatial resolution. However, low resolution significantly suppresses turbulence-driven fragmentation, which is a crucial mechanism for forming dense, compact clumps from thermally unstable gas (Scannapieco & Brüggen 2015; Sparre et al. 2019). High-resolution simulations of isolated CGM patches have shown that resolving small-scale interactions between cold and hot phases is critical for accurately capturing cloud survival and dynamics (Cooper et al. 2009; McCourt et al. 2015; Gronke & Oh 2018; Sparre et al. 2020). Without sufficient spatial resolution, simulations produce cold gas clouds that are overly smooth, diffuse, and unrealistically large when compared to observations. Fielding et al. (2020) simulated radiative cooling layers and demonstrated that lower resolution simulations can exhibit significant pressure dips at intermediate entropy levels where cooling rates are highest. These pressure dips vanish with increased resolution and are caused by numerical diffusion, which artificially broadens the cooling layer. This artificial mixing of cool and hot phases creates intermediate-temperature gas instead of distinct cold clouds (Hummels et al. 2019). This impacts predictions for the column densities of low-ionisation species that are essential for tracing cold CGM gas in observations (Peeples et al. 2019; van de Voort et al. 2019).

7.2. Effect of resolution on cosmic ray-cloud interactions

Resolution is equally crucial for accurately modelling CR transport because CR pressure gradients and interactions with gas and magnetic fields depend sensitively on small-scale structures. Low-resolution simulations often smooth out these gradients, which distorts the effective CR transport velocities and the forces CR exerted on the thermal gas. Resolving these gradients is essential for understanding how CRs influence the formation, evolution, and morphology of cold gas clouds in the CGM.

As discussed in Sect. 6, the size of simulated clouds significantly influences the role of CRs in shaping internal cloud dynamics. The timescale for CR escape from a cloud depends quadratically on its radius, $t_{\text{cr}}\propto r_{\text{cloud}}^2$ (see Eq. (17)). Larger clouds confine CRs for extended periods, allowing their pressure forces to act on the cloud over longer timescales. This prolonged interaction amplifies the impact of CRs on the thermodynamic evolution of the cloud, potentially enhancing heating by CRs, altering cooling rates, and modifying condensation. Consequently, variations in cloud sizes can directly affect how CRs influence the stability, structure, and overall dynamical behaviour of cold CGM gas.

Figure 16 illustrates the impact of resolution on resulting cloud sizes. We show projections of the hydrogen column density, N_H, for simulations employing two-moment CR transport that are initialised with X_{cr, 0}=3. From left to right, the initial resolution of the simulations continuously increases. Higher resolution simulations show a higher number of more compact clouds or cloud complexes with visibly smaller associated cloud radii.

Fig. 16. Projection showing the hydrogen column density, NH, after the simulations run for 8 tcool. From left to right, we plot the results from simulations with increasing resolution where each simulation employs two-moment CR transport and is initialised with Xcr, 0=3. The sizes and number of the cold clouds significantly depend on resolution.

A possible mechanism for this trend is that CRs can escape more rapidly from small cloud structures in the high-resolution simulations while large clouds, present in the low-resolution simulations, confine CRs for longer times. We verify this expectation in Fig. 17, which displays the ratio of CR escape time to cloud collapse time, t_cr/t_collapse (left panel), and the relative CR pressure, X_cr (right panel) for clouds of varying radii. Each cloud is colour-coded based on the initial resolution of the corresponding simulation. We show 1D probability densities for the respective variables above each axis. The highest-resolution run produces the smallest clouds, which facilitates rapid CR escape. In these simulations, the majority of clouds exhibit escape times shorter than their collapse times. Conversely, at lower resolutions, clouds are larger, which slows CR escape because the escape times become longer than the collapse times. This trend is particularly evident in the right panel where we observe that the CR pressure support increases for larger clouds formed in lower-resolution simulations, which amounts to X_cr values ranging from approximately ten to 22. By contrast, the high-resolution runs produce smaller, more compact clumps, resulting in lower CR pressure support and values of X_cr ranging from approximately four to eight. This comparison emphasises the significant influence of resolution on the interplay between CR dynamics and cloud formation in the CGM.

Fig. 17. Same as Fig. 13 but for different initial resolutions. All simulations employ two-moment CR transport and are initialised with Xcr, 0=3. A lower resolution produces larger clouds, inhibiting rapid CR escape and leading to enhanced CR pressure support within cold clouds.

8. Discussion and limitations

8.1. Comparison to previous work

In the pioneering study by Butsky et al. (2020), the first comprehensive 3D MHD simulations investigating the role of CRs in the CGM were conducted using the Enzo code (Bryan et al. 2014). Their work explores the influence of CRs on the formation and evolution of cold gas driven by TI within the CGM. The simulation setup consists of a stratified box in hydrostatic equilibrium with an isocooling density profile. The initial magnetic field is oriented uniformly along the x-axis and has a magnetic-to-thermal pressure ratio of X_mag=0.01. To initiate gas condensation, density perturbations are manually introduced in the initial conditions. Our setup is similar to that described by Butsky et al. (2020) but with some revisions. We incorporate turbulent velocity and magnetic fields in the initial conditions, such that fluctuations in these two fields do not need to be seeded from the dynamics started by density fluctuations. These initial turbulent fields naturally introduce density perturbations, negating the need for manually imposed density variations and better capturing the interplay between turbulence, magnetic fields, and CR dynamics in shaping cold gas structures. The turbulent magnetic fields also point in random directions starting from the beginning of the simulations, which leads to anisotropic CR transport that is not artificially correlated with one of the coordinate axes of the simulation box or the symmetry axis of the density and pressure profiles.

Butsky et al. (2020) find that CRs provide non-thermal pressure support, preventing cold gas from compressing as it cools, which results in lower densities in comparison to the pure MHD case. This effect allows cold gas to cool almost isochorically when CR pressure is dominant. Additionally, the presence of CR pressure allows cold gas clouds to reach larger sizes. In cases where CR pressure is relatively low, the cold gas forms as a ‘mist’ of small, pressure-supported cloudlets, while cloud sizes can increase by orders of magnitude in CR-pressure-dominated halos with inefficient CR transport compared to purely thermally driven models. Our work supports these findings for the limiting case of purely advective CR transport. By incorporating the streaming and diffusion CR transport mechanisms, Butsky et al. (2020) demonstrate that CRs effectively redistribute pressure from dense cold clouds to the surrounding hot medium, leading to reduced cloud sizes and increased cold gas density and mass flux relative to adiabatic CR models. However, they also conclude that actively transported CRs have a substantial influence on cloud dynamics – a finding that is not supported by our results, which suggests that cloud formation and evolution is governed by TI while CRs have at best a minor impact.

In Sect. 7, we discussed that insufficient resolution significantly affects the cloud-CRs interaction. This is because cloud radii are naturally limited by the numerical resolution, leading to larger CR escape times, which allows them to counteract TI-induced collapse for longer times. In their setup, Butsky et al. (2020) employ a static grid with a spatial resolution of 0.685 kpc per computational cell. Estimating that the smallest cold clouds span roughly three cells in diameter, the minimum cloud radius in their simulations is approximately r_cloud∼1 kpc. With their maximum CR diffusion coefficient of κ_cr=7.9×10²⁹ cm² s⁻¹, the minimum CR escape time is t_cr∼4×10⁻⁴ Gyr. Using their cooling function and the temperature probability density function for the cold, dense phase at T∼10^4.75 K, we estimate the cooling timescale to be t_cool∼2×10⁻⁴ Gyr. Based on these parameters, we find that the ratio of the CR escape time to the collapse time is t_cr/t_cool∼2, indicating that CRs cannot escape collapsing regions quickly enough. This delayed escape maximises the non-thermal pressure support of CRs on the clouds. Even in their highest-resolution simulations, this ratio can be reduced by a maximum factor of two and is thus still above the parameter regime where CRs could escape faster than clouds are collapsing due to TI. From these considerations, we conclude that the limited cloud-growth observed by Butsky et al. (2020) can be explained by a CR-induced delay of cloud collapse, which is artificially promoted by insufficient resolution.

Tsung et al. (2023) investigate the impact of CR heating on TI in a CGM-like gravitationally stratified medium with vertical magnetic fields and two-moment CR transport utilising the approach of Jiang & Oh (2018) in various 2D simulations with the Athena++ code (Stone et al. 2020). In the linear phase, they confirm that CR heating causes entropy modes to propagate at a velocity proportional to the Alfvén speed, as predicted by theory (Kempski & Quataert 2020). However, this propagation, which could theoretically suppress TI, is limited to a narrow range of conditions where the ratio of gas cooling timescale and CR heating timescale approaches unity. The study's main findings lie in the non-linear phase where TI leads to significant mass condensation and a cold disc formation in the mid-plane, altering CGM properties and phase structure.

8.2. Limitations

Our simulations are intentionally highly idealised, which is advantageous for isolating the specific influence of CRs on TI and their interactions with cold clouds. However, this approach lacks a fully realistic physical context that could affect cold gas formation and survival in the CGM. The heating model we use is a simplified representation of an isolated, dynamically stable, long-lived CGM, where we balance total cooling and heating globally on average. In a more physically realistic setup, local processes such as gas accretion, star formation, stellar and AGN feedback, and galactic mergers would naturally contribute to this equilibrium. This is particularly relevant in simulations evolved over several hundreds of Myr. On these timescales, the assumption of an isolated CGM becomes increasingly invalid as a consequence of the above processes (Tumlinson et al. 2017; Hummels et al. 2019; Nelson et al. 2020).

The initial density and temperature profiles in our simulations follow an isocooling relation ensuring a uniform initial cooling time across the computational domain. While this approach results in profiles that closely resemble isothermal ones, we do not intend to replicate realistic halo configurations with this setup. Instead, the primary motivation for using these profiles is to maintain consistency and comparability with previous studies that employed similar initial conditions and to allow for a more direct comparison of results. To maintain the global stability of the initial atmosphere while including additional CR pressure, the total mass in the system increases.

In our simulations, we adopt a fiducial ratio of τ=t_cool/t_ff=0.3, representative of regions in the CGM where rapid cooling promotes the condensation of cold gas. However, in reality, the CGM exhibits a wide range of thermodynamic conditions, with local cooling times varying significantly due to differences in density, thermal pressure, and metallicity. Even for fast CR transport models, there are likely to be dense, rapidly cooling regions of the CGM where t_cr/t_cool≳1, leading to effective CR trapping. Thus, while our simulations focus on a specific cooling regime that favours the formation of multiphase gas, we expect a broad spectrum of CR–gas coupling efficiencies across the CGM, depending sensitively on the local cooling environment and corresponding collapse timescales.

The CR transport velocity in our simulation is governed by interactions between CRs and gyroresonant Alfvén waves alongside the processes that drive the growth and damping of these waves. Specifically, we consider the gyroresonant interactions of CRs and the non-linear Landau damping of Alfvén waves (Kulsrud 2005) because these are generally regarded as the primary influences on CR transport. However, additional damping mechanisms including ion-neutral damping (Kulsrud & Pearce 1969), linear Landau damping (Wiener et al. 2018), and turbulent damping (Farmer & Goldreich 2004; Lazarian 2016; Holguin et al. 2019) may act in realistic systems to limit Alfvén wave amplification to reduce the CR scattering rate. This reduction, however, results in less confined CR populations with effectively higher CR transport velocities and positions the results of this paper as a lower bound on the CR escape timescale from collapsing regions and an upper bound on their pressure support in the cold phase.

Our simulations are incapable of fully resolving the cooling length, l_cool=c_s t_cool as demonstrated in Appendix A. While this limitation introduces some uncertainty in capturing the smallest scales of TI and cloud fragmentation (McCourt et al. 2018), we believe our key findings remain robust, particularly regarding the CR escape timescale from cooling clouds. Resolving l_cool is important for accurately capturing small-scale structures in the cold gas. However, our results consistently show that the dominant processes governing CR transport and their interaction with cooling clouds are largely determined by larger-scale dynamics and the interplay between the CR escape timescale and the cloud collapse timescale. This suggests that our conclusions regarding CR escape and its effects on cold gas dynamics are not strongly dependent on resolving l_cool.

In addition, our simulations do not fully resolve the small-scale structure of the magnetic field within cold clouds and their turbulent boundary layers. Capturing these scales could impact CR transport. The turbulent boundary layers between the cold clouds and the surrounding hot, diffuse medium generate complex density and pressure gradients, which enhance CR scattering and thereby reduce their effective transport speed (Yan & Lazarian 2004). These boundary layers can therefore act as bottlenecks for CRs (Wiener et al. 2017b), limiting their escape from the cloud interior and causing them to accumulate near the cloud edges. This accumulation increases the local CR pressure, which in turn exerts an outward force on the cloud outskirts. Moreover, magnetic fields within these turbulent layers are amplified and become more organised within cold gas, often forming filamentary structures (Zhao & Bai 2023; Das & Gronke 2024). While such local coherence may promote CR streaming along field lines, the overall magnetic topology across the multiphase medium remains complex. The suppression of turbulent mixing and increased magnetic pressure reduce the connectivity between different gas phases, potentially confining CRs within cold regions and impeding their large-scale transport.

9. Conclusions

In this paper, we have explored the impact of CRs on the formation and evolution of cold clouds in the CGM. Utilising the moving-mesh code Arepo, we performed 3D CRMHD simulations of idealised CGM environments. These simulations started from initial conditions characterised by an isocooling gas stratification, a turbulent velocity field, a turbulent magnetic field, a realistic gas cooling function, and a global heating function to replicate heating from galactic feedback. We incorporated CRs as a constant fraction of the thermal pressure, investigated various CR transport mechanisms, and assessed their significance. Our key findings are outlined as follows:

• 1.

  Purely advective CRs significantly alter the morphology of the CGM. The additional non-thermal pressure support from CRs either inhibits the collapse of thermally unstable regions or delays their collapse in comparison to purely thermal scenarios (cf. Fig. 4). The strength of this effect is proportional to the initial CR pressure in the simulation. However, purely advective CRs represent an academic scenario because realistic CRs are transported down their pressure gradient along magnetic field lines via streaming and diffusion in realistic astrophysical environments.

• 2.

  Including CR streaming and diffusion fundamentally alters the evolution of a collapsing cold cloud. As CRs are no longer fully coupled to the thermal gas movement, their dynamics are predominantly dictated by the CR pressure gradient along the magnetic field lines. When a cold cloud collapses, not only is its gas compressed, but the CRs that are tied to the magnetic field lines are also compressed. This increases the CR pressure within the collapsing regions and creates a pressure gradient between the cloud interior and the surrounding medium. The magnetic field lines that penetrate the cloud serve as pathways for the CRs to escape the collapsing regions on very short timescales (cf. Fig. 14). This counteracts the build-up of CRs inside the cloud and impedes the stabilising effect of the CR pressure. Consequently, CRs cannot prevent the cloud's collapse, and the morphology of the CGM nearly resembles that of a purely thermal setup. This conclusion is mostly independent of the initial CR pressure (cf. Fig. 7).

• 3.

  The effective CR diffusion coefficient, κ_eff, exhibits a wide dynamic range, spanning several orders of magnitude from 10²⁵ to 10³⁴ cm² s⁻¹ (see Fig. 10), and it is generally dependent on X_cr. Given the significant variation in X_cr expected in realistic CGMs, this broad range underscores the importance of employing the two-moment CR transport method. This approach is able to capture the anisotropic and dynamic nature of CR propagation in these diverse conditions.

• 4.

  The CR escape timescale from collapsing regions is dictated by the effective CR diffusion coefficient, κ_eff (cf. Fig. 11), rather than by the CR pressure. When κ_eff is close to or smaller than its canonical value in the ISM of the Galaxy (≲3×10²⁸ cm² s⁻¹), CRs can sustain their pressure support inside the thermally unstable regions for a longer period, thereby delaying the onset of collapse. For higher values of κ_eff (≳3×10²⁹), CRs escape quickly from the cold clouds (cf. Fig. 13).

• 5.

  Resolution plays a pivotal role in accurately assessing the influence of CRs on cold gas in the CGM. The size of condensing clouds is strongly influenced by the resolution of computational cells, with low-resolution simulations yielding clouds up to ten times larger than those formed in high-resolution runs (cf. Fig. 16). Since the CR escape timescale depends quadratically on cloud size, $t_{\text{cr}}=r_{\text{cloud}}^2/{\kappa }_{\text{eff}}$, larger clouds retain CRs for significantly longer durations. This extended confinement amplifies the CR impact on cloud dynamics (cf. Fig. 17).

These findings collectively demonstrate that the dynamical role of CRs in the CGM is highly sensitive to their transport properties and to the numerical resolution of the simulations. Accurate modelling of CR transport, particularly through two-moment methods, is therefore essential for capturing their true impact on cold gas evolution and the resulting CGM morphology.

Acknowledgments

We thank the anonymous referee for helpful comments that improved the clarity of this manuscript. We acknowledge support by the European Research Council under ERC-AdG grant PICOGAL-101019746. This work was supported by the North-German Supercomputing Alliance (HLRN) under project bbp00070. The projections and slices presented in this work were generated using Paicos (Berlok et al. 2024). Additionally, all Python-based analysis scripts relied on matplotlib (Hunter 2007), seaborn (Waskom 2021) and cmocean (Thyng et al. 2016) for visualisation including colour maps, and numpy (Harris et al. 2020) for numerical computations.

References

Appendix A: Resolution of cooling length

To highlight the impact of numerical resolution in our fiducial simulation with τ=0.3 and X_{cr, 0}=3, we show the relationship between the cooling length, l_cloud=c_s t_cool, and the diameter of the computational cell, d_cell=(6V_cell/π)^1/3, in Fig. A.1. It shows a 2D histogram of the cooling length versus the cell size. Each bin is colour-coded by its average temperature, providing a detailed view of how resolution varies with the thermodynamic state of the gas. The green dashed line marks the thermodynamic states where the corresponding cooling length is equal to the cell size, delineating the boundary where thermal processes are adequately resolved. Our results show that gas with temperatures T≳10⁵ K lies to the right of this line, indicating that our simulations properly resolve the cooling length for this temperature regime. Conversely, gas with T≲10⁵ K falls to the left of this line, suggesting that cooling processes in this regime may be under-resolved.

Fig. A.1. Two-dimensional histogram displaying the cooling length across various computational cell sizes of our fiducial simulation with τ=0.3 and Xcr, 0=3. Each bin is colour-coded according to its average temperature. The green dashed line indicates where the cooling length equals the cell size, demonstrating that gas with T≳105 K is well-resolved in our simulations.

Appendix B: Thermal instability without cosmic rays

In this appendix, we examine the onset of TI in the absence of CR physics. In their seminal work, McCourt et al. (2012) highlighted the critical importance of the cooling-to-free-fall time ratio, τ=t_cool/t_ff. They found that gas condensation occurs when τ≲1. Since t_cool is governed by the cooling function, changes in τ are achieved by modifying $t_{\text{ff}}=\sqrt{2z/g(z)}$, where z represents the vertical distance from the simulation box centre, and g(z) denotes the gravitational acceleration, as defined in Eq. (5). Specifically, we adjust the scale height, H, to increase or decrease t_ff, thereby altering τ.

Figure B.1 shows hydrogen column density projections, N_H, from simulations with identical initial density and temperature profiles but varying τ. In simulations with smaller τ, cold gas condenses rapidly from the background medium on timescales much shorter than t_ff. As a result, a significant amount of cold gas remains suspended at high altitudes, even after multiple cooling cycles. Conversely, as τ increases, gravitational acceleration becomes the dominant process. Under these conditions, condensed gas clumps increasingly precipitate toward the gravitational centre, where they accumulate. In the case with τ=10, we do not observe any TI. Notably, in contrast to Butsky et al. (2020), we do not find an increase in cloud sizes with increasing τ. We attribute this discrepancy to the significantly higher resolution of our simulations, which prevents artificial cloud coagulation caused by low-resolution-induced pressure gradients (Gronke & Oh 2020; Fielding et al. 2020). Overall, our findings are consistent with previous studies by McCourt et al. (2012), Ji et al. (2018), and Butsky et al. (2020).

Fig. B.1. Projections showing the hydrogen column density, NH, after the simulations run for 6 tcool without CRs. From left to right, we plot the results of the simulations employing different ratios of the cooling-to-free-fall time, τ. From the figure, it is evident that this ratio plays a crucial role in determining the onset of TI.

Figure B.2 presents the cold gas metrics used in this study: density fluctuations, cold mass and volume fractions, and cold mass flux (from left to right). After the initial turbulence decays, density perturbations begin to grow in simulations with τ≲3 at around t∼2 t_cool. This leads to the condensation of cold gas from the hot ambient medium, increasing both the cold mass and the volume it occupies. For runs with τ≳1, these values decrease at later times due to cold gas precipitation out of the analysed domain. The simulation with τ=3 produces significantly less cold gas because the shorter free-fall timescale allows cooling gas parcels to descend into higher-pressure regions near the gravitational centre earlier, leading to enhanced compressive heating of the cooling gas. Meanwhile, the cold mass flux increases with higher τ due to the faster velocities of cold gas driven by the shorter t_ff.

Fig. B.2. Time evolution of cold gas metrics for simulations without CRs. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. The different colours correspond to different values of τ. We evaluate each computational cell in the region 0.25 H≤|z|≤2.75 H. The dashed vertical line denotes the time of the analyses shown in Fig. B.1.

All Tables

All Figures

Fig. 1. Simulation domain overview. The figure shows a slice through the x–z plane at y=0 of our tall-box setup. The box extends from −3H to +3H in the z-direction and from −0.5H to +0.5H in the x- and y-directions, where H=30.11 kpc. Periodic boundary conditions are applied in the x- and y-directions, while the z-direction features an outflow boundary condition. The simulation is conducted without CRs and assumes a target mass of mtarget=22 M⊙. The computational domain is initialised with 256×256×1536 equally spaced mesh-generating points. The left elongated panel illustrates the resulting mass density, and the right elongated panel displays the temperature. The five smaller panels provide magnified views of various structures of different sizes, shapes, densities, and temperatures formed due to TI.

In the text

	Fig. 2. Projection in x–z of the hydrogen column density, NH. Red circles represent the radii of clouds detected by the cloud-identification algorithm and red dots indicate their respective centroids. We note that overlapping circles result from the projection of clouds located at different y positions and do not necessarily represent physical overlaps in three-dimensional space.
In the text

Fig. 3. Initial density profiles of the simulation box for various Xcr, 0. The solid lines show the implemented isocooling density profile, the dashed lines show the corresponding isothermal density profile with T0=106 K for comparison. For a higher initial Xcr, 0, more mass is contained in the atmosphere because the additional CR pressure supports the HSE without contributing to the gravitational force.

In the text

Fig. 4. Projections of the hydrogen column density, NH, employing purely advective CR transport. The panels show the x−z plane of a section extending from z=0 to z=2H and the depth of the projection is 1 H. The snapshot is at t=7 tcool≈2 tff. All simulations were initialised with different relative CR pressure, Xcr, 0. The density of the clouds decreases with increasing CR pressure, while the cloud size increases. In the CR pressure dominated run with Xcr, 0=3, cold gas condensation is entirely suppressed. This illustrates that advective CRs alter the morphology of the CGM significantly and suppress the formation of cold clouds.

In the text

	Fig. 5. Time evolution of cold gas metrics for simulations employing advective CR transport with varying Xcr, 0. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. We evaluated each computational cell in the region 0.25 H≤|z|≤2.75 H. The dashed vertical line marks the simulation time corresponding to the results shown in Fig. 4.
In the text

	Fig. 6. Gallery showing 2D slices of various quantities from the simulations employing two-moment CR transport with an initial CR pressure fraction of Xcr, 0=3. All snapshots are at t=7 tcool. The top row illustrates quantities of the thermal gas, while the bottom row presents variables related to CR transport.
In the text

	Fig. 7. Slices showing the mass density in the x−z plane centred at y=0 of simulations employing two-moment CR transport and different relative CR pressures, Xcr, 0. All snapshots are at t=8 tcool. The initial CR pressure barely influences the morphology of the cold gas in terms of cloud size or density.
In the text

Fig. 8. Time evolution of cold gas metrics for simulations employing two-moment CR transport with varying Xcr, 0. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. We evaluate each computational cell in the region 0.25 H≤|z|≤2.75 H. All cold gas metrics show a very similar evolution and are almost independent of the initial CR pressure.

In the text

Fig. 9. Intrinsic CR diffusion coefficient, κcr, for simulations utilising two-moment CR transport with varying initial values of Xcr, 0. Scatter points represent the distribution of CRs as a function of κcr and Xcr at t=6 tcool, with darker colours indicating higher CR-energy-weighted probability density in a logarithmic colour scheme. The corresponding 1D probability density of κcr is displayed with a linear scaling on the right-hand side. In the regime of efficient CR scattering (small κcr), larger Xcr, 0 values increase the overall CR energy density in the CGM allowing for more CR energy to be converted into gyroresonant Alfvén waves and larger CR scattering rates, thereby reducing κcr (see Eq. (13)). In regions of a smooth CR distribution, the driving of Alfvén waves is inefficient in comparison to wave damping processes so that CRs are poorly coupled to the plasma, leading to a very large value of κcr.

In the text

	Fig. 10. Same as Fig. 9 but for the effective CR diffusion coefficient, κeff.
In the text

	Fig. 11. Projections of the hydrogen column density, NH, for simulations utilising different CR transport models. All runs were initialised with Xcr, 0=3, and the snapshots are at t=7 tcool. The figure demonstrates that the onset of collapse is governed by the CR transport velocity rather than the CR pressure.
In the text

	Fig. 12. Time evolution of cold gas metrics for simulations with different CR transport models and Xcr, 0=3. Metrics are evaluated for computational cells within 0.25 H≤|z|≤2.75 H. The dashed vertical line marks the time corresponding to the snapshot shown in Fig. 11. Lower CR transport velocities delay the onset of cloud collapse and reduce the amount of condensed cold gas.
In the text

Fig. 13. Ratios of the CR transport timescale to the cloud collapse timescale (left) and the relative CR pressure within cold clouds (right), depicted for various cloud radii. The 1D probability density of the corresponding variable is displayed above each axis with a linear scaling. We analyse cold gas clouds within the region −0.25H≤x, y≤0.25H and 0.25H≤z≤1.75H (We focus on this reduced region to avoid double counting clouds across periodic boundaries and to ensure consistent analysis across resolutions (see Sect. 7) because the memory-intensive cloud finder makes analysing the full domain computationally very expensive at high resolutions.), considering eight consecutive snapshots following the identification of the first cloud by the cloud finder. Different colours represent different CR transport models (as detailed in the figure legend). All simulations are initialised with Xcr, 0=3. The dashed vertical line in the left panel separates regions where clouds collapse is slower (tcr/tcollaspe<1) or faster (tcr/tcollaspe>1) than CRs escape these clouds. Fast transport enables CRs to escape collapsing clouds on timescales shorter than the cloud collapse time, reducing their relative pressure. In contrast, slow transport traps CRs within collapsing regions, increasing their pressure support.

In the text

Fig. 14. Visual impression of the magnetic field topology in and around one of the clouds. This cloud has formed through TI and is currently falling down. The volume rendering highlights dense gas and the magnetic field lines are coloured according to the local magnetic field strength. The shown magnetic field lines connect the inside of the cloud to the ambient gas above the cloud. The magnetic field inside the cloud is stronger (achieving B∼1 μG) compared to the magnetic field in the ambient gas.

In the text

Fig. 15. Mass-weighted nH–T diagrams from our simulation suite. Columns represent different CR transport methods, while rows correspond to varying initial CR pressure fractions, Xcr, 0. The 1D probability densities of the respective variables are shown above each axis: columns reflect the resulting nH distributions for a fixed CR transport method across different Xcr, 0 and rows show the resulting T distributions for a fixed Xcr, 0 across different CR transport methods. The phase diagrams are constructed using 70 bins of equal (logarithmic) width for both variables. All simulations are analysed at t=6 tcool.

In the text

	Fig. 16. Projection showing the hydrogen column density, NH, after the simulations run for 8 tcool. From left to right, we plot the results from simulations with increasing resolution where each simulation employs two-moment CR transport and is initialised with Xcr, 0=3. The sizes and number of the cold clouds significantly depend on resolution.
In the text

	Fig. 17. Same as Fig. 13 but for different initial resolutions. All simulations employ two-moment CR transport and are initialised with Xcr, 0=3. A lower resolution produces larger clouds, inhibiting rapid CR escape and leading to enhanced CR pressure support within cold clouds.
In the text

	Fig. A.1. Two-dimensional histogram displaying the cooling length across various computational cell sizes of our fiducial simulation with τ=0.3 and Xcr, 0=3. Each bin is colour-coded according to its average temperature. The green dashed line indicates where the cooling length equals the cell size, demonstrating that gas with T≳105 K is well-resolved in our simulations.
In the text

	Fig. B.1. Projections showing the hydrogen column density, NH, after the simulations run for 6 tcool without CRs. From left to right, we plot the results of the simulations employing different ratios of the cooling-to-free-fall time, τ. From the figure, it is evident that this ratio plays a crucial role in determining the onset of TI.
In the text

	Fig. B.2. Time evolution of cold gas metrics for simulations without CRs. From left to right, we show density fluctuations, cold mass fraction, cold volume fraction, and cold mass flux. The different colours correspond to different values of τ. We evaluate each computational cell in the region 0.25 H≤|z|≤2.75 H. The dashed vertical line denotes the time of the analyses shown in Fig. B.1.
In the text