Electron-Scale Reconnection in Three-Dimensional Shock Turbulence.

Magnetic reconnection has been observed in the transition region of quasi-parallel shocks. In this work, the particle-in-cell method is used to simulate three-dimensional reconnection in a quasi-parallel shock. The shock transition region is turbulent, leading to the formation of reconnecting current sheets with various orientations. Two reconnection sites with weak and strong guide fields are studied, and it is shown that reconnection is fast and transient. Reconnection sites are characterized using diagnostics including electron flows and magnetic flux transport. In contrast to two-dimensional simulations, weak guide field reconnection is realized. Furthermore, the current sheets in these events form in a direction almost perpendicular to those found in two-dimensional simulations, where the reconnection geometry is constrained.

Introduction

Magnetic reconnection is a process in which magnetic field lines in a plasma change their topology, often accompanied by the conversion of stored magnetic energy to kinetic energy of accelerated particles (Priest & Forbes, 2000; Yamada, 2011). Reconnection plays an important role in laboratory and space plasma processes, including sawtooth crashes in tokamaks (von Goeler et al., 1974), magnetic substorms in the Earth's magnetosphere (Angelopoulos et al., 2008) and solar flares (Sweet, 1969).

A recent development in the study of reconnection has been the observation of purely electron‐scale reconnection regions, in which ions do not participate in the reconnection process (Gingell et al., 2019, 2020, 2021; Phan et al., 2018; Stawarz et al., 2022; Wang et al., 2019). This differs from the standard picture of collisionless reconnection (Ishizawa et al., 2004; Malyshkin, 2008; Yamada, 2011), where the reconnection region consists of an electron‐scale layer within a wider ion‐scale diffusion region. These electron‐only reconnection regions have been observed in the Earth's magnetosheath and shock transition region (Gingell et al., 2019, 2020, 2021; Phan et al., 2018; Stawarz et al., 2022; Wang et al., 2019), foreshock (Liu et al., 2020; Wang et al., 2020), magnetotail (S. Lu et al., 2020) and laboratory experiments (Shi et al., 2022).

In this work we focus on reconnection in the turbulent environment downstream and in the transition region of quasi‐parallel shocks. Numerous electron‐scale current sheets are observed in these regions, leading to a favorable environment for electron‐scale reconnection (Gingell et al., 2017, 2019, 2020, 2021; Phan et al., 2018; Sharma Pyakurel et al., 2019; Stawarz et al., 2022; Wang et al., 2019). While existing kinetic simulation studies in the shock environment have shown reconnection at both electron and ion scales, the simulations have been two‐dimensional (Bessho et al., 2019; Bessho et al., 2020, 2022; Q. Lu et al., 2021), restricting reconnection to a single plane. This is usually the plane containing the upstream flow and magnetic field vectors. These studies find strong guide‐field reconnection, in which there is a magnetic field component parallel to the current. In contrast, observations show a wide range of guide fields (Gingell et al., 2020), with statistics showing that weak guide‐field events are slightly less favored, but still observed. Three‐dimensional effects have also been shown to enhance the reconnection rate in laminar electron‐scale current sheets (Pyakurel et al., 2021). It is therefore important to study shock‐driven reconnection in three dimensions.

We perform fully‐kinetic, three‐dimensional simulations of a high Mach number quasi‐parallel shock. We find numerous reconnecting and non‐reconnecting current sheets, and demonstrate reconnection in two events with strong and weak guide fields.

Simulation

We perform three‐dimensional simulations of a quasi‐parallel shock using the fully‐kinetic particle‐in‐cell code VPIC (Bowers, Albright, Bergen, et al., 2008; Bowers, Albright, Yin, et al., 2008). The initial condition consists of a uniform plasma and electromagnetic fields, with B  = B 0 cosθ, B  = B 0 sinθ, and E  = V flow B 0 sinθ. The initial plasma moves in the negative x direction with velocity −V flow. The lower x boundary uses conducting walls for fields and reflecting walls for particles, while plasma and the z‐component of the electric field are injected at the upper x boundary with the initial field and flow values. The y and z boundaries are periodic. The simulation domain is 2000 × 500 × 200 (d e)3 covered by 4,000 × 1,000 × 400 cells, and is initialized with 150 particles per species per cell (note that typical 2D simulations take place in what we define as the x‐y plane). Physical parameters used in the simulation are ω pe/Ω ce = 4, m /m  = 100, , θ = 30° and M  = 10. Here ω pe is the electron plasma frequency, Ωce the electron cyclotron frequency, β the ratio between thermal pressure and magnetic pressure for either species and M  = V flow/v the Alfvén Mach number of the injected plasma. As the simulation develops, the shock front propagates from the lower x boundary in the positive x direction. Unless otherwise mentioned in the text, length scales in the paper are normalized to d , and velocities to c, and number densities to the initial upstream density.

Results

The shock propagates in the positive x direction at 1.9 v , giving rise to an Alfvén Mach number of 11.9 in the laboratory frame. An overview of the shock transition region at tΩci = 17.5 is shown in Figure 1. In the quasi‐parallel shock geometry, the interaction between incident and reflected particles in the foreshock leads to strong electromagnetic fluctuations and plasma turbulence, resulting in the formation of numerous current sheets in the transition region and downstream of the shock. Examples of current sheets can be found in Figure 1a, in which regions of enhanced positive and negative J can be seen. As can be seen in Figure 1a, where there is structure in both the x‐y and x‐z planes, the three‐dimensional geometry allows current sheets, and hence reconnection, to have different orientations.

Figure 1

(a) Three‐dimensional overview of the shock transition region at tΩci = 17.5. The upstream region is on the right (+x), and plasma flows in the −x direction, while the shock propagates in the +x direction. Multiple regions with intense J are visible in the x‐y and x‐z planes. (b–d) Evolution of a current sheet. Colored panels show the current density in the out‐of‐plane direction, while the arrows point out the structure that evolves into the reconnecting current sheet. The red, yellow (out‐of‐plane) and green lines on the figure show the L, M and N directions where L is the direction of the reconnecting field, M is in the direction of the current and N is normal to the current sheet. These lines have length 10 d . The simulation axes are shown by a glyph at the bottom.

Similar to prior two‐dimensional simulations of quasi‐parallel shocks (Bessho et al., 2019, 2020), both reconnecting and non‐reconnecting current sheets are found in the simulation. Here we focus on two active reconnection sites, the first of which is shown in Figures 1b–1d and 2. In Figure 1, we show the evolution of the current sheet as it develops in the turbulent plasma. The panels show the out‐of‐plane current density in a local coordinate system constructed using the method of Denton et al. (2018), which will be described in more detail later. Unlike two‐dimensional simulations where the current sheet forms along the z (out‐of‐plane) direction, this current sheet is primarily oriented along the y direction, as shown by the orientation of the planes in (b–d) relative to the simulation axis glyph. At tΩci = 17, there is a region with small positive J just below and to the left of the intersection of the axes marked by the arrow. The current density intensifies at tΩci = 17.25 as it moves upward, before becoming its most intense in the final panel at tΩci = 17.5. We also note there is some out‐of‐plane motion of the current sheet that is not captured by this figure, but the overall increase of the intensity of J in this current sheet holds within the volume around this plane. This process takes place over 0.5/Ωci using the initial conditions, or approximately 2.5/Ωci,local using the local magnetic field, meaning the event is transient.

Figure 2

Structure of the reconnection region from various angles at tΩci = 17.5. (a) Magnetic field line configuration with colors showing the different topologies associated with reconnection. Black and yellow are associated with the inflow, while pink and green are associated with the outflow regions. (b) Electron velocity u in the L‐N plane. The red, yellow and green lines in the figure show the L, M and N directions respectively. L and N are labeled and M is approximately out of plane. The L, M and N lines have length 10 d . (c) Three‐dimensional volume rendering of current density J showing the reconnection region and field lines. (d) Energy conversion from fields to particles E′ ⋅ J. (e and f) Magnetic flux transport (MFT) (see text) diagnostic showing active reconnection. The simulation axes are shown by the glyphs, with (c) having a different orientation from the other five plots.

In order to characterize the reconnection region, we find a local coordinate system for the current sheet using a hybrid minimum variance/maximum directional derivative method (Denton et al., 2018). In this system, L is the direction of the reconnecting magnetic field, N is normal to the current sheet and M completes the orthogonal triad. For this event, , and when evaluated at tΩci = 17.5. Due to the motion of the current sheet, the LMN coordinate system does not remain the same with time, though the directions remain similar. For instance, after accounting for the motion of the current sheet, at tΩci = 17.25, which shows the largest change, is rotated approximately 18° from at tΩci = 17.5. The coordinate system is illustrated in panels (b) and (c) of Figure 2, in which the L, M and N axes are red, yellow and green respectively, and the original Cartesian system is shown by the glyphs. In terms of the simulation coordinates, the current is mainly flowing in the y direction, in contrast to two‐dimensional simulations where the reconnection current is restricted to flow in the z direction (Bessho et al., 2019, 2020).

The structure of the reconnection region and its signatures can be seen in Figure 2, which shows two‐dimensional plots in the L‐N plane and three‐dimensional views from different angles. The colored lines in Figures 2a–2c are magnetic field lines traced from points close to the current sheet, and are associated with four distinct regions. The black and yellow lines are associated with the reconnection inflows, while the pink and green lines are associated with the outflows. Panel (b) shows the electron outflow velocity in the simulation frame, and it can be seen that there is a flow reversal close to the region where the four different types of field lines meet, indicative of a reconnection region. Panel (c) shows a volume rendering of the enhancement of current density in the M direction, with values of J below 0.28 made transparent. This figure illustrates that the current sheet has a finite extent in the M direction, and we find that its width at the half maximum of J is approximately 14 d .

Further evidence of active reconnection is shown in Figures 2d–2f. Panel (d) shows E′ ⋅ J where E′ = E + v × B, indicating the conversion of energy from fields to particles which can occur in reconnection regions. In panels (e) and (f), we use a recently developed diagnostic based on magnetic flux transport (MFT), which has been used to diagnose reconnection sites in turbulent regions in both simulations and observations even when flow patterns are unclear (Li et al., 2021; Qi et al., 2022). For this reconnection site, the MFT velocity is calculated as where E is the electric field in the M direction and B is the in‐plane (L‐N plane) magnetic field. This quantity can be applied to 3D simulations (Li et al., 2021) and shows the transport of in‐plane magnetic flux. In Figures 2e and 2f, the converging inflows of magnetic flux can be seen in the U signature, while the U signature exhibits diverging outflows. We also note that the guide field for this event is approximately −0.12, approximately 9% when compared to the upstream magnetic field in the L direction.

To calculate the reconnection rate, we study the asymmetry of the reconnection region. Figure 3c shows the density, magnetic field and current density along a cut in the N direction across the current sheet indicated by the white line in panel (b). The current sheet is asymmetric, with B 1 = 1.8, n 1 = 4.3, B 2 = 0.95 and n 2 = 4.3, where n and B are measured at the location where B is maximized on either side of the current sheet. The density in the center of the current sheet is higher, and we note that there is density variation in the L direction as well, ranging from approximately three on the left outflow to 7.5 on the right outflow. The electron and ion L (outflow) velocities are plotted in Figure 3d, where the primed quantities indicate the velocities are measured in the frame moving with the B , B reversals. Due to the curvature of the current sheet, the outflows are not exactly aligned, as indicated by the lines labeled L 1 and L 2 in Figure 3b. As seen in Figure 3d, the L 2 outflow with negative is stronger. The participation of the ions in this reconnection event is an interesting question. As seen in Figures 3a and 3c, while bulk ion flow in the L direction is negative in the moving frame, there is a (white) region of less negative ion velocity just above the current sheet. However, this is associated with a region of enhanced ion density outside the current sheet away from the electron flow reversal, and there is no strong ion flow along the L 2 cut corresponding to the larger electron outflow.

Figure 3

(a) Ion velocity in the frame moving with the reconnection site. (b) Current density in the L‐N plane. The white lines show the cuts along which quantities are evaluated in (c and d). (c) Density, magnetic field and current density along the N direction. (d) Ion and electron velocities along the L 1 (top left) and L 2 (bottom right) cuts.

From asymmetric reconnection scaling laws, the theoretical electron outflow velocity is given by , where ρ is the electron mass density rather than the ion density used for standard asymmetric reconnection (Cassak & Shay, 2007). Using upstream parameters, the theoretical outflow velocity is 0.64, whereas the peak outflow velocity in this plane is approximately 0.22. If we modify the prediction to account for the pressure variation in the outflow direction (Murphy et al., 2010), the predicted velocity is reduced to 0.53. The reconnection electric field in the region is approximately E  = 0.1, so that the reconnection rate calculated using the theoretical outflow speed E /B v theory = 0.15 where B  = 2B 1 B 2/(B 1 + B 2). When using the measured outflow speed instead, the rate is 0.36. These values are reasonable when compared with reconnection in two‐dimensional shock simulations (Bessho et al., 2019, 2020, 2022).

To illustrate the different configurations of reconnection, we show a second event in Figure 4. In this case, there is an intermediate B of approximately 0.7 compared to the nominal upstream magnetic field B , indicating that this region exhibits guide‐field reconnection. The L, M, and N axes are illustrated in panels (a) and (b), and the axis directions remain similar over the preceding 0.5/Ωci, similar to the first event. Because of the presence of the guide field, the change in field‐line topology is not as drastic as in the previous event, as all the field lines are traced out in the M direction. Instead, one can observe the magnetic shear from one side of the sheet to the other, and the spreading of field lines traced from different regions in the vicinity of the B ‐B reversal. The current sheet is extended in the M direction with a length of approximately 30 d as shown in panel (a). In the L‐N plane, the current density J is shown in panel (b), where the black line is a contour of J  = 0.5J to serve as a guide in subsequent panels. Signatures of reconnection can be seen in panels (c) and (d), where the U shows inflows and outflows of magnetic flux. The structure of the electron flows is shown in panel (e), and there are oppositely directed flows in the top‐right and bottom‐left of the region inside the contour. Again, the jets are not collinear or symmetric, with the negative jet reaching a peak of approximately and the positive jet only reaching in the x‐line frame. Similar to the previous event, the peak velocity is lower than the theoretical velocity for this event, which is approximately 0.6. As seen in panel (f), there are negative L‐directed ion flows in the region in and around the current sheet, with the flows most negative in the bottom‐left, and least negative toward the top‐right, indicatating they may be associated with the reconnection event, though the region is too small for ions to couple fully to the reconnection process (Sharma Pyakurel et al., 2019). These ion flows also do not show the same small‐scale structure as the electron flows. This event shows some electron heating, as seen in panel (g) where there is an increase in electron temperature toward the upper‐left of the current layer where density is lower, similar to experimental findings (Yoo et al., 2017). The ion temperature shows variation in both L and N directions at larger spatial scales and may not be related to the reconnection event. The reconnection electric field shows spatial variation and is stronger upstream. Within the current sheet, when normalized to B v out, the reconnection rate is approximately 0.05, lower than the first event but within the range of values seen in Bessho et al. (2022).

Figure 4

Guide‐field reconnection site. (a) Current density J and magnetic field lines. (b) Current density J in the L‐N plane. L and N axes are red and green respectively, and the lines have length 10 d for scale. (c and d) magnetic flux transport (MFT) velocity showing the in‐ and outflow of magnetic flux, indicating active reconnection. (e) Electron L velocity in the frame moving with the x‐line. (f) Ion velocity L in the frame moving with the x‐line. (g) Electron temperature. (h) Ion temperature.

For these results, we note that a reduced mass ratio of 100 has been used due to computational limitations. This reduces the scale separation between electrons and ions, which is important for the realization of electron‐only reconnection (Sharma Pyakurel et al., 2019). With a reduced scale separation, the electron‐scale current sheets may not be much smaller than the ion scales, which could lead to slower reconnection rates, and may account for the partial coupling of ions in the second event. Current sheets in this simulation may also be less intense than with a realistic mass ratio—if we consider a domain with the same size in terms of ion inertial lengths, current sheets would be thinner due to the smaller electron scale, and if the typical magnitude of the downstream magnetic field fluctuations remains the same, this would lead to stronger field gradients across the current sheet.

Discussion

The three dimensional reconnection sites show similarities and differences compared to earlier two‐dimensional simulations of reconnection in quasi‐parallel shocks. Similar to two‐dimensional simulations, we find multiple reconnecting and non‐reconnecting current sheets, and the structure of reconnection regions is different from laminar reconnection regions such as those in the magnetotail, with asymmetric inflows and outflows (Bessho et al., 2019, 2020).

On the other hand, there were differences in the occurrence and regime of reconnection in the 3D simulation. In an equivalent two‐dimensional simulation (in the x‐y plane) we performed with the same physical parameters, we did not find reconnecting current sheets at the same stage of evolution. The additional degree of freedom allows reconnection to occur in different planes. This has more consistency with observations as they show that reconnecting current sheets are observed for a wide range of shock orientations (Gingell et al., 2020).

For the two events discussed in this paper, the orientation of the current sheets is such that J is mainly in the y direction, which would not be possible in a 2D simulation in which the reconnecting current is constrained to be in the z direction. In the x‐y plane, we see long wavelength waves in the upstream region propagating oblique to the magnetic field with wavelength ≈2.5d , similar to Bessho et al. (2020), where oblique waves with wavelength ≈3d are seen. In the 2D case, these fluctuations cause the bending of field lines which contributes to the occurrence of reconnection. In the x‐z plane in 3D simulations, the projection of the wavevector of the fluctuations is almost parallel to the in‐plane magnetic field, which has an angle of approximately 20° to the shock normal. The parallel fluctuations bring alternating B to the shock, which may account for the existence of reconnection regions oriented in the x‐z plane.

Another major difference between two‐ and three‐dimensional simulations is the ability to study weak guide field reconnection. In two dimensional simulations (e.g., Bessho et al., 2019; Q. Lu et al., 2021), the evolution of the system leads to the generation of strong out‐of‐plane magnetic fields. Reconnection events observed in such simulations are mainly in the strong guide field regime. In this three‐dimensional simulation, we find reconnection events with a range of guide fields, with the first event we focus on having a weak guide field since the orientation is such that the strong B contributes to the reconnecting component of the magnetic field. Again, this is relevant to observations, as a statistical study of reconnection in the shock transition region has indicated that reconnection sites have a wide range of guide fields, with stronger guide fields slightly favored (Gingell et al., 2020). The presence of weak guide field reconnection sites also has implications for the conversion of energy in 2D and 3D quasi‐parallel shock simulations. In reconnection regions, acceleration is most efficient when the guide field is comparable to the upstream field or weaker (Dahlin et al., 2017), suggesting that energy conversion could be more efficient in 3D simulations. However, it should be noted that Dahlin et al. (2017) studied standard reconnection, and the results may differ in the electron‐scale case. Nevertheless, it is important to understand whether acceleration mechanisms differ in 2D and 3D simulations, the total volume of the reconnection regions, and their impact on the overall energy conversion at quasi‐parallel shocks. These are open questions for future investigations.

The two events studied have reconnection rates on the order of 0.1 Bv Ae, in agreement with other simulation studies of reconnection in shock turbulence (Bessho et al., 2022). No particular enhancement of the reconnection rate for these events due to three‐dimensional effects is seen as in Pyakurel et al. (2021), but a further statistical study will need to be conducted to determine if the reconnection rate in shock‐driven reconnection is modified by 3D effects.