text:
# EXACT EXPRESSIONS FOR THE CRITICAL MACH NUMBERS IN THE TWO-FLUID MODEL OF COSMIC-RAY MODIFIED SHOCKS

## 1 INTRODUCTION

It is now widely accepted that acceleration in supernova-driven shock waves plays an important role in the production of the observed cosmic ray spectrum up to energies of $`10^{15}`$eV (Heavens 1984a; Ko 1995a), and it is plausible that acceleration in shocks near stellar-size compact objects can produce most of the cosmic radiation observed at higher energies (Jones & Ellison 1991). In the generic shock acceleration model, cosmic rays scatter elastically with magnetic irregularities (MHD waves) that are frozen into the background (thermal) gas (Gleeson & Axford 1967; Skilling 1975). In crossing the shock, these waves experience the same compression and deceleration as the background gas, if the speed of the waves with respect to the gas (roughly the Alfvén speed) is negligible compared with the flow velocity (Achterberg 1987). The convergence of the scattering centers in the shock creates a situation where the cosmic rays gain energy systematically each time they cross the shock. Since the cosmic rays are able to diffuse spatially, they can cross the shock many times. In this process, an exponentially small fraction of the cosmic rays experience an exponentially large increase in their momentum due to repeated shock crossings. The characteristic spectrum resulting from this first-order Fermi process is a power-law in momentum (Krymskii 1977; Bell 1978a, b; Blandford & Ostriker 1978).

It was recognized early on that if the cosmic rays in the downstream region carry away a significant fraction of the momentum flux supplied by the incident (upstream) gas, then the dynamical effect of the cosmic-ray pressure must be included in order to obtain an accurate description of the shock structure (Axford, Leer, & Skadron 1977). In this scenario, the coupled nonlinear problem of the gas dynamics and the energization of the cosmic rays must be treated in a self-consistent manner. A great deal of attention has been focused on the “two-fluid” model for diffusive shock acceleration as a possible description for the self-consistent cosmic-ray modified shock problem. In this steady-state theory, first analyzed in detail by Drury & Völk (1981, hereafter DV), the cosmic rays and the background gas are modeled as interacting fluids with constant specific heat ratios $`\gamma _c`$ and $`\gamma _g`$, respectively. The coupling between the cosmic rays and the gas is provided by MHD waves, which serve as scattering centers but are otherwise ignored. The cosmic rays are treated as massless particles, and second-order Fermi acceleration due to the stochastic wave propagation is ignored. Within the context of the two-fluid model, DV were able to demonstrate the existence of multiple (up to 3) distinct dynamical solutions for certain upstream boundary conditions. The solutions include flows that are smooth everywhere as well as flows that contain discontinuous, gas-mediated “subshocks.” Subsequently, multiple solutions have also been obtained in modified two-fluid models that incorporate a source term representing the injection of seed cosmic rays (Ko, Chan, & Webb 1997; Zank, Webb, & Donohue 1993). Only one solution can be realized in a given flow, but without incorporating additional physics one cannot determine which solution it will be. The two-fluid model has been extended to incorporate a quantitative treatment of the MHD wave field by McKenzie & Völk (1982) and Völk, Drury, & McKenzie (1984) using a “three-fluid” approach.

During the intervening decades, a great deal of effort has been expended on analyzing the structure and the stability of cosmic-ray modified shocks (see Jones & Ellison 1991 and Ko 1995b for reviews). Much of this work has focused on the time-dependent behavior of the two-fluid model, which is known to be unstable to the development of acoustic waves (Drury & Falle 1986; Kang, Jones, & Ryu 1992) and magnetosonic waves (Zank, Axford, & McKenzie 1990). Ryu, Kang, & Jones (1993) extended the analysis of acoustic modes to include a secondary, Rayleigh-Taylor instability. In most cases, it is found that the cosmic-ray pressure distribution is not substantially modified by the instabilities.

The two-fluid theory of DV suffers from a “closure problem” in the sense that there is not enough information to compute the adiabatic indices $`\gamma _g`$ and $`\gamma _c`$ self-consistently, and therefore they must be treated as free parameters (Achterberg, Blandford, & Periwal 1984; Duffy, Drury, & Völk 1994). This has motivated the subsequent development of more complex theories that utilize a diffusion-convection transport equation to solve for the cosmic-ray momentum distribution along with the flow structure self-consistently. In these models, “seed” cosmic rays are either advected into the shock region from far upstream, or injected into the gas within the shock itself. Interestingly, Kang & Jones (1990), Achterberg, Blandford, & Periwal (1984), and Malkov (1997a; 1997b) found that diffusion-convection theories can also yield multiple dynamical solutions for certain values of the upstream parameters, in general agreement with the two-fluid model. Frank, Jones, & Ryu (1995) have obtained numerical solutions to the time-dependent diffusion-convection equation for oblique cosmic-ray modified shocks that are in agreement with the predictions of the steady-state two-fluid model. These studies suggest that, despite its shortcomings, the two-fluid theory remains one of the most useful tools available for analyzing the acceleration of cosmic rays in shocks waves (Ko 1995b).

In their approach to modeling the diffusive acceleration of cosmic rays, DV stated the upstream boundary conditions for the incident flow in terms of the total Mach number,

$$M\frac{v}{a},a\sqrt{\frac{\gamma _gP_g}{\rho }+\frac{\gamma _cP_c}{\rho }},$$
$`(1.1)`$
and the ratio of the cosmic-ray pressure to the total pressure,

$$Q\frac{P_c}{P},PP_c+P_g,$$
$`(1.2)`$
where $`v`$, $`a`$, $`\rho `$, $`P`$, $`P_c`$, and $`P_g`$ denote the flow velocity of the background gas, the total sound speed, the gas density, and the total, cosmic-ray, and gas pressures, respectively. We will use the subscripts “0” and “1” to denote quantities associated with the far upstream and downstream regions, respectively. DV described the incident flow conditions by selecting values for $`M_0`$ and $`Q_0`$. Once these parameters have been specified, the determination of the flow structure (and in particular the number of possible solutions) in the simplest form of the two-fluid model requires several stages of root finding. Since the analysis is inherently numerical in nature, the results are usually stated only for specific upstream conditions.

The characterization of the upstream conditions in terms of $`M_0`$ and $`Q_0`$ employed by DV turns out to be an inconvenient choice from the point of view of finding exact critical relations describing the number of possible flow solutions for given upstream conditions. As an alternative approach, it is possible to work in terms of the individual gas and cosmic-ray Mach numbers, defined respectively by

$$M_g\frac{v}{a_g},M_c\frac{v}{a_c},$$
$`(1.3)`$
where

$$a_g\sqrt{\frac{\gamma _gP_g}{\rho }},a_c\sqrt{\frac{\gamma _cP_c}{\rho }},$$
$`(1.4)`$
denote the gas and cosmic-ray sound speeds, respectively. According to equation (1.1), $`a^2=a_g^2+a_c^2`$, and therefore the Mach numbers $`M_g`$ and $`M_c`$ are related to the total Mach number $`M`$ and the pressure ratio $`Q`$ via

$$M^2=M_g^2+M_c^2,Q=\frac{\gamma _gM_g^2}{\gamma _cM_c^2+\gamma _gM_g^2},$$
$`(1.5)`$
or, equivalently,

$$M_g=\left(1+\frac{\gamma _c}{\gamma _g}\frac{Q}{1Q}\right)^{1/2}M,M_c=\left(1+\frac{\gamma _g}{\gamma _c}\frac{1Q}{Q}\right)^{1/2}M.$$
$`(1.6)`$
Since these equations apply everywhere in the flow, the boundary conditions in the two-fluid model can evidently be expressed by selecting values for any two of the four upstream parameters $`(M_0,Q_0,M_{g0},M_{c0})`$. In their work, DV described the upstream conditions using $`(M_0,Q_0)`$, whereas Ko, Chan, & Webb (1997) and Axford, Leer, & McKenzie (1982) used $`(M_{g0},Q_0)`$. Another alternative, which apparently has not been considered before, is to use the parameters $`(M_{g0},M_{c0})`$. Although these choices are all equivalent from a physical point of view, we demonstrate below that the set $`(M_{g0},M_{c0})`$ is the most advantageous mathematically because it allows us to derive exact constraint curves that clearly delineate the regions of various possible behavior in the parameter space of the two-fluid model. This approach exploits the formal symmetry between the cosmic-ray quantities and the gas quantities as they appear in the expressions describing the asymptotic states of the flow.

The remainder of the paper is organized as follows. In § 2 we discuss the transport equation for the cosmic rays and derive the associated moment equation for the variation of the cosmic-ray energy density. In § 3 we employ momentum and energy conservation to obtain an exact result for the critical upstream cosmic-ray Mach number that determines whether smooth flow is possible for a given value of the upstream gas Mach number. In § 4 we derive exact critical conditions for the existence of multiple solutions containing a discontinuous, gas-mediated subshock. The resulting curves are plotted and analyzed for various values of the adiabatic indices $`\gamma _g`$ and $`\gamma _c`$. In § 5 we present specific examples of multiple-solution flows that verify the predictions made using our analytical critical conditions. We conclude in § 6 with a general discussion of our results and their significance for the theory of diffusive cosmic-ray acceleration.

## 2 GOVERNING EQUATIONS

The two-fluid model is developed by treating the cosmic rays as a fluid with energy density comparable to that of the background gas, but possessing negligible mass. In this section we review the basic equations relevant for the two-fluid model. For integrity and clarity of presentation, we also include re-derivations of a few of the published results concerning the overall shock structure and the nature of the transonic flow.

### 2.1 Lagrangian Equations

The diffusive acceleration of energetic cosmic rays due to the convergence of scattering centers in a one-dimensional, plane-parallel flow is described by the transport equation (Skilling 1971; 1975)

$$\frac{Df}{Dt}=\frac{p}{3}\frac{f}{p}\frac{v}{x}+\frac{}{x}\left(\kappa \frac{f}{x}\right),$$
$`(2.1)`$
where $`p`$ is the particle momentum, $`v(x,t)`$ is the flow velocity of the background gas (taken to be positive in the direction of increasing $`x`$), $`\kappa (x,p,t)`$ is the spatial diffusion coefficient, and the operator

$$\frac{D}{Dt}\frac{}{t}+v\frac{}{x}$$
$`(2.2)`$
expresses the comoving (Lagrangian) time derivative in the frame of the gas. Equation (2.1) describes the effects of Fermi acceleration, bulk advection, and spatial diffusion on the direction-integrated (isotropic) cosmic-ray momentum distribution $`f(x,p,t)`$, which is normalized so that the total number density of the cosmic rays is given by

$$n_c(x,t)=_0^{\mathrm{}}4\pi p^2f𝑑p.$$
$`(2.3)`$
Note that equation (2.1) neglects the second-order Fermi acceleration of the cosmic rays that occurs due to stochastic wave propagation, which is valid provided the Alfvén speed $`v_A=B/\sqrt{4\pi \rho }`$ is much less than the flow velocity $`v`$, where $`B`$ is the magnetic field strength. Furthermore, equation (2.1) does not include a particle collision term, and therefore it is not applicable to the background gas, which is assumed to have a thermal distribution. The momentum, mass, and energy conservation equations for the gas can be written in the comoving frame as

$$\frac{Dv}{Dt}=\frac{1}{\rho }\frac{P}{x},\frac{D\rho }{Dt}=\rho \frac{v}{x},\frac{DU_g}{Dt}=\gamma _g\frac{U_g}{\rho }\frac{D\rho }{Dt},$$
$`(2.4)`$
respectively, where $`U_g=P_g/(\gamma _g1)`$ is the internal energy density of the gas. The expression for $`DU_g/Dt`$ in equation (2.4) implies a purely adiabatic variation of $`U_g`$, and therefore it neglects any heating or cooling of the gas due to wave generation or damping. This adiabatic equation must be replaced with the appropriate Rankine-Hugoniot jump conditions at a discontinuous, gas-mediated subshock, should one occur in the flow. In the case of a relativistic subshock, the momentum conservation equation for the gas must be modified to reflect the anisotropy of the pressure distribution (e.g., Kirk & Webb 1988).

### 2.2 Cosmic-Ray Energy Equation

The pressure $`P_c`$ and the energy density $`U_c`$ associated with the isotropic cosmic-ray momentum distribution $`f`$ are given by (Duffy, Drury, & Völk 1994)

$$P_c(x,t)=_0^{\mathrm{}}\frac{4\pi }{3}p^3Vf𝑑p,U_c(x,t)=_0^{\mathrm{}}4\pi p^2\epsilon f𝑑p,$$
$`(2.5)`$
where

$$\epsilon =(\gamma 1)mc^2,V=\frac{p}{\gamma m},\gamma =\sqrt{\frac{p^2}{m^2c^2}+1},$$
$`(2.6)`$
denote respectively the kinetic energy, the speed, and the Lorentz factor of a cosmic ray with momentum $`p`$ and mass $`m`$. Although the lower bound of integration is formally taken to be $`p=0`$, in practice the cosmic rays are highly relativistic particles, and therefore $`f`$ vanishes for $`p\stackrel{<}{}mc`$. If the distribution has the power-law form $`fp^q`$, then we must have $`4<q<5`$ in order to avoid divergence in the integrals for $`P_c`$ and $`U_c`$ (Achterberg, Blandford, & Periwal 1984), although this restriction can be lifted if cutoffs are imposed at high and/or low momentum (Kang & Jones 1990).

We can obtain a conservation equation for the cosmic-ray energy density $`U_c`$ by operating on the transport equation (2.1) with $`_0^{\mathrm{}}4\pi p^2T𝑑p`$, yielding

$$\frac{DU_c}{Dt}=\gamma _cU_c\frac{v}{x}+\frac{}{x}\left(\overline{\kappa }\frac{U_c}{x}\right),$$
$`(2.7)`$
where the mean diffusion coefficient $`\overline{\kappa }`$ is defined by (Duffy, Drury, & Völk 1994)

$$\overline{\kappa }(x,t)\frac{_0^{\mathrm{}}p^2T\kappa (f/x)𝑑p}{_0^{\mathrm{}}p^2T(f/x)𝑑p},$$
$`(2.8)`$
and the cosmic-ray adiabatic index $`\gamma _c`$ is defined by (Malkov & Völk 1996)

$$\gamma _c(x,t)\frac{4}{3}+\frac{1}{3}\frac{_0^{\mathrm{}}p^2Tf\gamma ^1𝑑p}{_0^{\mathrm{}}p^2Tf𝑑p}=1+\frac{P_c}{U_c}.$$
$`(2.9)`$
Note that in deriving equation (2.7), we have dropped an extra term that arises via integration by parts because it must vanish in order to obtain finite values for $`P_c`$ and $`U_c`$. The integral expression in equation (2.9) indicates that $`\gamma _c`$ must lie in the range $`4/3\gamma _c5/3`$. It also demonstrates that $`\gamma _c`$ will evolve in response to changes in the shape of the momentum distribution $`f`$. The closure problem in the two-fluid model arises because $`f`$ is not calculated at all, and therefore $`\gamma _c`$ must be imposed rather than computed self-consistently.

### 2.3 Eulerian Equations

The conservation equations can be rewritten in standard Eulerian form as

$$\frac{\rho }{t}=\frac{J}{x},$$
$`(2.10)`$
$$\frac{}{t}\left(\rho v\right)=\frac{I}{x},$$
$`(2.11)`$
$$\frac{}{t}\left(\frac{1}{2}\rho v^2+U_g+U_c\right)=\frac{E}{x},$$
$`(2.12)`$
where the fluxes of mass, momentum, and total energy are given respectively by

$$J\rho v,$$
$`(2.13)`$
$$I\rho v^2+P_g+P_c,$$
$`(2.14)`$
$$E\frac{1}{2}\rho v^3+v\left(P_g+U_g\right)+v\left(P_c+U_c\right)\overline{\kappa }\frac{U_c}{x}.$$
$`(2.15)`$
The momentum and energy fluxes can be expressed in dimensionless form as

$$\frac{I}{Jv_0}=u+𝒫_g+𝒫_c,$$
$`(2.16)`$
$$\frac{E}{Jv_0^2}=\frac{1}{2}u^2+\frac{\gamma _g}{\gamma _g1}u𝒫_g+\frac{\gamma _c}{\gamma _c1}u𝒫_c\frac{\overline{\kappa }}{Jv_0^2}\frac{}{x}\left(\frac{Jv_0𝒫_c}{\gamma _c1}\right),$$
$`(2.17)`$
where $`v_0`$ is the asymptotic upstream flow velocity and the dimensionless quantities $`u`$, $`𝒫_g`$, and $`𝒫_c`$ are defined respectively by

$$u\frac{v}{v_0},𝒫_g\frac{P_g}{Jv_0},𝒫_c\frac{P_c}{Jv_0}.$$
$`(2.18)`$
Note that the definition of $`u`$ implies that the incident flow has $`u_0=1`$. These relations can be used to rewrite equations (1.3) for the gas and cosmic-ray Mach numbers as

$$M_g^2=\frac{u}{\gamma _g𝒫_g},M_c^2=\frac{u}{\gamma _c𝒫_c}.$$
$`(2.19)`$
### 2.4 The Dynamical Equation

In this paper we shall adopt the two-fluid approximation in the form used by DV, and therefore we assume that the adiabatic indices $`\gamma _g`$ and $`\gamma _c`$ are each constant throughout the flow. The assumption of constant $`\gamma _g`$ is probably reasonable since the background gas is expected to remain thermal and nonrelativistic at all locations. The assumption of constant $`\gamma _c`$ is more problematic, since we expect the cosmic ray distribution to evolve throughout the flow in response to Fermi acceleration, but it is justifiable if the “seed” cosmic rays are already relativistic in the far upstream region. We also assume that a steady state prevails, so that the fluxes $`J`$, $`I`$, and $`E`$ are all conserved. In this case the quantities $`𝒫_g`$ and $`𝒫_c`$ express the pressures of the two species relative to the upstream ram pressure of the gas $`\rho _0v_0^2=Jv_0`$, where $`\rho _0`$ is the asymptotic upstream mass density. The Eulerian frame in which we are working is necessarily the frame of the shock, since that is the only frame in which the flow can appear stationary (Becker 1998). In a steady state, the adiabatic variation of $`𝒫_g`$ implied by equation (2.4) indicates that along any smooth section of the flow, the gas pressure can be calculated in terms of the velocity using

$$𝒫_g=𝒫_g\left(\frac{u}{u_{}}\right)^{\gamma _g},$$
$`(2.20)`$
where $`𝒫_g`$ and $`u_{}`$ denote fiducial quantities measured at an arbitrary, fixed location within the section of interest. According to equation (2.19), the associated variation of the gas Mach number along the smooth section of the flow is given by

$$M_g^2=M_g^2\left(\frac{u}{u_{}}\right)^{1+\gamma _g},$$
$`(2.21)`$
where $`M_g`$ denotes the gas Mach number at the fiducial location. Substituting for $`𝒫_g`$ in equation (2.16) using equation (2.20) and differentiating the result with respect to $`x`$ yields the dynamical equation (Achterberg 1987; Ko, Chan, & Webb 1997)

$$\frac{du}{dx}=\frac{d𝒫_c/dx}{M_g^21}.$$
$`(2.22)`$
Critical points occur where the numerator and denominator vanish simultaneously. The vanishing of the denominator implies that $`M_g=1`$ at the critical point, and therefore the critical point is also a gas sonic point (Axford, Leer, & McKenzie 1982). The vanishing of the numerator implies that $`d𝒫_c/dx=0`$ at the gas sonic point.

### 2.5 Transonic Flow Structure

We can rewrite the dimensionless momentum flux $``$ by using equations (2.19) to substitute for $`𝒫_g`$ and $`𝒫_c`$ in equation (2.16), yielding

$$=u\left(1+\frac{M_g^2}{\gamma _g}+\frac{M_c^2}{\gamma _c}\right).$$
$`(2.23)`$
Similarly, equation (2.17) for the dimensionless energy flux $``$ becomes

$$=u^2\left(\frac{1}{2}+\frac{M_g^2}{\gamma _g1}+\frac{M_c^2}{\gamma _c1}\right)\frac{1}{\gamma _c1}\frac{\overline{\kappa }}{v_0}\frac{d𝒫_c}{dx}.$$
$`(2.24)`$
Using equation (2.23) to eliminate $`M_c`$ in equation (2.24) yields for the gradient of the cosmic-ray pressure

$$g(u)\frac{1}{\gamma _c1}\frac{\overline{\kappa }}{v_0}\frac{d𝒫_c}{dx}=\left(\frac{1}{2}\mathrm{\Gamma }_c\right)u^2+\mathrm{\Gamma }_cu+\left(\mathrm{\Gamma }_g\mathrm{\Gamma }_c\right)\frac{u^2}{\gamma _gM_g^2},$$
$`(2.25)`$
where

$$\mathrm{\Gamma }_g\frac{\gamma _g}{\gamma _g1},\mathrm{\Gamma }_c\frac{\gamma _c}{\gamma _c1}.$$
$`(2.26)`$
Along any smooth section of the flow, $`g`$ depends only on $`u`$ by virtue of equation (2.21), which gives $`M_g`$ as a function of $`u`$. In the two-fluid model, the flow is assumed to become gradient-free asymptotically, so that $`d𝒫_c/dx0`$ in the far upstream and downstream regions (DV; Ko, Chan, & Webb 1997). The function $`g`$ must therefore vanish as $`|x|\mathrm{}`$, and consequently we can express the values of $``$ and $``$ in terms of the upstream Mach numbers $`M_{g0}`$ and $`M_{c0}`$ using

$$=1+\frac{M_{g0}^2}{\gamma _g}+\frac{M_{c0}^2}{\gamma _c},$$
$`(2.27)`$
and

$$=\frac{1}{2}+\frac{M_{g0}^2}{\gamma _g1}+\frac{M_{c0}^2}{\gamma _c1},$$
$`(2.28)`$
where we have also employed the boundary condition $`u_0=1`$.

The critical nature of the dynamical equation (2.22) implies that $`g=0`$ at the gas sonic point. Hence if the flow is to pass smoothly through a gas sonic point, then $`g`$ must vanish at three locations. It follows that one of the key questions concerning the flow structure is the determination of the number of points at which $`g=0`$. We can address this question by differentiating equation (2.25) with respect to $`u`$, which yields

$$\frac{dg}{du}=\frac{u}{\gamma _c1}\left(M^21\right),$$
$`(2.29)`$
where $`M=(M_g^2+M_c^2)^{1/2}`$ is the total Mach number, and we have used the result

$$\frac{dM_g^2}{du}=\frac{M_g^2}{u}\left(1+\gamma _g\right)$$
$`(2.30)`$
implied by equation (2.21). For the second derivative of $`g`$ we obtain

$$\frac{d^2g}{du^2}=\frac{1\gamma _cM_g^2(\gamma _g\gamma _c)}{\gamma _c1}.$$
$`(2.31)`$
Since the cosmic rays have a higher average Lorentz factor than the thermal background gas, $`\gamma _g>\gamma _c`$ (cf. eq. 2.9), and therefore $`d^2g/du^2<0`$, implying that $`g(u)`$ is concave down as indicated in Figure 1. Hence there are exactly two roots for $`u`$ that yield $`g=0`$, one given by the upstream velocity $`u=u_0=1`$ and the other given by the downstream velocity, denoted by $`u=u_1`$. We therefore conclude that if the flow includes a gas sonic point, then the velocity at that point must be either $`u_0`$ or $`u_1`$. Consequently the flow cannot pass smoothly through a gas sonic point, as first pointed out by DV.

The flows envisioned here are decelerating, and therefore the high-velocity root $`u=u_0`$ is associated with the incident flow. In this case, the fact that $`g(u)`$ is concave-down implies that $`dg/du<0`$ in the upstream region, and therefore based on equation (2.29) we conclude that $`M>1`$ in the upstream region and $`M<1`$ in the downstream region, with $`M=1`$ at the peak of the curve where $`dg/du=0`$. Hence the flow must contain a sonic transition with respect to the total sound speed $`a`$. In this sense, the flow is a “shock” whether or not it contains an actual discontinuity. For the flow to decelerate through a shock transition, the total upstream Mach number must therefore satisfy the condition $`M_0>1`$. This constraint also implies that the upstream flow must be supersonic with respect to both the gas and cosmic-ray sound speeds (i.e., $`M_{g0}>1`$ and $`M_{c0}>1`$), since $`M_0^2=M_{g0}^2+M_{c0}^2`$. Furthermore, for a given value of $`M_{g0}`$, the upstream cosmic-ray Mach number $`M_{c0}`$ must exceed the minimum value

$$M_{c,\mathrm{min}}\left(1M_{g0}^2\right)^{1/2},$$
$`(2.32)`$
corresponding to the limit $`M_0=1`$. The requirement that $`M_{g0}>1`$ forces us to conclude that if a gas sonic point exists in the flow, then it must be identical to the gradient-free asymptotic downstream state. Consequently the flow must either remain supersonic everywhere with respect to the gas, or it must cross a discontinuous, gas-mediated subshock. If $`M_g>1`$ everywhere, then the flow is completely smooth and the gas sonic point is “virtual,” meaning that it exists in the parameter space, but it does not lie along the flow trajectory. In this case, the gas pressure evolves in a purely adiabatic fashion according to equation (2.20), although the total entropy of the combined system (gas plus cosmic rays) must increase as the flow crosses the shock, despite the fact that it is smooth. In § 3 we derive the critical value for the upstream cosmic-ray Mach number $`M_{c0}`$ that determines whether or not smooth flow is possible for a given value of the upstream gas Mach number $`M_{g0}`$.

## 3 CRITICAL CONDITIONS FOR SMOOTH FLOW

The overall structure of a cosmic-ray modified shock governed by the dynamical equation (2.22) can display a variety of qualitatively different behaviors, as first pointed out by DV. Depending on the upstream parameters, up to three distinct steady-state solutions are possible, although only one of these can be realized in a given situation. The possibilities include globally smooth solutions as well as solutions containing a discontinuous, gas-mediated subshock. Smooth flow is expected when the upstream cosmic-ray pressure $`P_{c0}`$ is sufficiently large since in this case cosmic ray diffusion is able to smooth out the discontinuity. In this section we utilize the critical nature of the dynamical equation to derive an analytic expression for the critical condition that determines when smooth flow is possible, as a function of the upstream (incident) Mach numbers $`M_{c0}`$ and $`M_{g0}`$.

### 3.1 Critical Cosmic-Ray Mach Number

Whether or not the flow contains a discontinuous, gas-mediated subshock, it must be smooth in the upstream region (preceding the subshock if one exists). We can therefore apply equation (2.21) for the variation of $`M_g`$ in the upstream region, where it is convenient to use the incident parameters $`u_0=1`$ and $`M_{g0}`$ as the fiducial quantities. The requirement that $`M_g=1`$ at the gas sonic point implies that the velocity there is given by

$$u_sM_{g0}^{2/(1+\gamma _g)},$$
$`(3.1)`$
which we refer to as the “critical velocity.” If a sonic point exists in the flow, then $`u_s`$ must correspond to the downstream root of the equation $`g(u)=0`$, i.e., $`u_s=u_1`$. Note that the value of $`u_s`$ depends only on $`M_{g0}`$, and consequently it is independent of $`M_{c0}`$.

The asymptotic states of the flow are assumed to be gradient-free, and the critical conditions associated with the dynamical equation (2.22) imply that $`d𝒫_c/dx=0`$ at the gas sonic point. The constancy of $``$ and $``$ therefore allows us to link upstream quantities to quantities at the gas sonic point by using equations (2.23), (2.24), (2.27), and (2.28) to write

$$=u_s\left(1+\frac{1}{\gamma _g}+\frac{M_{cs}^2}{\gamma _c}\right)=1+\frac{M_{g0}^2}{\gamma _g}+\frac{M_{c0}^2}{\gamma _c},$$
$`(3.2)`$
and

$$=u_s^2\left(\frac{1}{2}+\frac{1}{\gamma _g1}+\frac{M_{cs}^2}{\gamma _c1}\right)=\frac{1}{2}+\frac{M_{g0}^2}{\gamma _g1}+\frac{M_{c0}^2}{\gamma _c1},$$
$`(3.3)`$
respectively, where $`M_{cs}`$ denotes the value of the cosmic-ray Mach number at the gas sonic point. If we substitute for $`u_s`$ using equation (3.1) and eliminate $`M_{cs}`$ by combining equations (3.2) and (3.3), we can solve for $`M_{c0}`$ to obtain an exact expression for the critical upstream cosmic-ray Mach number required for the existence of a sonic point in the asymptotic downstream region,

$$M_{c0}=M_{cA}\left\{\frac{\left[\left(R_01\frac{1}{\gamma _g}\right)M_{g0}^2+\frac{R_0}{\gamma _g}\right]\gamma _c\left[\left(\frac{R_0^2}{2}\frac{1}{2}\frac{1}{\gamma _g1}\right)M_{g0}^2+\frac{R_0^2}{\gamma _g1}\right]\left(\gamma _c1\right)}{R_0(R_01)M_{g0}^2}\right\}^{1/2},$$
$`(3.4)`$
where

$$R_0M_{g0}^{2/(1+\gamma _g)}.$$
$`(3.5)`$
Note that $`M_{cA}`$ is an explicit function of the upstream gas Mach number $`M_{g0}`$. The interpretation is that if $`M_{c0}=M_{cA}`$ for a given value of $`M_{g0}`$, then the flow is everywhere supersonic with respect to the gas sound speed $`a_g`$ except in the far-downstream region, where it asymptotically approaches the gas sonic point. Surprisingly, this simple solution for $`M_{cA}`$ has apparently never before appeared in the literature, probably due to the fact that the analytical form is lost when one works in terms of the alternative parameters $`M_0`$ and $`Q_0`$ employed by DV. This can be clearly demonstrated by using the expressions

$$M_{g0}=\left(1+\frac{\gamma _c}{\gamma _g}\frac{Q_0}{1Q_0}\right)^{1/2}M_0,M_{c0}=\left(1+\frac{\gamma _g}{\gamma _c}\frac{1Q_0}{Q_0}\right)^{1/2}M_0,$$
$`(3.6)`$
to substitute for $`M_{g0}`$ and $`M_{c0}`$ in equations (3.4) and (3.5) and then attempting to solve the resulting equation for either $`Q_0`$ or $`M_0`$. It is easy to convince oneself that it is not possible to express either of these quantities explicitly in terms of the other. In Figure 2 we depict the curve $`M_{c0}=M_{cA}`$ in the $`(M_{g0},M_{c0})`$ parameter space using equations (3.4) and (3.5) for various values of $`\gamma _g`$ and $`\gamma _c`$.

### 3.2 Smooth Flow Criterion

We have determined that smooth flow into an asymptotic downstream gas sonic point is possible if $`M_{c0}=M_{cA}`$. However, in order to obtain a complete understanding of the significance of the critical upstream cosmic-ray Mach number $`M_{cA}`$, we must determine the nature of the flow when $`M_{c0}M_{cA}`$. The resulting flow structure can be analyzed by perturbing around the state $`M_{c0}=M_{cA}`$ by taking the derivative of the asymptotic downstream velocity $`u_1`$ with respect to $`M_{c0}`$, holding $`M_{g0}`$ constant. The fact that $`M_{g0}`$ is held fixed implies that the critical velocity $`u_s`$ also remains unchanged by virtue of equation (3.1). Upon differentiating, we obtain

$$\left(\frac{u_1}{M_{c0}}\right)_{M_{g0}}=\frac{2}{M_{c0}^3}\left[\frac{\gamma _c+1}{2}+\frac{(\gamma _g\gamma _c)(\gamma _g1+u_1^{\gamma _g}\gamma _gu_1)}{\gamma _g(\gamma _g1)M_{g0}^2(1u_1)^2u_1^{\gamma _g}}\right]^1,$$
$`(3.7)`$
which is always negative since the flow decelerates, and therefore $`u_1<1`$. This indicates that if $`M_{c0}`$ is decreased from the value $`M_{c0}=M_{cA}`$ for fixed $`M_{g0}`$, then the downstream velocity $`u_1`$ increases above the critical velocity $`u_s`$, and therefore the flow is everywhere supersonic with respect to the gas sound speed. In this case there is no gas sonic point in the flow, and consequently a globally smooth solution is possible. Conversely, when $`M_{c0}>M_{cA}`$, a gas sonic point exists in the flow, and therefore the flow cannot be globally smooth because that would require smooth passage through a gas sonic point, which we have proven to be impossible. In this case, the flow must pass through a discontinuous, gas-mediated subshock. We conclude that globally smooth flow is possible in the region below each of the critical curves plotted in Figure 2.

Note that in each case there is a critical value for $`M_{g0}`$, denoted by $`M_{gA}`$, to the right of which smooth flow is always possible for any value of $`M_{c0}`$. This critical value is the solution to the equation

$$\left[\left(R_A1\frac{1}{\gamma _g}\right)M_{gA}^2+\frac{R_A}{\gamma _g}\right]\gamma _c\left[\left(\frac{R_A^2}{2}\frac{1}{2}\frac{1}{\gamma _g1}\right)M_{gA}^2+\frac{R_A^2}{\gamma _g1}\right]\left(\gamma _c1\right)=0,$$
$`(3.8)`$
where

$$R_AM_{gA}^{2/(1+\gamma _g)},$$
$`(3.9)`$
corresponding to the limit $`M_{cA}\mathrm{}`$ in equation (3.4). We plot $`M_{gA}`$ as a function of $`\gamma _g`$ and $`\gamma _c`$ in Figure 3. When $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$, we find that $`M_{gA}=12.28`$, in agreement with the numerical results of DV and Heavens (1984b).

## 4 CRITICAL CONDITIONS FOR MULTIPLE SOLUTIONS

DV discovered that two new discontinuous solutions become available in addition to either a smooth solution or another discontinuous solution when the upstream total Mach number $`M_0`$ is sufficiently large. Subsequent authors have confirmed the existence of multiple dynamical solutions within the context of the two-fluid model (Achterberg, Blandford, & Periwal 1984; Axford, Leer, & McKenzie 1982; Kang & Jones 1990; Ko, Chan, & Webb 1997; Zank, Webb, & Donohue 1993). However, most of this work utilized numerical root-finding procedures and therefore it fails to provide much insight into the structure of the critical conditions that determine when multiple solutions are possible. We revisit the problem in this section by recasting the upstream boundary conditions using the same approach employed in § 3. In particular, we show that when the incident flow conditions are stated in terms of the upstream gas and cosmic-ray Mach numbers $`M_{g0}`$ and $`M_{c0}`$, respectively, it is possible to obtain exact, analytical formulae for the critical curves that form the border of the region of multiple solutions in the parameter space.

### 4.1 Post-Subshock Flow

The existence of multiple dynamical solutions is connected with the presence in the flow of a discontinuous subshock mediated by the pressure of the gas. We can therefore derive critical criteria related to the multiple-solution phenomenon by focusing on the nature of the flow in the post-subshock region, assuming that a subshock exists in the flow. As we demonstrate in the Appendix, the energy, momentum, and particle fluxes for the cosmic rays and the background gas are independently conserved as the flow crosses the subshock. This implies that the quantities associated with the gas satisfy the usual Rankine-Hugoniot jump conditions, as pointed out by DV. Far downstream from the subshock, the flow must certainly relax into a gradient-free condition if a steady state prevails. In fact, it is possible to demonstrate that the entire post-subshock region is gradient-free, so that $`u=`$constant downstream from the subshock. This has already been shown by DV using a geometrical approach, but it can also be easily established using the following simple mathematical argument. First we combine the dynamical equation (2.22) with the definition of $`g(u)`$ given by equation (2.25) to obtain the alternative form

$$\frac{du}{dx}=(\gamma _c1)\frac{v_0}{\overline{\kappa }}\frac{g(u)}{M_g^21}.$$
$`(4.1)`$
In the post-subshock gas, $`M_g<1`$ and therefore $`M<1`$ regardless of the value of $`M_c`$, since $`M^2=M_g^2+M_c^2`$. It follows from equation (2.29) that $`dg/du>0`$ downstream from the subshock. Referring to Figure 4, we wish to prove that the subshock transition must take the velocity $`u`$ directly to the gradient-free asymptotic downstream root denoted by $`u_1`$, so that $`g=0`$ on the immediate downstream side of the subshock. To develop the proof, let us suppose instead that $`g>0`$ in the post-subshock gas, corresponding to the post-subshock velocity $`u_a>u_1`$ in Figure 4. In this case, equation (4.1) implies that $`du/dx>0`$ in the downstream region, so that $`u`$ increases along the flow direction, evolving away from the gradient-free root $`u_1`$ in the post-subshock flow. Conversely, if $`g<0`$ in the post-subshock gas (corresponding to the velocity $`u_b<u_1`$ in Figure 4), then $`du/dx<0`$ in the downstream region and consequently $`u`$ decreases, again evolving away from the gradient-free root $`u_1`$. Hence if the flow is to be stationary, then $`u`$ must jump directly to the value $`u_1`$, and the entire post-subshock region must therefore be gradient-free. This conclusion is valid within the context of the “standard” two-fluid model studied by DV, but Zank, Webb, & Donohue (1993) suggest that it may be violated in models including injection.

### 4.2 Global Flow Structure

We can derive an expression for $`g(u)`$ suitable for use in the post-subshock region by employing equation (2.21) to eliminate $`M_g`$ in equation (2.25). This yields

$$g(u)=\left(\frac{1}{2}\mathrm{\Gamma }_c\right)u^2+\mathrm{\Gamma }_cu+\left(\mathrm{\Gamma }_g\mathrm{\Gamma }_c\right)\frac{u_+^{1+\gamma _g}u^{1\gamma _g}}{\gamma _gM_{g+}^2},$$
$`(4.2)`$
where we have adopted the immediate post-subshock values $`u_+`$ and $`M_{g+}`$ as the fiducial quantities in the smooth section of the flow downstream from the subshock. Since we have already established that the post-subshock flow is gradient-free, we can obtain an equation satisfied by the post-subshock velocity $`u_+`$ by writing

$$g(u_+)=0=\left(\frac{1}{2}\mathrm{\Gamma }_c\right)u_+^2+\mathrm{\Gamma }_cu_++\left(\mathrm{\Gamma }_g\mathrm{\Gamma }_c\right)\frac{u_+^2}{\gamma _gM_{g+}^2}.$$
$`(4.3)`$
The gradient-free nature of the post-subshock flow also trivially implies

$$u_1=u_+,M_{g1}=M_{g+},$$
$`(4.4)`$
where $`M_{g1}`$ is the asymptotic downstream gas Mach number associated with the downstream velocity $`u_1`$. Equation (4.3) can be interpreted as a relation for the immediate pre-subshock velocity $`u_{}`$ by utilizing the standard Rankine-Hugoniot jump conditions (Landau & Lifshitz 1975)

$$\frac{u_+}{u_{}}=\frac{2+(\gamma _g1)M_g^2}{(\gamma _g+1)M_g^2},M_{g+}^2=\frac{2+(\gamma _g1)M_g^2}{1\gamma _g+2\gamma _gM_g^2},$$
$`(4.5)`$
where the pre-subshock gas Mach number $`M_g`$ is related to $`u_{}`$ and $`M_{g0}`$ via

$$M_g^2=M_{g0}^2u_{}^{1+\gamma _g},$$
$`(4.6)`$
which follows from equation (2.21).

By using equations (4.5) and (4.6) to eliminate $`u_+`$, $`M_{g+}`$, and $`M_g`$, we can transform equation (4.3) into a new equation for the pre-subshock velocity $`u_{}`$, which we write symbolically as

$$h(u_{})=0,$$
$`(4.7)`$
where

$$h(u_{})[(\frac{1}{2}\frac{\gamma _g+1}{\gamma _g1}\mathrm{\Gamma }_r\frac{\mathrm{\Gamma }_r}{\gamma _g}M_{g0}^2u_{}^{1\gamma _g})u_{}^2$$

$$+\mathrm{\Gamma }_ru_{}\left]\right(\frac{\gamma _g1+2M_{g0}^2u_{}^{1\gamma _g}}{\gamma _g+1}).$$
$`(4.8)`$
Recall that the constants $``$ and $``$ are functions of $`M_{c0}`$ and $`M_{g0}`$ by virtue of equations (2.27) and (2.28). In addition to satisfying the condition $`h=0`$, acceptable roots for the pre-subshock velocity $`u_{}`$ must also exceed the critical velocity $`u_s`$. This is because the flow must be supersonic with respect to the gas before crossing the subshock, if one exists. Equation (4.7) allows us to solve for the pre-subshock velocity $`u_{}`$ as a function of the upstream cosmic-ray and gas Mach numbers $`M_{c0}`$ and $`M_{g0}`$, respectively. In general, $`u_{}`$ is a multi-valued function of $`M_{c0}`$ and $`M_{g0}`$, and this results in the possibility of several distinct subshock solutions in certain regions of the parameter space. Fortunately, additional information is also available that can be utilized to calculate the shape of the critical curve in $`(M_{g0},M_{c0})`$ space bordering the region of multiple subshock solutions. The nature of this information becomes clear when one examines the topology of the function $`h(u_{})`$ as depicted in Figures 5 and 6 for $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$. We consider a sequence of situations with $`M_{g0}`$ held fixed and $`M_{c0}`$ gradually increasing from the minimum value $`M_{c,\mathrm{min}}`$ given by equation (2.32), corresponding to the limit $`M_0=1`$. Note that since $`M_{g0}`$ is held constant, the critical velocity $`u_s`$ also remains constant according to equation (3.1). Two qualitatively different behaviors are observed depending on whether or not $`M_{g0}`$ exceeds $`M_{gA}`$, where $`M_{gA}=12.28`$ is the critical upstream gas Mach number for smooth flow calculated using equation (3.8) with $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$.

In Figure 5a we plot $`h`$ as a function of $`u_{}`$ and $`M_{c0}`$ for $`M_{g0}=8`$, which yields for the critical velocity $`u_s=0.210`$. In this case the minimum upstream cosmic-ray Mach number is $`M_{c,\mathrm{min}}=1.008`$. Note that initially, for small $`M_{c0}`$, there is one root for $`u_{}`$ corresponding to the single crossing of the line $`h=0`$. Since this root does not exceed $`u_s`$, a subshock solution is impossible and instead the flow must be globally smooth as discussed in § 3. The choice $`M_{g0}=8`$ satisfies the condition $`M_{g0}<M_{gA}`$, and therefore as $`M_{c0}`$ increases, the root for $`u_{}`$ eventually equals $`u_s`$, which occurs when $`M_{c0}=M_{cA}=4.25`$. In this case the subshock is located at the asymptotic downstream limit of the flow, and therefore it is identical to the gas sonic point. Equations (3.1) and (4.6) indicate that the pre-subshock gas Mach number $`M_g=1`$ as expected.

As we continue to increase $`M_{c0}`$ beyond the critical value $`M_{cA}`$ in Figure 5a, the root for $`u_{}`$ exceeds $`u_s`$, and therefore the smooth solution is replaced by a subshock solution. We refer to this solution as the “primary” subshock solution. Since $`M_g>1`$ in this region of the parameter space, the subshock plays a significant role in modifying the flow. If $`M_{c0}`$ is increased further, the primary root for $`u_{}`$ increases slowly, the concave-down shape changes, and a new peak begins to emerge at large $`u_{}`$. The peak touches the line $`h=0`$ when $`M_{c0}=41.65`$, and therefore at this point a new subshock root appears with the value $`u_{}=0.944`$. The development of this new root can be clearly seen in Figure 5b, where we replot $`h`$ at a much smaller scale. As $`M_{c0}`$ continues to increase, the peak continues to rise, and the new $`u_{}`$ root bifurcates into two roots. Since the two new roots for $`u_{}`$ are larger than the primary root, the corresponding pre-subshock gas Mach numbers are also larger, and therefore the two new subshocks are stronger than the primary subshock. We conclude that in this region of the $`(M_{g0},M_{c0})`$ parameter space, three discontinuous subshock solutions are possible, although only one can occur in a given situation.

Another example is considered in Figure 6a, where we plot $`h`$ as a function of $`u_{}`$ and $`M_{c0}`$ for $`M_{g0}=13`$, which yields $`u_s=0.146`$ and $`M_{c,\mathrm{min}}=1.003`$. In this case, $`M_{g0}>M_{gA}`$, and consequently there is always a root for $`u_{}`$ below $`u_s`$, indicating that globally smooth flow is possible for all values of $`M_{c0}`$. Hence the “primary” subshock solution never appears in this example. However, for sufficiently large $`M_{c0}`$, a peak develops in $`h`$ just as in Figure 5b. This peak rises with increasing $`M_{c0}`$ and eventually crosses the line $`h=0`$ at $`u_{}=0.981`$ when $`M_{c0}=116`$, corresponding to the appearance of a new physically acceptable subshock root for $`u_{}`$. This new root bifurcates into two roots as $`M_{c0}`$ is increased, which can be clearly seen in Figure 6b, where $`h`$ is replotted on a much smaller scale. It follows that in this region of $`(M_{g0},M_{c0})`$ space, two discontinuous solutions are possible in addition to a single globally smooth solution. It is apparent from Figures 5 and 6 that the onset of multiple solutions is connected with the vanishing of $`h`$ coupled with the additional, simultaneous condition

$$\left(\frac{h}{u_{}}\right)_{M_{g0},M_{c0}}=0,$$
$`(4.9)`$
which supplements equation (4.7).

### 4.3 Critical Mach Numbers for Multiple Solutions

Equations (4.7) and (4.9) can be manipulated to obtain explicit expressions for the critical upstream gas and cosmic-ray Mach numbers corresponding to the onset of multiple subshock solutions as functions of the pre-subshock velocity $`u_{}`$. These critical Mach numbers are denoted by $`M_{gB}`$ and $`M_{cB}`$, respectively. The logical procedure for obtaining the relations is straightforward, although the algebra required is somewhat tedious. The first step in the process is to solve equations (4.7) and (4.9) individually to derive two separate expressions for $`M_{cB}`$. Equation (4.7) yields

$$M_{cB}=M_{gB}\left(\frac{F_1+F_2M_{gB}^2}{F_3+F_4M_{gB}^2+F_5M_{gB}^4}\right)^{1/2},$$
$`(4.10\mathrm{a})`$
where

$$F_14\gamma _g(\gamma _g1)u_{}^{\gamma _g},$$
$`(4.10\mathrm{b})`$
$$F_22\gamma _g(\gamma _g1)\left[1\gamma _g+(\gamma _g1)u_{}\right]u_{}^{2\gamma _g},$$
$`(4.10\mathrm{c})`$
$$F_34\gamma _c(\gamma _g1)(u_{}^{\gamma _g}1),$$
$`(4.10\mathrm{d})`$
$$F_42u_{}^{\gamma _g}[\gamma _cu_{}(u_{}^{\gamma _g}1)+\gamma _g(\gamma _g+1)(u_{}^{\gamma _g}u_{})$$

$$+\gamma _c\gamma _g^2(u_{}1)(u_{}^{\gamma _g}2)\gamma _g\gamma _c(25u_{}+u_{}^{\gamma _g}+2u_{}^{1+\gamma _g})],$$
$`(4.10\mathrm{e})`$
$$F_5\gamma _g(\gamma _g1)(u_{}1)u_{}^{2\gamma _g}\left[(\gamma _g+1)(\gamma _c1)+(1+\gamma _g3\gamma _c+\gamma _g\gamma _c)u_{}\right].$$
$`(4.10\mathrm{f})`$
The solution to equation (4.9) is given by

$$M_{cB}=M_{gB}\left(\frac{G_1+G_2M_{gB}^2}{G_3+G_4M_{gB}^2+G_5M_{gB}^4}\right)^{1/2},$$
$`(4.11\mathrm{a})`$
where

$$G_12\gamma _g^2u_{}^{\gamma _g},$$
$`(4.11\mathrm{b})`$
$$G_2\gamma _g(\gamma _g1)u_{}^{1+2\gamma _g},$$
$`(4.11\mathrm{c})`$
$$G_32\gamma _g\gamma _c(2+u_{}^{\gamma _g}),$$
$`(4.11\mathrm{d})`$
$$G_4u_{}^{\gamma _g}\left[2\gamma _g^2\gamma _c+(\gamma _g+\gamma _c5\gamma _g\gamma _c+\gamma _g^2+2\gamma _g^2\gamma _c)u_{}+\gamma _c(\gamma _g1)u_{}^{1+\gamma _g}\right],$$
$`(4.11\mathrm{e})`$
$$G_5\gamma _gu_{}^{1+2\gamma _g}\left[\gamma _c(\gamma _g1)+(1+\gamma _g3\gamma _c+\gamma _g\gamma _c)u_{}\right].$$
$`(4.11\mathrm{f})`$
The similarity of the dependences on $`M_{gB}`$ in equations (4.10) and (4.11) for $`M_{cB}`$ suggests that we can derive an exact solution for $`M_{gB}`$ as a function of $`u_{}`$ by equating these two expressions. The result obtained is

$$M_{gB}=\left(\frac{2T_12^{1/3}T_5T_6^{1/3}+2^{1/3}T_6^{1/3}}{3T_2}\right)^{1/2},$$
$`(4.12\mathrm{a})`$
where

$$T_1(\gamma _g1)u_{}^{2\gamma _g}\{\gamma _g^2(1+\gamma _g)(1+\gamma _c)+\gamma _g(2+2\gamma _g7\gamma _c+4\gamma _g\gamma _c\gamma _g^2\gamma _c)u_{}^2$$

$$+(\gamma _g+1)[(\gamma _g\gamma _c)u_{}^{\gamma _g}+\gamma _g+\gamma _c5\gamma _g\gamma _c+\gamma _g^2+2\gamma _g^2\gamma _c]u_{}\},$$
$`(4.12\mathrm{b})`$
$$T_2\gamma _g(\gamma _g1)u_{}^{1+3\gamma _g}[(\gamma _c+1)(\gamma _g^21)2(\gamma _g+1)(1+\gamma _g3\gamma _c+\gamma _g\gamma _c)u_{}$$

$$+(\gamma _g1)(1+\gamma _g3\gamma _c+\gamma _g\gamma _c)u_{}^2],$$
$`(4.12\mathrm{d})`$
$$T_32\gamma _c(1\gamma _g^2)+(5\gamma _c1\gamma _g5\gamma _g\gamma _c+2\gamma _g^2\gamma _c)u_{}+(\gamma _g\gamma _c\gamma _g\gamma _c+\gamma _g^2)u_{}^{\gamma _g},$$
$`(4.12\mathrm{c})`$
$$T_4216\gamma _g\gamma _c(\gamma _g1)T_2^216T_1^372\gamma _gT_1T_2T_3u_{}^{\gamma _g},$$
$`(4.12\mathrm{e})`$
$$T_54T_1^212\gamma _gT_2T_3u_{}^{\gamma _g},$$
$`(4.12\mathrm{f})`$
$$T_6T_4+\left(T_4^2+4T_5^3\right)^{1/2}.$$
$`(4.12\mathrm{g})`$
Equation (4.12) can be used to evaluate the critical upstream gas Mach number $`M_{gB}`$ as a function of the pre-subshock velocity $`u_{}`$. Once $`M_{gB}`$ is determined, we can calculate the corresponding critical upstream cosmic-ray Mach number $`M_{cB}`$ using either equation (4.10) or equation (4.11), which both yield the same result. Hence equations (4.10), (4.11), and (4.12) provide a direct means for calculating $`M_{gB}`$ and $`M_{cB}`$ as exact functions of $`u_{}`$. As an example, setting $`\gamma _g=5/3`$, $`\gamma _c=4/3`$, and $`u_{}=0.944`$ yields $`M_{gB}=8`$ and $`M_{cB}=41.65`$, in agreement with the results plotted in Figure 5. Likewise, setting $`\gamma _g=5/3`$, $`\gamma _c=4/3`$, and $`u_{}=0.981`$ yields $`M_{gB}=13`$ and $`M_{cB}=116`$, in agreement with Figure 6.

Although our expressions for $`M_{gB}`$ and $`M_{cB}`$ simplify considerably in the special case $`\gamma _g=5/3`$, $`\gamma _c=4/3`$, we have chosen to derive results valid for general (but constant) values of $`\gamma _g`$ and $`\gamma _c`$ in order to obtain a full understanding of the sensitivity of the critical Mach numbers to variations in the adiabatic indices. While admittedly somewhat complex, these equations can nonetheless be evaluated using a hand calculator, and replace the requirement of utilizing the root-finding techniques employed in most previous investigations of the multiple-solution phenomenon in the two-fluid model.

In Figure 7 we plot the critical curves for the occurrence of multiple solutions using various values of $`\gamma _g`$ and $`\gamma _c`$. This is accomplished by evaluating $`M_{gB}`$ and $`M_{cB}`$ as parametric functions of the pre-subshock velocity $`u_{}`$ using equation (4.12) along with either equation (4.10) or (4.11). The critical curves denote the boundary of the wedge-shaped multiple-solution region. Inside this region, two new subshock solutions are possible, along with either the primary subshock solution or the globally smooth solution. Outside this region, the flow is either described by the primary subshock solution or else it is globally smooth. The lower-left-hand corner of the multiple-solution region curves to the left and culminates in a sharp cusp. The presence of the cusp implies that there is a minimum value of $`M_{g0}`$, below which multiple solutions are never possible for any value of $`M_{c0}`$. Note that the multiple-solution region becomes very narrow as the values of $`\gamma _g`$ and $`\gamma _c`$ approach each other, suggesting that the physically allowed solutions converge on a single solution as the thermodynamic properties of the two populations of particles (background gas and cosmic rays) become more similar to each other. In the limit $`\gamma _g=\gamma _c`$, multiple solutions are not allowed at all, and the single available solution is either smooth or discontinuous depending on the values of $`M_{g0}`$ and $`M_{c0}`$.

In Figure 8 we plot all of the various critical curves for the physically important case of an ultrarelativistic cosmic-ray distribution ($`\gamma _c=4/3`$) combined with a nonrelativistic background gas ($`\gamma _g=5/3`$). This is the most fully self-consistent example of the two-fluid model, since in this case we expect that the adiabatic indices will remain constant throughout the flow, as is assumed in the model (see eq. \[2.9\]). The area of overlap between the multiple-solution region and the smooth-solution region implies the existence of four distinctly different domains within the $`(M_{g0},M_{c0})`$ parameter space. Within Domain I, which lies outside both the multiple-solution region and the smooth-solution region, the flow must be discontinuous, with exactly one (i.e., the primary subshock) solution available. Domain II lies inside the multiple-solution region and outside the smooth-solution region, and therefore within this area of the parameter space, three distinct subshock solutions are possible, while globally smooth flow is impossible. Domain III is formed by the intersection of the multiple-solution region and the smooth-solution region, and therefore two discontinuous subshock solutions are available as well as one globally smooth solution. Finally, Domain IV lies outside the multiple-solution region and within the smooth-flow region, and therefore one globally smooth flow solution is possible. The existence of Domain IV is consistent with our expectation that smooth flow will occur for sufficiently large values of the upstream cosmic-ray pressure $`P_{c0}`$ due to diffusion of the cosmic rays.

If we consider a trajectory through the $`(M_{g0},M_{c0})`$ parameter space that crosses into the multiple-solution region, then the sequence of appearance of the subshock roots for $`u_{}`$ depends on whether the boundary is crossed through the “top” or “bottom” arcs of the wedge. We illustrate this phenomenon for $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$ in Figure 8, where we consider two possible paths approaching the point $`P`$ inside the multiple-solution region, starting outside at either points $`Q`$ or $`R`$. The two paths cross the multiple-solution boundary on different sides of the wedge. In Figure 9 we plot $`h(u_{})`$ along the segment $`RP`$ (with $`M_{g0}=6.5`$), and in Figure 10 we plot $`h(u_{})`$ along the segment $`QP`$ (with $`M_{c0}=45`$). Note that the primary subshock root for $`u_{}`$ already exists at points $`Q`$ and $`R`$ since they both lie inside Domain I. When the lower section of the multiple-solution boundary is crossed along segment $`RP`$ ($`M_{c0}=M_{cB}=25.9`$ in Fig. 9), two new roots for $`u_{}`$ appear at larger values than the primary root, which is the same sequence observed in Figure 5. However, when the upper section of the boundary is crossed along segment $`QP`$ ($`M_{g0}=M_{gB}=5.83`$ in Fig. 10), the two new roots for $`u_{}`$ appear at smaller values than the primary root. Hence the order of appearance of the subshock roots for $`u_{}`$ is different along each path. Despite this path dependence, the actual set of roots obtained at point $`P`$ is the same regardless of the approach path taken, and therefore there is no ambiguity regarding the possible subshock solutions available at any point in the $`(M_{g0},M_{c0})`$ parameter space.

In Figure 11 we present summary plots that include all of the various critical curves derived in §§ 3 and 4 for several different values of $`\gamma _g`$ and $`\gamma _c`$. Note that the area of overlap between the multiple-solution region and the smooth-solution region observed when $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$ rapidly disappears when the difference between $`\gamma _g`$ and $`\gamma _c`$ is reduced, due to the increasing similarity between the thermodynamic properties of the cosmic rays and the background gas. Ko, Chan, & Webb (1997) and Bulanov & Sokolov (1984) have obtained similar parameter space plots depicting the critical curves for smooth flow and for the onset of multiple solutions. The parameters used to describe the incident flow conditions are $`(M_{g0},Q_0)`$ in the case of Ko, Chan, & Webb (1997) and $`(M_0,Q_0)`$ for Bulanov & Sokolov (1984), where $`Q_0`$ and $`M_0`$ refer to the incident pressure ratio and total Mach number, respectively (see eqs. \[1.1\] and \[1.2\]). When these parameters are employed directly instead of the quantities $`(M_{g0},M_{c0})`$ utilized in our work, the critical curves must be determined by numerical root-finding. However, if desired, equations (1.5) can be used to transform our exact solutions for the critical curves from the $`(M_{g0},M_{c0})`$ space to the $`(M_{g0},Q_0)`$ and $`(M_0,Q_0)`$ spaces. This procedure was used to obtain the results depicted in Figure 12, which are identical to those presented by Ko, Chan, & Webb (1997) and Bulanov & Sokolov (1984).

## 5 EXAMPLES

Our interpretation of the various critical Mach numbers derived in §§ 3 and 4 can be verified by performing specific calculations of the flow dynamics for a few different values of the upstream Mach numbers $`M_{g0}`$ and $`M_{c0}`$. We shall focus here on cases with $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$, but this is not essential. Along smooth sections of the flow, we can calculate the velocity profile by integrating the dynamical equation obtained by combining equations (2.21) and (4.1), which yields

$$\frac{du}{d\tau }=(\gamma _c1)g(u)\left[M_g^2\left(\frac{u}{u_{}}\right)^{1\gamma _g}1\right]^1,$$
$`(5.1)`$
where

$$g(u)=\left(\frac{1}{2}\mathrm{\Gamma }_c\right)u^2+\mathrm{\Gamma }_cu+\left(\mathrm{\Gamma }_g\mathrm{\Gamma }_c\right)\frac{u_{}^{1+\gamma _g}u^{1\gamma _g}}{\gamma _gM_g^2},$$
$`(5.2)`$
and we have introduced the new spatial variable

$$\tau (x)v_0_0^x\frac{dx^{}}{\overline{\kappa }(x^{})}.$$
$`(5.3)`$
The quantities $`M_g`$ and $`u_{}`$ are fiducial parameters measured at an arbitrary location along the smooth section of interest, and the constants $``$ and $``$ are functions of the incident Mach numbers $`M_{g0}`$ and $`M_{c0}`$ via equations (2.27) and (2.28). As discussed in § 3, if $`M_{c0}<M_{cA}`$, then the flow is everywhere supersonic with respect to the gas, and therefore we can use equation (5.1) to describe the structure of the entire flow by setting the fiducial parameters $`M_g`$ and $`u_{}`$ equal to the asymptotic upstream quantities $`M_{g0}`$ and $`u_0=1`$, respectively. Conversely, if $`M_{c0}>M_{cA}`$, then the flow must contain a discontinuous, gas-mediated subshock. In this case the upstream region is governed by equation (5.1) with $`M_g=M_{g0}`$ and $`u_{}=u_0=1`$, and the downstream region is governed by equation (5.1) with $`M_g=M_{g+}`$ and $`u_{}=u_+`$. The pre-subshock and post-subshock regions are linked by the jump conditions given by equations (4.5), which determine $`M_{g+}`$ and $`u_+`$ and yield $`g=0`$ in the entire downstream region as expressed by equation (4.3). In order to illustrate the various possible behaviors predicted by our critical conditions, we shall present one example from each of the four domains discussed in § IV (see Fig. 8) calculated using $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$.

In Figure 13 we plot the velocity profile $`u(\tau )`$ obtained by integrating the dynamical equation (5.1) with $`M_{g0}=4`$ and $`M_{c0}=4`$. According to Figure 8, this point lies within Domain I, and therefore we expect to find a single subshock solution, and no globally smooth solution. Equation (2.25) confirms that in this case the downstream root for $`g(u)`$ is $`u=0.27`$, and therefore smooth flow is impossible since this does not exceed the critical velocity $`u_s=0.35`$. Analysis of equation (4.7) yields a single acceptable value for the pre-subshock velocity root $`u_{}=0.50`$, with an associated pre-subshock gas Mach number $`M_g=1.60`$. Figure 13 also includes plots of the gas Mach number $`M_g(\tau )`$ obtained using equation (2.21), the dimensionless gas pressure $`𝒫_g(\tau )`$ obtained using equation (2.20), and the dimensionless cosmic-ray pressure $`𝒫_c(\tau )`$ obtained using equation (2.16). The first three quantities exhibit a jump at the subshock, whereas the cosmic-ray pressure is continuous since a jump in this quantity would imply an infinite energy flux. The overall compression ratio is $`1/u_1=3.66,`$ and the pressures increase by the factors $`𝒫_{g1}/𝒫_{g0}=9.27`$ and $`𝒫_{c1}/𝒫_{c0}=9.89`$ between the asymptotic upstream and downstream regions. Recall that $`𝒫_g`$ and $`𝒫_c`$ express the pressures of the two species divided by the upstream ram pressure of the gas. Hence in this example the two species each absorb comparable fractions of the upstream ram pressure. Note that the increase in the cosmic-ray pressure is entirely due to the smooth part of the transition upstream from the discontinuous subshock, while the gas pressure experiences most of its increase in crossing the subshock.

In Figure 14 we set $`M_{g0}=6`$ and $`M_{c0}=60`$, corresponding to Domain II in Figure 8. We therefore expect to find three distinct, physically acceptable subshock solutions associated with these parameter values. Analysis of equation (2.25) verifies that no smooth solutions are possible in this case because the downstream root $`u=0.18`$ is less than the critical velocity $`u_s=0.26`$. Equation (4.7) indicates that the three acceptable roots for the pre-subshock velocity are $`u_{}=0.53,0.72,0.99,`$ with associated pre-subshock gas Mach numbers $`M_g=2.54,3.84,5.94,`$ respectively. The gas and cosmic-ray pressures increase by the factors $`𝒫_{g1}/𝒫_{g0}=23,32,44,`$ and $`𝒫_{c1}/𝒫_{c0}=2126,1306,37,`$ and the overall compression ratios are $`1/u_1=5.20,4.64,3.72`$. Note that the solution with the largest cosmic-ray pressure has the largest overall compression ratio and the weakest subshock. This is due to the fact that the cosmic-ray pressure is amplified in the smooth-flow region upstream from the subshock, and this region is most extended when the flow does not encounter a subshock until the smallest possible value of $`M_g`$. Conversely, the solution with the largest gas pressure is the one with the strongest subshock and the smallest overall compression ratio.

In Figure 15 we use the upstream parameters $`M_{g0}=12.5`$ and $`M_{c0}=200`$. This point lies within Domain III in Figure 8, and therefore we expect to find two distinct subshock solutions along with one globally smooth solution. Equation (2.25) confirms that in this case a smooth solution is possible since the downstream root $`u=0.152`$ exceeds the critical velocity $`u_s=0.150`$. Analysis of equation (4.7) yields two acceptable values for the pre-subshock velocity given by $`u_{}=0.962,0.997`$. The corresponding pre-subshock gas Mach numbers are $`M_g=11.9,12.4,`$ respectively. In the discontinuous (subshock) solutions the gas and cosmic-ray pressures increase by the factors $`𝒫_{g1}/𝒫_{g0}=188,194,`$ and $`𝒫_{c1}/𝒫_{c0}=2011,170,`$ and the overall compression ratios are $`1/u_1=4.07,3.94`$. The subshocks in this example are strong, and therefore in the solutions containing a subshock most of the deceleration is due to the buildup of the pressure of the background gas, rather than the cosmic-ray pressure. Hence both of the subshock solutions are gas-dominated. In the globally smooth solution, which is cosmic-ray dominated, the pressures increase by the factors $`𝒫_{g1}/𝒫_{g0}=23`$ and $`𝒫_{c1}/𝒫_{c0}=40702`$, and the compression ratio is $`1/u_1=6.57`$. In this case the cosmic rays absorb almost all of the ram pressure of the upstream gas.

Finally, in Figure 16 we set $`M_{g0}=14`$ and $`M_{c0}=20`$. According to Figure 8, this point lies within Domain IV, and therefore we expect that only a single, globally-smooth solution exists for these upstream parameters. This prediction is verified by equation (4.7), which confirms that no acceptable subshock roots for $`u_{}`$ exists. Furthermore, equation (2.25) indicates that smooth flow is possible since the downstream root $`u=0.15`$ exceeds the critical velocity $`u_s=0.14`$. Hence the only acceptable solution is globally smooth, with the pressure increases given by $`𝒫_{g1}/𝒫_{g0}=23`$ and $`𝒫_{c1}/𝒫_{c0}=417`$. The overall compression ratio is $`1/u_1=6.56`$ and the flow is cosmic-ray dominated.

## 6 DISCUSSION

In this paper we have obtained a number of new analytical results describing the critical behavior of the two-fluid model for cosmic-ray modified shocks. It is well known that in this model, up to three distinct solutions are possible for a given set of upstream boundary conditions. The behaviors of the various solutions can be quite diverse, including flows that are smooth everywhere as well as flows that contain a discontinuous, gas-mediated subshock. The traditional approach to the problem of determining the types of possible solutions, employed by Ko, Chan, & Webb (1997) and Bulanov & Sokolov (1984), is based on stating the upstream boundary conditions in terms of the incident total Mach number $`M_0`$ and the incident pressure ratio $`Q_0`$ (see eqs. \[1.1\] and \[1.2\]). In this approach the determination of the available solution types requires several steps of root-finding, and there is no possibility of obtaining analytical expressions for the critical relationships.

The analysis presented here utilizes a fresh approach based upon a new parameterization of the boundary conditions in terms of the upstream gas and cosmic-ray Mach numbers $`M_{g0}`$ and $`M_{c0}`$, respectively. The analytical results we have obtained in §§ 3 and 4 for the critical upstream Mach numbers expressed by equations (3.4), (4.11), and (4.12) provide for the first time a systematic classification of the entire parameter space for the two-fluid model, which remains one of the most powerful and practical means available for studying the problem of cosmic-ray modified shocks. These expressions eliminate the need for complex root-finding procedures in order to understand the possible flow dynamics for a given set of upstream boundary conditions, and are made possible by the symmetry between the gas and cosmic-ray parameters as they appear in the expressions describing the asymptotic upstream and downstream states of the flow. We have compared our quantitative results with those of Ko, Chan, & Webb (1997) and Bulanov & Sokolov (1984), and they are found to be consistent. The results are valid for arbitrary (but constant) values of the gas and cosmic-ray adiabatic indices $`\gamma _g`$ and $`\gamma _c`$, respectively. In § 5 we have presented numerical examples of flow structures obtained in each of the four parameter space domains defined in Figure 8. These examples verify the predictions made using our new expressions for the critical Mach numbers, and confirm that the largest overall compression ratios are obtained in the globally-smooth, cosmic-ray dominated cases.

The existence of multiple distinct solutions for a single set of upstream boundary conditions demands that we include additional physics in order to determine which solution is actually realized in a given situation. This question has been addressed by numerous authors using various forms of stability analysis as well as fully time-dependent calculations. DV speculated that when three distinct solutions are allowed (in Domains II and III of the parameter space plotted in Fig. 8), the solution with the intermediate value of the cosmic-ray pressure $`𝒫_c`$ will be unstable. The argument is based on the idea that if the cosmic-ray pressure were to increase slightly due to a small perturbation, then the gas would suffer additional deceleration, leading to a further increase in the cosmic-ray pressure. This nonlinear process would drive the flow towards the steady-state solution with the largest value of $`𝒫_c`$. Conversely, a small decrease in the cosmic-ray pressure would decrease the deceleration, leading to a smaller value for the cosmic-ray pressure. In this case the flow would be driven towards the steady-state solution with the smallest value of $`𝒫_c`$. Recently, Mond & Drury (1998) have suggested that this type of behavior may be realized as a consequence of a corrugational instability. Other authors (e.g., Drury & Falle 1986; Kang, Jones, & Ryu 1992; Zank, Axford, & McKenzie 1990; Ryu, Kang, & Jones 1993) have argued that the globally smooth, cosmic-ray dominated solutions are unstable to the evolution of MHD waves in certain situations. Jones & Ellison (1991) suggest that even when formally stable, the smooth solutions may not be realizable in nature. On the other hand, Donohue, Zank, & Webb (1994) report time-dependent simulations which seem to indicate that the smooth, cosmic-ray dominated solution is indeed the preferred steady-state solution in certain regions of the parameter space. Hence, despite the fact that much effort has been expended in analyzing the stability properties of cosmic-ray modified shocks, there is still no clear consensus regarding which of the steady-state solutions (if any) is stable and therefore physically observable for an arbitrary set of upstream conditions.

In light of the rather contradictory state of affairs regarding the stability properties of the various possible dynamical solutions, we propose a new form of entropy-based stability analysis. In this method, the entropy of the cosmic-ray distribution is calculated by first solving the transport equation (2.1) for the cosmic-ray distribution $`f`$ and then integrating to obtain the Boltzmann entropy per cosmic ray,

$$\mathrm{\Sigma }_ck_0^{\mathrm{}}4\pi p^2F\mathrm{ln}Fdpk\mathrm{ln}\mathrm{}^3+k\mathrm{ln}V,$$
$`(6.1)`$
where $`k`$ is Boltzmann’s constant, $`n_c`$ is the cosmic-ray number density, $`\mathrm{}`$ is Planck’s constant, $`V`$ is the system volume, and $`Ff/n_c`$. According to equation (2.3), $`4\pi p^2F(p,x)dp`$ gives the probability that a randomly selected cosmic ray at location $`x`$ has momentum in the range between $`p`$ and $`p+dp`$. The cosmic-ray entropy density $`S_c`$ is computed using

$$S_c=n_c\mathrm{\Sigma }_ckn_c\mathrm{ln}\left(n_cV\right)+kn_c,$$
$`(6.2)`$
where the final two terms stem from the fundamental indistinguishability of the cosmic ray particles (of like composition), and is necessary in order to avoid the Gibbs paradox (Reif 1965). The inconvenient reference to the system volume $`V`$ can be removed by combining equations (6.1) and (6.2) to obtain the equivalent expression

$$S_ckn_c_0^{\mathrm{}}4\pi p^2F\mathrm{ln}Fdpkn_c\mathrm{ln}\left(n_c\mathrm{}^3\right)+kn_c.$$
$`(6.3)`$
The total entropy per particle $`\mathrm{\Sigma }_{\mathrm{tot}}`$ for the combined gas-particle system is calculated using

$$\mathrm{\Sigma }_{\mathrm{tot}}=\frac{S_c+S_g}{n_c+n_g},$$
$`(6.4)`$
where $`S_g`$ and $`n_g`$ denote the entropy density and the number density of the background gas, respectively. One may reasonably hypothesize that the state with the largest value for $`\mathrm{\Sigma }_{\mathrm{tot}}`$ will be the preferred state in nature. This criterion may prove useful for identifying the most stable solution when multiple solutions are available. We plan to pursue this line of inquiry in future work. The results may shed new light on the structure of cosmic-ray modified shocks.

The authors acknowledge several stimulating conversations with Frank Jones and Truong Le. The authors are also grateful to the anonymous referee for several useful comments. PAB also acknowledges support and hospitality from NASA via the NASA/ASEE Summer Faculty Fellowship program at Goddard Space Flight Center.

APPENDIX

In this section we demonstrate that the cosmic-ray energy flux must be continuous across a velocity discontinuity (subshock), if one is present in the flow. This is turn implies that the subshock must be mediated entirely by the pressure of the background gas, and therefore the discontinuity is governed by the classical Rankine-Hugoniot jump conditions. The cosmic-ray energy flux in the $`\widehat{x}`$ direction is given by

$$F_c=\overline{\kappa }\frac{dU_c}{dx}+\gamma _cvU_c.$$
$`(\mathrm{A1})`$
In a steady-state, the cosmic-ray energy equation (2.7) reduces to

$$v\frac{dU_c}{dx}=\gamma _cU_c\frac{dv}{dx}+\frac{d}{dx}\left(\overline{\kappa }\frac{dU_c}{dx}\right),$$
$`(\mathrm{A2})`$
which can be combined with equation (A1) to express the derivative of $`F_c`$ as

$$\frac{dF_c}{dx}=(\gamma _c1)v\frac{dU_c}{dx}.$$
$`(\mathrm{A3})`$
Integration in the vicinity of the subshock located at $`x=x_0`$ yields

$$\underset{ϵ0}{lim}_{x_0ϵ}^{x_0+ϵ}\frac{dF_c}{dx}𝑑x=\underset{ϵ0}{lim}_{x_0ϵ}^{x_0+ϵ}(\gamma _c1)v\frac{dU_c}{dx}𝑑x.$$
$`(\mathrm{A4})`$
In order to avoid unphysical divergence of $`F_c`$ at the subshock, $`U_c`$ must be continuous, and therefore the integrand on the right-hand side of equation (A4) is no more singular than a step function. This implies that in the limit $`ϵ0`$, the right-hand side vanishes, leaving

$$\mathrm{\Delta }F_c\underset{ϵ0}{lim}F_c(x_0+ϵ)F_c(x_0ϵ)=0.$$
$`(\mathrm{A5})`$
Hence $`F_c`$ remains constant across the subshock. The energy, momentum, and particle fluxes of the gas and the cosmic rays are therefore independently conserved across the discontinuity. This allows us to use the standard Rankine-Hugoniot jump conditions to describe the subshock transition in equations (4.5).

FIGURE CAPTIONS

Fig. 1. – Topology of the function $`g(u)`$ given by eq. \[2.25\]. The quantities $`u_0`$ and $`u_1`$ respectively denote the upstream and downstream roots for the velocity in a globally smooth solution.

Fig. 2. – Critical upstream cosmic-ray Mach number $`M_{cA}`$ (eq. \[3.4\]) for smooth flow plotted as a function of the upstream gas Mach number $`M_{g0}`$ in the $`(M_{g0},M_{c0})`$ parameter space for (a) $`\gamma _g=5/3`$ and various values of $`\gamma _c`$ as indicated; (b) $`\gamma _c=4/3`$ and various values of $`\gamma _g`$ as indicated. Smooth flow is not possible in the region above each curve.

Fig. 3. – Critical upstream gas Mach number for smooth flow $`M_{gA}`$ (eq. \[3.8\]) plotted as a function of the gas and cosmic-ray adiabatic indices $`\gamma _g`$ and $`\gamma _c`$, respectively. When $`M_{g0}>M_{gA}`$, smooth flow is possible for any value of $`M_{c0}`$.

Fig. 4. – Schematic depiction of the function $`g(u)`$ (eq. \[2.25\]). If the flow contains a discontinuous, gas-mediated subshock, then the velocity must jump directly to the final asymptotic value $`u_1`$ in crossing the shock. Otherwise the flow is unstable (see the discussion in the text).

Fig. 5. – Function $`h(u_{})`$ (eq. \[4.8\]) is plotted for the parameters (a) $`M_{g0}=8`$, $`\gamma _g=5/3`$, $`\gamma _c=4/3`$. The value of $`M_{c0}`$ is indicated for each curve. In this example, $`M_{g0}<M_{gA}=12.28`$, and therefore the primary subshock solution appears when the low-velocity root $`u_{}>u_s`$, which occurs when $`M_{c0}>M_{cA}=4.25`$. The same function is plotted on a smaller scale in (b), where we see that two new subshock roots for $`u_{}`$ appear when $`M_{c0}>41.65`$.

Fig. 6. – Function $`h(u_{})`$ (eq. \[4.8\]) is plotted for the parameters (a) $`M_{g0}=13`$, $`\gamma _g=5/3`$, $`\gamma _c=4/3`$. The value of $`M_{c0}`$ is indicated for each curve. In this example, $`M_{g0}>M_{gA}=12.28`$, and therefore the primary subshock solution never appears. Hence smooth flow is possible for all values of $`M_{c0}`$. The same function is plotted on a smaller scale in (b), where we see that two new subshock roots for $`u_{}`$ appear when $`M_{c0}>116`$.

Fig. 7. – Critical Mach numbers $`M_{gB}`$ and $`M_{cB}`$ for the onset of multiple solutions (eqs. \[4.11\] and \[4.12\]) are plotted as parametric functions of the pre-subshock velocity $`u_{}`$ in the $`(M_{g0},M_{c0})`$ parameter space for (a) $`\gamma _g=5/3`$ and the indicated values of $`\gamma _c`$; (b) $`\gamma _c=4/3`$ and the indicated values of $`\gamma _g`$. The interior of each wedge is the multiple-solution region for the associated parameters.

Fig. 8. – Critical upstream Mach numbers for the occurrence of multiple solutions (eqs. \[4.11\] and \[4.12\]; solid line) and for smooth flow (eq. \[3.4\]; dashed line) are plotted together in the $`(M_{g0},M_{c0})`$ parameter space for the case $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$. The minimum upstream cosmic-ray Mach number required for decelerating flow is also shown (eq. \[2.32\]; dotted line). There are four distinct domains in the parameter space as discussed in the text.

Fig. 9. – Function $`h(u_{})`$ (eq. \[4.8\]) is plotted for the parameters (a) $`M_{g0}=6.5`$, $`\gamma _g=5/3`$, $`\gamma _c=4/3`$ along the segment $`RP`$ in Fig. 8. The values of the upstream cosmic-ray Mach number are $`M_{c0}=45`$ (solid line), $`M_{c0}=26`$ (dashed line), $`M_{c0}=15`$ (dotted line). When $`M_{c0}>M_{cB}=26`$, there are three distinct subshock solutions available. In panel (b) the same function is plotted on a smaller scale.

Fig. 10. – Function $`h(u_{})`$ (eq. \[4.8\]) is plotted for the parameters (a) $`M_{c0}=45`$, $`\gamma _g=5/3`$, $`\gamma _c=4/3`$ along the segment $`QP`$ in Fig. 8. The values of the upstream gas Mach number are $`M_{g0}=6.5`$ (solid line), $`M_{g0}=5.83`$ (dashed line), $`M_{g0}=4.8`$ (dotted line). When $`M_{g0}>M_{gB}=5.83`$, there are three distinct subshock solutions available. In panel (b) the same function is plotted on a smaller scale.

Fig. 11. – Critical upstream Mach numbers for the occurrence of multiple solutions (eqs. \[4.11\] and \[4.12\]; solid line) and for smooth flow (eq. \[3.4\]; dashed line) are plotted together in the $`(M_{g0},M_{c0})`$ parameter space for (a) $`\gamma _g=5/3`$, $`\gamma _c=4/3`$; (b) $`\gamma _g=5/3`$, $`\gamma _c=1.35`$; (c) $`\gamma _g=1.6`$, $`\gamma _c=4/3`$; (d) $`\gamma _g=1.6`$, $`\gamma _c=1.35`$. Also indicated is the minimum value of $`M_{c0}`$ required for decelerating flow (eq. \[2.32\]; dotted line).

Fig. 12. – Our analytical results for the critical curves generated using equations (3.4), (4.11), and (4.12) are combined with equations (1.5) to create corresponding curves in the alternative parameter spaces $`(M_0,Q_0)`$ and $`(M_{g0},Q_0)`$ employed by Bulanov & Sokolov (1984) and Ko, Chan, & Webb (1997), respectively. Panel (a), with $`\gamma _g=5/3`$ and $`\gamma _c=4/3`$, is identical to Fig. 4 of Bulanov & Sokolov (1984). Panel (b), with $`\gamma _g=2`$ and $`\gamma _c=4/3`$, is identical to Fig. 1(a) of Ko, Chan, & Webb (1997). Note that these authors generated their curves using root-finding procedures. The interpretation of the line styles is the same as in Fig. 11.

Fig. 13. – Numerical solutions for (a) $`u`$, (b) $`M_g`$, (c) $`𝒫_g`$, and (d) $`𝒫_c`$ are plotted as functions of $`\tau `$ (see eq. \[5.3\]). The solutions were obtained by integrating the dynamical eq. \[5.1\] with $`M_{g0}=4`$ and $`M_{c0}=4`$, which corresponds to Domain I in Fig. 8. In this case one discontinuous solution is possible, and smooth flow is impossible.

Fig. 14. – Same as Fig. 13, except $`M_{g0}=6`$ and $`M_{c0}=60`$, which corresponds to Domain II in Fig. 8. In this case three distinct discontinuous solutions are possible, and smooth flow is impossible. The values of the pre-subshock gas Mach number are $`M_g=2.54`$ (solid line), $`M_g=3.84`$ (dashed line), $`M_g=5.94`$ (dotted line).

Fig. 15. – Same as Fig. 13, except $`M_{g0}=12.5`$ and $`M_{c0}=200`$, which corresponds to Domain III in Fig. 8. In this case two distinct discontinuous solutions are possible in addition to one globally smooth solution (solid line). The values of the pre-subshock gas Mach number are $`M_g=11.9`$ (dashed line), $`M_g=12.4`$ (dotted line).

Fig. 16. – Same as Fig. 13, except $`M_{g0}=14`$ and $`M_{c0}=20`$, which corresponds to Domain IV in Fig. 8. In this case one globally smooth solution solution is possible, and discontinuous flow is impossible.


format: text

text:
# Weighted Evolving Networks

## Abstract

Many biological, ecological and economic systems are best described by weighted networks, as the nodes interact with each other with varying strength. However, most network models studied so far are binary, the link strength being either 0 or 1. In this paper we introduce and investigate the scaling properties of a class of models which assign weights to the links as the network evolves. The combined numerical and analytical approach indicates that asymptotically the total weight distribution converges to the scaling behavior of the connectivity distribution, but this convergence is hampered by strong logarithmic corrections.

PACS number: 84.35.+i, 05.40.-a, 02.50.Cw, 87.23.Ge

With the increased availability of detailed topological data on complex networks describing communication, social or biological systems emerges the need to develop tools to understand in general terms the origin and evolution of these complex evolving systems. The recognition that real networks are fundamentally different from the random models that dominated the mathematical literature in the past forty years lead to a surge of activity in addressing the statistical properties of these systems. Despite rapid advances in uncovering the topology of complex networks, in one aspect most models are incomplete when compared with real systems: current network models assume that all links are equivalent. But in many fields it is well known that the interaction strengths can vary widely, such variations being essential to the network’s ability to carry on its basic functions. Sociologists have repeatedly argued about the importance of assigning strengths to social links, finding that the weak links people have outside their close circle of friends play a key role in keeping the social system together. Recently, Newman has showed that assigning weights to the links between scientists allows for a better characterization of the scientific collaboration web. Similarly, there is an ongoing discussion about the importance of weak links between species in guaranteeing the stability of an ecosystem; in addition viewing the economy as a complex network of companies, the monetary value of the various business transactions linking them is a key component of the systemic characterization. Finally, many transportation networks nature designed, ranging from cardiovascular to respiratory networks, have well defined weights or flow rates assigned to the links, whose magnitude is intimately determined by the network’s topology. Despite the known importance of interaction strengths in various well known systems, there have been no attempts to model networks other than binary nets, whose links have weights $`0`$ or $`1`$.

In this paper we take a first step in the direction of a systematic study of networks with non-binary connectivities. We introduce and investigate two models that assign weights to new links as they are dynamically created, providing a prototype of a weighted evolving network. While we choose the simplest possible models, in which the weights are driven by the network connectivity only, numerical simulations indicate that the distribution of the total weight scales differently from the total connectivity. However, an analytical solution reveals that the different scaling behavior can be explained by strong logarithmic correction, and asymptotically the investigated weighted networks belong to the same universality class as their unweighted counterparts.

Weighted Scale-Free (WSF) Model : Starting from a small number ($`m_0`$) of vertices, at each time step we add a new node which links to $`m`$ existing nodes in the system. The probability that a new node $`j`$ will connect to a existing node $`i`$ is

$`\mathrm{\Pi }_i={\displaystyle \frac{k_i}{_jk_j}},`$ (1)
where $`k_i`$ is the total number of links that the node $`i`$ has. In assigning a weight to the newly established link $`ji`$, we assume that the weight $`w_{ji}(=w_{ij})`$ is proportional to $`k_i`$, i.e., more connected (and therefore more ”powerful”) nodes gain more weight. Also, one can assume that all new nodes have fairly uniform total ‘resources’ for linking to other nodes in the system, we therefore require that each new node has a fixed total weight, i.e. we normalize $`w_{ij}`$ such that the sum of the weights for the $`m`$ new links is $`_{\{i^{}\}}w_{ji^{}}=1`$, where $`\{i^{}\}`$ represents a sum over the $`m`$ existing nodes to which the new node $`j`$ is connected. As a result of the two assumptions, each link $`ij`$ of the newly added node $`j`$ is assigned a weight as

$$w_{ji}=\frac{k_i}{_{\{i^{}\}}k_i^{}}.$$
(2)
Weighted Exponential (WE) Model : The model is inspired by model A discussed in Refs. , and is defined as follows : at every time step we add a new node with $`m(m_0)`$ links, connected with equal probability to the nodes present in the system. The weights of the links are assigned again by using (2).

The difference between the WSF and WE models comes in preferential attachment, which is known to fundamentally alter the topology : The WSF model generates a scale-free network whose connectivity distribution follows $`P(k)k^3`$, while the network generated by the WE model is exponential with the connectivity distribution following $`P(k)=\frac{e}{m}e^{k/m}`$. Since the weights of the links are driven by the connectivity, this difference is expected to lead to significant changes in the distribution of the link strengths as well.

We start by investigating the weight distribution of the two models. As Fig 1 (a) and (b) shows, both the WE and the WSF models lead to a peaked and skewed weight distribution, whose tails decay exponentially (or faster) for large $`w_{ij}`$. The boundedness of $`P(w_{ij})`$ is due to the normalization condition, which does not allow individual weights to be larger than 1. Most important, however, we find that the distribution is stationary, i.e. $`P(w_{ij})`$ is independent of time (and system size).

While the individual weights assigned to links, $`w_{ij}`$, are bounded, we get a very different picture when we study the total weight associated with a selected node. In binary networks the node’s importance is characterized by the total number of links it has, $`k_i`$. Similarly, in a weighted network the importance of a node $`i`$ can be measured by its total weight, obtained by summing the weights of the links that connect to it, $`w_i=_{\{j\}}w_{ij}`$.

Due to the normalization condition (2) a new node has $`w_i=1`$, but $`w_i`$ increases in time every time when a subsequently added nodes link to $`i`$. Since in both models the weights are determined by the network connectivity, we expect that $`P(w)`$ closely follows $`P(k)`$. In contrast, the numerical results summarized in Fig. 2 indicate striking differences between $`P(k)`$ and $`P(w)`$. As Fig 2a shows, while for the WE model $`P(k)`$ decays exponentially, $`P(w)`$ systematically deviates from a simple exponential behavior. This difference is even more evident in the network dynamics: while both $`k_i(t)`$ and $`w_i(t)`$ appear to increase logarithmically in time, they can be fitted with a different slope on a $`\mathrm{log}`$-linear plot (Fig. 2b). Similar systematic discrepancies are observed for the WSF model as well: as Fig. 2c indicates, while $`P(w)`$ can be fitted by a power law, $`P(w)w^\sigma `$, it appears that the exponent $`\sigma `$ is different from $`\gamma =3`$. Furthermore, we find that $`\sigma `$ depends strongly on $`m`$ (Fig. 2c). Again, this difference is reflected in the dynamical behavior of $`k_i(t)`$ and $`w_i(t)`$: as Fig. 2e indicates, $`w_i(t)t^\beta `$ with $`\beta >1/2`$, in contrast with $`k_i(t)t^{1/2}`$ predicted by the binary scale-free model.

To understand the different behaviors of $`w_i`$ and $`k_i`$ uncovered by the numerical simulations, we resort to analytical method in determining the averaged behavior of $`w_i(t)`$ for the discussed model. To simplify the discussion in the following we assume $`m=2`$, however, the calculations can be generalized for arbitrary $`m`$. The total weight of node $`i`$ at time $`t`$ can be written as

$`w_i(t)`$ $`=`$ $`1+{\displaystyle \underset{\{j\}}{}}w_{ij}=1+{\displaystyle _{t_i^0}^t}\stackrel{~}{P}_i(t^{})w_{ij}(t^{})𝑑t^{},`$ (3)
where $`\stackrel{~}{P}_i(t)`$ is the probability that node $`i`$ is selected to be connected to a new node $`j`$ at time $`t`$ and $`t_i^0`$ is the time at which the node $`i`$ has been added to the system. $`w_{ij}`$ is the average weight of link $`ij`$ once the link is established. When a new node $`j`$ and the list of $`m`$ nodes $`\{i^{}\}`$ to which it connects are selected, the weights of the links, $`w_{ji^{}}`$ are assigned according to (2). These weights depend on the number of links the selected nodes have, i.e. $`\{k_i^{}\}`$. If we assume that node $`j`$ is connected to nodes $`i`$ and $`l`$ ($`m=2`$), we have

$$w_{ij}(t)=_m^{\mathrm{}}w_{ji}(l)𝒫(k_l)𝑑k_l$$
(4)
where $`w_{ji}(l)`$ is the weight between the $`j`$ and $`i`$ nodes, $`𝒫(k_l)`$ is the probability distribution of $`k_l`$, the total link number of node $`l`$. Substituting (4) into (3), we obtain

$$w_i(t)=1+_{t_i^0}^t_m^{\mathrm{}}\stackrel{~}{P}_i(t^{})w_{ji}(l)𝒫(k_l)𝑑k_l𝑑t^{}.$$
(5)
According to (2) for $`m=2`$, the weight $`w_{ji}(l)`$ is given by

$$w_{ji}(l)=\frac{k_i}{k_i+k_l},$$
(6)
thus Eq.(5) becomes

$$w_i(t)=1+_{t_i^0}^t_m^{\mathrm{}}\stackrel{~}{P}_i(t^{})\frac{k_i}{k_i+k_l}𝒫(k_l)𝑑k_l𝑑t^{}.$$
(7)
Eq. (7) represents a general expression for calculating $`w_i(t)`$ for $`m=2`$. To apply it to the WE and WSF models, we need to calculate explicitly $`\stackrel{~}{P}(t)`$ and $`𝒫(k_l)`$.

WE model: In the WE model the nodes to which a new node connects to are selected uniformly among all existing nodes, thus the probability that node $`i`$ will be picked is independent of this node’s connectivity and is given by

$`\stackrel{~}{P}_i(t)`$ $`={\displaystyle \frac{m}{t+m_0}}.`$ (8)
Similarly, the connectivity distribution and the dynamical behavior of a single node are given by

$`𝒫(k)`$ $`=Ae^{k/m}={\displaystyle \frac{e}{m}}e^{k/m},`$ (9)
$`k_i(t)`$ $`=m\left[\mathrm{ln}(m_0+t1)\mathrm{ln}(m_0+t_i^01)+1\right]`$ (11)
$`=m\left[\mathrm{ln}(at+b)+1\right],`$
where $`a=\frac{1}{m_0+t_i^01}`$, $`b=\frac{m_01}{m_0+t_i^01}`$ and the normalization condition is $`1=_m^{\mathrm{}}𝒫(k)𝑑k`$.

Substituting (9) into (7), we obtain

$`w_i(t)`$ $`=`$ $`1+e{\displaystyle _{t_i^0}^t}{\displaystyle _m^{\mathrm{}}}{\displaystyle \frac{1}{t^{}+m_0}}{\displaystyle \frac{k_i(t^{})}{k_i(t^{})+k_l}}e^{k_l/m}𝑑k_l𝑑t^{}.`$ (12)
After performing the integration and inserting $`k_i(t)`$ from (9), for large $`t`$ we obtain

$$w_i(t)m\mathrm{ln}(at+b)m\mathrm{ln}(\mathrm{ln}(at+b)+2)+C,$$
(13)
where $`C`$ is an integration constant independent of $`t`$. Therefore the relation between $`w_i(t)`$ and $`k_i(t)`$ for large $`t`$ follows

$$w_i(t)k_i(t)m\mathrm{ln}\mathrm{ln}t+C.$$
(14)
The prediction (14) is fully supported by numerical simulations: in Fig. 3a we plot the difference $`w_i(t)k_i(t)`$ as function of $`\mathrm{ln}\mathrm{ln}(t)`$, showing that the difference indeed follows a double logarithmic law. This result is very interesting since it indicates that the different slopes observed in Fig. 2b for $`k_i(t)`$ and $`w_i(t)`$ do not represent distinct power law scaling behaviors, but are the result of logarithmic corrections.

WSF model: In the scale-free model the probability distributions and $`k_i(t)`$ are given by

$`\stackrel{~}{P}_i(t)`$ $`=m{\displaystyle \frac{k_i(t)}{_j^tk_j}}=m{\displaystyle \frac{k_i(t)}{2mt}}={\displaystyle \frac{k_i(t)}{2t}},`$ (15)
$`𝒫(k)`$ $`=mk^2(kP(k)),`$ (16)
$`k_i(t)`$ $`={\displaystyle \frac{m}{\sqrt{t_i^0}}}\sqrt{t}.`$ (17)
Substituting (15) into (7), and performing the integrals we obtain

$`w_i(t)`$ $``$ $`k_i(t){\displaystyle \frac{m}{8}}\left(\mathrm{ln}{\displaystyle \frac{m^2t}{t_i^0}}\right)^2+{\displaystyle \frac{m}{2}}\mathrm{ln}m\mathrm{ln}{\displaystyle \frac{t}{t_i^0}}+C^{},`$ (18)
indicating that despite a different scaling behavior suggested by the numerical simulations (Fig. 2e), we are dealing with strong logarithmic corrections and asymptotically we have $`\beta ^{}=\beta `$. Again, the analytical prediction (18) is confirmed by more detailed numerical simulations shown in Fig. 3b.

Our ability to calculate analytically $`w_{ij}`$ for the discussed models is based on the fact that the weights are driven by the connectivity distribution. To address the generality of our results we investigated several extensions of these two models, that we discuss in the following.

Weight driven weight- In general one could expect that in some systems the quantity determining the weight is not the connectivity, but are the weights themselves. To investigate this possibility we replaced (2) with

$$w_{ji}=\frac{w_i}{_i^{}w_i^{}}.$$
(19)
i.e. the weight of the newly added links are determined by the total weight of the nodes. While we cannot solve this model analytically, the numerical results are similar to those observed for the WE and WSF models: an apparently different scaling behavior for $`k`$ and $`w`$ can be attributed to slow corrections to scaling.

Weight driven connectivity- In some systems the topology could be driven by the total weights, and not by the connectivity. Thus we assume that the probability (1) that a new node is connected to a node $`j`$ is

$`\mathrm{\Pi }_i={\displaystyle \frac{w_i}{_jw_j}},`$ (20)
where $`w_i`$ is the weight of node $`i`$. The weights are then assigned following (2). We find that the scaling of this network is identical to that of the scale-free model, and the evolution of the weights also follows the paradigm established for the WSF model.

Discussion- Weighted links are a common feature of real networks, thus addressing their scaling behavior is of primarily importance if we are to understand complex networks in general. Here we take a first step in this direction by investigating the scaling properties of several simple models that incorporate mechanisms to assign weights to the links. The weight distribution in the discussed models can be determined analytically and by using numerical simulations, allowing for a thorough analysis of the model’s scaling properties. A comparison between the analytical and numerical results brings deeper understanding of the observed scaling behavior for the distribution functions. First, extensive simulations of networks whose size is comparable to the real networks that are currently available indicate the emergence of new scaling exponents for the behavior of the total weights. However, the analytical solutions reveal that the results are affected by strong logarithmic corrections, and asymptotically the scaling behaviors of the weighted and unweighted models are identical. This result raises important questions regarding our ability to uncover the correct scaling behavior of real weighted networks, should such data become available in the near future: the real exponents could be easily shadowed by corrections to scaling similar to that encountered in the investigated models here.

An important feature of the studied models is the fact that the distribution of the individual weights, $`w_{ij}`$, is bounded and stationary. This implies that the scaling behavior of the total weight distribution can not be explained by the distribution of individual weights. Indeed it is the correlation between the total weight and the linkage probability of the new node that drives the distribution towards power law behavior. In general, one can imagine systems and models where the link distribution, $`P(w_{ij})`$, could also have nontrivial distribution, such as a power law. Investigating such models could be of major future interest.

The results presented in this paper represent only the starting point towards understanding weighted networks. In some real systems, diverse dynamical rules can govern the assignment of weights to links, which could result in statistical properties of the network that are different from that discussed here. In particular, we assumed that once a weight has been assigned to a link, it stays unchanged, which is often not the case in more realistic networks: weights can evolve dynamically just as the network topology does. For example, acquaintance can turn into friendship by strengthening a previously weak link. Determining the generic behavior of such complex evolving systems is a real challenge for future research. Despite these limitations, the investigated models give a glimpse into the complex behavior we are facing as we attempt to make network modeling more realistic by incorporating weights. Finally, the ultimate understanding of weighted networks will be determined by the available data on real systems. While currently such data is rare, we believe that the increasing interest in network modeling and creative data collection methods will soon lead to the development of such data sets, offering further guidance for modeling these complex systems.

Research at Notre Dame supported by NSF, PHY-9988674 and CAREER DMR97-01998.