PARTICLE ACCELERATION BY COLLISIONLESS SHOCKS CONTAINING LARGE-SCALE MAGNETIC-FIELD VARIATIONS

Collisionless shocks in space and other astrophysical environments are efficient accelerators of energetic charged particles. Diffusive shock acceleration (DSA; Krymsky 1977; Axford et al. 1977; Bell 1978; Blandford & Ostriker 1978) is the most popular theory for charged-particle acceleration. It naturally predicts a universal power-law distribution f ∝ p^−γ with γ ∼ 4.0 for strong shocks, where f is the phase-space distribution function, close to cosmic rays observed in many different regions of space. The basic conclusions of DSA can be drawn from the Parker transport equation (Parker 1965) by considering the shock to be a compressive discontinuity in an infinite one-dimensional and time steady system. DSA is thought to be the mechanism that accelerates anomalous cosmic rays (ACRs) in the Heliospheric termination shock and also galactic cosmic rays (GCRs) with energy up to at least 10¹⁵ eV in supernova blast waves. However, recent in situ observations in the termination shock and the Heliosheath by Voyager 1 (Stone et al. 2005) found that the intensity of ACRs is not peaked at the termination shock and the energy spectrum is still unfolding after entering the Heliosheath, which strongly indicates that the simple planar shock model is inadequate to interpret the acceleration of ACRs. Numerical and analytical studies suggest that the possible solution can be made by considering the temporary and/or spatial variation (Florinski & Zank 2006; McComas & Schwadron 2006; Jokipii & Kóta 2008; Kóta & Jokipii 2008; Schwadron et al. 2008). In particular, McComas & Schwadron (2006) discussed the importance of the magnetic geometry of a blunt shock on particle acceleration. They argued that the missing ACRs at the nose of the Heliospheric termination shock are due to particle energization occurring primarily back along the flanks of the shock where magnetic field lines have had a longer connection time and higher injection efficiency. Kóta & Jokipii (2008) presented a more sophisticated simulation which gives results similar to that described by McComas & Schwadron (2006). Schwadron et al. (2008) also developed a three-dimensional analytic model for particle acceleration in a blunt shock, including perpendicular diffusion and drift motion due to large-scale shock structure.

Large-scale magnetic field line meandering is ubiquitous in the heliosphere and other astrophysical environments (Jokipii 1966; Jokipii & Parker 1969; Parker 1979). The acceleration of charged particles in collisionless shocks has been shown to be strongly affected by magnetic-field turbulence at different scales (Giacalone 2005a, 2005b; Giacalone & Neugebauer 2008; Guo & Giacalone 2010). The large-scale magnetic-field variation will have important effects on the shock acceleration since the transport of charged particles is different in the directions parallel and perpendicular to the magnetic field, as shown in early work (Jokipii 1982, 1987). The blunt shocks and shocks with fluctuating front (Li & Zank 2006), which have similar geometry, are also relevant to this problem. In this study, we analyze the effect of the large-scale spatial variation of magnetic field on DSA by considering a simple system which captures the basic physical ideas.

The DSA can be studied by solving the Parker transport equation (Parker 1965), which describes the evolution of the quasi-isotropic distribution function f(x_i, p, t) of energetic particles with momentum p dependent on the position x_i and time t including the effects of diffusion, convection, drift, acceleration, and source particles:

$$\begin{equation} \frac{\partial f}{\partial t} = \frac{\partial }{\partial x_i}\left[\kappa _{ij}\frac{\partial f}{\partial x_j}\right] - U_{i}\frac{\partial f}{\partial x_i}-V_{d,i}\frac{\partial f}{\partial x_i}+ \frac{1}{3}\frac{\partial U_{i}}{\partial x_i} \left[\frac{\partial f}{\partial \ln p} \right] + Q. \end{equation} \tag{ 1 }$$

Here, κ_ij is the diffusion coefficient tensor, U_i is the convection velocity, and Q is a local source. The drift velocity is given by $\textbf {V}_d = (pcw/3q)\nabla \times (\textbf {B}/B^2)$, where w is the velocity of the particle, c is the speed of light, and q is the electric charge of the particle. Parker's transport equation does not include the pre-acceleration process known as the "injection problem," which has been considered by a number of authors (e.g., Kucharek & Scholer 1995; Chalov & Fahr 1996; Giacalone 2005b; Guo & Giacalone 2010). At present there is no consensus on this issue.

Kóta & Jokipii (2008) and Kóta (2010) considered analytically a model in which the upstream magnetic field is in a plane (say x, y), with average direction in the y-direction. The x-component of the magnetic field was composed of uniform sections (straight field lines) alternating in sign, which were periodic in y. They find "hot spots" and spectral effects which illustrate the effect of an upstream meandering in the magnetic field.

Here, we consider a two-dimensional (x, z) system with a planar shock at x = 0 and a sinusoid magnetic field $\textbf {B} = \textbf {B}_0 + \sin (k z)\delta \textbf {B}$. For most of this paper, we discuss the case shown in Figure 1. In this figure, the magnetic lines of force are illustrated by blue lines. The shock is denoted by the red dashed line. In the system of interest here, the magnetic field is small enough that its effects are very small (dynamic pressure/magnetic pressure ∼ρV²_w/(B²/8π) ∼ 100 in the solar wind at 1 AU). The system is periodic in the z-direction, with the magnetic field convecting from upstream (x < 0) to downstream (x>0). In the shock frame, the particles will be subjected to convection and diffusion due to the flow velocities U₁ (upstream) and U₂ (downstream), and the diffusion coefficients parallel and perpendicular to the large-scale magnetic field (κ_1(∥,⊥) and κ_2(∥,⊥)), respectively. The gradient and curvature drifts in this case are only in the direction out of the x–z plane and thus irrelevant to this study. Because of the steady velocity difference between upstream and downstream, charged particles which travel through the shock layer will be accelerated. However, since we consider the large-scale magnetic-field variation, transport of energetic particles in the fluctuating magnetic field becomes important. The diffusion coefficient in the x–z system can be expressed as

$$\begin{equation} \kappa _{ij} = \kappa _\perp \delta _{ij} - \frac{(\kappa _\perp -\kappa _\parallel)B_iB_j}{B^2}. \end{equation} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

The normalization units chosen in this study are the spatial scale X₀ = 10 AU, the upstream velocity U₁ = 500 km s⁻¹, the timescale T₀ = 3 × 10⁶ s, and the diffusion coefficients are in units of $\kappa _0 = 7.5 \times 10^{21} \rm cm^{2} \;s^{-1}$. The shock compression ratio r = U₁/U₂ is taken to be 4.0. The shock layer is considered to be a sharp variation U_x = (U₁ + U₂)/2 − (U₁ − U₂) tan h(x/th)/2 with thickness th = 1 × 10⁻³, which is required to be less than κ_xx1/U₁ everywhere in the upstream simulation domain, where κ_xx1 is the upstream diffusion coefficient normal to the shock surface. The simulation domain is taken to be [ − 2.0 < x < 2.0, − π < z < π]. The parallel diffusion coefficients upstream and downstream are assumed to be the same and taken to be κ_∥1 = κ_∥2 = 0.1 at p = p₀. The ratio between the parallel diffusion coefficient and the perpendicular diffusion coefficient is taken to be κ_⊥/κ_∥ = 0.05, which is consistent with that determined by integrating the trajectories of test particles in magnetic turbulence models (Giacalone & Jokipii 1999). The momentum dependence of the diffusion coefficient is taken to be κ ∝ p^4/3, corresponding to non-relativistic particles in a Kolmogorov turbulence spectrum (Jokipii 1971). We use a stochastic integration method, which is described in the Appendix, to obtain the numerical solution of the transport Equation (1). In Equation (1), the source function Q = Q₀δ(p − p₀)δ(x) is represented by injecting pseudo-particles at the shock x = 0 with initial momentum p = p₀. The trajectories of pseudo-particles are integrated at each time step to obtain the numerical solution. The pseudo-particles will be accelerated if they cross the shock as predicted by DSA. The time step is 1 × 10⁻⁷, which is small enough to resolve the variation of U_x at the finite shock layer. Particles which move past the upstream or downstream boundaries will be removed from the simulation. The system is periodic in the z-direction, so a pseudo-particle crossing the boundaries in z will re-appear at the opposite boundary and continue to be followed. A particle splitting technique similar to Giacalone (2005a) is used in order to improve the statistics. Although we use very approximate parameters, we note that the results are insensitive to the precise numbers. Our results are also qualitatively unchanged after allowing the injection rate to vary as a function of shock normal angle and different downstream diffusion coefficients. We also note that our model is simplified to illustrate the physics of the magnetic field variation. Other effects such as changes in plasma properties are small for problems of current interest.

3.1. A Shock Propagating Perpendicular to the Average Magnetic Field

Consider first the case where the average magnetic field is in the z-direction and the fluctuating magnetic field is δB = B₀. As shown in Figure 1, the magnetic field is convected through the shock front and is compressed in the x-direction; thus, B_z2 = rB_z1. For the sinusoid magnetic field considered in this paper, the local angle between upstream magnetic field and shock normal, θ_Bn, will vary along the shock surface. As a magnetic field line passes through the shock surface, its connection points (the points where the field lines intersect the surface of the shock) will be moving apart in the middle of the plane (z = 0) and approaching each other on both sides of the system (|z| = π). Since κ_∥ ≫ κ_⊥, the particles tend to remain on the magnetic field lines. Because the acceleration only occurs at the shock front, as the magnetic lines of force convect downstream, the particles will be trapped and accelerated at places where the connection points converge toward each other, leading to further acceleration. For the regions where the field lines separate from each other, the particles are swept away from those regions. Figure 2 displays the spatial distribution contours of accelerated particles in three energy ranges: 3.0 < p/p₀ < 4.0 (top), 8.0 < p/p₀ < 10.0 (middle), and 15.0 < p/p₀ < 30.0 (bottom). The density is represented by the number of particles in simulation and its unit is arbitrary. It can be seen that "hot spots" form in the regions where the connection points are approaching each other at all energy ranges, with lobes extending along the magnetic field lines. The density of the accelerated particles at the connection-point separating region (in the middle of the plane) is clearly much smaller, although there is still a concentration of low-energy accelerated particles there since the acceleration of low-energy particles is rapid and efficient at perpendicular shocks. At higher energy ranges (middle and bottom), the lack of accelerated particles may be interpreted as due to the fact that the acceleration to high energies takes time.

Zoom In Zoom Out Reset image size

Figure 3 illustrates the profiles of the density of accelerated particles for different energy ranges at z = 0 (top) and z = π (bottom). In each panel, the black solid lines show the density of low-energy particles (3.0 < p/p₀ < 4.0), the blue dashed lines show the density of intermediate-energy particles (8.0 < p/p₀ < 10.0), and the red dot-dashed lines show the density of particles with high energies (15.0 < p/p₀ < 30.0). In connection-point separating regions z = 0 (top), it can be seen that while the downstream distribution of low-energy particles is roughly a constant, the density of particles with higher energies increases as a function of distance downstream from the shock. These particles are not accelerated at the shock layer in the center of the plane but in the "hot spots." At z = π (bottom), the density of particles of all energies decreases as a function of distance, which indicates that the accelerated particles diffuse away from the "hot spots." Since the high-energy particles have larger diffusion coefficients than the particles with low energy, it is easier for them to transport to the middle of the plane. The profile at z = 0 is similar to Voyager's observation of ACRs at the termination shock and the Heliosheath (Stone et al. 2005; Cummings et al. 2008) which shows that the intensity of the ACRs is still increasing and the energy spectrum is unfolding over a large distance after entering the Heliosheath. The same physics has been discussed by Jokipii & Kóta (2008), where a "hot spot" of energetic particles is produced by the spatial variation of the injection of the source particles. In this work, the concentration of energetic particles is a consequence of particles accelerated in a shock containing large-scale magnetic-field variation.

Zoom In Zoom Out Reset image size

The top panel in Figure 4 represents the positions in the z-direction and the times as soon as the particles reached a certain momentum p_c = 3p₀. We show that particles are accelerated mainly at the connection-point converging region. There are also particles accelerated at the middle of the plane because the particles can gain energy rapidly at perpendicular or highly oblique shock due to the smallness of the perpendicular diffusion coefficient (Jokipii 1987); however, further acceleration is suppressed by the effect that the charged particles travel away from the connection-point separating region; see also the top panel in Figure 5. It is clear that since the particles tend to follow the magnetic field lines, when the field line connection points separate from each other as the field convects through the shock, the particles travel mainly along the magnetic field and away from the middle of the plane. The characteristic time for a field line convecting from upstream to downstream is τ_c = D/U₁ ∼ 1; therefore, there is no significant acceleration in the middle of the plane after t = τ_c. Some of the particles can get more acceleration traveling from other regions to "hot spots." Figure 4 (bottom) shows the distance particles traveled in the z-direction from their original places |z − z₀| versus the time when the particles get accelerated at a certain energy (p_c = 3p₀). It shows that many particles are accelerated close to their original position, which is related to the acceleration in the "hot spot." Nevertheless, there are also a number of particles that travel from the connection-point separating region to the "hot spot" and get further accelerated, which is represented by the particles that travel a large distance in the z-direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the same plot as Figure 4, except here the critical momentum is p_c = 10.0p₀. It is shown again in Figure 5 (top) that most of the particles accelerated to high energies are in the hot spot. However, in contrast to Figure 4 (top), there are very few particles accelerated at the center of the plane since energetic particles transport away from the middle region and the time available is not long enough. For a quick estimate, the acceleration time is approximately

$$\begin{eqnarray} \tau _{\rm acc} &=& \frac{3}{U_1-U_2}\int ^p_{p_0}\left(\frac{\kappa _{xx1}}{U_1}+\frac{\kappa _{xx2}}{U_2}\right) d\ln p > \frac{3}{U_1-U_2}\nonumber\\ &&\times\,\int ^{10p_0}_{p_0}\left(\frac{\kappa _\perp }{U_1}+\frac{\kappa _\perp }{U_2}\right) d\ln p > \tau _c. \end{eqnarray} \tag{ 3 }$$

Therefore, most of the particles do not have sufficient time to be accelerated to high energies at the center. A number of particles accelerated at the center will travel to the "hot spot" and get more acceleration, as shown in Figure 5 (bottom).

Clearly, this presents different pictures of particle acceleration by the shock containing two-dimensional spatial magnetic field variations, indicates that the resulting distribution function is spatially dependent. In Figure 6, we show the steady-state energy spectra obtained in the following regions. Top: [0.1 < x < 0.3, π − 0.2 < z < π] (black solid line) and [0.1 < x < 0.3, −0.1 < z < 0.1] (green dashed line); bottom: [0.8 < x < 1.0, π − 0.2 < z < π] (black solid line) and [0.8 < x < 1.0, −0.1 < z < 0.1] (green dashed line). It is shown that the spatial difference among distribution functions at different locations caused by large-scale magnetic-field variation is considerable. The black lines in both the top and bottom plot, which correspond to the "hot spots," show power-law-like distributions except at high energies. At high energies, the particles will leave the simulation domain before gaining enough energy which causes the roll over in the distribution function; this roll over is mainly caused by a finite distance to upstream boundary. For other locations, the two-dimensional effect we discussed will produce the modification in distribution functions. The most pronounced effect can be found at the nose of the shock (green lines), the distribution of which in the top panel shows a suppression of acceleration at all the energies. This insufficient acceleration is most prominent in the range of 6–12p₀. At these energies, the acceleration timescales are longer than the time for the field line convection and swipe the particles away from the connection-point separating region, as we discussed above. The bottom plot shows that deep downstream the spectrum of accelerated particles is similar at high energies because of the mobility of these particles.

Zoom In Zoom Out Reset image size

3.2. An Oblique Shock

The previous discussion has established the effect of a spatially varying upstream magnetic field on the acceleration of fast charged particles at a shock which is propagating normal to the average upstream magnetic field. We next consider the case where the shock propagation direction is not normal to the magnetic field.

Clearly, if the varying direction of the upstream magnetic field is such that at some places the local angle of the magnetic field relative to the average field direction exceeds the angle of the average magnetic field to the shock plane, we will have situations similar to that discussed in the previous sections. There will be places where the connection points of the magnetic field to the shock move further apart or closer together. Hence we expect that the same physics can be applicable. An example is given in Figure 7. In this case the ratio δB/B₀ is taken to be 0.5, the averaged shock normal angle θ_Bn = 70°. It can be seen from this plot that the connection points can still move toward each other in some regions. Figure 8 shows the density contours of accelerated particles the same as Figure 2, but for the case of the oblique shock. We find that in this case the process we discussed in the previous section is still persistent, even for an oblique shock and relatively smaller δB/B₀. The "hot spot" forms correspond to the converging magnetic connection points, and particle acceleration is suppressed in the region where connection points separate from each other. We may conclude that, for a shock which is oblique, if some magnetic field lines can intersect the shock multiple times, we have "hot spots" of accelerated particles where the connection points are converging together.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The acceleration of charged particles in shock waves is one of the most important unsolved problems in space physics and astrophysics. The charged particle transport in turbulent magnetic field and acceleration in shock region are two inseparable problems. In this paper we illustrate the effect of a large-scale sinusoidal magnetic-field variation. This simple model allows a detailed examination of the physical effects. As the magnetic field lines pass through the shock, the connection points between field lines on the shock surface will move accordingly. We find that the region where connection points are approaching each other will trap and preferentially accelerate particles to high energies and form "hot spots" along the shock surface, somewhat in analogy to the "hot spots" postulated by Jokipii & Kóta (2008). The shock acceleration will be suppressed at places where the connection points move apart from each other. Some of the particles injected in those regions will transport to the "hot spots" and get further accelerated. The resulting distribution function is highly spatial dependent at the energies we studied, which could give a possible explanation for the Voyager observation of ACRs. Although we have discussed a simplified, illustrative model, the resulting spectra and radial distributions show qualitative similarity with the in situ Voyager observations. Thus, the intensities do not, in general, peak at the shock, and the energy spectra are not power laws. We show that this process is robust even in the case of oblique shocks with relatively small magnetic field variations. Large-scale magnetic-field variation, which could be due to magnetic structures like magnetic clouds, or the ubiquitous large-scale field line random walk will strongly modify the simple planar shock solution. This effect could work in a number of situations for large-scale shock acceleration including magnetic variations, for example, the solar wind termination shock and supernova blast waves.

F.G. thanks Joe Giacalone for the tutorial and discussion on the stochastic integration method. F.G. also thanks Chunsheng Pei and Erica McEvoy for valuable discussion on stochastic differential equations. We acknowledge partial support from NASA under grants NNX08AH55G, NNX08AQ14G, and NNX09AB13G, and NSF under grant ATM 0641600.

In the two-dimensional system considered, the Parker transport Equation (1) can be written in a Fokker–Planck form as

$$\begin{eqnarray} \frac{\partial f}{\partial t}&=& \frac{\partial ^2}{\partial x^2}\left(\kappa _{xx} f \right) + \frac{\partial ^2}{\partial z^2}\left(\kappa _{zz} f\right) + \frac{\partial ^2}{\partial x \partial z}\left(2\kappa _{xz}f \right) \nonumber \\ & &-\frac{\partial }{\partial x}\left[\left(U_{x} + \frac{\partial \kappa _{xx}}{\partial x} + \frac{\partial \kappa _{xz}}{\partial z} \right)f\right]\nonumber \\ && -\frac{\partial }{\partial z}\left[\left(\frac{\partial \kappa _{zz}}{\partial z} + \frac{\partial \kappa _{xz}}{\partial x} \right)f\right] +\frac{1}{3}\frac{\partial U_{x}}{\partial x} \left[\frac{\partial f}{\partial \ln p} \right] + Q,\nonumber\\ \end{eqnarray} \tag{ A1 }$$

where κ_xx = Δx²/2Δt, κ_zz = Δz²/2Δt, and κ_xz = ΔxΔz/2Δt. Following the usual approach in stochastic integration (Jokipii & Owens 1975; Jokipii & Levy 1977), the solution can be calculated by successively integrating trajectories of pseudo-particles using:

$$\begin{eqnarray} \Delta x &=& r_1 (2\kappa _\perp \Delta t)^{1/2} + r_3 (2(\kappa _\parallel -\kappa _\perp)\Delta t)^{1/2}\frac{B_x}{B} + U_x\Delta t\nonumber\\ && + \left(\frac{\partial \kappa _{xx}}{\partial x} + \frac{\partial \kappa _{xz}}{\partial z}\right) \Delta t\vspace*{-10pt} \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} \Delta z &=& r_2 (2\kappa _\perp \Delta t)^{1/2} + r_3 (2(\kappa _\parallel -\kappa _\perp)\Delta t)^{1/2}\frac{B_z}{B}\nonumber\\ && + \left(\frac{\partial \kappa _{zz}}{\partial z} + \frac{\partial \kappa _{xz}}{\partial x}\right) \Delta t\vspace*{-10pt} \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} &&\Delta p = -\frac{p}{3} \frac{\partial U_x}{\partial x} \Delta t, \end{eqnarray} \tag{ A4 }$$

where r₁, r₂, and r₃ are different sets of random numbers which satisfy 〈r_i〉 = 0 and 〈r²_i〉 = 1. It can be easily demonstrated that the ensemble average of stochastic differential Equations (A3), (A4), and (A5) is the solution of transport Equation (A1). In order to study DSA, we approximate the shock layer as a sharp variation U_x = (U₁ + U₂)/2 − (U₁ − U₂) tan h(x/th)/2 with a thickness th much smaller compared with the characteristic length of diffusion acceleration κ_xx1/U₁. At the same time, we have to make sure that the time step Δt is small enough to resolve the motion in the shock layer.