Diffusive shock acceleration (Fermi acceleration) of cosmic rays in a system containing two shock waves is investi- gated using test particle simulations and analysis of the diffusion convection equation. We assume that the cosmic ray acceleration timescale is much less than the approaching timescale of the two shocks. Low-energy cosmic rays are primarily accelerated as they interact with one of the two shocks. However, as they are energized, and their mean free path becomes comparable to the separation distance of the two shocks, they start to accelerate as they interact with both shocks. As a result, a double power-law-type spectrum for the cosmic ray distribution is produced. The numerical result is well explained by a model based on the diffusion convection equation, which accounts for the dependence of the diffusion coefficient on the cosmic ray energy. The break point of the spectrum can be estimated by consider - ing transport of cosmic rays from one shock to the other. Keywords: Cosmic rays, Diffusive shock acceleration, Multiple shocks authors found that ions were accelerated by first-order Introduction Fermi acceleration through the two shocks, and their One of the most widely accepted mechanisms for pro- energy was higher than expected from theory based on ducing cosmic rays (energetic particles) is diffusive shock the single shock DSA model. When there are two shocks acceleration (DSA). This is also known as first-order located in close proximity, the DSA model using a single Fermi acceleration, in which the cosmic rays gain energy shock may not explain the acceleration process. as they are scattered by turbulence upstream and down- Particle acceleration in multiple shocks has been stream of the shock, repeatedly cross the shock, and treated theoretically by Melrose and Pope (1993). They are effectively compressed by shock flow Blandford and considered a system where “seed” particles initially given Eichler (1987). Although the DSA of the cosmic rays has at the location of the shock waves are continuously accel- been studied extensively, most past studies have exam- erated by a sequence of identical shock waves. Further- ined the case wherein cosmic rays interact with a single more, they concluded that the energy spectrum becomes shock wave. The DSA in a system with multiple shocks is −1 n(p) ∝ p for infinite number of shocks, where n(p) not well understood. is the density distribution function. The spectrum for Shock waves are ubiquitous in space. For instance, dual shocks is harder than that of a single shock whose shock waves can be generated in front of coronal mass compression ratio is 4. Recently, Zhao and Li (2014) ejections, in association with the co-rotational interac- performed test particle simulations for a pair of paral- tion regions, and in front of planetary magnetospheres, lel shocks approaching each other. They included mag - including that of the earth. These shocks frequently netic turbulence and calculated the equation of motion approach and collide with each other. The collision and pitch-angle diffusion of particles. Their calculations of an interplanetary shock with the earth’s bow shock showed that the power-law spectrum index for a pair of was reported and analyzed in Hietala et al. (2011). The shocks is harder than that of a single shock, and the max- imum energy is greater than that of a single shock. *Correspondence: firstname.lastname@example.org In the present study, we account for accelerated parti- Interdisciplinary Graduate School of Engineering Sciences, Kyushu cles from one shock traveling to the other shock against University, 6-1 Kasuga-Koen, Kasuga, Fukuoka 816-8580, Japan Full list of author information is available at the end of the article a plasma flow, an effect that is not considered in Melrose © The Author(s) 2018. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Nakanotani et al. Earth, Planets and Space (2018) 70:33 Page 2 of 6 and Pope (1993). We discuss the particle acceleration in and “Region 3” (L < x), which correspond to upstream of two shock waves by performing test particle simulations the two shocks, downstream of Shock 1 and upstream of and evaluating the convection diffusion equation. Shock 2, and downstream of both shocks, respectively. This paper is organized as follows. In “Introduction ” Next, p and p are the cosmic ray momentum vectors section, we describe our model and assumptions used defined in the shock rest and plasma rest frames, respec - in our test particle simulation and theoretical analy- tively. These are related through the Lorentz transforma - sis. The results of the test particle simulations are pre - tion using the local plasma flow velocity, u(x)e , where sented in “Test particle simulation” section. To explain u(x) = u (= const .) for the three regions, i = 1, 2, 3, and the simulation results, we develop a theoretical analysis e is the unit vector in the x direction. In the present in “Theoretical analysis” section. In “Discussion ” section, model, the cosmic ray travels freely with a constant speed we discuss the breaking points of the double power-law- along a straight line for the time period, τ, until it collides type spectrum by examining the dependence of the break with a scatterer, i.e., MHD turbulence, which is assumed point energy on other parameters. Finally, a summary is to move with the local plasma flow. provided in “Summary” section. Because cosmic ray energy is conserved in the plasma rest frame, the cosmic ray momentum after the collision ′ ′ ′ ′ Test particle simulation can be written as, p = p e , where p =|p | is the magni- Let us consider a system in which there are two hydro- tude of the momentum in the plasma rest frame, and e is dynamic planar shocks, “Shock 1” and “Shock 2,” located a unit vector representing the momentum direction after at x = 0 and x = L in Cartesian coordinates, respec- the collision. tively (Fig. 1). For simplicity, we ignore the background For simplicity, the scattering time, τ, is constant. One magnetic field. Because we consider this problem as can specify τ according to some given probability distri- one-dimensional, we do not have to include adiabatic bution function, but within the context of the present decompression, which is treated in Melrose and Pope study, the generalization for τ does not alter the results (1993). The shocks are assumed to be infinitesimally thin, qualitatively, as long as the central limit theorem is sat- and both shock normals are pointing toward the x direc- isfied for the convolution of the probability distribution tion. The assumption of zero scale length of the shock function. transition layer is made due to its simplicity, but it may We also assume that cosmic ray scattering due to tur- be violated when a shock precursor with a cosmic ray bulence is isotropic in the plasma rest frame. (e is chosen scattering length scale is generated in front of the shock randomly in all the directions in three dimensions.) This Drury and Voelk (1981). Here, we simply assume that the corresponds to the assumption of large-angle scattering, energy density of the cosmic rays is small enough. which may represent the scattering due to a strongly tur- We further assume that the cosmic ray speed is much bulent magnetic field, or an approximation of successive faster than the shock speed, and that the timescale of par- small-angle scattering, such as the pitch-angle scattering ticle acceleration is sufficiently shorter than the timescale by small-amplitude waves. The difference between the for the two shocks to approach each other. Specifically, large-angle scattering and the pitch-angle scattering mod- we assume that shocks are at rest and time stationary. els is discussed in Naito and Takahara (1995). We treat the The two shocks divide the entire space into three regions motion of particles in electromagnetic fields as random. as in Fig. 1: “Region 1” (x < 0), “Region 2” (0 < x < L), First, we inject seed particles at Shock 1 and Shock 2. The initial particle momentums are all equal to p , but their directions are random. They then travel for τ and are scattered by turbulence, and a new momentum is computed after collision; they then travel again, and the cycle continues. After repeating this cycle for a certain length of time, we collect particle data and compute the cosmic ray spectrum. The thick solid line in the large panel shown in Fig. 2 shows the result from the test particle simulation. Test particles are injected at both Shock 1 and Shock 2. Parameters used for this run are: u = 4, u = 2, 1 2 u = 1 (normalized to u ), L = 125 (normalized to τu ), 3 3 3 p =|p |= 7.5 (normalized to m u , where m is the rest 0 0 0 3 0 mass of the particle), and the number of particles used is Fig. 1 Schematic showing shock positions, regions, and flow speeds 2 × 10 . Nakanotani et al. Earth, Planets and Space (2018) 70:33 Page 3 of 6 The calculation time is (normalized to tau). 2 × 10 τ The horizontal axis is the logarithm of the cosmic ray momentum in the shock rest frame, p, and the verti- -4 cal axis is the logarithm of the cosmic ray density dis- tribution function, n(p). Specifically, n(x, p)�p�x corresponds to the number of particles with its position within the range (x, x + �x) and its momentum within 4 -4 -2 the range (p, p + �p). In using these ranges, we ensure p that the numerically obtained cosmic ray distribution is 4 nearly isotropic. We collect all the particles in Region 3 to evaluate the distribution function. Because Region 3 cor- 0 0.5 1 1.5 2 responds to the downstream of Shock 2 , n(x, p) must be 0 0.5 1 1.5 2 independent of x. log (p/p ) 10 0 The thick solid line in the small panel in Fig. 2 shows Fig. 2 Large panel: comparison of density distribution obtained from the result of the test particle simulation for one single the test particle simulation for two shocks (thick solid line), theoreti- cal model based on the diffusion convection equation with constant shock. The conditions are the same with the above simu - diffusion coefficient (dotted line), and the same model that takes into lation except for omitting Shock 2 (u = 2) and injecting account the dependence of the diffusion coefficient on the cosmic particles only from Shock 1. It is well known that the DSA ray energy (thin solid line). The two blue lines with indices − 4 and −Ŵ process produces the power-law spectrum, n(p) ∼ p , − 2 are indicated for reference. Test particles are injected at both with Ŵ = (r + 2)/(r − 1), where r is the shock compres- “Shock 1” and “Shock 2.” Small panel: density distribution obtained by the test particle simulation for one single shock sion ratio. For Shock 1, r = u /u = 2 and Ŵ = 4. We 1 2 confirm that the spectrum obtained from the test particle simulation corresponds to the index from the theoretical (i = 1, 2, 3), with u and k constant. Then, the time sta - analysis. i i tionary solution of (1) for each region can be written as As shown in the large panel, the spectrum in the low- energy regime (0 < log(p/p )< 0.78) and the high- C (p) u x i i f (x, p) = + B (p) exp , energy regime (0.78 < log(p/p )< 2.5) is approximated 0 i i (3) u k i i by power-law distributions with different indices. The where C (p) and B (p) are arbitrary functions of p. From low-energy part is comprised of cosmic rays accelerated i i the boundary conditions at x → ±∞ and continuity of by one of the two shocks, and the high-energy part is f (x, p) across the two shocks we find comprised of particles that travel long enough distances i to interact with both shocks. For Shock 1 and Shock 2, C (p) C (p) = B (p) = 0, B (p) = + B (p) , 1 3 1 2 r = u /u = u /u = 2 and Ŵ = 4; for the two shocks 1 2 2 3 together the compression ratio is r = u /u = 4, which 1 3 (4) C (p) C (p) u L 3 2 2 gives Ŵ = 2. These values explain the slope of the spec - = + B (p)exp . u u k 3 2 2 trum for the low- and high-energy regimes in Fig. 2. Integration of (1) across the first shock at x = 0 yields Theoretical analysis ∂f p ∂f −(k + u) = 0, In this section we describe a theoretical model based on (5) ∂x 3 ∂p the diffusion convection equation (e.g., Skilling 1975). If where the bracket denotes the difference across the we assume that the cosmic ray distribution is isotropic, shock. Substituting (3) into (5), we have its time evolution can be described by p dB u C p d C 2 1 2 2 ∂f ∂f p du ∂f ∂ ∂f B + + + + u − = k + Q(x, p, t), (1) 3 dp u − u u 3 dp u 1 2 2 2 ∂t ∂x 3 dx ∂p ∂x ∂x (6) where f(x, p, t) is the cosmic ray distribution func- = δ(p − p ). u − u 1 2 tion defined in the shock rest frame, p =|p|, u(x) is the plasma velocity, k(x) is the spatial diffusion coefficient, Similarly, across the second shock at x = L, we obtain and Q(x, p, t) is a source term, u B p dB pB dξ p d C 2 2 2 2 2 + + e + Q(x) = A δ(x)δ(p − p ) + A δ(x − L)δ(p − p ), 1 0 2 0 (2) u − u 3 dp 3 dp 3 dp u 2 3 2 where A and A are constant. Within the three regions 1 2 A (7) = δ(p − p ), we defined before, we let u(x) = u and k(x) = k i i u − u 2 3 log n 10 Nakanotani et al. Earth, Planets and Space (2018) 70:33 Page 4 of 6 where ξ = u L/k (p) is a dimensionless parameter repre- After substituting the above into (6) and (9), and seeking 2 2 senting the relative magnitude of the convection and dif- a non-trivial solution for α and β, we find fusion for cosmic rays in Region 2. If ξ >> 1, convection s(s − 3) = (s − s )(s − s )e , (12) 1 2 is more dominant than diffusion, and cosmic rays start - ing from Shock 2 cannot reach Shock 1. Therefore, the where s = 3u /(u − u ) and s = 3u /(u − u ) are the 1 1 1 2 2 2 2 3 two shocks are recognized as two well-separated shocks power-law indices of the density distribution function for particles seeded from Shock 2. On the other hand, if due to the DSA by Shock 1 and Shock 2, respectively. ξ << 1, diffusion is more dominant than convection, and In this calculation, we have neglected the second term the two shocks can essentially be viewed as a single shock in (9) for simplicity. When ξ is very large, (12) gives two with a velocity jump from u to u . roots, s = s and s = s , suggesting that the two shocks 1 3 1 2 The value of ξ depends on the momentum of the parti- are independent. The DSA calculated with the two shocks cle through the spatial diffusion coefficient, separately produces distributions specified by their com - pression ratios. In contrast, in the limit of ξ → 0, (12) 2 2 2 2 (v τ) v τ τc p becomes a linear equation for s, and the solution is k =< >= = , 2 (8) 2 2 2τ 6 6(m c + p ) 0 s s 3u 1 2 1 s = = . (13) s + s − 3 u − u 1 2 1 3 where <> denotes the ensemble average, and v is the x component of the velocity vector in three dimensions. As expected, the spectrum index from the DSA due to the Then (7) is rewritten as two shocks combined is obtained. Equation (13) implies that the spectrum of the two shocks can be harder than u 2ξ p dB p d C 2 2 2 the strongest single shock. If there are two strong shocks − B + e + u − u 3γ 3 dp 3 dp u 2 3 2 with compression ratios of 4, Eq. (13) becomes 3.2. When considering a series of multiple shocks, the maximum = δ(p − p ), 0 (9) u − u 2 3 spectrum index can be obtained from Eq. (13), apply- ing the upstream flow of the first shock as u and the where γ = 1 + . downstream flow of the last shock as u . It is easily seen 2 2 m c that the maximum index approaches 3 as the number of Equations (6) and (9) are simultaneous differential shocks increases because the downstream flow of the last equations for B (p) and C (p), and the solutions can be 2 2 shock approaches 0, in agreement with Melrose and Pope obtained numerically. Then, the distribution function is (1993). given as Discussion C (p) u L 2 2 2 2 n(p) = 4πp f (p) = 4πp + B (p) exp . 3 2 We have seen that the diffusion convection equation well u k 2 2 explains the test particle simulation results presented (10) in the previous section. In this section, we discuss the For the DSA with a single shock, the cosmic ray spectrum dependence of spectrum breaking using Eqs. (6) and (9). does not depend on the diffusion coefficient, which is not Below we discuss determining the momentum at the the case for the present problem. The dotted line in Fig. 2 breaking point, p , by examining its dependence on shows (10), in which k is given by (8) with p = p . (p is 2 0 0 shock parameters, L, u , u , and u . 1 2 3 the initial momentum.) In this case, k is a constant. All Figure 3a shows p as a function of L (the crosses), other parameters are chosen according to the test parti- keeping other parameters unchanged at u = 4, u = 2 1 2 cle simulation run. and u = 1. This plot is made by numerically evaluat - The result does not agree with the simulation. In par - ing the first derivative of log n(p) with respect to log p. ticular, the model does not have a double power-law dis- We define the break point as the point at which the first tribution produced by the test particle simulation. derivative has the minimum value. Figure 3b suggests The thin solid line in Fig. 2 shows (9) with the p that p becomes larger as L increases. This relationship dependence correctly kept in the definition of k in (8). is due to the increase in L when a particle starting from Now, the result agrees well with the test particle simula- Shock 2 must travel a long distance to reach Shock 1. tion. In particular, the breaking of the power spectrum is Figure 3b shows p as u (circles) varies with param- b 2 included. eters L = 125, u = 4, and u = 1. When u increases, 1 3 2 To seek an asymptotic solution for (6) and (9), let us p also increases, because a particle starting from Shock write 2 must overcome the fast convection of the background −s −s flow to reach Shock 1. C (p) = αp ; B (p) = βp . 2 2 (11) Nakanotani et al. Earth, Planets and Space (2018) 70:33 Page 5 of 6 The break point is not strongly affected by variations The spatial diffusion coefficient approaches for τ c /6 in u and u (not shown here) because u and u do not highly relativistic particles. In this case, is written 1 3 1 3 max influence the particle transport between Shock 1 and as (c/4u )τc. Therefore, considering the time station - Shock 2. ary solution, the condition L << (c/4u )τc should be We also estimate the value of p by comparing the dif- applied. fusion and convection in Region 2. Suppose a particle with momentum p, located at Shock 2, starts a random Summary walk at time t = 0 in a convecting flow with a speed u . In summary, we examined the DSA process in a system The motion of the particle may then be given as a super - including two shocks. We assume that an acceleration position of the random walk and convection, timescale in a shock is much less than the approaching timescale of the two shocks. These two shocks can be (t) = 6k t − u t, (14) 2 2 treated as stationary under this assumption. Although where is the travel distance along the x direction meas- we have to consider finite acceleration time for a more ured from the location of Shock 2 at x = L. Because k realist case, which means two shocks moving at differ - is an increasing function of the particle momentum p, ent speeds, we believe that it does not change the result faster particles can travel longer distances compared with significantly. slower particles. Moreover, due to convection, there is a Test particle simulations are performed, and a double maximum distance a particle can travel, = 3k /2u . power-law-type spectrum is obtained. A simple model max 2 2 We estimate that the momentum at which = L gives based on the diffusion convection equation explains the max p . Namely, particles with momentum larger than p can result well. In the model, it is important to include the b b travel from one shock to the other, potentially accelerat- momentum dependence on the spatial diffusion coeffi - ing from both shocks. From this argument, we estimate cient, unlike the DSA model with a single shock. We discussed the dependence of the break point, p , p δ b on other shock parameters, and showed that the value of = , (15) m u p can be estimated by balancing the convection and dif- 0 3 1 − δ b fusion in a region between the two shocks. Our results where δ = 4(L/τu )(u /u ) and c is the speed of light. indicate that an acceleration in two shocks is more effi - 3 2 3 In Fig. 3 we plot (15) as the solid line for L (a) and the cient than a single shock, which is consistent with previ- dotted line for u (b). This model well explains the func - ous results. tional dependence of p on L and u . This dependence of The DSA in a two-shock system is similar to the DSA b 2 p on L becomes different between simulation and theo - from cosmic ray-mediated shocks, in that the back- retical results for large L because we implicitly assume ground flow velocity difference used to accelerate parti - that the distribution function can spread out over the dis- cles depends on particle energy, which forms the double tance , which a particle can reach when we consider power-law spectrum. Similarly, spectral hardening has max time stationary solution (4). been observed at solar flare termination shocks Li et al. a b 1.4 1.2 0.8 0.6 0 250 500 750 1000 1 2 3 4 Fig. 3 The breaking point p as a function of L (crosses) and u (circles), numerically evaluated using the first derivative of log n(p) with respect to b 2 log p. Solid lines are Eq. (15) log (p /p ) 10 b 0 Nakanotani et al. Earth, Planets and Space (2018) 70:33 Page 6 of 6 References (2013) and Kong et al. (2013). In these previous studies, Blandford R, Eichler D (1987) Particle acceleration at astrophysical shocks: a the finite thickness of the termination shock was consid - theory of cosmic ray origin. Phys Rep 154:1–75 ered, and it was concluded that spectral hardening occurs Drury LO, Voelk JH (1981) Hydromagnetic shock structure in the presence of cosmic rays. Astrophys J 248:344–351 because high-energy electrons can see the full width, Hietala H, Agueda N, AndréEová K, Vainio R, Nylund S, Kilpua EKJ, Koskinen while low-energy particles accelerate for only part of the HEJ (2011) In situ observations of particle acceleration in shock–shock shock width. interaction. J Geophys Res (Space Physics) 116:10105 Kong X, Li G, Chen Y (2013) A statistical study of the spectral hardening of Authors’ contributions continuum emission in solar flares. Astrophys J 774:140 MN designed the test particle simulation algorithm and drafted the manu- Li G, Kong X, Zank G, Chen Y (2013) On the spectral hardening at GSIM 300 kev script. TH and SM helped analyze and interpret the results and draft the in solar flares. Astrophys J 769:22 manuscript. All authors read and approved the final manuscript. Melrose DB, Pope MH (1993) Diffusive shock acceleration by multiple shocks. Proc Astronom Soc Aust 10:222 Author details Naito T, Takahara F (1995) Monte carlo simulation of diffusive particle accelera- Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, tion in shock waves with oblique magnetic fields. Mon Not R Astronom 6-1 Kasuga-Koen, Kasuga, Fukuoka 816-8580, Japan. Faculty of Engineering Soc 275:1077–1092 Sciences, Kyushu University, 6-1 Kasuga-Koen, Kasuga, Fukuoka 816-8580, Skilling J (1975) Cosmic ray streaming. I—effect of alfven waves on particles. Japan. International Center for Space Weather Science and Education, Mon Not R Astronom Soc 172:557–566 Kyushu University, 744 Motooka, Nishi-Ku, Fukuoka 819-0395, Japan. Zhao L, Li G (2014) Particle acceleration at a pair of parallel shocks near the sun. J Geophys Res (Space Physics) 119:6106–6119 Acknowlegements We thank C. Mazelle for carefully reading the manuscript. Competing interests The authors declare that they have no competing interests. Ethics approval and consent to participate Not applicable. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub- lished maps and institutional affiliations. Received: 8 November 2017 Accepted: 6 February 2018
Earth, Planets and Space – Springer Journals
Published: Feb 21, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud