A numerical approach to the non-uniqueness problem of cosmic ray two-fluid equations at shocks

TL;DR

This paper examines the non-uniqueness issues in modeling cosmic ray two-fluid equations at shocks and proposes a robust numerical method that reduces dependence on discretization details, emphasizing the importance of subgrid closure for accurate shock physics.

Contribution

It introduces a numerical approach that mitigates non-uniqueness in cosmic ray two-fluid models at shocks and highlights the significance of subgrid closure based on kinetic theory.

Findings

Results depend on numerical methods without subgrid closure.

A robust method reduces sensitivity to numerical details.

Subgrid closure improves shock microphysics modeling.

Abstract

Cosmic rays (CRs) are frequently modeled as an additional fluid in hydrodynamic (HD) and magnetohydrodynamic (MHD) simulations of astrophysical flows. The standard CR two-fluid model is described in terms of three conservation laws (expressing conservation of mass, momentum and total energy) and one additional equation (for the CR pressure) that cannot be cast in a satisfactory conservative form. The presence of non-conservative terms with spatial derivatives in the model equations prevents a unique weak solution behind a shock. We investigate a number of methods for the numerical solution of the two-fluid equations and find that, in the presence of shock waves, the results generally depend on the numerical details (spatial reconstruction, time stepping, the CFL number, and the adopted discretization). All methods converge to a unique result if the energy partition between the thermal and non-thermal fluids at the shock is prescribed using a subgrid prescription. This highlights the non-uniqueness problem of the two-fluid equations at shocks. From our numerical investigations, we report a robust method for which the solutions are insensitive to the numerical details even in absence of a subgrid prescription, although we recommend a subgrid closure at shocks using results from kinetic theory. The subgrid closure is crucial for a reliable post-shock solution and also its impact on large scale flows because the shock microphysics that determines CR acceleration is not accurately captured in a fluid approximation. Critical test problems, limitations of fluid modeling, and future directions are discussed.

Tables

Table 1. Table 1: Various methods for implementing the non-conservative term of CR-HD equations † .

Assumptions used		Coupling terms	Nomenclature & References e.g.
in numerical implementation
		A.  p dv	A1. Eg+Ecr (OpSplit-pdv)
			Pfrommer et al. 2006; Sharma et al. 2009;
		$+p_{\mathrm{cr}}\nabla \cdot 𝐯$	Salem & Bryan 2014; Butsky & Quinn 2018
	Option 1: $Eg+Ecr$	$-p_{\mathrm{cr}}\nabla \cdot 𝐯$	A2. Eg+Ecr (Unsplit-pdv)
	Gas energy + CR energy eqs		This work
		B.  v dp	B1. Eg+Ecr (OpSplit-vdp)
		$-𝐯\cdot \nabla p_{\mathrm{cr}}$
		$+𝐯\cdot \nabla p_{\mathrm{cr}}$	B2. Eg+Ecr (Unsplit-vdp)
Two-fluid CR-HD		A.  p dv	A1. Et+Ecr (OpSplit-pdv)
		$0$	Yang et al. 2012;
	Option 2: $Et+Ecr$	$-p_{\mathrm{cr}}\nabla \cdot 𝐯$	A2. Et+Ecr (Unsplit-pdv)
	Total energy + CR energy eqs		This work
		B.  v dp	B1. Et+Ecr (OpSplit-vdp)
		$0$
		$+𝐯\cdot \nabla p_{\mathrm{cr}}$	B2. Et+Ecr (Unsplit-vdp)
	Option  3: $Et+Scr$		Et+Scr
	Total energy + CR entropy eqs		Kudoh & Hanawa (2016)

Table 2. Table 2: Test Problems

Test	Initial Conditions	$t_{\mathrm{stop}}$	$N_x$	Figure
	$\{\rho ,v_x,p_g,p_{\mathrm{cr}}\}$
$1$. Pressure balance	$\begin{array}{cc}L:\{1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.9}\}& \\ R:\{1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}0.9},\mathrm{\hspace{0.17em}0.1}\}& \end{array}$	$1.0$	$200$	3
$2$. Shock tube A	$\begin{array}{cc}L:\{1,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}1}\}& \\ R:\{0.2,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0.02},\mathrm{\hspace{0.17em}0.1}\}& \end{array}$	$0.1$	$1000$	4, 6, 6
$3$. Shock tube B	$\begin{array}{cc}L:\{1,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}6.7}\times {10}^4,\mathrm{\hspace{0.17em}1.3}\times {10}^5\}& \\ R:\{0.2,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}2.4}\times {10}^2,\mathrm{\hspace{0.17em}2.4}\times {10}^2\}& \end{array}$	${10}^{-4}$	$1000$	7
$4$. Blast wave	see section 5.4	section 5.4	section 5.4	8

Equations

$$\begin{array}[]{l}\displaystyle\frac{\partial f}{\partial t}+\left({\bf v}+{\bf v}_{\rm st}\right)\cdot{\nabla}f=\frac{\mathcal{P}}{3}\frac{\partial f}{\partial\mathcal{P}}\nabla\cdot\left({\bf v}+{\bf v}_{\rm st}\right)\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle+{\nabla}\cdot\left[D_{\rm s}\hat{b}\left(\hat{b}\cdot{\nabla}f\right)\right]+\frac{1}{\mathcal{P}^{2}}\frac{\partial}{\partial\mathcal{P}}\left(\mathcal{P}^{2}D_{\rm\mathcal{P}}\frac{\partial f}{\partial\mathcal{P}}\right).\end{array}$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}+\left(\mathbf{v}+{\mathbf{v}}_{\mathrm{st}}\right)\cdot \nabla e_{\mathrm{cr}}$$

$$\frac{\partial \mathcal{U}}{\partial t}+\nabla \cdot \mathsf{F}=\mathcal{S},$$

$$\mathsf{F}={\left[\begin{array}{c}\rho \mathbf{v}\\ \\ \rho \mathbf{v}\mathbf{v}+p_{\mathrm{t}}\mathsf{I}\\ \\ \left(e_{\mathrm{t}}+p_{\mathrm{t}}\right)\mathbf{v}-{\kappa }_{\mathrm{cr}}\nabla e_{\mathrm{cr}}\\ \\ e_{\mathrm{cr}}\mathbf{v}-{\kappa }_{\mathrm{cr}}\nabla e_{\mathrm{cr}}\end{array}\right]}^{⊺},\mathcal{S}=\left[\begin{array}{c}0\\ \\ 0\\ \\ 0\\ \\ -p_{\mathrm{cr}}\nabla \cdot \mathbf{v}\end{array}\right]$$

$$e_{\mathrm{t}}=e_{\mathrm{g}}+e_{\mathrm{cr}}=\left(\frac{1}{2}\rho |\mathbf{v}{|}^2+\frac{p_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}\right)+\frac{p_{\mathrm{cr}}}{{\gamma }_{\mathrm{cr}}-1},$$

$$p_{\mathrm{t}}=p_{\mathrm{g}}+p_{\mathrm{cr}}$$

$$\frac{\partial \mathcal{V}}{\partial t}+{\mathsf{A}}_{\mathrm{p}}(\mathcal{V})\frac{\partial \mathcal{V}}{\partial x}=0,$$

$$\mathcal{V}=\left[\begin{array}{c}\rho \\ v_{\mathrm{x}}\\ v_{\mathrm{y}}\\ v_{\mathrm{z}}\\ p_{\mathrm{g}}\\ p_{\mathrm{cr}}\end{array}\right],{\mathsf{A}}_{\mathrm{p}}=\left[\begin{array}{cccccc}{v}_{\mathrm{x}}& \rho & 0& 0& 0& 0\\ 0& v_{\mathrm{x}}& 0& 0& \frac{1}{\rho }& \frac{1}{\rho }\\ 0& 0& v_{\mathrm{x}}& 0& 0& 0\\ 0& 0& 0& v_{\mathrm{x}}& 0& 0\\ 0& \rho a_{\mathrm{g}}^2& 0& 0& v_{\mathrm{x}}& 0\\ 0& \rho a_{\mathrm{cr}}^2& 0& 0& 0& v_{\mathrm{x}}\end{array}\right].$$

$${\lambda }_{\mathrm{p}}=v_{\mathrm{x}}-a,v_{\mathrm{x}},v_{\mathrm{x}},v_{\mathrm{x}},v_{\mathrm{x}},v_{\mathrm{x}}+a.$$

$${\mathsf{R}}_{\mathrm{p}}=\{r_{\mathrm{p}}^{mn}\}$$

$$\left[\rho v\right]$$

$$\left[\rho v^2+p_{\mathrm{g}}+p_{\mathrm{cr}}\right]$$

$$\left[\left(\frac{1}{2}\rho v^2+\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}{p}_{\mathrm{g}}+\frac{{\gamma }_{\mathrm{cr}}}{{\gamma }_{\mathrm{cr}}-1}{p}_{\mathrm{cr}}\right)v\right]$$

$${\mathcal{M}}_{\mathrm{g},1}=\frac{v_1}{({\gamma }_{\mathrm{g}}{p}_{\mathrm{g},1}/{\rho }_1{)}^{\frac{1}{2}}},{\mathcal{M}}_{\mathrm{cr},1}=\frac{v_1}{({\gamma }_{\mathrm{cr}}{p}_{\mathrm{cr},1}/{\rho }_1{)}^{\frac{1}{2}}},$$

$$\begin{array}{r}{\mathcal{M}}_1^2+\frac{1+{\mathcal{P}}_{\mathrm{cr},1}}{{\gamma }_{\mathrm{g}}+{\gamma }_{\mathrm{cr}}{\mathcal{P}}_{\mathrm{cr},1}}=\frac{{\mathcal{M}}_1^2}{\mathcal{R}}+\frac{{\mathcal{P}}_{\mathrm{g},2}+{\mathcal{P}}_{\mathrm{cr},2}}{{\gamma }_{\mathrm{g}}+{\gamma }_{\mathrm{cr}}{\mathcal{P}}_{\mathrm{cr},1}}\end{array}$$

$$\begin{array}{r}\mathcal{R}\left\{\frac{{\mathcal{M}}_1^2}{2}+\frac{1}{{\gamma }_{\mathrm{g}}+{\gamma }_{\mathrm{cr}}{\mathcal{P}}_{\mathrm{cr},1}}\left(\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}+\frac{{\gamma }_{\mathrm{cr}}}{{\gamma }_{\mathrm{cr}}-1}{\mathcal{P}}_{\mathrm{cr},1}\right)\right\}=\\ \left\{\frac{{\mathcal{M}}_1^2}{2\mathcal{R}}+\frac{1}{{\gamma }_{\mathrm{g}}+{\gamma }_{\mathrm{cr}}{\mathcal{P}}_{\mathrm{cr},1}}\left(\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{g}}-1}{\mathcal{P}}_{\mathrm{g},2}+\frac{{\gamma }_{\mathrm{cr}}}{{\gamma }_{\mathrm{cr}}-1}{\mathcal{P}}_{\mathrm{cr},2}\right)\right\},\end{array}$$

$${\mathcal{P}}_{\mathrm{cr},1}=\frac{p_{\mathrm{cr},1}}{p_{\mathrm{g},1}}=\left(\frac{{\gamma }_{\mathrm{g}}}{{\gamma }_{\mathrm{cr}}}\right){\left(\frac{{\mathcal{M}}_{\mathrm{g},1}}{{\mathcal{M}}_{\mathrm{cr},1}}\right)}^2.$$

$$w_{\mathrm{cr}}=\frac{p_{\mathrm{cr},2}}{p_{\mathrm{g},2}+p_{\mathrm{cr},2}}=\frac{{\mathcal{P}}_{\mathrm{cr},2}}{({\mathcal{P}}_{\mathrm{g},2}+{\mathcal{P}}_{\mathrm{cr},2})},$$

$${ϵ}_{\mathrm{cr}}$$

$$p_{\mathrm{cr},2}=p_{\mathrm{cr},1}{\left(\frac{{\rho }_2}{{\rho }_1}\right)}^{{\gamma }_{\mathrm{cr}}}\to {\mathcal{P}}_{\mathrm{cr},2}={\mathcal{P}}_{\mathrm{cr},1}{\mathcal{R}}^{{\gamma }_{\mathrm{cr}}}.$$

$$\begin{array}[]{l}\displaystyle\quad{\cal U}^{*}={\cal L}_{h}({\cal U}^{n})\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\quad{\cal U}^{n+1}={\cal L}_{s}({\cal U}^{*})\,,\end{array}$$

$${\mathcal{U}}^{n+1}=({\mathcal{L}}_h+{\mathcal{L}}_s)({\mathcal{U}}^n).$$

$$\text{noindent}\frac{\partial e_{\mathrm{g}}}{\partial t}$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}$$

$$e_{\mathrm{cr},i}^{\ast }$$

$$e_{\mathrm{cr},i}^{n+1}$$

$$v_{i+\frac{1}{2}}^{\ast }=\frac{1}{2}(v_i^{\ast }+v_{i+1}^{\ast }).$$

$$e_{\mathrm{cr},i}^{n+1}=e_{\mathrm{cr},i}^n-\Delta t\left[\left(\frac{f_{i+\frac{1}{2}}^n-f_{i-\frac{1}{2}}^n}{\Delta x_i}\right)+p_{\mathrm{cr},i}^n\left(\frac{v_{i+\frac{1}{2}}^n-v_{i-\frac{1}{2}}^n}{\Delta x_i}\right)\right].$$

$$v_{\mathrm{i}+\frac{1}{2}}^n=\frac{{\mathcal{U}}_{[m_x],i+\frac{1}{2}}^n}{{\mathcal{U}}_{[\rho ],i+\frac{1}{2}}^n},$$

$$p_{\mathrm{cr},\mathrm{i}}^n=\frac{{\gamma }_{\mathrm{cr}}-1}{2}\left(e_{\mathrm{cr},i-\frac{1}{2}}^n+e_{\mathrm{cr},i+\frac{1}{2}}^n\right),$$

$$\frac{\partial e_{\rm cr}}{\partial t}=-{\nabla}\cdot\left[(e_{\rm cr}+p_{\rm cr}){\bf v}\right]+\color[rgb]{0,0,0}{\bf v\cdot\nabla}p_{\rm cr}\,,$$

$$\frac{\partial e_{\mathrm{t}}}{\partial t}=-\nabla \cdot \left[(e_{\mathrm{t}}+p_{\mathrm{t}})\mathbf{v}\right],$$

Full text

††pagerange: A numerical approach to the non-uniqueness problem of cosmic ray two-fluid equations at shocks–C

1 Introduction

Macroscopic extension of the cosmic ray (CR) transport equation in the form of the two-fluid CR-HD/MHD equations (Skilling 1971; Drury & Voelk 1981) is very convenient to study the effects of non-thermal processes in astrophysical systems. A two fluid-model provides important insights about the dynamical effects of CRs on large scales (such as blast waves, wind bubbles, galaxies; see e.g., Chevalier 1983; Salem & Bryan 2014; Gupta et al. 2018a), which are expensive to capture from the CR kinetic theory (see Drury 1983 for a review). In many astrophysical systems, the CR energy density is comparable to the thermal/magnetic energy density, and the CR pressure cannot be ignored.

A fluid description of CRs is justified since the Larmor radius of energy density-dominating CRs ($\sim 10^{-5}\,E_{\rm 10GeV}/B_{\rm\mu G}\,{\rm pc}$) is much smaller than the scales of interest and they are confined along the direction of magnetic fields by self-generated magnetic fluctuations at this scale (Kulsrud & Pearce 1969). Therefore, the two-fluid model is applicable in a variety of astrophysical systems, ranging from a star-forming cloud to clusters of galaxies.

Many astrophysical phenomena on large scales are studied with the two-fluid model of CRs. CRs represent an attractive agent for feedback heating as their energy loss time scale is much longer than the cooling time of the thermal gas. CRs retain energy in the cloud for a long time and provide extra pressure leading to a moderate star formation rate. For a similar reason, CRs help to launch galactic winds (e.g., Salem & Bryan 2014; Wiener et al. 2017). Investigations with two-fluid models also found that CRs can reduce the temperature of the circumgalactic medium and account for the observed absorption lines of various elements (e.g., Butsky & Quinn 2018). CR heating can be an efficient heating mechanism in galaxy clusters (e.g., Guo & Oh 2008) and can affect buoyancy instabilities in the intracluster medium (e.g., Sharma et al. 2009). These important conclusions are based on numerical simulations of the two-fluid model and therefore it is necessary to closely examine the properties of these equations and their numerical implementation.

The two-fluid CR-HD model contains three conservation laws (for mass, momentum, and total energy) plus one additional equation, which governs the evolution of CR energy density. Since the CR particle density is usually negligible compared to the thermal particle density, the evolution of the CR particle density is not considered. However, a CR particle carries significant kinetic energy, which makes the energy density contribution of CRs (integrated over entire spectrum) comparable to that of the thermal plasma.

Similar two-fluid equations are also used in other contexts. For example, the two-fluid MHD equations (with pressures along and perpendicular to field lines and equations governing their evolution, and the generalization to include electrons and ions separately) are used to study pressure anisotropy in astrophysical and fusion plasmas (e.g., Sharma et al. 2006; Jardin et al. 2007). Both CR-HD/MHD and MHD with anisotropic pressure have two internal energy equations but they do not have a separate density/velocity equation for the second fluid. In the CR two-fluid model, the equation that determines the evolution of CR energy density cannot be cast in a conservative form owing to the presence of a source (or coupling) term containing spatial derivatives (either $p_{\rm cr}\nabla\cdot{\bf v}$ or ${\bf v}\cdot\nabla p_{\rm cr}$). While this term does not represent an issue for smooth flows, it poses a serious challenge in determining the weak form of the equations. Therefore, the shock jump conditions are non-unique and the solutions across a shock depend on numerical details. Note that there are other two-fluid systems that do not have the spatial derivatives in the source terms and for these the solution across a shock is unique (e.g., Balsara et al. 2016).

The numerical discretization of the non-conservative terms in the two-fluid CR-HD/MHD system can be done in different ways. However, these implementations may not produce an identical solution. This crucial point was highlighted by Kudoh & Hanawa (2016). They suggested that if the CR energy density is estimated from the advection of a passive scalar, namely $\chi=p^{1/\gamma_{\rm cr}}_{\rm cr}/\rho$ (where $\gamma_{\rm cr}$ is adiabatic index for CR fluid, $p_{\rm cr}$ is CR pressure, $\rho$ is gas density), then the results are consistent with the underlying mathematical formalism. This is equivalent to positing that the CR entropy $p_{\rm cr}/\rho^{\gamma_{\rm cr}}$ is conserved across a shock. One advantage of this method is that the two-fluid equations apparently become conservative, which makes numerical application of Godunov-type shock-capturing schemes straightforward. However, the assumption of a constant CR entropy across a shock is inconsistent with the idea that most CRs are accelerated at shocks. This motivates us to explore more physically motivated strategies.

In this paper, we discuss several numerical discretizations of the two-fluid CR-HD equations by implementing the equations in the PLUTO code (Mignone et al. 2007). We find that for most of the common methods, the numerical solutions depend on the choice of spatial reconstruction, time stepping, and even the CFL number. We suggest a method which gives robust numerical results (i.e., solutions are less sensitive to these choices compared to other methods). A physically reliable two-fluid implementation of CRs must be calibrated with the results from kinetic (e.g., Particle-In-Cell) simulations. Our results are crucial for simulating the impact of CRs on large scales (e.g., galaxy), where the CR acceleration microphysics is difficult to resolve.

The non-uniqueness problem that we discuss in this paper applies to both the two-fluid hydro and MHD equations but we exclusively consider the former to focus on this problem and its possible solution in a simpler setting. Generalization to MHD is straightforward.

We organize this paper as follows. After presenting the basic equations in section 2, we discuss the closure problems in section 3. Different methods for implementing two-fluid equations are given in section 4. Section 5 presents results of various test problems. Section 6 discusses broader implications of our work. Our main results are summarized in section 7.

2 Governing Equations

The two-fluid CR HD/MHD equations are obtained from the Fokker-Planck CR transport equation (Skilling 1975), which is given by

[TABLE]

Here $f({\bf x},{\bf\mathcal{P}},t)$ is the CR distribution function assumed to be isotropic in momentum space, $\mathcal{P}$ is the CR momentum, $D_{\rm s}$ and $D_{\rm\mathcal{P}}$ are the diffusion coefficients in spatial and momentum space, ${\bf v}$ is the velocity of thermal plasma. The term ${\bf v}_{\rm st}=(v_{\rm x}^{st},v_{\rm y}^{st},v_{\rm z}^{st})^{\intercal}$ represents the bulk velocity of the scattering centers of CR particles w.r.t. the background plasma, known as the streaming velocity, which is along the direction of the magnetic field ($\hat{b}$) but down the gradient of CR pressure (for a brief discussion see Appendix A in Pfrommer et al. 2017). The first term on the RHS represents the CR convection term, the second term is spatial diffusion while the third term represents CR diffusion in the momentum space (Skilling 1975).

Drury & Voelk (1981) suggested that Eq. (1) can be simplified if one rewrites it in terms of macroscopic variables. This is done by taking the energy moment of Eq. (1), which yields

[TABLE]

Here $\kappa_{\rm cr}$ is the CR diffusion coefficient integrated over the CR distribution function (see equation $7$ in Drury & Voelk 1981). CR streaming and anisotropic diffusion along field lines cannot be captured in hydrodynamics, therefore, we take $|{\bf v}|/|{\bf v}_{\rm st}|\gg 1$ and $\kappa_{\rm cr}\hat{b}(\hat{b}\cdot{\nabla}e_{\rm cr})\approx\kappa_{\rm cr}{\nabla}e_{\rm cr}$. It is also assumed that CR diffusion in momentum space is negligible, so that $\Gamma_{\rm m}\rightarrow 0$. This leads to the two-fluid CR-HD equations, written as

[TABLE]

where ${\cal U}=\{\rho,\,\rho\,{\bf v},\,e_{\rm t},\,e_{\rm cr}\}^{\intercal}$ is the conservative variable array, while the flux tensor and the source term are given, respectively, by

[TABLE]

where

[TABLE]

and

[TABLE]

are, respectively, the total energy density (the sum of kinetic energy density, gas thermal energy density, and CR energy density $e_{\rm cr}$) and the total pressure (the sum of gas pressure $p_{\rm g}$ and CR pressure $p_{\rm cr}$). The adiabatic constants of the thermal and CR fluids are defined as $\gamma_{\rm g}=1+p_{\rm g}/e_{\rm th}$ (where thermal energy density $e_{\rm th}=e_{\rm g}-\rho\,|{\bf v}|^{2}/2$) and $\gamma_{\rm cr}=1+p_{\rm cr}/e_{\rm cr}$, and their values are taken to be $5/3$ and $4/3$ respectively. Note that the last element of ${\cal U}$ ($e_{\rm cr}$) does not obey a conservative equation. CRs can lose/gain energy due to the term $-p_{\rm cr}{\nabla\cdot{\bf v}}$ representing coupling between thermal and CR fluids. This term, involving a derivative, gives a non-zero and non-unique contribution across a shock which depends on numerical implementation. Later we explore various implementations of this term and the associated numerical challenges.

2.1 Characteristic structure

In order to understand the characteristic structure of the two-fluid equations, we neglect the CR diffusion term (i.e., the term $\kappa_{\rm cr}{\nabla}e_{\rm cr}$ of Eq. 4)111Diffusion term is usually separately implemented using the standard parabolic schemes (e.g., Vaidya et al. 2017). and focus only on the hyperbolic structure of the equations. Considering 1D Cartesian coordinate with three velocity components, the primitive form of Eq. (3) then becomes

[TABLE]

where

[TABLE]

Here $a_{\rm g}=(\gamma_{\rm g}p_{\rm g}/\rho)^{\frac{1}{2}}$ and $a_{\rm cr}=(\gamma_{\rm cr}p_{\rm cr}/\rho)^{\frac{1}{2}}$. Defining the effective sound speed as $a=(a^{2}_{\rm g}+a^{2}_{\rm cr})^{\frac{1}{2}}$, the eigenvalues of the characteristic matrix $A_{\rm p}$ are found to be (also see Fig. 1)

[TABLE]

Right eigenvectors of the characteristic matrix $\mathsf{A}_{\rm p}$ can be written as,

[TABLE]

In the absence of CR fluid, the $5^{\rm th}$ column and the $6^{\rm th}$ row of the matrix $\mathsf{R}_{\rm p}$ are absent, i.e., the right-going eigenvalue is $\lambda=v_{\rm x}+a_{\rm g}$ and the corresponding eigenvector $r_{\rm 5}=\{1,a_{\rm g}/\rho,0,0,a_{\rm g}^{2}\}^{\intercal}$. This implies that the ratios between the fluctuations in density, velocity, and pressure scale as $1:a_{\rm g}/\rho:0:0:a_{\rm g}^{2}$. Therefore, any change in the x-velocity in a right going sound wave causes a simultaneous change in the density and the pressure, which increase/decrease without any ambiguity. However, the presence of CRs can introduce ambiguity in the fluctuations, which is illustrated as follows.

As shown in Eq. (10), for the right-going sound waves in the two-fluid CR equations the ratios between the fluctuations in density, velocity, thermal pressure, and CR pressure scale as $1:a/\rho:0:0:a_{\rm g}^{2}:a_{\rm cr}^{2}$, i.e., the fluctuation of the x-velocity can be satisfied for multiple values of $a_{\rm g}$ and $a_{\rm cr}$ [because $a=(a^{2}_{\rm g}+a^{2}_{\rm cr})^{\frac{1}{2}}$]. This ambiguity affects the fluctuations of thermal and CR pressures. Moreover, the $5^{\rm th}$ column of the right eigenvector shows that at the contact discontinuity (represented by $r_{2}$ to $r_{5}$), the ratio between the fluctuations of the thermal pressure and CR pressure is $-1$. This means that, although the total pressure ($p_{\rm t}$, see Eq. 6) on both side of the contact discontinuity is constant, the thermal and CR pressures may not remain the same across a contact discontinuity. The fluctuation in thermal pressure across the contact discontinuity is compensated by a gain/loss of the CR pressure or vice versa. This subtlety in the characteristic structure of the two fluid CR equations leads to numerical issues addressed in the next sections.

3 Non-uniqueness problem at shocks

Related to the ambiguities discussed in the previous section, there are serious problems regarding the uniqueness of the solutions of the two-fluid equations at shocks. In two-fluid CR-HD, CR diffusion coefficient ($\kappa_{\rm cr}$) and CR adiabatic index ($\gamma_{\rm cr}$) need to be specified. These are generally known as the closure parameters, as they must be prescribed to close the system of equations (Eq 3). Although these choices are crucial to quantify the overall impact of CRs, their role is significantly important at shocks. For example, if one considers the time-evolution of CR distribution function (see e.g., Duffy et al. 1994), then a self-consistent treatment for $\gamma_{\rm cr}$ and $\kappa_{\rm cr}$ in the two-fluid model is necessary. However, the fluid equations are not sufficient to estimate these parameters (e.g., Achterberg et al. 1984). Here, we show that an additional equation of state is required to obtain a unique post-shock solution, even when a time-independent $\gamma_{\rm cr}$ and $\kappa_{\rm cr}$ are used (as is common in the two-fluid treatment). This EOS can be regarded as a closure across shocks when the length scale of interest is much larger than the shock transition scale. Since the galactic scales are orders of magnitude larger than the CR Larmor radius, below we neglect CR diffusion and focus on the far- upstream and downstream of a shock. In section 6, we further discuss the role of cosmic ray diffusion in determining the small-scale structure of a two-fluid shock.

In the shock rest frame, the mass, momentum, and energy conservation equations for the two-fluid CR-HD model obey the Rankine-Hugoniot jump conditions

[TABLE]

respectively, where $[\mathcal{A}]=\mathcal{A}_{\rm 1}-\mathcal{A}_{\rm 2}$ and the subscripts $1$ and $2$ represent the far-upstream and downstream fluid variables. These equations show that there are three conservation laws (for mass, momentum, and total energy) but four unknowns ($\rho$, ${v}$, $p_{\rm g}$, $p_{\rm cr}$). This constitutes a closure problem across shocks.

In order to highlight this problem, we introduce the compression ratio $\mathcal{R}=\rho_{\rm 2}/\rho_{\rm 1}=v_{\rm 1}/v_{\rm 2}$, and the upstream gas and CR Mach numbers, are respectively defined as

[TABLE]

while the effective upstream Mach number of the composite fluid is $\mathcal{M}_{\rm 1}=(\mathcal{M}_{\rm g,1}^{-2}+\mathcal{M}_{\rm cr,1}^{-2})^{-\frac{1}{2}}$. We then normalize the gas and CR pressures relative to the upstream gas pressure and denote them by $\mathcal{P}_{\rm g,i}=p_{\rm g,i}/p_{\rm g,1}$ and $\mathcal{P}_{\rm cr,i}=p_{\rm cr,i}/p_{\rm g,1}$ respectively, where $i\in 1,2$. The shock jump conditions can then be rewritten as

[TABLE]

[TABLE]

where the normalized upstream CR pressure, $\mathcal{P}_{\rm cr,1}$ is

[TABLE]

Eqs (15) and (16) must be solved for the downstream state ($\mathcal{P}_{\rm g,2}$, $\mathcal{P}_{\rm cr,2}$, and $\mathcal{R}$), given the upstream conditions. This is an underdetermined system since one has more unknowns than equations.

To get a unique solution, we need an additional equation and this is the closure problem. We can fix the downstream CR pressure ($\mathcal{P}_{\rm cr,2}$) in several ways. Here we briefly discuss three possible equation-of-states (EoSs; see Appendix A for details).

• •

$w_{\rm cr}$-EoS: In this case, the downstream CR pressure is set to a fraction $w_{\rm cr}$ of the downstream total pressure, i.e.,

[TABLE]

where $0\leq w_{\rm cr}\leq 1$. The post-shock solution depends on the prescribed value of $w_{\rm cr}$ (Appendix A.1).

• •

$\epsilon_{\rm cr}$-EoS: For this EoS, the downstream CR enthalpy flux is set as a fraction $\epsilon_{\rm cr}$ of the upstream total energy flux, i.e.,

[TABLE]

where $0\leq\epsilon_{\rm cr}\leq 1$ ensures energy conservation. In Appendix A.2.1, we show that wcr-EOS and $\epsilon_{\rm cr}$-EOS are interchangeable and therefore the shock solutions are expected to be similar in these two cases.

• •

Adiabatic-EoS: In this case, the CR entropy across a shock is assumed to be constant, i.e., the downstream CR pressure can be set to

[TABLE]

A method for obtaining post-shock solution for this EoS is discussed in Appendix A.3. This EoS is sometimes used to check consistency between analytic solution and test problem in numerical simulation (see e.g., Pfrommer et al. 2006; Kudoh & Hanawa 2016). However, as we show here, this is only one of the possible EoSs.

Despite the closure problem discussed here, the CR two-fluid equations are mostly used without an EoS imposed across shocks. The impact of the closure problem in numerical implementation, as well as in the resulting solutions, are presented in the next sections.

4 Numerical framework

In this section, we explore a variety of different discretization strategies for the numerical solution of the two-fluid equations. All methods differ essentially in the representation of the coupling term (the non-conservative term in Eq. 4). Let ${\cal L}_{h}()$ and ${\cal L}_{s}()$ be the discrete operators corresponding to the evolution of the homogeneous part of the equations and to the source term alone, respectively. The update step can then be achieved via operator splitting,

[TABLE]

or in a fully unsplit fashion:

[TABLE]

We label these two approaches as OpSplit (Eq. 21) and Unsplit (Eq. 22), respectively. Notice that the operator split method presented here is only first-order accurate in time, but it can be made formally second order accurate using Strang (or alternate) splitting (Strang 1968).

The rationale for choosing ${\cal L}_{h}$ and ${\cal L}_{s}$ depends on the implementation of fluid energy densities and for this reason we classify them into three different options: Option $1$ - “Eg+Ecr”, $2$ - “Et+Ecr”, and $3$ - “Et+Scr”. The differences between various methods are briefly illustrated in Table 1 and the nomenclature of different indices is shown in Fig. 2. Note that since the implementation of mass and momentum equations remain the same in all of the methods, we focus solely on the energy equations.

4.1 Option 1: Eg+Ecr

Here the energies of thermal and CR fluids, coupled by the source term, are evolved separately; i.e., we solve the following equations

[TABLE]

By suitably rewriting the RHS of Eqs. (23) and (24), the coupling term may be implemented either as $p_{\rm cr}{\nabla\cdot\bf{v}}$ (pdv method) or ${\bf v\cdot\nabla}p_{\rm cr}$ (vdp method). Although the former choice (pdv method) is found to be numerically robust and more common, we discuss both implementations separately in the following sections. Note that total energy conservation is ensured by adopting the same discretization for the source term in Eq. (23) and (24).

A. pdv method

In this method, the CR energy flux is defined, following Eq. (24) by ${\bf F}_{\rm cr}=e_{\rm cr}{\bf v}$. Since the source term is $p_{\rm cr}{\nabla}\cdot{\bf v}$, we refer to this method as the pdv method.

A1. OpSplit-pdv method: In many studies, the $p_{\rm cr}{\nabla\cdot\bf{v}}$ term has been implemented using an operator splitting method (e.g., Pfrommer et al. 2006; Sharma et al. 2009; Salem & Bryan 2014; Gupta et al. 2018a). Following Eq. (21), $e_{\rm cr}$ is first evolved along with the other hydrodynamic variables to yield $e^{*}_{\rm cr}$. This value is then used to calculate $p_{\rm cr}{\nabla}\cdot{\bf v}$ at the next step. For a $1^{\rm st}$-order time-stepping scheme this amounts to

[TABLE]

where, in a Godunov-type finite volume scheme, $f^{n}_{i+\frac{1}{2}}$ represents an upwind flux computed with a Riemann solver. Here we employ the HLL Riemann solver (see Chapter $10$ in Toro 2009) to evaluate the fluxes of the conservative terms, while the interface velocity (see Fig. 1) used for the non-conservative term in Eq. (26) is defined as

[TABLE]

In other words, the quantity ${\nabla}\cdot{\bf v}$ is estimated using a cell centered method 222Note that at the time of updating Eqs. (23) and (24) without $p_{\rm cr}{\nabla\cdot\bf{v}}$, the Riemann solver does not have any information about the term $p_{\rm cr}{\nabla}\cdot{\bf v}$. We find that such splitting modifies the effective sound speed to $a=[(a^{2}_{\rm g}+(\gamma_{\rm g}/\gamma_{\rm cr})a^{2}_{\rm cr}]^{\frac{1}{2}}$, instead of $a=[a^{2}_{\rm g}+a^{2}_{\rm cr}]^{\frac{1}{2}}$ (see Appendix B for the derivation). The effective sound speed is required for computing the signal speed/time step of the Riemann solver.. The gas energy equation (23) is updated using a similar discretization as Eqs. (25) and (26). Eq. (25) is readily generalized to second-order temporal acuracy by using a Runge-Kutta method already available in the code while direct extension to curvilinear geometries is thoroughly described in the work by Mignone (2014).

A2. Unsplit-pdv method: In the fully unsplit scheme (Eq. 22), the $p_{\rm cr}{\nabla\cdot\bf{v}}$ term is directly added to the RHS of Eqs. (23) and (24). This yields, for the CR energy equation, the following update:

[TABLE]

Although several options are possible for choosing $p^{n}_{{\rm cr,}i}$ and $v^{n}_{i+\frac{1}{2}}$ (see e.g., Appendix C), one can take advantage of the Riemann solver information to estimate $p^{n}_{{\rm cr,}i}$ and $v^{n}_{\rm i+\frac{1}{2}}$. We find that a robust choice for the cell interface velocity is

[TABLE]

where ${\cal U}_{[m_{x}]}$ and ${\cal U}_{[\rho]}$ are the momentum- and density- state variables obtained from the HLL Riemann solver (Toro 2009).

The CR pressure in Eq. (28) can also be chosen in various ways. We find that a robust selection is

[TABLE]

where $e^{n}_{\rm cr,i+\frac{1}{2}}$ is the state-variable obtained from the HLL Riemann solver, without including the $p_{\rm cr}\nabla\cdot{\bf v}$ term. We label this method as “Eg+Ecr (Unsplit-pdv)".

One of the advantages of choosing the HLL solver is that it does not distinguish between the left- and right- states of the contact wave. This helps avoid a selection ambiguity in the $5^{\rm th}$ eigenvector as discussed in section 2. Thus the choice of HLL Riemann solver is motivated by simplicity and consistency in treating the fluxes and the source terms.

B. vdp method

Alternatively, the coupling term can be implemented in ${\bf v\cdot\nabla}p_{\rm cr}$ form. The evolution of CR energy density is then obtained by rewriting the CR energy equation (Eq. 24) as

[TABLE]

which, at the discrete level, becomes similar to Eq. (25) and (28) by exchanging $\bf v$ and $p_{\rm cr}$.

Although, Eqs. (24) and (31) are mathematically equivalent, the numerical discretizations are not the same. Since the coupling term now is ${\bf v\cdot\nabla}p_{\rm cr}$, we call this method as vdp method. As before, this coupling term can be included via an operator splitting method (OpSplit-vdp) or an unsplit method (Unsplit-vdp).

One may now choose a cell-centered discretization for ${\bf v}_{i}$ and the arithmetic average between the left and right interface values for $p_{{\rm cr},i+\frac{1}{2}}$. However, as it will be shown in section 5, this choice leads to results that depend on the type and order of the spatial reconstruction algorithm. This effect is considerably less pronounced with the previous ($pdv$) method.

4.2 Option 2: Et+Ecr

In this case, we replace Eq. (23) with the total energy equation

[TABLE]

where $e_{\rm t}$ and $p_{\rm t}$ are given by Eqs. (5) and (6).

The total energy density now directly obeys a conservative equation. However, the CR energy equation still contains the coupling term, which can be implemented using the pdv or vdp method as described earlier. After extensive numerical testing, we have found that the numerical results are consistent for various reconstruction schemes only with the $pdv$ method in a fully unsplit fashion, and $v$ and $p_{\rm cr}$ are chosen as given in Eqs. (29) and (30). We thus label this method as “Et+Ecr (Unsplit-pdv)". It is worth mentioning that both “Et+Ecr (Unsplit-pdv)" and “Eg+Ecr (Unsplit-pdv)" give an identical result when Eqs. (29) and (30) are used.

4.3 Option 3: Et+Scr

Kudoh & Hanawa (2016) suggested that the difficulties related to the presence of the coupling term may be avoided if the evolution of CR energy density is obtained from the advection of a passive scalar. They redefined the CR pressure as $p_{\rm cr}\equiv\rho_{\rm cr}^{\gamma_{\rm cr}}$ and used a passive scalar: $\chi=\rho_{\rm cr}/\rho\equiv p_{\rm cr}^{1/\gamma_{\rm cr}}/\rho$ to update the CR energy equation. Under this formalism, Eq. (24) or (31) reads

[TABLE]

The remaining equations maintain the same form as for Option $2$. The advantage of this method is that the two-fluid equations become manifestly conservative, which makes the application of Godunov-type formalism straightforward.

Since Eq. (33) is a tracer equation, the CR entropy $S_{\rm cr}\equiv p_{\rm cr}/\rho^{\gamma_{\rm cr}}$ does not experience a jump across a shock front333Entropy conservation also holds for the Euler equation in smooth regions but the conservative/integral form of the equations must be evolved at shocks for a physical solution.. However, we argue that CRs, like the thermal plasma, should not be adiabatic across a shock. Indeed, CRs are accelerated with non-negligible efficiencies across strong shocks, but this possibility is not allowed by the strict isentropic evolution imposed by the above scalar equation. Isentropic evolution of CRs across shocks is discussed as one of the equation-of-states in section 3 (Eq. 20). The results from this method are labeled by “Et+Scr”.

4.4 Shock detection

In previous sections, we have discussed various possible numerical schemes for implementing CR energy equation. Since the coupling term in this equation involves a spatial derivative, both the numerical and physical two-fluid shock jump condition can be non-unique. While a shock is generally defined as a surface where a sharp transition between upstream and downstream occurs, in numerical simulations, a shock is often broadened over several grid zones due to non-negligible numerical viscosity. To impose an EoS across a shock, it is important to discuss how one can detect shocked zones. For shock detection, we use the following three conditions:

Compressibility: $\nabla\cdot{\bf v}<0$, 2. 2.

Pressure jump: $\left(\nabla p_{\rm t}\cdot\Delta{\bf x}\right)/p_{\rm t,min}\geq\delta_{\rm threshold}$, and 3. 3.

Bypassing spurious waves at contact discontinuity which may fulfil conditions (i) and (ii): $\nabla T\cdot\nabla\rho>0$,

where $T$ is the fluid temperature (see e.g., section $3.1.1$ in Pfrommer et al. 2017; also see Appendix B in Mignone et al. 2012). We find that the choice $\delta_{\rm threshold}=0.5-1$ works quite well for the test problems discussed in this work (see e.g., the top panel of Fig. 8).

5 Test problems & results

In this section we compare the previously presented numerical methods for the solution of selected one- and multi-dimensional benchmarks. We briefly outline the used methods in the following:

• •

Eg+Ecr (OpSplit-pdv) [Eqs 25 - 27]

• •

Et+Ecr (Unsplit-vdp) [Eq. 31],

• •

Et+Ecr (Unsplit-pdv) [Eqs. 28 - 30, and 32], and

• •

Et+Scr (Eq. 33).

As we shall show, the “Et+Ecr (Unsplit-pdv)” method appears to be less sensitive to numerical details (spatial reconstruction, time-stepping, and CFL number). The initial conditions for the selected test problems are listed in Table 2. Computations are performed using either a $1^{\rm st}$-order scheme (Euler time stepping with flat reconstruction), a $2^{\rm nd}$-order scheme (RK2 time stepping with linear reconstruction) or a $3^{\rm rd}$-order scheme (RK3 time stepping with WENO reconstruction). The fiducial CFL number is set to $C_{\rm a}=0.6$ unless otherwise stated.

5.1 Pressure balance mode

This is an important test problem designed to check the evolution of a pressure balance mode. At $t=0$, the left and right states are defined as $(\rho,v_{\rm x,}p_{\rm g},p_{\rm cr})_{\rm L}=(1,1,0.1,0.9)$ and $(\rho,v_{\rm x,},p_{\rm g},p_{\rm cr})_{\rm R}=(1,1,0.9,0.1)$. Since, the total pressure is the same across the interface ($p_{\rm g}+p_{\rm cr}=1$), the profiles should not change in time. The pressure balance test is applicable to any multi-fluid system because of its trivial analytical solution.

To obtain numerical solutions, we use the $2^{\rm nd}$-order scheme and impose periodic boundary conditions. The snapshots of fluid quantities at $t=1$ (i.e., after one advection time across the domain) are shown in Fig. 3. OpSplit-pdv (blue) and Unsplit-pdv (green) curves do not show spurious oscillations. However, Unsplit-vdp (red)444For the Unsplit-vdp method, one can find a suitable combination of $p_{\rm cr}$ and $v$, which does not show these spurious oscillations. However, we find that the same choice does not provide a good solution for other test problems. and Et+Scr (magenta) schemes fail to maintain the pressure balance mode (for both HLL and HLLC Riemann solvers).

Kudoh & Hanawa (2016) showed that the amplitude of the spurious waves in Fig. 3 (for the Et+Scr method) can be reduced by increasing the spatial resolution or including additional numerical diffusion. However, both options are computationally expensive.

5.2 Shock tube A

The initial condition for this test consists of two constant states separated by a discontinuity placed at $x=0$. Left and right states are defined as $(\rho,v_{\rm x,}p_{\rm g},p_{\rm cr})_{\rm L}=(1,0,2,1)$ and $(\rho,v_{\rm x,},p_{\rm g},p_{\rm cr})_{\rm R}=(0.2,0,0.02,0.1)$ respectively. In all cases we employ $1000$ equally spaced zones. The set-up for this problem is identical to Kudoh & Hanawa (2016).

Fig. 4 shows fluid variable profiles at $t=0$ (black dashed curves) and $t=0.1$ (circles) where different colours correspond to the four methods discussed in section 4 (also see section 5.1). The solutions, obtained using $1^{\rm st}$-order reconstruction, are not identical. Blue and green curves - representing the Eg+Ecr (OpSplit-pdv) and Et+Ecr (Unsplit-pdv) methods respectively - look similar although, as we shall see in Fig. 6, the solution obtained with the OpSplit method depends on the details of spatial reconstruction and time-stepping.

We now turn our attention to the top rightmost panel of Fig. 4 showing the CR entropy profile ($p_{\rm cr}/\rho^{4/3}$). In the two-fluid CR-HD model, it is sometimes assumed that CRs are adiabatically compressed in the post-shock region (Pfrommer et al. 2006). This assumption also helps to obtain an analytical solution of the shock tube problem. In numerical simulations this assumption is not automatically fulfilled. In fact, we find that CRs are adiabatically compressed only for the Et+Scr method which is constructed to do so in the first place. In general, the post-shock solution depends on the method one uses.

As noted before, the computed results (for all methods except unsplit-pdv and Et+Scr) depend on numerical details such as reconstruction and the CFL number. While the Et+Scr method is expected to yield robust results at shocks, a constant CR entropy across shocks is not physically justified as discussed earlier. Moreover, this method is unable to maintain pressure balance across a contact discontinuity as shown in section 5.1.

Fig. 6 shows the zoomed-in CR entropy profile from the shock tube problem using four different methods. In each panel, black, green and blue curves correspond to solutions obtained using $1^{\rm st}$, $2^{\rm nd}$ and $3^{\rm rd}$-order schemes. The top left panels indicate that the results are not unique: the profile obtained using the Et+Ecr (Unsplit-vdp, first panel) and Eg+Ecr (OpSplit-pdv, second panel) methods depend on the choice of spatial reconstruction. On the contrary, the Et+Scr and Et+Ecr (Unsplit-pdv) schemes (shown in the $3^{\rm rd}$ and $4^{\rm th}$ columns) are in better agreement555We have also tested the HLLI solver of Dumbser & Balsara (2016) on this problem and find non-uniqueness of the solution..

We now consider the bottom panels of Fig. 6 which, on the other hand, exhibit very similar profiles and have been obtained by fixing the CR pressure behind the shock, as explained in the following. First, a shock detection algorithm (section 4.4) is employed to identify shocked zones. Then we redistribute the shocked zone thermal and CR energies using the parameter

[TABLE]

where

[TABLE]

The redistribution of energy among CRs and the thermal fluid does not change the total energy (in particular $e_{\rm cr}+e_{\rm th}$) of the cell. In the bottom panels of Fig. 6, we have used $w_{\rm cr}=0.5$ to obtain the solutions. While using a parameter $w_{\rm cr}$ to fix the post-shock CR pressure, we notice that the post-shock conditions mainly depend on the total pressure of the shocked-zone nearest to the downstream. Thus, we are able to obtain identical numerical solutions irrespective of the method by imposing a sub-grid closure for CRs at shocks (see e.g., Caprioli & Spitkovsky 2014, Mignone et al. 2018). In section 3, we have discussed some other possible equation-of-states (EoS) which can be used. Although a sub-grid closure is recommended, one would still prefer an implementation that is insensitive to the numerical details. Therefore, we recommend the ‘Et+Ecr (Unsplit-pdv)’ method.

Fig. 6 shows the robustness of the ‘Et+Ecr (Unsplit-pdv)’ method for different spatial reconstructions and time stepping (Appendix C shows insensitivity of this method on the CFL number). As a reference solution, we have performed a simulation using an extremely high spatial resolution. The snapshot of various profiles at $t=0.1$ are shown. Four different solid curves, i.e., black/grey ($1^{\rm st}$ order), green ($2^{\rm nd}$ order), and blue ($3^{\rm rd}$ order) curves show the solutions with different spatial reconstruction. Since the profiles are identical, our implementation of the ‘Et+Ecr (Unsplit-pdv)’ method is numerically robust.

Summarizing our results, we conclude therefore that: (i) the solution of the two-fluid equations in the post-shock region depends on the choice of numerical methods, (ii) all methods give identical results only when the CR pressure behind the shock is imposed ‘by hand’ (e.g., using a subgrid closure at shocks as given in Eq. 34), (iii) the results from our method Et+Ecr (Unsplit-pdv) (Eqns. 29 and 30) numerically robust even in absence of subgrid closure at shocks. For these reasons we discuss, from the next section onwards, the results from our best performing method, “Et + Ecr (Unsplit-pdv)”.

5.3 Shock tube B

We now want to assess the sensitivity of the left- and right- going wave speed estimates in our method. Such estimates are typically required to calculate both the Riemann flux and the cell interface state of the conservative variables for HLL-type solvers. In order to check this, we set up the problem $3$ in Table 2, which leads to the formation of a right-going shock and left-going rarefaction waves.

In the literature, various choices are available for the left- and right- going wave speeds. (see e.g., Chapter $10$ in Toro 2009). The simplest choice is the direct wave speed estimates suggested by Davis (1988), yielding (using the same notations of Fig. 2):

[TABLE]

The present implementation of the Et+Ecr (Unsplit-pdv) shows a spurious feature near the opening of the rarefaction wave (red curves in Fig. 7). We address this issue by adjusting the spatial width of the Riemann fan through a redefinition of the right- and left- wave speeds as

[TABLE]

respectively, and choose ${\phi}=1.1$. Note that, when ${\phi}=1$ then Eq. (37) reduces to Eq. (36). This trick basically increases the width of the mixed state, which becomes $2\,\phi\,a$ at $t+\Delta t$, instead of $2\,a$ (here $a=(a_{\rm g}^{2}+a_{\rm cr}^{2})^{\frac{1}{2}}$ is the sound speed of composite fluid). Increasing the width of the mixed state is justified by the fact that during the solution of the Riemann problem we do not include the coupling term which changes the effective signal speed (see e.g., Appendix B). Although in the Unsplit-pdv method we update the CR energy equation (Eq. 28) in a single step, the required sound wave speed may lie between $[a_{\rm g}^{2}+a_{\rm cr}^{2}]^{\frac{1}{2}}$ and $[a_{\rm g}^{2}+(\gamma_{\rm g}/\gamma_{\rm cr})a_{\rm cr}^{2}]^{\frac{1}{2}}$. Note that, the factor $\phi$ introduced here is obtained from our numerical experiments, with no rigorous proof.

Fig. 7 shows different pressure profiles at $t=10^{-4}$ obtained with $\phi=1$ (red) and $\phi=1.1$ (green). For both solutions, in the post-shock region, the ratio $p_{\rm cr}/(p_{\rm cr}+p_{\rm g})\simeq 0.03$. Solutions are smooth only for $\phi=1.1$, implying the necessity of a slightly higher wave speed estimate in order to obtain a robust solution.

5.4 Sedov-Taylor Blast Wave

We have extended our implementation of CR-HD equations to multi-dimensions and different geometries in the PLUTO code. We perform the standard blast wave problem in $3$D Cartesian and $1$D spherical geometries. The numerical set-up is discussed below.

At $t=0$, we create high pressure in a small region by setting $p_{\rm g}=(\gamma_{\rm g}-1)E/\Delta V$ where $\Delta V=4/3\pi\,r_{\rm inj}^{3}$ is the small volume, $r_{\rm inj}=0.01$ (in code units), and $E=10^{51}$ erg. In the rest of the computational domain, density, velocity, and pressures are set to $1$, [math], and $60$ respectively. We normalize density, velocity, and length to $\rho_{u}=m{\rm{}_{H}\,cm^{-3}}$ ($1.67\times 10^{-24}$ g cm*-3*), $v_{u}=10^{5}\,{\rm cm\,s^{-1}}$ and $L_{u}=3.086\times 10^{18}\,{\rm cm}$ ($1$ parsec) respectively. For $1$D spherical geometry, the inner radial boundary at $r=0.001$ is set to be reflective while for $3$D Cartesian, we set all boundaries to outflow. For the 3D run, we employ $200$ equally spaced zones along $x,y$, and $z$ directions, i.e., $N_{\rm x}N_{\rm y}N_{\rm z}=200^{3}$ (but we also perform a low resolution run using $N_{\rm x}N_{\rm y}N_{\rm z}=100^{3}$ to check convergence), and for the $1$D run we use $200$ equally spaced zones. Simulations are performed using a CFL number $C_{\rm a}=0.2$ for both geometries. Since we do not inject mass initially, the Sedov-Taylor phase starts right at the beginning of the shock evolution. The initial injection is purely thermal, but we inject CRs at the shock as follows.

First, we identify shocked zones using the detection algorithm (see section 4.4) and then inject CRs in these zones using the parameter $w_{\rm cr}$ (Eq. 34) as done previously in 1D. We set $w_{\rm cr}=0.5$, i.e., equipartition between thermal and CR pressures in the post-shock zones. The results from $3$D (Cartesian) and $1$D (spherical) runs are shown in Fig. 8. In the top panels we show the snapshots of density, thermal and CR pressures in the $z=0$ plane at $t=2\times 10^{-6}$ (code unit) obtained from our $3$D run. The second panel (top row) shows the zoomed-in view of the shocked zones where CR energy is injected. In the bottom panels we compare the results between the $3$D case and the $1$D runs, showing good agreement. We, therefore, conclude that our implementation (i.e., Et+Ecr (Unsplit-pdv) method) is well suited for multi-dimensional calculations.

6 Discussion

The results from the previous section have shown that the numerical solution of the two-fluid equations depends on the discretization of the coupling term and that the results become insensitive to the choice of the method only when the post-shock CR pressure is fixed by the parameter $w_{\rm cr}$ (Fig. 6), which provides a possible closure to the problem discussed in section 3. Below we further discuss to what extent the CR diffusion term that we have neglected until now can affect our conclusions, and also present the limitations and broader implications of this work.

6.1 Shock structure with CR diffusion

Several works on two-fluid shocks (Drury & Voelk 1981; Becker & Kazanas 2001) have considered the impact of CR diffusion on the shock structure. A non-negligible CR diffusion implies that the CR pressure (unlike the gas pressure) remains continuous across the shock transition (length scale $\sim\kappa_{\rm cr}/v_{\rm 1}$). Since the spatial width of the shock is now broadened, some additional equations can be formulated between far- upstream and downstream of a shock, which are known as precursor-EoSs (Drury & Voelk 1981; Voelk et al. 1984; Jun et al. 1994). The precursor-EoS connects the upstream and downstream and thereby allows to investigate the shock structure, and it can be seen as a closure to two-fluid equations. Although this method has been promising to study shock modification due to CR acceleration, some features are worth highlighting.

In CR-HD, two distinct shock structures are possible: a discontinuous, gas-mediated sub-shock with CRs diffusing ahead of the shock; and a CR-dominated smooth “shock” (across which the gas properties vary smoothly) for large Mach numbers, i.e., $\gtrsim 12$, see Fig. 2 in Becker & Kazanas 2001, or the right panels in Fig. 5 of Gupta et al. 2018a for an astrophysical realization in a wind-driven shock. In the latter case, most of the upstream kinetic energy flux goes into CR acceleration rather than heating the thermal plasma, and the astrophysical implications of this are enormous (see e.g., Gupta et al. 2018b). However, the physical existence of such shocks is yet to be established by kinetic plasma simulations. In some parameter space, multiple solutions are predicted for the same set of upstream conditions and the results are found to be sensitive to the upstream variables close to the shock surface (for a brief discussion, see section $6$ in Becker & Kazanas 2001).

Apart from the multiplicity (non-unique) of the post-shock solutions, the choice of the CR diffusion coefficient is crucial. Near a shock, magnetic fluctuations due the plasma instabilities (Kulsrud & Pearce 1969; Bell 2004) can significanly modify $\kappa_{\rm cr}$ from those used in large scale flows (Schroer et al. 2020). Therefore although in theory, CR diffusion allows to find shock structure, the results on the shock structure are inconclusive and further investigation is needed. The results from kinetic studies, irrespective of whether a globally smooth shock forms or not, can always be implemented as a subgrid EoS.

6.2 Limitations and broader implications

Non-unique numerical results at the shocks that have been discussed in this work is ultimately connected to the missing microphysics in fluid approximation (Hoyle 1960). Both fluid and kinetic studies showed that CRs drive a zoo of plasma instabilities that are responsible for CR acceleration and eventualy build-up the post-shock CR pressure (Dorfi & Drury 1985; Drury & Falle 1986; Ryu et al. 1993; for a review see Marcowith et al. 2016). Since the majority of these instabilities occurring on the CR Larmor radius scale depend on the CR energy flux, neglecting them far away from the shock is justified as long as their effects are incorporated in the subgrid EoS connecting the far-upstream and downstream regions. However, to obtain a detailed shock structure in fluid theory, we must include a self-consistent (derived from kinetic considerations) CR diffusion.

Note that the CR injection model presented in this work (see e.g., Fig. 8) is different from those suggesting injection at the gas sub-shocks (see e.g., Falle & Giddings 1987; Kang & Jones 1990; Saito et al. 2013). The sub-shock injection still relies on the fluid approximation and the results can be affected by the choice of the diffusion coefficient and other missing kinetic processes at shocks. Moreover, injection at the gas sub-shock requires resolving the shock transition region (Hin Navin Tsung et al. 2020) which is computationally infeasible in galaxy-scale simulations. From this perspective, our method may seem a simplification, but we believe that in absence of a self-consistent model that takes into account the necessary shock microphysics, the proposed method can give better control on shock solutions and their impact on the large scale flow.

Although we recommend using such closures at shocks, it is worth mentioning that this requires shock identification on the computational grid. We suggest the “Et+Ecr (Unsplit-pdv)” method (section 4.2) for a robust numerical solution, irrespective of the use of a physically motivated subgrid closure at shocks.

The main limitation of our method is that it does not allow investigating the time evolution of shock structure self-consistently. This is not the aim of studies large scale impact of CRs (length scale $\gg$ shock transition width). By defining the CR injection as a function of time, the shock evolution may be inspected. A physically reliable value of $w_{\rm cr}$ (or an equivalent parameter; e.g., see section 3 depends on upstream parameters (e.g., Mach number and magnetic field orientation w.r.t. shock normal) and it needs to be prescribed using Particle-In-Cell simulations (see e.g., Caprioli & Spitkovsky 2014). In MHD, the injection of CRs will also depend on the nature of the shock (e.g., fast versus slow) and detailed calibration with kinetic plasma simulations will be needed (Bret 2020).

7 Summary

In this work we have shown that the two-fluid CR-HD model suffers from a non-uniqueness problem at shocks. It contains three conservation laws and one additional equation (for the CR pressure) which cannot be cast in a satisfactory conservative form and causes difficulties in its numerical implementation. A unique shock jump condition is possible only if one makes an additional assumption. The steady-state shock structures can be predicted by assuming a suitable downstream CR pressure/energy (as discussed in section 3). There is a degeneracy between the gas and CR pressures (see last paragraph of section 2) because of which the solutions may depend on numerics. Without fixing the fraction of upstream energy transformed into CR energy and simply relying on the numerical discretization of the non-conservative exchange term involving derivatives $p_{\rm cr}\nabla\cdot{\bf v}$ (or equivalently, ${\bf v}\cdot\nabla p_{\rm cr}$), makes the solutions of the two-fluid equations across shocks depend on the details of numerics (e.g., spatial reconstruction, time-stepping, and even the CFL number). In this work, we have investigated numerical implementation of the two-fluid CR-HD equations as applied to large scale simulations. Our findings are summarized as follows:

• •

Numerical solution of the two-fluid equations depends on implementation of the coupling term ($p_{\rm cr}{\nabla\cdot\bf{v}}$ or ${\bf v\cdot\nabla}p_{\rm cr}$). We show that the different discretizations do not show an identical solution (Figs. 4 and 6). This is because of the non-negligible and non-unique contribution of the source term involving a derivative at the shock.

• •

We suggest a method (“Et+Ecr (Unsplit-pdv)”, section 4.2) for which the solutions are robust to the choice of spatial reconstruction, time-stepping, and the CFL number (see Fig. 6). In order to ensure that the characteristic speed remains within the fastest signal propagation speed, this method demands a slightly higher signal speed than the standard estimate of the two-fluid sound speed (Fig. 7).

• •

We show that all methods give an identical solution only when the CR pressure across the shock is fixed by an imposed equation of state. This can be done, for example, by specifying CR pressure in the shocked zones (as done in Fig. 6; also see section 3). A physically realistic implementation of the CR pressure across a shock is possible by calibrating with kinetic simulations using different upstream parameters. In this approach, a shock must be properly identified. We suggest the “Et+Ecr (Unsplit-pdv)” method for a robust solution irrespective of whether a subgrid closure is used at shocks.

In summary, this work highlights the critical aspects of the two-fluid CR-HD equations. Although here we have not discussed the CR-MHD system, these problems are also present there. The implementation of CR-MHD will be discussed in a future work.

Acknowledgements

We thank Dinshaw Balsara, Damiano Caprioli, Michael Dumbser, Yuki Kudoh, and Christoph Pfrommer for helpful discussions. SG acknowledges CSIR, India for the SPM Fellowship. SG thanks the postdoctoral research program of the University of Chicago, where some of this work was done to address referee’s suggestions. PS thanks the Humboldt Foundation for supporting his sabbatical at MPA where some of this work was done. PS acknowledges the Department of Science and Technology for a Swarnajayanti Fellowship (DST/SJF/PSA-03/2016-17). We thank our reviewer, Luke Drury, for constructive comments on our manuscript.

Data Availability

The data and code associated with the paper can be made available on request to the corresponding author.

Appendix A Details of various EoSs

Here we present the detailed calculations for finding shock solution for a given EOS outlined in section 3.

A.1 wcr-EoS

Using Eq. (18), we substitute $\mathcal{P}_{\rm g,2}$ in Eqs. (15) and (16). After combining them, we obtain a quadratic equation:

[TABLE]

[TABLE]

Note that $w_{\rm cr}$ is a parameter specified in the downstream, and $\mathcal{P}_{\rm cr,1}$ (Eq. 17) and $\mathcal{M}_{\rm 1}(\mathcal{M}_{\rm cr,1},\mathcal{M}_{\rm cr,1})$ are upstream parameters. The solutions of Eq. (38) is given by

[TABLE]

Since $\mathcal{R}_{\rm-}\leq 1$, $\mathcal{R}_{\rm+}$ represents the physical solution and its values as a function of upsteam Mach numbers for four different $w_{\rm cr}$ are shown in Fig 10. The figure indicates that the compression ratio increases with $w_{\rm cr}$. It can also be shown analytically:

[TABLE]

Taking $\gamma_{\rm g}=5/3$ and $\gamma_{\rm cr}=4/3$, we obtain $\mathcal{R}_{\rm+}=4+3w_{\rm cr}$, which shows that the compression ratio approaches $7$ when $w_{\rm cr}\rightarrow 1$. Note that, CR diffusion can increase the compression ratio even for a smaller $w_{\rm cr}$ due to diffusive loss of CR energy (section $4.2.1$ in Gupta et al. 2018a).

A.2 $\epsilon_{\rm cr}$-EoS

In this case, we follow a similar approach as done in Appendix A.1. We first normalize the gas and CR pressures of Eq. (19) w.r.t. the upstream gas pressure and find

[TABLE]

Then we substitute $\mathcal{P}_{\rm cr,2}$ (Eq. 44) in Eqs. (15) and (16), which results in a quadratic equation identical to Eq (38), where

[TABLE]

Therefore, the solutions can be found by solving Eq. (42). We find that the qualitative nature of the solutions remains similar to Fig. 10, i.e., the compression ratio increases with $\epsilon_{\rm cr}$.

A.2.1 Relation between wcr and $\epsilon_{\rm cr}$

The relation between wcr-EoS and $\epsilon_{\rm cr}$-EoS can be obtained as follows. From Eq. (18), we note that the value of $w_{\rm cr}$ depends on the downstream CR and gas pressures. In $\epsilon_{\rm cr}$-EOS, the downstream CR pressure ($\mathcal{P}_{\rm cr,2}$) is already parameterized by $\epsilon_{\rm cr}$, which is given in Eq. (44). The remaining quantity, $\mathcal{P}_{\rm g,2}$, can be obtained using the momentum conservation equation, Eq. (15), which gives

[TABLE]

It is worth mentioning that the only constrain on the above equation is $\mathcal{P}_{\rm g,2}>0$, implying that the downstream CR pressure $\mathcal{P}_{\rm cr,2}<\left\{\left(\gamma_{\rm g}+\gamma_{\rm cr}\mathcal{P}_{\rm cr,1}\right)\mathcal{M}_{\rm 1}^{2}\left(1-\frac{1}{\mathcal{R}}\right)+1+\mathcal{P}_{\rm cr,1}\right\}$. Since we now have $\mathcal{P}_{\rm g,2}$ and $\mathcal{P}_{\rm cr,2}$, we can obtain $w_{\rm cr}$ as a function of $\mathcal{R}$, $\mathcal{M}_{\rm 1}$, $\mathcal{P}_{\rm cr,1}$, and $\epsilon_{\rm cr}$ where $\mathcal{R}(\mathcal{M}_{\rm 1},\mathcal{P}_{\rm cr,1},\epsilon_{\rm cr})$ can be calculated from the solution of Eq. (42). Therefore, the parameter $w_{\rm cr}$ becomes a function of $\mathcal{M}_{\rm g,1}$, $\mathcal{M}_{\rm cr,1}$, and $\epsilon_{\rm cr}$. The value of $w_{\rm cr}$ for various upstream Mach numbers and $\epsilon_{\rm cr}$ are shown in Fig. 10. In the limit $\mathcal{M}_{\rm 1}\rightarrow\infty$, it can be shown that

[TABLE]

where we have taken $\mathcal{R}_{\rm+}$ from Eq. (43) and the constrain $0\leq w_{\rm cr}\leq 1$ implies that $w_{\rm cr,+}$ is the physical solution.

A.3 Adiabatic-EoS

To find the compression ratio, we substitute $\mathcal{P}_{\rm cr,2}$ from Eq. (20) to Eq. (15) and obtain $\mathcal{P}_{\rm g,2}$ as a function of $\mathcal{R}$, $\mathcal{P}_{\rm cr,1}$, and $\mathcal{M}_{\rm 1}$. Replacing this $\mathcal{P}_{\rm g,2}$ in Eq. (16) we obtain

[TABLE]

[TABLE]

Eq. (50) can be solved using a standard root-finding scheme (e.g., the Newton-Raphson method). After obtaining the solution for $\mathcal{R}$, we can calculate $w_{\rm cr}$ (Eq. 18) and the results are shown in Fig. 11. The left panel of Fig. 11 displays the compression ratio as a function of $\mathcal{M}_{\rm g,1}$ and $\mathcal{M}_{\rm cr,1}$, which shows that the compression ratio $\leq 4$. The right panel shows $w_{\rm cr}$ as a function of $\mathcal{M}_{\rm g,1}$ and $\mathcal{M}_{\rm cr,1}$. It shows that if the upstream CR pressure is very small compared to the upstream ram pressure (i.e., when $\mathcal{M}_{\rm cr,1}\gg 1$) then the post-shock CR pressure is almost negligible compared to the gas pressure. For a large upstream CR pressure (${\cal M}_{\rm cr,1}\sim 1$), the value of $w_{\rm cr}\rightarrow 1$, i.e., the downstream is mostly dominated by CR pressure. However, even in this regime, the compression ratio $\leq 4$, which can be seen by comparing the left and right panels. These results are easy to understand intuitively because the effective Mach number $\mathcal{M_{\rm 1}}\sim 1$; i.e., the shock is too weak to show a noticeable change in CR pressure.

Appendix B Eigenvalues without the coupling term

To obtain the HLL flux at left- and right-interface of the computational cell, we need to provide an estimate of the signal speed to the solver. Since in the operator splitting method the solver does not have any information of the $p_{\rm cr}\nabla\cdot{\bf v}$ (or ${\bf v}\cdot\nabla p_{\rm cr}$) coupling term, the signal speed can be different from the actual speed of the complete system. This is illustrated as follows. Consider that the source term in Eq. (3), $p_{\rm cr}\nabla\cdot{\bf v}$, is absent. In 1D Cartesian geometry, the Jacobian matrix without the coupling and diffusion terms is found to be

[TABLE]

Eigenvalues of this matrix are $\lambda=v-a_{\rm eff},v,v,v+a_{\rm eff}$ where $a_{\rm eff}=[a^{2}_{\rm g}+(\gamma_{\rm g}/\gamma_{\rm cr})a^{2}_{\rm cr}]^{\frac{1}{2}}$. The methods involving the $pdv$ source term give the same set of eigenvalues. As $\gamma_{\rm g}>\gamma_{\rm cr}$, we note that the effective propagation speed $a_{\rm eff}$ is slightly larger than the actual sound speed of the composite fluid, which is $[a^{2}_{\rm g}+a^{2}_{\rm cr}]^{\frac{1}{2}}$. Without using the effective sound speed as an estimate of the maximum signal propagation speed, one may find spurious oscillations in the solution. Such an example has been discussed in section 5.3. The effective sound speed ($a_{\rm eff}$) with vdp splitting is also different from $[a_{\rm g}^{2}+a_{\rm cr}^{2}]^{1/2}$, and we advice caution even in this case for the signal speed estimate.

Appendix C Dependences on the CFL number

Other than for the Et+Scr method, the coupling term has to be implemented in the two-fluid equations. In order to implement this term, we have to choose $p_{\rm cr}$ and $v$ at the cell center/interface, which can be chosen in various ways. We find that all of these choices do not produce a unique/consistent result. To explicitly show this, we present the result from the shock tube A (problem $2$ in Table 2) using three possible combinations of $p_{\rm cr}$ and $v$ in Eq. (28). These choices are given below.

[TABLE]

Note that the possible choices are not limited to the above three combinations. Despite the similarity of these choices, we find that the results are not identical and also depend on the CFL number. Fig. 12 shows the solutions for two different CFL numbers $0.1$ (solid curves) and $0.6$ (dotted curves) respectively. Each colour represents different possible combinations of $p_{\rm cr}$ and $v$ given above. We have tried several other possibilities and found that all give different result. This experiment shows that our choice of $v$ and $p_{\rm cr}$ (i.e., Eqs. 29 and 30) provides a robust result (solid/dotted black curves), i.e., solutions are robust to the choice of spatial reconstructions (Fig. 6), time stepping (Fig. 6), and the CFL number (Fig. 12).