Cosmic Ray Acceleration of Cool Clouds in the Circumgalactic Medium

TL;DR

This study explores how cosmic rays can accelerate cool clouds in the circumgalactic medium through magnetohydrodynamic simulations, revealing that CR pressure buildup can propel clouds to high velocities, with implications for understanding high velocity clouds.

Contribution

First 2D MHD simulations of cosmic ray bottlenecks in the CGM, demonstrating cloud acceleration and the effects of radiative cooling on cloud stability.

Findings

CR pressure buildup can accelerate cool clouds to hundreds of km/s.

Radiative cooling helps keep clouds intact against CR wave heating.

Cloud acceleration scales with CR flux, roughly linear at low fluxes.

Abstract

We investigate a mechanism for accelerating cool (10$^4$ K) clouds in the circumgalactic medium (CGM) with cosmic rays (CRs), possibly explaining some characteristics of observed high velocity clouds (HVCs). Enforcing CRs to stream down their pressure gradient into a region of slow streaming speed results in significant buildup of CR pressure which can accelerate the CGM. We present the results of the first two-dimensional magnetohydrodynamic (MHD) simulations of such `CR bottlenecks,' expanding on simpler simulations in 1D from \cite{wiener17a}. Although much more investigation is required, we find two main results. First, radiative cooling in the interfaces of these clouds is sufficient to keep the cloud intact to CR wave heating. Second, cloud acceleration depends almost linearly with the injected CR flux at low values (comparable to that expected from a Milky Way-like star formation rate), but scales sublinearly at higher CR fluxes in 1D simulations. 2D simulations show hints of sublinear dependence at high CR fluxes but are consistent with pure linear dependence up to the CR fluxes tested. It is therefore plausible to accelerate cool clouds in the CGM to speeds of hundreds of km s$^{-1}$.

Figures

Equations

$${\mathcal{H}}_c=-{\mathbf{v}}_{\mathbf{A}}\cdot \nabla P_c=v_A|\hat{\mathbf{b}}\cdot \nabla P_c|.$$

$$\hat{\mathbf{b}}\cdot \nabla P_c{\rho }^{-{\gamma }_c/2}=0.$$

$$P_c=\text{constant}.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot \left(\rho \mathbf{uu}-\frac{\mathbf{BB}}{4\pi }\right)+\nabla P_{\text{tot}}=0$$

$$\frac{\partial \mathbf{B}}{\partial t}-\nabla \times (\mathbf{u}\times \mathbf{B})=0$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[(E+P_{\text{tot}})\mathbf{u}-\frac{\mathbf{B}(\mathbf{B}\cdot \mathbf{u})}{4\pi }\right]=-\nabla \cdot {\mathbf{F}}_c-\mathcal{C}$$

$$\frac{\partial E_c}{\partial t}+\nabla \cdot (E_c\mathbf{u})=-P_c\nabla \cdot \mathbf{u}-{\mathcal{H}}_c-\nabla \cdot {\mathbf{F}}_c+q_c$$

$$P_{\text{tot}}=P_g+P_c+\frac{B^2}{8\pi }$$

$$E=E_g+E_c+\frac{1}{2}\rho u^2+\frac{B^2}{8\pi }$$

$${\mathbf{F}}_c=(E_c+P_c){\mathbf{v}}_s$$

$$v_s=v_A\tanh \left(\frac{|\hat{\mathbf{b}}\cdot \nabla P_c|}{P_c/L}\right)$$

$$\rho (r)={\rho }_{\text{hot}}+({\rho }_{\text{cold}}-{\rho }_{\text{hot}})\left(\frac{1}{2}-\frac{1}{2}\tanh \left(\frac{r-r_c}{t_c}\right)\right).$$

$$q_c(t)=q_0\left(1-e^{-t/t_{\text{ramp}}}\right)\left(1-e^{(t-t_{\text{end}})/t_{\text{ramp}}}\right)$$

$$\mathcal{C}=n_H^2\Lambda (T)-n_H\Gamma$$

$$\Lambda (T)=1.1\times 10^{-21}\ast 10^{\Theta (\log (T/10^5\text{K}))}{\text{ erg cm}}^3{\text{ s}}^{-1}$$

$$\Theta (x)=0.4x-3+\frac{5.2}{e^{x+0.08}+e^{-1.5(x+0.08)}}.$$

$$\begin{gathered} I_{cr}=SF{R}_{\ast }{R}_{\text{SN}}(\xi 10^{51}{\text{ erg SN}}^{-1}) \\ =\left(\frac{SF{R}_{\ast }}{1M_{\odot }{\text{ yr}}^{-1}}\right)\left(\frac{R_{\text{SN}}}{1\text{ SN}/200M_{\odot }}\right)5\times 10^{47}{\text{ erg yr}}^{-1} \end{gathered}$$

$$\begin{gathered} {\Phi }_{cr}=\frac{I_{cr}}{2\pi r_{\ast }^2}=3\times 10^{-6}{\text{ erg cm}}^{-2}{\text{s}}^{-1} \\ \times \left(\frac{SF{R}_{\ast }}{1M_{\odot }{\text{yr}}^{-1}}\right)\left(\frac{R_{\text{SN}}}{1\text{ SN}/200M_{\odot }}\right){\left(\frac{r_{\ast }}{10\text{ kpc}}\right)}^{-2} \end{gathered}$$

$$\Delta tto0.0pt\text{lower}2.5pt\sim \text{hss}\text{raise}1.5pt<\frac{\Delta x}{v_s}\frac{\Delta x}{L}$$

$$\Delta tto0.0pt\text{lower}2.5pt\sim \text{hss}\text{raise}1.5pt<\frac{\Delta x}{v_s}$$

Full text

1 Introduction

The circumgalactic medium (CGM), a halo of gas surrounding the Galactic disk, is an important component of any model of galaxy evolution. The CGM contains most of the baryons in the Galaxy and serves as a potential reservoir of fresh material to form new stars. At the same time, the chemical makeup of the CGM holds imprints of the Galaxy’s star formation history, as supernova eject material into the CGM in large-scale Galactic fountains or outflows. The makeup and structure of the CGM is thus of great interest (Bregman et al., 2018; Zhang, 2018).

One of many probes into the CGM is absorption line spectroscopy of bright background (extraglactic) sources. Wakker et al. (2012) analyze dozens of sight lines through the Galactic halo and use absorption line features to measure the column densities of various ions and their velocity structure. Among other findings, they found large column densities of OVI, suggesting the existence of “transition temperature” gas at 105.5 K. Additionally, Doppler shifts in the absorption lines suggest bulk velocities of up to 100 km s*-1* for some components. The multiphase nature of the CGM of other galaxies is revealed in observations of sight-lines through galaxy halos, such as the COS-Halos survey. Werk et al. (2014) examine low-ionization state absorption lines to characterize the cool 104 K component of the CGM and find that the electron density of this cool phase is 100 times lower than what is expected from thermal equilibrium with the hot 106 K phase. Werk et al. (2016) study the high-ionization state component using OVI absorption lines and find that while many OVI absorption features are cospatial with low-ionization absorption, some are not, suggesting the presence of a 105.5 K phase. Low-ionization features are also observed in galactic outflows (see Chisholm et al. (2016)), suggesting that galactic winds also contain multiple gas phases.

From these and other observations, astronomers infer the presence of cool or warm (104 K) clouds moving very quickly within the otherwise hot (106 K) CGM, with thick intermediate temperature (105.5 K) interfaces. These are referred to as high velocity clouds (HVCs). The properties of HVCs introduce some puzzling questions about the dynamics of the CGM. Were they accelerated along with the launching of a Galactic wind, or did they condense in situ out of the hot phase of a previously launched wind? We consider here only the former possibility, addressing the issue that Galactic wind models using thermal pressure alone have trouble accelerating cool clouds to the 100 km s*-1* speeds measured without destroying them (for example, see Scannapieco & Brüggen (2015)). Additionally, few models are able to generate high OVI column densities, or even the right ion ratios (Wakker et al., 2012).

In Wiener et al. (2017) (hereafter W17) we investigated the influence of cosmic ray (CR) transport physics in this environment. Using simple 1D models, we showed that due to a bottleneck effect arising from the spatially-varying CR streaming speed, large CR pressure gradients build up as CRs exit a cool cloud. This gradient has dynamic effects, accelerating the cloud to high velocities consistent with observations, as well as energetic effects, heating the interface at the location of the gradient in such a way as to nearly reproduce observed ion ratios.

These bottlenecks may also significantly impact the spatial distribution of CR pressure in a galactic wind, critically affecting the nature of CR-driven winds. CR-driven winds are advantageous over thermally-driven winds since CRs are less subject to radiative losses and maintain their energy over longer distances. Generally speaking, CR-driven wind models involve weak CR pressure gradients pushing on the gas gently over these long distances, resulting in wind velocities that increase with height above the disk. However, as we found in W17, CR bottlenecks redistribute the CR pressure in space, resulting in sharp CR gradients at one edge of the cloud, and eliminating CR gradients upwind of this point. As such the nature of a CR driven wind through a medium with many such bottlenecks could be radically different from one without them.

In this work we expand upon these simple 1D hydrodynamic (HD) simulations with a set of 2D magnetohydrodynamic (MHD) simulations using the FLASH code, modified to include CR streaming physics (Ruszkowski et al. (2017)). Unlike in W17, these simulations handle radiative cooling in a self-consistent way, and account for the potential bending of magnetic field lines as the system evolves. This is important in regards to the bottleneck effect, since the amount of magnetic flux penetrating the cloud controls how many CRs go into the bottleneck. A magnetically-isolated cloud, for example, would have no CR bottleneck at all. By performing simulations in 2D (and possibly in future work, 3D) we can examine the impact of these lateral effects.

In §2 we briefly review the CR streaming physics that result in the bottleneck effect. In §3 we describe the FLASH code and the simulation setups. In §4 we present the results of our simulations and discuss their implications, focusing on the difference with our simple 1D simulations. We conclude in §5. Effects of our boundary conditions and resolution tests are presented in Appendices A and B respectively.

2 Cosmic Ray Streaming

All the scales in this problem are much larger than the scattering mean free path of individual CRs, so we will describe the CRs in a simple fluid approximation. The CRs are treated as a single fluid with pressure $P_{c}$ and adiabatic index $\gamma_{c}$, assumed to be 4/3. We evolve this pressure in time along with a background gas described by density $\rho$, velocity $u$, and internal energy $E_{g}$, as well as a magnetic field $B$.

We use the self-confinement theory of CR transport, wherein the bulk flow of CRs is regulated by the streaming instability (Kulsrud & Pearce (1969); Skilling (1971); Kulsrud (2005); Zweibel (2017)). In this framework, Alfvén waves are amplified by gyroresonant cosmic rays wherever the bulk CR flow speed along the magnetic field $v_{s}$ exceeds the local Alfvén speed $v_{A}$. These waves then scatter CRs in pitch angle until $v_{s}$ falls below $v_{A}$ and the instability shuts off. Moreover, in the limit of short mean free path, the streaming speed is proportional to the CR pressure gradient. The CR streaming speed $v_{s}$ is thus determined by the marginal stability condition for the streaming instability, which, in the case where Alfvén wave damping processes are negligible, is simply described by $v_{s}=v_{A}$ with the constraint that the direction of net CR flow be along the local magnetic field, and down the CR pressure gradient. We then say that CRs are ‘locked to the wave frame’ wherever there is a non-zero CR pressure gradient along the local magnetic field.

One consequence of the self-confinement picture is an energy exchange between CRs and the thermal plasma, mediated by the Alfvén waves. One can show that streaming CRs will heat the gas at a rate equal to

[TABLE]

This was first done by Wentzel (1971). Here, $\mathbf{\hat{b}}$ is a unit vector pointing along the magnetic field. This heating has potential consequences on the energetics of a system, especially at the locations of sharp CR pressure gradients. Such gradients may develop because of a second consequence of the self-confinement picture, the CR bottleneck. For a detailed description of the CR bottleneck effect, which was first predicted by Skilling (1971), we refer the reader to W17, but we summarize the main points here. Consider a one-dimensional flow of CRs along a single magnetic flux tube. As $B$ or $\rho$ change along the tube, the streaming speed $v_{A}$ changes. If $v_{A}$ increases along the direction of the CR flow, the CR pressure $P_{c}$ must drop according to

[TABLE]

If $\rho$ increases along the direction of the CR flow, we might expect $P_{c}$ to increase accordingly. But this would violate the constraint that CRs only stream down their gradient. Therefore, the only steady state solution in this case is

[TABLE]

A maximum in $\rho$ along a flux tube therefore produces a bottleneck. In the steady state, $P_{c}$ ‘upstream’ of the maximum must be constant. ‘Downstream’ of the maximum $P_{c}$ follows (2). This effect will be present in cool clouds in the CGM provided the magnetic field lines penetrate the cloud. We would then expect $\rho$ to have a maximum in the center of the cloud, causing such a CR bottleneck, and building up a strong pressure gradient on one end of the cloud.

This picture is subject to a number of caveats. We do not include any wave damping in our analysis or simulations, but inside the cloud it may be cold enough for Alfvén waves to suffer significant ion-neutral damping. At late times we expect the bottleneck to erase CR pressure gradients in this region (we will see this in the simulations in §4). Waves are not excited in the first place, so the presence or absence of ion-neutral damping has no effect on the dynamics at this point. However, at early times ion-neutral damping may have a transient effect by allowing CRs to “fill up” the cloud region more rapidly. More importantly, bottlenecks can be completely erased if magnetic field lines do not penetrate the cloud at all. The field line topology of a cool cloud will depend strongly on how it is formed. We refer the reader to W17 for further discussion of these issues.

3 Simulation Setup

3.1 Basic Equations

We use the magnetohydrodynamics (MHD) code FLASH in 2D to solve the ideal MHD equations with modifications to include CRs and CR streaming, according to Ruszkowski et al. (2017) as well as thermal gas heating and cooling terms. The equations we solve are as in Ruszkowski et al. (2017), but without gravity ($\mathbf{g}=0$), explicit CR diffusion ($\mathbf{\kappa}=0$), or supernova feedback ($\dot{p}_{\textrm{SN}}=\mathcal{H}_{\textrm{SN}}=0$):

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the above $\rho$ and $\mathbf{u}$ are the gas mass density and velocity, $\mathbf{B}$ is the magnetic field, and $P_{\textrm{tot}}$ is the total pressure

[TABLE]

consisting of thermal pressure $P_{g}=\rho k_{B}T/(\mu m_{p})$, CR pressure $P_{c}$, and magnetic pressure. The total energy density $E$ is

[TABLE]

and includes the thermal energy density $E_{g}=P_{g}/(\gamma_{g}-1)$, the CR energy density $E_{c}=P_{c}/(\gamma_{c}-1)$, the kinetic energy density and the magnetic energy density. The CR flux $\mathbf{F}_{c}$ is written in our model as a simple advective flux

[TABLE]

where the streaming speed $\mathbf{v}_{s}$ has magnitude up to $v_{A}$ (see below) and points along the local magnetic field. A CR source term $q_{c}$ is included, which we use to add CRs to the domain (see §3.2). Lastly, $\mathcal{C}$ represents radiative cooling and a heating function which maintains thermal equilibrium in the absence of cosmic rays, discussed below, and $\mathcal{H}_{c}$ is CR heating, given by (1).

As a brief aside, let us elaborate on the CR streaming velocity $\mathbf{v_{s}}$, both in the self-confinement model and in the FLASH code implementation. In the model, which neglects wave damping, CRs will stream along the background magnetic field at the Alfvén speed $v_{A}$ anywhere there is a CR pressure gradient along the field to drive the flow. If the CR pressure is constant along the field, $\mathbf{v_{s}}$ will still be along the field, but will have magnitude between 0 and $v_{A}$, determined by other factors such as a steady-state requirement. For instance, as discussed in §2, the steady state solution upstream of a bottleneck is $P_{c}=$constant, $|v_{s}|=v_{A,\textrm{min}}<v_{A,\textrm{local}}$.

In simulation, the requirement that CRs only flow down their gradient presents a computational challenge which in the CR module in FLASH (see Ruszkowski et al. (2017)) is handled with a regularization method described in Sharma et al. (2010). In this method a scale length $L$ is chosen, and the streaming speed is determined by

[TABLE]

The length scale sets a characteristic CR pressure gradient $P_{c}/L$ - gradients much steeper than this will have $v_{s}=v_{A}$, while smaller gradients will have $v_{s}<v_{A}$. In this way the streaming speed is consistently calculated everywhere, and it exhibits the behavior required by the self-confinement model: upstream of a bottleneck, CR pressure builds up, approaching a flat profile. The reduced CR gradient in this region results in the streaming speed $v_{s}$ dropping according to (12) above. This emulates the theoretical steady state picture in an approximate way which depends on the chosen length scale $L$. The larger $L$ is, the more accurate the simulation, but the more restrictive the time step. There are other ways to treat CR streaming - Jiang & Oh (2018) introduce a two-moment method where the CR flux $\mathbf{F}_{c}$ is treated as an independent variable with its own evolution equation. We briefly discuss this in §4.6.

3.2 Specifics of Model

The background plasma is initialized at uniform pressure $P_{g}=3.23\times 10^{-13}$ dyne cm*-2* and mass density $\rho_{\textrm{hot}}=2.1\times 10^{-27}$ g cm*-3*. We assume hydrogen and helium are fully ionized. We use an average mass per particle of $\mu=0.60m_{p}$, and average mass per Hydrogen nucleus of $\mu_{H}=1.43m_{p}$. This mass density and thermal pressure therefore correspond to a Hydrogen number density of $n_{H}=9.0\times 10^{-4}$ cm*-3* and a temperature $T=1.1\times 10^{6}$ K. In this background we insert a circular, fully-ionized, cooler cloud at the same pressure but with mass density $\rho_{\textrm{cold}}=2.35\times 10^{-25}$ g cm*-3*, corresponding to a density $n_{H}=0.10$ cm*-3* and a temperature $T=9.9\times 10^{3}$ K111Please note that while we quote three significant figures for some quantities in the spirit of transparency, the derived quantities are only held to two significant figures.. The density follows a tanh profile of radius $r_{c}=50$ pc and interface thickness $t_{c}=0.25$ pc according to

[TABLE]

This initial density profile is a 2D generalization of the one used in W17. The cloud is centered on coordinate $(x,y)=(1\ \textrm{kpc},0)$. The magnetic field is initialized to a uniform value of 1 $\mu$G throughout the domain, pointing in the (horizontal) $x$-direction. With these parameters, $v_{A}=62$ km s*-1* in the hot phase. The Alfvén travel time from the left boundary to the cloud is $\sim 15$ Myr, while the adiabatic acoustic travel time is $\sim 6$ Myr. These timescales are important for understanding the evolution of the system.

The initial CR distribution is set to a very low, uniform pressure. Each time step, CRs are injected into all the leftmost active cells of the domain via the source term $q_{c}$. The CR flux through the left boundary is forced to zero to ensure that all injected CRs stay in the domain and flow towards the cloud (this also prevents CRs from advecting into the domain with the gas, so that the amount of CR energy added to the domain is never more than what is added via the source term). The injection is time-dependent, emulating a sudden burst of star formation, according to

[TABLE]

with $t_{\textrm{ramp}}=3$ Myr and $t_{\textrm{end}}=30$ Myr. At $t$ greater than $t_{\textrm{end}}$, the CR source is set to zero. The duration of the burst exceeds both the acoustic and Alfvén travel times, but not by large factors; thus, a steady state is never achieved. The source strength is set to $q_{0}=5.0\times 10^{-25}$ erg cm*-3* s*-1* for our fiducial run. We increase this by factors of 3 and 10 in comparison runs. The region of CR injection extends from the left boundary at $x=0$ to $x=D=1$ pc corresponding to a CR injection rate per unit area of $q_{0}D=1.54\times 10^{-6}$ erg cm*-2s-1* for the fiducial run. In a 1D situation, we expect the cosmic ray pressure between the cloud and the source to build up to a value of order $q_{0}D/(3v_{A})\sim 8.2\times 10^{-14}$ dyne cm*-2*. This is borne out by Fig. 2, and is only about 25% of the thermal pressure, suggesting that the there will be little lateral magnetic field line distortion or spreading, and that the dynamics will remain 1D.

As in W17, the onset of CR injection produces an initial sound wave which propagates outward and tends to reflect off the right boundary of the domain. To prevent the non-physical reflected wave from interfering with our results we place the right boundary at a very large distance (20 kpc). To save computation time, the extent of the domain in the (vertical) $y$-direction is limited to $\pm$125 pc, and we make use of the mesh refinement in FLASH to reduce resolution far to the right of the cloud. The mesh is fixed in time, and is forced to its maximum resolution of 1.95 pc in the region $0<x<3$ kpc.

The simulation includes radiative cooling proportional to $n_{H}^{2}$, with a supplemental heating rate proportional to $n_{H}$:

[TABLE]

The supplemental heating is added to ensure two thermal equilibrium phases at about 104 and 106 K. Its form is consistent with the heating expected from photoelectric expulsion of electrons from dust grains. We use a heating constant of $\Gamma=10^{-25}$ erg s*-1*, which is similar to values used in other works (Inoue et al. (2006) for example). The cooling function is approximated with the analytic function

[TABLE]

where the exponent is

[TABLE]

This form is taken from Bustard et al. (2016), modified by Chad Bustard to fit CLOUDY data from Wiersma et al. (2009) (Chad Bustard, private communication). The comparison of this function to the data is shown in figure 1. We overplot the quantity $\Gamma/n_{H}$ at the initial thermal pressure given above to show the presence of the two thermal equilibria. We turn off radiative cooling for temperatures below 104 K. With this setup we have, at the chosen thermal pressure, one thermal equilibrium at exactly 104 K and one just above 106 K. The equilibrium at 104 K is stable. The equilibrium near 106 K is formally unstable, however the time scales of the growth of deviations from the equilibrium are long compared to the simulation times, so the hot phase of our simulated medium does not change. Although thermal conduction is not included explicitly in our models, numerical diffusion causes the cloud boundary to evolve. This is a small effect on the timescales considered here.

4 Results

4.1 Fiducial Run

The fiducial simulation is run out to 120 Myr, four times the duration of the injected CR pulse. A time series of the density and CR pressure in the vicinity of the cloud is shown in figure 2. The qualitative picture is the same as the 1D simulations in W17: the onset of CRs injected at the left run into the cloud, building up a bottleneck. The resulting CR pressure gradient pushes the cloud to the right, reducing the thermal pressure inside until equilibrium is regained.

The influence of the CRs in this setup turns out to be small. The cloud is only accelerated to speeds of 5-10 km s*-1*, the field lines are only slightly perturbed from their initial configuration, and the CR pressure inside the cloud peaks at about 30% of the total pressure. This is shown in figure 3 which plots various quantities along a horizontal slice through the center of the cloud ($y=0$).

The top panel of figure 3 shows horizontal profiles of the mass density $\rho$ along the line $y=0$, which goes through the center of the cloud. The initial reaction of the cloud to the push from the CRs is to stretch horizontally, reducing the mass density and thermal pressure inside the cloud. This is the same behavior as in the simple 1D simulations in W17. At later times however the density slowly increases. This may be due to radiative cooling causing condensation onto the cloud.

One of the significant results of this work is the finding that the inclusion of radiative cooling prevents the cloud from being destroyed by CR heating as we saw in the 1D simulations without radiative cooling in W17, wherein we found that without radiative cooling, CR heating will destroy a cloud in around 100 Myr. This lack of any erosion of the cloud in these simulations suggests it is possible to accelerate clouds to high velocities with this mechanism without destroying them.

The middle panel of figure 3 shows the horizontal velocity profiles of the gas for different times in the simulation. Other than the initial profile which is uniformly zero, the velocity profiles all show two sharp dropoffs corresponding to the two edges of the cloud. The cloud velocity can be read off as the value of the “plateau” between these two edges - it has risen to 4 km s*-1* by 40 Myr, and accelerates further to about 7 km s*-1* by the last time step at 120 Myr. These velocities fall far short of the observed 100 km s*-1* velocities of HVCs (and also the velocities attained in the 1D models). As the bottom panel shows, this fiducial run has a very weak CR source, resulting in the thermal pressure dominating over the CR pressure everywhere. Since this results in very little field line bending, the dynamics are basically one-dimensional. This is evidenced by the qualitative agreement of these profiles with the 1D results from W17. As to a quantitative comparison, we note that the $\sim 100$ km s*-1* cloud velocities reported in W17 are reached using a CR source 20 times larger than in this 2D simulation. There are a number of other differences between the 1D and 2D simulations, such as magnetic field and initial cloud position (these factors affect at least the timing of the acceleration), but the difference in source strength is likely the largest factor in the difference in cloud velocities. We discuss this comparison further in §4.2.

The bottom panel of figure 3 shows the thermal and CR pressure profiles at a single time slice at 40 Myr, chosen when the CR bottleneck is around its steepest. As in W17 we can clearly see the bottleneck where the CRs exit the cloud at its right edge, near $x=1.15$ kpc. To the right (downstream) of the bottleneck, CR pressure gradients drive the streaming instability and $P_{c}$ drops as $v_{A}$ rises. To the left (upstream) of the bottleneck, the net CR flow is sub-Alfvénic, the streaming instability is shut off, and waves are not generated.

Figure 3 also reveals the limitations of this set of simulations. The 1.95 pc cell size is sufficient to resolve the cloud as a whole, but falls far short of resolving the cloud interfaces. In the interfaces the relevant length scale is set by the balance of radiative cooling with heating mechanisms. Even at the right edge, where the bottleneck provides some CR heating, these length scales are of order a few pc. At the left edge there is no CR pressure gradient and no CR heating, and so the equilibrium length scale is formally zero. Incorporating thermal conduction would alleviate this, but typical Field lengths in this type of system are fractions of a pc at the longest. As Figure 9 of W17 shows, cosmic ray mediated thermal fronts have structure on 0.1 pc scales under the conditions studied here.

The consequences of an underresolved thermal interface are evident in the sharp dip in thermal pressure at the left edge of the cloud shown in the bottom panel of figure 3. Because the temperature drops from 106 K to 104 K in just a few grid cells, the cooling rate can be very spatially sporadic, resulting in the spiky drops in thermal pressure seen here. Presumably there is a non-physical overabundance of material at temperatures around 105.5 K, at the peak of the cooling curve, resulting in more energy lost to radiative cooling than in a fully resolved simulation. The large scale effects of this are not clear, but we speculate that the strange nature of the swept-up “wings” at the top and bottom of the cloud at late times (visible in the left panel of figure 2) may be due to this.

4.2 Stronger CR sources

We ran similar simulations with identical initial conditions, but with the rate of CR energy injection increased by factors of 3 and 10, in an attempt to accelerate the simulated cloud to higher velocities. Due to computational limitations these were not run for as long in simulated time, but instead were run to 59 and 40 Myr respectively. Time series of density and CR energy with magnetic field streamlines for each are shown in figures 4 and 5.

Qualitatively, these runs are similar to the fiducial run. In all cases the onset of the CR source sends a pressure wave through the domain, followed by the CRs themselves. Both of these push on the cloud, accelerating it to some top speed. At these higher source strengths the cloud is accelerated to higher speeds (about 11 km s*-1* and 25 km s*-1* at the end of each run respectively).

With higher CR pressures the 2D effects start to arise. The initial sound wave causes significant flow around the cloud. This causes the magnetic field to roll around the top and bottom of the cloud in a hairpin shape. However, due to flux freezing, the set of field lines which penetrate the cloud stays the same (neglecting the small amount of evaporation or condensation which may occur). As such, the magnetic field topology does not change, and the resulting CR pressure gradients should be roughly the same as the 1D case except that the 2D runs may take longer to build up since CRs have to travel along the longer, rolled-up field lines to reach the cloud.

We do in fact see that the CR pressure upstream of the bottleneck increases significantly as the CR source is increased, as shown in the bottom panels of figures 6 and 7. The 3x run reaches rough equipartition, while in the 10x run CR pressure exceeds the thermal pressure. We also start to see significant CR pressure downstream of the bottleneck. These qualitatively match the behavior of the 1D steady-state solutions in W17, as well as the 1D time-dependent simulation. The 1D simulation used a source strength 20 times that of the fiducial 2D run, equal to $3.08\times 10^{-5}$ erg cm*-3s-1*.

Before discussing the cloud speeds reached in our simulations, we take a moment to address a numerical peculiarity that emerged in the 2D simulation with the strongest CR source. Figure 5, right side, shows the CR energy densities at different times for this run. At 20 Myr and beyond, CR pressure gradients have mostly vanished except at the bottleneck itself. This is in contrast to the other two runs, where lateral CR pressure gradients remain throughout the simulation. At first glance this appears to imply that CRs are somehow flowing from the bottlenecked region to the left of the cloud to the areas above and below the cloud. This would itself imply that the magnetic topology of the cloud has changed, that some field lines which started inside the cloud have left it. But this would imply flux freezing was violated, and having checked the field carefully we do not find this to be the case. The vanishing of lateral pressure gradients instead seems to stem from a numerical effect of our outflow boundary conditions on the top and bottom of the domain. If a bottlenecked field line, that is, a field line which enters the cloud, is bent enough that it leaves the domain, this results in extra CRs being added by streaming into the domain boundary nearby. The details of this non-physical mechanism are laid out in appendix A, but basically it should not affect any of the dynamics along bottlenecked field lines, and so the acceleration of the cloud should not be significantly affected.

The higher pressures are accompanied by faster cloud speeds, shown in the middle panels. Comparing speeds at the same simulation time of 40 Myr between all three simulations, we see that the increase is sub-linear. For CR sources in the ratio 1:3:10, we reach cloud speeds at 40 Myr of around 4, 9, and 25 km s*-1* respectively. These ratios are quite close to the ratios of maximum cosmic ray pressures achieved between cloud and source. This is further evidence that the transverse spreading that occurs in 2D deprives the cosmic ray piston of some of its force. In addition, a certain amount of magnetic field reconnection takes place due to numerical diffusion. This reduces the cosmic ray flux into the cloud, and hence the cosmic ray momentum available for cloud acceleration. While it is unlikely that numerical diffusion faithfully replicates reconnection physics, it does highlight the sensitivity of the cloud - cosmic ray interaction to magnetic field topology.

Although the 1D runs of W17 are quite different from these 2D runs in many ways, we can put them all on a plot of CR source strength versus cloud velocity. For ease of comparison, we choose a time of 30 Myr to measure the cloud velocity - this is the latest time shown in the 1D run, and the CR time-dependences of the CR sources differ strongly between the 1D and 2D runs after this time (the 1D runs have a steady source throughout rather than the pulsed source described in section 3). For each of the 2D runs we take a horizontal line through the cloud center $y=0$ at 30 Myr, and take the minimum and maximum velocities in the cloud, defined as anywhere with temperature less than $2\times 10^{4}$ K (the resulting velocity range is not sensitive to this cutoff). These velocity ranges are plotted against injected CR flux in blue in figure 8. The 1D run in W17 used a source strength twice that of the strongest source used in the 2D runs, equal to $3.08\times 10^{-5}$ erg cm*-3s-1*, and resulted in a cloud speed of about 90 km s*-1* after 30 Myr. This is shown in cyan, shifted to the right slightly to avoid overlap with other points. We emphasize that spatial resolution, initial cloud position, cloud edge width, magnetic field strength and other properties for this particular run differ from the 2D runs.

For a more apples-to-apples comparison we have also performed new 1D simulations in ZEUS which more closely reproduce the spatial resolution and initial conditions of the 2D FLASH runs. These runs use the same initial density and temperature profiles, the same 1 $\mu$G magnetic field, the same time-dependent CR source, and a spatial resolution of 2 pc. However, as with the 1D run from W17, these runs do not include radiative cooling and instead throw away all energy from CR wave heating. The velocity ranges for these are shown as the red lines in figure 8. We have also run one 1D simulation with a higher magnetic field of 3 $\mu$G, shown in orange. A higher magnetic field results in a lower CR pressure for the same injected CR flux (since the Alfvén speed at the bottleneck is higher). However, we see that increasing the field by a factor of 3 decreases the cloud velocity by only a little. This suggests that the injected CR flux is more important than CR pressure gradients in determining cloud velocities.

We see a clear upward trend of cloud velocity with CR source strength for both the 1D and 2D runs. The velocity range centers for the 2D runs are about 2.5, 7, and 20 km s*-1* respectively, which are very close to the 1:3:10 ratio of the CR source strengths, consistent with a linear relationship. This result is promising, and suggests that higher velocities could be reached with proportionally higher CR sources. This conclusion is partially born out by the results of the 1D simulations, which extend to slightly higher CR fluxes. In the 1D runs, we see that the cloud velocity trends linearly with CR source at low values, but falls off somewhat at higher values. The trend among the 2D runs is linear throughout (though there are hints of sublinearity at high CR sources), and all reach lower velocities than their 1D counterparts. This may partially be because radiative cooling increases the mass of the cloud with time relative to the 1D runs (see figure 11 and discussion in §4.3). Since these mass differences are small, the 2D geometry itself seems to be the likelier cause. Although flux freezing implies that ‘strength’ of the CR bottleneck is unaffected by lateral motions (see discussion in §4.2), it is possible that lateral motions reduce the effective ram pressure seen by the cloud from the accelerated ambient medium. In 1D, all rightward gas flow impacts the cloud, but in 2D some is redirected to its sides.

The response of the cloud density to varying CR source strength is far less clear. The three times higher CR source run has a central density profile which drops significantly as the cloud stretches out horizontally, as shown in the top panel of figure 6. The magnitude of this effect lies somewhere between that of the fiducial 2D run in the top panel of figure 3 and the 1D simulation with high CR source in W17, in which the cloud density drops by as much as a factor of 10. However, in 2D increasing the CR source does not further decrease the cloud density. The density profiles for the ten times higher CR source are shown in the top panel of 7, and we can see the central cloud densities at 40 Myr have not significantly decreased from their initial values. 2D effects are likely involved, as the shape of the cloud is quite different in figures 4 and 5. The proper inclusion of radiative cooling, as opposed to the strategy of throwing away all CR wave heating used in the 1D simulations, is likely a major factor as well.

There are further numerical limitations in these simulations in addition to those discussed in §4.1. These runs do not extend as long as the fiducial run in simulation time. This was due to limited computing resources, but is compounded by the fact that at higher CR source strengths, the cloud moves further in the same amount of time, requiring a greater amount of the domain to be at maximum resolution. The extent of the domain in the vertical direction is also limited. In the fiducial run this was of little consequence since the evolution was almost one-dimensional, but at higher CR source strengths lateral motions are significant and we would ideally want to put the lateral domain boundaries farther away.

4.3 The Role of Radiative Cooling

In W17 we found that without any radiative cooling to offest the CR wave heating at the bottleneck, clouds can be eroded on the scales of 10s of Myr. We speculated that radiative cooling would efficiently get rid of the heat deposited by the CRs in this way and simply removed this energy in subsequent 1D runs. One of the goals of this work is to determine whether a proper treatment of radiative cooling supports this assumption. The fact that the cloud’s edges stay sharp222As mentioned in §4.1, the expected equilibrium length scales are at the pc or sub-pc level and are not completely resolved in these simulations. as seen in figures 2-5 implies that radiative cooling is sufficient to prevent cloud destruction by CR heating at bottlenecks.

To isolate the role that radiative cooling plays we also ran the fiducial simulation with radiative cooling turned off. The results are shown in figures 9 and 10. Comparing with the equivalent results from the fiducial run in figures 2 and 3, we see the response of the cloud to the onset of CRs is similar. However, without radiative cooling the cloud does not maintain a sharp boundary, and slowly erodes due to CR heating. As a result the total cloud mass drops slightly with time, in contrast to the fiducial cloud which gains mass due to condensation (see figure 11). This reinforces the idea that radiative cooling can keep clouds intact even when accelerated by CR pressure gradients, and is one of the significant conclusions of this work.

The cloud without cooling reaches velocities of 2.1 - 4.6 km s*-1* in 30 Myr, roughly the same values as the fiducial cloud but slightly higher. This may simply be due to the slightly lower mass of the cloud as it evolves. Peculiarly, unlike the fiducial cloud, the velocity profiles show little spatial variation. While some of this may be attributed to the lack of sharp density gradients, there is little velocity difference between the inside and outside of the cloud, in stark contrast with the fiducial run (see figure 3 middle panel).

4.4 Equivalent Star Formation Rates

We have discussed the trend of cloud velocity with increasing CR flux, but what CR fluxes might we expect in a Milky Way-like galaxy, or in a starburst galaxy and how do they compare to the values used in this work? To find out, let us roughly convert the Milky Way’s star formation rate (SFR) to a CR flux near the disk. The basic picture is as follows. As stars form in the Galactic disk, some fraction of them explode as core-collapse supernovae. Each supernova injects a roughly fixed amount of energy in CRs into the disk. These CRs then flow away from the disk in something approaching plane-parallel symmetry.

Denote the Galactic SFR by $SFR_{*}$, known to be a few $M_{\odot}$ yr*-1*, and denote the core-collapse supernova rate per unit stars formed as $R_{\textrm{SN}}$. For typical initial mass functions (IMF) of star formation, we expect approximately 1 core-collapse supernova for every 100-200 $M_{\odot}$ of stars formed. Each supernova releases about 1051 erg of kinetic energy, some fraction $\xi$ of which goes into CRs. Observations of the CR content in the Milky Way suggest $\xi$ is around 0.1. The total injection rate of energy in CRs throughout the Galaxy $I_{cr}$ is then

[TABLE]

Let us make the gross simplification that this energy is injected uniformly throughout the Galactic disk, which is itself approximated as a thin disk with radius $r_{*}$ of about 10 kpc. The resulting CR energy flux $\Phi_{cr}$ leaving the disk on both sides is then

[TABLE]

This CR flux, which represents the relatively low SFR of the Milky Way, lies in the lower-middle of the range of values tested in these simulations (see figure 8). This suggests that significantly higher SFRs such as those in starburst galaxies could drive clouds of the same mass to much faster speeds. M82 for instance has a higher SFR of about 10 $M_{\odot}$ yr*-1* (see Gao & Solomon (2004)), contained in a much smaller region of radius 200 pc (Förster Schreiber et al. (2003)). Converting this to a CR flux according to (18) and extrapolating figure 8 to this value, we may expect CRs expelled from the starburst region of M82 to accelerate clouds to speeds of  104-105 km s*-1*. This is an optimistic value as it is under the assumption that the CRs drive winds perpendicular to the disk and don’t spread out, and also ignores the possibility of clouds being destroyed by flow instabilities. Furthermore, at the higher densities in starburst galaxies, hadronic losses of CRs, which are ignored in this work, would become non-negligible. Nevertheless, equation (18) and figure 8 suggest that even galaxies with modest SFRs could inject enough CRs to drive a 104 K cloud of this size to speeds of  100 km s*-1*.

4.5 Other Discussion

Figure 8 neatly summarizes the quantitative result of this work - CR bottlenecks can accelerate clouds to speeds which are roughly proportional to the injected CR flux. But there are clearly other factors which determine the cloud’s speed at any given time. We briefly discuss some of these below.

The magnetic field in the vicinity of the cloud will set much of what occurs in its evolution in these simulations. With all other quantities fixed, a higher magnetic field results in a higher streaming speed at the bottleneck. This would lower the steady-state CR pressure for a given CR injection rate, and may result in a slower cloud. A full study of the effects of magnetic field strength would require another set of simulations and is beyond the scope of this paper. However we do have some points of reference for the 1D simulations: the main simulations for comparison to the 2D simulations, which have a 1 $\mu$G field; an equivalent 1D simulation with a 3 $\mu$G field; and the old 1D simulation from W17 which has a 3.26 $\mu$G field (as well as a number of other differences).

We see from figure 8 that the two runs with higher magnetic field do indeed lie somewhat lower on the plot, i.e. their clouds reach slower speeds than the fiducial 1 $\mu$G field runs. However, the reduction in speed is small, much less than the factor of 3 difference in bottleneck streaming speeds. This suggests that cloud velocity is set more by the incoming CR flux than by the steady-state CR pressure at the bottleneck. On the other hand, with different streaming speeds the timing of events is different, and it may be incorrect to compare runs with different magnetic fields at the same time slices.

Another significant factor for which we have only 1 data point is the total cloud mass. All simulations here have the same column density through the center of the cloud of 7.2$\times 10^{-5}$ g cm*-2* (the old 1D run from W17 has a slightly smaller column density of 6.1$\times 10^{-5}$ g cm*-2*). We therefore cannot make any concrete statements about the dependence of cloud velocity on cloud mass. In multiple dimensions the shape of the cloud is also important, as its cross section determines how wide the CR piston is while its volume is tied to its mass. A flat cloud which is parallel to the disk might reach higher speeds than a spherical cloud of the same mass. We leave the study of these questions to future work.

4.6 Future Work

There is ample room for extension of these simulations. The greatest computational burden arises from the time constraint required from modeling CR streaming. Using the tanh regularization method of Sharma et al. (2010), the simulation time step must be very small,

[TABLE]

where $\Delta x$ is the grid cell size, $v_{s}$ is the CR streaming speed, and $L$ is a smoothing length used in the regularization. The longer $L$ is, the more accurate the simulation, but the harsher the time constraint. For reasonable results we typically require $L$ to be about the size of the domain. For all simulations shown here we use $L=10$ kpc.

However, Jiang & Oh (2018) have recently introduced a new numerical method to handle CR streaming which is based on a radiative transfer method. This method has no smoothing length $L$, and so the associated time constraint is a simple CFL condition

[TABLE]

which is longer than (19) by orders of magnitude. Implementing this method in the future will allow us to address many of the numerical issues present here, such as limited spatial resolution, vertical domain extent, and total simulated time.

5 Conclusion

We presented a series of two-dimensional MHD simulations exploring the dynamic effects of a CR bottleneck on a cool 104 K cloud embedded in the hot 106 K CGM. In earlier 1D simulations (W17), the CR pressure gradient resulting from the bottleneck was found to accelerate clouds to high velocities ($\sim 100$ km s*-1*) and significantly reduce the density inside the cloud. We explored how these results changed with the inclusion of the possibility of lateral motion around the cloud.

We found that at low CR source strengths, lateral pressure gradients are not enough to bend magnetic field lines significantly, and the evolution is largely one-dimensional, and our results resemble those of W17. We further found that the inclusion of radiative cooling prevents the cloud from being destroyed by CR heating.

At higher CR strengths approaching those required to reach cloud speeds of 100 km s*-1* in the 1D runs, lateral effects become very important. High pressures induce flow around the cloud, wrapping field lines up around its sides. Despite this, cloud velocity increases for higher CR sources in a nearly-linear correlation up to the source strengths tested here, suggesting it is possible for CR bottlenecks to accelerate clouds to the speeds of observed HVCs. A Milky Way sized galaxy with a star formation rate of 20 - 30 $M_{\odot}$ yr*-1* would generate enough CRs to accelerate clouds to 100 km s*-1*. However, limited computational resources make further testing prohibitive until the development of CR bottleneck simulations which take advantage of new numerical methods for streaming.

Acknowledgements

We acknowledge Marcus Brüggen for useful discussions. We acknowledge Karen Yang for her support with using the FLASH code. We acknowledge Chad Bustard for additional code support as well as supplying the analytic fit to the cooling curve. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. JW and EGZ acknowledge support by NSF Grant AST-1616037, the WARF Foundation, and the Vilas Trust. MR acknowledges support by NSF Grant AST-1715140.

Appendix A Boundary-Related Issues

In §4.2 we alluded to a numerical issue with our boundary conditions which resulted in the erasure of lateral CR pressure gradients in our last 2D run shown in figure 5. We describe here the details of this issue, which occurs wherever there is a lateral CR pressure gradient across adjacent field lines which intersect the top or bottom domain boundary at an angle other than 90 degrees.

First we present evidence that field lines leaving the domain are the cause of the loss of lateral CR pressure gradients seen in 5. Figure 12 shows a gas density plot at 20 Myr for each of the three 2D runs presented in this work. Overlaid on top of the density are streamlines which follow the magnetic field and are colored according to CR pressure. Yellow implies high CR pressure, while purple indicates low CR pressure.

In the first two runs, the lateral CR pressure gradients clearly remain, as seen by the color contrast between field lines. As seen from the plots, the border between the high CR pressure in yellow and the lower CR pressure in green is formed by the top- and bottom-most magnetic field lines which still intersect the cloud. This is consistent with the idea that CRs streaming along field lines which intersect the cloud are effectively bottlenecked, while those streaming along other field lines can stream faster. Note that none of the bottlenecked field lines leave the domain out the top or bottom.

In the third run, however, there is no lateral color contrast seen, indicating there are no CR pressure gradients perpendicular to the field lines. Note that the field lines where we would expect such a gradient to be present, which again are the top- and bottom-most lines which intersect the cloud, also leave the domain out the top and bottom.

This removal of cross-field CR gradients appeared in one of our very early tests of the FLASH code, where we injected CRs along the left boundary and had them stream down a titled B-field toward a dense cloud. In these tests we turned off the MHD and only evolved the CRs in time. The top and bottom boundaries were periodic, and the left boundary was reflecting. We expected the CRs to build up in pressure on field lines which entered the cloud. Instead we saw roughly the same CR pressure everywhere, except for the reduced pressure past the bottleneck (see figure 13, left). This did not fit with our theory.

We found that our results were consistent with our theory only if we negated any CR flux through the left boundary (see figure 13, right). This suggested there was some unintended CR flux through the boundary. A schematic for how the boundary conditions cause this is shown in figure 14. This figure illustrates the issue when reflecting boundary conditions are used, but the same problem arises when using outflow boundary conditions. (a) Consider two adjacent field lines, one which enters the cloud (shown in blue) and one which doesn’t (shown in red). (b) CR pressure will build up on one line because of the bottleneck effect, but not on the other. Once this buildup reaches the boundary, the boundary conditions fill in the guard zones (c) according to the condition chosen. Outflow conditions will copy values from the first active cell to all the guard cells. Reflecting conditions will mirror the values of the first active cells, with normal the gas velocity component flipping sign (normal magnetic field is not flipped). Guard cells across the boundary from a bottlenecked field line will be given large CR pressures. But this induces a CR gradient along the adjacent field line in red, causing CRs to flow into the domain along this line via streaming (d). This raises the CR pressure on a field line even though it doesn’t pass through the cloud. This unphysically high CR pressure then propagates to further field lines (e,f).

The same effect can actually cause CRs to leave the domain unphysically as well. This is shown in figure 15, which considers two adjacent field lines on the opposite end of the cloud. Here the boundary conditions (either outflow or reflecting) introduce a boundary cell with low CR pressure along one of the bottlenecked field lines. The resulting pressure gradient causes CRs to stream out of the domain, reducing the pressure in the bottleneck. This explains why the CR pressure in the bottlenecked region of figure 13 left is less than that of figure 13 right despite using the same setup.

We applied the fix of negating all CR streaming flux through the left boundary to avoid this problem. We also ended up negating advective flux of CRs through this boundary for a different reason - to control the amount of CRs being added to the domain. However, as we wanted CRs to be able to stream out of the top and bottom of the domain, the streaming flux was not negated there. Unexpectedly, in the third 2D simulation field lines connected to the cloud left the domain through the top and bottom, so this issue reappeared in a different form. This simulation would properly be redone with a domain which is wider in the vertical direction to accommodate the spreading of the field lines.

According to the geometry of the magnetic field lines leaving the domain in our cloud simulations, e.g. figure 12 bottom panel, we only expect CRs to be spuriously added (as in figure 14) to the domain along non-bottlenecked field lines. We don’t expect any CRs to be spuriously removed (as in figure 15) anywhere. This implies that the CR dynamics along bottlenecked field lines are not significantly impacted by these boundary issues. Hence, we do not expect the acceleration of the cloud to be affected much.

Appendix B Resolution Test

We unfortunately did not have enough computing resources to do a proper resolution test. Instead we present a qualitative test wherein we ran the fiducial run at half the spatial resolution. Density slices for this run are shown in figure 16 alongside the equivalent time slices of the fiducial run. Similarly, CR energy densities with magnetic field lines overlaid are shown in 17. While there are some clear differences, the overall behavior is largely the same between the two resolutions.

We can make a more direct comparison by plotting quantities along the horizontal line through the center of the cloud as in figure 3. This is shown in figure 18. The top two panels show the density and horizontal velocity profiles respectively. The solid lines show the fiducial run, and the dashed lines show the lower resolution run. There is very good agreement between the density profiles. Agreement is not as good between the velocity profiles, but they are close to each other in the cloud region.

The bottom panel of figure 18 shows the thermal pressure profiles (in solid lines) and CR pressure profiles (in dashed lines) at 40 Myr, with the fiducial run in orange and the low resolution run in red. These also show very good agreement.