{'fraction_non_alphanumeric': 0.039308015149964175, 'fraction_numerical': 0.009253761899887398, 'mean_word_length': 4.251693366304698, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Glimm-Godunov's Method for Cosmic-ray-hydrodynamics

A numerical method for integrating the equations describing a dynamically coupled system made of a fluid and cosmic-rays is developed. In smooth flows the effect of CR pressure is accounted for by modification of the characteristic equations and the energy exchange between cosmic-rays and the fluid, due to diffusive processes in configuration and momentum space, is modeled with a flux conserving method. Provided the shock acceleration efficiency as a function of the upstream conditions and shock Mach number, we show that the Riemann solver can be modified to take into account the cosmic-ray mediation without having to resolve the cosmic-ray induced substructure. Shocks are advanced with Glimm's method which preserves their discontinuous character without any smearing, thus allowing to maintain self-consistency in the shock solutions. In smooth flows either Glimm's or a higher order Godunov's method can be applied, with the latter producing better results when approximations are introduced in the Riemann solver.


Introduction
We wish to formulate a numerical method to solve a system of equations characterizing a fluid that is dynamically coupled to suprathermal particles through the exchange of momentum and energy. Such conditions occur commonly in astrophysical plasmas. The fluid system is an ordinary nonrelativistic gas described by the following modified equations of hydrodynamics : where (ρ, u, P g , e g ) indicate the gas density, velocity, pressure and specific energy respectively; i, d index the spatial components and summation over repeated indexes is assumed; δ id is the Kronecker's delta. The gas total specific energy is, e g = u 2 /2 + e th , and a γ-law equation of state is assumed so that the gas specific internal energy is related to the gas pressure through e th = P g /ρ(γ g − 1). The inhomogeneous terms proportional to ∂P c /∂x on the right hand side of Eq. (2)-(3) account for the effects of the suprathermal pressure. Σ is a source term describing the transfer of energy between the fluid and the suprathermal component. This may be due to, e.g., particle acceleration processes at the expenses of the fluid energy or, conversely, energy losses from the suprathermal particles that end up heating the fluid.
As for the suprathermal component we consider cosmic-ray (heretofore CR) particles described by a distribution function, f (x, p, t), which depends upon a spatial, a momentum and a temporal coordinate. In what follows p is in units of 'm c c', with m c the CR particle mass, and the normalization of f is such that the number density of particles with momentum between p and p + dp is dn c = 4π p 2 f dp. In addition, f is assumed to be isotropic in momentum space and evolves according to the following diffusion-convection equation [35] ∂f ∂t The second and third term on the left hand side of the above equation represent, respectively, spatial advection and diffusion with a coefficient, κ(x, p). The first term on the right hand side accounts for adiabatic effects and, b ℓ (p) ≡ −(dp/dt) loss , describes the particle momentum change due to energy losses associated with mechanical and radiative processes. In addition, D p (p) is the diffusion coefficient in momentum space. The CR pressure in Eq. (2) and (3) is then defined through the distribution function, f , as where p min , p max are the minimum and maximum CR momenta, respectively. More specifically, the momentum p min marks (somewhat loosely) the transition between the thermal and nonthermal components and p max is the maximum momentum the particles can achieve and still be confined inside the system.
The CR energy and adiabatic index are given by The evolution of the CR energy is obtained from Eq. (4) and reads The first line of Eq. (8) refers to the effects of advection, diffusion (with an energy averaged diffusion coefficient κ ) and adiabatic compression, and the second line to energy losses/gains introduced above. The surface terms on the third line describe changes in the CR energy due to the flux of particles across the low and high boundaries in momentum space. The first of these two surface terms is typically negligible whereas the second is important in case of efficient shock acceleration and can cause significant energy losses in the system. Finally the source term in Eq. (p 2 + 1) 1 2 dp −4πm c c 2 p 2 f 1 3 where b mℓ (p) now includes mechanical losses only (i.e. radiative losses are excluded).
The spatial diffusion coefficient introduces a physical scale characterizing the particles mean free path due to diffusion, i.e. λ mf p ∼ κ(p)/v(p), where v(p) ∼ c, the velocity of a particle of momentum p, is of order the speed of light. In the following we shall distinguish two different regimes of application of the equations (1)-(4), namely smooth flows and shock waves. The reason for doing so is that for astrophysical systems, λ mf p ≪ λ system , and the entire dynamic range of scales cannot be resolved with currently available computers. However, while on the scales that can be resolved by simulations diffusion in smooth flows can be safely assumed to become either slow or negligible, this is not the case for shocks.
In smooth flows (and on large enough scales, λ ≫ λ mf p ) the presence of CRs enhances the propagation speed of sound waves but simultaneously causes a damping of their amplitude due to CR diffusion [32,34]. In addition energy is exchanged non adiabatically between the thermal and nonthermal components according to Eq. (9). These effects arise from diffusive processes and as long as the relevant transport coefficients, κ and D p , are defined correctly, they can be properly modeled numerically with schemes available in the literature.
Around shocks, however, the situation is more complicated because the diffusion process gives rise to an efficient mechanism for transferring energy from the flow to the particles. This topic is discussed in detail in several review articles [9,5]. Here we emphasize two basic and related points relevant for the present discussion. Firstly, the dissipation of energy into CRs changes the value of the total pressure generated by the shock dissipation mechanism due to the different thermodynamic properties of gas and CRs. In addition, as illustrated above (cf. Eq. [8]) the escape of high energy particles upstream of the shock allows for energy to be removed from the system. This can reduce the pressure support in the downstream region, allowing for compression ratios higher than the hydrodynamic limit. Finally, the CR pressure gradient produced by the CR particles diffusing upstream decelerates the flow approaching the shock front. As a result the velocity structure is not a sharp transition anymore but is broadened up to scales of order the diffusive scale length of the most energetic CR particles. This is λ κ (p max ) = κ(p max )/u shock , where u shock is the shock speed. This effect creates the so called shock precursor where the upstream gas is adiabatically compressed before being shocked. Thus, even though in a numerical calculation the CRs do not diffuse out of the resolution element during a timestep, shock acceleration can modify significantly the shock jump conditions with respect to the simple fluid case (see Section 3).
There are, therefore, at least two different limits of interest for solving the equations (1)-(4) in the case of shocks. One which focuses on the study of the diffusive shock acceleration process itself. In this case one requires: (i) a kinetic approach in which the evolution of the distribution function in momentum space given by Eq. (4) is followed accurately [2,23]; (ii) enough spatial resolution to properly resolve the full range of relevant scales that enter the problem, from the thickness of the shock, ℓ shock , to the diffusive scale length of the highest energy CR particles, λ κ (p max ). A number of codes, with different levels of sophistication, employing various numerical methods and devoted to this type of approach have been developed (e.g. [13,1,19,11,4,14,22,16]). Sometime and to various extents they also include the processes that regulate the diffusive properties of the medium (e.g. wave amplification and damping), a key factor for the efficiency of the acceleration mechanism. By necessity, they focus on a very narrow region around the shock, of order λ κ (p max ). Complemented by analytic studies [12,25,23,24], among their ultimate goals is to investigate, as a function of the upstream gas conditions, U − , and the shock Mach number, M, the downstream CR distribution function, f + (p), and the efficiency of the shock acceleration process, η. Here and in the following this quantity is defined as the fraction of total pressure upstream of the shock that is converted into downstream CR pressure, P + c : where P + c is related to f + through Eq. (5).
There is then the opposite limit to the one described above, which focuses on the dynamical effects of CRs in smooth flows and shocks, for systems of size λ system ≫ λ κ (p max ). We can refer to it as the astrophysical limit. This approach is more application oriented [17,31,30,28,29,36,18] and does not aim at studying f + (p; U − , M) or η(U − , M). In this paper we attempt to design a numerical method that serves this purpose, without having to resolve scales of order ℓ shock . In doing so, we seek to eliminate all but the essential information about the CR distribution function in momentum space, so that a fluid-like description is approached. Information at the kinetic level must, however, be preserved in two parts of the formulation: (a) when computing shock solutions and (b) when computing the time-evolution of the CR pressure. The first requirement is set because, as already pointed out, a correct shock solution can only be obtained with a fully kinetic approach [2,23]. This means that the only way to meaningfully include the effects of CRs acceleration in a fluid-like approach is to assume that f + (p; U − , M) and the shock acceleration efficiency, η(U − , M), are provided independently (e.g. from kinetic models) as part of the input. When the acceleration efficiency is high and the contribution from the highest energy particles dominate the CR energy and pressure, the flux of CR particles across p max must also be specified [26]. The second requirement stems from the fact that dP c /dt is basically the result of energy losses/gains and diffusion of the CR particles. Since these processes are strongly momentum dependent, dP c /dt is bound to be different for different shapes of the CR distribution function and one ought to be able to account for this. Note that the two-fluid approximation alone, in which Eq. (4) is integrated in momentum space to derive an equation for the time evolution of P c , is not sufficient for these two purposes.
In order to properly evolve P c , we divide momentum space in a set of N p (∼ 10) coarse kinetic volumes or momentum bins and integrate Eq. (4) within the boundaries of each of them. This provides an equation for the evolution of the number density of CRs within each bin and is a cost-effective way of following the change of shape of the CR distribution function resulting from the momentum dependent CR processes mentioned above. This essentially works because the CR distribution, f , is typically a smooth power law with a slowly varying slope as a function of momentum. This fact provides a natural way for describing f (p) within each bin [27,17], which compensates for the coarseness of the discretization of momentum space.
In addition, we show that once the shock acceleration efficiency, η, is specified it is possible to account for the modifications induced by the CRs on the hydrodynamic shock solution, even though the structure of the shocks remains completely unresolved. This is achieved by solving a slightly more complicated Riemann problem, after proper modification of the definition of the nonlinear waves that appear in it. This shock solution can still be thought of in terms of a two-fluid model description [10] but with the important difference that, among the family of admissible solutions [2], we select the one demanded by the (explicitly) adopted shock acceleration model. This ensures a self consistent description of the CR-hydrodynamic system.
In order for this to work, however, it is essential that the shock discontinuity does not spread as a result of numerical dissipation. This is because in general the dissipation of CR energy at a shock is a nonlinear function of the full jump conditions. Therefore, if the shock is artificially spread over a few zones the sum of the CR energy generated at each numerical subshock will in general not be the same as that predicted by the model for the full jump conditions.
A suitable hydrodynamic method for our purpose is the one originally proposed by Glimm [15]. Introduced as part of a constructive proof of existence of solution to nonlinear hyperbolic conservation laws it was turned into an effective numerical scheme for hydrodynamics by Chorin [6,7]. Glimm's method maintains unsmoothed all the sharp features that are present in the flow. In particular, shocks remain unsmoothed jumps as they propagate across the grid. This allows us to maintain self-consistency in the shock solution. The limitation in using Glimm's method is that at the moment its multidimensional extension does not work properly at shocks [8]. Thus here we focus on a one dimensional algorithm and leave its generalization to more than one dimension for future work.
In smooth flows either Glimm's or Godunov's method can be applied. Therefore it is possible to define a hybrid scheme where Glimm's method is applied at shocks and Godunov's method in smooth parts of the flow. In either case, (in smooth flows) the effects due to the CR pressure are included by proper modification of the characteristic analysis.
Note that recently a method has been proposed in [18,33] to include the dynamical effect of CR pressure on the hydrodynamics. That approach consists effectively of a two fluid model in which the generation of CRs at shocks is treated as an explicit source term, similar to the scheme in [27]. However, a hydro scheme that includes self-consistently the kinetic effects on both the shock substructure and the CR pressure evolution, in the sense mentioned in (a) and (b) above, is still lacking.
This paper is structured as follows. In Section 2 we describe the discretization of momentum space and compute the fluxes due to energy losses and diffusion in that dimension. In Section 3 we discuss the effects of CRs on the structure of the flow and define a modified Riemann problem to include the effects of CRs both in smooth flows and at shocks. In Section 4.1 and 4.2 we describe the implementations of Glimm's method and of a hybrid Glimm-Godunov's method, respectively. Tests follow in Section 5 and conclusions are presented in Section 6.

Diffusion-Convection Equation
In order to formulate a finite-kinetic-volume description of the diffusion-convection equation (4), we divide momentum space into N p logarithmically spaced bins, each with boundaries ), is taken as constant (although this is not necessary). We then follow the evolution of the CR number density associated with each bin, namely For a piecewise power-law distribution function we have where f 0j and q j are the normalization and logarithmic slope for f in the  th momentum bin. Once the set, {n p j ; 0 ≤ j < N p }, is defined and the boundary conditions for the slopes q −1 and q Np are provided, we can compute the set of slopes, {q j ; 0 ≤ j < N p }, and normalizations, {f 0j ; 0 ≤ j < N p }, as follows [17]. For each bin, i, we use Eq. (13) with j = i, i ± 1 and further assume: (a) that the spectral curvature, ∂q/∂ ln p, remains constant across adjacent bins; (b) that f (p), as given in Eq. (12), is continuous across the bins boundaries, p i− 1 2 and p i+ 1

2
. This provides six equations for six variables that can be efficiently solved with an iterative method. This procedure is applied for each bin, and allows to reconstruct the set {f j (p)} from {n p j }. In the following we will use both n p j and f j (p) assuming that we can reconstruct the latter from the former through the procedure just described 1 .
The equation describing the evolution of, n p j , is obtained by integrating Eq. (4), multiplied by a factor 4π p 2 , within the interval, . This leads to where∇ p is the undivided (one-dimensional) gradient in momentum space and we have introduced In writing Eq. (14)- (16) we have emphasized that both advection and diffusion terms along the momentum coordinate can be cast in conservative form, in analogy to the corresponding terms in configuration space. This allows us to adopt a Godunov-like scheme for the numerical integration of those terms. However, we place them on the right hand side of the Eq. (14) because they will effectively be treated as source terms of n p j . Finally, on the right hand side of Eq. (14), we have added a source term, J j , which represents the rate of production of CR particles due to the diffusive shock acceleration mechanism. We do this only pro forma: because we always use Glimm's method to advance shocks, the Riemann solver effectively subsumes the role of the term J j . In other words, J j will not be treated as an explicit source term but as an implicit part of the shock solution computed through the Riemann solver. Thus, it will be sufficient for our purposes to only formulate a prescription for the postshock values of the set {n p j }, which we do in Section 3.3.
1 When spectral curvature is important, the following alternative approach first proposed in [27] can be used instead: for each bin, in addition to the number density of CR particles, one also follows the energy density, ǫ p j . For each bin the definitions of n p j and ǫ p j provide two equations which can be readily solved for the two unknowns (f 0j , q j ). While the equation for n p j is derived below (cf Eq. 14), the one for ǫ p j is obtained analogously by multiplying Eq. (4) by a factor, 4π p 2 [(p 2 + 1) . See also [16] for the effectiveness of this method. Here we take the simpler approach, however, in which we follow n p j only. This is because we wish to focus on the novelty of the method which is related to the fluid aspect of the solutions. In fact, the method for the evolution of f (p) in phase space was extensively studied in [27].

Fluxes in Momentum Space
Time integration of Eq. (14) due to the fluxes in momentum space is done following the method proposed in [27]. We can also retain the diffusive term although here we limit the discussion to the case of a small diffusion coefficient, D p , which can be treated explicitly, i.e.
Here τ Dp is the characteristic diffusion time in momentum space, ∆p = p j+ 1 , the momentum bin size, and ∆t the time step. Then, the terms on the right-hand-side of Eq. (14) produce a change in n p j such that where the superscript n + 1 2 indicates time centering and the last term on the right hand side will be specified in Sections 4.1 and 4.2. We can start estimating the time-averaged flux in momentum space at time t = n∆t as [27] F n which is obtained by time integrating Eq. (16) and by changing integration variable from time to momentum [27]. In Eq. (21) f n ju is the upstream distribution function at time t defined at the grid point in momentum space and p u is the upstream momentum, solution of the integral equation Note that the denominator of the above integrand function has typically a polynomial form so that the integral can be computed in closed form [27].
Finally, a time centered estimate of the term,∇F , is obtained by taking the average of the fluxes, F p , computed at t and t + ∆t, as usually done for nonstiff source terms. Time centering is relevant because the function in Eq. (23) in general depends on the local properties of the fluid, such as density and temperature, which change with time. Note that there are two ways to compute the divergence of the velocity in theṗ term: with a cell centered scheme (∇u) i = (u i+1 − u i−1 )/2∆x or with a staggered scheme, is computed as part of the solution to the Riemann problem that one has to solve in order to advance the fluid equations and estimate the spatial terms appearing in Eq. (20). This is described in the next section. In the following we always use the staggered scheme except in the Godunov's predictor step.

Riemann Problem for Cosmic-ray Hydrodynamics
In this section we describe the modifications to the Riemann problem in the presence of CRs. Without loss of generality we restrict to the one-dimensional case. In addition, for the sake of clarity, in the following discussion we shall neglect the spatial diffusion term, except for the fact that it is implicitly assumed to be at work at shocks, where it causes CR particles to be accelerated. In smooth flows this term is assumed to be slow and is taken into account with an explicit conservative formulation.
We begin with rewriting our system of equations in conservative form: where the entry n p j is repeated for all momentum bins, e.g. for j = 0, N p − 1, and we have introduced the total pressure and specific energy as .
Losses in the total energy of the system are due to escape of energetic particles and, in principle, radiative losses. Thus we writė with b r referring to radiative losses only. Note that the first term is different from zero only ifṗ > 0 because we assume f (p) = 0 for p > p max . In addition, just like the term J j discussed in the previous section, the source termĖ loss will not be treated explicitly but will be part of the shock solution computed through the Riemann solver.
The conservative character of the Eq. (24) for the system in Eq. (25) suggests that the jump conditions across a shock wave can be written in the usual way, provided that the total, thermal plus cosmic-ray, energy and pressure are used. However, there is an additional complication related to the way in which the total pressure and energy is partitioned between CRs and thermal components. This is further addressed in Section 3.2.
For the sake of clarity in the ensuing discussion, we now outline the Riemann problem. Suppose that at t = 0 the gas is described by The solution at t > 0 is in general characterized by two waves: a backward moving wave separating the states (U l , U * l ) and a forward moving wave separating the states (U r , U * r ). Each wave will be either a shock wave or a rarefaction wave. The central states U * l , U * r are separated by a slip line across which the velocity and total pressure (u * , P * ) remain constant, but the density and individual pressure components, (ρ * o , P * o g , P * o c ), o = l, r, in general will change. The value of these three quantities in each intermediate state, U * o , will be reconstructed from U o and the type of wave connecting the two states. In the following two subsections we describe the structure of such waves.

Characteristic Analysis
Transforming into primitive variable space we obtain where we have introduced the CR concentrations, y p j ≡ n p j m p /ρ. These quantities are followed as passive scalars and, for simplicity, will be omitted in the following analysis. Note, however, that the source term ensures that there is consistency between their evolution and that of P c , to which they are related through Eq. (12), (13) and (5).
The system of equations for the primitive variables, V , is obtained with the following transformation where Here c g and c c correspond to the speed of sound associated with the gas and CR pressure respectively. (In principle, a coefficient (∂ ln P c /∂ ln ρ) s should be used instead of γ c in the definition of c g , but the difference is negligible for the purpose of this paper.) Solving for Det(A − λI) = 0, the eigenvalues of A are found to be λ 0 = u − c s , λ 1 = u, λ 2 = u, λ 3 = u + c s . The associated left and right eigenvalues are, respectively, and According to this analysis 'simple waves', including rarefaction waves, are described by the following equations where C # = ρc # , # = g, c, s, is the Lagrangian sound speed or mass flux across the wave. The first two equations (37)-(38) are the usual relation for hydrodynamics with the sound speed modified to account for the CR pressure. The last equation describes the change of CR pressure as a function of the gas pressure during an adiabatic process.

Cosmic-ray Mediated Numerical Shocks
In the following we consider the structure of a shock modified by the presence of accelerated CR particles. As in the case of hydrodynamics, we assume that steady state conditions have been reached. When CR acceleration operates this takes of order the timescale to accelerate thermal particles to relativistic energies. While substantially longer than for a pure hydrodynamic case, in our astrophysical limit this is typically still much shorter than the size of a time step (cf. [21] for further details on this).
If we label quantities far upstream and far downstream of the shock with, − and +, respectively, the Rankine-Hugoniot conditions read where W is the Lagrangian speed of propagation of nonlinear waves. These do not yet specify the amount of CR energy and pressure generated by the dissipation mechanism. In fact, P c and γ c are not and cannot be specified by the above equations alone. To do that we use the following facts without proof (but see [10,3,2]). In the presence of CRs the shock discontinuity is replaced by a precursor, where the gas is compressed adiabatically by the upstreaming CRs and P g ∝ ρ γ ∝ u −γ ; the precursor is immediately followed by a viscous subshock where entropy is generated but both the CR pressure and CR energy flux remain continuous. This means that the structure of the subshock is purely hydrodynamical. If we use the label, 0, to indicate the quantities just prior to the subshock, the compression at the precursor is given by Define the shock Mach number as The gas pressure jump across the total shock transition, being the result of the precursor adiabatic compression and subshock compression 2 , reads [2] An analogous relation can be obtained for the conservation of mass equation, namely Using these results with Euler equation gives [2] P + Recalling that, M = W/C, Eq. (45) and (47) then lead to the following modified definition of the Lagrangian speed of nonlinear waves where C − g = (γ g ρ − P − g ) 1 2 , P = P g + P c . Finally, the value of the CR adiabatic index downstream of the shock, γ + c , is determined by the energy equation through the following relations: whereĖ loss is given in Eq. (28). The term Q loss ≤ 0 describes the energy losses occurring at the precursor and shock front. It becomes important when the CR acceleration efficiency is very high. In this case, this term must also be specified consistently with the acceleration efficiency by the kinetic solution.

Riemann Solver Procedure
Having specified the form of the rarefaction and compression waves modified by the CRs, we can now define the procedure for solving the Riemann problem. First note that, provided the shock acceleration efficiency, η(U − , M), as a function of the upstream conditions (U − ) and the shock Mach number (M), we solve Eq. (47) to derive a similar function for the compression at the precursor, Given the left and right states in Eq. (29), we then want to compute the intersection point, (u * , P * ), of the two wave curves passing through U l , U r in the P −u plane. For this we use the iterative technique proposed in [37], with the two shock approximations [8] and additional modifications which we describe next.
In the absence of shocks we change the Lagrangian speed of nonlinear waves as follows in order to account for the CR pressure W = C l,r s 1 + γ g + 1 2γ g P * − P l,r P l,r with C l,r s = ρ l,r c l,r s . If a shock is present we instead use W given in Eq. (48) with P + replaced by P * . However, unlike the pure hydrodynamic case, W now also depends on r p (M). Thus, using M = W/C g , we have the implicit equation which also needs to be solved iteratively. In addition the tangent slopes to the wave curves in the P −u plane, which are used in the iterative procedure to find P * , are modified according to To summarize, the iteration procedure is now given by : ν = 0, P * 0 = [C l s P r + C r s P l − C l s C r s (u r − u l )]/(C l s + C r s ) while not converged do ν++ if P * ν−1 > P l,r then W l,r ν = W l,r ν (P * ν−1 , W l,r ν /C l,r g ) M l,r = W l,r ν /C l,r g r l,r p = r p (U l,r , M l,r ), Y l,r = Y (r l,r p ), C l,r = C l,r The criterion for convergence can be set to be, |P * ν − P * ν−1 |/P * ν < ǫ, where ǫ is a parameter that sets the error tolerance. Note that at the end of the above procedure, in addition to P * and u * , we have also solved for W l,r and, therefore, M l,r and r l,r p .
Once the left and right moving waves have been determined, we proceed as follows [8]. In searching for the solution at a given point, ξ ≡ x/t, we set σ = sign(ξ − u * ) and define We then complete the definition of the intermediate state U * o . If the latter is separated from U o through a rarefaction wave, knowing P * , P o g and P o c we can use Eq. (39) to estimate P * o g , P * o c . This amounts to solving for P * o g the nonlinear equation and then setting In Eq. (60), γ o c is the CR adiabatic index of the U o state which remains unchanged during an adiabatic process. The density is then estimated through the polytropic law We then evaluate c * s through Eq. (34) and define the wave speedŝ If ξ lies ahead or behind the o-wave we can set the solution to ρ, u, P g , P c , However, ifλ * < ξ <λ o , we have to evaluate the solution inside a rarefaction wave. This requires integration of the system (37)-(39), which cannot be done in closed form and can be expensive. An alternative method is to linearly interpolate between the states U o and U * o as This works just fine for Godunov's method. However, we find that for Glimm's method, when strong rarefactions are involved, it is important to use the exact approach in order to avoid spurious effects.
Before concluding the description of the Riemann solver we point out the presence of a slight inconsistency in the formulation. In fact lim rp→1 W = C s , i.e. the speed of weakly nonlinear waves does not tend to the speed of sound waves. This is a consequence of the large gap in physical scales between the sound waves driven by the total CR and thermal pressure and the CR mediated shock waves. The former are in fact long wavelength perturbations on which scale the diffusion is slow and unimportant. Such scales are much larger than those characterizing the structure of a shock. Thus the conflict in trying to reconcile the two solutions is due to the impossibility of following the intermediate scales.
Since the nonlinear effects due to the process of CR shock acceleration are only important for strong shocks, a natural way of solving the above conflict is to assume that a shock solution is adopted only if the shock propagation speed exceeds that of the sound speed given in Eq. (54). This leads to the following condition

Implementations
The method described in the previous section for the solution of the Riemann problem can be employed to construct time dependent solutions to the system of equations (1)-(4) on a grid.

Glimm's Method
The first implementation we describe is based on Glimm's method. Following [8], here we simply outline the main procedure that we use and refer the reader to the original references [15,6,7,8] for a detailed description.
Consider a piecewise constant approximate solution at time t n = n∆t where ∆x is the mesh size, ∆t the timestep and D the computational domain. We seek to advance the solution by one timestep to t = (n + 1)∆t. To do that we solve the Riemann problem at each cell interface, i − 1 2 , with left and right states given by U n i−1 and U n i . Denote the solution with where we have made explicit use of its property of self-similarity. If the choice of the timestep is sufficiently small, say the wave solutions of the Riemann problems at each cell interface will not interact with each other. Then the set of Riemann solutions, {R i− 1 2 ,n , i ∈ D}, each covering a region, (i − 1)∆x ≤ x < i∆x, defines an exact solution, U e,n (x, t), to the initial value problem in (70) for the time interval, t n < t ≤ t n+1 . The solution at each grid point, i, and time, t n+1 , is obtained by random sampling U e,n as follows: evaluate the solution at the point, x = (i − 1 + a n+1 )∆x, within the region covered by R i− 1 2 ,n , where, a n+1 , is a randomly chosen number, a n+1 ∈ [0, 1). Note that x ∈ i if a n+1 ≥ 1 2 and x ∈ (i − 1) if a n+1 < 1 2 . We can then define (73) Following [8], we use a sampling procedure that is based on van der Corput's pseudo-random sequence, so that a n is the n th element of that sequence.

Hybrid Glimm-Godunov's Method
Shock waves are the only features that need to be propagated without numerical smearing. Therefore, we have also implemented a hybrid scheme which uses Glimm's method to advance shock fronts and Godunov's method for smooth parts of the flow. Having described Glimm's method in the previous section, here we briefly outline a scheme based on the higher order Godunov's method.
In Godunov's method the solution is updated with a conservative scheme where the source term has been described in Section 2.
The fluxes at the cell faces are given by where V is obtained by solving the Riemann problem (discussed extensively in Section 3) with left and right states (V i,+ , V i+1,− ). These states correspond to up-wind time averages, which allow to achieve second order accuracy. They are reconstructed from the cell center, taking into account the effects of spatial gradients and the source term, as follows. At each grid point, i, we compute centered and one-sided slopes and use van Leer's limiter to make the final choice about the local slope, ∆V i . Then, the up-wind, time averaged left (−) and right (+) states at cell faces are Here A is given in Eq. (33), I is the identity operator and n indicates the time-step corresponding to time t. In addition projects out from the state V the components carried by characteristics that propagate away from the cell interface (l j , r j are the left and right eigenstates respectively and λ j is the corresponding eigenvalue, described in Section 3.1).
In conclusion, our hybrid scheme can be summarized as We say that the cell, i, is shocked if a shock is going to cross it during the next timestep. In order for this to happen at least a shock moving with speed u s with respect to the grid must be present at the interface, i − sign(u s ) 1 2 . The criterion for deciding whether or not a wave across an interface qualifies as a shock will be based on the strength of the pressure jump across it, |P i+1 − P i |/ min(P i+1 , P i ), and shall take the condition (68) into account.

Tests
We now present a few tests illustrating the performance of the methods described in the previous sections. The tests consist of a set of Riemann problems with initial conditions (ρ, u, P g , P c , γ c , f (p))[x, t = 0] =      (ρ l , u l , P l g , P l c , γ l c , f l (p)) if x ≤ 0.5, (ρ r , u r , P r g , P r c , γ r c , f r (p)) if x > 0.5, for which we compare the numerical and the 'exact' solutions. The 'exact' solution is obtained by solving the Riemann problem as outlined in Section 3, numerically but without any of the approximations involved in the Riemann solvers for the numerical methods. In particular, no two shock approximation is made, and the exact expression for the speed of rarefaction waves is used.
In order to allow for an easier comparison with the 'exact' solution we only retain the adiabatic terms in the diffusion-convection equation, i.e. we neglect the terms D p in Eq. (16) and b ℓ (p) in Eq. (18). Similarly we let momentum space range over 15 orders of magnitude with p min = 10 −5 and p max = 10 10 . This choice, while unrealistic, is made in order to minimize energy losses due to fluxes across boundaries in momentum space when studying rarefaction waves.
For simplicity the initial conditions for the CR distribution function are specified as an unbroken power-law, with p min ≤ p ≤ p max . CR particles return to the thermal pool for p ≤ p min and escape the system for p ≥ p max . Due to the lack of energy losses the evolution of f (p), followed with the scheme presented in Section 2, becomes trivial. However, the accuracy of that scheme has already been extensively tested in [27], so that here we focus solely on the quality of the hydrodynamic solutions.
In solving the Riemann test problems with the numerical methods presented in this paper we always employ a grid of 128 mesh points on a domain of size unity (so that ∆x = 1/128). We use γ g = 5/3. In addition we use N p = 16 momentum bins. Throughout the Section time is expressed in adimensional code units.
As already mentioned, when evaluating the solution inside a rarefaction wave with the Riemann solver, we have a choice of either integrating directly the Eq. (37)-(39) or use the approximate Eq. (66). When using Glimm's method we test both approaches and compare the results, whereas when using Godunov's method we only employ the approximate approach which turns out sufficiently accurate.
For shock waves, we adopt the following simple prescription defining the shock acceleration efficiency in which the fraction of total momentum impinging on the shock and dissipated into CR pressure depends solely on the shock Mach number. In the above expression, M min and M s are a threshold and scale parameter, respectively. While clearly a simplification, the functional form of η, with a sharp rise for M > M min , followed by a flattening for M > M s , is partially inspired by thermal leakage models and the numerical results described in [20]. We take A = 0.8, M s = 5.77 and M min = 1.5. In this simplified model we use values of Q loss ≤ 0 that, based on the prescribed acceleration efficiency and the energy equation, allow for γ c ≥ 4/3. With the above choices, the resulting shock solutions always admit a subshock, i.e. completely smooth shock transitions do not appear.
The accelerated CR distribution function is assumed to be an unbroken powerlaw as in Eq. (81). Using Eq. (5)- (7), we can thus write the following relation between the slope and the CR adiabatic index in which α(q) ≥ 0 is a function of q through P c (q) and f (p). The above relations imply that q takes a value in the interval (3, ∞) as γ c ranges between 4/3 and 5/3.
It should be pointed out that in a realistic shock the slope of the distribution function, q, in general depends on p and is determined self-consistently with the velocity profile in the precursor and subshock. However, the simplification made here about the distribution function, as well as other naive assumptions made earlier in this section, are solely for the sake of simplicity or easy comparison with exact solutions. Nothing prevents the use of more sophisticated and realistic kinetic models for η(M), p min , p max , q(p) with the numerical method presented in this paper.

Shocks
We begin with three shock problems. Two shocks with mild upstream ratio of CR to thermal pressure, one moving leftward with M = 20 and the other moving rightward with M = 5. And a shock with upstream ratio of CR to thermal pressure equal to one, and moving leftward with M = 10. The problems are specified by the following parameters: for the third problem. The initial CR distribution functions are specified by Eq. (81), with slope q l,r determined respectively by γ l,r c through Eq. (83). In addition to the initial left/right states of the Riemann problems we have also specified the parameter Q loss defined in Eq. (52) and its ratio to the upstream energy flux, F − (ρe). So, in the above three examples about 46%, 0% and 22%, respectively, of the energy flux through the shock front is carried away by escaping CR particles with p > p max .
In these tests the role of Godunov's method in the hybrid formulation is trivial. Therefore we only show the results obtained with Glimm's method. The left and right moving shocks described by Eq. (85)-(86) are presented in the left and right panels of Fig. 1, respectively, while Fig. 2 refers to the case described by Eq. (87). All plots correspond to a solution time, t = 0.05. For each plot the four panels show, from top to bottom, gas density, velocity, gas pressure and CR pressure. The numerical solution reproduces the 'exact' solution very well, without oscillations or artifacts, despite the fact that the CR pressure is comparable or significantly higher than the thermal pressure. Note that the front of the right moving shock is displaced with respect to the 'exact' solution by one cell. This is a characteristic of Glimm's method: as it advances the shock front in discrete steps of size ∆x, it may inevitably place the shock on a grid position that is offset with respect to the 'true' position. By using van der Corput sequence, however, the offset is at most one zone and it eventually becomes negligible when compared to the distance traveled by the shock [8].

Rarefactions
We now turn to the following initial value problem:  representing a rarefaction wave in the λ + characteristic family. The initial CR distribution functions are power-laws, Eq. (81), with the same slope q determined by γ l c = γ r c through Eq. (83). Note that in this case Q loss = 0. Fig. 3 shows the solution at time t = 0.15, obtained with Glimm's method using the approximate (left) and exact (right) rarefaction solution in the Riemann solver, respectively. Again, from top to bottom, each panel shows gas density, velocity, gas pressure and CR pressure.
Glimm's solutions appear slightly ragged. The raggedness in the right panel is purely due to the sampling procedure and does not correspond to an oscillatory behavior. In particular, despite their appearance, the points on the rarefaction curve remain connected through the correct wave solution and the ragged character may or may not be there depending on the specific sample that is being drawn. We show an example of this in the next test, where Glimm's solution of a shock tube problem is characterized by a smooth rarefaction wave.
On the other hand, the solution on the left panel shows additional irregularity which is due to the approximations in Riemann solver. In particular the approximations involved in evaluating the solution inside a rarefaction wave add spurious structure to the Riemann solution which is picked up by Glimm's method. This is undesirable because some of the spurious features may be amplified and even become unstable.
Finally, the left panel of Fig. 4 shows the solution obtained with Godunov's method for the same rarefaction wave problem. Overall the numerical solution reproduces accurately the evolution of the rarefaction wave. Compared with Glimm's method, the solution is now smooth across the wave, although it appears slightly less sharp at the head of the wave.

Shock Tubes
We conclude the set of tests with a shock tube problem with the following initial conditions ρ l = 1.0, u l = 0.0, P l g = 1.0, P l c = 0.0, (γ l c ), ρ r = 9.0, u r = 0.0, P r g = 10, P r c = 6.0, γ r c = 1.33433.
The initial CR distribution function for the right state is a power-law, Eq. (81), with slope q r determined by γ r c through Eq. (83). Note that given the null value of the CR pressure on the left state, we need not specify γ l c nor f l (p). In the shock tube problem in general a shock, a rarefaction and a contact discontinuity develop. In this case the shock is weak and Q loss = 0.
The numerical results at t = 0.15 are shown in the right panel of Fig. 4 for the hybrid Glimm-Godunov's method and in Fig. 5 for Glimm's method. As in the previous test, the latter figure shows both the results obtained by employing an approximate (left) and exact (right) rarefaction solution in the Riemann solver. Note that from top to bottom, each panel now shows gas density, velocity, total (gas+CR) pressure and CR pressure, respectively.
The left panel of Fig. 5 shows additional evidence that when using the approximate rarefaction solution with Glimm's method, spurious structure appears in the numerical solution. This affects not only the rarefaction wave but also the intermediate state. When the exact rarefaction solution is employed however, Glimm's method produces a highly accurate result. Unlike the earlier test, the rarefaction curve is now smooth, confirming how the raggedness in the previous test was due to the sampling procedure. Notice how Glimm's method preserves the sharpness of both the shock and contact discontinuity. As already pointed out, however, their position may be displaced with respect to the exact solution by at most one grid cell.
The right panel of Fig. 4 shows that the numerical solution produced by the hybrid method is also highly accurate. As a sanity check we notice that in switching between Glimm and Godunov's formulations no spurious effect is introduced. In addition, the shock position is the same as in Glimm's method solution. The rarefaction part of the solution is well captured, although now the foot of the rarefaction is not as sharp as in Glimm's method. Finally, the contact discontinuity spreads over a few cells, which is characteristic of Godunov's method. Note that the contact discontinuity appears both in the density and in the pressure components, P g and P c . This, however, does not affect the total pressure P g + P c , which remains perfectly constant in the intermediate state between the rarefaction and the shock, guaranteeing a correct velocity profile as well.

Conclusions
In this paper we have developed a method to include the dynamical effects of CR pressure on a hydrodynamical system.
In smooth flows this is achieved by modification of the characteristic equations which define the waves that carry information in the CR-hydro fluid. The exchange of energy between the two CR particles and the fluid component as a result of diffusive processes both in configuration and momentum space, is modeled with a flux conserving method.
Regarding the solution at shock waves, we have shown that once the acceleration efficiency has been specified as a function of the upstream conditions and shock Mach number, the shock CR mediation and the induced substructure can be correctly taken into account in the fluid solution by modifying the procedure for the Riemann solver.
We have implemented two numerical schemes for obtaining time dependent solutions on a computational grid. One based on Glimm's method and another based on a hybrid Glimm-Godunov method. In both approaches we exploit the ability of Glimm's method to preserve the discontinuous character of shocks. This is useful because when combined with the aforementioned modified Riemann solver it provides a natural scheme for advancing the shock solution at the correct speed, meaning that the CR dynamical effects are taken into account without resolving the shock substructure.
In smooth flows Glimm's method is not the only possible choice and Godunov's method can also be employed. Our tests show that Glimm's method is rather sensitive to the approximations assumed in the Riemann solver procedure. In particular, when evaluating the solution inside a rarefaction wave it is important to solve the exact equations (37)-(39) or else spurious features will appear. This is not the case with Godunov's method. Compared to the version of Glimm's method with the exact rarefaction solution in the Riemann solver, Godunov's produces smoother solutions and is only slightly less sharp in capturing the head of rarefaction waves. So the hybrid method is also a viable option.
The proposed method can be readily employed for the study of one-dimensional models of astrophysical systems, such as simplified radially symmetric descriptions of, e.g., Supernova Remnants or Galaxy Clusters. The potential benefit of this scheme, however, lies in the possibility of coupling it with three-dimensional shock tracking algorithms or extending Glimm's method to three dimensions. In this case the dynamical role of CRs in astrophysical systems can be studied without geometrical restrictions. Since numerical dissipation is necessary for stable hydrodynamic shocks in three dimensions [8], but spoils the correct CR-hydrodynamic shock solution, shock tracking algorithms may offer the only way to self-consistently study CR-hydrodynamics in multi-dimensions.