The Role of Cosmic-ray Transport in Shaping the Simulated Circumgalactic Medium

As an Eulerian code, Enzo models gas as a fluid moving through grid cells that are fixed in space. At every time step, Enzo advances the state of the fluid in the simulation by numerically approximating a solution to the equations below. These equations encompass the conservation of mass (Equation (1)), conservation of momentum (Equation (2)), the induction equation (Equation (3)), and the energy equation (Equation (4)). In simplified terms, these equations assert that the change over time in the value of a given conserved quantity in one cell (the $\tfrac{\partial }{\partial t}$ term ) is equal to the flux of that conserved quantity through the boundaries of that cell (the ∇ · () term). Source terms that encompass both energy gains and losses (e.g., an injection of CR energy after a supernova explosion) appear on the right-hand side of the equality. Traditionally, the evolution of thermal gas is encompassed in Equations (1)–(4) (setting the CR pressure term P_c to zero). We model the evolution of CRs as an additional ultra-relativistic proton fluid with adiabatic index γ_c = 4/3 (Jun et al. 1994; Drury & Falle 1986) in Equation (5),

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot (\rho {{\boldsymbol{vv}}}^{T}+{P}_{g}+{P}_{c})=-\rho {\rm{\nabla }}{\boldsymbol{\Phi }}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{B}}{v}^{{\rm{T}}}-{{\boldsymbol{vB}}}^{{\rm{T}}})=0\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\varepsilon }_{g}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{v}}{\varepsilon }_{g})=-{P}_{g}{\rm{\nabla }}\cdot {\boldsymbol{v}}+{H}_{c}+{{\rm{\Gamma }}}_{g}+{{\rm{\Lambda }}}_{g}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\varepsilon }_{c}}{\partial t}+{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{c}}}=-{P}_{{\rm{c}}}{\rm{\nabla }}\cdot {\boldsymbol{v}}-{H}_{{\rm{c}}}+{{\rm{\Gamma }}}_{{\rm{c}}}+{{\rm{\Lambda }}}_{{\rm{c}}}.\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{{\rm{c}}}={\boldsymbol{v}}{\varepsilon }_{{\rm{c}}}+{{\boldsymbol{v}}}_{{\rm{s}}}({\varepsilon }_{{\rm{c}}}+{P}_{{\rm{c}}})-{\kappa }_{\varepsilon }{\boldsymbol{b}}({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}{\varepsilon }_{{\rm{c}}}).\end{eqnarray} \tag{ 6 }$$

Here, we define Φ to be the gravitational potential (where ${{\boldsymbol{\nabla }}}^{2}{\rm{\Phi }}=4\pi {\boldsymbol{G}}{\rho }_{\mathrm{tot}}$), ρ_tot to be the total density, and ρ, ${\boldsymbol{v}}$ to be the gas density and velocity. ${\boldsymbol{B}}$ is the magnetic field strength, and ${\boldsymbol{b}}={\boldsymbol{B}}/| {\boldsymbol{B}}|$ is the magnetic field direction. The superscript T denotes the vector transpose. The internal gas pressure, P_g is related to the internal thermal energy density ε_g = (γ_g − 1)P_g, where γ_g = 5/3. Similarly, the CR pressure is related to the CR energy density, ε_c = (γ_c − 1)P_c, with γ_c = 4/3. From here on, we use the subscript "g" to refer to properties of the thermal gas and the subscript "c" to refer to properties of the CR fluid. We model the diffusion coefficient κ_ε as a constant, which is observationally constrained to be on the order of κ_ε ≃ 10²⁸ cm² s⁻¹ (Strong & Moskalenko 1998; Ptuskin et al. 2006; Tabatabaei et al. 2013). Although we neglect CR transport perpendicular to the magnetic field, it could have observable effects (Kumar & Eichler 2014). Γ and Λ are energy source and sink terms, respectively. In our simulations, the source of CR energy is supernova events. We do not implement hadronic CR energy losses.

Equation (6) encompasses the three modes of CR transport that we have implemented. The first term on the left represents advection, in which the CR fluid moves with the bulk velocity of the thermal gas. This term is solved explicitly in the Riemann solvers of Enzo. The next two terms describe CR streaming and diffusion, respectively. These are approximations to CR motion relative to the gas and are therefore solved separately from the advection equation. The implementation of these transport methods is described in detail in Sections 2.2 and 2.3.

In the streaming approximation, CRs move relative to the gas with a velocity given by

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\rm{s}}}=-\mathrm{sgn}({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}{\varepsilon }_{{\rm{c}}}){f}_{s}{{\boldsymbol{v}}}_{A}.\end{eqnarray} \tag{ 7 }$$

Here, f_s is the constant streaming factor described in Ruszkowski et al. (2017). The Alfvén velocity ${{\boldsymbol{v}}}_{A}={\boldsymbol{B}}/\sqrt{4\pi \rho }$ represents the transverse waves propagating along the magnetic field lines in a plasma. The function sgn returns the sign of the enclosed expression. In this limit, CRs also transfer momentum to the thermal gas through the heating term H_c where

$$\begin{eqnarray}&&{H}_{c}=| {{\boldsymbol{v}}}_{A}\cdot {\rm{\nabla }}{P}_{c}| .\end{eqnarray} \tag{ 8 }$$

Although this term appears only in the streaming approximation, we follow the example of Wiener et al. (2017) and include it in some simulations with CR diffusion to isolate the underlying differences in the different transport mechanisms (see Table 1).

Table 1.  Simulation Parameters

Run ID	Min. Grid Size (pc)	CR Parameters
		fc	Transport Mode	κε [×1028 cm2 s−1]	fs	Hc
ncr	160	⋯	No CRs	⋯	⋯	no
	320	⋯	No CRs	⋯	⋯	no
adv	160	0.1	Advection only	0	⋯	no
	320	0.01	Advection only	0	⋯	no
	320	0.1	Advection only	0	⋯	no
	320	0.3	Advection only	0	⋯	no
isod	160	0.1	Isotropic diffusion	3	⋯	no
	320	0.01	Isotropic diffusion	3	⋯	no
	320	0.1	Isotropic diffusion	3	⋯	no
	320	0.3	Isotropic diffusion	3	⋯	no
	320	0.1	Isotropic diffusion	1	⋯	no
aanisd	160	0.1	Anisotropic diffusion	3	⋯	no
	320	0.01	Anisotropic diffusion	3	⋯	no
	320	0.1	Anisotropic diffusion	3	⋯	no
	320	0.1	Anisotropic diffusion	10	⋯	no
	320	0.3	Anisotropic diffusion	3	⋯	no
	320	0.1	Anisotropic diffusion	1	⋯	no
anisdh	160	0.1	Anisotropic diffusion	3	⋯	yes
	320	0.1	Anisotropic diffusion	3	⋯	yes
stream	160	0.1	Streaming	⋯	4	yes
	160	0.1	Streaming	⋯	1	yes
	320	0.01	Streaming	⋯	4	yes
	320	0.1	Streaming	⋯	1	yes
	320	0.1	Streaming	⋯	2	yes
	320	0.1	Streaming	⋯	4	yes
	320	0.3	Streaming	⋯	4	yes

Note.

^aA summary of the initial conditions of the simulation. The run ID and the corresponding boldface entries describe the fiducial runs discussed in depth throughout this paper. The fiducial runs differ from each other only by the invoked CR transport mechanism. We deviate from the fiducial runs by changing the resolution, the fraction of supernova energy injected as CRs (${f}_{{\rm{c}}}$), the diffusion coefficient (κ_ε), the streaming factor (f_s), and by including the CR heating term (H_c). Figure 11 compares the effect of these variables on the CR pressure distribution in the CGM.

Download table as:  ASCIITypeset image

Both CR streaming and diffusion are approximations to CRs scattering off of Alfvén waves. The difference lies in the source that is generating these waves. In the streaming approximation, CRs are assumed to drive the growth of Alfvén waves through the streaming instability (Kulsrud & Pearce 1969). This is often referred to as the "self-confinement" case (Zweibel 2017).

We isolate the streaming behavior from the general CR transport equation (Equation (6)) using

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\varepsilon }_{c}}{\partial t}-{\rm{\nabla }}\cdot [{{\boldsymbol{v}}}_{{\rm{s}}}({\varepsilon }_{{\rm{c}}}+{P}_{{\rm{c}}})]=0.\end{eqnarray} \tag{ 9 }$$

At each simulation time step, we calculate CR streaming by updating the value of the CR energy density in each cell. The evolution of the CR energy density in cell i is given by

$$\begin{eqnarray}&&{\varepsilon }_{c,i}^{n+1}={\varepsilon }_{c,i}^{n}-{\rm{\Delta }}t\displaystyle \sum _{j}-\mathrm{sgn}({{\boldsymbol{b}}}_{{ij}}\cdot {\boldsymbol{\nabla }}{\varepsilon }_{{\rm{c}},{ij}}^{n}){{\boldsymbol{v}}}_{A,{ij}}\cdot {{\boldsymbol{n}}}_{{ij}}{({\rm{\Delta }}{x}_{{ij}})}^{-1},\end{eqnarray} \tag{ 10 }$$

where ${\varepsilon }_{c,i}^{n}$ is the value of the CR energy density in cell i before streaming is applied. Δt is the time step, and the terms ${{\boldsymbol{b}}}_{{ij}}$ and ∇ε_c,ij are the direction of the magnetic field and the gradient of CR energy density computed at cell face j. ${{\boldsymbol{n}}}_{{ij}}$ describes the plane parallel to face j, and Δx_ij is the length of the cell axis that is perpendicular to face j. If the cells in the simulation are constructed as cubes, then Δx can be treated as a constant and moved out of the summation term.

The streaming time step is set by the bulk motion of the gas and the Alfvén velocity, so that ${t}_{\mathrm{stream}}\lt \tfrac{{\rm{\Delta }}x}{{{\boldsymbol{v}}}_{g}+{f}_{s}{{\boldsymbol{v}}}_{a}}$. We avoid instabilities near local extrema of CR energy density by employing the regularization described in Sharma et al. (2009) and Ruszkowski et al. (2017), where

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{c}=-({e}_{c}+{p}_{c}){{\boldsymbol{u}}}_{A}\tanh ({h}_{c}\hat{{\boldsymbol{v}}}\cdot {\rm{\nabla }}{e}_{c}/{e}_{c}).\end{eqnarray} \tag{ 11 }$$

We follow the recommendation of Ruszkowski et al. (2017) and set the free regularization parameter h_c to 10 kpc.

In the "extrinsic turbulence" case, Alfvén waves are excited by turbulent motions in the thermal gas. Since the gyroradius of CRs around magnetic field lines is significantly smaller than the best-resolved cells in our simulation, CR transport in this regime is modeled as diffusion parallel to the direction of magnetic field lines. At each time step, in addition to solving the conservation equations, Enzo computes CR diffusion according to

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\varepsilon }_{c}}{\partial t}-{\rm{\nabla }}\cdot [{\kappa }_{\varepsilon }{\boldsymbol{b}}({\boldsymbol{b}}\cdot {\rm{\nabla }}{\varepsilon }_{c})]=0.\end{eqnarray} \tag{ 12 }$$

With diffusion, the updated CR energy density in a grid cell i follows the prescription

$$\begin{eqnarray}&&{\varepsilon }_{c,i}^{n+1}={\varepsilon }_{c,i}^{n}+{\rm{\Delta }}t\displaystyle \sum _{j}{\kappa }_{\varepsilon }({{\boldsymbol{b}}}_{{ij}}\cdot {\rm{\nabla }}{\varepsilon }_{c,{ij}}^{n}){{\boldsymbol{b}}}_{{ij}}\cdot {{\boldsymbol{n}}}_{{ij}}{({\rm{\Delta }}{x}_{{ij}})}^{-1},\end{eqnarray} \tag{ 13 }$$

where Δt is the diffusion time step.

In the case of isotropic diffusion, Equation (13) becomes

$$\begin{eqnarray}&&{\varepsilon }_{c,i}^{n+1}={\varepsilon }_{c,i}^{n}+{\rm{\Delta }}t\displaystyle \sum _{j}{\kappa }_{\varepsilon }{\rm{\nabla }}{\varepsilon }_{c,{ij}}^{n}\cdot {{\boldsymbol{n}}}_{{ij}}{({\rm{\Delta }}{x}_{{ij}})}^{-1}.\end{eqnarray} \tag{ 14 }$$

We employ an explicit diffusion scheme because we find that our simulation time step is not limited by the diffusion time step constraint of ${t}_{\mathrm{diff}}\lt \tfrac{1}{2N}\tfrac{{\rm{\Delta }}{x}^{2}}{{\kappa }_{\varepsilon }}$, where N is the dimensionality of our simulation. For a detailed discussion on numerically modeling semi-implicit anisotropic CR diffusion, see Pakmor et al. (2016a).

For some geometrical configurations, using a simple estimate of the CR energy gradient when calculating anisotropic diffusion violates the second law of thermodynamics and leads to an unphysical flow of CR energy from cells with lower energy density to cells with higher energy density. In more extreme cases, this leads to some cells developing unphysical (negative) values of CR energy. To solve this problem, we employ a limited gradient as described in van Leer (1977).

The simple gradient only considers the flux through a given cell face, so that

$$\begin{eqnarray}&&{\rm{\nabla }}{\varepsilon }_{c}=\displaystyle \frac{{\varepsilon }_{c,i+1}-{\varepsilon }_{c,i}}{{\rm{\Delta }}{x}_{i}},\end{eqnarray} \tag{ 15 }$$

where i and i + 1 are the indices of neighboring cells.

Where the simple gradient produces an unphysical flux, we estimate the limited CR energy density gradient on the interface (Sharma & Hammett 2007; Pakmor et al. 2016a) to be

$$\begin{eqnarray}&&{\rm{\nabla }}{\varepsilon }_{c}=\displaystyle \frac{4}{{\displaystyle \sum }_{n}{({\rm{\nabla }}{\varepsilon }_{\mathrm{cr}}^{n})}^{-1}}.\end{eqnarray} \tag{ 16 }$$

This approximation takes into account the component of the gradient that is parallel to the cell face, n.

Star formation in our simulations follows the prescription described in Cen & Ostriker (1992) with minor modifications. Stellar particles are only formed in grids at the maximum level of refinement that meets predetermined density, mass, and minimum dynamical time thresholds: ${\rho }_{\mathrm{cell}}\gt {\rho }_{\mathrm{thres}}$, ${M}_{\mathrm{cell}}\gt {M}_{\mathrm{thres}}$, and ${t}_{\mathrm{cool},\mathrm{cell}}\lt {t}_{\mathrm{dyn},\min }$. The exact values for ${\rho }_{\mathrm{thres}},{M}_{\mathrm{thres}},{t}_{\mathrm{dyn},\min }$ depend on the resolution of the simulation, because a highly resolved cell may not be capable of holding sufficient mass to satisfy the formation criteria. In our simulations, we choose ${\rho }_{\mathrm{thres}}=3.0\times {10}^{-26}$ g cm⁻³, ${M}_{\mathrm{thres}}=3.0\times {10}^{5}{M}_{\odot }$, and ${t}_{\mathrm{dyn},\min }=10\,\mathrm{Myr}$. Additionally, the gas in the parent cell must be collapsing (determined by measuring negative velocity divergence) and Jeans unstable. If all conditions are met, the parent cell produces a star cluster particle with 10% mass efficiency so that the initial mass is given by ${M}_{\mathrm{SC}}^{\mathrm{init}}=0.1{\rho }_{g}{\rm{\Delta }}{x}^{3}$.

Over the course of 120 Myr, 25% of that mass is ejected into the cell in which the star cluster particle resides, modeling the effects of SNe II. We inject 10⁵¹ erg of energy for every $42{M}_{\odot }$, so that a star cluster particle of M = 3.0 × 10⁵M_⊙ expels a total of E_tot = 5.4 × 10⁵⁴ erg of energy. In our model, this total energy is comprised of thermal, magnetic, and CR energy so that ${E}_{\mathrm{th}}={f}_{\mathrm{th}}{E}_{\mathrm{tot}}$, ${E}_{B}={f}_{B}{E}_{\mathrm{tot}}$, and ${E}_{\mathrm{cr}}={f}_{\mathrm{cr}}{E}_{\mathrm{tot}}$, where ${f}_{\mathrm{th}}+{f}_{B}+{f}_{\mathrm{cr}}=1$. We adopt a conservative estimate of ${f}_{B}=0.01$, ${f}_{\mathrm{cr}}=0.1$ (Wefel 1987; Ellison et al. 2010) and assume 2% of the ejecta to be metals.

In our simulations, we explicitly track all ion species of hydrogen. All other species are tracked together in a metallicity variable. We calculate cooling using the GRACKLE chemistry library (Bryan et al. 2014; Smith et al. 2017), which is integrated with Enzo. It uses pre-computed tables of metal cooling rates as functions of the gas density and temperature generated by CLOUDY (Ferland et al. 2013). In addition, we consider uniform photoelectric heating of 8.5 × 10⁻²⁶ erg s⁻¹ cm⁻³ without self-shielding (Tasker & Bryan 2008). In our idealized isolated disk setup, we do not consider the ultraviolet background radiation from distant quasars and galaxies.

We use Trident (Hummels et al. 2017), a python-based tool integrated with yt (Turk et al. 2012) to construct ion densities and generate synthetic spectra from our simulated data sets. For ions that are not explicitly tracked by the simulation, we can determine their number densities, n_X, in post-processing with

$$\begin{eqnarray}&&{n}_{X}={n}_{{\rm{H}}}Z{\left(\displaystyle \frac{{n}_{X}}{{n}_{{\rm{H}}}}\right)}_{\odot },\end{eqnarray} \tag{ 17 }$$

where n_H is the total hydrogen number density, Z is the metallicity (a value that is kept track of throughout the simulation), and (n_X/n_H) is the solar number abundance for any element that is not explicitly tracked.

Trident computes relative ion abundances in a simulation cell by considering both photoionization from an extragalactic UV background (Haardt & Madau 2012) and collisional ionization within the gas. The collisional ionization rate is determined by the cell temperature, density, and metallicity, using an extensive set of lookup tables precalculated by CLOUDY, which assume ionization abundances as predicted by collisional ionization equilibrium (CIE). CIE holds even when CR pressure dominates thermal pressure because the underlying assumption of our model is that the CR fluid is collisionless. Any deviations from CIE (which is more likely in the dense regions of the galactic disk than in the CGM) would likely increase the ionization rate in low-temperature gas (Oppenheimer & Schaye 2013).

We use this functionality in Trident to plot column densities of different ions as a function of impact parameter in Figure 9. We calculate the column densities by defining sightlines through the simulation box. The ion number density is then calculated in each length element along this projected ray. The column density along the line of sight is then given by the summation ${n}_{\mathrm{col}}=\sum {dl}\cdot n$.

We simulate a suite of idealized isolated disk galaxies with initial conditions described by the AGORA collaboration (Kim et al. 2014). The fixed dark matter halo of mass M₂₀₀ = 1.074 × 10¹²M_⊙ follows the Navarro–Frenk–White (NFW; Navarro et al. 1997) profile and is situated in a hot (10⁶ K), stationary, uniform density box of (1.31 Mpc)³. The concentration and spin parameters are defined to be c = 10 and λ = 0.04, respectively. The gaseous disk has a total mass of M_d = 4.297 × 10¹⁰M_⊙ and follows the analytic exponential profile

$$\begin{eqnarray*}&&\rho (r,z)={\rho }_{0}{e}^{-r/{r}_{d}}{e}^{-| z| /{z}_{d}},\end{eqnarray*}$$

where ${r}_{d}=3.432\,\mathrm{kpc}$, z_d = 343.2 pc, and ${\rho }_{0}={M}_{d}{f}_{\mathrm{gas}}/4\pi {r}_{d}^{2}{z}_{d}$. We define ${f}_{\mathrm{gas}}=0.2$ to be the gas mass fraction of the disk. The rest of the disk mass (80%) is contained in 10⁵ stellar particles. The stellar bulge has a stellar mass of 4.297 × 10⁹M_⊙ and follows the Hernquist density profile (Hernquist 1990).

We include an initial toroidal magnetic field of strength ${B}_{0}=1\,\mu {\rm{G}}$ in the disk and a toroidal field of B₀ = 10⁻¹⁵ G in the halo. The strong initial field in the disk follows the example of Ruszkowski et al. (2017) to achieve sufficiently fast CR streaming velocities at early simulation times. The strength of the initial halo field lies in the accepted theoretical range of primordial magnetic fields, ${10}^{-20}\mbox{--}{10}^{-9}$ G (Cheng et al. 1994; Durrer & Neronov 2013). To avoid numerical interpolation errors, galaxy models that include CR physics are initialized with an isotropic background CR energy density of 0.1 erg cm⁻³. This background CR pressure is 15 orders of magnitude weaker than the initial conditions of the gas pressure and does not alter the thermal gas in any significant way. Because we are interested in tracing the cycling process of metals, we set the initial metallicities of the disk and halo to be 0.3Z_⊙ and 10⁻³ Z_⊙, respectively.

Our simulation suite is designed to isolate the effect of different implementations of CR transport. All of the galaxy models share the initial conditions described above. The fiducial models differ from each other only in their CR transport prescription and are described below:

• 1.  
  Model ncr does not include CR physics and serves as our control model.
• 2.  
  Model adv only simulates CR advection with the bulk motion of the gas.
• 3.  
  Model isod assumes isotropic CR diffusion using the algorithm described in Salem & Bryan (2014) with a constant diffusion coefficient of ${\kappa }_{\varepsilon }=3\times {10}^{28}$ cm² s⁻¹.
• 4.  
  Model anisd assumes anisotropic CR diffusion along magnetic field lines with a constant diffusion coefficient of κ_ε = 3 × 10²⁸ cm² s⁻¹.
• 5.  
  Model anisdh builds on model anisd with an added CR heating term H_c.
• 6.  
  Model stream assumes CR streaming with a streaming factor of f_s = 4 (Ruszkowski et al. 2017).

All of our fiducial models that include CR physics (in bold font in Table 1) inject 10% of their supernova ejecta in the form of CR energy. See Table 1 for a summary of the different models.

The CGM in the isolated galaxy models is enriched solely by the outflows that expel gas from the disk. For this reason, we first turn our attention to analyzing the CR-driven winds and their dependence on the CR transport mechanism.

Generally speaking, a CR-driven wind begins when CRs move down their gradient, out of the midplane of the galaxy. CR pressure support then lifts low-entropy gas out of the galactic potential well, triggering outflows. Figure 1 displays the relationship between gas entropy (top row), the vertical component of the velocity of the galactic winds (middle row), and the ratio of CR to thermal pressures (bottom row) after 2 Gyr of evolution. From left to right, the columns show models with no CR physics (ncr), CR advection (adv) only, isotropic CR diffusion (isod), and anisotropic CR diffusion (anisd, and CR streaming (stream)).

Zoom In Zoom Out Reset image size

Although models ncr and adv experienced a brief period of weak gas expulsion, by t = 2 Gyr these galaxies have lost all signs of outflows and show no signs of strong inflows. This is reinforced by the gas entropy profiles, which resemble the initial conditions. The main difference between models ncr and adv at this time is near the midplane of the galaxy, where the CR pressure of model adv has created a thicker vertical disk profile. Confined to move solely through advection, the CRs in model adv have no efficient mechanism for escaping the disk.

Models with CR transport relative to the gas (isod, anisd, and stream) all drive strong outflows with velocities reaching 10² km s⁻¹, consistent with previous works (e.g., Pakmor et al. 2016b; Salem et al. 2016; Wiener et al. 2017).

Model isod drives relatively thin and uniform conical outflows. The gas entropy profile is relatively unaltered, reaching values near 5 × 10² cm² keV near the midplane of the disk. CR pressure dominates thermal pressure, tracing the shape of the outflows out to a cylindrical radius of 10 kpc. Outside of the active outflow region, the CR pressure is below one tenth of the gas pressure.

The outflows generated with anisotropic CR diffusion in model anisd have a thinner radial profile and reach higher velocities compared to those driven by isotropic diffusion. For z > 20 kpc outside of the midplane, there are signs of infalling gas surrounding the outflow. The CR pressure dominates thermal pressure for nearly all radii in the 100 × 100 kpc projection. The gas entropy reaches values of 50 cm² keV just outside the disk, roughly an order of magnitude lower than that in model isod.

Compared to both of the diffusion models, the outflows generated by CR streaming start higher above the midplane of the disk and have a wider horizontal extent, reaching radii of nearly 40 kpc at a vertical height of 50 kpc. Near the galactic center, there are several filaments of inflowing gas. Compared to model anisd, model stream has lower gas entropy immediately outside the disk and at larger radii, maintaining values around 5 × 10² cm² keV at radii 50 kpc. Unlike either of the diffusion models, the distribution of the CR pressure ratio in model stream is patchy, ranging from 0.2 to 50 times that of the thermal pressure.

Figure 2 follows the time evolution of the density-weighted outflow velocity as a function of the height above the galactic disk for models with different CR transport prescriptions. At early times, gas in our control model, ncr, is inflowing at all radii within 200 kpc of the galactic center. After 2 Gyr, weak (${v}_{{\rm{z}}}=10$ km s⁻¹) outflows develop at heights above 30 kpc. As we show in Figure 6, these weak outflows fail to enrich the CGM with metals. Instead, we turn our attention to models isod, anisd, and stream, for which 10% of supernova feedback was injected as CR energy.

Zoom In Zoom Out Reset image size

Model isod has the strongest outflows at t = 1 Gyr, with radially averaged velocities surpassing 100 km s⁻¹. By t = 4 Gyr, these outflows weaken significantly, with average velocities hovering around 10 km s⁻¹. After 8 Gyr, notable inflows develop at radii above 30 kpc and persist for the rest of the simulation time.

This galactic wind is consistent with expectations of isotropic CR diffusion. The initial burst of supernova feedback in model isod creates a steep gradient of CR energy density. With isotropic diffusion, CRs drag gas out of the galaxy as they travel down their own gradient, uninhibited by the presence of magnetic fields. This process continues until the CR gradient flattens, slowing down the diffusion process. At later times, the star formation has decreased enough that the newly injected CRs can no longer create a gradient steep enough to drive a strong wind. As the CRs in the CGM continue to stream away from the galaxy, they cannot provide enough pressure support to the gas they expelled, and it begins to collapse back down onto the galaxy.

Model anisd retains consistent outflow velocities within 50 kpc of the disk throughout its evolution. Beyond 50 kpc, the winds are sensitive to the star formation history. For example, the double-valley feature at t = 8 Gyr traces the suppressed star formation at t = 7 Gyr (see Figure 5).

The steady nature of the galactic outflows in model anisd is fueled by anisotropic CR diffusion. In this approximation, the velocity with which CRs can escape the disk is regulated by the complex geometry of magnetic field lines. Therefore, the galaxy releases its built-up CR energy slowly over time. Although the galaxy is quenched after t = 10 Gyr, weak outflows in model anisd persist out to t = 13 Gyr.

Like in model anisd, the key to sustained outflows in model stream is its anisotropic CR transport. Unlike any other fiducial run, model stream shows signs of inflow near the disk while simultaneously hosting outflows at larger radii.

The integrated mass-loading factor (IMLF), defined as the ratio of the total injected gas mass over the total stellar mass, quantifies the expelled gas content. Figure 3 shows the evolution of the IMLF in models isod, anisd, anisdh, and stream. We exclude models nocr and adv, which did not drive outflows, from this analysis.

Zoom In Zoom Out Reset image size

For the first 1.5 Gyr, the IMLF is nearly indistinguishable between models isod, anisd, and anisdh. The low IMLF at early times in model stream is most likely due to the later onset of the galactic wind as well as to the wind location above the galactic midplane.

At later times, the IMLF decreases in model isod, yet increases in the other three models. Model anisd reaches the highest IMLF due to the accumulated reservoir of CR energy near the galactic disk that continues to drive outflows even after star formation has ceased. Models anisdh and stream have weaker IMLFs than model ansid due to higher SFRs and losses in CR energy through the heating term.

The relationship between CR transport and a galaxy's SFR can dramatically influence the galaxy morphology. The outflows generated by different CR transport prescriptions remove gas from the disk that would otherwise have been available for star formation. In addition, the presence of strong CR pressure in the disk can prevent star formation by stabilizing the thermal gas against collapse. In this section, we compare the morphologies of our fiducial galaxy models after evolving the simulations for 13 Gyr. Because a galaxy's morphology is so intricately tied to its star formation history, we describe Figures 4 and 5 simultaneously.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 displays the density-weighted face-on and edge-on projections of the gas density in our fiducial galaxy models after 13 Gyr of evolution. In the left panel, model ncr serves as the benchmark example of an Milky Way-type disk galaxy. This galaxy has a clear spiral structure, with the disk extending to a radius of roughly 25 kpc and a vertical height of roughly 3 kpc. The SFR of model ncr begins with one solar mass per year and slowly decreases over time, hovering around a few tenths of a solar mass at later times. With a modest star formation history and no galactic winds, model ncr retains much of its gas in its disk, with average density values around 3 × 10⁻²⁵ g cm⁻³.

Although model adv has had a similar outflow history as ncr, its star formation history has dramatically altered its morphology. Model adv quenched only 200 Myr after its initial burst of star formation (see Figure 5). Confined to advection-only transport, the CRs in this galaxy are trapped inside the galactic disk. After the first episode of star formation induces CR feedback, CR pressure dominates gas pressure (see the bottom-left panel of Figure 1 and the discussion in Section 4.1). This additional pressure stabilizes the gas against against collapse, halting future star formation. Without galactic outflows to carry CRs out of the disk and without any new stars to create supernovae that could potentially drive outflows, this galaxy will remain quenched. Although model adv has a substantial reservoir of gas in its disk, CR pressure expands its vertical profile, keeping the gas at lower average densities.

The star formation history of model isod most closely follows that of model ncr. Episodes of star formation that expel CRs into the disk keep the SFR of model isod consistently below that of the control. However, because isotropic diffusion is efficient at removing CRs from the galactic disk, star formation is only weakly suppressed in this galaxy. The morphology of model isod is qualitatively similar to that of model ncr out to a cylindrical radius of 15 kpc. Having expelled much of its ISM through outflows, model isod has a shorter radial and vertical extent to its gas.

Anisotropic CR transport in model anisd retains CR pressure in the disk, which suppresses star formation at early times. As galactic outflows develop, the CR pressure is lifted out of the disk, allowing for star formation to resume. After roughly 2 Gyr of evolution, the star formation in model anisd briefly matches that of the control model. The injected CR pressure following this star formation episode decreases future star formation, ultimately quenching the galaxy after 10 Gyr. CR heating in model anisdh relieves some of the CR pressure in the disk, leading to a higher SFR than that in model anisd. Although at the time of the snapshot in Figure 4 model anisd has been quenched for 3 Gyr, this galaxy model retains some spiral structure and high gas densities within a 15 kpc radius of its center. The lingering CR pressure from past star formation creates an extended low-density gas profile around the disk.

Both models anisd and stream have extended, thick gaseous profiles that are most apparent in the edge-on view. Model anisd has a bimodal distribution of gas, with a core around 2 × 10²⁵ g cm⁻³ extending to 15 kpc, surrounded by a less dense gaseous halo extending to 25 kpc. Although some spiral structure is present, it is significantly less pronounced than in model isod. Model stream has relatively low gas density in the disk with a spiral arm structure that is thinner than that of models ncr or isod.

Because of the toroidal geometry of the initial magnetic field, the CRs in model stream suppress star formation for the first 500 Myr. Once the magnetic field develops a sufficiently strong vertical component, the CRs escape the disk, dragging thermal gas along with them. The outflows relieve the ISM of the CR pressure, allowing star formation to resume. Filaments of inflowing gas near the midplane of the disk (see Figure 1) supply the additional gas necessary for extended episodes of star formation. Compared to the other galaxy models, model stream has the lowest gas densities and the most extended vertical and radial disk profile.

We now turn our attention to the different ionization structures within the CGM of the simulated galaxy models. Our control model, ncr, has limited amounts of metals outside of the galactic disk. This is an unrealistic result, and we will not dwell on its implications here. Our goal is to isolate the effects of CR transport on the CGM temperature and ionization structure, and so the following figures and discussions focus on models isod, anisd, and stream. Where relevant, we discuss model anisdh, which includes the heating term H_c, which is present in the CR streaming prescription. Although the heating term is an unphysical addition to CR diffusion, its presence helps isolate which aspect of the CR streaming or diffusion mechanism is responsible for observed differences in temperature and ionization states (Wiener et al. 2017).

Figure 6 shows the density-weighted projections of gas metallicity and temperature and the unweighted projections of H i and O vi column densities after 13 Gyr of evolution. Each column holds a different mode of CR transport. Starting from the left, the columns show results for galaxy models ncr, isod, anisd, anisdh, and stream.

Zoom In Zoom Out Reset image size

The outflows in models with CR diffusion or streaming have populated their CGM with metals from the disk out to radii surpassing 200 kpc. Model isod has a relatively uniform distribution of metallicity in its CGM, with a density-weighted average value around 0.2 Z_⊙. Similarly, its temperature and H i and O vi column densities are also spatially uniform. CR pressure support creates a cooler temperature profile (roughly 3 × 10⁵ K) compared to the control. The cooler temperatures result in stronger column densities of H i and O vi compared to model ncr.

Anisotropic diffusion in model anisd produces a steep temperature and H i column density gradient. The temperature of the gas is coolest near the disk, where CR pressure is strongest, and decreases radially outwards. The H i column density follows the shape of the temperature profile, extending out to where the density-weighted temperature is 10⁶ K. The O vi column density is sensitive to metallicity, so the column density profile traces the radial extent of the outflows. Compared to the other fiducial runs, model anisd has the strongest column densities of H i and O vi. We point out that although this model has been quenched for 3 Gyr, the gradual release of CR pressure results in the accumulation of cold gas just outside the disk.

The metallicity distribution of model anisdh is qualitatively similar to that of model anisd. However, the added heating term dramatically decreases the impact of CR pressure, resulting in warmer temperatures and weaker column densities of H i and O vi than those in model anisd.

The outflows in model stream reach larger radii than those in model anisd. The temperature profile is dominated by clumps of cool and warm gas extending out to radii of 150 kpc. Model stream has a patchy distribution of both H i and O vi column densities, with both coexisting in relative abundance at radii above 100 kpc from the galactic center.

To better understand the CGM structure, we analyze the abundance of O vi as a function of temperature and density in Figure 7. Since we are primarily interested in the gas in the CGM, we exclude data points contained within a cylindrical radius of 25 kpc and a vertical height of 3 kpc of the galactic center.

Zoom In Zoom Out Reset image size

The CGM of models ncr and adv has lower masses of O vi and smaller amounts of low-density gas compared to the other models. Since low-density values trace large radii in the CGM, the artificially high density floors are due to the lack of outflows. Suppressed star formation in model adv is responsible for the low metallicities in its disk that ultimately result in low masses of O vi.

Models with additional CR transport drive strong outflows that enrich the CGM with metals (see Figure 6). This enrichment creates higher column densities of O vi at large radii. However, the different prescriptions for CR transport result in varied phase structures within the CGM.

One such difference is the shape of the temperature-density phase diagram. Models with weak CR pressure support (models isod and anisdh) have less cool gas at low densities than models anisd and stream. The strong CR pressure support in models anisd and stream supports a wide range of temperatures at each density. This feature is most pronounced at low densities in model stream.

CR transport also affects the abundance of O vi and the temperature of the gas that produces it. Both models stream and anisd predict a reservoir of O vi ionized with the help of CR pressure support in low-density gas. However, only model stream predicts an abundance of O vi photoionized with gas temperatures around 3 × 10⁵ K, consistent with predictions in McQuinn & Werk (2018).

Figure 8 shows the total enclosed mass as a function of spherical radius for six different galaxy models. The colors of the bars denote the fractional distribution of the gas temperature at different radial shells. We consider three temperature bins of gas: cold gas ($T\lt {10}^{5}\,{\rm{K}}$) in dark purple, warm gas (${10}^{5}\,{\rm{K}}\lt T\lt {10}^{6}\,{\rm{K}}$) in medium purple, and hot gas ($T\gt {10}^{6}\,{\rm{K}}$) in light purple. The colored gas fraction is not cumulative since the cold mass contribution from the disk would dominate the distribution in some cases.

Zoom In Zoom Out Reset image size

Having driven weak outflows, model ncr has trace amounts of cold or warm gas outside of the disk. In contrast, although model adv had similarly weak winds, the presence of CR pressure keeps the gas immediately around the disk below 10⁶ K. The influence of CR pressure drops abruptly after a radius of 50 kpc. Because the CR pressure in the disk for this model stopped star formation almost immediately, model adv has more gas available both inside and outside the disk of this model.

Models in which CRs can move relative to the gas (via diffusion or streaming) have an altered temperature structure in their CGM. Model isod is dominated by warm gas out to large radii. Since isotropic diffusion depends solely on the direction and magnitude of the CR energy density gradient, the CR pressure distribution is nearly uniform throughout the halo. Over time, with more stellar feedback releasing CRs into the disk, the CR pressure accumulates in the halo, providing pressure support to the thermal gas.

Models with anisotropic CR transport exhibit a multiphase temperature structure. Models anisd and anisdh are both dominated by cold gas near the disk. In both models, the cold gas extends out to radii of 125 kpc and the warm gas extends to radii of 200 kpc.

Comparing this result with Figure 6, we see that the cool gas, traced by H i column densities, decreases radially away from the center. The added heating term in model anisdh converts some CR pressure, which is responsible for supporting the warm and cold gas, into heating the gas, consistent with results in Wiener et al. (2017). Even with the heating term, the distribution of warm gas in model anisdh is still similar to that of anisd.

Model stream shows signs of a true multiphase medium. This model has the lowest cumulative gas mass in its disk, likely due to its low gas densities and recent episodes of star formation (see Figure 5). Although warm gas dominates its CGM out to radii of 150 kpc, cold gas survives in abundance 200 kpc away from the galactic center. Using Figure 6, we see that the temperature structure in model stream does indeed result in a patchy, multiphase medium.

Figure 9 holds the column densities of H i, Si iv, C iii, and O vi as a function of the spherical radius from the center of the galaxy. Each scatter point in the plot represents one column density measurement at that radius. In each panel, we include points from both the fiducial model and its low-resolution counterpart. The solid black lines show the average column densities for the fiducial run, while the dashed black lines show the average column densities for the lower resolution run with the same physics. For details on how the column densities were constructed, see Section 3.3.

Zoom In Zoom Out Reset image size

To sufficiently sample our simulation space, we calculate the column densities at randomly oriented sightlines passing through the CGM. To generate a sightline, we first select a random point on the sphere of radius r ∈ [10, 200] kpc from the galaxy center. We then chose a sightline by selecting a random angle in the plane tangent to the sphere at that point.

The control model, ncr, underpredicts the column densities of all four ions. This is both due to the lack of metals in its CGM and a deficit of cool gas.

Of the remaining models, model isod has the lowest column densities with a flat distribution across the impact parameter. This picture is consistent with weak and spatially uniform CR pressure in the CGM. Model anisd has higher column densities with a wider spread at large radii. This is likely due to the uneven metallicity distribution (see Figure 6). Model stream predicts column densities that are similar in strength to model anisd. However, the column density measurements in model stream are more tightly clustered around the average value, which is qualitatively similar to observations.

For all models, a lower resolution predicts higher column densities, possibly because the complicated structure of the magnetic field lines is not fully resolved. The column densities in model stream are the least sensitive to changes in resolution.

The column densities presented here are a qualitative example of the differences in ionization structure between different models of CR transport. In reality, the exact values of the ion column densities would be influenced by the metallicity and inflows from the IGM. Therefore, in order to better match observations, we would need to simulate galaxies in a cosmological context.

CR pressure in the CGM impacts the temperature and ionization structure of the gas. For this reason, we investigate the role of CR transport mechanisms in altering this CR pressure distribution.

In Figure 10 we explore the distribution of CR pressure as a function of spherical radius for models isod, anisd, anisdh, and stream. The pixels in the phase plot are colored by the density of the gas. Like in Figure 7, we cut out the disk from our data sample.

Zoom In Zoom Out Reset image size

Model isod has a flat and narrow profile of CR pressure as a function of radius, with the ratio ${P}_{{\rm{c}}}/{P}_{{\rm{g}}}$ ranging between 0.5 and 5. At larger radii, the CR pressure ratio widens considerably, revealing regions where CR pressure is sparse.

Model anisd has a much wider range of CR pressure ratios across all radii. Compared to isotropic diffusion, anisotropic diffusion creates CR pressure ratios that are nearly an order of magnitude higher across all radii. In this model, the CR pressure dominates thermal pressure near the galactic disk and at high gas densities. The heating term in model anisdh lowers the CR pressure ratio at all radii. However, this model still retains the same qualitative shape of the CR pressure distribution that is present in model anisd.

Although model stream also has a wide spread of CR to thermal gas pressure ratios, the distribution of that spread is qualitatively different from that in model anisd. In this model, CR pressure in and near the disk is roughly 1.5 dex lower, spanning nearly 2 dex in strength. Unlike in model anisd, the CR pressure is low near the galactic disk and dominates thermal pressure at large radii (and low gas densities).

Figure 11 shows the density-averaged ratio of CR to thermal pressure as a function of impact parameter, measured after 3 Gyr of evolution. The solid black line represents the low-resolution runs of our fiducial models. All other lines differ from this run only in the parameter specified in the label. Blue and green lines represent galaxy models in which supernova feedback injected 1% and 30% of its energy as CRs, respectively. Purple lines represent slower CR transport velocities with either a diffusion coefficient of ${\kappa }_{\varepsilon }={10}^{27}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ in models isod and anisd, or a streaming factor of f_s = 1.0 in model stream. To better compare models with anisotropic diffusion to those with isotropic diffusion, we explore faster CR transport (${\kappa }_{\varepsilon }={10}^{29}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$) in the middle panel, depicted by the red line. Finally, the black dotted line is the high-resolution fiducial run. For a summary of the parameters described above, see Table 1.

Zoom In Zoom Out Reset image size

Varying the injected CR fraction does not correspond with a linear change in the ratio of CR to thermal pressure in the CGM. For example, increasing f_CR to 0.3 only marginally increases CR pressure in model isod. Counterintuitively, increasing f_CR in models anisd and stream actually decreases the CR pressure ratio in the CGM by suppressing star formation in the disk and thus preventing the injection of additional CRs. On the flip side, decreasing the injected CR fraction by a factor of 10 from the fiducial value does decrease the CR pressure in model isod by roughly an order of magnitude at radii above 100 kpc. However, this decrease in f_CR does not significantly alter the CR pressure ratio for models anisd and stream.

The CR pressure distribution in all three models depends strongly on their transport velocities. Decreasing the diffusion coefficient by a factor of 10 limits the reach of outflows to 100 kpc in model isod and 60 kpc in model anisd. Similarly, decreasing the streaming factor by four limits the reach of CR pressure to 40 kpc in 3 Gyr. In addition to limiting the radial extent of CR pressure, the decrease in transport velocity in model isod increases the CR pressure ratio by a factor of 10 at small radii. Increasing the transport velocity in model anisd decreases the CR pressure within 75 kpc of the galactic center, but has no discernible effects at larger radii.

Doubling the resolution decreases the CR pressure ratio by roughly an order of magnitude for models isod and anisd. In model stream, the CR pressure ratio of the higher resolution run only deviates from the fiducial run above 150 kpc.

Model isod is the most sensitive to changing parameters. In this approximation of CR transport, CRs can only move relative to the gas by diffusing down their own gradient. This motion is very sensitive to the injection fraction of CRs, which defines the strength of the CR gradient. Anisotropic diffusion is moderated by transport along magnetic field lines. The complicated geometry of the magnetic field lines slows anisotropic diffusion near the disk, making model anisd less sensitive to the CR injection fraction and star formation history than model isod. The most robust model is stream, in which the CR transport depends on the CR gradient, magnetic field strength, and Alfvén wave velocity. The CR pressure ratio in all three models decreases with slower CR transport and increased resolution. Breaking this degeneracy will require additional parameter studies.

In this work, we explored the qualitative impact of the choice of CR transport on the simulated CGM. However, before simulations with CR physics can hold predictive power, several improvements must be made to constrain the details of this transport. We discuss these factors in greater detail below.

5.1.1. Magnetic Fields

5.1.2. Improved CR Physics

5.1.3. Cosmological Context

The Sod shock-tube, first described in Sod (1978), is used to test the behavior of gas in the presence of strong shocks in numerical simulations. Because the original Sod shock-tube does not include CRs, we use a modified version first described in Pfrommer et al. (2006). The initial conditions are described in Table 2 below.

Table 2.  Initial Conditions for the Modified Sod Shock-tube

	ρ	Pg	εc	v
Left	1.0	$2/3\times {10}^{5}$	4.0 × 105	0.0a
Right	0.2	267.2	801.6	0.0

Note.

^aThe initial density, thermal pressure, CR energy density, and velocity for the modified Sod shock-tube. The left and right regions of the shock-tube are defined as $(0\lt x\lt 250\,\mathrm{cm})$ and (250 < x < 500 cm), respectively.

Download table as:  ASCIITypeset image

In Figure 12 we test our implementation of the non-diffusive (κ_ε = 0.0) two-fluid CR model and compare it to the existing implementation in the Zeus hydro scheme in Enzo (Salem & Bryan 2014). From left to right, the panels show the density, velocity, and pressures of this shock-tube after t = 0.31 code units of evolution. Where applicable, the analytic solution is depicted as a black dashed line. Compared to Salem & Bryan (2014), our implementation has a sharper velocity profile at x ≃ 250 cm, but slightly underpredicts the density value around x = 400 cm. Overall, the two implementations of CR advection agree well with each other and with the analytic solution.

Zoom In Zoom Out Reset image size

The Brio-Wu shock-tube (Brio & Wu 1988) tests the advection properties of gas interactions with magnetic fields in the presence of an extreme shock. In the modified Brio-Wu shock-tube, we include CR pressure and evolve the system for 0.2 code units. The initial conditions of the gas are described in Table 3. In Figure 13, we show the time evolution of the Brio-Wu shock-tube using our CR implementation.

Table 3.  Modified Brio-Wu Shock-tube

	ρ	Pg	${\varepsilon }_{c}$	${v}_{x}$	${v}_{y}$	Bx	By
Left	1.0	0.6	1.2	0.0	0.0	1.0	1.0a
Right	0.125	0.06	0.12	0.0	0.0	1.0	−1.0

Note.

^aThe initial density, thermal pressure, CR energy density, velocity and magnetic field strengths for the modified Brio-Wu shock-tube. The left and right regions of the shock-tube are defined by $(-0.5\lt x\lt 0)$ and $(0\lt x\lt 0.5)$, respectively.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

In the anisotropic diffusion approximation, CRs propagate down their gradient, along magnetic field lines. In Figure 14 we test our implementation with the anisotropic ring problem initially described by Parrish & Stone (2005) and Sharma & Hammett (2007) and explored in detail in Pakmor et al. (2016a). The experiment uses a uniform density box in with no gas (or a gas that is fixed in space such that there is no CR advection) in a domain of [−1, 1]².

Zoom In Zoom Out Reset image size

Following the setup in Pakmor et al. (2016a), we set initial conditions for the CR energy density to be

$$\begin{eqnarray}{\epsilon }_{c}(x,y)=\left\{\begin{array}{ll}12 & \mathrm{if}\,0.5\lt r\lt 0.7\,\mathrm{and}\,| \phi | \lt \displaystyle \frac{\pi }{12}\\ 10 & \mathrm{else},\end{array}\right.\end{eqnarray} \tag{ 18 }$$

where the radial coordinate $r=\sqrt{{x}^{2}+{y}^{2}}$, $\phi =\mathrm{atan}2(y,x)$, and the diffusion coefficient is set as κ_εκ_c = 0.01.

In Figure 14 we compare the performance of our anisotropic diffusion implementation against the analytic solution given by

$$\begin{eqnarray}&&\begin{array}{l}{\epsilon }_{c}(x,y)\\ =\left\{10+\mathrm{erfc}\left[\left(\phi +\displaystyle \frac{\pi }{12}\right)\displaystyle \frac{r}{D}\right]-\mathrm{erfc}\left[\left(\phi -\displaystyle \frac{\pi }{12}\right)\displaystyle \frac{r}{D}\right]\right.,\end{array}\end{eqnarray} \tag{ 19 }$$

where $D=\sqrt{4{\kappa }_{\epsilon }t}$. ${\epsilon }_{c}(x,y)=10$ everywhere else.

After evolving this simulation for 10 code time units, we see the initial wedge of uniform CR energy density moving down its gradient around the magnetic fields lines. The lower resolution runs have lower peak values of CR energy density and show signs of CR transport perpendicular to the magnetic field direction. Both of these properties improve with increased resolution. The best-resolved simulation (800 × 800, fixed grid) matches the analytic solution well.

At a later time (t = 100 code units), the CR energy density has reached the opposite side of the circle traced by the toroidal magnetic field lines. The energy density approaches a steady state and is nearly evenly distributed around the annulus. In low-resolution runs, the average final CR energy is lower due to conservation of CR energy in the thicker rings.

We test the behavior of our CR streaming implementation with a 1D simulation of an initial Gaussian profile of CR energy density (see Figure 15). The initial CR energy density profile is set by

$$\begin{eqnarray}&&\varepsilon ={\varepsilon }_{0}{e}^{-{x}^{2}/2D},\end{eqnarray} \tag{ 20 }$$

where x is the spatial coordinate. In our example, we chose constants ε₀ = 100 and D = 0.05. We include a magnetic field in the $\hat{x}$ direction of a strength of 0.25 in code units. We isolate the effects of CR streaming and diffusion by fixing the gas so that no advection is taking place.

Zoom In Zoom Out Reset image size

For CR diffusion, the analytic solution for the evolution of the CR energy density over time is given by

$$\begin{eqnarray}&&\varepsilon ={\varepsilon }_{0}\sqrt{\displaystyle \frac{D}{D+2{\kappa }_{\varepsilon }t}}\exp \left(\displaystyle \frac{-{x}^{2}}{2(D+2{\kappa }_{\varepsilon }t)}\right).\end{eqnarray} \tag{ 21 }$$

The CR diffusion implementation follows the analytic solution well. Although there is no analytic solution for the CR streaming case, our results are qualitatively comparable to previous studies (e.g., Uhlig et al. 2012; Wiener et al. 2017).