A&A 452, 763-771 (2006)
DOI: 10.1051/0004-6361:20064885

Two-fluid models of cosmic ray modified radiative shocks

A. Y. Wagner 1 - S. A. E. G. Falle 2 - T. W. Hartquist 1 - J. M. Pittard 1


1 - School of Physics and Astronomy, University of Leeds LS2 9JT, UK
2 - Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK

Received 20 January 2006 / Accepted 5 March 2006

Abstract
Aims. We demonstrate that cosmic rays suppress the thermal overstability of radiative shocks for a wide range of upstream parameters and initial conditions.
Methods. We perform one-dimensional time-dependent two-fluid simulations, introducing cosmic rays into overstable radiative shocks.
Results. When a cosmic ray fluid is introduced, the oscillations of the shock front are damped and a smooth cosmic ray modified shock develops and propagates away from the cold dense layer. A large increase in cosmic ray pressure occurs for lower Mach numbers in radiative cosmic ray modified shocks than in adiabatic cosmic ray modified shocks. The compression ratio is limited to 7 and steady solutions do not exhibit a cold dense layer like those in radiative shocks that are not modified by cosmic rays. If the diffusion length is much larger than the cooling length, the shock is nearly isothermal. The structure of such shocks obtained with the two-fluid code agree well with analytically derived solutions for strictly isothermal flow containing cosmic rays.

Key words: shock waves - ISM: cosmic rays - hydrodynamics - instabilities - ISM: supernova remnants

  
1 Introduction

Models of radiative shocks have been developed for filaments in the Vela (Slavin et al. 2004), Cygnus loop (Raymond et al. 2001), and other SNRs (Blair et al. 2000,1995). The radiative shocks are often assumed to be steady, although they may display time-dependent behaviour. For instance, multiple shock formation occurs when a remnant undergoes a transition from the adiabatic phase to the radiative phase (Falle 1981) and oscillatory thermal instabilities (or overstability) develop when the remnant is in the radiative phase. This overstability, discovered by Langer et al. (1981) in simulations of accretion flows onto white dwarfs, has been extensively examined analytically (Chevalier & Imamura 1982; Mignone 2005; Bertschinger 1986; Blondin & Cioffi 1989) and numerically (Pittard et al. 2005; Mignone 2005; Strickland & Blondin 1995; Sutherland et al. 2003; Gaetz et al. 1988; Imamura et al. 1984; Innes et al. 1987; Walder & Folini 1996) for a wide range of boundary conditions. If the radiative cooling is approximated by a power law in temperature then there is a critical value of the exponent below which the overstability appears. The validity of the use of steady shock models to interpret emission lines from SNRs has been questioned for some time (Benvenuti et al. 1980; Raymond 1984), and some observed line intensities are not reproduced with steady shock models.

At around the same time as the thermal overstability was discovered by Langer et al. (1981), Drury & Völk (1981) and Axford et al. (1982) analytically studied the effect of diffusive particle acceleration on shock structure within a two-fluid description. The solutions to the two-fluid equations show two notable non-linear features. One is the conversion of a large proportion of the upstream ram pressure into particle energy. The other is the possible existence of up to three different downstream states for a given upstream state. The former naturally encouraged the idea that diffusive shock acceleration was a very efficient method of particle acceleration (Drury 1983). The latter prompted further numerical hydrodynamic investigations to test whether all solutions really exist and are time asymptotically stable (Zank et al. 1993; Drury & Falle 1986; Mond & O'C. Drury 1998; Falle & Giddings 1987; Drury et al. 1995). Jones & Ellison (1991) and Ko (1995) have written extensive reviews on the subject. Becker & Kazanas (2001) have analytically classified the different types of steady solutions.

We note that the use of a two-fluid model is likely to give the maximum cosmic ray pressure production efficiency and a corresponding minimum postshock density enhancement. We address this point more completely in Sects. 4 and 5, where mechanisms for transferring energy from the cosmic rays to the thermal gas are mentioned.

A linear stability analysis of a two-fluid medium including radiative cooling has revealed that cosmic rays suppress the growth of small amplitude perturbations (Wagner et al. 2005). To investigate whether they can damp the thermal overstability of radiative shocks, we performed time-dependent hydrodynamic simulations of a two-fluid medium with radiative cooling. Furthermore, we systematically investigate steady cosmic ray modified shock structures in a wide domain of parameter space.

In the next section we give the full set of equations governing the flow and identify the parameters that characterise the flow. In Sect. 3 we describe the effects of cosmic rays on radiative shocks and in Sect. 4 we concentrate on steady solutions for radiative shocks and compare their structures with those of strictly isothermal shocks for which we can find analytic solutions. Finally, we discuss the results and conclude the paper in Sect. 5.

  
2 Equations and parameter space

The one-dimensional equations that describe the cosmic ray fluid immersed in an inviscid radiative gas are:
  
    $\displaystyle \frac{\partial\rho}{\partial{}t}+\frac{\partial\rho{}u}{\partial{}x}=0,$ (1)
    $\displaystyle \frac{\partial\rho{}u}{\partial{}t}+\frac{\partial\rho{}u^2}{\par...
...ac{\partial{}P_{\rm G}}{\partial{}x}+\frac{\partial{}P_{\rm C}}{\partial{}x}=0,$ (2)


 
$\displaystyle \frac{\partial}{\partial{}t}\left(\frac{\rho{u^2}}{2}+\frac{P_{\r...
...rm C}u\right)
-P_{\rm C}\frac{\partial{}u}{\partial{}x}+\rho^2\mathcal{L}(T)=0,$     (3)


 \begin{displaymath}\frac{\partial{}P_{\rm C}}{\partial{}t}+\frac{\partial{}P_{\r...
...tial{}x}-\kappa\frac{\partial{}^2P_{\rm C}}{\partial{}x^2}=0
.
\end{displaymath} (4)

$\rho, P_{\rm G}, T$ and $\gamma_{\rm G}$ are the density, pressure and temperature and adiabatic index of the thermal gas respectively. $P_{\rm C}, \kappa$ and $\gamma_{\rm C}$ are the pressure exerted by the cosmic rays, the diffusion coefficient and the cosmic ray adiabatic index. x and t are the space and time coordinates, and u is the speed of the flow.

$\rho^2\mathcal{L}(T)$, the cooling rate per unit volume, is taken to be

 \begin{displaymath}
\rho^2\mathcal{L}(T)=\rho^2\Lambda{}T^\alpha,
\end{displaymath} (5)

where $\Lambda$ and $\alpha$ are constants. The gas is allowed to cool to the upstream temperature, towards which $\mathcal{L}\rightarrow0$. $\mathcal{L}<0$ if T is less than the upstream temperature.

In our runs, $\gamma_{\rm G}$ and $\gamma_{\rm C}$ have fixed values of 5/3 and 4/3, respectively, and the diffusion coefficient is also assumed to be independent of time and space. Apart from the cooling term, they are the standard two fluid equations (e.g. Drury & Völk 1981).

Equations (1)-(3) are integrated numerically with an upwind second order Eulerian Godunov scheme. Equation (4) was put into a Crank-Nicholson form and $P_{\rm C}$ solved implicitly. Free flowing boundary conditions are applied at either side of the grid. The initial conditions of the time-dependent runs are given in Sect. 3 and the methods to find steady solutions are explained in Sect. 4. The scheme is stable under the Courant condition based on the total sound speed (Drury & Falle 1986)

\begin{displaymath}\Delta{}t<\frac{\Delta{}x}{\sqrt{\rho(\gamma_{\rm C}P_{\rm C}+\gamma_{\rm G}P_{\rm G})}}\;
\end{displaymath} (6)

where $\Delta{}t$ and $\Delta{}x$ are the time step and the cell width, respectively.

The dimensionless parameters that determine the steady solution are: the upstream Mach number with respect thermal sound speed

\begin{displaymath}M_0=\sqrt{\frac{u_0^2\rho_0}{\gamma_{\rm G}P_{{\rm G}0}}};
\end{displaymath} (7)

the ratio of the upstream cosmic ray pressure to the upstream ram pressure

\begin{displaymath}\Phi_0=\frac{P_{{\rm C}0}}{\rho_0u_0^2};
\end{displaymath} (8)

the exponent $\alpha$ in the cooling law; and the ratio of the cosmic ray diffusion length, $l_{\rm d}$, to the cooling length, $l_{\rm c}$,

\begin{displaymath}q=\frac{l_{\rm d}}{l_{\rm c}}\cdot
\end{displaymath} (9)

Here

 \begin{displaymath}
l_{\rm d}=\frac{4\kappa}{u_{\rm s}}\quad\textrm{and}\quad{}l...
...frac{u_{\rm s}\mu{}kT_{\rm s}}{\rho\Lambda_0T^\alpha_{\rm s}},
\end{displaymath} (10)

where $u_{\rm s}$ is the shock speed, $\rho_{\rm s}$ and $T_{\rm s}$ are the immediate post shock density and temperature given by the Rankine-Hugoniot jump conditions. $\mu$ and k are the average mass per gas particle and Boltzmann's constant. Throughout the paper, the subscripts 0 and 1 denote asymptotic upstream and downstream values respectively. In the examples presented in Sect. 3 we note down the ratio $N_0=\gamma_{\rm C}P_{{\rm C}0}/\gamma_{\rm G}P_{{\rm G}0}$ in brackets to highlight the relative pressure contributions by the cosmic ray and gas upstream. In general, $\Phi=P_{\rm C}/\rho_0u_0^2$ and $N=\gamma_{\rm C}P_{\rm C}/\gamma_{\rm G}P_{{\rm G}0}$.

For a wide part of the $M_0, \Phi_0, q, \alpha $ parameter space, we find little diversity in the qualitative nature of the shock structures. In the following section we will concentrate mainly on M0=3, 5, 10 and 15 shocks. A distinction between high and low Mach number adiabatic cosmic ray modified shocks can be drawn at $M_0\approx6$. For $M_0\ga6$ the two-fluid equations (with $\gamma_{\rm G}=5/3$ and $\gamma_{\rm C}=4/3$) allow multiple solutions. The solution with the largest increase in cosmic ray pressure may either be smooth or contain a subshock. This solution exhibits a large increase in cosmic ray pressure that is "saturated'' with respect to changes in the upstream cosmic ray pressure, in the sense that $P_{{\rm C}1}$ is close to $\rho_0u_0^2$ even if $P_{{\rm C}0}$ approaches 0. Such shocks are said to be "efficient'' in accelerating cosmic rays. Below $M_0\la6$, a single solution exists for a given set of M0 and $\Phi _0$. For such solutions $P_{{\rm C}1}\rightarrow0$ as $P_{{\rm C}0}\rightarrow0$. We will see that this distinction does not occur when radiative cooling is included.

We have studied cases for which $-1.5\leq\alpha\leq0.7$ and employed free flowing boundary conditions. Radiative shocks that are not modified by cosmic rays exhibit thermal overstability if $\alpha$ is less than some $\alpha_{{\rm crit}}$ which is sensitive to the upstream Mach number and the postshock equilibrium temperature (Pittard et al. 2005). M0=3, 5, 10, and 15 shocks are unstable to the fundamental mode of the overstability at $\alpha$ below $\alpha_{\rm crit}\approx-0.7, -0.4, -0.1$, and 0.1 respectively.

We take the parameter q to range from 10 to 0.1. We find that the flow is nearly isothermal when $q\ga5$.

  
3 Suppression of the thermal overstability


  \begin{figure}
\par\includegraphics[width=16.3cm,clip]{4885f01.eps}\end{figure} Figure 1: Evolution of a M0=10 shock, $\alpha =0.7$. At t<0 the radiative shock is stable and is not modified by cosmic rays. The cosmic rays are switched on at time t=0 with $\Phi=0.0075\:(N=1)$ and q=10. a) Density profiles. b) Density space-time diagram. Darker shading represents higher density and the scale is logarithmic. c) Cosmic ray pressure profiles. d) Space-time diagram of $\Phi $. Darker shading represents higher cosmic ray pressure and the scale is logarithmic.
Open with DEXTER

We now present time-dependent solutions to Eqs. (1)-(4). The figures display, in the reference frame of the initial shock front, the evolution of a radiative shock that is not modified by cosmic rays to a cosmic ray dominated shock. The flow enters the grid from the left and leaves the grid on the right. The cosmic rays are introduced into the thermal gas at time t=0 with a pressure homogeneous over the whole grid.


  \begin{figure}
\par\includegraphics[width=16.4cm,clip]{4885f02.eps}\end{figure} Figure 2: Evolution of a M0=10 shock, $\alpha =-1$. At t<0 the radiative shock is overstable and is not modified by cosmic rays. The cosmic rays are switched on at time t=0, $\Phi=0.0075\;(N=1)$, and q=2. The overstability is suppressed. a) Density profiles. b) Density space-time diagram. Darker shading represents higher density and the scale is logarithmic. c) Density space-time diagram: close up of the damping of the overstability. d) Cosmic ray pressure profiles. e) Space-time diagram of $\Phi $. Darker shading represents higher cosmic ray pressure and the scale is logarithmic.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=16.4cm,clip]{4885f03.eps}\end{figure} Figure 3: Evolution of a M0=15 shock, $\alpha =-0.3$. At t<0 the radiative shock is overstable and is not modified by cosmic rays. The cosmic rays are switched on at time t=0, $\Phi=0.0075\;(N=2.25)$, and q=0.1. The overstability is suppressed. a) Density profiles. b) Density space-time diagram. Darker shading represents higher density and the scale is logarithmic. c) Density space-time diagram: close up of the damping of the overstability.
Open with DEXTER

Figure 1a shows the density profiles at several stages during the transition from an initially stable $M_0=10, \alpha=0.7$ radiative shock without cosmic rays to a cosmic ray dominated shock. The cosmic rays are switched on at time t=0 with $\Phi=0.0075\:(N=1)$, and q=10. The abscissa is in units of the diffusion length and the times on the plot are given in units of the diffusion time, $t_{\rm d}=l_{\rm d}/u_{\rm s}$. The sharp jump associated with the beginning of the cold dense layer is smoothed out as the shock becomes modified by the particle acceleration and moves outwards due to the additional cosmic ray pressure downstream.


  \begin{figure}
\par\includegraphics[width=16.2cm,clip]{4885f04.eps}\end{figure} Figure 4: Evolution of a M0=5 shock, $\alpha =-0.3$. At t<0 the radiative shock is overstable and is not modified by cosmic rays. The cosmic rays are switched on at time t=0, $\Phi=0.0075\;(N=0.25)$, and q=10. The overstability is suppressed. a) Density profiles. b) Density space-time diagram. Darker shading represents higher density and the scale is logarithmic.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=15.5cm,clip]{4885f05.eps}\end{figure} Figure 5: Evolution of a M0=3 shock, $\alpha =-1.2$. At t<0 the radiative shock is overstable and is not modified by cosmic rays. The cosmic rays are switched on at time t=0, $\Phi=0.00083\;(N=0.01)$, and q=10. The overstability is suppressed. a) Density profiles. b) Density space-time diagram. Darker shading represents higher density and the scale is logarithmic. c) Cosmic ray pressure profiles. d) Space-time diagram of $\Phi $. Darker shading represents higher cosmic ray pressure and the scale is logarithmic.
Open with DEXTER

The density profile evolution can also be seen in the space-time diagram shown in Fig. 1b. The abcissa is the time in units of $t_{\rm d}$ and the ordinate is the grid axis labeled in units of $l_{\rm d}$, and the flow enters the grid from the top. Darker shading corresponds to higher density and the scale is logarithmic. Here we can follow how the transition occurs smoothly in a couple of diffusion timescales as the cosmic ray dominated shock front moves away from the diminishing cold dense layer. The evolution of the initially homogenous cosmic ray pressure distribution is shown in Figs. 1c and d.

We now show how the cosmic ray pressure suppresses the overstability of a radiative shock. Figure 2 shows, in the same manner as Fig. 1, the density and cosmic ray pressure profiles as well as the corresponding space-time diagrams during the transition from an initially overstable radiative shock to a cosmic ray dominated shock. At times t<0, the M0=10 radiative shock with $\alpha =-1$ exhibits strong oscillations of the shock height and strong variability of the post-shock flow. The cosmic rays are switched on at time t=0, with a homogenous pressure distribution of $\Phi=0.0075\;(N=1)$, and q=2 over the extent of the grid. The periodic displacement of the shock front is damped after a couple of oscillations within a timescale of the order of several $t_{\rm d}$. Figure 2c shows a close up of the initial damping of the overstability. The downstream flow is smoothed and the overstability suppressed as the cosmic ray dominated shock develops and moves away from the slowly decaying cold dense layer. The damping of an overstable M0=15 shock with $\alpha =-0.3$ is shown in Fig. 3.

We now look at lower Mach number shocks. Figure 4 demonstrates the stabilisation of an overstable M0=5 shock with $\alpha =-1.5$. A sequence of density profiles and a space-time diagram for the density are shown. The cosmic rays appear at time t=0, q=10 and $\Phi=0.0075\;(N=0.25)$. The periodic variation in shock height and the large disturbances in the cold dense layer are suppressed, just as in the high Mach number cases.

Figures 5a and b show the effect of cosmic rays on an overstable M0=3 shock with $\alpha=-1.2, \Phi=0.00083\;(N=0.01)$, and q=10. Because q is large, the oscillations of the shock front continue for many periods before the shock is sufficiently modified by the cosmic rays. Nonetheless, the perturbations are smoothed.

  
4 Steady cosmic ray modified shocks


  \begin{figure}
\par\includegraphics[width=6.8cm,clip]{4885f06.eps}\hspace*{1.2cm}
\includegraphics[width=7cm,clip]{4885f07.eps}\end{figure} Figure 6: The steady structure of the cosmic ray modified shock front at time $t>58t_{\rm d}$ during the evolution shown in Fig. 1. In the frame of the shock, $M_0=12.4, \Phi_0=0.005\;(N_0=0.005)$ and q=5.8. $\alpha =-1$.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=6.8cm,clip]{4885f08.eps}\hspace*{1.2cm}
\includegraphics[width=7cm,clip]{4885f09.eps}\end{figure} Figure 7: Steady structure of a M0=10 shock, $\Phi_0=0.0075\;(N_0=1), \alpha=-1$, and q=10.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=6.8cm,clip]{4885f10.eps}\hspace*{1.2cm}
\includegraphics[width=7cm,clip]{4885f11.eps}\end{figure} Figure 8: Steady structure of a M0=5 shock, $\Phi_0=0.0003\;(N_0=0.01), \alpha=-1.5,$ and q=10.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=6.8cm,clip]{4885f12.eps}\hspace*{1.2cm}
\includegraphics[width=7cm,clip]{4885f13.eps}\end{figure} Figure 9: Steady structure of a M0=5 shock, $\Phi_0=0.03\;(N_0=1), \alpha=-1.5$, and q=0.01.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=6.8cm,clip]{4885f14.eps}\hspace*{1.2cm}
\includegraphics[width=7cm,clip]{4885f15.eps}\end{figure} Figure 10: The steady structure of the cosmic ray modified shock front at time $t>50t_{\rm d}$ during the evolution shown in Fig. 5. In the frame of the shock, M0=3.7, $\Phi_0=0.0007\;(N_0=0.8), q=3.5$, and $\alpha =-1.2$.
Open with DEXTER

The steady cosmic ray modified shock structures presented here are found in one of two ways. One is to transform into the frame of the cosmic ray modified shock that is advancing from the cold dense layer in the time-dependent calculations described in Sect. 3. The other way is to start with an adiabatic cosmic ray modified shock and let the gas cool. Unless we refer to a time-dependent run from Sect. 3, the steady solutions were found by the latter method. Either method can be used, however, to obtain any of the steady shock structures in this section.

The flow profile for the steady solutions are shown in the frame of the shock. The upstream and downstream gradients in the flow variables vanish. The upstream state is always to the left and the downstream state to the right of the grid. The quantities shown on the top left of the panels are the upstream values of the flow parameters.

At time $t>58t_{\rm d}$ in the time-dependent run shown in Fig. 2, $M_0=12.4, \Phi_0=0.005\;(N_0\approx1)$, and q=5.8 in the reference frame of the cosmic ray modified shock. The steady profiles of the flow parameters are depicted in Fig. 6.

In Fig. 7 we show the steady structure of a M0=10 shock with $\Phi_0=0.0075\;(N_0=1), \alpha=-1$, and q=10. This nearly isothermal shock is smooth in structure and globally stable, even though a radiative M0=10 shock that is not modified by cosmic rays would be overstable at $\alpha =-1$.

In Fig. 8 we show the structure of a steady $M_0=5, \Phi_0=0.0003\;(N_0=0.01), \alpha=-1.5, q=10$ shock. The shock is stable even though a radiative shock that is not modified by cosmic rays with the same upstream conditions would be overstable. The jump in cosmic ray pressure is about 2300. This is remarkable because the jump in cosmic ray pressure in an adiabatic shock with M0=5 and $\Phi_0=0.0003$ is only 25. The large jump is present even for q=0.1. Clearly, radiative cooling extends the possibility of efficient non-linear acceleration to low Mach numbers. Figure 9 depicts a $M_0=5, \Phi_0=0.03\;(N_0=1), \alpha=-1.5, q=0.01$ shock which contains a subshock. Cosmic ray modified radiative shocks contain a subshock when q is sufficiently low.

For the resulting shock in the time-dependent run shown in Fig. 5, $M_0=3.7, \Phi_0=0.0007\;(N_0=0.8)$ and q=3.5. $\alpha =-1.2$. The profiles are shown in Fig. 10. The downstream to upstream cosmic ray pressure ratio is 900. In an adiabatic $M_0=3.7, \Phi_0=0.0007$ shock, the cosmic ray pressure ratio is about 14. In our search of parameter space, we have not found an example of an overstable cosmic ray modified shock. Strong and weak shocks appear to be stable down to $\alpha =-1.5$ for all values of q and $\Phi _0$.

In Figs. 6-10, the profiles and the jump conditions, except those of the temperature, resemble adiabatic cosmic ray modified shocks rather than radiative shocks without cosmic rays. Some slight differences are seen. For instance, the radiative M0=10 shock for which results are shown in Fig. 7 has globally a smooth structure, whereas an adiabatic shock with the same upstream conditions contains a small subshock. If q is decreased, a subshock appears. The jump in cosmic ray pressure in radiative cosmic ray modified shocks is always higher than in adiabatic cosmic ray modified shocks.

In the radiative cosmic ray modified shock, the gas is prevented from being compressed to a cold dense layer. The overall compression ratio is limited to 7. Shocks for which $q\ga5$ are nearly isothermal. The gas cools significantly in the long cosmic ray precursor leading to the temperature jump becoming very small and little postshock cooling being needed for the gas to drop back to the upstream temperature. Although a smaller q leads to a larger rise in temperature, the peak temperature of the shock is always smaller than that in the corresponding adiabatic case, which is itself much smaller than that given by the single-fluid Rankine-Hugoniot jump condition.


  \begin{figure}
\par\includegraphics[width=15cm,clip]{4885f16.eps}\end{figure} Figure 11: The jump conditions for isothermal cosmic ray modified shocks. a) Compression ratio as a function of $\Phi _0$ for various M0. b) $\Phi _1$ versus $\Phi _0$ for various M0.
Open with DEXTER


  \begin{figure}
\par\includegraphics[width=14.8cm,clip]{4885f17.eps}\end{figure} Figure 12: In all panels the points are from hydrodynamic simulations in which the flow is nearly isothermal and the solid line shows the ideal isothermal solution obtained by integrating the steady state fluid equations. a) Cosmic ray pressure profile of a M0=10 shock. $q=10, \Phi _0=0.0075,$ and $\alpha =0.7$. b) The density profile. c) Cosmic ray pressure profile of a M0=5 shock. $q=10, \Phi _0=0.0003$, and $\alpha =-1.5$. d) The density profile.
Open with DEXTER

The steady state versions of Eqs. (1) and (2) can be used to study jump conditions for isothermal cosmic ray modified shocks. They are

  
    $\displaystyle \rho_0u_0=\rho_1u_1,$ (11)
    $\displaystyle \rho_0u_0^2+P_{{\rm G}0}+P_{{\rm C}0}=\rho_1u_1^2+P_{{\rm G}1}+P_{{\rm C}1}.$ (12)

The far downstream and upstream states can be related by the energy equation for a flow that is isothermal everywhere and contains cosmic rays:
 
$\displaystyle \frac{1}{2}\rho_0u_0^3+P_{{\rm G}0}u_0\ln{u_0}+\frac{\gamma_{\rm ...
...P_{{\rm G}1}u_1\ln{u_1}+\frac{\gamma_{\rm C}}{\gamma_{\rm C}-1}u_1P_{{\rm C}1},$     (13)

where

 \begin{displaymath}\left.\frac{\partial{}P_{\rm C}}{\partial{}x}\right\vert _0=\left.\frac{\partial{}P_{\rm C}}{\partial{}x}\right\vert _1=0.
\end{displaymath} (14)

We also have the isothermal equation of state, which implies that

 \begin{displaymath}P_{{\rm G}0}u_0=P_{{\rm G}1}u_{1}.
\end{displaymath} (15)

Equations (11)-(13) and (15) form a closed set of algebraic equations in four unknowns. One can combine these to obtain an equation for the compression ratio $y=\rho_1/\rho_0$
 
    $\displaystyle \frac{1}{y^2}\left(\frac{1}{2}-\frac{\gamma_{\rm C}}{\gamma_{\rm ...
...eft(1+\frac{1}{\gamma_{\rm G}M_0^2}+\frac{1}{\gamma_{\rm C}M_{\rm C0}^2}\right)$  
    $\displaystyle ~~~~+\frac{1}{\gamma_{\rm C}M_0^2}\ln{\frac{1}{y}}-\frac{1}{2}-\f...
...ft(\frac{1}{\gamma_{\rm G}M_0^2}+\frac{1}{\gamma_{\rm C}M_{\rm C0}^2}\right)=0.$ (16)

$M_{\rm C}\equiv\sqrt{\rho{}u^2/\gamma_{\rm C}P_{\rm C}}$ is the cosmic ray Mach number. In the strong shock limit ( $M_0\rightarrow{}\infty$, $M_{\rm C0}\rightarrow{}\infty$), the solution to Eq. (16) is

\begin{displaymath}y=\frac{\gamma_{\rm C}+1}{\gamma_{\rm C}-1}=7,\quad{}\gamma_{\rm C}=\frac{4}{3}\cdot
\end{displaymath} (17)

Figure 11a shows y as a function of $\Phi _0$ and M0. In Fig. 11b we have plotted the downstream cosmic ray pressure as a function of the upstream cosmic ray pressure and M0. The compression ratio converges to a value of 1 as $\Phi _0$ approaches unity. The analytic results for the jump conditions agree well with those found in the hydrodynamic simulations for radiative cosmic ray modified shocks for large q.

One can solve Eqs. (11), (12) and (15), which are valid throughout an isothermal flow, for u and substitute into the time independent version of Eq. (4). Equation (4) can be numerically integrated to allow a comparison of an isothermal shock profile with a profile from a time asymptotic two-fluid result for a shock in which the flow is not exactly isothermal. Comparisons are shown in Fig. 12 for M0=10, $\Phi_0=0.0075\;(N_0=1),q=10$, and $\alpha =0.7$ as well as $M_0=5, \Phi_0=0.0003\;(N_0=1), q=10$, and $\alpha =-1.5$.

In general, the downstream jump conditions of radiative cosmic ray modified shocks depend on the upstream Mach number and cosmic ray pressure. We found no strong dependence on $\alpha$ or q. A large gain in cosmic ray pressure in a radiative cosmic ray modified shock is found for small $\Phi _0$ even at low values of M0, just as found from the algebraic calculations of the isothermal jump conditions (see Fig. 11b). In fact, the upstream cosmic ray pressure in the steady solutions does not influence the qualitative structure and stability of high Mach number shocks. The quantitative influence is also small.

The two fluid model for steady cosmic ray modified shocks implies a maximum compression ratio of 7, which corresponds to a cosmic ray dominated flow downstream ( $\gamma_{\rm C}=4/3$). Radiative cooling increases the overall compression towards this limit, and, consequently, efficient solutions exist down to lower Mach numbers than in adiabatic shocks. This is directly related to our finding that the thermal overstability is suppressed even for weak shocks: the overstability is suppressed if the cold dense layer is absent and the presence of cosmic rays act to ensure its absence. Without a cold dense layer that reflects sound waves effectively (Pittard et al. 2005), perturbations are not confined within the cooling region and are damped by the large post-shock cosmic ray pressure (Wagner et al. 2005).

In adiabatic cosmic ray modified shocks, multiple solutions may occur when $M_0\ga6$ and the upstream cosmic ray pressure is sufficiently small. If three steady solutions exist for one set of upstream conditions, up to two of them may be corrugationally stable and stable against spontaneous accoustic emission (Mond & O'C. Drury 1998). If we start with these stable solutions in our code and let the gas cool, only one solution readjusts itself to a stable radiative shock. We have not encountered any hint that multiple solutions exist in cosmic ray modified radiative shocks. In the isothermal limit they cannot exist.

Even though the two fluid model allows cosmic rays to be accelerated from small upstream pressures to a "saturated'' value far downstream, the result should be viewed with caution. The rise in energy density in these "efficient'' shocks is due to a diminishing number of particles being trapped for a long time in the vicinity of the shock front and perpetually gaining energy (Drury 1983). Fully kinematic calculations in which $\kappa$ has a power law dependence on particle momentum, $\kappa\propto{}p^\beta$, demonstrate this (Duffy 1992; Falle & Giddings 1987). Drury et al. (1995) explained that if $\beta<1$ in plane parallel shocks (or $\beta<1/3$ in spherical symmetry) the particle energy density can keep increasing with acceleration time to a value at which the increase is balanced by downstream advection. Unless there is an additional loss mechanism for particles with high momenta, the cosmic ray pressure enhancement can be very large. The existence of "efficient'' solutions in cosmic ray modified radiative shocks may be explained in the same manner.

On the level of a fluid description of cosmic rays, the strong rise in cosmic ray pressure may be moderated by invoking a process that transfers energy from the cosmic ray fluid to the gas. The model of Völk et al. (1984) includes the heating of the thermal gas by the dissipation of Alfvén waves generated by cosmic ray streaming. The work done by the cosmic rays is $V(\partial{}P_{\rm C}/\partial{}x)$, where V is the Alfvén speed, and the energy goes directly into the gas. This process, however, is not efficient if $u\gg{}V$. Although the two-fluid model we have worked with here does not contain a mechanism to limit the cosmic ray energy and will therefore always yield steady solutions with maximum acceleration efficiencies, our work serves as the first step in the study of the cosmic ray modification of radiative shocks within a wide parameter range.

  
5 Conclusion

The principal result of our study of cosmic ray modified radiative shocks is that there is no region of $M_0, \alpha, \Phi$, and q parameter space that we explored in which the inclusion of cosmic rays does not suppress the radiative thermal overstability. We have demonstrated the damping of the thermal overstabilty in time-dependent calculations in which the initial state corresponds to an overstable radiative shock without cosmic rays, and subsequently introduce the cosmic rays with a homogeneous pressure over the whole grid.

The steady shock structures shown in this paper can all be reached from initial states in which the shocks are unaffected by cosmic rays but experience the subsequent homogeneous injection of cosmic rays. They can also be reached if one starts with a steady adiabatic cosmic ray modified shock structure with the desired upstream conditions and spontaneously switches on the cooling of the thermal gas. In general, the steady shocks resemble cosmic ray modified adiabatic shocks more than one-fluid radiative shocks. In particular a maximum compression ratio of 7 and an absence of a cold dense layer are seen. For $q\ga5$ a shock is nearly isothermal and the jump conditions and flow profile can be calculated from the steady state equations. We have also found that the strong rise in cosmic ray pressure can occur at low Mach numbers, whereas in adiabatic cosmic ray modified shocks it occurs only for $M_0\ga6$. This is a direct consequence of radiative cooling in the two-fluid model. Radiative cooling enhances compression but the compression ratio cannot exceed 7 as the downstream flow becomes cosmic ray dominated.

In the absence of a cold dense layer there is no sharp boundary in the postshock flow with high reflectivity that would contain the sound wave energy of the overstable oscillations inside the cooling region. Instead, any perturbations compressing the downstream flow quickly lead to a conversion of the sound wave energy to an increase in the internal energy of the cosmic rays. The perturbations are quickly damped as the cosmic rays carry away this energy as they escape downstream.

The results presented here encourage one to re-address the assumptions on the nature of shocks when performing SNR shock diagnostics. The steady shock structures reported here are so radically different from structures used in diagnostic studies that they may not be compatible with observations. Disagreements between observational and model spectra obtained using cosmic ray modified shocks may be avoided if the shocks are younger than the cosmic ray diffusion time. Another possibility is that an instability, such as that studied by Drury & Falle (1986), leads to a conversion of cosmic ray energy to thermal energy of the non-cosmic ray fluid.

References

 

Copyright ESO 2006