© The Authors 2025

1. Introduction

Cosmological simulations have become powerful tools for understanding the physics of cosmic structure formation. Some of these numerical calculations only include gravity, thereby simulating the evolution of mock universes that are solely comprised of (collisionless) dark matter particles. However, recent simulations now commonly include baryons and many of the associated physical processes that are considered important for galaxy formation and evolution (for a recent overview of cosmological simulations, see Vogelsberger et al. 2020).

For instance, it is common practice to include models for the radiative cooling of gas, either with a simple primordial network (e.g., Katz et al. 1996) or incorporating additional processes such as metal line cooling (e.g., Wiersma et al. 2009). This can be paired with a model for a time-dependent UV background (e.g Faucher-Giguère et al. 2009; Haardt & Madau 2012) with corrections for self-shielding (e.g., Rahmati et al. 2013). As gas reaches sufficiently high densities, star formation can occur, creating stellar particles that can represent individual stars (e.g., Gutcke et al. 2021; Smith et al. 2021) or entire populations (e.g., Springel & Hernquist 2003).

Galactic feedback processes are an important mechanism to regulate the star formation rate (SFR) and the amount of stellar mass growth. In low-mass galaxies (M_⋆ ≲ 10^10.5 M_⊙), outflows driven by stellar and supernova feedback play an important role (e.g., Fielding et al. 2017). In addition to keeping the SFR in check, they also inject significant amounts of energy and metals into the surrounding halo gas (e.g., Nelson et al. 2019a; Pandya et al. 2021). Higher-mass galaxies are instead dominated by the activity of the central supermassive black hole (SMBH), which may release energy in thermal, kinetic, or other physical channels (e.g., Schaye et al. 2015; Weinberger et al. 2017; Davé et al. 2019). These feedback episodes launch gas from the galaxies at high velocities (e.g., Oppenheimer et al. 2020; Ramesh et al. 2023a), alter the morphology and kinematics of the circumgalactic medium (CGM; Nelson et al. 2015; Kauffmann et al. 2019), increase the cooling time of gas in the surrounding halo (Zinger et al. 2020), and hinder future gas accretion and star formation (e.g., Davies et al. 2020).

In addition, recent simulations now start to include important physics arising from nonthermal processes, most notably, magnetic fields (e.g., Pakmor et al. 2017; Marinacci et al. 2018) and their impact on the propagation of galactic outflows (Steinwandel et al. 2020) as well as halo-scale gas properties (Pakmor et al. 2020). In addition, theory suggests that magnetic fields may play an important role in the evolution of small-scale CGM clouds, either by suppressing fluid instabilities (Cottle et al. 2020; Sparre et al. 2020; Ramesh et al. 2025), which provides additional pressure support (Nelson et al. 2020; Ramesh et al. 2023b; Fielding et al. 2023), or, in the presence of a draped layer (Dursi & Pfrommer 2008; Pfrommer & Dursi 2010), by boosting the drag force they experience (McCourt et al. 2015; Ramesh et al. 2024). These cool clouds may play a part in replenishing the galactic cold gas supply (Lepine & Duvert 1994). This unavoidably links their survival and evolution to the star formation activity of galaxies.

The second nonthermal physical component of significant interest for galaxy evolution is cosmic rays (CRs). The energy density of CRs in galaxies is non-negligible and in rough equipartition with magnetic, thermal, and turbulent energies (Cox 2005; Naab & Ostriker 2017). Depending on their momenta, CRs can influence gas in various ways. Cosmic-ray particles in the ≲100 MeV c⁻¹ range of the spectrum can ionize their surrounding gas (e.g., Ivlev et al. 2018; Padovani et al. 2020), while CRs with ∼GeV c⁻¹ can drive powerful galaxy-scale outflows (e.g., Zirakashvili et al. 1996; Girichidis et al. 2016; Chan et al. 2022) through gradients in their pressure structure (Pfrommer et al. 2017; Butsky et al. 2020). Much like magnetic fields, this added nonthermal pressure support may also provide additional stability to small-scale cold gas structures (Brüggen & Scannapieco 2020; Butsky et al. 2022).

More energetic CRs (≳TeV c⁻¹) are not dynamically relevant because they are subdominant in terms of their contribution to the total energy budget. However, pion production via hadronic interactions leads to a subsequent decay into γ-rays, which provides an important and observable nonthermal radiative signature (Kotera & Olinto 2011; Werhahn et al. 2021). While these high-energy CRs are less impacted by streaming losses (Girichidis et al. 2024), their lower-energy counterparts excite Alfvén waves through the CR streaming instability when they propagate faster than the Alfvén velocity (Kulsrud & Pearce 1969). This effectively heats the gas because these waves damp on short timescales (Wiener et al. 2017; Ruszkowski et al. 2017; Thomas et al. 2023).

Cosmic rays may thus play an important role in the evolution of galaxies and their gaseous halos. They may affect the accretion of gas onto galaxies by altering the flow structure of CGM gas, thereby modifying structural properties such as the sizes of galactic disks (Buck et al. 2020) and time-integrated quantities such as the star formation rates within the ISM (Jubelgas et al. 2008; Farcy et al. 2022). Outflows driven by CRs can be smoother, cooler, and denser than those driven by supernovae (Girichidis et al. 2018; Armillotta et al. 2024), which in turn may affect the escape of Lyα photons from galaxies (Gronke et al. 2018). Their added nonthermal pressure may result in halos that are significantly cooler (Ji et al. 2020). Several studies also suggested that CRs may play a role in magnetizing gas through the so-called galactic dynamo (Parker 1992; Hanasz et al. 2009; Pfrommer et al. 2022).

The impact of CRs has been studied so far in isolated (idealized) galaxy setups and in cosmological zoom-in simulations of single halos. Their impact on a cosmologically representative galaxy population remains unknown, however: They have never been studied in a large-volume cosmological (magneto)hydrodynamical simulation. This is partly due to the small timesteps that are required to track their evolution (Hopkins 2017; Chan et al. 2019), as well as to the complexity involved in representing CRs over a broad spectral distribution (e.g., Girichidis et al. 2020; Hopkins 2023). This makes it challenging to include this component in already computationally expensive simulations (see also Jiang & Oh 2018; Gupta et al. 2021).

To overcome this limitation, Hopkins et al. (2023) recently proposed a subgrid toy recipe aimed at capturing the leading first-order impact of CRs on galaxy growth. The model is parameterized by two simple quantities, and the implementation is similar to the gravity-tree algorithm (Barnes & Hut 1986). It allows the user to estimate the CR contribution at any point in space with negligible computational effort.

We implemented this simple CR scheme into the AREPO moving-mesh code (Springel 2010) and coupled it to the IllustrisTNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018). We then ran a suite of (25 Mpc h⁻¹)³ cosmological magnetohydrodynamical boxes at the resolution of TNG100, and we studied the impact of CRs on the galaxy and halo properties. The paper is organized as follows: In Section 2 we provide an overview of the simulation suite and methods. The results are presented in Section 3, and they are summarized in Section 4.

2. Methods

2.1. Initial conditions and overview of the setup

The results we presented were derived from a suite of (25 Mpc h⁻¹)³ ∼ (37 Mpc)³ cosmological magnetohydrodynamical boxes run with the moving mesh code AREPO (Springel 2010). The initial conditions, generated using N-GenIC (Springel et al. 2005), are identical to the large number of TNG variation boxes presented in Pillepich et al. (2018). In particular, we use the L25n512 box containing 512³ dark matter particles and 512³ initial gas cells, giving a dark matter (average baryonic) mass of ∼10⁷ (10⁶) M_⊙. This corresponds to a spatial resolution of ∼350 pc (the average star-forming gas cell size) to ∼740 pc (the gravitational softening of stars and dark matter), equivalent to the TNG100-1 simulation (Nelson et al. 2019b).

Due to the relatively small size of the simulation domain, massive objects (M_200c ≫ 10¹² M_⊙) are limited in number. At the other end of the mass spectrum (M_200c ≪ 10^10.5 M_⊙), structures are typically sampled by ≲O(1000) baryonic elements. In order to maximize statistics and resolution, we focus our main analysis on halos in the mass range ∼10¹¹–10¹² M_⊙, identified via the friends-of-friends algorithm (Davis et al. 1985). In addition, unless otherwise stated, we restrict the analysis to central galaxies of these halos, as identified by SUBFIND (Springel et al. 2001). Throughout this work, the term halo mass refers to M_200c, and (galaxy) stellar mass to the cumulative mass of all stars within twice the stellar half mass radius.

We adopt a cosmology consistent with the Planck 2016 analysis (Planck Collaboration XIII 2016): Ω_Λ = 0.6911, Ω_m = 0.3089, Ω_b = 0.0486 and h = 0.6764 [100 km s⁻¹ Mpc⁻¹].

2.2. Galaxy formation model

We used the IllustrisTNG model (hereafter, TNG model) for the physics of galaxy formation and evolution (Weinberger et al. 2017; Pillepich et al. 2018). TNG includes models for primordial and metal line cooling in the presence of a time-dependent metagalactic radiation field, formation of stars and supermassive black holes (SMBH), stellar evolution and enrichment, feedback processes driven by both supernovae (SNe) and SMBHs, and other related processes. The TNG model also includes ideal magnetohydrodynamics (Pakmor et al. 2011, 2014) with an initial (z ∼ 127) primordial seed field of ∼10⁻¹⁴ cG, while the Powell et al. 1999 eight-wave divergence cleaning scheme is used to maintain ∇ ⋅ B ∼ 0.

We froze every aspect of the TNG model but one: Stellar winds, that is, the deposition of mass, energy and momentum transfer from SNe feedback (Springel & Hernquist 2003). These are given 90% of their original (kinetic) energy, while the remaining 10% is assumed to accelerate CRs (e.g., Helder et al. 2012; Ackermann et al. 2013). At all other positions in space, we estimate the energy density of CRs (e_cr) using the Hopkins et al. (2023) subgrid model. This adopts a number of approximations, including but not limited to: assuming spherical symmetry around point-like sources, constant streaming/diffusion coefficients in space and time, and that the CR energy equation is in a steady-state. Ultimately, e_cr at the location of a gas cell in the simulation i is estimated as the sum over all (j) star forming gas cells: $e_{\mathrm{cr},\mathrm{i}}=Q_i^{\mathrm{atn}}{\displaystyle \sum _j}{\dot{E}}_{\mathrm{cr},\mathrm{j}}^{\mathrm{atn}}F(r_{\mathit{ij}})$, where $Q_i^{\mathrm{atn}}=e^{-\mathrm{\Delta }{\tau }_{\mathrm{cr},\mathrm{i}}}$ accounts for the attenuation of CRs. Here,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\tau }_{\mathrm{cr},\mathrm{i}}=\frac{{\psi }_{\mathrm{loss}}^i}{2}{[\mathrm{\Delta }{x}_i^2+{\left(\frac{{\rho }_{\mathrm{gas},\mathrm{i}}}{|\mathrm{\nabla }{\rho }_{\mathrm{gas},\mathrm{i}}|}\right)}^2]}^{1/2}\end{array}$$(1)

incorporates the local gas properties. ψ_loss ∼ 10⁻¹⁶ (6.4 n_n + 3.1 n_e + 1.8 n_HI)¹ cm³ s⁻¹ represents a cumulative loss term for various dissipation processes (Mannheim & Schlickeiser 1994; Guo & Oh 2008), while Δx, ρ and ∇ρ correspond to the size, density and density gradient of the corresponding gas cell.

The energy input rate is ${\dot{E}}_{\mathrm{cr},\mathrm{j}}^{\mathrm{atn}}=Q_j^{\mathrm{atn}}{⟨{\dot{E}}_{\mathrm{cr}}⟩}_j$, where the (attenuated) contribution by a star-forming gas cell j is given by ${⟨{\dot{E}}_{\mathrm{cr}}⟩}_j={ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}{E}^{\mathrm{SNe}}{\dot{M}}_j$. Here, ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ is the fraction of SNe energy transferred to CRs, which we fix at 10%. An energy input of E^SNe is available for every unit stellar mass formed, which for the initial mass function and SNe progenitor mass we employ is ∼ 5.2 × 10⁴⁹ erg M_⊙⁻¹, and Ṁ is the star formation rate. In effect, this uses an approach of the type realized by the locally extincted background radiation in optically thin networks (LEBRON) (Hopkins et al. 2018), where the attenuation over the entire trajectory from the source to the target is approximated by the local terms at these locations.

Last, $F(r_{\mathit{ij}})={[4\pi r_{\mathit{ij}}({\kappa }_{\mathrm{eff}}+v_{\mathrm{st},\mathrm{eff}}{r}_{\mathit{ij}})]}^{-1}{e}^{-r_{\mathit{ij}}^2/r_{\max }^2}$ modulates the e_cr contribution based on the distance (r_ij) between the source and the target. κ_eff and v_st, eff are the two free parameters of the model that correspond to effective diffusion- and streaming-like terms. The scale r_max = $(t_{\max }{v}_{\mathrm{st},\mathrm{eff}}/2)\{1+[1+16{\kappa }_{\mathrm{eff}}/{(v_{\mathrm{st},\mathrm{eff}}^2{t}_{\max )}]}^{1/2}\}$ is a proxy for the maximum distance that CRs may possibly propagate within a time interval t_max, which at each timestep we set to the time since the start of the simulation (see Hopkins et al. 2023 for a more detailed discussion).

At every timestep, e_cr was computed for all active gas cells using a tree-walk similar to the tree-walk that is frequently used for gravity (Barnes & Hut 1986). Following the calculation, the corresponding (nonthermal) pressure contribution [P_cr = (γ_cr − 1)e_cr = e_cr/3] is added under the ultra-relativistic approximation (γ_cr = 4/3). In practice, this not only involves updating the total pressure as P_tot = P_non, cr + P_cr while solving the Riemann problem, but also modifying the velocity of the fast magneto-acoustic wave that is required by the HLLD solver (Miyoshi & Kusano 2005), as

$$\begin{array}{c}\hfill {\mathrm{c}}_{\mathrm{f}}={\left[\frac{{\gamma }_{\mathrm{eff}}\mathrm{P}+{|\mathbf{B}|}^2+{[({\gamma }_{\mathrm{eff}}\mathrm{P}+{|\mathbf{B}|}^2{)}^2-4{\gamma }_{\mathrm{eff}}\mathrm{P}{\mathrm{B}}_{\mathrm{n}}^2]}^{1/2}}{2\rho }\right]}^{1/2},\end{array}$$(2)

where γ_eff P = γ_th P_th + γ_cr P_cr, ρ is the density of gas, B the magnetic field vector, and B_n the component of the field perpendicular to the face of the cell. In addition, we include a heating term [${\dot{e}}_{\mathrm{th},\mathrm{gas},\mathrm{cr}}$ = e_cr (0.9 n_n + 1.6 n_e)×10⁻¹⁶ cm³ s⁻¹] that accounts for Coulombic and hadronic losses (Mannheim & Schlickeiser 1994; Guo & Oh 2008), and an ionization rate [Γ_cr = 7.5 × 10⁻¹⁸ s⁻¹ (e_cr/eV cm⁻³)] due to low-energy nonrelativistic CRs (Indriolo et al. 2015). Effectively, we follow a similar numerical implementation of the coupling of CRs to gas as other explicit CR modules in AREPO (e.g., Pfrommer et al. 2017), with the exception that the time-evolution of e_cr is not tracked, that is, we solved a coupled system of conservation equations

$$\begin{array}{c}\hfill \frac{\partial \mathit{U}}{\partial t}+\mathbf{\nabla }·\mathit{F}=\mathit{S}\end{array}$$(3)

for a state of conserved variables U, a flux function F, and the associated source terms S of the form

$$\begin{array}{c}\hfill \mathit{U}=\left(\begin{array}{c}\rho \\ \rho \mathit{v}\\ \rho e\\ \mathbf{B}\end{array}\right),\mathit{F}=\left(\begin{array}{c}\rho \mathit{v}\\ \rho \mathit{v}{\mathit{v}}^T+{\mathrm{P}}_{\mathrm{tot}}-\mathbf{B}{\mathbf{B}}^T\\ (\rho e+{\mathrm{P}}_{\mathrm{tot}})\mathit{v}-\mathbf{B}(\mathit{v}·\mathbf{B})\\ \mathbf{B}{\mathit{v}}^T-\mathit{v}{\mathbf{B}}^T\end{array}\right),\mathit{S}=\left(\begin{array}{c}0\\ \mathbf{0}\\ {\mathit{S}}_{\mathbf{e}}\\ \mathbf{0}\end{array}\right),\end{array}$$(4)

where S_e = P_cr ∇  ⋅  v² − v_st  ⋅  ∇ P_cr³ + Λ_net is the energy source term, e = u + $\frac{{|\mathit{v}|}^2}{2}$ + $\frac{{|\mathbf{B}|}^2}{2\rho }$ is the total energy per unit mass, with u being the thermal energy per unit mass, and Λ_net is the net cooling rate, which includes the ė_{th, gas, cr} term described above. In practice, note that the above equations are solved in comoving coordinates, for which a number of variables are transformed and redefined, as extensively discussed in Pfrommer et al. (2017).

Throughout this work, we only consider SNe (in practice, star-forming gas cells) as CR sources. While one could additionally add analogous source terms arising, for instance, from relativistic jets driven by AGN, or acceleration in cosmic shocks, we forego these contributions in this preliminary study. Such contributions will be important in massive elliptical galaxies where black holes drive powerful jets (e.g., Omma et al. 2004; Teyssier et al. 2011), and can particularly dominate the CR injection fraction in the absence of young populations of massive stars. For this reason, we restrict our current analysis to Milky Way-mass halos and below, with M_200c ≲ 10¹² M_⊙.

In the main body of work we also intentionally neglect CR cooling terms due to the excitation of Alfvén waves (Wiener et al. 2017; Ruszkowski et al. 2017; Buck et al. 2020), that is, Alfvén heating of the gas. Appendix A demonstrates that this process is negligible in our simulations.

3. Results

3.1. Selecting a fiducial set of parameters

We begin with a brief exploration of v_st, eff and κ_eff – the effective streaming velocity and diffusion coefficient – that are the two free parameters of the cosmic ray transport model. Fig. 1 explores this parameter space in terms of the resulting stellar-to-halo mass relation at z = 0. Black curves show the median relation from the realization with the fiducial TNG model, which is in reasonable agreement with observational/empirical constraints (Behroozi et al. 2013; Moster et al. 2013). Analogous median relations from the various CR runs are shown through the various colored curves. The bottom panels show the ratio of the various CR cases to the TNG model, that is, the deviation of the (median) stellar mass, at fixed halo mass.

Fig. 1. Impact of changing vst, eff and κeff, the two free parameters of the CR model, on the stellar-to-halo mass relation at z = 0. We varied κeff at fixed vst, eff (left panels), and vice versa (right panels). The black curve shows the relation from the realization with the fiducial TNG model (no CRs), and the various colored curves are simulations with CRs. The bottom panels show the ratio of the various CR cases to the TNG model, where the dashed black (gray) horizontal lines show values of one (0.25 and 0.1). In most models with CRs, the stellar mass growth is significantly suppressed. Varying κeff at fixed vst, eff has a minimal impact, while larger values of vst, eff at fixed κeff lead to higher stellar masses, closer to TNG values. Throughout the paper, we adopt values of 100 km s−1 (5 × 1028 cm2 s−1) vst, eff (κeff) as our fiducial set of transport parameters (red curve, labeled TNG-CR). For this case, we also contrast against a simulation without SNe/SMBH feedback (NF; dashed orange curve). Finally, for comparison in subsequent plots, we also consider a weaker model (the blue curve, labeled TNG-Weak CR).

To choose values for these two CR parameters, we first draw inspiration from the calibration of Hopkins et al. (2023). By contrasting a set of properties against simulations run using an explicit CR transport scheme with a constant anisotropic diffusivity κ_∥ = 3  ×  10²⁹ cm² s⁻¹ (Hopkins et al. 2020), they found that values of (v_st, eff, κ_eff) = (20 km s⁻¹ ,  5  ×  10²⁸ cm² s⁻¹) offer a good fit. We therefore explore variations about this starting point, since they are in agreement with observationally favored values (see discussion below).

In the left panels, we fixed v_st, eff = 20 km s⁻¹ and varied κ_eff by one order of magnitude in each direction. The median curves corresponding to κ_eff = 5  ×  [10²⁷, 10²⁸, 10²⁹] cm² s⁻¹ are shown in green, teal and violet, respectively. In all cases, a clear suppression in stellar mass is seen at fixed halo mass. Starting from a value of ∼ 0.8 at M_200c ∼ 10¹⁰ M_⊙, the ratio (lower panel) reaches a minima of ∼ 0.1 at the Milky Way-mass (M_200c ∼ 10¹² M_⊙). A similar qualitative behavior was reported by Hopkins et al. (2020), that is, for simulations that were run with a more advanced CR model, in which stellar masses were only suppressed in halos with a mass higher than M_200c ≳ 10¹¹ M_⊙.

Unexpectedly, the three curves are similar over almost the entire range of M_200c, that is, varying κ_eff by two orders of magnitude has minimal impact on the stellar-to-halo mass relation, except for the high-mass end (M_200c ≳ 10¹³ M_⊙) where the red curve begins to separate out from the others. While Hopkins et al. (2020) also noted that the diffusion coefficient has a minimal impact on stellar mass growth at masses M_200c ≲ 5  ×  10¹¹ M_⊙, they found that larger values of κ_∥ led to greater suppression in stellar mass for more massive halos (M_200c ≳ 5  ×  10¹¹ M_⊙). The difference must arise from the simplicity of the CR model we use and/or differences in other baryonic physics models or numerical schemes between the simulations.

In the right panel of Fig. 1, we fix κ_eff = 5  ×  10²⁸ cm² s⁻¹, and vary v_st, eff = [20,  100,  500] km s⁻¹, shown in teal, red and blue, respectively. Although stellar masses (at fixed halo mass) are once again displaced with respect to the black curve, a monotonic trend is visible: the level of suppression is reduced toward higher values of v_st, eff, particularly at the low-mass end. The red curve (v_st, eff = 100 km s⁻¹) returns stellar masses within a factor of ≲ 4 of the fiducial TNG result at all masses, with the difference maximal at the Milky Way-mass range. An even larger v_st, eff of 500 km s⁻¹ brings the median to within a factor of ≲ 1.5 at the low mass end (M_200c ≲ 10¹² M_⊙), and in fact over-estimates the stellar mass in more massive halos. As we discuss below, this is due to delayed SMBH growth, making the AGN feedback model of TNG relatively less effective above the Milky Way-mass range.

It is encouraging that larger values of v_st, eff produce stellar-to-halo mass relations similar to TNG, as such a model can likewise be consistent with observational constraints. It is unclear, however, whether these CR transport parameters are physically plausible. Hopkins et al. (2021a) suggest that the narrow range in parameter space defined by ${\kappa }_{\mathrm{eff}}+{\mathit{v}}_{\mathrm{st},\mathrm{eff}}{r}_{\odot }^{\mathrm{gal}}\sim O({10}^{29}{\mathrm{cm}}^2{\mathrm{s}}^{-1})$, where $r_{\odot }^{\mathrm{gal}}\sim 8$kpc is the galactocentric distance of the sun, is needed to be consistent with observations in the solar neighborhood (see also Korsmeier & Cuoco 2022).

Irrespective of the value assumed for κ_eff, an effective streaming speed v_st, eff = 500 km s⁻¹ is on the higher side (v_st, eff $r_{\odot }^{\mathrm{gal}}$ ∼ 10³⁰ cm² s⁻¹). A value of v_st, eff = 100 km s⁻¹, paired with either κ_eff = 5  ×  10²⁷ or 5  ×  10²⁸ cm² s⁻¹, satisfies the above constraint. Both values of κ_eff yield largely similar results, and we adopt the latter throughout the rest of this work for computational convenience, that is, we fix κ_eff = 5  ×  10²⁸ cm² s⁻¹ and v_st, eff = 100 km s⁻¹, and label this TNG-CR.

In order to offer a comparison with a case where the stellar-to-halo mass relation is more realistic, we also include a subset of results in a number of plots from the run with v_st, eff = 500 km s⁻¹, which we label TNG-Weak CR. We again stress, however, that this case differs from observational CR constraints (Hopkins et al. 2021a), and is only meant as a proxy for a case where the impact of CRs is relatively weak.

The final comparison we make throughout this paper is with a case without feedback. Here, we adopt the fiducial CR values of κ_eff = 5  ×  10²⁸ cm² s⁻¹ and v_st, eff = 100 km s⁻¹, but completely disable all stellar and SMBH feedback processes. We label this TNG-CR-NF (dashed orange line, upper right panel of Fig. 1). As seen, this differs from TNG-CR (red curve) only for relatively small halos (M_200c ≲ 10¹¹ M_⊙), beyond which the two are largely similar. This suggests that the suppression of stellar masses (at fixed halo mass) is dominated by the effect of CRs, and furthermore that no recalibration that is, weakening of the TNG feedback model could bring the results into agreement with data. We discuss the implications of this further below.

Before expanding our view on the statistics of galaxy properties across the population, we consider the spatial distribution of cosmic rays on cosmological scales. Fig. 2 shows a visualization⁴ of the large-scale structure of the CR energy density (e_cr) across the simulated volume of the TNG-CR realization. Brighter colors correspond to larger e_cr, as shown by the colorbar. CR energy densities are as high as ≳10⁵ K cm⁻³ at the centers of galaxies, corresponding to sites of active star formation. CR energy decays as distance increases, in a roughly spherical symmetric fashion owing to the functional form of the distance-modulation factor (F(r_ij); Section 2.2). This gives rise to cosmic ray halos that permeate throughout, and extend beyond, dark matter halos.

Fig. 2. Visualization of the projected CR energy density through the simulated volume (at z = 0) of the realization with (vst, eff, κeff) = (100 km s−1 ,  5  ×  1028 cm2 s−1), i.e., the fiducial set of parameters we use. The image extends 25 h−1 ∼ 37 Mpc along the x- and y-axes of the plane, as well as in the projection direction. CR energy densities are as high as ≳105 K cm−3 at the centers of galaxies, i.e., at the sites of star formation, and decay with increasing distance. Halos with multiple such sources, i.e., contributions from satellite galaxies in addition to the central, have multiple bright dots, leading to a certain degree of asymmetry in the CR energy density profile; the halo toward the lower left corner of the image is a representative example. The corresponding nonthermal CR pressure, as well as secondary effects due to associated heating plus ionization terms, has a noticeable impact on the evolution of galaxies, as we shall explore throughout the rest of this work.

Halos with multiple sites of star formation, that is, those that have star-forming satellites in addition to a star-forming central, have multiple such halo features. The superposition of these gives rise to a certain degree of asymmetry in the e_cr structure. The halo toward the lower left corner of the image is a representative example. In addition, complexity in the gas distribution in and around galaxies, that is, in ρ and ∇ρ, also lead to deviations from perfectly spherical e_cr profiles (Δτ_cr, i; Section 2.2). Such features are most clearly visible toward the outer edge of halos, sometimes as linear or filamentary-like streaks between halos. The corresponding nonthermal pressure, as well as secondary effects due to associated heating and ionization terms, has a noticeable impact on the evolution of galaxies, as we shall explore throughout the rest of this work.

3.2. Integrated and large-scale halo properties

Fig. 3 shows six properties and relations of the galaxy population: the star formation rate density, stellar mass function, stellar mass to halo mass relation, black hole mass versus stellar mass, stellar sizes, and halo gas fractions. Black curves always show the fiducial TNG model, contrasting against TNG-CR (red) and TNG-Weak CR (blue). Gray points and shaded bands show various observational and empirical constraints. Note that these are exactly the observables used in the calibration of the TNG model, undertaken in the relative sense with respect to the outcome of the original Illustris simulation (Pillepich et al. 2018). We therefore reflect on how TNG, in its calibration space, is impacted by the inclusion of cosmic rays. Except for the top left panel, all curves show the median, and all give z = 0 results.

Fig. 3. Summary of galaxy- and halo-scale integrated properties of the galaxy population, all at z = 0, except for the top left panels. In all panels, the black curves correspond to the fiducial TNG model, and the red (blue) curves correspond to TNG-CR (TNG-Weak CR). Scatter points show various observational and empirical constraints: Baldry et al. (2008), Giodini et al. (2009), Baldry et al. (2012), Bernardi et al. (2013), Behroozi et al. (2013), McConnell & Ma (2013), Moster et al. (2013), D’Souza et al. (2015), Oesch et al. (2015), Lovisari et al. (2015), Bland-Hawthorn & Gerhard (2016). The fiducial TNG model is in reasonable agreement with these observational constraints, as this is the set of observables used for calibration. The TNG-CR simulation however produces a reduced number of galaxies at M⋆ ≳ 108.5 M⊙ (center left) as a result of lower star formation rate densities over cosmic epochs (top left). Despite the suppressed stellar mass at fixed halo mass (lower left), and lower black hole masses in general (top right), the galaxy size-stellar mass relation (center right) and halo gas fractions (lower right) are within a factor of ∼ few between TNG-CR and the fiducial TNG case, as well as the different data points (see main text for a more detailed discussion).

The top left panel shows the cosmic star formation rate density (SFRD) as a function of redshift. At z ≳ 2, the SFRD in the TNG-CR case (red) is suppressed by a factor ∼4 − 6 with respect to TNG (black), with the disparity between the two growing smaller at z ≲ 2. The inclusion of CRs shifts the peak of the SFRD to lower redshift, while the shape at high redshift is similar. Coincidentally, the z = 0 SFRD of TNG-CR reaches the same value of TNG. Neglecting all possible complications including obscuration effects, both CR models shown are in significant tension with data (Behroozi et al. 2013; Oesch et al. 2015).

The reduced global star formation rates naturally lead to galaxies being typically under-massive in their stellar content, which we show through the galaxy stellar mass function (Φ, middle left panel). Solid curves show Φ for central galaxies, that is, the galaxies we primarily focus on in this work, and dashed curves for all galaxies, that is, centrals and satellites. In both these sets of populations, the number density of M_⋆ ≳ 10^{8.0 − 8.5} M_⊙ galaxies is lower in TNG-CR (red) versus TNG (black). The TNG-Weak CR case typically lies in between these two, except for z ≲ 1 where the SFRD is larger. An imprint of this enhanced star formation activity is visible as a larger number density of M_⋆ ≳ 10^11.0 M_⊙ galaxies in the center left panel. Broadly speaking, the dashed all galaxies lines should be compared against the observational constraints (Baldry et al. 2008; Bernardi et al. 2013; D’Souza et al. 2015). While all models are roughly consistent with data in the dwarf regime, both CR models are ruled out by the observed abundance of galaxies at the knee, that is, Milky Way-mass systems.

For a given halo mass, the inclusion of CRs has a noticeable impact on the stellar mass content of galaxies. The lower left panel shows the stellar-to-halo mass relation at z = 0, now normalizing stellar mass by halo mass. The TNG-CR simulation clearly has a different behavior when compared to TNG (black), both in terms of shape and normalization. A peak, if it exists, is marginal and moved to roughly one order of magnitude higher in halo mass. In contrast, TNG-Weak CR (blue) is qualitatively similar to TNG, but horizontally offset by ∼0.5 dex. Neither are consistent with semi-empirical models (Behroozi et al. 2013; Moster et al. 2013).

In addition to the reduction of stellar mass, the growth of supermassive black holes is also hindered. The black hole versus halo mass relation is shown in the top right panel with the dashed curves outlined in orange, with x-axis values (i.e., halo mass) given by the top (orange) x-axis. At M_200c ≳ 10¹¹ M_⊙, TNG-CR is offset by ∼0.5 dex below TNG, signifying smaller black holes at fixed halo mass. Interestingly, the median relation in TNG-CR falls closer to the mass of SgrA* in the Milky Way (Bland-Hawthorn & Gerhard 2016), shown through the gray star with an orange outline, which is otherwise difficult to reproduce in MW-mass galaxies (Pillepich et al. 2024). Once again, the TNG-Weak CR case lies in between the fiducial TNG and TNG-CR outcomes, as expected.

The upper right panel also shows the black hole-to-stellar mass relation, with the solid curves and lower (black) x-axis. In this case, the TNG and TNG-CR models vary only by a factor of a ∼few, and are both in broad agreement with the observed relation reported by McConnell & Ma (2013). The TNG-Weak CR case is largely similar to TNG-CR in this case, that is, at fixed stellar mass, the mass of the black hole is largely independent of v_st, eff. This is expected from the functional form of F(r_ij) (Section 2), since the distance modulation factor is largely independent of v_st, eff at small r_ij, and thus has a minimal impact on the accretion rate of black holes, which is rather modulated by the density of the local gas reservoir.

The center right panel shows the galaxy size-mass relation, with the former defined as the stellar half mass radius. All three models are largely similar, and all broadly consistent with observations (Baldry et al. 2012). At these masses, the population is dominated by blue galaxies, and so best compared to the circular gray data points. Although galaxies are under-massive for their parent halo mass in both TNG-CR and TNG-Weak CR, the size comparison demonstrates that some galaxy properties are realistic, that is, largely unchanged, for a given galaxy stellar mass.

Last, we show the halo gas fraction as a function of the halo mass in the lower right panel. To guide the eye, the dashed gray line is drawn at the cosmic baryon fraction, that is, Ω_b/Ω_m ∼ 16%. While the three curves vary by a few percent at halo masses M_200c ≲ 10¹³ M_⊙, they are in good agreement at higher masses. For simplicity we plot the FoF gas mass normalized by the total FoF mass, which can only roughly be compared to observational constraints (Giodini et al. 2009; Lovisari et al. 2015), against which there is reasonably agreement in all models. Analogous observations at lower mass scales will play an important role in discriminating between various galaxy formation models in the future (Zhang et al. 2024).

Having focused on integral galaxy population statistics and scaling relations, we now transition to focus on the properties of gas in and around galaxies, to better understand the physical impact of CRs. We begin with a visual motivation in Fig. 4, which contrasts projections of four quantities in the fiducial TNG case (upper panels) versus TNG-CR (center panels) and TNG-CR-NF (no feedback; lower panels). From left to right, we show temperature, magnetic field strength, (total) pressure, that is, thermal plus nonthermal, and metallicity.

Fig. 4. Visual impression of a subset of gas quantities through the simulated volumes for the fiducial TNG (upper panels), TNG-CR (centre panels) and TNG-CR-NF (lower panels) runs. From left to right, we show temperature, magnetic field strength, (total) pressure, and metallicity. All images extend ∼37 Mpc along the plane and the perpendicular direction, i.e., through the entire (z = 0) simulated volume, and scale bars in all panels correspond to 15 Mpc. The large-scale structure of temperature and pressure is similar between the three runs, but for the lack of hot, rarified bubbles around massive galaxies in TNG-CR. The runs differ significantly in the other quantities, with lower levels of magnetization and metal enrichment of halo gas in the two CR runs as compared to TNG. The similarity in total pressure despite weaker field strengths, implying lower levels of magnetic pressure, shows that the nonthermal support by CRs is significant.

On these large scales, the distributions of temperature and total pressure are largely similar between the three runs, except for the lack of hot rarified bubbles around massive halos (M_200c ≳ 10¹² M_⊙) in TNG-CR. In the (fiducial) TNG runs, these are understood to be produced by kinetic mode AGN feedback (Pillepich et al. 2021), the dominant channel for massive black holes (M_BH ≳ 10⁸ M_⊙; Weinberger et al. 2017). Given that black holes are less massive at fixed halo mass in the CR runs, these features are therefore suppressed in the lower panels.

We remark on an interesting qualitative result. Although stellar and black hole feedback processes are turned off in TNG-CR-NF, and thus it contains no AGN driven outflows, bubbles similar to TNG are visible in this case. We speculate that this reflects outflows that are purely driven by CR pressure gradients. This phenomenon was previously seen in idealized simulations and spatial scales that are smaller by one order of magnitude (Hanasz et al. 2013; Salem & Bryan 2014; Girichidis et al. 2016; Huang & Davis 2022), and in the FIRE-2 cosmological simulations on Mpc scales (Hopkins et al. 2021b).

The similarity of the (total) pressure images is surprising, since the large-scale magnetic field strengths are notably different between the two runs. While the centers of halos have similar levels of magnetization, halo gas surrounding these galaxies has much weaker field strengths. Despite the significant drop in magnetic pressure (P_B ∼ |B|²), the total pressure is nearly unchanged. This occurs because cosmic rays instead provide a similar level of nonthermal support, that is, cosmic ray pressure naturally compensates for the missing magnetic pressure.

Gas in the circumgalactic medium of galaxies also has much lower metallicity in the CR model runs. This is a natural consequence of the lower stellar masses of galaxies, and thus overall less efficient production of metals across cosmic time. As we discuss below, lower feedback-driven outflow velocities in the CR simulations also prevent metals from being distributed to as large of spatial scales into the intergalactic medium, instead being confined closer to their production sites.

In Fig. 5, we zoom-in to smaller scales around a single Milky Way-mass halo, cross-matched across simulations. As before, we contrast the fiducial TNG run (top row) to TNG-CR (middle row) and TNG-CR-NF (no feedback; lower row). Images span ±R_200c from edge to edge, and the white scale bars always indicate 100 kpc. From left to right, we show temperature, density, the ratio of cosmic ray to magnetic pressure, and the ratio of nonthermal to thermal pressure.

Fig. 5. Visual impression of a subset of gas quantities around a Milky Way-mass halo in the fiducial TNG (upper panels), TNG-CR (center panels), and TNG-CR-NF (lower panels) runs, cross-matched across the various realizations. Images extend ±R200c in all cases, along the plane and in the perpendicular direction, and scale bars correspond to a 100 kpc. From left to right, we show temperature, density, ratio of PCR to Pmag, and ratio of Pnon, th(= PCR + Pmag) to Pth. While profiles of temperature are largely similar across the three runs, the density of gas is lower in the TNG-CR-NF run, particularly in the halo. In the two CR runs, PCR dominates over Pmag in all regions except the central disk, where the two are in rough equipartition. In all three runs, only the central region of the halo is dominated by nonthermal pressure support. The fraction of nonthermal pressure in outer regions of the halo is relatively higher in the CR runs than in Fid. TNG, however.

Similar to Fig. 4, the spatial structure of gas temperature is broadly similar across the three runs: the disk is dominated by a cool ∼10⁴ K component, and is embedded in a halo of hot ≳10^5.5 K gas. Hotter ∼10⁶ K gas is evident in the TNG model, and somewhat suppressed in the CR runs. The extent and abundance of the cool gas disk is more prominent in the CR runs, although we attribute this primarily to a timing offset between the simulations, and not to any underlying physical cause.

On the contrary, the density structure of gas clearly differs between the runs – most notably, gas densities are lower in the TNG-CR-NF run, particularly in the halo gas surrounding the central disk. While densities in TNG-CR are also lower in comparison to the fiducial TNG run, the level of disparity is smaller.

In the CR runs, magnetic pressure is in rough equipartition with CR pressure in the disk, as indicated by the black color. In all other regions of the halo, owever, the magnetic fields are clearly subdominant to CRs by up to a factor of ∼10³. We note that this is in qualitative agreement with results from the FIRE-2 simulations (Hopkins et al. 2020). Last, the fourth panel of Fig. 5 shows that only the innermost regions of the halo are dominated by nonthermal pressure support, in all three runs. As one transitions farther out into the halo, thermal pressure begins to dominate, although to a much smaller degree in the CR runs as compared to fiducial TNG.

We again note the emergence of a feature visible in TNG-CR and even in TNG-CR-NF: bubble-like outflows that are nonthermal pressure dominated. These are clearly evident in both magnetic and CR pressure, although not particularly so in temperature. The spatial scales of these bubbles is ∼10 − 50 kpc, much smaller than the features previous discussed in Fig. 4, and reminiscent of the Fermi/eROSITA-like bubbles seen in TNG50 (Pillepich et al. 2021). Because they are present in the NF (no feedback) case and thus are clearly sourced by CR pressure gradients, however, they are of a different origin. We defer a detailed study of these bubbles and their properties to future work.

Fig. 6 provides a more quantitative assessment of these trends, through spherically averaged radial profiles of multiple gas properties. In all panels, black curves correspond to the fiducial TNG run, and red to TNG-CR. We select (central) galaxies whose parent halo masses are M_200c ∼ 10¹¹, 10^11.5 and 10¹² M_⊙, shown through dotted, dashed and solid lines, respectively. For the fiducial TNG case, note that we only select those galaxies in which the black hole in the quasar mode, thereby excluding galaxies with strong AGN feedback from our analysis. This enables us to focus on interpreting the impact of CRs with respect to star formation and stellar feedback, also facilitating comparison with previous work. All curves show median values across galaxies in the respective bin. For the Milky Way-mass bin alone (M_200c ∼ 10¹² M_⊙), we also show median curves from TNG-CR-NF (orange) and TNG-Weak CR (blue). Insets in all panels show individual profiles for Milky Way-mass halos in the TNG-CR run. Curves are colored by the virial mass of the halo, with the darkest (brightest) colors corresponding to M_200c ∼ 10^11.8 (10^12.2) M_⊙. For lower halo masses the temperatures in the CGM are noticeably colder, which is in agreement with models of isolated halos (e.g., Girichidis et al. 2024, see also Thomas et al. 2025).

Fig. 6. Spherically averaged stacked radial profiles of six gas properties. Clockwise from the top left, we show: temperature, density, nonthermal to thermal pressure ratio, metallicity, magnetic field strength, and total pressure. Colors contrast TNG-CR (red) versus the fiducial TNG run (black), while the line styles correspond to different halo mass bins. For the 1012 M⊙ mass bin in each panel, we also show profiles from the CRs but no feedback run (TNG-CR-NF; orange) and TNG-Weak CR (blue). The inclusion of CRs impacts most properties: halo gas is less magnetized and metal enriched in the TNG-CR and TNG-CR-NF runs. Profiles in TNG-Weak CR are always closer to those of fiducial TNG. Most notably, the nonthermal pressure support of CGM gas (≳0.2 R200c) is relatively larger in the CR cases as compared to the TNG run, although profiles in the galaxy (≲0.2 R200c) are largely similar across the various runs. Insets in all panels show individual profiles of Milky Way-mass halos in TNG-CR, with curves colored by halo mass (see text).

We begin with temperature in the top left panel. In all three mass bins, profiles are similar between TNG and TNG-CR. At the Milky Way mass, as well as the other two bins (not shown), TNG-CR-NF yields a similar radial profile, suggesting that these temperature profiles are largely modulated by gravity, and only to a secondary extent by feedback processes. We note, however, that kinetic mode feedback in TNG can produce a noticeable impact on temperature structure by driving outflows of hot, super-virial gas (Pillepich et al. 2021; Ramesh et al. 2023a), which we have intentionally excluded from these curves to avoid complication of interpretation of results (see above). The TNG-Weak CR curve is also similar, albeit offset to higher temperatures by ∼0.05 − 0.1 dex.

The radial density profiles are shown in the top right panel. TNG and TNG-CR are similar in the inner halo (≲0.2 R_200c), although densities are typically lower by ∼0.1 − 0.2 dex in the latter at larger galactocentric distances. While the blue (TNG-Weak CR) solid curve lies intermediate of the black and red, the orange (TNG-CR-NF) is vertically offset further to lower values, by 0.1 − 0.5 dex depending on distance. In all the CR cases, our interpretation is that the added CR pressure inhibits gas accretion into the halo, thereby reducing total mass and overall gas densities. In the TNG-CR-NF case, star formation rates of low-mass halos (M_200c ≲ 10¹¹ M_⊙; Fig. 1) are enhanced in the absence of feedback processes, giving rise to larger CR pressures in the halo, and eventually to lower densities as they evolve into their z = 0 descendants. We note that the trend of reduced halo gas densities in the presence of CRs is in qualitative agreement with Hopkins et al. 2020.

In the center panels, we focus on quantities related to the pressure of gas. In the center left panel, the total pressure profiles are similar for the TNG and TNG-CR runs, as well as TNG-Weak CR, albeit with minor ∼O(0.1) dex variations, much like those in temperature and density. In TNG-CR-NF (orange), the total pressure is lower by 0.2 − 0.4 dex, depending on distance. This is primarily due to lower densities giving rise to lower thermal pressures, but also in part because of weaker halo magnetic fields, as we elaborate further below.

The type of pressure support, however, differs significantly between TNG and the CR runs, as shown by the ratio of nonthermal to thermal pressure (center right panel). While the profiles of TNG and TNG-CR are largely similar in the inner halo (≲0.2 R_200c), where both are dominated by nonthermal pressure to similar extents, profiles begin to separate out at larger distances. In particular, the TNG ratio of P_non, th/P_th monotonically drops as one transitions to the outer regions of the halo, to as low as ∼0.03 at ∼0.6 R_200c. On the contrary, the TNG-CR (and TNG-Weak CR) profiles flatten at ∼0.3, roughly independent of halo mass. Thus, while the total pressure support is similar between TNG and TNG-CR, thermal pressure dominates to a smaller extent in the CR runs. We mention that results from the FIRE-2 simulations are significantly different, both qualitatively and quantitatively. In FIRE-2, nonthermal pressure support dominates the CGM in halos of ∼10¹² M_⊙ (Hopkins et al. 2020), which is not the case in our CR runs. We make a more quantitative comparison below.

Previous observational studies have inferred that the CR energy density within the Milky Way ISM may be in rough equipartition with analogous magnetic and thermal reservoirs (Cox 2005; Naab & Ostriker 2017). For our CR runs, this would imply a pressure ratio P_non, th/P_th of about 2⁵. In all mass bins in TNG-CR, the ratio is ≳10 at the very centers of galactic disks (≲0.05 R_200c). In the outer disk and the inner CGM (0.05 ≲ r / R_200c ≲ 0.3), however, the ratio is bracketed between ∼ [0.5, 10], that is, an approximate equipartition is attained between thermal energy and the total nonthermal component at these intermediate distances. Further out in the halo (≳0.1 − 0.2 R_200c), magnetic field strengths are weak, and no longer in equipartition with other components (see discussion below). In such a case, P_non, th (∼ P_cr)/P_th is expected to be close to 0.5 if the thermal and CR energies are in equipartition. Intriguingly, profiles of all mass bins flatten close to this value at distances ≳0.3 R_200c. We suspect that this near perfect balance between the two energy reservoirs may be linked to the assumption of a steady-state e_cr solution by the Hopkins et al. (2023) model (Section 2.2).

As discussed previously, halo magnetic fields are typically weaker in the CR runs, which we quantify further in the lower left panel of Fig. 6. Comparing the TNG and TNG-CR profiles, it is apparent that field strengths are weaker in the CR runs by as large a factor as ≳30 in the outer reaches of the halo. The dashed and dotted lines show a similar qualitative behavior, indicating its universality with halo mass, as to the TNG-Weak CR and TNG-CR-NF profiles, suggesting that any level of CR pressure tends to lower CGM magnetic field strengths.

In all three mass bins, the B field strength is largely unaffected at the very centers of galaxies, however. We suggest that this is because fields in the innermost regions of halos are primarily amplified by small-scale dynamos (Pakmor & Springel 2013). This exponential growth phase proceeds until magnetic and turbulent kinetic energy reservoirs reach rough equipartition, after which a linear growth phase takes place in the outer-disk through differential galactic rotation (Pakmor et al. 2017). This explains why the TNG and TNG-CR profiles begin to differ at ∼0.05 − 0.1 R_200c; since galaxies are smaller at fixed halo mass in TNG-CR, this linear growth phase is effectively damped, yielding weaker fields. Finally, magnetic fields in the halo are closely linked to the transport of magnetized gas from the galaxy into the CGM (Pakmor et al. 2020; Ramesh et al. 2023c), suggesting that the lower level of magnetization in the outer halo in TNG-CR is likely a result of reduced mass-outflow rates.

A similar trend is seen in the case of metallicity profiles (lower right panel), where the metal enrichment of the gas in the halo outskirts is lower by a factor of ≳3 − 5 in TNG-CR as compared to the fiducial run. This also indicates suppressed outflow rates, as the transport of metal-enriched gas from the galaxy and into halo plays an important role in enriching the CGM with metals (Péroux et al. 2020; Pandya et al. 2021). Note however that the metallicity profiles are also offset at the halo centers, albeit to a lower level than the CGM, a direct impact of the tight relation of the stellar mass to the gas-phase metallicity (Tremonti et al. 2004; Ma et al. 2016; Torrey et al. 2019).

Last, the insets in all panels show that profiles of individual halos are diverse. A variation of ≲1 dex at fixed distance is seen in most cases, often correlating with the mass of the parent halo. For instance, more massive halos typically have the hottest CGMs, most pressure support and strongest magnetic fields. In cases like metallicity and pressure ratio, the scatter is largely random with respect to halo mass, suggesting that these variations are driven by other factors.

In Fig. 7, we dig deeper into the outflow dynamics of gas. The top panel shows outflow velocity, while the lower shows mass-outflow rate. At a galactocentric distance r, the latter is defined as ${\dot{M}}_{\mathrm{out}}(r)=1/\mathrm{\Delta }r\sum m_i{\mathit{v}}_{\mathrm{rad},\mathrm{i}}$ where the sum is over all gas cells i with mass m_i and radial velocity v_rad, i that have positive radial velocities (v_rad, i > 0, i.e., outflowing) and reside within a shell of thickness Δr at a radius r (i.e., |r_i − r|≤Δr/2). The outflow velocity is computed as the mean radial velocity of the same subset of gas cells. We fix Δr to be 0.05 R_200c, and note that results are qualitatively similar for other choices of Δr.

Fig. 7. Mass-outflow rates (lower panel) and outflow velocity (upper panel) as a function of distance. Colors show different runs, and line styles different halo mass bins. In all mass bins, the outflow rate and velocities are suppressed throughout the halo in TNG-CR with respect to the fiducial TNG case. In the Milky Way-mass range, the outflow rate is suppressed further in TNG-CR-NF, and the radial trend reverses, leading to larger velocities and outflow rates in the outer versus inner halo (≳0.5 R200c). TNG-Weak CR is typically intermediate between TNG and TNG-CR. We note, however, that the stellar masses are different in the CR runs as compared to fiducial TNG (see main text).

As before, we focus on central galaxies of halos in bins of M_200c ∼ 10¹¹, 10^11.5 and 10¹² M_⊙, shown in dotted, dashed and solid curves, respectively. Curves in black correspond to fiducial TNG, red to TNG-CR, orange to TNG-CR-NF and blue to TNG-Weak CR. Although the latter two are only shown for the Milky Way-mass bin, we mention that they portray qualitatively similar trends for the other two mass bins as well. For the Fid. TNG case, we once again restrict the analysis to those halos in which the black hole is active in the quasar mode.

In all three mass bins, the mass-outflow rates are suppressed in TNG-CR as compared to TNG. The offset is roughly 0.1 − 0.2 dex in the galaxy, but rises to ∼0.3 − 0.5 dex as one transitions further out into the halo. The TNG-Weak CR result (blue) lies between the TNG and TNG-CR values and shows a qualitatively similar behavior, that is, increasing from the center of the galaxy to a distance of ∼0.1 R_200c (typical extent of the gas disk), and dropping at larger distances. The orange curve, however, is qualitatively different, both in shape and normalization: while the curve also has a peak at ∼0.1 R_200c, thereafter it rises toward R_200c and the outer halo.

In the absence of supernova and black hole feedback, this outflow is purely driven by gradients in CR pressure (e.g., Quataert et al. 2022; Modak et al. 2023; Martin-Alvarez et al. 2023), and is different in nature from outflows driven by other physical mechanisms (Simpson et al. 2016; Thomas et al. 2025). Most importantly, the mass-outflow rate is lower at all distances than in the other curves, that is, at a fixed halo mass, outflows that are only driven by CR pressure gradients are overall weaker than other feedback processes. It is important to note, however, that stellar masses (and star formation rates) are not constant for the different physics variations, and mass-loading factors portray additional behaviors (see Fig. B.1).

Outflow velocities follow similar trends: in all three mass bins, at most distances, velocities are smaller in the TNG-CR run as compared to TNG. Only at the center of halos are the two similar, a direct result of the TNG model, where the wind velocity at injection is set by the local velocity dispersion of dark matter (Pillepich et al. 2018), and is thus primarily set by the halo mass, modulo any local disturbances. Interestingly, outflow velocities in the TNG-CR-NF run (orange) are in fact larger than the red curve at distance ≳0.4 R_200c, likely as a result of smaller ambient gas densities which provide lower levels of resistance as the outflow propagates outwars. We speculate that this enhanced outflow velocity, which increases toward larger distances, is the primary driver of large-scale bubbles seen in Fig. 4, which are absent in TNG-CR since velocities are typically lower and decay toward the halo outskirts.

As outflow rates and velocities in TNG-CR are clearly suppressed with respect to TNG, this will impact magnetic field strengths and metallicites in the CGM. Fig. 8 shows two dimensional plots of radial velocity (y-axis) as functions of magnetic field strength (x-axis) with background color indicating mean gas-phase metallicity, contrasting TNG (left) versus TNG-CR (right). We focus exclusively on outflowing gas (v_rad ≥ 0) within R_200c of the centrals of Milky Way-mass halos, that is, with satellite gas excised. As before, halos in which the black hole is in the kinetic mode are excluded. Dashed white curves show contours drawn at the [1, 10, 50]% levels⁶, increasing in level toward lower velocities. For comparison purposes, the dotted black curves on the right replicate the contours from the left panel.

Fig. 8. Two-dimensional plots of the radial velocity (y-axis) as functions of the magnetic field strength (x-axis). The background pixels are colored by the mean gas-phase metallicity to compare the Fid. TNG run (left) with TNG-CR (right). We focus exclusively on outflowing gas (vrad ≥ 0) within R200c of the centrals of Milky Way-mass halos, i.e., without satellite gas. The dashed white curves show contours drawn at the [1, 10, 50] % levels, increasing in level toward lower velocities. For ease of visualization, the dotted black curves on the right overplot the contours from the left panel. A greater fraction of outflows gas in Fid. TNG extends to velocities as large as O(100) km s−1, carrying metal enriched gas outward from the galaxy. The Fid. TNG contours are shifted toward slightly larger field strengths, implying the transport of relatively higher magnetized gas by these outflows.

Contrasting equi-level contours between the TNG and TNG-CR simulations reveals that a greater fraction of gas in the fiducial TNG model case is outflowing at large velocities, as high as ∼O(100) km s⁻¹. This is in agreement with the upper panel of Fig. 7. Background colors further indicate that the slower outflows in TNG-CR carry gas that is less metal enriched, at the level of 0.01 − 0.1Z_sun. We reiterate the complication, however, that the metallicites of the galaxies themselves differ across the two runs (lower right panel of Fig. 6).

Contrasting contours also shows that outflows in TNG-CR are shifted toward weaker magnetic field strengths. The slower outflows in TNG-CR thus also transport relatively weakly magnetized gas, yielding overall smaller field strengths in the halo (lower left panel of Fig. 6). There is also a clear correlation between weakly magnetized and metal pristine gas, suggesting that the increasing mass-outflow rates at larger halocentric distances in the CR models accelerate in situ halo gas, rather than ejecting ISM gas, resulting in a different mass loading trend with radius (see also Mitchell et al. 2020).

In Fig. 9, we conclude this subsection with density-temperature phase diagrams of gas within Milky Way-mass halos (M_200c ∼ 10¹² M_⊙). The panel on the left corresponds to the fiducial TNG run, and the right to TNG-CR. As before, halos with kinetic-mode, that is, strong SMBHs are excluded. All halos in the mass bin are stacked, and all gas within R_200c is considered, with satellite gas excised. Color shows the relative mass of gas in each bin, that is, normalized by the bin with the maximum mass. The dashed black contours are drawn at the 10% level, and the two crosses correspond to (local) maxima of cold (< 10^4.5 K) and hot (> 10^5.5 K) gas. For ease of comparison, the dotted curves and gray crosses on the right copy the corresponding features from the left panel.

Fig. 9. Phase diagrams of gas within R200c of Milky Way-mass halos in the Fid. TNG run (left) vs. TNG-CR (right). The crosses in both panels correspond to (local) maxima of the cold (< 104.5 K) and hot (> 105.5 K) phases, while dashed curves show contours drawn at the 10% level. For ease of visualization, the dotted curves and gray crosses on the right overplot the corresponding features from the left panel. While the peak density of the cold phase is larger by ∼0.3 dex in Fid. TNG with respect to TNG-CR, the hot phase peak is only offset by ∼0.1 dex in density, and ∼0.05 dex in temperature. Furthermore, contours are more broad in Fid. TNG as compared to TNG-CR, implying larger variation in these gas properties at fixed halo mass.

In density-temperature space, contours enclosing fixed fractions of mass are broader in TNG than in TNG-CR, that is, the distribution of densities and temperatures is narrower in TNG-CR. Furthermore, there is a systematic shift of the local maxima (i.e., mean values) toward lower densities in TNG-CR. This is more pronounced in the cold phase (∼0.3 dex) as opposed to the hot phase (∼0.1 dex). While not shown here, we mention that the peak densities of star forming gas are also lower in TNG-CR (∼0.5 dex). This provides a compelling explanation for the suppressed star formation rates and stellar masses in the CR models, as lower gas density translates directly into lower SFR (Kennicutt 1998). We attribute the overall lower densities primarily to a direct impact of the added CR pressure support suppressing the amount of gas accretion onto halos and galaxies, resulting in lower densities (Fig. B.2; see also Buck et al. 2020).

A minor offset in the location of the temperature maxima for the hot phase (of the CGM) is visible – temperatures are slightly higher in TNG versus TNG-CR (∼0.05 dex), in agreement with the top left panel of Fig. 6. This is likely a combination of effects: increased star formation activity in TNG results in more supernovae and thus overall more energy injection from stellar feedback, while the presence of nonvanishing CR pressure at constant total pressure implies somewhat lower thermal pressures. Second order differences are also visible: larger fractions of warm ∼10⁵ K gas, as well as more super-virial ≳10⁶ K gas, in TNG than in TNG-CR.

Overall, while properties like the pressure ratio, magnetic field strength and metallicity are significantly impacted by the inclusion of CRs, either directly or indirectly, the temperature, density and total pressure structure of gas are modified to a smaller extent. It is notable, however, that removing both SNe and SMBH feedback yields halos that are underpressurized and slightly less dense (Fig. 6), and thus some of these results may also be sensitive to the physics employed for other galactic processes. Future projects with additional physics variations will help explore such degeneracies further.

3.3. Comparison with other CR galaxy simulations

So far we have qualitatively discussed our results in the context of previous work from the literature. We now contrast our simulations against other projects in a more quantitative manner, in particular with those that employ more complex and sophisticated models for the transport physics of cosmic rays.

We begin with a face-value comparison between the results of our simulations versus those from a set of isolated galaxy simulations presented in Girichidis et al. (2024). In particular, we contrast our findings against their runs which employ a state-of-the-art spectral CR scheme (Girichidis et al. 2020) coupled to MHD (Girichidis et al. 2022). These simulations also include radiative cooling and star formation using the Springel & Hernquist (2003) model, but no strong galactic-scale feedback.

Fig. 10 shows spherically averaged radial profiles of cosmic ray pressure (P_cr = e_cr/3) out to a distance of ∼0.25 R_200c⁷, that is, focused in on the galaxy, the disk-halo interface, and the inner CGM. Solid, dashed and dotted curves correspond to M_200c ∼ 10¹², 10¹¹ and 10¹⁰ M_⊙ halos, respectively. In the Girichidis et al. (2024) simulations, the galaxies at the center of these halos have stellar masses M_⋆ ∼ 10^10.5, 10^9.5 and 10^8.5 M_⊙, shown through the different black curves. Red curves correspond to median TNG-CR profiles of galaxies matched with respect to stellar mass, while profiles of the halo mass matched TNG-CR sample are shown in orange. For the solid red curve alone, we additionally plot the 16^th-84^th variation across the sample.

Fig. 10. Comparison of spherically averaged radial profiles of the CR pressure (Pcr = ecr/3) between TNG-CR and an advanced spectral CR scheme (Girichidis et al. 2024). Solid, dashed and dotted curves correspond to M200c ∼ 1012, 1011 and 1010 M⊙ halos, respectively. In the Girichidis et al. (2024) simulations, the galaxies at the center of these halos have stellar masses M⋆ ∼ 1010.5, 109.5 and 108.5 M⊙, shown through the different black curves. Red curves correspond to TNG-CR profiles of galaxies matched with respect to stellar mass, while profiles of the halo-matched TNG-CR sample is shown in orange. Only for the solid red curve, we additionally plot the 16th–84th variation across the sample. The stellar mass matched sample, except for the least massive bin, roughly matches the black curves (within a factor of ∼1.5). The halo mass matched sample however is offset by much larger factors, owing to suppressed stellar masses at fixed halo mass.

Fig. 11. Comparison of spherically averaged radial profiles of the CR pressure (Pcr; lower panel) and the ratio of nonthermal-to-thermal pressures (upper panel) between TNG-CR and other cosmological simulations: FIRE-2 and Auriga. Gray and black curves correspond to halos from FIRE-2, coral to multiple CR physics variations of Auriga halo 6 (Au6), and red (stellar mass matched) and orange (halo-mass matched) to those from TNG-CR. While the Pcr profiles are in broad agreement between TNG-CR and FIRE-2, the type of pressure support in these halos is significantly different: in TNG-CR, halos are dominated by nonthermal (thermal) support in the inner (outer) regions, largely independent of halo mass. In contrast, low (high) mass halos in FIRE-2 are almost always dominated by thermal (nonthermal) support. Results of Au6 depend on the CR physics, but profiles of pressure ratios are typically more qualitatively similar to TNG-CR than to FIRE-2.

Starting from a P_cr/k_B of ∼10^4.1 K cm⁻³ at the center of the galaxy, the Girichidis et al. (2024) profile for their Milky Way-mass halo (solid black) steadily declines by over an order or magnitude to ∼10^2.5 K cm⁻³ at ∼0.25 R_200c. The M^⋆-matched TNG-CR (solid red) profile shows the same qualitative behavior, but is vertically offset downward ∼0.2 dex at all distances. At a fixed stellar mass of 10^10.5 M_⊙, the CR pressure support in our TNG-CR runs (with the transport parameters we select) is thus underestimated by a factor of ∼1.5. Note however that the percentile regions are relatively broad, ∼± 0.4 dex at all distances, and the P_cr profiles in some of these galaxies are in excellent agreement with Girichidis et al. (2024).

A similar trend is seen between the dashed curves, that is, for M_⋆ 10^9.5 M_⊙ galaxies. At even lower stellar masses (M_⋆ 10^8.5 M_⊙; dotted curves), the offset is larger (a factor of ∼10). This difference arises because CR transport at this galaxy stellar mass is dominated by advection for prolonged periods of time in the Girichidis et al. (2024) run, thereby giving rise to greater CR pressure support in the halo. In more massive halos, increased gravity limits advective transport to shorter intervals of time, beyond which diffusion dominates. Assuming a universal constant v_st, eff at all halo masses and at all cosmic epochs, as is required by the Hopkins et al. (2023) model and adopted in our CR runs, neglects such effects. We note however that P_cr in TNG-CR at M_⋆ ∼ 10^8.5 M_⊙ may also be under-estimated as a result of numerical resolution (see Appendix C).

Trends are qualitatively similar with the halo matched sample, although quantitatively different. For instance, the M_halo-matched (solid orange) is offset low by ∼0.6 dex at all distances, that is, by a factor of ∼4 with respect to Girichidis et al. (2024) (solid black), as a result of suppressed stellar masses at fixed halo mass (Fig. 1). The offset is larger for lower halo masses (up to a factor of ∼20), but this is likely due to the stellar mass of the Girichidis et al. (2024) run being larger than the empirically value motivated value (black curve in Fig. 1), due to the idealized nature of their setup. Note that the dotted orange curve, which corresponds to galaxies that are very close to unresolved at our resolution (M_⋆ ≲ 10^7.0 M_⊙), is absent in this plot.

Overall we conclude that, at fixed stellar mass, our simple CR model and assumed transport parameters result in a reasonable match with the P_cr profiles of Girichidis et al. (2024), based on much more sophisticated CR physics. This agreement is reassuring and demonstrates that the net impact of CRs in our runs is plausible, at least from an energetics point of view. As stellar masses are suppressed in TNG-CR, the halo-mass matched sample is naturally offset in terms of pressure support.

In Fig. 11, we contrast our results against those from the FIRE-2 project (Hopkins et al. 2020, 2023). These simulations use a single energy, that is, gray approximation for CRs, they are fully cosmological in nature and hence make it possible to extend the above analysis to larger, CGM scales. The FIRE-2 runs we contrast against employ an explicit CR transport scheme with a constant anisotropic diffusivity κ_∥ = 3  ×  10²⁹ cm² s⁻¹, and include ideal MHD, supernova feedback, metal line cooling, amongst other processes (see Hopkins et al. 2018, 2020).

These runs are shown through the various gray (m11b, m11i) and black (m12i, m12f, m12m) curves. The former correspond to halos of mass M_200c ≲ 10¹¹ M_⊙, while the latter to ∼ 10¹² M_⊙. Profiles of the halo-mass matched sample from TNG-CR is shown in orange, and stellar mass matched in red. Note that the latter is shown only for the Milky Way mass, since the less massive FIRE-2 halos have stellar masses M_⋆ ∼ 10⁸ M_⊙, at which scale we are limited by numerical resolution (Appendix C).

The lower panel focuses on radial profiles of cosmic ray pressure, P_cr. Remarkably, we find that the stellar mass matched sample is in excellent agreement with the corresponding FIRE-2 sample. The halo-mass matched sample is offset vertically by ∼ 0.4 dex at all distances, a consequence of suppressed stellar masses at fixed halo mass (Fig. 1). The gray curves portray a similar shape as the dashed orange, but are offset vertically below, since these halos are in fact slightly less massive than 10¹¹ M_⊙. That is, a one-to-one comparison is not possible in terms of halo mass simply due to limited numerical resolution at the low mass end (Appendix C). Overall, we conclude that the simple transport model and parameters in TNG-CR give a decent match with expected nonthermal CR pressure support predicted by the FIRE-2 + CRs models.

The upper panel of Fig. 11 compares profiles of the ratio of nonthermal to thermal pressure. To guide the eye, the dotted horizontal violet line demarcates the regime of nonthermal pressure dominance (above) from that of thermal support (below). Note that only the m11b and m12i curves from FIRE-2 are available, and so are shown here.

The behavior of these pressure ratios is markedly different from that of the total CR pressure profiles. While halos in TNG-CR are dominated by nonthermal pressure support in the inner halo (≲0.2 R_200c), gas at larger distances is instead supported by thermal pressure, roughly independent of the two halo masses considered here. In stark contrast, the ∼10¹² M_⊙ FIRE-2 halo is dominated by thermal support only at the very center of the galaxy, and transitions into a nonthermal pressure support regime throughout the rest of the halo. The ∼10¹¹ M_⊙ FIRE-2 halo is thermal pressure dominated at all distances, that is, the halo-mass dependence of this finding also differs.

The reason for this inversion of pressure ratio across mass scales in FIRE-2, as well as the difference with respect to TNG-CR, is unclear. As Ji et al. (2020) discuss, the physical properties of gas in the CGM of FIRE-2 galaxies with CRs are substantially altered. Not only does the mean temperature decrease by levels on some order of magnitude, but the trend of temperature with density is inverted. Overall, these simulations feature a cooler, CR-dominant CGM, substantially altering predicted observables, for instance, for metal ion absorption as well as X-ray emission, and likely representing an extreme case for the impact of CRs which is in tension with such data (Ji et al. 2020).

There are few additional cosmological, that is, nonidealized galaxy formation simulations incorporating CR physics. Notably, Buck et al. (2020) present resimulations of two Auriga Milky Way-mass halos with variations on CR physics. In Fig. 11 we also compare to the Au6 run from that work, contrasting simulations assuming pure advection based CR transport (Adv; coral dashed), with added anisotropic diffusion using a diffusion coefficient κ_∥ = 10²⁸ cm² s⁻¹ (Diff; coral dotted), and further with the inclusion of Alfvén cooling of CRs (DiffAlfven; coral dot dashed). Overall, the CR pressure profiles are steeper than TNG-CR, being larger in the inner halo and lower at ∼R_200c/2. With respect to the pressure ratio, the Buck et al. (2020) results depend on the CR physics. Some models yield trends more qualitatively similar to TNG-CR than to FIRE-2, although the Diff case produces a nonthermal pressure dominated inner CGM extending to scales intermediate between the TNG-CR and FIRE-2 result.

In addition to the difference in CR-schemes used, note that the diffusion coefficient in both the FIRE-2 and Auriga simulations we compare against do not perfectly match the value set in TNG-CR. This likely also contributes to the differences between the various curves, further complicating this comparison. In Appendix E, we offer a simple comparison across the suite of TNG-CR variation runs, giving a zeroth-order idea of how transport parameters may impact pressure profiles, at least within the realm of the simple CR-scheme that we use. Last, differences in baryonic feedback models as well as the underlying hydrodynamical scheme, may all give rise to significant differences. The impact of CRs on galaxy and CGM properties clearly eludes current consensus.

3.4. Discussion and outlook

In TNG-CR, the inclusion of CRs leads to a suppression of stellar masses at fixed halo mass (Fig. 1), a direct consequence of reduced star formation rates over cosmic epochs (Fig. 3). There are at least four possible mechanisms: (1) the added nonthermal CR pressure support may reduce accretion rates of CGM gas onto the galaxy, or the accretion rate of gas into the halo from larger scales (Hopkins et al. 2020; Ruszkowski & Pfrommer 2023), (2) the increased ionization and (3) heating terms may reduce the amount and/or density of cold gas in the galaxy or the CGM, or (4) outflows driven by CR pressure gradients (Pakmor et al. 2016; Girichidis et al. 2018) may suppress galactic star formation rates. These reflect both preventative and ejective feedback mechanisms.

In our models, options (2) and (3) do not appear to cause a significant impact, as the magnitude of these contributions is subdominant to other relevant ionization/cooling terms. As discussed in Fig. 7, our CR driven outflows are not generally as powerful as other feedback-driven outflows. This is likely because the star formation rates are suppressed, however, which effectively reduces the CR energy density. If stellar masses were not underestimated for a given halo mass (Fig. 1), and CR pressure gradients were stronger, it is possible that corresponding outflows could compete with other mechanisms in terms of energetics, as has been proposed in the past (Pfrommer et al. 2017; Girichidis et al. 2018; Chan et al. 2022).

In TNG-CR, we therefore suggest that the added CR pressure support in the CGM leads to the suppression of galactic star formation. This is a result of inhibited gas accretion onto halos leading to lower gas fractions (lower right panel of Fig. 3), ultimately giving rise to lower gas densities in the ISM (Figure 9), that is, CRs do not directly impact star-forming gas cells in our simulations (see below).

Although it is possible that our results are biased by the limitations of the simple model that we employ, similar qualitative behaviors have been found with more advanced transport schemes (Hopkins et al. 2020; Farcy et al. 2022; Montero et al. 2024). We refer the reader to Hopkins et al. (2023) for a discussion on the possible aspects in which this model is likely to fail quantitatively. In particular, it is least realistic at the injection scale, that is, within the ISM. In addition, the Springel & Hernquist (2003) model that we adopt implies that the thermodynamic properties of star-forming ISM gas is set, that is, overwritten, by the effective equation of state (eEOS), and thus CRs have no direct impact on its phase structure. Future extensions of TNG-CR that include more sophisticated, that is, resolved, treatments for the ISM can test the impact that CRs may directly have on star-forming gas, in which case it would be important to incorporate the physics of ion-neutral damping, which may contribute to the decoupling of CRs from neutral gas (Armillotta et al. 2021; Sike et al. 2024).

As motivated in Section 3.1, we have selected transport parameters that reasonably reproduce a subset of realistic galactic-scale properties, while also being in agreement with the narrow range in parameter space proposed by Hopkins et al. (2021a). These parameters also offer a decent match to CR pressure support predicted by more advanced models (Section 3.3). Certain properties of galaxies in the TNG-CR run as well as TNG-Weak CR and other variants are clearly ruled out by available observational constraints, however. These include the cosmic SFRD, the stellar mass function, and the stellar-to-halo-mass relation (SMHM; Fig. 3). Moreover, the reduced magnetic field strengths and metallicities will alter observable rotation measures (RM) and ionic column density and covering fractions, where observational datasets and constraints are rapidly improving (e.g., Wu et al. 2025; Anderson et al. 2024).

Stellar masses are also suppressed even in the total absence of supernova and black hole feedback from the TNG model (dotted curve in the right panel of Fig. 1), particularly at the high mass end (M_200c ≳ 10¹¹ M_⊙). This suggests that, independent of feedback physics, the TNG-CR model is too strong in its impact, and furthermore, that no recalibration, that is, weakening of the TNG feedback model can bring TNG-CR into agreement with observational constraints.

How can the galaxy properties of TNG-CR be reconciled with observational constraints? For a given set of transport parameters, one approach to reduce the nonthermal CR pressure support is to adopt a lower value for the fraction of CR energy per SNe (${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$). While ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ was set to 10% for all the runs presented in the main text, Appendix D explores the impact of varying this fraction. For example, considering a value of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ = 1 % returns a SMHM relation much closer to emperical constraints. While Murase et al. (2019) infer a constraint of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ ≲ 5 − 10 %, this is based on observations of a single SNe, and is also model dependent. More observations of this kind are needed to better constrain ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ and its possible dependencies.

Within the realm of the Hopkins et al. (2023) model, an alternate approach may be to relax some subset of the underlying assumptions. In particular, by requiring that CR transport parameters are fixed in space and time, one is likely to miss CR dynamics that depend on halo mass. For instance, Girichidis et al. (2024) suggest that CR transport in low-mass halos (M_200c ∼ 10^{10 − 11} M_⊙) is dominated by advection, at least temporarily, while more massive halos are instead dominated by diffusion. This would correspond to a time- (and space-) dependent v_st, eff that increases (locally) during an advection dominated period, and reduces to smaller values when diffusion is instead the primary transport mechanism. More generally, the transport mode (diffusion/streaming) and the effective transport speed need to be computed based on local hydrodynamic properties, which naturally lead to locally (strongly) varying quantities (see, e.g., Jiang & Oh 2018; Thomas & Pfrommer 2019, 2022). Exploring such a model remains the topic of future work.

Even though they are computational expensive, more sophisticated, that is, explicit, CR models must also be explored in large-volume cosmological simulations. This includes more realistic transport, and spectral schemes that can capture the energy-dependent effects of the CR population. It remains to be seen if, and how, realistic galaxy populations can be obtained with explicit and physically plausible CR scenarios.

4. Summary and conclusions

We have used a suite of (25 Mpc h⁻¹)³ cosmological magnetohydrodynamical boxes to explore the impact of CRs on the evolution of galaxies. To do this, we coupled the IllustrisTNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018) with a simple subgrid model for CR transport (Hopkins et al. 2023). Our main findings are listed below.

1. We compared a fiducial TNG simulation with multiple models including CRs: a fiducial but strong TNG-CR parameterization, a TNG-Weak CR alternative, and a control case with CRs, but without other feedback. In general, the quantitative impact of CRs on the galaxy properties depends on the values of our two CR transport parameters, that is, on the effective diffusion coefficient and on the streaming velocity.

2. The inclusion of CRs reduces the star formation rates and suppresses galaxy stellar masses at a fixed halo mass. This effect can reach a factor of about ten in the range of the Milky Way mass, and it decreases to a factor of about two for M^⋆ ∼ 10⁹ M_⊙ dwarf galaxies (Fig. 1).

3. The results of our fiducial TNG-CR model differ strongly in several summary statistics and scaling relations from observations. The differences include the SFRD, the galaxy SMF, and the relation of the stellar mass to the halo mass. The relations of the galaxy size to mass and of the halo scale to the gas fractions are not qualitatively changed and remain in broad agreement with data. The growth of SMBHs is reduced because the densities in the central star-forming ISM are lower. This flattens the relations of the black hole mass (Fig. 3).

4. The inclusion of CRs does not significantly impact the temperature, density, or (total) pressure structure of the CGM compared to the TNG model. Halo gas, however, is supported by nonthermal pressure to a larger extent, and the magnetization and metal enrichment are lower with CRs (Figs. 4, 5, and 6).

5. The primary physical driver of the differences between TNG and TNG-CR is that (i) CRs add appreciable nonthermal pressure support to galactic halos. This suppresses gas inflow rates into both the halo and galaxy. Simultaneously, (ii) CRs reduce the mass-outflow rates and velocities of feedback-driven winds (Figs. 7 and 8).

6. An appreciable CR energy density permeates galaxies and galactic halos and extends into the IGM. We identified a variety of CR-driven outflows that are supported by nonthermal pressure and bubble-like features at a variety of scales, from 10 − 50 kpc galactic outflows to about megaparsec-sized halo-scale structures (Figs. 2, 4, and 5).

7. The phase structure of gas within the virial radii of Milky Way-like halos is qualitatively unchanged. The presence of CRs lowers the mean density (by ∼0.1 dex) and the temperature (∼0.05 dex) of the CGM (Fig. 9).

8. The assumed transport parameters mach the CR pressure profiles predicted by more advanced CR models well for isolated galaxy setups with a spectral scheme (Girichidis et al. 2024) and for the FIRE-2 cosmological simulations. The ratio of nonthermal to thermal pressure and the overall impact on the physical properties of the CGM is starkly different in our simulations from FIRE-2, however (Figs. 10 and 11).

This work is an initial study of the CR impact on different galaxy populations. It represents our first attempt at illuminating the IllustrisTNG universe with cosmic rays, and it is a first step toward including CR physics in large-volume cosmological simulations, which is an important future direction for numerical models of galaxy formation and evolution. A key and unresolved caveat is that the galaxy population of our TNG-CR model differs significantly more from observations than the original TNG model without CRs. Future work will advance our current modeling by incorporating more advanced treatments of CR production and transport physics.

Data availability

The original IllustrisTNG simulations are publicly available and accessible at www.tng-project.org/data (Nelson et al. 2019b), where the variations presented here including cosmic ray physics will be made public in the future. Data directly related to this publication is available on request from the corresponding authors. This work has benefitted from the scida analysis library (Byrohl & Nelson 2024).

Acknowledgments

We thank Annalisa Pillepich for insightful discussions and sharing relevant data, and Matthew C. Smith for valuable inputs on the tree walk implementation. RR and DN acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) through an Emmy Noether Research Group (grant number NE 2441/1-1). RR is a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). PG acknowledges funding by the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130). This work has made use of the VERA supercomputer of the Max Planck Institute for Astronomy (MPIA) operated by the Max Planck Computational Data Facility (MPCDF), and of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Impact of Alfvén heating

Cosmic rays that propagate faster than the Alfvén-velocity excite Alfvén waves through the CR streaming instability (Kulsrud & Pearce 1969), effectively heating up gas as these waves dampen on short timescales (Wiener et al. 2017; Ruszkowski et al. 2017). For the approximations made in the derivation of the CR subgrid model, the corresponding heating term is given by Λ_alf = |v_alf ⋅ ∇P_cr|/3 (Hopkins et al. 2023), where |v_alf| = |B|/(4πρ)^1/2 is the magnitude of the Alfvén-velocity.

As mentioned in Section 2.2, this process was not included in the results presented in the main body of the paper. In Fig. A.1, we briefly explore the impact of adding this heating term. The upper panel shows the stellar to halo-mass relation, while radial profiles of magnetic field strength for three different halo mass bins are shown in the bottom. Colors correspond to different variation runs, and line styles in the to lower panel to different mass bins.

In terms of both stellar mass and magnetic field strength there is minimal difference between the two runs, except for the high halo mass end in the former (M_200c ≳ 10¹³ M_⊙). This is also qualitatively the case for all other properties considered in this work (not shown). This suggests that the inclusion of Alfvén heating has a largely subdominant impact under the framework of the model we employ. We speculate that this may be due to the relative simplicity of the model, that is, the computation of P_cr through a distance-modulated locally attenuated approach results in typically small values of ∇P_cr.

Even simulations run with more sophisticated CR transport models, however, suggest that the Alfvén cooling time is typically longer than its hadronic counterpart in most regions except the very center of galaxies (Girichidis et al. 2024). This would imply that the effect of this added heating term is most pronounced on ≲ kpc scales. This regime is both poorly resolved and modulated by our subgrid ISM pressurization model. It is also where the Hopkins et al. 2023 CR model is severely limited. More sophisticated galaxy formation, ISM, and CR models can revisit this question.

Fig. A.1. The impact of including Alfvén heating. The top panel shows the stellar to halo mass relation, while the lower panel shows radial profiles of magnetic field strength. Colors correspond to different variation runs, and line styles in the to lower panel to different mass bins. Minimal difference is seen between the two runs, suggesting that the inclusion of Alfvén heating has a minimal impact under the framework of this simple model.

Appendix B: Mass loading factors and inflow rates

As discussed in Fig. 7, at fixed halo mass, gas outflow rates are suppressed in TNG-CR as compared to the fiducial TNG run. However, stellar masses (and star formation rates) vary across the different physics variations, meaning that the amount of energy available for stellar feedback driven outflows varies substantially between simulations. This makes interpretation nontrivial.

To better unravel this effect, Fig. B.1 instead shows the mass loading factor, that is, the mass-outflow rate normalized by the star formation rate of the galaxy. Halos in bins of mass M_200c ∼ [10^11.0, 10^11.5, 10^12.0] M_⊙ are shown through dotted, dashed and solid curves, respectively. Results from TNG are shown in black, and TNG-CR in red. Only for the M_200c ∼ 10^11.5 M_⊙ bin, we additionally show curves from TNG-CR-NF (orange) and TNG-Weak CR (blue), and also from a run which uses the TNG model but with no feedback, that is, no stellar and black hole feedback, as well as no CRs (Fid. TNG-NF; gray).

Contrasting the various black versus red curves, higher mass loading factors are seen in the TNG-CR runs, in all three halo mass bins. Outflows driven by the CR pressure gradients, added to the already existing SNe and black hole driven feedback in TNG, thus produce higher mass loading factors. Due to smaller stellar masses and star formation rates, this however produces a smaller outflow rate (Fig. 7).

In the TNG-CR-NF run, mass loading factors are smaller in the inner halo (≲ 0.5 R_200c) as compared to TNG and TNG-CR, and values are comparable to the TNG-CR run only at the very outskirts of the halo. With our simple CR model and assumed transport parameters, outflows driven solely by CR pressure gradients are typically weaker than other feedback processes in TNG (see also Hanasz et al. 2013; Dashyan & Dubois 2020). Similar to our general finding in the main text, the TNG-Weak CR case (blue curve) is intermediate of between TNG and TNG-CR. Last, mass-loading factors in the Fid. TNG-NF run (no feedback of any kind; gray) are by far the smallest, as expected.

As a final point of discussion on gas flows, Fig. B.2 measures inflow rates of gas through the halo. Analogous to our measurement of mass-outflow rates, at a galactocentric distance r, the inflow rate is defined as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{in}}(r)=\frac{1}{\mathrm{\Delta }r}{\displaystyle \sum _{\begin{array}{c}i=0\\ |r_i-r|\le \mathrm{\Delta }r/2\\ v_{\mathrm{rad},i}<0\end{array}}^n}-m_i\times v_{\mathrm{rad},i}\end{array}$$(B.1)

where notations carry the same definitions as earlier. As with Fig. B.1, halos in bins of mass M_200c ∼ [10^11.0, 10^11.5, 10^12.0] M_⊙ are shown through dotted, dashed and solid curves, respectively. Results from TNG are shown in black, and TNG-CR in red. Only for the M_200c ∼ 10^12.0 M_⊙ bin, we additionally show curves from TNG-CR-NF (orange) and TNG-Weak CR (blue), and also from Fid. TNG-NF (gray).

Similar to Fig. 7, in all halo mass bins and at all distances, inflow rates are smaller in TNG-CR as compared to TNG. On one hand, this is linked to the reduced outflow rates, that then result in lower levels of circulation and fountain flows (Oppenheimer & Davé 2008; Péroux & Howk 2020). However, contrasting the gray (no feedback whatsoever) versus orange (only CRs, no SNe or AGN feedback) curves suggests that the added nonthermal CR pressure support may be an important driver of this trend, that is, even in the total absence of SNe, black hole or CR driven outflows, and thus a complete lack of gas circulation between the galaxy and the CGM, the Fid. TNG-NF simulation produces larger inflow rates than the TNG-CR-NF case. Similar to discussions in the main text, this suppression of gas inflow is due to the added CR nonthermal pressure support.

Appendix C: Assessing convergence of numerical resolution

As mentioned in Section 2, all results presented in the main text are based on the L25n512 suite of simulations, that is, those run at an average baryonic mass resolution of ∼ 10⁶ M_⊙ (equivalent to TNG100-1). We here assess numerical convergence by running analogous TNG-CR boxes at increasingly coarser resolutions.

In Fig. C.1, we show CR pressure (P_cr) profiles of central galaxies from three different mass bins: M_200c ∼ [10^11.0, 10^11.5, 10^12.0] M_⊙ in dotted, dashed and solid curves, respectively. Results from L25n512 are shown in black, while lower resolution counterparts L25n256 (m_b ∼ 10⁷ M_⊙) and L25n128 (m_b ∼ 10⁸ M_⊙) in red and orange, respectively. Note that a subset of curves are absent, as these produce profiles of P_cr ∼ 0 as a result of limited resolution.

At a given halo mass, values of P_cr are typically smaller at coarser resolution. For instance, the red solid curve is offset by ∼ 0.15 dex with respect to the black solid curve throughout most of the halo, with a slightly larger difference of ∼ 0.4 dex at the very center of the galaxy. We understand this to be a direct impact of lower star formation rates at coarser levels of resolution (Pillepich et al. 2018).

A similar trend is seen while contrasting black and red dotted curves, that is, the M_200c ∼ 10^11.5 M_⊙ bin. At an even lower halo mass of M_200c ∼ 10^11.0 M_⊙, the red (dotted) curve is absent, implying that these halos are the smallest ones that one could explore in the L25n512 run, which is the cutoff that we use throughout the main body of this work (Section 2).

Fig. C.1. Similar to the bottom panel of Fig. 7, but normalized by the star formation rate of the galaxy, i.e., the mass loading factor. In all halo mass bins, values are larger in TNG-CR as compared to Fid. TNG. The added outflow driven by the presence of CRs thus boosts the mass loading factor, albeit with smaller mass-outflow rates since star formation rates as diminished.

Fig. C.2. Similar to the bottom panel of Fig. 7, but the inflow rate of gas. In all halo mass bins, and at all galactocentric distances, inflow rates are smaller in TNG-CR as compared to Fid. TNG, signifying lower levels of gas infall through the halo and into the galaxy.

Fig. C.3. Comparison of CR pressure (Pcr) profiles across simulations run at varying numerical resolutions. For a given halo mass bin, values of Pcr are typically smaller at coarser resolution, a direct impact of suppressed star formation rates. Note that a subset of curves are absent, as these produce profiles of Pcr ∼ 0 as a result of limited resolution.

Appendix D: Impact of varying ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$

Results presented throughout the main text are based on simulations where ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$, the CR energy fraction available from SNe, is set to 10 %. We here explore the impact of varying this fraction.

In Fig. D.1, we compare the (z = 0) stellar mass-(upper panel) and black hole mass-(lower panel) halo mass relations from several simulations. The black curve shows results from the TNG run, while the TNG-CR is shown in red. Orange and blue correspond to new runs with ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ = 5 and 1 %, respectively.

Reducing the value of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ yields a SMHM relation closer to the black curve, that is, to empirical constraints. For instance, the ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ = 1 % curve is in excellent agreement with the black curve up to a halo mass of 10¹² M_⊙. We understand this to be a result of lower amounts of CR pressure support in the halo (not shown), thus allowing a greater amount of gas to accrete on to galaxies, yielding larger star formation rates.

However, the most massive halos above 10¹² M_⊙ have even larger stellar masses than TNG for the lower-energy fraction cases (blue and orange). This is a result of diminished black hole masses (lower panel), thus delaying the onset of kinetic mode BH feedback, that is the main quenching mechanism in massive halos in TNG. While CR pressure support in the halo is reduced, the added nonthermal component at the center of galaxies is still large enough that SMBH accretion rates (via a Bondi model; Weinberger et al. 2017) are appreciably reduced, thereby limiting their growth.

While Murase et al. 2019 propose a constraint of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ ≲ 5 − 10 %, this is based on observations of a single SNe, and also model dependent. More observations of this kind are needed to assess whether the ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$ = 1 % case is indeed realistic. In such a model, one would additionally need to adjust the models for SMBH growth and feedback in order to reproduce TNG galaxy properties at the high-mass end.

Fig. D.1. A comparison of the stellar mass-(upper panel) and black hole mass-(lower panel) halo mass relation between simulations with varying values of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$. Lower values effectively result in weaker CR pressure support in the halo, yielding a more realistic SMHM relation. However, the growth of SMBHs is still suppressed in all cases, due to high CR pressures in the very center regions of star-forming galaxies.

Appendix E: TNG-CR variations – Pressure profiles

As discussed in Fig. 11, a number of differences are present between profiles of cosmic ray pressure (P_cr) and ratio of nonthermal to thermal pressures (β⁻¹) of TNG-CR galaxies versus those from the FIRE-2 and Auriga simulations. While the differences in cosmic ray models, and in general other numerical methods, between these different projects most certainly play a role, the transport parameters between all the runs also vary, thereby further complicating this comparison.

In order to better explore the impact of transport parameters, in Fig. E.1, we show a comparison of profiles of P_cr (lower panel) and β⁻¹ (upper panel) between a number of TNG-CR variations runs for M_⋆ ∼ 10^10.5 M_⊙ galaxies. The black curve shows the fiducial TNG-CR run, that is, [κ_eff, v_st, eff] = [5 × 10²⁸, 100] (cm² s⁻¹, km s⁻¹). The teal curve corresponds to [5 × 10²⁷, 100], that is, a variation of κ_eff at fixed v_st, eff. The red and blue curves vary v_st, eff to 20 and 500 km s⁻¹, respectively, at fixed κ_eff = 5 × 10²⁸ m² s⁻¹.

Contrasting the black versus teal curves shows that varying κ_eff at fixed v_st, eff has a minimal impact on P_cr profiles, albeit for marginal differences toward the very central regions of the halo. Profiles of β⁻¹ are rather similar in the inner halo (≲ 0.2 R_200c), and separate out by ∼ 0.15 dex at larger distances.

On the other hand, in both P_cr and β⁻¹ profiles, varying v_st, eff at fixed κ_eff has an impact at all distances, although to different degrees in the two panels based on the value of v_st, eff employed. Although not shown explicitly, results are qualitatively similar for other stellar and/or halo mass bins.

Fig. E.1. A comparison of profiles of cosmic ray pressure (Pcr; lower panel) and ratio of nonthermal to thermal pressures (β−1; upper panel) between a number of TNG-CR variations runs for M⋆ ∼ 1010.5 M⊙ galaxies. Contrasting the black versus teal curves shows that varying κeff at fixed vst, eff has a minimal impact on Pcr profiles, albeit for marginal differences toward the very central regions of the halo. Profiles of β−1 are rather similar in the inner halo (≲ 0.2 R200c), and separate out by ∼ 0.15 dex at larger distances. On the other hand, in both Pcr and β−1 profiles, varying vst, eff at fixed κeff has an impact at all distances, although to different degrees in the two panels based on the value of vst, eff employed.

All Figures

Fig. 1. Impact of changing vst, eff and κeff, the two free parameters of the CR model, on the stellar-to-halo mass relation at z = 0. We varied κeff at fixed vst, eff (left panels), and vice versa (right panels). The black curve shows the relation from the realization with the fiducial TNG model (no CRs), and the various colored curves are simulations with CRs. The bottom panels show the ratio of the various CR cases to the TNG model, where the dashed black (gray) horizontal lines show values of one (0.25 and 0.1). In most models with CRs, the stellar mass growth is significantly suppressed. Varying κeff at fixed vst, eff has a minimal impact, while larger values of vst, eff at fixed κeff lead to higher stellar masses, closer to TNG values. Throughout the paper, we adopt values of 100 km s−1 (5 × 1028 cm2 s−1) vst, eff (κeff) as our fiducial set of transport parameters (red curve, labeled TNG-CR). For this case, we also contrast against a simulation without SNe/SMBH feedback (NF; dashed orange curve). Finally, for comparison in subsequent plots, we also consider a weaker model (the blue curve, labeled TNG-Weak CR).

In the text

Fig. 2. Visualization of the projected CR energy density through the simulated volume (at z = 0) of the realization with (vst, eff, κeff) = (100 km s−1 ,  5  ×  1028 cm2 s−1), i.e., the fiducial set of parameters we use. The image extends 25 h−1 ∼ 37 Mpc along the x- and y-axes of the plane, as well as in the projection direction. CR energy densities are as high as ≳105 K cm−3 at the centers of galaxies, i.e., at the sites of star formation, and decay with increasing distance. Halos with multiple such sources, i.e., contributions from satellite galaxies in addition to the central, have multiple bright dots, leading to a certain degree of asymmetry in the CR energy density profile; the halo toward the lower left corner of the image is a representative example. The corresponding nonthermal CR pressure, as well as secondary effects due to associated heating plus ionization terms, has a noticeable impact on the evolution of galaxies, as we shall explore throughout the rest of this work.

In the text

Fig. 3. Summary of galaxy- and halo-scale integrated properties of the galaxy population, all at z = 0, except for the top left panels. In all panels, the black curves correspond to the fiducial TNG model, and the red (blue) curves correspond to TNG-CR (TNG-Weak CR). Scatter points show various observational and empirical constraints: Baldry et al. (2008), Giodini et al. (2009), Baldry et al. (2012), Bernardi et al. (2013), Behroozi et al. (2013), McConnell & Ma (2013), Moster et al. (2013), D’Souza et al. (2015), Oesch et al. (2015), Lovisari et al. (2015), Bland-Hawthorn & Gerhard (2016). The fiducial TNG model is in reasonable agreement with these observational constraints, as this is the set of observables used for calibration. The TNG-CR simulation however produces a reduced number of galaxies at M⋆ ≳ 108.5 M⊙ (center left) as a result of lower star formation rate densities over cosmic epochs (top left). Despite the suppressed stellar mass at fixed halo mass (lower left), and lower black hole masses in general (top right), the galaxy size-stellar mass relation (center right) and halo gas fractions (lower right) are within a factor of ∼ few between TNG-CR and the fiducial TNG case, as well as the different data points (see main text for a more detailed discussion).

In the text

Fig. 4. Visual impression of a subset of gas quantities through the simulated volumes for the fiducial TNG (upper panels), TNG-CR (centre panels) and TNG-CR-NF (lower panels) runs. From left to right, we show temperature, magnetic field strength, (total) pressure, and metallicity. All images extend ∼37 Mpc along the plane and the perpendicular direction, i.e., through the entire (z = 0) simulated volume, and scale bars in all panels correspond to 15 Mpc. The large-scale structure of temperature and pressure is similar between the three runs, but for the lack of hot, rarified bubbles around massive galaxies in TNG-CR. The runs differ significantly in the other quantities, with lower levels of magnetization and metal enrichment of halo gas in the two CR runs as compared to TNG. The similarity in total pressure despite weaker field strengths, implying lower levels of magnetic pressure, shows that the nonthermal support by CRs is significant.

In the text

Fig. 5. Visual impression of a subset of gas quantities around a Milky Way-mass halo in the fiducial TNG (upper panels), TNG-CR (center panels), and TNG-CR-NF (lower panels) runs, cross-matched across the various realizations. Images extend ±R200c in all cases, along the plane and in the perpendicular direction, and scale bars correspond to a 100 kpc. From left to right, we show temperature, density, ratio of PCR to Pmag, and ratio of Pnon, th(= PCR + Pmag) to Pth. While profiles of temperature are largely similar across the three runs, the density of gas is lower in the TNG-CR-NF run, particularly in the halo. In the two CR runs, PCR dominates over Pmag in all regions except the central disk, where the two are in rough equipartition. In all three runs, only the central region of the halo is dominated by nonthermal pressure support. The fraction of nonthermal pressure in outer regions of the halo is relatively higher in the CR runs than in Fid. TNG, however.

In the text

Fig. 6. Spherically averaged stacked radial profiles of six gas properties. Clockwise from the top left, we show: temperature, density, nonthermal to thermal pressure ratio, metallicity, magnetic field strength, and total pressure. Colors contrast TNG-CR (red) versus the fiducial TNG run (black), while the line styles correspond to different halo mass bins. For the 1012 M⊙ mass bin in each panel, we also show profiles from the CRs but no feedback run (TNG-CR-NF; orange) and TNG-Weak CR (blue). The inclusion of CRs impacts most properties: halo gas is less magnetized and metal enriched in the TNG-CR and TNG-CR-NF runs. Profiles in TNG-Weak CR are always closer to those of fiducial TNG. Most notably, the nonthermal pressure support of CGM gas (≳0.2 R200c) is relatively larger in the CR cases as compared to the TNG run, although profiles in the galaxy (≲0.2 R200c) are largely similar across the various runs. Insets in all panels show individual profiles of Milky Way-mass halos in TNG-CR, with curves colored by halo mass (see text).

In the text

Fig. 7. Mass-outflow rates (lower panel) and outflow velocity (upper panel) as a function of distance. Colors show different runs, and line styles different halo mass bins. In all mass bins, the outflow rate and velocities are suppressed throughout the halo in TNG-CR with respect to the fiducial TNG case. In the Milky Way-mass range, the outflow rate is suppressed further in TNG-CR-NF, and the radial trend reverses, leading to larger velocities and outflow rates in the outer versus inner halo (≳0.5 R200c). TNG-Weak CR is typically intermediate between TNG and TNG-CR. We note, however, that the stellar masses are different in the CR runs as compared to fiducial TNG (see main text).

In the text

Fig. 8. Two-dimensional plots of the radial velocity (y-axis) as functions of the magnetic field strength (x-axis). The background pixels are colored by the mean gas-phase metallicity to compare the Fid. TNG run (left) with TNG-CR (right). We focus exclusively on outflowing gas (vrad ≥ 0) within R200c of the centrals of Milky Way-mass halos, i.e., without satellite gas. The dashed white curves show contours drawn at the [1, 10, 50] % levels, increasing in level toward lower velocities. For ease of visualization, the dotted black curves on the right overplot the contours from the left panel. A greater fraction of outflows gas in Fid. TNG extends to velocities as large as O(100) km s−1, carrying metal enriched gas outward from the galaxy. The Fid. TNG contours are shifted toward slightly larger field strengths, implying the transport of relatively higher magnetized gas by these outflows.

In the text

Fig. 9. Phase diagrams of gas within R200c of Milky Way-mass halos in the Fid. TNG run (left) vs. TNG-CR (right). The crosses in both panels correspond to (local) maxima of the cold (< 104.5 K) and hot (> 105.5 K) phases, while dashed curves show contours drawn at the 10% level. For ease of visualization, the dotted curves and gray crosses on the right overplot the corresponding features from the left panel. While the peak density of the cold phase is larger by ∼0.3 dex in Fid. TNG with respect to TNG-CR, the hot phase peak is only offset by ∼0.1 dex in density, and ∼0.05 dex in temperature. Furthermore, contours are more broad in Fid. TNG as compared to TNG-CR, implying larger variation in these gas properties at fixed halo mass.

In the text

Fig. 10. Comparison of spherically averaged radial profiles of the CR pressure (Pcr = ecr/3) between TNG-CR and an advanced spectral CR scheme (Girichidis et al. 2024). Solid, dashed and dotted curves correspond to M200c ∼ 1012, 1011 and 1010 M⊙ halos, respectively. In the Girichidis et al. (2024) simulations, the galaxies at the center of these halos have stellar masses M⋆ ∼ 1010.5, 109.5 and 108.5 M⊙, shown through the different black curves. Red curves correspond to TNG-CR profiles of galaxies matched with respect to stellar mass, while profiles of the halo-matched TNG-CR sample is shown in orange. Only for the solid red curve, we additionally plot the 16th–84th variation across the sample. The stellar mass matched sample, except for the least massive bin, roughly matches the black curves (within a factor of ∼1.5). The halo mass matched sample however is offset by much larger factors, owing to suppressed stellar masses at fixed halo mass.

In the text

Fig. 11. Comparison of spherically averaged radial profiles of the CR pressure (Pcr; lower panel) and the ratio of nonthermal-to-thermal pressures (upper panel) between TNG-CR and other cosmological simulations: FIRE-2 and Auriga. Gray and black curves correspond to halos from FIRE-2, coral to multiple CR physics variations of Auriga halo 6 (Au6), and red (stellar mass matched) and orange (halo-mass matched) to those from TNG-CR. While the Pcr profiles are in broad agreement between TNG-CR and FIRE-2, the type of pressure support in these halos is significantly different: in TNG-CR, halos are dominated by nonthermal (thermal) support in the inner (outer) regions, largely independent of halo mass. In contrast, low (high) mass halos in FIRE-2 are almost always dominated by thermal (nonthermal) support. Results of Au6 depend on the CR physics, but profiles of pressure ratios are typically more qualitatively similar to TNG-CR than to FIRE-2.

In the text

Fig. A.1. The impact of including Alfvén heating. The top panel shows the stellar to halo mass relation, while the lower panel shows radial profiles of magnetic field strength. Colors correspond to different variation runs, and line styles in the to lower panel to different mass bins. Minimal difference is seen between the two runs, suggesting that the inclusion of Alfvén heating has a minimal impact under the framework of this simple model.

In the text

	Fig. C.1. Similar to the bottom panel of Fig. 7, but normalized by the star formation rate of the galaxy, i.e., the mass loading factor. In all halo mass bins, values are larger in TNG-CR as compared to Fid. TNG. The added outflow driven by the presence of CRs thus boosts the mass loading factor, albeit with smaller mass-outflow rates since star formation rates as diminished.
In the text

	Fig. C.2. Similar to the bottom panel of Fig. 7, but the inflow rate of gas. In all halo mass bins, and at all galactocentric distances, inflow rates are smaller in TNG-CR as compared to Fid. TNG, signifying lower levels of gas infall through the halo and into the galaxy.
In the text

	Fig. C.3. Comparison of CR pressure (Pcr) profiles across simulations run at varying numerical resolutions. For a given halo mass bin, values of Pcr are typically smaller at coarser resolution, a direct impact of suppressed star formation rates. Note that a subset of curves are absent, as these produce profiles of Pcr ∼ 0 as a result of limited resolution.
In the text

Fig. D.1. A comparison of the stellar mass-(upper panel) and black hole mass-(lower panel) halo mass relation between simulations with varying values of ${ϵ}_{\mathrm{cr}}^{\mathrm{SNe}}$. Lower values effectively result in weaker CR pressure support in the halo, yielding a more realistic SMHM relation. However, the growth of SMBHs is still suppressed in all cases, due to high CR pressures in the very center regions of star-forming galaxies.

In the text

Fig. E.1. A comparison of profiles of cosmic ray pressure (Pcr; lower panel) and ratio of nonthermal to thermal pressures (β−1; upper panel) between a number of TNG-CR variations runs for M⋆ ∼ 1010.5 M⊙ galaxies. Contrasting the black versus teal curves shows that varying κeff at fixed vst, eff has a minimal impact on Pcr profiles, albeit for marginal differences toward the very central regions of the halo. Profiles of β−1 are rather similar in the inner halo (≲ 0.2 R200c), and separate out by ∼ 0.15 dex at larger distances. On the other hand, in both Pcr and β−1 profiles, varying vst, eff at fixed κeff has an impact at all distances, although to different degrees in the two panels based on the value of vst, eff employed.

In the text