Particle-in-cell simulation study of the scaling of asymmetric magnetic reconnection with in-plane flow shear

Magnetic reconnection is the fundamental plasma process where a change in magnetic topology facilitates the conversion of magnetic energy to plasma kinetic energy and heat. It plays a fundamental role in causing eruptions in the coronae of the sun and other stars, in the interaction between the solar wind and the magnetospheres of Earth and other planets, for confinement in toroidal fusion devices, and in a large collection of astrophysical settings.¹ 

There has been increased interest in the properties of reconnection in realistic systems going beyond the simplifying assumptions of the Sweet-Parker model.^2,3 In this classical model, the magnetic fields, densities, and temperatures are equal on either side of the reconnection site, and the upstream plasmas has no bulk flow other than the inflow in the reference frame of the reconnection site. One example of going beyond this model is to allow asymmetries in the reconnecting magnetic fields, densities, and temperatures on either side of the reconnecting current sheet. A second example is including the effect of a bulk flow in the upstream plasma, whether in the direction of the reconnecting magnetic field or out of the reconnection plane.

Understanding how these effects impact the reconnection process, both quantitatively and qualitatively, is often of great importance for applying our understanding of reconnection to real systems. One example is reconnection at the dayside magnetopauses of Earth and other planets. The plasmas at the magnetopauses of Earth⁴ and Mercury⁵ differ on the two sides and the solar wind drives a bulk flow in Earth's magnetosheath⁶ and undoubtedly does at Mercury's, as well. When the interplanetary magnetic field (IMF) is northward, the magnetosheath flow is parallel/anti-parallel to the reconnecting magnetic field in the polar regions. When the IMF is southward, magnetosheath flow at the flanks is mostly out of the reconnection plane. The effect of upstream bulk flow is even more dramatic at the magnetospheres of Jupiter and Saturn, where rotation of the magnetosphere is much stronger of an effect than at Earth.⁷ 

A second example where upstream asymmetries and bulk flow are important is in tokamaks. The density and temperature profiles are peaked in the plasma core with a spatially varying magnetic field, so the plasma profiles at the reconnection site are non-uniform.⁸ Further, there are often bulk flows causing the toroidal and poloidal rotation of the plasma,⁹ especially those driven by neutral beam injection. Therefore, both asymmetries and flows effects are present and are important to the dynamics in magnetospheres and tokamaks.

While the effect of asymmetries and flow shear have separately received much attention,¹⁰ only a few studies have treated systems that simultaneously contain both effects. Studies of the shock structure far downstream of the reconnection site were carried out analytically^11–13 and using magnetohydrodynamic (MHD) modeling.¹⁴ Particle-in-cell (PIC) simulations were used to study systems simultaneously including asymmetries, flow shear, and an out-of-plane (guide) magnetic field.¹⁵ It was shown that the flow shear and diamagnetic drifts set up with the pressure gradient and the guide field can either reinforce or counteract each other.

More recently, a scaling analysis for systems including both asymmetries and upstream flow in the reconnection plane was performed.¹⁶ It was argued that the reconnection site (the X-line) typically convects in the outflow direction. The convection speed of the X-line and the rate of reconnection were predicted as a function of arbitrary upstream plasma parameters for isolated systems; the results will be reviewed in Sec. II. In symmetric reconnection with a flow shear, reconnection does not occur if the flow is super-Alfvénic because the tension in the reconnecting magnetic field cannot overcome the energy of the flow.¹⁷ There is also a critical flow speed above which reconnection does not occur for asymmetric reconnection; a generalization of the symmetric result for the asymmetric case was also derived.¹⁶ These predictions were successfully tested with two-dimensional numerical simulations using the two-fluid model (MHD with the Hall term and electron inertia). However, it is known that the fluid model is not well-suited to describe systems with asymmetric density and temperature as the fluids do not mix in the absence of thermal conduction;^18–20 even if conduction is present, the fluid model may not faithfully describe mixing in a nearly collisionless plasma as is the case in many applications. These shortcomings are not present in kinetic simulations, such as the particle-in-cell numerical technique²¹ where macro-particles are evolved in time and plasma mixing naturally occurs. Thus, it is important to investigate the scaling of the reconnection rate and the drift speed of isolated X-lines within a fully kinetic model.

In this study, we perform a systematic numerical study of magnetic reconnection with asymmetries and in-plane upstream flow using the PIC technique. We measure relevant quantities in independent simulations in which all quantities are held fixed other than the upstream flow. We find that the theoretical predictions previously tested in fluid simulations¹⁶ are consistent with the results of the PIC simulations.

In Sec. II, we review the predictions for the convection speed of isolated X-lines and the reconnection rate in terms of upstream parameters. Section III discusses the simulations we perform as well as our methodology for analyzing the simulation data. Section IV presents the simulation results and compares them to the predictions. We summarize our results and discuss applications of the results in Sec. V.

Scaling laws for the dissipation region's convection speed and the reconnection rate were derived¹⁶ for isolated configurations including asymmetric (but anti-parallel) magnetic fields and asymmetric densities and temperatures, along with arbitrary in-plane upstream flow. We define the magnetic field strengths of the two upstream regions as $BL,1$ and $BL,2$⁠, the plasma mass densities as ρ₁ and ρ₂, and the upstream flow speeds as $vL,1$ and $vL,2$⁠, where the subscript L is borrowed from boundary normal coordinates in magnetospheric applications to denote the reconnecting component, and the subscripts 1 and 2 denote the two upstream sides of the reconnection site. The reconnecting magnetic fields are treated as positive quantities, and the speeds are defined as positive if in the direction of $BL,1$ and negative in the direction of $BL,2$⁠.

The convection speed $vdrift$ of the X-line along the current sheet, the reconnection rate $Eshear,asym$⁠, and the upstream flow speed $vshear,crit$ required for steady state reconnection to be prevented scale (in cgs units) as

and

In writing these expressions, the asymmetric Alfvén speed $cA,asym$⁠, the asymmetric reconnection rate in the absence of a flow shear $Easym,0$⁠, and the velocity shear $vshear$ are

and

where δ and L_d are the half-thickness in the normal direction and half-length in the outflow direction of the dissipation region and c is the speed of light.

Equation (1) was derived using conservation of momentum in the L direction into and out of the dissipation region. Equation (2) follows from treating the energetics of the release of magnetic tension on the outflow jets while including the upstream flow, keeping track of the fact that the X-line and stagnation point are not collocated in asymmetric reconnection.²² Equation (3) is the condition that makes $Eshear,asym=0$ in Eq. (2), which is the condition for when reconnection shuts off. Equations (4) and (5) follow from the analysis of asymmetric reconnection with no upstream flow.²² Equation (6) is merely shorthand for the quantity of import for the reconnection rate in Eq. (2).

We emphasize a few important assumptions made in this analysis. It was assumed no upstream out-of-plane (guide) magnetic field. It was also assumed that the outflow speeds in the two downstream directions are equal and opposite in the reference frame of the moving X-line. As a scaling analysis, it assumes a single characteristic value represents each quantity in the appropriate regions.

We perform 2D kinetic particle-in-cell simulations using the P3D code²³ to test the predictions. Particles are stepped forward using a relativistic Boris algorithm, while electromagnetic fields are updated with a second order trapezoidal leapfrog. Magnetic field strengths are normalized to an arbitrary strength B₀, and plasma number densities are normalized to an arbitrary density n₀. Values of length and speed are normalized to the ion inertial length $di0=(mic2/4πn0e2)1/2$ and the Alfvén speed $cA0=B0/(4πmin0)1/2,$ respectively, where e is the ion charge and m_i is the ion mass. The unit of time is therefore $t0=di0/cA0=Ωci0−1$⁠.

As in the fluid simulations,¹⁶ the boundary conditions are doubly periodic and the magnetic field profile is initialized as a double Harris sheet

where $w0=1.0 di0$ is the initial current sheet width and L_y is the domain size in the inflow direction. (The x direction in the simulations corresponds to the L direction in boundary normal coordinates.) There is initially no out-of-plane guide field. The temperature profile of species j, which can denote electrons e or ions i, is initialized as

where $Tj1$ and $Tj2$ are selected asymptotic initial temperatures. We use $Ti1/Te1=Ti2/Te2=2$ for all simulations. Initial electron and ion densities are chosen to be equal with asymptotic values of n₀₁ and n₀₂. The density profiles initially enforce pressure balance across the current sheet in the fluid sense. There is no known general asymmetric kinetic equilibrium^24,25 (although there are approximations²⁶ that are not employed here). As in many previous studies, our system rings at early times but settles to a steady state by the time of interest for this study. By this time, the initial kinks in the initial magnetic field and bulk velocity profiles at the current sheets also smooth out and therefore do not present any problems.

The ion and electron bulk flow speeds are initialized as a double tanh profile

and there is no out-of-plane component of the flow. This is accomplished by loading particles with a Harris-type drifting Maxwellian distribution function with a non-zero $vj,x$ contribution given by Eq. (9), which is equivalent to a previously used approach.²⁷ The electron and ion bulk flow speeds in the x direction are identical.

For all simulations, the speed of light is $c=15 cA0$ and the electron mass is $me=mi/25$⁠. The time step for particles is $dt=0.006 Ωci0−1$ and the electromagnetic fields have a time step half as much. The grid scale is $dx=dy=0.05 di0$⁠. The simulations performed for this study are summarized in Table I. The set labeled “B” employs asymmetric magnetic fields with symmetric density, the set labeled “n” has symmetric magnetic fields with asymmetric density, and the set labeled “ms” is for representative magnetospheric conditions.²⁸ The domain size is $Lx×Ly=204.8×102.4 di0$ for the B simulations and $Lx×Ly=102.4×51.2 di0$ for the n and ms simulations. The initial number of particles-per-grid cell is 1000 for the B and n simulations and 500 for the ms simulations. The upstream flow speeds are varied for each set; Table II shows the values used in the present study.

To reach the nonlinear phase of reconnection more rapidly, the simulations are initialized using a coherent divergence-free sinusoidal perturbation to the magnetic fields of amplitude 0.1 with one full wavelength of the perturbation in the x direction and two full wavelengths in the y direction. Each simulation is evolved until magnetic reconnection reaches a steady-state. Since the reconnection rate differs for different upstream flow speeds, the steady state is reached at different times for different simulations. Consequently, we use the half-width of the primary magnetic island as a common indicator across the simulations; a range of island widths of 4–6 $di0$ is used to identify comparable times. If needed, the interval is slightly adjusted to ensure that the system is in a steady state.

For each time step and each current sheet, the X- and O-line are found using standard techniques by calculating the flux function ψ as $B=ẑ×∇ψ$⁠. The saddle point of ψ is the X-point, and the extremum is the O-point. The convection velocity of the reconnection site is measured as the time derivative of the X-line position. The reconnection rate is the time rate of change of the magnetic flux difference between the X- and O-line. These values are averaged over the steady state interval to provide a representative value for that simulation.

We begin by showing an overview of the plasma parameters in a simulation with representative magnetospheric conditions, specifically the ms simulation with $vL,1=1$ and $vL,2=0$⁠. Figure 1 contains (a) the out-of-plane current density J_z with magnetic field lines overplotted and (b) the out of plane magnetic field B_z at a time of 67.5 when the reconnection rate has reached a steady state, with the coordinate system shifted so that the X-line is at the origin. Only a fraction of the total computational domain is plotted. Interestingly, the results are quite similar to standard systems without flow shear. In particular, the Hall magnetic field in (b) is mostly bipolar, as is the norm in strongly asymmetric systems.^25,29 The normal electric field E_y (not shown) shows the typical asymmetric Hall electric field dominated by a positive E_y on the strong (magnetospheric) field side of the dissipation region and a negative Larmor electric field upstream of it.³⁰ Panel (c) shows the reconnection rate E as a function of time from this simulation, showing that the system reaches a steady-state by approximately t = 60.

Figure 2 shows cuts 4 $di0$ downstream of the X-line from the same simulation. The ion density n_i and ion bulk flow velocity $vi,x$ to the left and right of the X-line are shown in panels (a) and (b), respectively. They reveal the negative and positive deflections from the background flow profile due to the reconnection exhausts. Panel (c) shows E_y on the left axis and B_z on the right axis in a cut 4 $di0$ to the left of the X-line. The dashed horizontal lines mark zero for the two plots. The Hall magnetic (y = 0.5), Hall electric (⁠$y=−0.25$⁠), and Larmor electric (y = −1.5) fields are present. These results suggest that the upstream bulk flow largely does not alter the kinetic signatures of reconnection for typical magnetospheric conditions.

Another feature of reconnection with an upstream flow shear is the tilting of the current sheet near the X-line.^14,17,31 We see the current sheet tilt in the present simulations, as well. The tilting in the n simulations with asymmetric density is more pronounced than similar upstream flows from the B simulations with asymmetric magnetic field (not shown). The tilt is more pronounced for higher flow shear, as is to be expected.

Next, we test the predictions for the X-line drift speed, reconnection rate, and cutoff speed. Raw measured and predicted values for all simulations are given in Table II. We first consider the B simulations containing asymmetries in magnetic field strength. From Eq. (3), the predicted cutoff speed is about $2.2 cA0$⁠. We vary the upstream flow speed up to this cessation condition, with $vL,1=−vL,2$ for simplicity, and measure the drift speed in each simulation. The results are shown in the top plot in Fig. 3. Here and throughout, blue triangles and red squares are for the two current sheets in the initial setup. Error bars are determined as the standard deviation during the steady state intervals. A linear trend in drift speed for increasing shear flow speed is observed, qualitatively consistent with Eq. (1). The dashed black line is the prediction of Eq. (1), so the quantitative agreement is good as well.

We carry out the same analysis on the asymmetric density simulations, with results shown in the bottom plot of Fig. 3. The predicted cessation condition for these simulations is approximately $1.8 cA0$⁠. Again, the trend and absolute agreement is quite good. Note that some of the measurements from the simulation are slightly below the prediction for both the B and n simulations; this could be due to the inertia of moving larger primary islands in this finite sized and periodic domain.

Next, we test the reconnection rates obtained in both sets of simulations. The results are plotted in Fig. 4 for (top) asymmetric magnetic field and (bottom) asymmetric density sets, respectively. The dashed lines denote the predicted values; in calculating the predictions, we use the average reconnection rate from the two sheets in the zero flow shear case as $Easym,0$ in Eq. (2). We find excellent agreement between these measurements and the predictions.

A similar analysis is carried out for the representative magnetospheric simulations. From the prediction, the convection speed should double and the reconnection rate should drop by 5% as $vL,1$ is increased from 1.0 to 2.0. The measured convection speed is 0.88 for $vL,1=1.0$ and 1.64 for $vL,2=2.0$⁠, which differ by a factor of 1.86; this agrees well with the prediction. The reconnection rate drops from 0.143 to 0.123, a decrease of 13%. While ostensibly greater than the prediction, the uncertainty in these measurements is enough where the difference is not expected to be significant. The key is that the significant increase in upstream flow speed only slightly impacts the reconnection rate, as predicted for reconnection in systems with a strongly asymmetric density.¹⁶ 

As a test of the prediction for the cutoff flow shear speed to shut off asymmetric reconnection, we report results from simulations beyond the predicted cessation conditions. We run simulations with v_shear = 2.4 and 2.8 for the asymmetric field and v_shear = 2.0 and 2.4 for the asymmetric density simulations. The current sheets in these systems tend to contort and become sinusoidal rather than flat, as if beginning to go Kelvin-Helmholtz unstable. This is qualitatively different than the simulations below the threshold, where reconnection clearly starts from the beginning. We also point out that some of them nonlinearly experience reconnection in the strongly bent fields. We argue this form of reconnection is different than the robust form of reconnection discussed earlier, so we assert that the simulation results are consistent with Eq. (3).

Finally, we comment on whether the predictions remain valid when the upstream flow on the two sides is in the same direction, i.e., both $vL,1$ and $vL,2$ are positive. We do an asymmetric magnetic field simulation with $vL,1=vL,2=2$⁠. The predictions are that the X-line will drift with the common upstream drift speed and the reconnection rate is the same as if there was no flow. As shown in Table II, this is borne out in the simulations. In summary, the simulation results agree quite well with the predictions discussed in Sec. II.

We use particle-in-cell simulations to study the scaling of 2D asymmetric anti-parallel reconnection with an in-plane upstream flow. The particle-in-cell approach is necessary to faithfully capture the effect of plasma mixing between the two disparate plasmas in the exhaust region, which is not well-described in the fluid description. We find very good agreement with the scaling predictions as a function of upstream plasma parameters for the drift speed of the X-line, the reconnection rate, and the critical upstream flow speed necessary to suppress reconnection found in a recent study.¹⁶ 

One area in which this study goes beyond the previous fluid simulation study is by testing the theory for systems in which the flow is in the same direction on each upstream side, which has been studied for solar wind applications.^32–34 The results confirm the theory¹⁶ works in this case as well, and in particular, confirms that the figure of merit is $vshear∝vL,1−vL,2$⁠. Physically, for a case with $vL,1$ and $vL,2$ in the same direction, the upstream plasma on both sides enter the diffusion region with momentum in the direction of the reconnecting field, and this momentum makes the X-line convect at the weighted average of the two speeds. In the special case of equal flows on the two sides, the upstream plasmas are stationary in the reference frame of the X-line, so there is no effect on the reconnection rate. This reveals that reconnection in the solar wind should not be suppressed by flow shear, which is consistent with the observation of active reconnection in the solar wind.^35–38 

An important consideration before applying the results here to naturally occurring reconnection or reconnection in the laboratory is that the X-lines in question must be “isolated” in the sense that they are free to convect in the external flow for the theory to apply. This is essentially satisfied in the solar wind, for example, and may also be the case for neoclassical tearing modes (NTMs) in tokamaks. However, for the magnetopause of Earth and other planets, one should proceed with caution. A primary X-line is undoubtedly controlled by global considerations such as being line-tied to the ionosphere. Therefore, it is not clear if a single X-line would follow the predictions of the theory tested here. However, a flux rope or flux transfer event (FTE)³⁹ could be considered isolated, so the theory may apply. There are differences between the predictions tested here and the leading model of open flux (i.e., FTE) motion,⁴⁰ so future work on this topic would be interesting.

An area of potential interest for future study is the properties of the system in regimes where the flow is higher than the cessation condition. In our simulations, the current sheet contorts as if beginning to undergo Kelvin-Helmholtz instability. Then, reconnection as a secondary process on the bent current sheets begins. The current sheets flatten and reconnect robustly. It is not clear whether this is physical or only a function of the finite system size in the simulations, so it is worthy of future study. We point out that a limit with such strong flows is not likely to apply at the magnetosphere except possibly when the reconnection site interacts with the dense, cold plasmas in plasmaspheric drainage plumes,^41–43 but may be important in tokamaks.

One limitation of this study is that it is in 2D. Any 3D dynamics, including drift waves set up with wave vector in the out-of-plane direction due to the in-plane pressure asymmetry being normal to the reconnecting magnetic field, are artificially suppressed in 2D. It is not expected that drift waves will change the bulk properties of the reconnection, but this is worth further consideration.

Another limitation of the present study is that it does not take into account an out-of-plane component of bulk flow velocity or an out-of-plane (guide) magnetic field. It is known that diamagnetic effects arise in systems with a guide field and a pressure asymmetry,²⁴ and that the effects of flow shear and diamagnetic drifts compound.¹⁵ The effect of out-of-plane flow has been studied.^44–47 Incorporating diamagnetic effects into the theoretical predictions^48,49 should be the subject of future work.

Support from West Virginia University's Summer Undergraduate Research Experience (SURE) program (CED), the NASA West Virginia Space Grant Consortium (CED), NSF Grant AGS-0953463 (PAC), and NASA Grants NNX16AF75G (PAC) and NNX16AG76G (PAC and MS) are gratefully acknowledged. This research uses resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank K. Malakit for use of his color table in some of the plots, and we thank M. T. Beidler, C. M. Komar, and Y. H. Liu for stimulating discussions. The simulation data used to produce the results of this paper are available from the authors by request.