The Astrophysical Journal, 511:L53–L56, 1999 January 20 q 1999. The American Astronomical Society. All rights reserved. Printed in U.S.A.
ASYMPTOTIC PARTICLE SPECTRA AND PLASMA FLOWS AT STRONG SHOCKS M. A. Malkov Max-Planck-Institut fu¨r Kernphysik, D-69029, Heidelberg, Germany;
[email protected] Received 1998 May 5; accepted 1998 November 18; published 1998 December 15
ABSTRACT In a test-particle (linear) acceleration theory, the shock compression controls the particle spectrum. We present a nonlinear self-similar solution in which the accelerated particles change the flow structure near the shock so strongly that the total shock compression r may become very high. Despite this, the energy spectrum remains close to E23/2 independent of r. This result is valid for the particle diffusivity k(p) growing with momentum faster than p1/2, which is believed to be the case. Subject headings: acceleration of particles — cosmic rays — diffusion — hydrodynamics — shock waves — supernovae: general 1. INTRODUCTION
2. EXACT SOLUTION OF DIFFUSION-CONVECTION EQUATION
The momentum distribution of particles accelerated in shocks by the Fermi process is believed to be a power law f ∝ p2q, with the index q controlled by the shock compression r according to the relation q 5 3r/(r 2 1) (Krymsky 1977; Axford, Leer, & Skadron 1977; Bell 1978; Blandford & Ostriker 1978). This is certainly true if the back-reaction of accelerated particles on the shock structure is negligible and the shock transition is infinitely thin. If, however, the shock accelerates efficiently in the sense of its energy conversion into energetic particles, the accelerated particles (cosmic rays) tend to smear out the shock itself and the above formula is no longer valid. The broadening of the shock transition softens the spectrum (see, e.g., Drury 1983). On the other hand, particles accelerated to ultrarelativistic energies drive the specific heat ratio down to its limiting value, g 5 4/3. This and (particularly) the escape of cosmic rays (CRs) increase the total shock compression, thus hardening the spectrum. The major difficulty, however, is caused by the momentum-dependent particle diffusivity. It allows only the particles with highest momenta to sample the total flow compression. Therefore, the above formula for the index q becomes meaningless, since the majority of particles simply do not know what the total compression is. None of the above effects is small in coupled equations describing the kinetics of the CR gas and hydrodynamics of the thermal plasma in the efficient acceleration regime. Under these circumstances, an exact solution seems to be the only reliable way to find the net result. In view of a significant age of the problem, a complete exact solution is hardly possible. What is possible, however, is to find a solution that tends to the exact one if system parameters tend to their extremes, remaining physically quite realistic. That is, if the maximum energy is very high and, hence, the shock transition is very broad, the solution well inside the shock transition approaches a self-similar form that can be found exactly. In the next section, we first find this solution in a special, one-parameter flow profile. Then, in § 3 we determine the value of this parameter, balancing the CR pressure with the shock ram pressure and thus making our solution self-consistent. This balance tends to an exact one deep inside the shock transition and is approximate at its periphery.
FOR ARBITRARY k(p)
The standard formulation of the problem includes the diffusion-convection equation constrained by conservation of the fluxes of mass and momentum (see, e.g., Drury 1983). We assume that a strong CR-modified shock propagates in the positive x-direction. In the shock frame, the steady mass flow profile is defined as U(x) 5 2u(x), x ≥ 0 and U(x) 5 2u 2 , x ! 0, where u 2 1 0 is the (constant) downstream mass velocity, u(01) 5 u 0 1 u 2, and u(`) 5 u1 1 u 0. The equations read g 1 du g ug 1 k( p) 5 p , x x 3 dx p
[
]
ru 5 r1 u1 , Pc 1 ru 2 5 r1 u12,
(1) (2)
x 1 0.
(3)
Here the number density of CRs is normalized to 4pgdp/p, the particle momentum p is normalized to mc, r(x) is the mass density, r1 5 r(`), and Pc is the CR pressure Pc (x) 5
4p 2 mc 3
E
p1
p0
pd p
Îp 2 1 1
g( p, x).
(4)
The upper limit p1 stands for a boundary in the momentum space (cutoff) beyond which particles are assumed to leave the system instantaneously (g { 0 , p 1 p1 ). In the downstream medium, x ! 0, the only bounded solution is g 5 G( p) { g( p, x 5 0). As indicated, equation (3) is written in the region x 1 0, where we have neglected the contribution of the cold gas (i.e., particles with 0 ! p ! p 0), confining our consideration to sufficiently strong shocks with M 2 { r1 u12/gPg1 k (u1 /u 0 )g, where g is the specific heat ratio of the plasma and Pg1 is the gas kinetic pressure at x 5 ` (for a detailed discussion of this approximation, see Malkov 1997a, hereafter M97a). The subshock strength can be obtained from the Rankine-Hugoniot condition for the gas: rs {
L53
u0 g11 5 . u 2 g 2 1 1 2(u1 /u 0 )g11M 22
(5)
L54
ASYMPTOTIC PARTICLE SPECTRA
Introducing the flow potential f, such as u 5 df/dx, we seek the solution of equation (1) in the form
{
g 5 g0 ( p) exp 2
11b f(x) , x 1 0, k( p)
}
(6)
where b( p) { 2(1/3)d ln g0 /d ln p.
(7)
The substitution given by equation (6) is a central one. It was first found in M97a to reduce equations (1)–(3) to one integral equation. It balances exactly all three terms in equation (1) inside the shock transition where they all are of the same order and, as we shall see in the next section, for a physically meaningful self-consistent velocity field u(x). We know of no other solution of this kind with a momentum-dependent k. Considering u as u(f) and substituting equation (6) in equation (1), we obtain du/df 5 lu/f,
p
(8)
(
)
db d ln k 3 5 (1 1 b) 2 b , dp d ln p l
(10)
where f 0 5 f(0) is another constant (which will be determined below). It is straightforward to verify that the following expression is the first integral of the system made up of equations (7) and (9): g0 ( p)k (1 1 b)
5 const.
3
() p p0
b 0 1 1 23/l #113 p lk 0 0
[
E
2l
p
p0
0
k( p )p
E
s0
u(f) 1 z
s1
ds l21 2fs p 2 (s) s e 5 u1. Î1 1 p 2 (s) b(s)
03/l21
dp
]
0
.
(12)
s 5 (1 1 b)/k
4p mc 2 g ( p ), 3 r1 u12 0 0
(16)
and the function p(s) in equation (14) should be determined from equation (15). The most plausible k(p) dependence is believed to be of the Bohm type, k( p) 5 Kp 2 (1 1 p 2 )21/2, i.e., the mean free path of a particle is proportional to its Larmor radius (here K is a reference diffusivity). Then, equation (14) can be rewritten as u(f) 5 u1 2
z K
E[ s0
s1
11
]
1 s l22e2fsds. b(s)
(17)
According to equations (13) and (15), b(s) is a very simple function, taking in the most part of its domain nearly constant and relatively close values, b . l/3 for Ks K 1 and b . 2l/3 for Ks k 1. It varies monotonically between these limiting values, where Ks ∼ 1. Differentiating equation (17) with respect to f, assuming 0 ! l ! 1 , and considering first the region 1/s 0 K f K 1/s 1, we may obviously replace the lower limit by zero and the upper one by infinity. From equation (17), we then obtain du z . df Kfl
g0 ( p) ∝ k2l ( p) and b . (l/3)d ln k/d ln p.
5
As we shall see, the parameter l depends on the scaling of k as well, and the most surprising consequence of this dependence is that the resulting slope of g0(p) is, in fact, independent of k(p). The region p & p 0 cannot be described within the present approach which produces two integration constants: g0(p0) and b0 in the solution of equation (1) given by equation (12). They serve as external parameters provided by the “injection” solution (see Malkov & Vo¨lk 1995; M97a; Malkov 1998).
(15)
and the limits s 0,1 5 s( p 0,1 ). We have also used equation (11), g0 ∝ s l. The parameter z 5 (l/3)hu1 s2l 0 , where the injection rate h is defined as
For p * p 0, more precisely for (k/k 0 )( p/p 0 ) 3/l k 1, the spectral slope is determined merely by k(p) and “forgets” its behavior at p . p 0: (13)
(14)
We have introduced a new variable s in place of p:
(11)
Denoting k 0 { k( p 0 ) and b 0 { b( p 0 ), for g0 we then have g0 ( p) 5 g0 ( p 0 )
What we have obtained so far is a one-parameter (l) family of exact solutions to equation (1) that require a rather special class of the flow profiles u(f). It is by no means guaranteed that any of these solutions satisfy the pressure balance equation (3). Fortunately, the dependence of Pc(f) in equation (3) can be made consistent with the dependence u(f) ∝ fl (eq. [8]) in a fairly large part of the shock transition, provided that the parameter l is chosen properly. To demonstrate this, we substitute equation (6) into equation (4). Using equation (2), equation (3) can be rewritten as
h5
f(x) 5 f2l/(12l) [(1 2 l)u 0 x 1 f 0 ]1/(12l), 0
2l
3. APPROXIMATE SELF-CONSISTENT SOLUTION
(9)
where l is a separation constant. Equation (8) may be readily integrated and yields for the flow potential
l
Vol. 511
E
`
[1 1 1/b(t/f)]t l21e2tdt
0
zG(l) ¯ [1 1 1/b(t/f)], Kfl
(18)
where G is the gamma function and b(t/f) is replaced in the ¯ with t¯ ∼ 1 . As we have last integral by its mean value at t/f ¯ seen already, the function b(t/f) varies slowly and it is close to l/3 for f 1 K and to 2l/3 for f ! K. Therefore, the f dependence of du/df is determined by the factor f2l and is indeed consistent with equation (8), i.e., with u ∝ fl provided that l 5 1/2. Equation (18) becomes invalid for f * 1/s 1, since the lower limit in equation (17) cannot be replaced by zero in
No. 1, 1999
MALKOV
Fig. 1.—Postshock spectrum G 5 g0(p) exp [2s(p)f0] (solid line) and g0(p) (dotted line; coincides with G, except at p ≈ p0) calculated using eq. (12) for p0 5 0.03, p1 5 1000, M 5 9 # 104, and h 5 0.84 and normalized to dp/p. The subshock compression and other parameters involved in this calculation are found to be rs ≈ 2.8, z ≈ 0.0025, b0 ≈ 1.8, and f0 /K ≈ 1.1 # 1024. The dashed lines show test-particle power-law spectra with the indices q 5 3/(r 2 1) corresponding to shocks of compression r 5 rs (qs ≈ 1.7); r 5 4 and r 5 7, respectively.
this case and the function du/df cuts off as (see eq. [17])
[
]
du z 1 . 11 e2s1f. df KÎs 1f b(s 1 )
(19)
The last asymptotic result obviously matches equation (18) at f ∼ 1/s 1, but equation (8) is no longer consistent with it. At the same time, this is a periphery of the shock transition where u(f) exponentially approaches its limit u1 at x * l { k( p1 )/u1. It should be added here that equation (6) must be slightly modified at the outer part of the precursor where u(x) r u1 (see M97a; Malkov 1997b, hereafter M97b). According to equation (10), the flow profile in the internal part of the precursor (u K u1) is simply linear (l 5 1/2): u(x) 5 u 0 1 (u 20 /2f 0 )x. Assuming K K f K s21 or, equiva1 lently, u 0 K u(x) K u1, from equation (18) we get for u(x) the following independent relation: u 5 98pK 22z 2 x. Equating linear parts of u(x), we find the constant f0: f 0 5 (9/49p)(Ku 0 /u1h) 2 s 0. The remaining unknown constant, the spectral index b0 in equation (12), derives from the test-particle formula d ln G/d ln p 5 23/(rs 2 1), approximately valid at p ≈ p 0 (M97a). For the strong shock modification (u1 k u 0), the parameter b 0 * 1 (see numerical example below), and we infer that f 0 (1 1 b 0 )/k 0 . l. Thus, for f0 we finally have f 0 . (3/7)(2p)21/2 u 0 K/u1h. The Rankine-Hugoniot relation for the precursor may be easily obtained by considering equation (17) at f 5 f 0. The main contribution to the integral comes from the lower limit, and after some simple algebra we have u 0 /u1 5 1 2 (8/p)1/4 (hp1 u 0 /u1 )1/2. For the case of interest, u1 k u 0, we have u1 /u 0 . (8/p)1/2 hp1. Note that this result is valid only in the strong shock limit, M k (hp1 ) 4/3, as a more general treatment in M97a shows. It is identifiable with the most efficient (out of three possible) solution studied there. Besides that, particle spectra of less efficient solutions differ from the spectrum in equation (12) and approach the standard
L55
Fig. 2.—Flow profile u(x) calculated using eq. (17) and df/dx 5 u(f) for the same parameters as in Fig. 1. The calculated precursor compression u1/u0 ≈ 2500. The particle distribution g(x, p) is drawn for different p’s and normalized for each p to its value at x 5 0.
test-particle result. However, these low-efficiency solutions can be realized only for sufficiently weak injection (M97a, M97b). To illustrate the asymptotic nature and limitations of the present solution, we choose a (moderate) ratio not quite favorable for this approximation: p1 /p 0 ≈ 3 # 10 4 . Note that for typical supernova shock conditions, one usually expects p1 /p 0 ∼ 10 8. To produce still a strong modification (hp1 k 1), we take the injection parameter h to be fairly (maybe somewhat unrealistically) high but, in order to prevent extremely strong subshock reduction and to keep the above simplified calculation of f0 and u1/u0 approximately correct, we set the Mach number M to a very high value as well. This does not affect the flow and the spectrum in the precursor where u(x) k u 0 . The case of moderate M can be found in M97a and M97b. Figure 1 shows equation (12) for particle spectrum with k(p) given above. We corrected the spectrum at higher momenta p ∼ p1 [i.e., we multiplied it by (1 1 p/p1 )1/2 ] according to the more consistent treatment of this energy range in M97a. As we mentioned before, at the low-energy end the solution liberates itself from the subshock control very rapidly. Already for ( p/p 0 ) 8 k 1 the spectrum becomes universal, G ∝ p21 , independent of rs. The next turn occurs at p ∼ 1, where the spectrum transforms to G ∝ p21/2, again independent of both the total and subshock compression. Figure 2 shows the flow profile u(x) (eq. [17]) along with the particle distribution g(p, x) at different values of p. One sees that the linear behavior of u(x) or, equivalently, Pc(x) in the internal part of the shock transition results from the sum of exponentially decaying partial contributions to the CR pressure coming from different momenta.
4. CONCLUDING RESULTS AND DISCUSSION
The downstream particle spectrum given by equation (11), with b from equation (13) being expressed in terms of kinetic energy E rather than momentum, exhibits a fairly uniform behavior throughout the entire energy range, relativistic and nonrelativistic. In a standard normalization, F(E)dE, this spectrum
L56
ASYMPTOTIC PARTICLE SPECTRA
has the form F ∝ E 23/2
Î
(E 1 1)(7E 2 1 14E 1 8) , (E 1 2) 3
where E is measured in mc2. Clearly, the overall spectral index is close to 1.5 everywhere except the injection energy (if rs ! 4), the cutoff energy (both ignored in the last formula), and the region E ∼ 1. We emphasize that this 1.5 value of the spectral index is not affected by any parameters involved in these calculations and is, at least in this sense, universal. One may ask, then, What does the power-law index depend on? It should be primarily the momentum dependence of the CR diffusivity k(p) that we have specified in our treatment above. To examine this surmise, we rescale k as k 0 5 k a , where a 1 1/2 and the rescaling of all primed variables below is induced by the above transformation of k. In particular, the spectral slope b is now to be replaced by b 0 5 0 (al/3)d ln k/d ln p (eq. [13]). Recalculating du/df in equation (18) with this rescaled spectrum and CR diffusivity k9, we 0 obtain du/df ∝ f1/a212l . Since the flow profile u(f) must still 0 obey equation (8), i.e., u ∝ fl , we deduce that the rescaled l 0 is l 5 1/2a { l/a. Consequently, the spectral slope b remains unchanged, b 0 5 b. We conclude that the spectral universality survives also the rescaling of k(p); in other words, it is insensitive to the spectrum of the underlying MHD turbulence. On the other hand, the
velocity profile does depend on a. From the above analysis we obtain u ∝ x 1/(2a21), which also explains the condition a 1 1/2.1 Of course, there is a small deviation from this scaling at the distances x corresponding to the diffusion length of particles with p ∼ 1, since b depends on f in equation (18) at f ∼ K. Summarizing the last results, when k(p) rescales, so does the flow profile u(x), but the index b remains invariant. It is not difficult to understand why this is so. As usual in the Fermi process, the spectral slope of course depends on the flow compression. But, since the flow is modified, a particle with momentum p, bound diffusively to the shock front, samples not the total compression but only a compression accessible to it. Now, the latter is determined by the relation f(x) ∝ k a (eq. [6]). As we have shown, u(f) ∝ f1/2a. Therefore, the flow compression u/u 2, as seen by this particle, scales as u/u 2 ∝ f1/2a ∝ [k( p)]1/2 . Since this is independent of a, the index b must also be. This work was done within the Sonderforschungsbereich 328 of the Deutsche Forschungsgemeinschaft (DFG). 1 Remarkably, a precisely opposite condition a ! 1/2 is required to produce a steady velocity jump without momentum cutoff and injection but with a secularly broadening CR precursor. This has been shown by Drury (1983). The fact that our strictly stationary solution that is essentially based on energy losses through the upper cutoff and particle injection at p 5 p0 appears immediately beyond a 5 1/2 is perhaps more than a coincidence. Also, the timesaving numerical solutions with a ≤ 1/2 might be inadequate for modeling the more realistic case of a 5 1.
REFERENCES Axford, W. I., Leer, E., & Skadron, G. 1977, Proc. 15th Int. Cosmic-Ray Conf. (Plovdiv), 11, 132 Bell, A. R. 1978, MNRAS, 182, 147 Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 Drury, L. O’C. 1983, Rep. Prog. Phys., 46, 973 Krymsky, G. F. 1977, Dokl. Akad. Nauk SSSR 234, 1306 (English transl. in Soviet Phys.—Dokl., 22, 327 [1977])
Vol. 511
Malkov, M. A. 1997a, ApJ, 485, 638 (M97a) ———. 1997b, ApJ, 491, 584 (M97b) ———. 1998, Phys. Rev. E, 58, 4911 Malkov, M. A., & Vo¨lk, H. J. 1995, A&A, 300, 605
