Kinetic simulations verifying reconnection rates measured in the laboratory, spanning the ion-coupled to near electron-only regimes

Magnetic reconnection¹ is the process through which the topology of magnetic field lines changes in the presence of a plasma, often resulting in an explosive release of magnetic energy. Well-known examples include solar flares² and auroral substorms in the Earth's magnetosphere.³ Reconnection is studied both in situ in the magnetosphere by satellites like NASA's MMS mission⁴ and in laboratory experiments here on Earth. One such experiment is the Terrestrial Reconnection EXperiment (TREX) at the University of Wisconsin-Madison,⁵ which is operated specifically to reach parameter regimes of magnetic reconnection relevant to those of collisionless space plasmas.⁶ Previous results from TREX have characterized the rate at which reconnection occurs in the experiment.⁷ The reconnection rate determines the speed at which plasma and magnetic field lines can move into and out of the reconnection region, effectively setting the timescale of the entire process.^8–10 These TREX results showed the importance of magnetic flux pileup and the formation of a shock preceding the reconnection layer in maintaining force balance and setting the normalized reconnection rate. Furthermore, the experimental rate has a dependence on the size of the system relative to the ion scale; smaller scale size produces higher rates, indicative of the transition from ion-coupled toward electron-only reconnection.^7,11,12

Similar to previous numerical simulations of the TREX geometry,¹³ in this Letter, we will apply fully kinetic simulations with the aim to confirm and reproduce the results of Ref. 7. This is part of an ongoing effort to synchronize data collection between the experimental and simulated TREX environments, using the VPIC code developed at Los Alamos National Laboratory.^13–16 After a brief introduction to the TREX experiment and the VPIC code, multiple simulations of TREX will be analyzed to verify that pressure balance exists across the reconnection shock front. Finally, further simulations of TREX will be evaluated to check if the reconnection rate results match the values measured in TREX and their dependence on normalized experimental system size.

An engineering schematic of TREX is presented in Fig. 1(a). The vacuum vessel, provided by the Wisconsin Plasma Physics Laboratory (WiPPL),⁵ is a 3 m diameter sphere that uses an array of permanent magnets embedded in the chamber wall to limit the plasma loss area to a very small fraction of the total surface area while keeping the bulk of the plasma unmagnetized. The setup includes a set of internal drive coils and an exterior Helmholtz coil that provides a near-uniform axial magnetic field with a magnitude up to $28 mT$⁠.^6,7 The current through the three internal drive coils (purple) ramps up to create a magnetic field that opposes and reconnects with the background Helmholtz field, resulting in an anti-parallel magnetic configuration (e.g., no significant guide field). The plasma source is a set of plasma guns located at the machine's pole (shown in yellow). This setup mimics the asymmetric conditions of the dayside magnetopause; the high-density, low-field inflow at low R (analogous to the solar wind) is opposed by the low-density, high-field inflow at high R (analogous to the Earth's magnetic field). TREX is typically operated in either hydrogen, deuterium, or helium plasmas.

In the planar cut of TREX shown in Fig. 1(b), the cyan lines illustrate the typical magnetic geometry of an experimental run. As the current through the drive coils ramps up, the reconnection region is pushed from underneath the drive coils radially inward [orange arrows in Fig. 1(b)]. With a typical reconnection layer speed of $50 km/s$⁠, the temporal resolution of our probes (⁠$10 MHz$⁠) translates to a high spatial resolution measurement of about $5 mm$⁠. This process, referred to as the “jogging method,” permits the magnetic structure of the entire reconnection geometry to be characterized in a single experimental shot. The various magnetic and temperature probes and their locations are represented by the blue and red rectangles in Fig. 1(b). In addition to the jogging method probes, a different array of three-axis $Ḃ$ probes can be moved between shots, allowing for the creation of multi-shot datasets. The coverage area of this probe is given by the light green rectangle in Fig. 1(b).

An example of data collected from a typical set of experimental shots is provided in Fig. 1(c), where data from 34 shots are combined into one picture; for each shot, the probe is at a different position within the green region in Fig. 1(b). The black lines are contours of the flux function $Ψ$ to illustrate the in-plane magnetic field lines. Typical plasma parameters near the reconnection region include $Ti≪Te ≃5–20 eV, ne ≃2×1018 m−3, and Brec ≃4 mT$⁠, yielding $βe ≃0.4$ and $S ≃104$⁠.

TREX is simulated using VPIC, a kinetic particle-in-cell code.^14–16 VPIC has previously been used to mimic the TREX setup and produce results comparable to experimental data. More information on the general usage of VPIC to simulate TREX may be found in Ref. 13. The number of grid-points in the 2D simulations described here is 756 × 1800, spanning a system size of about 4.5 × 9 ion skin depths (d_i) in the R and Z directions in our standard density case. The low R boundary is set at $R≈0.155di$ in the standard case and acts as a reflecting conductor; this is meant to replicate the effect of the cylindrical TREX current layer bouncing off of itself once it reaches $R=0 m$ in the experiment (cylindrical VPIC cannot operate at R = 0, so a relatively close value is chosen for the lower bound instead; the closer this value is to 0, the higher the computational load). The average number of macro-particles per cell is 500. In all simulations presented here, the ratio of the electron cyclotron frequency to the electron plasma frequency is 1, and the mass ratio $mi/me$ is 400. Sub-realistic mass ratios are typical of PIC simulations for computational tractability; as a consequence, the experimental size measured in d_e is slightly larger than the simulation domain. However, Ref. 7 ran the TREX experiment at different ion masses and verified that the reconnection rate is tied to d_i rather than d_e. A typical experimental reconnection drive is modeled by injecting a linearly increasing current through the simulated drive coils. The strength of this current drive is given as a fraction of a typical experimental measurement of the time rate of change of the current in the drive coils, $İ0≈6.28×108 A s−1$⁠.

In TREX and TREX simulations during a reconnection discharge, the upstream plasma below the current layer can be divided into several distinct regions, as illustrated in Fig. 2(a).⁷ Working in the reference frame of the reconnection layer (r, red), the far upstream (u, blue) moves toward the layer at a speed faster than the local Alfvén speed. This necessitates the formation of a region of magnetic flux pileup (p, purple); this region is separated from the far upstream by a sub-critical shock (s, green). Pressure balance between these regions was verified using TREX data in the shock's reference frame in Ref. 7; the assumptions and approximations involved with this calculation will be detailed below.

When combined with Ampère's law, the MHD momentum balance equation for the plasma in regions u and p is

where ρ is the plasma density, v is the plasma flow speed, B is the magnetic field, p is the plasma kinetic pressure, and $d/dt$ is the total convective derivative, $d/dt=∂/∂t+v·∇$⁠. By plugging in the full convective derivative and using a vector identity, Eq. (1) can be rewritten as

Due to the toroidal symmetry of our experiment and the periodic boundary conditions in $ϕ$ in our simulation, we assume that $∂a/∂ϕ=0$ for any quantity a. To evaluate the pressure balance between regions u and p, we will integrate Eq. (2) over a path $dr$ between arbitrary points g and h, where these points share the same value of the Z coordinate such that $dr=drr̂$⁠. The resulting momentum/pressure balance relation is

By taking g and h to be in regions u and p, respectively, we can evaluate the change in the different terms of this equation across the shock that separates the two regions; this analysis is shown in Fig. 2(b). Here, the equation has been split into distinct terms, where the magnetic pressure is in blue, the magnetic curvature is in cyan, the total convective acceleration is in red, the total kinetic pressure (given as the sum of the ion and electron pressures) is in green, with the total (the sum of all the terms) given in black. All terms are evaluated in the reference frame of the shock. The value for each term can be averaged over the points highlighted in regions u and p to give a single value for each, resulting in Fig. 2(c); in this plot, the values in s are simply the difference between the u and p regions. This averaging was done to mimic the limitations of the TREX experiment: the speeds of different regions cannot be measured simultaneously, so only a single value for each term can be calculated in each of regions u and p. Both u and p are taken to be outside the electron diffusion region, such that the electron contribution to the inertia term is neglected.⁷ As expected, the total momentum/pressure is constant across the shock layer.

When evaluating the pressure balance across the shock layer in the experiment, several approximations are needed to account for the limitations of data collection (including the region speed limitation detailed above). Most notably, the analysis in Ref. 7 assumed that in regions u and p, changes in the plasma's velocity with respect to time or spatial coordinate are both minor relative to the ram pressure term, $minivr,i2$⁠, and balanced by the change in the magnetic tension term. The experimental analysis also assumed that $Te≫Ti$⁠. Both of these assumptions are tested in Fig. 2(d), where the simulation data are used to recreate TREX's measurements. The total of the approximate pressure terms, shown in purple, is constant across the shock, as it is in the full momentum balance analysis detailed above and shown in Figs. 2(b) and 2(c). From this, we conclude that the assumptions that went into Ref. 7's are also consistent with the presented numerical results.

As described previously, reconnection in TREX is asymmetric in plasma density and magnetic field on the opposing sides of the reconnection layer. As such, the reconnection rate is appropriately normalized by the method derived in Ref. 17,

where α is the normalized reconnection rate, E_rec is the reconnection electric field, B_red is the reduced magnetic field, $Va,hyb$ is the hybrid Alfvén speed, and B and n are magnetic field and plasma density values at locations 2 and 3 as shown in Fig. 2(a), respectively. Experimental values of α were calculated in TREX, where the reconnection electric field was derived from the time rate of change of the magnetic flux function, $Ψ$⁠.⁷ Results from this evaluation are shown as the yellow points in Fig. 3(a); the normalized reconnection rate, α, is plotted against the normalized system size, $L/di$⁠, where d_i is the ion skin depth, $di=c/ωpi$⁠. Similar to other analyses of simulated reconnection,^18,19 the rate increases as the normalized system size decreases. This is consistent with the transition toward electron-only reconnection, where the ions become less coupled to the field lines on the scale of the reconnection region, allowing reconnection to proceed without being constrained by the inertia of the ion fluid.^11,12 The full analysis in Ref. 7 also showed that the effective reconnection rate is constant regardless of the applied current drive $(İ)$⁠; this is due to the interplay between particle density and magnetic field strength upstream of the layer. Even if the layer is forced down with a stronger drive, producing a larger E_rec value, the shock structure and flux pileup will develop in such a way to produce a similar increase in the product of B_red and $Va,hyb$⁠, resulting in a constant value for the scaled rate α.

The TREX experiment is typically operated between two different density settings with three different ion species (hydrogen, deuterium, and helium), yielding six experimental points shown in Fig. 3(a). To compare these results to simulated TREX setups in VPIC, we instead vary the value of initial plasma density. Examples of applied initial density profiles are shown in Fig. 3(b). Our standard density profile vs R is shown in red and labeled as n. To reach a range of values for the scaled system size, this standard density was varied up and down by a single factor multiplying the entire profile; for example, a profile of twice the standard density (2n) is shown in yellow, while another profile of half the standard density (⁠$0.5n$⁠) is shown in blue. Note that the profiles shown here have had their density gradient decreased from the actual profiles used in the simulations for the sake of clarity. The real density profiles are based on measurements in the experiment and show a much stronger dependence on the value of the R coordinate.

Simulations were run for density values as low as $n/5$ and as high as 4n, where n represents our standard density. Three runs of each density value with different random initialization seeds were compared to reduce the chances of an anomalous result skewing the conclusions. Within each individual run, multiple points in time corresponding to the reconnection region being in the range of R values most readily measured by the experiment (⁠$R=0.35–0.45 m$⁠) were selected, and data from each of these points in time were sampled in regions p and r. These data were then used to calculate the reconnection rate [following Eqs. (4)–(6)] and the local ion skin depth. The results from all these time points from each of the three repeated simulations of a given initial density setup were averaged together to produce the green points in Fig. 3(a); the error bars represent the uncertainty estimate obtained by propagating the standard deviation of the distribution of selected density and magnetic fields through the rate equations. In general, these points follow the same trend as the experimental points (yellow), with increasing rate as the system size decreases.

One point of interest in Fig. 3(a) is the dip in the green simulation rate results localized around $L/di≈2.2$⁠. This feature is a real aspect of the data trend, tied to subtle Alfvénic wave dynamics related to our cylindrical reconnection drive scenario. When the drive current begins to ramp up, the pressure balance of the initial configuration is suddenly violated, causing an Alfvénic perturbation to propagate from beneath the drive coils downward toward $R=0 m$⁠. Although this wavefront is propagating down, the wave itself corresponds to a radially outward expansion of the plasma. In the standard scenario with density scale n [Figs. 3(d) and 3(h)], this expansion persists throughout the reconnection layer formation and inward propagation, allowing the reconnection dynamics to adjust in a manner that keeps the effective rate consistent with expectation, as described earlier and in Ref. 7. However, in the lower density scenario [$n/3$⁠, Figs. 3(e) and 3(i)], the initial wavefront travels inward and then reflects off of the low R boundary while the reconnection layer is still evolving. On the tailing side of the reflected wave-front, the plasma expansion is significantly reduced, corresponding to a transient reduction in the drive as the front reaches the reconnection layer. The upstream conditions of the reconnection layer cannot instantly adjust to these effects, resulting in normalized reconnection drives that can be either enhanced or reduced.

So far, this feature has not been clearly observed in the experiment, as it exists in a parameter regime that is not reachable in TREX. TREX has reached values of $L/di≈2.2$⁠, but this was done with helium and deuterium plasmas rather than by going to lower density values. Furthermore, the effect of this Alfvénic feature may be influenced by our simulation's reflecting boundary condition at low R, which approximates TREX's behavior but may not be exactly analogous. This dip would not be expected in scenarios without some manner of reflection along one of the domain boundaries. Similarly, the results in Fig. 3(c) show variation relative to drive current ramp intensity due to the above feature's effect in extinguishing the upstream inflow. This feature only appears in the simulations, causing them to diverge from the results of Ref. 7 that showed experimentally that the drive intensity does not affect the scaled rate.

Experiments conducted in the Terrestrial Reconnection EXperiment (TREX) over a range of different scaled system sizes showed a range of reconnection rates, which increased as the system size decreased. As part of an ongoing effort to model the TREX experimental setup in a particle-in-cell simulation, VPIC was used to replicate TREX runs at a range of densities, many of which were outside TREX's normal operating parameters. Within the range of the TREX parameters, the numerical simulations confirmed the experimentally observed rates of reconnection, weakly dependent on the normalized size of the experiment with higher rates at smaller system size indicative of the transition toward electron-only reconnection. Additionally, the high detail of simulation data allows the full pressure balance equation across the reconnection shock front to be calculated and pressure balance to be confirmed. This calculation also allowed us to verify the accuracy of some of the assumptions that were needed in TREX's experimental pressure balance calculation. Together with previous results,¹³ these conclusions continue to verify VPIC's ability to accurately capture the full shock formation and reconnection dynamics observed in TREX.

We gratefully acknowledge DOE funds Nos. DE-SC0019153, DE-SC0013032, and DE-SC0010463 and NASA fund No. 80NSSC18K1231 for the support of the TREX experiment. In addition, the experimental work is supported through the WiPPL User Facility under DOE fund No. DE-SC0018266. Simulation work was supported by the DOE Basic Plasma Science program and by a fellowship from the Center for Space and Earth Science (CSES) at LANL. CSES is funded by LANL's Laboratory Directed Research and Development (LDRD) program under Project No. 20180475DR. This work used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001.

The authors have no conflicts to disclose.

Samuel Greess: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Software (supporting); Writing – original draft (lead); Writing – review & editing (equal). Jan Egedal: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Adam Stanier: Data curation (supporting); Methodology (equal); Software (lead); Supervision (equal); Writing – review & editing (equal). Joseph Olson: Data curation (equal); Formal analysis (supporting); Methodology (equal); Project administration (supporting); Writing – original draft (supporting); Writing – review & editing (equal). William Daughton: Funding acquisition (lead); Methodology (supporting); Project administration (lead); Resources (lead); Software (equal); Supervision (equal); Writing – review & editing (equal). Ari Le: Methodology (supporting); Software (equal); Writing – review & editing (equal). Alexander Millet-Ayala: Conceptualization (supporting); Data curation (equal); Methodology (supporting). Cameron Kuchta: Data curation (supporting); Formal analysis (supporting); Software (supporting); Writing – review & editing (supporting). Cary B. Forest: Conceptualization (supporting); Funding acquisition (lead); Project administration (lead); Resources (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.