The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas

There is growing evidence that magnetic reconnection and turbulence are intimately connected in magnetized plasmas (e.g., Refs. 1–5 and references therein). However, the link between shocks and these two processes has been less clear. We explore this connection using global kinetic simulations of the Earth's bow shock, and in the process address a number of puzzling observations. Due to the availability of in situ measurements,⁶ the Earth's bow shock is the most studied example of collisionless shocks. The bow shock is a curved shock formed due to interaction of the supermagnetosonic solar wind with the Earth's magnetic dipole which acts as a barrier. The shocked solar wind is called the magnetosheath and is bounded by the magnetopause which diverts the flow.

Observations, along with theory and simulations, have provided considerable insight into the physics of collisionless fast shocks. The curved nature of the bow shock means that at a given instant, there are regions with quasi-parallel $Q∥$ (⁠$0°≤θBN≤45°$ and quasi-perpendicular $Q⊥ (45°≤θBN≤90°)$ geometries. Here, $θBN$ is the angle between the shock normal and the incident magnetic field. In the quasi-parallel case, there exists a significant number of ions reflected from the shock that reach to large distances upstream. The relative streaming between these ions and the incoming plasma excites instabilities which generate low frequency waves that can grow to large amplitudes, giving rise to steepened fronts called magnetosonic shocklets and short large-amplitude magnetic structures (SLAMS) as well as other nonlinear structures such as cavitons (see review in Ref. 6). These structures are convected back into the shock and contribute to generation of turbulence in the magnetosheath.

There are also reflected ions associated with quasi-perpendicular shocks, but these ions remain confined to the foot of the shock, leading to shock reformation on the gyroscales of the ions for a certain range of parameters (e.g., Ref. 7). There can also be turbulence downstream of quasi-perpendicular shocks due to temperature anisotropy driven instabilities with magnetic fluctuation levels $Brms/Bo∼0.1$ as compared to $Brms/Bo∼1$ downstream of quasi-parallel shocks.⁸ 

Despite this progress, a number of important issues remain. One issue of considerable debate is whether bow shock turbulence has significant large scale consequences, especially since the turbulence is initiated on ion scales. In particular, observations have shown evidence of anomalous flows in the magnetosheath. Examples include observations of fast jets (see Refs. 9, 10, and references therein) and sunward flows (see Ref. 11 and references therein). The origin of such flows and their impact on the magnetosphere remain controversial. Since solar wind has embedded discontinuities that frequently impinge on the magnetosphere, models to explain these anomalous flows fall into two general categories: those that invoke the role of discontinuities and those that rely on the physics of the bow shock itself. Similarly, observations have reported reconnecting current sheets in the quasi-parallel magnetosheath.^12,13 These are difficult measurements and confirmation through simulations would be helpful but no evidence for such current sheets have yet been reported in the global simulations.

Simulations provide a global view of the magnetosphere and its dynamical evolution. In order to address the above issues, and to gain a better understanding of the turbulence associated with the bow shock, we use global hybrid and global fully kinetic simulations that capture the self-consistent formation of the bow shock and its associated turbulence. Most studies to date have focused on the upstream region of the $Q∥$ bow shock, called the ion foreshock. We use our simulations to compare and contrast the nature of fluctuations/turbulence in the ion foreshock, $Q∥$ magnetosheath, and $Q⊥$ magnetosheath and their large scale consequences. We have run the simulations to a factor of ∼9 times longer than most previous studies in order to gain insight into the evolution of turbulence in time. Through these simulations, we find that the $Q∥$ magnetosheath turbulence is much stronger than either that in the foreshock or in the $Q⊥$ magnetosheath. We demonstrate, for the first time, that this enhanced level of turbulence leads to generation of reconnecting current sheets and formation of magnetic islands. This finding establishes the link between shocks, turbulence, and reconnection. This result is consistent with observations of reconnection in the quasi-parallel magnetosheath,^12,13 where, however, no search for magnetic islands was performed. In 3D, magnetic islands correspond to flux ropes. In the remainder of this paper, we will use the term magnetic islands and flux ropes interchangeably. Since some of the flux ropes reach tens of ion inertial length, we expect a dedicated observational search for flux ropes in the magnetosheath to yield many examples. One possibility to search for such flux ropes in the data may be through the use of the Grad-Shafranov reconstruction technique.^14,15

The simulations also enable us to address large scale consequences of turbulence in the magnetosphere. We examined the possibility that foreshock turbulence, rather than interaction of solar discontinuities with the magnetosphere, could cause anomalous flows. We show that turbulence causes corrugation of the shock surface which in turn gives rise to collimated, high velocity plasma jets inside the magnetosheath. This is in agreement with the model originally proposed by Hietala et al.⁹ and the properties of the jets are found to be consistent with high velocity magnetosheath flows recently reported in the observations.^9,10 The jets have high dynamical pressure and we show that upon reaching the magnetopause they can at times trigger so-called flux transfer events.¹⁶ 

Next, we considered whether the jets, which penetrate into the magnetosheath, can trigger sunward flows. Such a connection was proposed by Shue et al.¹¹ who attributed the origin of the sunward flow to the interaction of a jet with the magnetopause, which in turn would cause the magnetopause to rebound, creating a sunward motion. However, we find sunward flows in the vicinity of many jets that terminated in the magnetosheath far from the magnetopause. A careful examination of the time history of the flow reveals a complex process that would have been difficult to uncover in observations. In this process, the jets “stir” the magnetosheath, creating large scale disturbances and pushing the ambient low density plasma and the embedded magnetic field. The plasma in the sunward flow regions is not the plasma in the jets that was deflected, rather it is the plasma accelerated by thermal and dynamical pressure that develop because of the jets. In other words, the energy of the jets is convected into sunward flow rather than the deflection of the jets themselves. This energy conversion can give rise to sunward flows that extend tens and even hundreds of ion inertial length into the solar wind, well past the nominal position of the bow shock. The strong turbulence in the magnetosheath also gives rise to regions where the magnetic field during anti-sunward interplanetary magnetic field (IMF) orientation can point sunward in the magnetosheath. This is another testable prediction for observations.

Another consequence of the jets is found to be the nonlinear formation of large scale vortices in the magnetosheath. These vortices are driven due to turbulence and are to be distinguished from the usual Kelvin-Helmholtz driven vortices along the flanks of the magnetopause (Refs. 17, 18, and references therein).

Although these results were obtained from simulations of the bow shock, which is curved, similar effects are also expected at planar shocks. For example, the current sheets formed downstream of the quasi-parallel shocks is due to the strength of the turbulence, rather than any physics tied to the curvature of the shock. Similarly, the rippling of the surface of the shock is due to interaction of the shock with the steepened fronts getting convected back into it, an effect that is also present in planar shocks. The paper concludes with the discussion of the implications of these findings for the general understanding of plasmas and draws the link between shocks, reconnection and turbulence.

Our hybrid simulation code is H3D,^19–21 a massively parallel hybrid code with a stretched mesh capability. Although we have conducted 3D global simulations of the magnetosphere,^20,21 here we will focus on 2D simulations for three reasons: (a) One of the key issues considered is whether the shock generated turbulence can lead to formation of reconnecting current sheets. The identification of reconnection in 3D remains a subject of active research and is still not completely resolved. (b) Confirmation of the formation of reconnection sites observed in our global hybrid simulations through equivalent fully kinetic global simulations can only be done in 2D since 3D global fully kinetic simulations will remain out of reach in the foreseeable future. (c) We have conducted an extensive parameter study to identify the most interesting regimes which serve as conditions for 3D simulations to be published elsewhere.

Distances are normalized to ion skin depth d_i, time is normalized to ion cyclotron frequency $Ωci$⁠, magnetic field is normalized to upstream magnetic field B_o, density is normalized to upstream density N_o, and velocity is normalized to the Alfvén speed $VA=Bo/(4πNomi)1/2$⁠, where m_i is the ion mass. The ratio of plasma to ion cyclotron frequency is set as 6000 and the ions and electrons are taken to have the same temperature T_i = T_e for these runs. The Alfvénic Mach number is given by M_A = V_sw/V_A, where V_sw is the solar wind speed. We use a polytropic electron model with γ_e of 5/3 and 7/6, with the latter being suggested as more appropriate for shock physics.^22,23 The de Hoffmann-Teller cross shock potential under polytropic electron closure is proportional to $γe/(γe−1)[kTe]$⁠. Thus using the more nearly correct γ_e for solar wind thermal electrons^22,23 of 7/6 leads to a larger de Hoffmann-Teller potential in hybrid simulations, raising the potential coefficient to 7 versus 5/2 for γ_e of 5/3.

We have conducted large number of hybrid simulations covering a range of D_p of 40 to the Earth's size of 600 d_i, plasma beta of 1 to 4, Alfvénic Mach number M_A of 4 to 12, IMF angle of 0° to 90°, and IMF B_z/B_o of 0 to 1. This range of parameters is not meant to be exhaustive but does represent a good coverage of the relevant parameter space. Here, D_p is the nose position of the magnetopause normalized to ion inertial length, which has been shown to be a good way to characterize the type of magnetospheric structure as a function of dipole strength.²⁴ We have also examined the effects of (i) spatial resolution by varying the grid from 0.25d_i to 1d_i, (ii) particle per cell by considering a range of 10 to 200, (iii) and resistivity by considering different models such as no imposed resistivity, and several models based on gradients of magnetic field. We have found no significant effect of the resistivity on the physical processes considered here. However, care is required in general since use of large uniform resistivity as is done in some studies can lead to excess heating of the magnetosheath. We will focus this paper around the runs listed in Table I that are from converged solutions.

The simulation setup is similar to that used in global MHD simulations and consists of a magnetic dipole embedded in a uniform solar wind plasma which is continually injected from the left boundary. Mirror dipole is used to ensure the magnetic field is just the IMF field at the injection boundary. This reduces magnetic field perturbations at the injection boundary. All other boundaries use open boundary conditions. As in global MHD simulations, an inner boundary is placed around the Earth to limit the Alfvén speed from reaching relativistic values. Two particle boundary conditions are implemented: reflecting and absorbing. The dipole can also be ramped up in time, although we have not found significant differences compared to the case where the dipole is initialized at full strength at the start of the simulation. As the solar wind plasma interacts with the dipole, the magnetosphere forms in time. In 2D, the dipole is a line dipole in order to have a divergence free field.

Figure 1 shows the intensity plot of magnetic field at three different times for run 1 where the IMF direction changes at the injection boundary from −45° to 10° at $Ωcit=90$⁠. This launches a rotational discontinuity (RD) that propagates towards the magnetosphere. The RD eventually runs into the reflected ions in the ion foreshock before reaching the bow shock. The time evolution of this interaction is shown in Figure 2 where we show a segment of the simulation zoomed in around the RD-reflected ion interaction region, along with 1D cuts of several parameters at y = 3100d_i. Figure 2 shows evidence for particles following the field lines across the RD into the upstream region. An observational evidence for this effect would be particles on magnetic field lines that do not point at the bow shock. Figure 2 also shows formation of at least three nonlinear structures: a fast shock, a cavity with hot plasma with depleted plasma density, and an anisotropic RD. The RD is seen to separate in time from the shockfront as expected.

The interaction of an RD with shock reflected ions was shown to lead to a compound structure referred to as the foreshock bubble (FB) which consists of a fast magnetosonic shock attached to a very hot, low density cavity behind it.^25,26 The solution found here differs from the example of FB in Ref. 25 in at least two important aspects. First, the structure in Fig. 2 is more consistent with a localized front such as a shocklet than a shock since the plasma state behind it goes back to the solar wind conditions as evidenced from the 1D cut of plasma variables (e.g., Figures 2(d) and 2(j)). Second, the cavity is seen to remain much smaller than the extent of the shockfront in y and does not grow much in time (Figures 2(a)–2(c)). In contrast, in the case of FB in Ref. 25, the cavity covers almost the entire extent of the area behind the shock. Given these differences, we have labeled this structure as “FB” to distinguish it from the previously reported FBs.

As it turns out, there is a spectrum of solutions between the example found here and that reported in Refs. 25 and 26. The resulting compound structure can take on several forms depending on the thickness of the RD, the IMF direction, among others. Particularly important is the level of density depression in the cavity formed near the shockfront. The larger the density depression, the hotter the temperature within the cavity. In moderate density depressions such as those in Figure 2 $(>0.1No)$⁠, the cavity does not increase in size significantly. Note, however, that the size of the shockfront here in both x and y directions can still reach several hundred ion inertial lengths as seen in Figure 2. For more extreme density depressions approaching $∼0.05No$⁠, the shockfront and the cavity become one structure and can grow to very large scales even reaching the size of the magnetosphere if given sufficient time. However, one has to be cautious regarding simulations of such structures in the hybrid code. Standard hybrid codes are not ideally suited for treatment of structures when the density depression reaches the density floor set in the simulations. And we have found a tendency to see weaker cavities with use of larger number of particles per cell and/or higher resolutions. We will examine the formation of foreshock bubble using fully kinetic simulations elsewhere. An observational signature of such a structure, assuming it can be realized in nature, would be a large region of the solar wind populated with very low density and hot plasma.

As is evident from Figure 1(a), the quasi-parallel magnetosheath appears more turbulent than the quasi-perpendicular magnetosheath. And after the IMF change, the new quasi-parallel magnetosheath in Figures 1(a)–1(c) appears more turbulent than the new quasi-perpendicular magnetosheath. This is confirmed by calculation of the magnetic field fluctuations which are about a factor of 4 larger in the $Q∥$ regions. This higher turbulence level is also confirmed by additional diagnostics in Figure 3. Figures 3(a) and 3(b) show the magnetic field spectrum and the kurtosis (fourth standardized moment) of magnetic field increments $δB(x,s)=B(x+s)−B(x)$ as a function of separation length s for the foreshock, $Q∥$ and $Q⊥$⁠, magnetosheath respectively. Kurtosis, or flatness, is a measure of deviations from Gaussian. The flatness of Gaussian is exactly three. The $Q∥$ magnetosheath shows clear turbulent features, i.e., an energy spectrum with close to −5/3 range and strong non-Gaussianity (i.e., a very large kurtosis which increases as the length scale decreases). Interestingly enough, the kurtosis values are in good agreement with in situ measurements of $Q∥$ magnetosheath in cases where evidence of turbulent reconnection was reported.¹³ In contrast, the kurtosis of the foreshock and the $Q⊥$ magnetosheath are much smaller. In particular, the kurtosis of $Q⊥$ is almost flat with values between 4 to 5, showing no or almost no signature of intermittency. Intermittency is a signature of the generation of coherent structures. As we will demonstrate shortly, turbulence generates coherent structures in the form of reconnecting current sheets and magnetic islands in the $Q∥$ but not in the $Q⊥$ magnetosheath, consistent with the kurtosis diagnostics. Figure 3(c) shows the comparison of the ion energy distribution function in the foreshock, $Q∥$⁠, and $Q⊥$⁠. The foreshock shows the presence of a counter streaming population as expected. The magnetosheath shows power-law distribution functions but with different power-law indices for the $Q∥$ and the $Q⊥$ cases. The thermal energy and the high energy particle production are significantly higher in the $Q∥$ magnetosheath. This is not surprising since higher levels of turbulence in the $Q∥$ region can lead to additional heating and acceleration of particles. Another mechanism for acceleration to high energies is the Fermi process which requires the particles to go back and forth across the shock several times. Since particles can move along the magnetic field more readily than across it, and because magnetic field lines move along the bow shock as the result of their advection with the solar wind, one might expect more energetic particles in the $Q∥$ magnetosheath since field lines there have been connected to the bow shock the longest. The relative importance between these competing processes is beyond the scope of the present work and will be explored elsewhere.

The power-law indices found here are not universal and can vary under different magnetosheath conditions. Previous 1D hybrid simulations of quasi-parallel shocks have also shown power-laws over a similar range^27,28 but they show flattening of the particle spectrum at higher energies where we do not have sufficient statistics. This flattening of the particle spectrum is also seen in 2D hybrid simulations of perpendicular shocks that include pre-existing, large-scale magnetic-field turbulence in the upstream region.²⁹ The details of the turbulence are clearly very different in our simulations as compared to 1D simulations where reconnection and magnetic islands are excluded. Direct comparison of particle energization in local shock simulations versus our global simulations will be presented elsewhere.

In the present study, the initial solar wind ions are taken to have a Maxwellian distribution function. However, nonthermal ion distribution functions are commonly observed in the solar wind. The presence of nonthermals can affect the instabilities in the foreshock which in turn can affect the details of turbulence and associated particle energization. This is an important extension of the present work that we are currently working on.

Next, we compare the evolution of the turbulence level as measured by the root mean square of the magnetic field fluctuation amplitude $Brms=(∑i=1n(Bi−⟨B⟩)2/n)1/2$ which is normalized to B_o. Here, ⟨B⟩ is the average magnetic field in the region where B_rms is calculated. We use run 2 where the IMF direction remains at 10° during the entire length of the simulation. We have run the simulation to longer times than previous hybrid simulations so that we can explore the asymptotic behavior of the turbulence. The result is shown in Figure 4(a). Although the foreshock turbulence is a major source of turbulence downstream of the quasi-parallel shocks, the correlation between the turbulence level in these two regions is more complex as evident in Figure 4(a). While on average the increase or decrease in the foreshock activity leads to a corresponding increase or decrease in the turbulence levels downstream, there is a temporal lag between changes in the foreshock and its transmission to magnetosheath. For example, there is no discernible increase in the foreshock turbulence levels from time of 200 to 500, whereas the $Q∥$ turbulence level has increased by a factor of nearly 2.4. However, in time the foreshock turbulence level starts to drop off, and the turbulence level downstream of the shock also starts to decrease.

As another way to see this, and also to get a handle on the timescale for decay of turbulence levels once the geometry changes from quasi-parallel to quasi-perpendicular (Fig. 4(b)), we examine time evolution of B_rms for run 3. This run, which is a higher resolution (0.5d_i) version of run 1, has identical parameters as run 2 except that the IMF direction is changed from −45° to 10° after $Ωcit=90$⁠. We used this run as part of our convergence study to ensure the results are not affected by cell resolution or number of particles. It is seen that the change to quasi-perpendicular geometry results in a rapid drop in the turbulence fluctuation levels to about $Brms∼0.3$ which is similar to that in the $Q⊥$ magnetosheath in Figure 4(a). It is interesting to note that the mean rate of solar wind IMF directional discontinuities is about 0.5−2.5 h (Ref. 30) or about $Ωcit∼150−750$⁠. This is competitive with decay times of turbulence in the magnetosheath and provides a natural low frequency “stirring” mechanism for the magnetosheath turbulence.

The significantly enhanced level of turbulence in the $Q∥$ magnetosheath has implications for the nature of turbulence as we now demonstrate.

In 2D, the presence of magnetic islands formed through reconnection can be demonstrated by plotting the contours of the out-of-the-plane vector potential or tracing of magnetic field lines. An effective way to visualize field lines is through the line-integral-convolution (LIC) technique (see the Appendix) which we first introduced in plasma physics in our study of shear driven turbulence.³ Figure 5 shows the magnetic field lines colored by V_x for a zoomed in area of the quasi-parallel shock for run 1 at $Ωci=694$⁠. The magnetic field intensity plot of a larger region of this simulation was shown in Figure 1(c). Figure 5 clearly shows a wide range of magnetic islands up to $∼30di$ in diameter. Some of the islands show evidence of outflowing jets (not shown).

These magnetic islands are observed in the quasi-parallel but not in the quasi-perpendicular magnetosheath (not shown). The absence of magnetic islands in the $Q⊥$ magnetosheath is the reason for differences in the kurtosis in the $Q∥$ and $Q⊥$ regions shown earlier (Fig. 3(b)). Also note that the islands can exist even in the vicinity of the $Q∥$ region of the bow shock. Due to the small size of these magnetic islands/flux ropes, they are not expected to be force free.

These results are consistent with Cluster observations of magnetosheath reconnection events^12,13 which also were only present in the quasi-parallel magnetosheath. Although postulated, no search for mangetic islands was conducted in these observational studies. Detection of small magnetic islands is challenging but given that some of the islands grow to as large as $∼30−40di$⁠, it may be possible to find such islands in the data.

We have found the reconnecting current sheets and magnetic islands to be a common feature of quasi-parallel shocks. Since the size of the magnetosheath is larger for larger D_p, the number of islands and their interactions increases at larger D_p. Figure 6 shows the density and the LIC of magnetic field colored by the strength of B for run 4 (D_p = 300). The magnetosheath is found to be almost completely volume filled with magnetic islands in this case.

The enhanced turbulence level in $Q∥$ magnetosheath also leads to stronger heating as compared to $Q⊥$ magnetosheath as evident in Figure 7. We have found evidence for this in the observations. Figure 8 shows the statistical map of THEMIS measurements between October 2007 and October 2013 when the probes occupied the magnetosheath. The quantity binned here is the mean value of the total ion temperatures during a Parker-Spiral IMF (⁠$Bx>0.2|B|$ and $By<0.2|B|$ or $Bx<0.2|B|$ and $By>0.2|B|$⁠), which is presented in the Magnetosheath InterPlanetary Medium (MIPM) reference frame. The MIPM frame uses upstream measurements to arrange within a normalised model magnetosheath thus accounting for boundary motion and planetary aberration. For a detailed description of the MIPM frame and subsequent data processing, we refer the reader to access the methodology paper which also compares the statistical mapping results against global MHD simulations produced by the CCMC BATS-R-US code.³¹ To briefly summarize, the MIPM frame arranges points as a function of the upstream IMF. Therefore, statistically it can be closely compared to an aberrated GSE system during parker-spiral IMF orientation where Y_MIPM > 0 represents a quasi-perpendicular shock and Y_MIPM < 0 quasi-parallel. The physical dimension of each bin is 0.5 × 0.5 Earth radii and the bin number density is typically a hundred data points per bin. Each data point is calculated from the mean averaged 3-min intervals of THEMIS data. Currently, the data are binned using a 20 min sliding averaged OMNI data around each 3 min THEMIS data interval. The points during which IMF cone angle varies more than 20° are removed. It is clear that the dayside MSH is hotter due to the shock compression and the hot plasma region extends more tail-ward at the quasi-parallel shock side of the magnetosheath, possibly due to kinetic processes. The quantitative calculation as a function of MIPM-frame zenith angle shows that the magnetosheath downstream of the quasi-parallel bow shock is on average about 10%–15% hotter. A more detailed temperature study using THEMIS data for various solar wind conditions and quantitative comparisons between dawn and dusk sectors will be discussed elsewhere.

Hybrid codes do not include electron kinetic effects. As such, magnetic reconnection in hybrid simulations occurs due to either numerical or imposed resistivity. To confirm that the magnetosheath reconnection events observed in our hybrid simulations are not due to numerical effects, we have also conducted a 2D global fully kinetic simulation for the radial IMF condition. The detailed comparison of this run with an equivalent hybrid simulation will be presented elsewhere. Here, our interest is only to establish that reconnection and magnetic islands also occur in the global fully kinetic simulation.

The simulation was performed using a general-purpose plasma simulation code VPIC.³² The simulation is of size 2000 × 1500d_i with 56,544 × 41,984 cells. The IMF magnetic field is in the x direction and the dipole is oriented at an angle of 80° with respect to the IMF. The incoming solar wind plasma has parameters M_A = 8, T_i = T_e, $β=8π(Te+Ti)no/Bo2=1$⁠, $ωpe/ωce=3$⁠. The solar wind is normal to the dipole. The ion-to-electron mass ratio is m_i/m_e= 50. The strength of the dipole corresponds to stand-off distance D_p = 50d_i. Reflecting boundary conditions for particles are prescribed on a circle of radius R = 0.3 D_p around the geometrical center of the dipole. The dipole is created by prescribing out-of-plane current in two circular regions of the simulation domain or radius 0.1R separated by distance 0.25R. The initial conditions correspond to a uniform plasma with solar wind parameters. In order to minimize the perturbations created by the initial conditions, the strength of the currents creating the dipole is ramped up in time as $j=j0[1−exp(−t/τ)]$⁠, where $τ=30Ωci−1$⁠.

The result is illustrated in Figure 9 where only a segment of the simulation around the quasi-parallel region is shown. The LIC of magnetic field is shown colored by ion density. As in the hybrid example, many magnetic islands are observed in the quasi-parallel magnetosheath. Our closer examination of the magnetic field structures (not shown) shows the hierarchy of islands extending down to the electron kinetic scales. However, the largest islands have similar size $(∼30di)$ as in the hybrid simulations.

It is well known that the relative streaming of ions reflected from the bow shock and the incoming solar wind leads to the generation of large amplitude ultralow frequency (ULF) waves.⁶ As these waves get convected back towards the shock, they can steepen, giving rise to shocklets and SLAMS. Here, we examine the evolution of such structures, their interaction with the bow shock, and their role in magnetosheath turbulence. For simplicity, here we refer to the steepened structures as wavefronts or shocklets. As seen from Figure 10, the wavefront-magnetosphere interaction falls into four stages. The IMF is nearly radial in this case (run 5). The part of the wavefront hitting the bow shock first gets compressed to shorter scales (Figs. 10(a), 10(c), and 10(d)), causes a rippling of the shockfront as well as local spike in the temperature (Fig. 10(b)). As the wavefront gets through the shock (Fig. 10(e)), its amplitude is amplified (Fig. 10(g)). The waves persist as wavefronts in the magnetosheath for tens of ion cyclotron times (Figs. 10(i)–10(l)) but then they start to break up and dissipate their energy. In this way, turbulence upstream gets amplified and contributes to the turbulence in the magnetosheath.

Interestingly enough, large amplitude waves have also been observed downstream of the quasi-parallel shocks at Venus. A recent study compared detailed wave diagnostics in both upstream and downstream regions of the Venusian bow shock and concluded that the downstream waves are transmitted ULF waves.³³ 

In cases where the steepened fronts are nearly parallel to the quasi-parallel region of the bow shock surface, they can significantly modify the shock structure and the magnetosheath thickness. This is shown in Fig. 11, from run 6, which shows the time evolution of density and B_y. The propagation of one wavefront, labeled as WF1, is marked in time. The steepened front is evident in Figures 11(a) and 11(e). Note that the steepened front has many properties similar to a shock, exhibiting localized density and magnetic pulses/enhancements over the background values. By the time $Ωcit=47.2$⁠, the front has gained strength (Figures 11(b) and 11(f)) and creates a large density jump ahead of the shock (Figure 11(b)). The thickness of the magnetosheath is much thinner in this region as compared to the upper region away from the shocklets. A rather large magnetic island of about 10d_i in width is evident as a circular loop in the density panels. WF1 grows in scale (Figures 11(c) and 11(g)) and eventually gets convected into the shock and moved downstream (Figures 11(d) and 11(h)). This process can repeat for several cycles. This is evident in Figures 11(d) and 11(h) which show the formation of a second wavefront (WF2) which will go through the same process as WF1. Figure 12 shows a blow up of B_y and an associated 1D cut across the wavefront at $Ωcit=47.2$⁠. The 1D cut shows the amplification of the wavefront amplitude as it crosses the shock into the magnetosheath.

Figures 13(a) and 13(b) illustrate the dynamics between the shocklet and the bow shock. The LIC of magnetic field colored by density and LIC of ion velocity colored by V_x are shown for a zoomed in region of Figure 11(c). The shocklet presents an obstacle to the flow, generating local spikes in the density (Fig. 13(a)) and the magnetic field, and decelerating the flow (Fig. 13(b)). The region between the shocklet and the magnetosheath can be quite turbulent. Several vortices and magnetic islands are clearly evident in this in-between region, similar to the magnetosheath. Note also the presence of a relatively large magnetic island spanning about 10d_i along the magnetopause surface. The vortices are apparent here since the flow speed is nearly zero in the shown regions. The crossing of such shocklet/bow shock structures may be mistaken as multiple shock crossings in the data, especially since there are regions with conditions similar to the solar wind in between them.

Another effect of ion foreshock that is perhaps not widely appreciated is that it reduces the magnetosonic Mach number as evident from the 1D cut in Figure 14. As such, the ion foreshock helps in the shock dissipation process and affecting the slowing of the plasma. The presence of the jets can also been as streaks of high magnetosonic Mach number within the magnetosheath (Figure 14). The modification of M_s due to the ion foreshock is consistent with an observational study where significant evolution of interplanetary shock parameters was observed during its motion along the bow shock.³⁴ 

In this section, we demonstrate the effects of turbulence on the magnetosheath flow. There are a number of diagnostics that are useful for visualizing flow. These include streamlines, pathlines, and streaklines. Streamlines represent instantaneous patterns of the flow. A true vortical pattern on scales larger than kinetic scales would wrap the magnetic field. However, the presence of such a vortex can be masked, using streamlines, in the presence of a large background flow as is the case in the magnetosheath. Instead, we use streaklines to study the mechanisms for formation of sunward flow and vortices in the magnetosheath. Streaklines provide a time history of the flow and can reveal the presence of vortical pattern even in case of large background flow. Each streakline is created by tagging and following the trace of fluid elements passing through a particular spatial point. Streaklines can be thought of the pattern traced by dye if one were to steadily inject dye into the fluid at a fixed point.

An interesting feature evident in Figure 5 is the presence of two ion jets, shown as red, that penetrate deep inside the magnetosheath. The magnetic field lines are seen to remain nearly radial as in the solar wind along these jets, embedded in an otherwise sea of magnetic islands. In this section, we examine these jets in more detail. Figure 15 shows plots of the dynamical pressure over a larger region of the simulation than in Figure 5. Presence of many jets is evident, seen as streaks of high dynamical pressure, with some reaching the magnetopause. As the jets slow down in the magnetosheath, they form a bow wave (Fig. 15) which can be as large as several hundred d_i. The interaction of the jets/bow wave with the magnetopause causes strong perturbations, giving rise to surface waves, and/or inner-magnetospheric compressional waves.^11,35,36 And in this case, this interaction has triggered flux transfer events (Fig. 15(b)), although the IMF is 10°. The jets also lead to significant inhomogeneities in the magnetosheath and serve as an additional driver of turbulence. Figure 15(b) shows the LIC of B for an area zoomed in around the large bow wave in Fig. 15(a). This leads to another source of reconnection and generation of magnetic islands in the magnetosheath.

In time, the jets get convected towards the tail while new jets are formed. Figure 16 shows the properties of the jets, which strongly deviate from the surrounding magnetosheath. In particular, jets exhibit enhanced magnetic field strength (Fig. 16(b)), lower temperature (Fig. 16(c)), and higher density (Fig. 16(d)) than the surrounding magnetosheath. Note also the strong kinking of the magnetic field (Fig. 16(b)) along the jets. In each case, we have found the wavelength of the kinked jet to be in good agreement with the wavelength of Kelvin-Helmholtz instability based on local measurements of velocity shear and the thickness of the shear layer.

Observations have revealed the presence of transient enhancements of plasma flow within the magnetosheath characterized by various measures such as high dynamic pressure and kinetic energy density.^{9,10,35,37–40} Hietala et al.⁹ reported high velocity jets during nearly radial conditions (IMF angle less than 20°) for M_A = 12. They found the plasma density to be compressed along the jets with density depletions around them, the velocity of the jets to be close to the upstream value, and the width of them on the order of 50 to 100 ion d_i. These features are fully consistent with the observed jets here. The velocities of the jets shown in Figure 5 are close to the solar wind speed of 8V_A, their widths are in the range of 20 – 60d_i, and there is density compression on the order of a factor of 6 along the jets (Fig. 16(d)) and large density depletions around them, with densities dropping below the solar wind density in some cases. The size and the level of the depletion area vary on time. The ion temperature (Fig. 16(c)) shows that the temperature does not change significantly along the jets. This is also consistent with properties of the observed jets.^9,10 This may also explain previous observations⁴¹ which had found cases where the solar wind passes through the shock layer without significant heating.

The origin of these jets has been controversial. Hietala et al.⁹ attributed their formation to ripples along the shock surface. Since the shock mainly decelerates the component of the upstream velocity parallel to the shock normal, ripples along the shock cause local changes to the angle between these two vectors. For large angles, the shock speed would remain close to its upstream value. While others have suggested the origin of the jets to be due to interaction of the magnetosphere with solar wind discontinuities (e.g., Ref. 42) and others (e.g., Ref. 39) questioned the ripple-based mechanism. In contrast, a more recent study¹⁰ using 4 years of subsolar magnetosheath observations concluded that 97% of the observed jets are consistent with the ripple idea. This finding is consistent with our simulation results where we find the bow shock surface to be highly rippled due to foreshock turbulence, and we find the jets even in cases where we keep the IMF direction fixed in time and there are no solar wind discontinuities.

In the absence of turbulence, the downstream flow must be directed away from the bow shock at all points as seen in global MHD simulations. However, THEMIS observations during radial IMF conditions have shown evidence of sunward flows. Such flows are also seen in our simulations. Figure 17 shows V_ix for the same time in the simulation as in Figure 16. The purple areas are regions in the magnetosheath with negative V_ix with peaks around 2V_A. Note that there are several regions along the bow shock where such sunward flows are present.

Several different mechanisms have been proposed to explain this anomalous flow. Shue et al.¹¹ attributed the origin of the sunward flow to the interaction of a fast anti-sunward flow of unknown origin with the magnetopause, which in turn causes the magnetopause to rebound, creating a sunward motion. As shown in Figure 16, the generation of fast anti-sunward flow is due to the rippling of the bow shock. Although some of the jets interact with the magnetopause, many terminate at significant differences away from the magnetopause (e.g., Figure 16(a)) and yet there are sunward flows in the vicinity of such jets (Fig. 17). This indicates that alternative mechanisms must be in play for the generation of sunward flow than that proposed by Shue et al.¹¹ The simulations enable us to study the generation mechanism of sunward flow in much more detail than is possible in spacecraft measurements.

We have seen that the jets snake around the highly inhomogeneous magnetosheath, some reaching the magnetopause while others are terminated earlier. As we will see, some of the jets encounter local magnetic field that is nearly perpendicular to them. This causes them to get deflected since it is difficult for plasma to traverse perpendicular to the magnetic field. This complex dynamics has a direct bearing on the formation of sunward flow. We use streaklines to follow the evolution of the flow. First, we examined the role of jets that reach the magnetopause in causing the sunward flow. Figure 18(a) shows the streaklines shortly after they were seeded in two regions in the foreshock. These two groups of seeds were carefully selected so as to follow both the plasma motion inside the two large jets that reach the magnetopause (painted in white) as well as the plasma motion around the jets and their bow waves (painted in blue). Several points are immediately clear in Figure 18. First, the sunward flow is not caused by reflection of the jets. Second, there is no evidence that magnetopause motion is creating a large scale sunward flow. Third, in Figure 18(d), one of the largest jets has started to wrap although it has not formed a complete vortex. Any wrapping clearly created a sunward flow locally. Fourth, it is evident in Figures 18(b)–18(d) that the large bow wave pushes on the plasma around it, seen as blue streaks, as it propagates downward. This is another way that sunward flows can be generated.

Next, we seeded the field lines to examine the origin of the sunward flows between the jets that extend into the solar wind as were seen in Figure 17. The evolution of the streaklines is shown in Figure 19. Here again we marked the seeds in two colors, those that interact with the two jets marked as J₁ and J₂, marked in white and those that are from the areas in the vicinity but outside of these jets, marked in blue. The streaklines at the mouth of the opening in the region between J₁ and J₂ are seen to bend sunward as the blue streaklines penetrate this area (Figs. 19(c)–19(f)). The blue streaklines represent plasma that was pushed downward due to the pressure of the jets on the ambient plasma. In this way, the jet energy causes deflection of the ambient sunward, rather than the plasma within the jets being deflected.

As another check on the formation mechanism of sunward flow, we show in Figure 20 a zoomed-in view of the region containing jets marked as J₁ and J₂ in Figures 19(d) and 19(e). The top two panels show the LIC of magnetic field colored by ion temperature. The jets are seen to have a temperature similar to the solar wind whereas the area between the jets, marked as R2, is much hotter. The magnetic field lines clearly show that the surface of the bow shock is highly corrugated. A few sample magnetic field lines are plotted in white, while density contours are plotted in black. Note the presence of a band of nearly laminar magnetic field that is nearly perpendicular to the jets as they enter the magnetosheath, marked as R1. As the nearly parallel propagating jets encounter this band of magnetic field, they cannot penetrate further into the magnetosheath as it is difficult for plasma to move across the magnetic field. The upper jet is seen to snake around in the magnetosheath and gets deflected when reaching this band of magnetic field. The lower jet also winds around but does not cross this band of magnetic field. The plasma in the region between the two jets R1 gets squeezed and follows the magnetic field lines downward. This plasma pushes outward the shocked solar wind between the two jets in region R2, causing sunward flow as evident from negative V_ix flows seen in Figures 20(c) and 20(d). This picture is consistent with the development of the flow seen in Figure 19. The blue streaklines in Figure 19 represent the flow that follows the band of magnetic field lines R1 and enters the region between the two jets R2, creating sunward flow.

These considerations have led us to the following picture for formation of sunward flows. The jets, which have high dynamical pressure, “stir” the magnetosheath, creating large scale disturbances and pushing the ambient lower density plasma and the embedded magnetic field. The plasma in the sunward flow regions is not the plasma in the jets that was deflected, rather it is the plasma that was already there that is accelerated by thermal and dynamical pressure that develops because of the jets. In other words, the energy of the jets is deflected into sunward flow rather than the deflection of the jets themselves. This energy deflection can give rise to sunward flow that extends tens and even hundreds of ion inertial length into the solar wind, well past the nominal position of the bow shock (Fig. 17).

The presence of jets and sunward flows in the $Q∥$ magnetosheath can drive velocity shears which can lead to formation of large vortices. Figure 21 (run 7) shows the formation of a large scale vortex which occurs in two stages. A structure akin to the Rayleigh-Taylor finger is formed (Fig. 21(a)) which then gets caught into the velocity shear in the magnetosheath (Fig. 21(b)), forming a large scale vortex (20c). The magnetic field lines are wrapped around, similar to a vortex generated by Kelvin-Helmholtz (KH) instability. Figure 22 illustrates the generation mechanism of the vortex using streaklines. A band of streaklines above the nominal stagnation line flow is seen to get deflected sharply and turn downward (pink) (Fig. 22(a)). This is due to the deflection of the energy of a jet which pushes the ambient plasma downward (not shown) similar to that for run 1 in Figure 19. However, unlike in Figure 19, here this plasma gets caught in the velocity shear in the magnetosheath (Fig. 22(b)) and turns into a vortex (Fig. 22(c)). The color coding of streakline is also useful for tracking the movement of the plasma in different regions of the magnetosheath as seen in Fig. 22(c). The structure in Figure 21(a) has the appearance of a Rayleigh-Taylor finger. The plasma in the jet is denser than the ambient plasma which in the presence of an accelerating force that would mimic gravity could give rise to the Rayleigh-Taylor instability. However, here the formation of this structure appears to be due to dynamics, i.e., deflection of the jet energy as it encounters a nearly perpendicular magnetic field. Although KH driven vortices are known to form on the flanks of the magnetopause (e.g., Ref. 18 and references therein) due to the velocity shear between the solar wind and the magnetopause, we believe this is the first report of a vortex driven due to turbulence in the magnetosheath.

The fundamental processes of magnetic reconnection, shocks, and turbulence have been demonstrated to be more intimately connected than previously realized. Laminar reconnecting current sheets can develop turbulence;⁴³ studies of turbulence show that formation of reconnection current sheets may be a generic feature of turbulence in magnetized plasmas (e.g., Refs. 1–5 and references therein). The link between shocks and the other two processes has, until this study, been less clear. For example, one early model of reconnection, the Petschek model, required four standing slow shocks as part of the reconnection exhaust.⁴⁴ However, the applicability of this model in collisionless plasmas remains in question.⁴⁵ Turbulence does not lead to shock formation, except possibly in shear driven turbulence,³ but even there turbulence is not the cause of shock formation.

Shocks of all types in collisionless plasmas share the property of generating reflected disturbances of electrons and especially ions that can get back upstream into the unshocked medium.^6,46 The angle of upstream magnetic field to the shock normal guides the reflected particles and is one of several parameters controlling the spatial volume to which they may reach before being swept back through the shock layer. In addition, the shock layer substructure may provide for the instability of the transmitted distributions, as is particularly well known for quasi-perpendicular shocks. The turbulence on the downstream side of the shocks reflects the geometry of B to the normal, with the quasi-parallel shocks possessing downstream fluctuations comparable to the mean field and the quasi-perpendicular downstream regions quieter. We have demonstrated for the first time that turbulence associated with the $Q∥$ shocks lead to generation of magnetic islands, with some islands even forming at the shock surface. In some cases (e.g., Fig. 13), the islands can even form ahead of the shock. We have also demonstrated that reflected upstream disturbances when convected back and transmitted through the shock layers eventually dissipate there. Thus, both “waves” and reconnection play a role in providing dissipation downstream of the quasi-parallel shocks. The question of the relative importance of these two processes for dissipation is beyond the scope of this paper and will be addressed elsewhere.

We have thus forged a link between shocks, turbulence, and magnetic reconnection in collisionless plasmas. In doing so, we have provided a theoretical argument that explains the location of reported reconnection events in the magnetosheath.^12,13 Previous studies have not reported the incidence of reconnection in simulations of fast shocks or global hybrid simulations possibly because (i) many fast shock studies were in 1-D where reconnection cannot occur; (ii) the very limited number of 2-D quasi-parallel shock simulations apparently did not check for incidence of reconnection; and (iii) the low number of particles per cell and or low grid resolution of these early simulations (e.g., Ref. 25) may have masked the detection of magnetic islands. Our findings here show that with increasing turbulence levels, the system exhibits greater structural organization in the form of coherent structures such as flux ropes, vortices, and fast jets. Whether such reorganization into coherent structures is a consequence of a general governing principle of all systems far from equilibrium that are driven remains an intriguing speculation.

Another interesting result here is in regards to the large scale consequences of the bow shock turbulence in the magnetosphere. This turbulence can, under certain conditions, enable formation of large scale jets that can interact with the magnetopause and lead to surface waves on it, or inner magnetospheric compressional waves.^11,35,36 The jets can also lead to formation of sunward flows and large scale vortices downstream of the shock.

Given that we found magnetic islands to be a common feature of quasi-parallel shocks with $Brms∼1$⁠, search for such structures in the data would be of interest. Ordinarily, flux transfer events are conceptually associated with southward IMF; our simulations have shown examples where flux transfer events can occur with northward IMF along the magnetopause (Fig. 15(b)) as well as inside the magnetosheath (e.g., Fig. 6(b)). While data identified flux transfer events can reach scales of hundreds of d_i, our results demonstrate that flux ropes due to shock “turbulence” appear to be smaller and at most reach tens of d_i in scale.

We have already started the extension of the present work to 3D. Our preliminary results indicate that the basic physical effects uncovered here, such as reconnection in the magnetosheath, remain intact in 3D. However, 3D introduces additional complexity such as fanning of the magnetic field lines in the ion foreshock and formation of Kelvin-Helmholtz instability along the flanks. These results will be reported elsewhere.

This work was partially supported by NASA (NNX12AD30G, Heliophysics Theory), NSF (0904734, 1104815), and DOE (DE-SC0004662). M.W. acknowledges support from NSF Grants No. AGS-1063439 and No. AGS-1156094 (SHINE), and from NASA Grant No. NNX11AJ44G. JS acknowledges support from NSF ATM 0802380 and NASA NNX07AF40G. Simulations were performed on Pleiades provided by NASA's HEC Program, and Blue Waters sustained-petascale computing project, which is supported by the NSF (OCI 07-25070) and the state of Illinois. We gratefully acknowledge useful discussions with J. Giacalone, J. Borovsky, and A. Retino. We also thank the referee for making useful suggestions that improved the paper.

The standard way to visualize streamlines of a vector field is to seed some points and integrate to trace curves that are instantaneously tangent to the velocity vector. Although useful, this approach can be cumbersome when interactively exploring a dataset. For example, the technique is inherently local and unless feature locations are know a priori it is difficult to select an effective set of seed points. An alternative technique called “line integral convolution” or LIC, which convolves noise with a vector field producing streaking patterns that follow vector field tangents, has the advantage that a very detailed view of the streamlines over the entire computational domain is found in one step without the need to explicitly specify a set of seed points.

The technique was initially developed for use on images,⁴⁷ but has since been extended to arbitrary surfaces,⁴⁸ and to volumes.⁴⁹ The algorithm has been ported to the GPU⁵⁰ and data-parallel implementations have been developed.⁵¹ We have extended this technique and have made it available through the latest release of ParaView,⁵² an open source, cross platform, tool for parallel interactive visualization of large datasets. In addition to traditional visualization algorithms, ParaView provides an MPI based data-parallel GP-GPU surface LIC algorithm which includes a number of enhancements to the basic LIC algorithm designed for interactive data exploration. ParaView's surface LIC can be used on massively parallel supercomputers without GPUs. It is described in detail in ParaView's online documentation.⁵³ To aid the reader in the interpretation of the figures presented here, a brief review of the LIC algorithm is presented here.

The LIC algorithm of a vector field defined on an image, I(x,y) is given by the following integral over streamline arcs computed from the center of each pixel location, x, y.

where L is the integration length, N generates the noise value at a given location, S_i is a position on a streamline arc centered at x, y, k(i) is an appropriate convolution kernel, and I(x,y) is the image pixel at x, y. In practice, streamline arcs are computed using an RK method over a fixed number of steps, and L can be a constant for all pixels in I, or it may be a function of the local vector field. When L is constant, the resulting visualization has a uniform look. When L varies as a function of the vector local field, the visualization accurately shows relative strengths of flow features. Because the LIC produces a dense representation of the flow field, features of interest can be quickly identified. LIC is often used in conjunction with scalar pseudocoloring. In this case, a specialized shader is used to combine the colors with the gray scale LIC image. Optimized implementations, especially data-parallel or GPU based, are very fast making the technique useful for interactive data exploration.

An example of image LIC of electron stream lines colored by the magnitude of electron velocity is shown in Figure 23(a). The data are from a 2D fully kinetic simulation of asymmetric current sheet. The formation of a vortex associated with a magnetic island (not shown) is clearly evident. The surface LIC technique we use is similar to image based LIC except first vectors are projected onto the surface then into image space where an image LIC algorithm is used to compute LIC.⁴⁸ Lit, pseudocolored surface geometry, is combined with the image LIC to produce a realistic rendering of the surface geometry. To show the utility of surface LIC in analysis of 3D data, we have applied it to data from a 3D fully kinetic simulation of reconnection.⁴³ Figure 23(b), using isocontour of density, shows the presence of multiple flux ropes as evident from their helical magnetic field lines.