PARTICLE ACCELERATION IN SHOCK–SHOCK INTERACTION: MODEL TO DATA COMPARISON

Collisionless shock waves are a key source of energetic particles in the heliosphere and beyond (Stone & Tsurutani 1985). While most studies address particle acceleration in a single shock, shock–shock interaction is an interesting topic as well. A geometry where two shocks collide head-on establishing a contracting magnetic bottle is an evident way to produce enhanced acceleration (Fermi 1949). Such configurations can form when a coronal mass ejection (CME) driven shock hits a planetary bow shock or a corotating interaction region reverse shock. Two CMEs lifting off at different angles may also produce such a trap in the solar corona. Naturally, two astrophysical shocks can collide head-on as well.

An early observational study of the shock–shock interaction effects on ions was the work by Scholer & Ipavich (1983). Later, particle acceleration in shock–shock collisions was investigated by Cargill et al. (1986; see also Cargill 1991) using a one-dimensional hybrid simulation. More recently, Lembège et al. (2010) analyzed collisions between quasi-perpendicular non-stationary shocks with a one-dimensional full particle-in-cell simulation.

Recently, Hietala et al. (2011) reported detailed in situ observations of particle acceleration in the interaction of an interplanetary (IP) shock with the bow shock of the Earth on 1998 August 9–10. They used energetic ion data from the Advanced Composition Explorer (ACE; Stone et al. 1998), Wind (Acuña et al. 1995), and Geotail (Nishida 1994) spacecraft during a quasi-radial interplanetary magnetic field (IMF) configuration, as illustrated in Figure 1. After crossing a tangential discontinuity (TD), the spacecraft moved into a flux tube that was filled with a seed population of energetic particles accelerated by the IP shock (phase 1 in Figure 1). Since ACE was magnetically connected to the IP shock but not to the bow shock, the seed population could be characterized by its measurements: a constant spectral index and flat flux intensity profiles during the 6–7 hr from the TD crossing to the IP shock crossing. No "shock-spike" could be identified at the crossing. Wind observed several particle bursts coming from the bow shock direction during the first part of the event (phase 2). Later, Wind became continuously connected to both shocks, and measured an increasing flux as well as two counter-streaming populations until the IP shock crossing (phase 3). Geotail was located closest to the Earth and continuously connected to both shocks. It recorded the highest intensity at the IP shock crossing. Furthermore, immediately after the crossing Geotail observed a burst of very high energy particles propagating sunward (phase 4). Based on the velocity dispersion of the burst and the analysis of the geometry of the two shocks, Hietala et al. (2011) deduced that these particles had been released from the magnetic trap between the shocks as they collided.

Zoom In Zoom Out Reset image size

The aim of the present study is to use a global 2.5D test-particle version of the Corsair code (http://corsair.fmi.fi) described below to further investigate particle acceleration in this event. We concentrate on the last hour of the shock–shock interaction, when the IP shock crossed the spacecraft and acceleration phases 3 and 4 were observed. In particular, we are interested whether the shock–shock interaction in this magnetic geometry can explain the observed shapes of the flux intensity profiles. In Section 2 we briefly describe our numerical model and its application to the event of 1998 August 9–10. In Section 3 we present the simulation results and compare them with the observations. Discussion and conclusions are given in Section 4.

The general scheme of our model is similar to Figure 1. The simulation box was a quasi-two-dimensional slice with periodic Z-boundaries. For the binning of the data the box was divided into 147 × 67 cells (Δx = 2 R_E) from which the magnetosheath region was then removed. The average solar wind properties were as described in Hietala et al. (2011). In the simulations presented here the angle between the IMF and the X-axis (cone-angle) was 17°.

At the beginning of the run, the box was filled with test-particle protons in the energy range 0.1–3 MeV, 60,000 per cell. We continued to inject this amount of particles to each boundary cell at every time step of 0.25 s. We used ACE observations of the event seed population to characterize the test-particles: a power-law energy distribution with the index σ = − 4 and a step-function pitch-angle distribution. As the particle densities in the test-particle model are arbitrary, they were normalized so that the density in the 0.22–0.47 MeV range matches ACE observations in the beginning of the run.

For the proton propagation we adopted the guiding center (GC) approximation with state variables $(\mathbf {R},V_{\parallel },M)$. Here $\mathbf {R}$ is the GC position, V_∥ is its (signed) speed along $\mathbf {B}$, M = mv²_gyro/2B is the (constant) magnetic moment, and $\mathbf {v}_{\mathrm{gyro}}$ is the velocity related to gyromotion. In the solar wind the GCs were propagated using a second-order Runge–Kutta scheme:

$$\begin{eqnarray} \frac{d\mathbf {R}}{dt} & = & \mathbf {V}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{dV_{\parallel }}{dt} & = & {-}\frac{M}{m}\nabla B_{\parallel } - \frac{\mathbf {V}_{\mathrm{E}}}{B}\cdot \left[\frac{\partial }{\partial t}+(\mathbf {V}_{\parallel }+\mathbf {V}_{\mathrm{E}})\cdot \nabla \right]\mathbf {B}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathbf {V} & = & V_{\parallel }\mathbf {e}_{\parallel }+\mathbf {V}_{\mathrm{D}}, \end{eqnarray} \tag{ 3 }$$

where $\mathbf {e}_{\parallel }$ is the unit vector to the direction of $\mathbf {B}$ and $\mathbf {V}_{\mathrm{D}}$ is the instantaneous drift velocity at GC position. Here only electric drift velocity is included, i.e., $\mathbf {V}_{\mathrm{D}}=\mathbf {V}_{\mathrm{E}}=\mathbf {E}\times \mathbf {B}/B^{2}$. In Equation (2) the second term on the right-hand side contains acceleration due to non-inertial effects, i.e., due to the curvature of $\mathbf {B}$. Without this term the curvature effects are not correctly included in shock drift acceleration when the GCs are crossing the IP shock. With the definitions given above the energy of the particle can be calculated as

$$\begin{equation} U = \frac{1}{2}\,m\,V_{\parallel }^2 + \frac{1}{2}\,m\,V_{D}^2 + M\,B. \end{equation} \tag{ 4 }$$

The IP shock and Earth's bow shock were modeled using the Rankine–Hugoniot jump conditions in an infinite Mach number approximation, V_shock ≫ V_A (see, e.g., Sandroos & Vainio 2006). The IP shock was modeled as a planar shock having a constant gas compression ratio R_IP = 2.25 and velocity V_IP = 350 km s⁻¹ as described in Sandroos & Vainio (2006). However, the fields near the planar shock were modified as follows in order to correctly calculate shock drift acceleration in the GC scheme. During IP shock crossing the GCs were first transformed into IP shock normal incidence frame (SNIF), where it is easier to evaluate Equation (2). Inside the IP shock SNIF the variables only vary with distance x* to the shock. The gas compression ratio was taken to increase linearly when −L/2 ⩽ x* ⩽ +L/2,

$$\begin{equation} R_{\mathrm{IP,w}}(x^{*})=1+(R_{\mathrm{IP}}-1)\frac{x^{*}+L/2}{L}, \end{equation} \tag{ 5 }$$

where the lower index w refers to the widened shock front. The magnetic field was modified accordingly,

$$\begin{equation} \mathbf {B}_{\mathrm{w}}(x^{*})=\mathbf {B}_{\mathrm{n,1}}+R_{\mathrm{IP,w}}(x^{*})\mathbf {B}_{\mathrm{t,1}}, \end{equation} \tag{ 6 }$$

where $\mathbf {B}_{1}$ is the (unmodified) magnetic field at the shock and the indices n and t refer to normal and tangential components. The parallel acceleration (2) can readily be calculated using Equations (5) and (6) as

$$\begin{eqnarray} \frac{dV_{\parallel }}{dt} & = & -\frac{M}{m}\frac{{B}_{\mathrm{w},x}}{B_{\mathrm{w}}}\frac{B_{\mathrm{w},y}CB_{\mathrm{1},y}+B_{\mathrm{w},z}CB_{\mathrm{1},z}}{B_{\mathrm{w}}}\nonumber \\ & & +\frac{C}{B_{\mathrm{w}}}(V_{\parallel ,x}+V_{\mathrm{E},x})\cdot (V_{\mathrm{E},y}B_{\mathrm{w},y}+V_{\mathrm{E},z}B_{\mathrm{w},z}), \end{eqnarray} \tag{ 7 }$$

where C = R'_{IP, w}(x*) = (R_IP − 1)/L. In Equation (7) we have assumed that $\partial _{x*}\mathbf {B}_{\mathrm{t,}1}\ll \mathbf {B}_{\mathrm{t,}1}\partial _{x*}R_{\mathrm{IP,w}}$, which is the case for constant $\mathbf {B}$. In the simulations presented here we used L = 10 km.

The shape of the bow shock was calculated using the empirical model by Merka et al. (2005). If particles were transmitted through the bow shock, i.e., particles were able to enter the magnetosheath, they were removed from the simulation. Thus we did not have to use a physical model for the bulk plasma in the magnetosheath. This also allowed us to use more physical models for the bow shock gas compression ratio than would have been otherwise possible.

When the GCs crossed Earth's bow shock from upstream to downstream, they were transformed into a local bow shock de Hoffman–Teller (HT) frame, and either a reflection or transmission operator was applied based on the pitch-angle cosine μ in the HT frame. The local transformation velocity $\mathbf {V}_{\mathrm{T}}=\mathbf {V}_{\mathrm{SNIF}}+\mathbf {V}_{\mathrm{HT}}$, where

$$\begin{eqnarray} \mathbf {V}_{\mathrm{SNIF}} & = & \mathbf {V}_{\mathrm{sw,t}}, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \mathbf {V}_{\mathrm{HT}} & = & {-}|\mathbf {V}_{\mathrm{sw,n}}||\tan \theta |\mathbf {B}_{\mathrm{t}}/|\mathbf {B}_{\mathrm{t}}|. \end{eqnarray} \tag{ 9 }$$

Here $\mathbf {V}_{\mathrm{SNIF}}$ is the transformation velocity from Earth's rest frame into the local bow shock SNIF, $\mathbf {V}_{\mathrm{HT}}$ is the transformation velocity between the local SNIF and HT frames, and tan θ = B_t/B_n is the tangent of the local shock obliquity angle in the upstream region. If the GC pitch angle was outside the loss cone $\cos \alpha _{\mathrm{HT}}<\cos \alpha _{\mathrm{LC}}=\sqrt{1-B_{2}/B_{1}}$, the GC was reflected back to the upstream, V_{HT, ∥} → V'_{HT, ∥} = −V_{HT, ∥}.

The gas compression ratio R_BS which influences the downstream magnetic field $\mathbf {B}_{2} = \mathbf {B}_{\mathrm{1n}} + R_{\mathrm{BS}}(r)\mathbf {B}_{\mathrm{1t}}$ was taken to depend on the position along the bow shock surface. In Figure 2, the red line shows the solution of the cubic equation for the compression ratio of an oblique MHD shock (Priest 1984, p. 202), where we have used purely radial solar wind velocity of 400 km s⁻¹ and β₁ = 1.4 measured by Wind together with the model for the bow shock shape. The black line depicts the estimated gas dynamic compression ratio. It was calculated using the scaling scheme:

$$\begin{eqnarray} M_{\mathrm{MS}}^2 & = & \frac{2 R_{\mathrm{nose}}n_{x}^2}{(\gamma +1)-(\gamma -1)R_{\mathrm{nose}}} \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} R_{\mathrm{BS}} & = & \frac{(\gamma +1)M_{\mathrm{MS}}^2}{2+(\gamma -1)M_{\mathrm{MS}}^2}, \end{eqnarray} \tag{ 11 }$$

where M_MS is the magnetosonic Mach number, R_nose = 3.78 is the compression ration at the bow shock subsolar point, n_x is the x-component of the bow shock normal vector, and γ = 5/3. When Equation (11) is multiplied by n^1/4.5_x, we obtain a good match with the MHD solution, as indicated by the blue dashed line in Figure 2. For the results presented here we have used this modified gas dynamic model as it is computationally less demanding than the numerical solution of the cubic equation.

Zoom In Zoom Out Reset image size

3.1. Spatial and Temporal Evolution

The animation available in the online version of this journal and the snapshots in Figure 3 show the evolution of particle density in the energy range 0.22–0.47 MeV during the simulation run. The locations of Wind and Geotail are indicated as well.

Zoom In Zoom Out Reset image size

During the first minutes of the run, the magnetic field lines connected to the dusk flank of the bow shock are filled with particles reflected and accelerated by the quasi-perpendicular part of the curved shock. Since there is no scattering in the model, the parallel part of the bow shock does not accelerate particles but lets them through. Thus the field lines connected to the pre-noon sector contain less particles than expected for the real foreshock region. However, in the present study we are not interested in this region as the spacecraft were located on the dusk side.

Figures 3(a) and (b) illustrate another transient feature of the simulation: the front formed by the first particles that have reflected from the IP shock proceeds through the box. The particles can be seen to bounce between the shocks gaining in energy while drifting to the −Y direction.

After these transients there is a steady temporal increase in the test-particle density between the two shocks (Figure 3(c)). On the other hand, the IP shock passage causes a spatial effect: the density of energetic particles is higher on its downstream side due to adiabatic heating. Furthermore, the orientation of the magnetic field lines containing the energetic particles changes across the shock. Thus when the IP shock crosses a particular point, e.g., the location of Wind, the field line passing through that point scans the bow shock dusk flank (Figure 3(d)). As the connection point moves to the +Y direction, the field line may eventually disconnect from the bow shock like it does for Wind (Figure 3(e)).

Around 2645 s (44 minutes) the shocks begin to pass through each other starting from the post-noon region. This results in a drop in the test-particle density on the magnetic field lines connected to this region, as there is no longer a magnetic trap between the two shocks. The animation and Figure 3(f) illustrate how the drop propagates downstream along the field lines at the speed of the last escaping particles.

The animation shows a stream of particles from the bow shock far flank into the region between Geotail and Wind until the IP shock reaches the top corner of the simulation box. This stream is probably not physical, or at least did not take place in the observed event, as it is related mainly to two factors of the simulation geometry. First, the bow shock is not compressed by the higher magnetic field and dynamic pressure of the IP shock downstream region. Second, in the particular run the downstream magnetic field direction is almost tangential to the model bow shock surface, while the observed cone-angle was larger. In addition, the simulated far flank compression ratio may be too large.

3.2. Comparison of Time Profiles with Observations

To compare the results of the model with the observations, we extracted time series of the magnetic field and the test-particle density in different energy channels from several virtual spacecraft placed in the simulation box. The virtual spacecraft placed on the field lines that were connected only to the IP shock record flat profiles throughout the run (not shown). The flatness is natural, as there is no scattering in model. Moreover, the profiles are like those observed by ACE that had a similar magnetic connection (Figure 4 in Hietala et al. 2011).

Figure 4 shows the simulated time series from locations corresponding to those of Wind and Geotail next to their observations. At the beginning of the run, both virtual spacecraft (Figures 4(a) and (c)) record two step-like increases related to the initialization: first the formation of the foreshock of the bow shock and then the arrival of the first particles that have reflected from the IP shock (at 7–8 minutes for Wind and at 10–12 minutes for Geotail). After these transients the virtual spacecraft record smoother increase until the IP shock crossing. Note that in the model the magnitude of the downstream magnetic field is approximately 2 nT lower than observed due to small changes in the IMF direction right before the shock crossing.

Zoom In Zoom Out Reset image size

Comparison of the shock crossing in Figures 4(a) and (b) shows that the overall profile is well produced by the model for Wind's location. For the last ∼5 minutes before the shock crossing, the densities at 0.47–0.62 MeV and 0.62–1.4 MeV channels match the observed ones quite well also quantitatively. The apparent density jump from upstream to downstream in the simulated profile is due to the scatter-free assumption of the model. The after-shock density decrease has a duration of 8–9 minutes similar to the observed decrease of ∼7 minutes. In addition, the bumpy structures in the decrease visible in the two lowest channels seem to match well.

At the location of Geotail (Figures 4(c) and (d)), the simulated densities are slightly lower than observed during the last 5–10 minutes before the shock crossing, but not significantly. The maximum density is reached after the IP shock crossing both in the model and in the observations. Although the peak produced by the model is smaller than observed, their timing—4 minutes after the shock in the simulation and ∼2 minutes in the data—is very similar.

Hietala et al. (2011) noted that Geotail's observations suggest a small velocity dispersion for the peak: the higher energy particles seem to peak before the lower energy particles. However, this dispersion was at the limit of the instrument's time resolution (48 s and 98 s at the energies shown in Figure 4(d)). From the results of the numerical model (Figure 3) we see that the energetic particle density near the bow shock nose suddenly drops as the shocks collide and the trap between them collapses. At the location similar to Geotail's, this drop ends the post-shock increase and results in the peak of the profile (Figure 4(c)). Furthermore, there is a 30–50 s velocity dispersion in the arrival time of the density drop to the virtual spacecraft, as the last high-energy particles reach it before the last low-energy particles.

We have tested the sensitivity of our results to the choice of the bow shock shape (different paraboloids, model by Merka et al. 2005), bow shock compression ratio (constant, linear, gas dynamic, modified gas dynamic), and the IMF cone-angle (22° and 17°). We find that profiles close to the bow shock are quite sensitive to bow shock parameters. Furthermore, the IMF cone-angle has a significant effect on the profiles recorded by both virtual spacecraft, especially on the rise before the IP shock crossing. According to the observations, the cone-angle increased slowly by 10°–15°during the last 1.5 hr before the shock crossing: at Wind from 10° to 23° and at Geotail from 18° to 27°. Moreover, there was a more rapid change during the 5 minutes preceding the shock crossing (Figures 4(b) and (d)), so that the observed θ^IP_Bn ∼ 55°–60°, while in the simulation with the 17° cone-angle θ^IP_Bn was 46°. Due to this change, the analytically calculated downstream magnetic field direction and magnitude differ from the observations. The downstream magnitude affects the energization of the particles, and the smaller compression may be one of the reasons why the simulated post-shock density peak is somewhat smaller than the observed peak.

The evident limitation of our model is the neglect of the third dimension: the particles trapped between the shocks drift to the Z direction but the bow shock does not curve away. Besides, the particles are reflected from the bow shock without drifting along its surface. However, the gradient drifts at the shocks would be in approximately opposite vertical directions and partly cancel out in consecutive bounces. Moreover, the scatter-free shock acceleration increases particles' parallel energy only (Sandroos & Vainio 2006). Hence, the accelerating particles will eventually enter the loss cone of one of the shocks. Another significant loss mechanism is the electric drift to the −Y direction and then escape through the parallel bow shock. Due to these factors the energy gain in the 2.5D model is not unrealistically large.

In the present study we have compared the 2.5D simulation results against the observations and found a very good agreement. The model produces the observed pre-shock increase in the energetic particle density at the locations of both Wind and Geotail. Based on the simulations it is evident that confinement by the IP shock is needed to get the high densities and energies. The simulations also exhibit the post-shock maximum at Geotail's location at approximately right time and with a velocity dispersion compatible with the observations. The model does not produce quite as high densities as were observed, but we have only simulated the end of a long process. The densities are probably also decreased by the lower-than-observed downstream magnetic field and the fact that the model does not contain any scattering, which would certainly be important at Geotail's location inside the turbulent foreshock.

In conclusion, we find that shock–shock interaction in this magnetic geometry can produce the main features of the observed profiles and, in particular, the post-shock maximum at Geotail's location. The next step in the analysis of this event would be to implement a three-dimensional version of the Corsair code. It is also possible to incorporate to the model MHD background fields as well as wave-observation-based mean free paths for scattering. A quantitative analysis of the possible particle energization in the actual shock merging, though intriguing, would require a self-consistent, microphysical plasma simulation.

We thank the CDAWeb at NSSDC for data access and the PIs for the data used in this study. The work of H.H. was funded by the Finnish Graduate School in Astronomy and Space Physics and of A.S. by the Academy of Finland grant 251797.