PARTICLE-IN-CELL SIMULATIONS OF PARTICLE ENERGIZATION VIA SHOCK DRIFT ACCELERATION FROM LOW MACH NUMBER QUASI-PERPENDICULAR SHOCKS IN SOLAR FLARES

Low Mach number (M), high plasma beta (β_p) fast mode shocks can occur in the magnetic reconnection outflows of solar flares. Hard X-ray data from Yohkoh and RHESSI have revealed that electrons are energized above flare loop tops and footpoints (e.g., Lin et al. 2003). Solar flares are diverse and the associated reconnection events are likely "acceleration environments" with potentially different mechanisms of particle acceleration operating on different scales. One potential source of particle acceleration, seen from analytic predictions (Blackman & Field 1994) and numerical simulations of reconnection configurations in which an obstacle was present (e.g., Forbes 1998; Workman et al. 2011), involves the presence of low Mach number fast shocks in reconnection outflows. These "termination" shocks may contribute to the high energy acceleration observed over some frequency range, studying the potential ways in which these shocks can accelerate particles is well motivated.

Mann et al. (2006, 2009) and Warmuth et al. (2009) suggested an electron energization via shock drift acceleration (hereafter SDA) in termination shocks. Guo & Giacalone (2010, 2012) performed hybrid simulations for termination shocks where test electrons were effectively energized via the interaction with pre-existing large-scale magnetic fluctuations. In our recent work (Park et al. 2012), we performed a full particle-in-cell (PIC) simulation for exactly perpendicular (i.e., magnetic field perpendicular to the shock normal), low M/high β_p shocks. We found that both electrons and ions participated in SDA.

Termination shocks may, however, deviate from exactly perpendicularity in solar flares, the extent to which is an open question. In the meantime, it is instructive to relax the constraint that the shocks are exactly perpendicular and consider the case of quasi-perpendicularity. The difference is significant because particles can cross back upstream for quasi-perpendicular shocks. There were analytical (Wu 1984; Krauss-Varban et al. 1989b) and hybrid simulation studies (Krauss-Varban et al. 1989a) of SDA in a nearly perpendicular bow shock for energetic electrons. Recently, Matsukiyo et al. (2011) performed one-dimensional (1D) PIC simulation for the electron SDA in low Mach number quasi-perpendicular shocks in galaxy clusters.

In this paper, we present the results of two-dimensional (2D) PIC simulations for quasi-perpendicular low M/high β_p shocks in solar flares. The upstream magnetic field makes angles of θ_B = 80°, 82°, and 835 respectively to the shock normal and the shocks satisfy the subluminal condition, V_sh/cos θ_B < c, where V_sh is the shock speed in the upstream rest frame and c is the speed of light. In such a subluminal shock, electrons can then be reflected at the shock front due to the magnetic mirror effect and gain energy (e.g., Ball & Melrose 2001; Mann et al. 2006, 2009; Warmuth et al. 2009).

One difficulty in performing PIC simulations for SDA in solar flares is that they require sufficiently high energy electrons above the threshold energy for SDA in the simulation to obtain statistically reliable results. The commonly used Maxwellian distribution has too few above the threshold electrons to be used in a simulation. Herein we use a kappa distribution with κ = 10 for both injected ions and electrons. This increases the number of high energy electrons but does not significantly alter the original shock structure, which is determined by the bulk of the thermal particles. Although the kappa distribution is used for computational convenience in this paper, the kappa distribution may be physically a relevant distribution in the outflow driven electron–ion jet (e.g., Yoon et al. 2006; Kas˘parová & Karlický 2009; Mann et al. 2006, 2009; Warmuth et al. 2009).

The goals of this paper are two-fold: (1) to study the formation and structure of such shocks, including the turbulent dissipation mechanism for collisionless shock sustenance and entropy creation, and (2) to study the particle acceleration mechanism relevant for the soft/hard X-ray flux observations in solar flares. We observe the same modified-two-stream instability (Krall & Liewer 1971) in the shock transition region as in the perpendicular shocks (Park et al. 2012) that can provide the turbulent dissipation in the downstream. Furthermore, a temperature anisotropy after the shock transition region is found to drive the electron whistler instability (e.g., Gary 1993).

As shown below, more SDA-accelerated electrons are found in the present simulations of the quasi-perpendicular case compared to the previously simulations of the perpendicular case (Park et al. 2012). In this context, we also extend the theoretical analysis in Mann et al. (2006, 2009) and Warmuth et al. (2009) to include the electric potential jump at the shock front (e.g., Ball & Melrose 2001) and generalize the electron energy spectrum. Both our theoretical analysis and our simulations show a transition energy, E_{trans, p}, between the thermal and non-thermal photon spectrum that is determined by the minimum angle of θ_B. We find that the transition energy E_{trans, p} and the spectral index δ from theory and our simulation are consistent with some of the X-ray spectra of solar flares from RHESSI (e.g., Altyntsev et al. 2012) in the energy range of E ≲ 40 ∼ 100 keV. Beyond this range, to maintain the power-law distribution up to E ∼ MeV, additional mechanisms beyond SDA are required.

The rest of the paper is organized as follows. The simulation setup is described in Section 2. The shock structure and particle acceleration are described in Section 3. We summarize in Section 4.

We use the fully-relativistic full PIC code OSIRIS (Fonseca et al. 2002) to study the formation of and particle energization in low-M, high-β_p shocks, where M is the Mach number and β_p is the ratio of thermal to magnetic pressure. To launch a shock, we use the moving wall boundary condition (Langdon et al. 1988; Park et al. 2012) at the right boundary of the 2D simulation box. The moving wall method generates a slowly propagating shock compared to the more standard fixed reflection boundary method, and allows for smaller box sizes and more efficient use of simulation time.

We adopt parameters typical of those found in solar flare reconnection outflows (Tsuneta 1996; Workman et al. 2011) as the upstream conditions for our shock. In particular, we use a plasma density n = 5 × 10⁹ cm⁻³, electron and ion temperatures T_e = T_i = 0.8 keV(= 9.27 × 10⁶ K), and the magnetic field strength B = 6 G with β_p ≡ 8πn(T_e + T_i)/B² = 8.93. The magnetic field is in the x-y plane $({\bf B}=B_x\hat{x}+B_y\hat{y})$ and has an angle of θ_B = 80°, 82°, and 835 from the shock normal (x-axis) in each simulation. We also performed another simulation where the upstream magnetic field is directed out of the simulation plane $(B=B_x\hat{x}+B_z\hat{z})$ and has an angle of θ_B = 80° to compare the electron energy spectrum with that in the "in-plane" B-field simulation. A reduced ion/electron mass ratio of m_i/m_e = 30 is used to reduce computational demands. The Alfvén Mach number is chosen to be $M_A\equiv V_1\sqrt{4\pi m_in}/B=6.62$, which equates to an upstream plasma flow velocity in the shock rest frame, V₁ = 0.032c, for m_i/m_e = 30. If the real mass ratio, m_i/m_e = 1836, is used for the same M_A, then the upstream flow velocity would be 1226 km s⁻¹(= 0.0041c). Effects of a realistic mass ratio are further discussed in Section 3.3. The fast-magnetosonic Mach number M satisfies $M\equiv M_A/\sqrt{1+(5/6)\beta _p}=2.28$. The ratio of the electron cyclotron frequency to the electron plasma frequency is Ω_ce/ω_pe = 0.0265. With these upstream values of M_A and β_p, the Rankine–Hugoniot relation (Tidman & Krall 1971) for the shocks with the angles used gives a compression ratio in the range of (2.06, 2.50), where the lower and the upper limits are calculated for 2D and 3D, respectively.

The simulation box is initialized with a kappa-distributed ion–electron plasma drifting with V_d = 0.0213c and T_e = T_i = 0.8 keV, where V_d is set to a smaller value than the upstream speed V₁(= 0.032c) in the shock rest frame anticipating that the shock will be traveling to the left. In the upstream rest frame, the energy distribution is a kappa distribution such as

$$\begin{eqnarray} f(E) &=& {2^{5/2}\over \sqrt{\pi }}{1\over (2\kappa -3)^{3/2}} {\Gamma (\kappa +1)\over \Gamma (\kappa -1/2)} {E^{1/2}\over T^{3/2}} \nonumber \\ && \times \left(1+{2\over 2\kappa -3}{E\over T}\right)^{-\kappa -1}, \end{eqnarray} \tag{ 1 }$$

where Γ is the gamma function and f(E) is normalized to 1, ∫^∞₀dEf(E) = 1. The kappa distribution approaches the Maxwellian distribution as κ goes to ∞. Figure 1 shows the initial energy distribution in the upstream rest frame from the simulation. The dashed line is a kappa distribution in Equation (1) with κ = 10 and T = 0.8 keV and the dotted line is a Maxwellian distribution with the same temperature. The bulk parts (E < 5 keV) of the two distributions are very similar. The implementation of the kappa distribution in OSIRIS is described in Appendix A.

Zoom In Zoom Out Reset image size

A uniform external E_z field is set up along the z-axis with E_z = −V_dB_y/c. A new plasma of the same kappa distribution is constantly injected from the left boundary (x = 0) throughout the simulation. The simulation box sizes are L_x × L_y = 250c/ω_pe × 50c/ω_pe. The grid size used is dx = dy = 0.08c/ω_pe and the time step used is dt = 0.056/ω_pe. For each particle species, 196 particles per cell are used. The particle number/cell is fairly large to reduce numerical collisions and maintain the kappa distribution for a sufficiently long time.

A linear current deposition scheme is used for all simulations in this paper. A periodic boundary condition is used in the y-direction for both particles and fields. For fields, an open boundary condition is used in the x-direction. Particles that reach x = 0 are re-injected into the box with the initial kappa distribution. At x = L_x, a moving wall boundary condition is adopted (Park et al. 2012) and the moving wall speed is set to V_wall = 0.004c.

3.1. Shock Structure

In this subsection, we present results for the shock structure for θ_B = 80° in-plane case only but other simulations for θ_B = 82° and 835 also show similar shock structures. The electron energization via SDA shows different energy spectra for the different angles of θ_B as will be seen in the next subsection.

In Figures 2(a) and (b), we plot for both species the upstream to downstream ratios of density n(x)/n₁, flow speeds V₁/V(x), and y-magnetic field, B_y(x)/B_y1. In Figures 2(c)–(l), we plot the phase–space distributions of p_x–x, p_y–x, and p_z–x, the y-averaged flow velocity profiles of V_x(x), V_y(x), and V_z(x), and the temperature profiles of T_∥(x), T_⊥(x), and T_z(x) (where ⊥ and ∥ are perpendicular and parallel to the magnetic field ${\bf B}(x)\equiv B_x(x)\hat{x}+B_y(x)\hat{y}$, respectively), for the electrons (the left column) and the ions (the right column) at t = 11200/ω_pe for the θ_B = 80° case.

Zoom In Zoom Out Reset image size

In Figure 2, the shock front is at x ≈ 140c/ω_pe and moves to the left with a speed of 0.01c in the simulation frame. In Figures 2(a) and (b), the compression ratio is r = 2.5 near the shock front and oscillates around 2.0 downstream, in reasonable agreement with the Rankine–Hugoniot relation (Tidman & Krall 1971). Figures 2(c)–(h) show that the electrons and ions are heated downstream but the ions are conspicuously less heated in the y-direction (Figure 2(f)). In Figures 2(i) and (j), the flow velocity, V_x = 0.021c, in the upstream decreases downstream to the moving wall speed, V_wall = 0.004c. The z-component of the velocity, V_z, is oscillating around zero downstream due to the E_x × B_y-drift while V_y ≃ 0. Figures 2(k) and (l) show that the downstream temperature is anisotropic, that is T_⊥ > T_∥ with T_e⊥ = T_ez ≈ 1.66 keV and T_e∥ ≈ 1.5 keV for the electrons, and T_i⊥ = T_iz ≈ 2.5 keV and T_i∥ ≈ 1.0 keV for the ions.

Figures 2(m)–(o) show B in the x–y space. The B_x and B_z fields show oscillatory patterns along the y-axis near the shock front (x ∼ 150c/ω_pe) and the B_y field shows a rippled surface along the y-axis. These variations are due to electron temperature anisotropy-driven instabilities, T_e⊥ > T_e∥.

Figure 3 shows the y-averaged E, B, and the potential energy $|e|\Phi (x)(=e\int _{x_1}^{x}E_xdx)$ at t = 11200/ω_pe, measured in the shock rest frame. The E_y, E_z, and B_x fields are approximately constant across the shock while E_x and B_z oscillate with a wave number k = 0.2ω_pe/c downstream. The potential energy jump |e|ΔΦ at the shock front is ∼3.5 keV (Figure 3(d)). Some ions in the low energy tail reflect at the shock front and can drive the modified two-stream instability.

Zoom In Zoom Out Reset image size

In Figures 4(a)–(d), we plot the Fourier spectra of the E_⊥, B_⊥, B_∥, and B_z fields, where ⊥ and ∥ are perpendicular and parallel to B₀, respectively. The ${\bf B}_0(=B_{0x}\hat{x}+B_{0y}\hat{y})$ is an averaged field over 140 < x < 170 (c/ω_pe) in Figure 2 and has an angle 86° from the x-axis. Two distinct dominant modes are excited by temperature anisotropy-driven instabilities (T_⊥ > T_∥), one is along the B₀ axis with k = 0.5ω_pe/c and another is sightly deviated from the B₀ axis with k = 0.8ω_pe/c. In Figure 4(d), the strong signals at k_x = ±0.2ω_pe/c are due to the oscillatory pattern of the B_z field along the x-axis as seen in Figure 3(f).

Zoom In Zoom Out Reset image size

To analyze the spectra from the simulation, we numerically solve the dispersion relation (e.g., Gary 1993) in Figures 4(e) and (f). Here we assume bi-Maxwellian electron and ion distributions and a uniform background magnetic field B₀. We use the parameters extracted from the simulation with B₀ = 13.5G, T_e∥ = 1.25 keV and T_e⊥ = 1.66 keV for the electrons, and T_i∥ = 0.86 keV and T_i⊥ = 2.55 keV for the ions, where ∥ and ⊥ are parallel and perpendicular directions to B₀, respectively.

In Figure 4(e), the electron whistler modes are centered on the B₀ axis with k = 0.45ω_pe/c, the ion cyclotron modes centered on the B₀ axis with k = 0.1ω_pe/c, and the ion mirror modes are obliquely off the B₀ axis with k = 0.12ω_pe/c. The maximum growth rates are ω_i = 5.6 × 10⁻⁴ω_pe, 2.5 × 10⁻⁴ω_pe, and 2.0 × 10⁻⁴ω_pe for the electron whistler, the ion cyclotron, and the ion mirror modes, respectively. In Figure 4(f), the real frequencies for the electron whistler and the ion cyclotron modes are around ω_r = 0.012ω_pe and 0.0012ω_pe, respectively, and the ion mirror modes have zero real frequency. These are consistent with the analytical results (e.g., Gary 1993), Ω_ci < ω_r < |Ω_ce| for the electron whistler modes and 0 < ω_r < Ω_ci for the ion cyclotron modes.

The linear theory based on bi-Maxwellian distributions with the parameters found in the PIC simulations show that electron mirror modes are marginally stable (e.g., Gary & Karimabadi 2006). The oblique modes found in Figures 4(a)–(d) could be due to some unknown modes. However, given the fact that the actual distribution is not exactly a bi-Maxwellian (e.g., with flows) and the uncertainty in the temperature anisotropy measurements in the simulation, the oblique modes with k = 0.8ω_pe/c seen in the simulation (Figures 4(a)–(d)) may possibly still be a electron mirror-type of modes for a non-bi-Maxwellian distribution. In addition, the ion cyclotron/mirror modes are not observed in this simulation because of the small box size L_y = 50c/ω_pe and/or their relatively weaker growth.

Figures 5(a) and (b) show the electron and ion distributions in the shock transition region in 135 < x < 155(c/ω_pe) at t = 11200/ω_pe from the simulation. We observe that 27% of the incoming ions are reflected at the shock front (Figure 5(a)) and the modified two-stream instability (MTSI; e.g., Matsukiyo & Scholer 2003, 2006; Umeda et al. 2010, 2012) can be excited. Here we solve the MTSI in the electrostatic limit with its wave vector along the x-axis (Appendix B). The electrons have a temperature of T_ex = 1.65 keV and drift with V_ex = 0.0064c in the simulation frame (Figure 5(b)). We fit the distributions using drifting Maxwellians. The electrons are magnetized with B = 10 G while the ions are assumed to be non-magnetized. In the electron rest frame, the drift velocities for the incoming and reflecting ions are V_xin = 0.0116c and V_xre = −0.0134c, respectively. Both incoming and reflecting ions have a temperature of T_xin = T_xre = 0.98 keV. In Figures 5(c) and (d), we solve the kinetic dispersion relation for the MTSI. The maximum growth rate is ω_i = 3.2 × 10⁻⁴ω_pe at k = 0.2c/ω_pe and the real frequency is ω_r = −0.0015ω_pe.

Zoom In Zoom Out Reset image size

Turbulent dissipation is needed to randomize the upstream flow to form collisionless shocks (e.g., Wu 1982). The macroscopic jump conditions across a shock are essentially independent of the source of microphysical turbulence, as long as there is such a source. Plasma instabilities in the shock transition region are natural sources for this needed turbulent dissipation and we consider three possible instabilities: (1) lower hybrid instability (the excited modes are in the B × ∇B direction; e.g., Zhou et al. 1983), (2) whistler instability, and (3) MTSI (e.g., Papadopoulos et al. 1971; Wagner et al. 1971). The lower hybrid instability is precluded in these 2D in-plane simulations where the B-field is in the x–y simulation plane. In addition, we notice that the shock structure (i.e., the compression ratio) in the simulation with an in-plane B field is the same as that with the out-of-plane B field which then precludes the whistler instability. We are led to conclude that the MTSI is the most likely candidate to provide the needed turbulence (e.g., Park et al. 2012).

3.2. Particle Heating via Shock Drift Acceleration

In Figure 6, we plot the electron and ion energy distributions, f(E)(= dN/dE) upstream (50 < x < 155c/ω_pe ((b) and (e))) and downstream (155 < x < 250c/ω_pe ((c) and (f))) at t = 8400/ω_pe for θ_B = 80°. We fit the thermal (bulk) part of the distributions using a kappa distribution with κ = 10 (dashed lines), yielding T_e = T_i = 0.8 keV upstream and T_e = 1.58 keV and T_i = 1.8 keV downstream. In Figures 6(a) and (d), we plot the electron and ion phase–space distributions in E − x, overlaid with B_y field to indicate the location of the shock front. Abundant non-thermal electrons are seen ahead of the shock at x = 155c/ω_pe, traveling as far back as x = 10c/ω_pe (Figure 6(a)).

Zoom In Zoom Out Reset image size

In Figure 6(b), the electron energy spectrum shows a deviation from a thermal distribution at E ∼ 12 keV and the spectrum becomes steeper at E ∼ 50 keV. The dotted line in Figure 6(b) is the theoretical energy distribution via SDA when the electric potential energy jump |e|ΔΦ ≈ 3.5(keV) at the shock front is considered (see Section 3.2.1). In Figure 6(e), the ion energy spectrum shows a deviation from a thermal distribution at E ∼ 15 keV and the spectrum becomes steeper at E ∼ 30 keV.

In Figures 7(a)–(d), we plot a typical track of an electron experiencing SDA, overlaid with the shock front and the subsequent compression peaks (Figure 7(a)). The electron is reflected at the shock front (Figure 7(a)) and gains energy from 8 keV to 50 keV (Figure 7(c)). After the reflection, the electron drifts along the upstream magnetic field lines (Figure 7(b)). In Figure 7(d), we plot the electron in the v_⊥ − v_∥ phase–space. The electrons in the region I are transferred to the region II after the reflection.

Zoom In Zoom Out Reset image size

In Figures 7(e) and (f), we plot a typical track for an ion experiencing SDA. Whether or not an ion gains energy via SDA depends on its incident speed and angle of incidence at the shock front (Kirk 1994). When the ion meets the shock front, it turns back toward the upstream with a larger gyro-radius (Figure 7(a)) and is accelerated by the E_z field. The kinetic energy of the ion increases from 5 keV up to 20 keV (Figure 7(b)). A detailed analysis for ions experiencing SDA in perpendicular shocks was described in our previous work (Park et al. 2012).

3.2.1. Electron Spectrum via SDA

In this subsection, we generalize the electron energy spectrum via SDA by Mann et al. (2006, 2009) and Warmuth et al. (2009) to include the electric potential energy, eΦ, at the shock front (e.g., Ball & Melrose 2001) and compare with the PIC simulation results.

First, we consider the de Hoffmann–Teller (dHT) frame (denoted by ') where the motional E_z(= −V₁/cB_y) field vanishes. The dHT frame is obtained by boosting with v_s = V_sh/cos θ_B along the magnetic field line in the upstream rest frame. (Here we consider the negative shock speed, V_sh < 0, for a shock traveling to the -$\hat{x}$ direction.) The maximum θ_B for the existence of the dHT frame is given by θ_Bmax = cos ⁻¹(V_sh/c)(< 90°).

In the dHT frame, the condition for an electron to reflect at the shock front is (e.g., Ball & Melrose 2001; Mann et al. 2006, 2009)

$$\begin{equation} \beta _{\perp }^{\prime }>\beta _{\parallel }^{\prime } \tan \alpha _0 \,\,\,{\rm and}\,\,\, \beta _{\parallel }^{\prime }>0, \end{equation} \tag{ 2 }$$

where ⊥ and ∥ are perpendicular and parallel to the upstream B₁, respectively, β'_{∥, ⊥} = v'_{∥, ⊥}/c, and $\alpha _0=\sin ^{-1}\sqrt{B_1/B_2}(\approx \sin ^{-1}\sqrt{1/r})$. If we include the potential energy eΦ(< 0) at the shock front, the reflection condition, Equation (2), can be written in the non-relativistic limit by (e.g., Ball & Melrose 2001)

$$\begin{equation} \beta _{\perp }^{\prime }>\sqrt{\left(\beta ^{\prime 2}_{\parallel }-{2e\Phi \over m_e c^2}\right)} \tan \alpha _0, \,\,\, \beta _{\parallel }^{\prime }>0. \end{equation} \tag{ 3 }$$

For relativistic electrons with $\beta _{\parallel }^{\prime }\gg \sqrt{2|e|\Phi /m_e c^2}(=0.117)$, the potential energy eΦ is negligible and Equation (3) is reduced to Equation (2).

Using the Lorentz transformation between the dHT frame (') and the upstream rest frame,

$$\begin{equation} \beta ^{\prime }_\parallel ={\beta _\parallel -\beta _s\over 1-\beta _\parallel \beta _s}, \,\,\, \beta ^{\prime }_\perp ={\beta _\perp \over \gamma _s(1-v_\parallel \beta _s)}, \end{equation} \tag{ 4 }$$

where β_s = v_s/c (= V_sh/c cos θ_B) and $\gamma _s=1/\sqrt{1-\beta _s^2}$, Equation (3) is transformed in the upstream rest frame into

$$\begin{eqnarray} \beta _{\perp } \; & {>} &\; \gamma _s\tan \alpha _0 \left[(\beta _{\parallel }-\beta _s)^2-{2e\Phi \over m_ec^2} (1-\beta _s\beta _\parallel)^2\right]^{1/2}. \nonumber \\ \beta _\parallel \; & {>} &\; \beta _s. \end{eqnarray} \tag{ 5 }$$

The velocity of the electron after reflection in the dHT frame is given by β'_r∥ = −β_i∥' and β'_r⊥ = β_i⊥', where the indices i and r represent the incoming and reflected electron, respectively. In the upstream rest frame, one gets (Mann et al. 2006, 2009)

$$\begin{eqnarray} &&\beta _{i\parallel }={2\beta _s-\beta _{r\parallel } \left(1+\beta _s^2 \right) \over 1-2\beta _{r\parallel }\beta _s+\beta _s^2}, \beta _{i\perp }={\beta _{r\perp }\over \gamma _s^2\left(1-2\beta _{r\parallel } \beta _s+\beta _s^2\right)}. \quad \end{eqnarray} \tag{ 6 }$$

In Figure 8, we plot the SDA condition in the β_⊥ − β_∥ space for β_s = −0.293 as an example (e.g., Mann et al. 2006, 2009). Equation (5) defines the region I (light shaded region) and the hyperbolic curve is given by the potential energy eΦ = −3.5 keV at the shock front. When eΦ goes to zero, the hyperbolic curve approaches the straight dotted lines as seen in Mann et al. (2006, 2009). Equation (6) implies that the incoming electrons in the region I are transferred to the region II (dark shaded region) after reflection. As θ_B increases to θ_Bmax(= cos ⁻¹(V_sh/c)), β_s goes to −1 in Figure 8 and the number of electrons satisfying the reflection condition in Equation (5) decrease to zero. Therefore, no electron reflects at the shock front in the superluminal shocks (where |β_s| ⩾ 1) or the perpendicular shocks.

Zoom In Zoom Out Reset image size

The threshold energy E_thres for SDA is given by the shortest distance from the origin to the hyperbolic curve, $\overline{oq}(=\beta _{{\rm thres}})$ in Figure 8, namely

$$\begin{eqnarray} E_{\rm thres} &=& \left(1/\sqrt{1-\beta _{\rm thres}^2}-1\right)m_ec^2, \,\, {\rm where} \nonumber \\ \beta _{{\rm thres}} &=& \sqrt{{D^2\left[(1-\beta _s)^2D^2P^2 \left\lbrace P- \left(1+\beta _s^2 \right) \right\rbrace +\beta _s^2-P\right]\over 1+D^2 \left(1-P\beta _s^2 \right)}}, \nonumber\\ \end{eqnarray} \tag{ 7 }$$

where D = γ_stan α₀ and P = 2eΦ/m_ec².

The transition energy E_{trans, e} between thermal and non-thermal electron populations is given by the distance from the origin to the point p, $\overline{op}(=\beta _{{\rm trans}})$, in Figure 8 such as

$$\begin{eqnarray} E_{{\rm trans},e}&=\left(1/\sqrt{1-\beta _{\rm trans}^2}-1\right)m_ec^2, \,\,{\rm where} \\ \beta _{{\rm trans}}&=\beta _s^2-\left(1-\beta _s^2 \right)\left(2e\Phi /m_ec^2 \right)\tan ^2\alpha _0. \nonumber \end{eqnarray} \tag{ 8 }$$

The energy point E_{2, e} where the maximum energy ratio of the reflected to the incoming electron occurs is given by Ball & Melrose (2001) in the non-relativistic approximation, namely

$$\begin{eqnarray} E_{2,e}&\approx \left(1/\sqrt{1-\beta _2^2}-1\right)m_ec^2, \nonumber \\ \beta _2&=2\beta _s\cos \,(\alpha _0/2). \end{eqnarray} \tag{ 9 }$$

Beyond E = E_{2, e}, the spectral index δ, defined as in f(E)∝E^−δ, increases. Figure 9 shows how the transition energy E_{trans, e} for eΦ = −3.5 keV (solid) and 0 keV (dashed), and E_{2, e} (dot-dashed) varies with θ_B when V_sh = 0.0041c, r = 2.5 and m_i/m_e = 1836 (M_A = 6.62 and β_p = 8.93).

Zoom In Zoom Out Reset image size

The reflected electron distribution in the upstream rest frame is written as

$$\begin{eqnarray} f_r(\boldsymbol{\beta }_r) &=& f_i[\boldsymbol{\beta }_i(\boldsymbol{\beta }_r)](d^3\beta _i/d^3\beta _r) \Theta (\beta _s-\beta _{r\parallel }) \nonumber \\ &&\times \Theta \left(\beta _{r\perp }-\gamma _s\tan \alpha _0 \left[(\beta _s-\beta _{r\parallel })^2 \nonumber -{2e\Phi \over m_ec^2} (1-\beta _s\beta _{r\parallel })^2\right]^{1/2}\right), \end{eqnarray} \tag{ 10 }$$

where Θ is the step function and the term d³β_i/d³β_r is given by the Jacobian determinant,

$$\begin{equation} {d^3\beta _i\over d^3\beta _r}={\beta _{i\perp }\over \beta _{r\perp }} \left|{\partial (\beta _{i\perp },\beta _{i\parallel })\over \partial (\beta _{r\perp },\beta _{r\parallel })}\right|= {(1-\beta _s)^4\over \left(1-2\beta _{r\parallel } \beta _s+\beta _s^2 \right)^{4}}. \end{equation} \tag{ 11 }$$

Here we let the incoming distribution in Equation (10) be a semi-relativistic kappa distribution (Mann et al. 2006, 2009),

$$\begin{equation} f_i(\boldsymbol{\beta }_i)=c_\kappa \left(1+{2(\gamma _i-1)\over 2\kappa -3}{m_ec^2\over T}\right)^{-\kappa -1}, \end{equation} \tag{ 12 }$$

where $\gamma _i=1/\sqrt{\vphantom{A^A}\smash{\hbox{${1-\boldsymbol{\beta }_i^2}$}}}$ and $c_\kappa =1/\int d^3\beta _i f_i(\boldsymbol{\beta }_i)$.

We calculate the upstream electron energy distributions, f(E) = f(β)dβ/dE, in the upstream rest frame using Equations (6)–(12) to obtain

$$\begin{eqnarray} &&f(\beta)=\left\lbrace \begin{array}{@{}cc}2\pi \beta ^2\int _{-1}^{1}dt f_i({\beta }) \\ \mbox{................. for}\,\,\beta <\beta _{\rm thres} \\ 2\pi \beta ^2\left(\int _{-1}^{t_3}dt+\int _{t_1}^{1}dt\right)f_i(\beta) \\ \mbox{................. for}\,\,\beta _{\rm thres}\le \beta <\beta _{\rm trans}, \\ 2\pi \beta ^2\int _{t_3}^{t_2}dt {d^3\beta _i\over d^3\beta _r} f_i[\boldsymbol{\beta }_i(\beta ,t)] \\ +2\pi \beta ^2\left(\int _{-1}^{t_2}dt+\int _{t_1}^{1}dt\right)f_i(\beta) \\ \mbox{................. for}\,\,\beta \ge \beta _{\rm trans} \end{array}\right. \end{eqnarray} \tag{ 13 }$$

where t = cos θ and θ is the electron's pitch angle, i.e., the angle between β and B₁. The boundaries t₁(= cos θ₁), t₂(= cos θ₂), and t₃(= cos θ₃) as shown in Figure 8 are given as t₂ = β_s/β, and roots of the equation,

$$\begin{eqnarray} &&t_{1,3}^2\left\lbrace 1+D^2\left(1-P\beta _s^2\right)\right\rbrace -2t_{1,3} D^2{\beta _s\over \beta }\left(1-P\right) \nonumber \\ &&+{D^2\over \beta ^2} \left(\beta _s^2-P\right)-1 =0. \end{eqnarray} \tag{ 14 }$$

In Figures 10(a)–(c), we plot the upstream electron distribution in the β_⊥ − β_∥ space of the upstream rest frame from the simulations for θ_B = 80, 82 and 835. Electrons in region I are transferred to region II. As θ_B increases, the number of electrons participating in SDA decreases since the energy threshold for SDA in Equation (7) increases.

Zoom In Zoom Out Reset image size

In Figures 10(d)–(f), we plot the normalized upstream electron energy distributions, f(E)(= (1/N)(dN/dE)), in the upstream rest frame for θ_B = 80, 82 and 835. We compare the results of Equation (13) for eΦ = 0 (dashed) and −3.5(keV) (dot-dashed) with the simulation result (solid). The dotted line is the incoming kappa distribution with T = 0.8 keV and κ = 10. Using Equation (8), the transition energy points are given by E_{trans, e} = 11.3 keV, 16.5 keV, and 24.2 keV for the angles, θ_B = 80°, 82°, and 835, respectively, when eΦ = −3.5 keV. Using Equation (9), the energy points E_{2, e}'s, beyond which the spectral index increases, are given by E_{2, e} = 34 keV, 56 keV, and 93 keV for the angles, θ_B = 80°, 82°, and 835, respectively. Here V_sh = −0.032c, r = 2.5 and m_i/m_e = 30 (M_A = 6.62 and β_p = 8.93).

In Figure 11, we plot the upstream electron energy distribution in the upstream rest frame from the out-of-plane upstream B-field simulation (solid) where ${\bf B}_1\!=\!B_{1x}\hat{x}+B_{1z}\hat{z}$ and  θ_B = 80°. The dashed and dot-dashed lines are the theoretical results for eΦ = 0 and eΦ = −3.5 keV at the shock front, respectively. In this simulation, the temperature anisotropy-driven instabilities seen in Figure 4 are precluded. By direct comparison with the electron energy distribution in Figure 10(d), the simulation result in Figure 11 is closer to the theoretical result for eΦ = −3.5 keV. That is, the high energy tail is steeper in Figure 11. This implies that the temperature anisotropy-driven instabilities near the shock front can contribute to further electron acceleration, possibly via the interaction between electrons and perturbed magnetic fields. There is a still discrepancy between the theory with eΦ = −3.5 keV and the simulation around the transition energy, E_{trans, e} = 11.3 keV. This is probably due to the effects of small-scale waves generated in the transition region as indicated by Matsukiyo et al. (2011).

Zoom In Zoom Out Reset image size

3.2.2. Bremsstrahlung Radiation

Electrons accelerated via SDA move along the magnetic field lines as can be seen from the electron tracking in Figure 7(b) and will collide with ions to emit bremsstrahlung radiation. In the low-frequency limit of ωb/γv ≪ 1, where ω is a photon angular frequency, b is the impact parameter, v is the electron speed, and $\gamma =1/\sqrt{1-(v/c)^2}$, the number of photons per unit frequency per unit volume per unit time produced by the bremsstrahlung process is (e.g., Rybicki & Lightman 1979)

$$\begin{equation} {dN\over d\omega dV dt}={16 Z^2e^6n_en_i\over 3c^3m_e^2\hbar ^2\omega } \int _{v_{\rm min}}^{c} {\gamma \over v}\ln \left({\gamma m_ev^2\over \hbar \omega }\right)f(v)dv, \end{equation} \tag{ 15 }$$

where f(v) is a normalized electron distribution and v_min is determined by the equation ℏω = (γ(v_min) − 1)mc².

Given the electron distributions in Figures 10(d)–(f), we calculate the number of photons per unit time (s) per unit energy (keV) using Equation (15) for the ion density of n_i = n_e and the volume of the region emitting the X-rays, V = 10²⁷ cm³. In Figures 12(a) and (b), we show the photon distribution for θ_B = 80 (solid), 82 (dashed), and 83.5 (dot-dashed) when the electron distribution is given by the theoretical distribution in Equation (13) with eΦ = −3.5 keV (Figure 12(a)) and by the simulation (Figure 12(b)). The dotted line is the photon distribution when the electron distribution is given by a kappa distribution with T = 0.8 keV and κ = 10. Here V_sh = −0.032c, r = 2.5 and m_i/m_e = 30 (M_A = 6.62 and β_p = 8.93).

Zoom In Zoom Out Reset image size

In Figure 12(c), we plot the averaged photon distribution of 80° ⩽ θ_B ⩽ θ_Bmax(= 8817) from the theoretical distribution in Equation (13) with eΦ = 0 (dashed) and eΦ = −3.5 keV (solid). In Figure 12(d), we plot the averaged photon distribution of θ_B = 80, 82, and 835 from the theoretical distribution in Equation (13) with eΦ = 0 (dashed) and eΦ = −3.5 keV (dot-dashed), and from the simulations (solid). In Figure 12(d), we notice that the theoretical result with eΦ = −3.5 keV (dot-dashed) is in good agreement with the simulation result (solid).

In Figure 12(d), the simulation result shows that the transition energy between the thermal and non-thermal photon spectrum is E_{trans, p} ≈ 10 keV and the energy point beyond which the photon spectrum becomes steeper is E_{2, p} ≈ 40 keV. The spectral index is δ = 3(simulation) and δ = 2.2(theory with eΦ = −3.5 keV) in 10 keV <E < 40 keV and δ = 7(simulation) in E > 40 keV. For emission from multiple shocks with different θ_B's, the transition energy E_{trans, p} and E_{2, p} would be dominated by the shock with the minimum θ_B. Note that the transition energy E_{trans, p} for the photon distribution is a bit smaller than E_{trans, e} for the electron distribution because bremsstrahlung photons are produced by electrons with higher energies than the photon energy (e.g., Holman 2003). The energy point E_{2, p} is approximately given by Equation (9).

RHESSI data for several solar flares (Table I in Altyntsev et al. 2012) show that the spectral index δ is in the range 2.5 < δ < 3 and the transition energy is in the range 12.1 keV <E_{trans, p} < 29.2 keV. Therefore, the electron energization via SDA well explains some of the RHESSI X-ray spectra for the energy regime, E < E_{2, p} and how the transition energy is related with the shock geometry, i.e., the minimum θ_B.

The observed RHESSI spectra do not show a steepening beyond E = E_{2, p}, and thus the theory herein by itself cannot account for the electron acceleration to produce those photons. This indicates that additional mechanisms, such as diffusive shock acceleration, are required for further electron energization to maintain the power-law spectrum up to E ∼ MeV. It is not unreasonable to expect that solar flares involve multiple acceleration mechanisms operating on a range of scales.

3.3. Effects of a Realistic Proton/Electron Mass Ratio on the Spectra

With the actual ion/electron mass ratio of m_i/m_e = 1836, the shock speed in the upstream rest frame is reduced by $\sqrt{30/1836}$ compared to our simulations when the Mach number M and the plasma β_p are fixed. The compression ratio r and the electric potential energy eΦ are unchanged for a fixed M and β_p (e.g., Hoshino 2001). Therefore, the shock structure are not expected to change for a realistic mass ratio simulation.

The electron energy spectrum in Equation (6) depends only on β_s, α₀, and eΦ. For a fixed M and β_p, $\alpha _0(\equiv \sin ^{-1}\sqrt{B_1/B_2})$ and eΦ are unchanged. From the definition of β_s(≡ V_sh/ccos θ_B), we only need to change the angle from θ_B into θ'_B to obtain the electron energy spectrum for the mass ratio m_i/m_e = 1836 from that with m_i/m_e = 30,

$$\begin{equation} \theta ^{\prime }_B=\cos ^{-1}\left[{\cos \theta _B\over \sqrt{1836/30}}\right]. \end{equation} \tag{ 16 }$$

For example, the averaged spectrum with 80° ⩽ θ_B ⩽ θ_Bmax(= 8817) for m_i/m_e = 30 corresponds to that with 8873 ⩽ θ_B ⩽ 8977 for m_i/m_e = 1836.

In Figure 13, we plot the averaged photon distribution of θ_Bmin ⩽ θ_B ⩽ θ_Bmax(= 8977) for a real ion/electron mass ratio, m_i/m_e = 1836, from the electron distribution given by Equation (13) with eΦ = −3.5 keV. Here V_sh = 0.0041c and r = 2.5 (M_A = 6.62 and β_p = 8.93). The different values of θ_Bmin, from θ_Bmin = 8873 to 8943, give the different transition energy points, from E_{trans, p} = 10 to 35 keV. The power indices of the photon spectrum are nearly the same as δ ∼ 2.2 in E_{trans, p} < E < E_{2, p}, where E_{2, p} runs from 40 to 150 keV.

Zoom In Zoom Out Reset image size

In summary, we studied quasi-perpendicular, low M/high β_p shocks with full PIC 2D simulations using a reduced ion/electron mass ratio m_i/m_e = 30. The shock compression ratio we found was in agreement with the Rankine–Hugoniot relation. Whistler instabilities driven by downstream temperature anisotropy were observed. A modified two-stream instability due to the incoming and reflecting ions in the shock transition region was also observed.

Abundant non-thermal electrons accelerated via SDA were observed upstream. We compared the electron energy distribution from the simulations with the distributions derived by extending a theoretical model (Mann et al. 2006, 2009; Warmuth et al. 2009) and found that they reasonably agree with each other.

In the perpendicular shocks, however, SDA can be achieved only by particles transmitting into the downstream, and their energy gains are smaller than those of the reflected ones in quasi-perpendicular shocks here (e.g., Ball & Melrose 2001). Therefore, such abundant non-thermal electrons observed in this paper were not seen in the perpendicular shocks (e.g., Park et al. 2012).

We calculated the photon flux via bremsstrahlung radiation from the electron distributions from both the theory and the simulations. We showed that a transition energy, E_{trans, p}, marking the transition from a thermal to a non-thermal part of the photon spectrum, is determined by the minimum θ_B from multiple shocks with different θ_B's in θ_Bmin ⩽ θ_B ⩽ 90°. Different solar flares have different θ_Bmin's in their termination shocks and therefore can show different transition energy points.

From the simulations, the averaged photon spectrum of θ_B = 80°, 82°, and 835 gives a spectral index δ ∼ 3 in 10 keV <E < 40 keV and the spectral index increases beyond E = 40 keV. The spectral index δ ∼ 3 as well as the transition energy, E_{trans, p} = 10 keV, well explains some of the RHESSI X-ray spectra in the energy regime, E < 40 keV. To account for the spectral index of the RHESSI X-ray spectrum up to E ∼ MeV, however, additional mechanisms other than SDA are required.

Note that although our simulations were performed using m_i/m_e = 30, we analytically scaled the results to the realistic ion/electron mass ratio, and found that the predicted photon spectra are indeed insensitive to this mass ratio.

This work was supported by NSF under Grant PHY-0903797, by DOE under Grant No. DE-FG02-06ER54879 and Cooperate Agreement No. DE-FC52-08NA28302, and by NSFC under Grant No. 11129503. We also thank the OSIRIS consortium for the use of OSIRIS. The research used resources of NERSC.

To generate a kappa distribution, we use the random number distribution by Leitner et al. (2011). We denote {N_i(0, 1)|i = 1, 2, ..., N} as a set of normally distributed random numbers with the mean of 0 and the deviation of 1, and {U_i(0, 1)|i = 1, 2, ..., N} as a set of uniformly distributed random numbers between [0,1]. Then a sequence of the random number x_i given by

$$\begin{equation} x_i=b_1 N_i(0,1)+b_2U_i(0,1)x_{i-1} \end{equation} \tag{ A1 }$$

generates a 1D kappa distribution with the mean of 0 and the deviation of 1, and the coefficients, b₁ and b₂ determine the index κ (Figure 8 in Leitner et al. 2011).

To implement the kappa generator in OSIRIS, we determine the particle's initial momentum p as

$$\begin{equation} {\bf p}_i={\bf p}_{th} x_i+{\bf p}_{d}, \end{equation} \tag{ A2 }$$

where x_i is given by Equation (A1), p_th and p_d are the thermal and the drift momentum, respectively. For κ = 10, we choose b₁ = 0.58 and b₂ = 1.15. Then a sequence of p_i in Equation (A2) generates a 3D kappa distribution in Figure 1.

The kinetic dispersion relation for electrostatic instabilities is

$$\begin{equation} 1+\sum _{s=e,i}K_s({\bf k},\omega)=0, \end{equation} \tag{ B1 }$$

where K_s(k, ω) is the susceptibility. For magnetized electrons with isotropic Maxwellian distributions, K_e(k, ω) is given by (Gary 1993)

$$\begin{equation} K_e({\bf k},\omega)={\omega _{{\rm pe}}^2\over v_{eth}^2 k^2} \left[1+\xi _e^0 e^{-\lambda _e} \sum _{m=-\infty }^{\infty }I_m(\lambda _e)Z(\xi _e^m)\right], \end{equation} \tag{ B2 }$$

where v_eth is the thermal velocity, I_m is the modified Bessel function of the first kind, λ_e = k²_⊥v²_eth/Ω_e², Ω_e = eB/(m_ec)(< 0), Z is the plasma dispersion function, and $\xi _e^m=(\omega -m\Omega _e)/(\sqrt{2}k_\parallel v_{eth})$. For unmagnetized ions with drifting Maxwellian distributions,

$$\begin{equation} K_i({\bf k},\omega)=-{\omega _{pi}^2\over 2 v_{ith}^2 k^2}Z^{\prime }(\xi _i), \end{equation} \tag{ B3 }$$

where $\xi _i=(\omega -{\bf k}\cdot {\bf V}_{id})/(\sqrt{2}k v_{ith})$ and ' is the derivative with respect to ξ_i.

Here we consider that ${\bf B}\approx B_y\hat{y}$, ${\bf k}=k_x\hat{x}$${\bf V}_{sd}=V_{sd}\hat{x}$. Then the term ξ⁰_eZ(ξ^m) in Equation (B2) becomes −ω/(ω − mΩ_e) as k_∥ goes to 0. The dispersion relation for the MTSI in the electrostatic limit becomes

$$\begin{eqnarray} 1 &+& {\omega _{{\rm pe}}^2\over v_{eth}^2 k^2} \left(1-e^{-\lambda _e}\sum _{m=-\infty }^{\infty }I_m(\lambda _e) {\omega \over \omega -m\Omega _e}\right) \nonumber\\ & -& \sum _{s=in,re}{\omega _{ps}^2\over 2 v_{sth}^2 k^2}Z^{\prime }(\xi _s)=0. \end{eqnarray} \tag{ B4 }$$