Kinetic-scale flux rope reconnection in periodic and line-tied geometries Open Access

Flux ropes or plasmoids are helical bundles of plasma and magnetic field.¹ They are a natural three-dimensional (3D) extension of the two-dimensional (2D) concept of a magnetic island. In both nature and the laboratory, individual flux ropes can interact with each other and undergo magnetic reconnection,² altering the magnetic topology and converting magnetic energy into particle thermal energy. A well-known example is the massive release of energy associated with solar flares.³ A significant portion of coronal mass ejections are found to have a clear flux rope structure,⁴ and measurements indicate that magnetic reconnection is important in flux rope formation.⁵ Within 2D models, plasmoids form naturally within all large-scale reconnection layers due to the properties of the plasmoid instability.⁶ They are thought to play a key role in permitting fast reconnection within the MHD regimes⁷ and in facilitating the transition between the MHD and kinetic regimes.⁸ Thus, there have been considerable research efforts devoted to the topic of plasmoids, mostly within 2D simulation models.

The physics of kinetic-scale flux ropes is of great interest in magnetospheric dynamics, as is partly evidenced by the recent results from the Magnetospheric Multiscale Mission.⁹ Flux ropes are produced naturally during reconnection within the ion-scale current layers along the magnetopause^10–12 and within the magnetotail,¹³ and fully kinetic simulations have demonstrated formation and interaction of flux ropes on the ion-inertial-scale due to the initial onset of reconnection within these thin current sheets.¹⁴ Such flux ropes undergo rapid development and complex interactions that may influence the large-scale dynamics. The prevalence of kinetic-scale flux ropes in magnetospheric dynamics has motivated a number of laboratory experiments. Measurements on the Relaxation Scaling Experiment¹⁵ demonstrate flux rope coalescence¹⁶ and kinking of the individual ropes.¹⁷ Similar experiments at the Large Plasma Research Device¹⁸ (LAPD) have measured partially line-tied ropes rotating at a frequency¹⁹ that closely matches theoretical predictions,²⁰ and quasi-periodic reconnection events are observed between two kink-unstable flux ropes.^21,22

Flux rope dynamics in both fluid and kinetic regimes have been extensively studied in 2D systems using initial conditions in which magnetic islands are allowed to coalesce. For these problems, magnetic reconnection occurs across well-defined topological separatrices,^23–25 but the extension of these ideas to realistic 3D systems remains challenging for a variety of reasons. First, applications in space and astrophysics are not periodic, and thus the influence of boundary conditions is a major concern for flux rope dynamics. In 2D, plasmoid interactions are limited to coalescence and bouncing, while in 3D, flux ropes have considerably more freedom to interact and can also undergo additional instabilities such as kinking. Consequently, it is more difficult to define topology and compute 3D reconnection rates. Numerical simulations of kinetic-scale reconnection in 2D have identified the importance of the electron pressure tensor and electron inertia in supporting the reconnecting electric field (see, for example, Refs. 26–29), but the extension of these results to large 3D systems remains an active area of research in both simulations^30,31 and with spacecraft observations.³² 

In this paper, we study the above issues using fully kinetic particle-in-cell (PIC) simulations. We employ parameters similar to recent LAPD experiments, although computational practicalities preclude whole-device modeling at present, and in this initial paper, we consider the collisionless limit (while collisions are important in LAPD). Specifically, we focus on the interaction of two initially parallel flux ropes that are brought into contact through a combination of the kink and coalescence instabilities. We consider both periodic boundary conditions and line-tied boundary conditions with particle absorption and injection. The 3D reconnection rate is computed, and the results are contrasted with 2D simulations using the same initial equilibrium. Some significant differences between the 2D and 3D reconnection rate are observed, depending upon the guide field and the boundary conditions. In addition, we identify which terms support the reconnecting electric field in both the 2D and 3D cases and compare with measurements on LAPD. In the collisionless limit considered in this paper, we find that the non-gyrotropic part of the electron pressure tensor is primarily responsible for supporting reconnection along field lines passing through the quasi-separatrix layer^33–35 (QSL), corresponding to the approximate topological boundary between the ropes (see Sec. II A for definition). For field lines adjacent to the QSL, electron inertia contributes significantly to the non-ideal electric field.

The remainder of the paper is organized as follows. A brief overview of the theory of 3D magnetic reconnection is presented first in Sec. II A, and this is followed by a description of the model equations, the initial equilibrium, and the boundary conditions employed here. Results from 2D simulations are presented in Sec. III to allow for comparisons between 2D and 3D reconnection. The 3D results are presented in Sec. IV and are subdivided as follows. In Sec. IV A, we present the global evolution of three different sets of simulations: 3D periodic with no kinking, 3D periodic with kinking, and 3D line-tied with kinking. The reconnection rate in these simulations is computed and discussed in Sec. IV B, and we examine which terms in the electric field are responsible for supporting reconnection in Sec. IV C. Comparisons with experimental measurements on LAPD are presented in Sec. IV D. Conclusions are presented in Sec. V.

Significant theoretical efforts have been undertaken to describe the magnetic reconnection and plasmoid interaction in 3D systems, and the theory of general magnetic reconnection^36–38 provides one approach for identifying under what conditions reconnection occurs. Reconnection is defined as a breakdown in magnetic connectivity as a result of a localized non-ideal response, and this is quantified via the quasi-potential

where the path is taken along magnetic field lines. Within 3D reconnection theory, it has been suggested that the maximum of Ξ provides a measure of the 3D reconnection rate.³⁹ A detailed review of the theory of general magnetic reconnection may be found in Ref. 40 and references therein.

A related concept is that of a quasi-separatrix layer (QSL), a region where gradients in the magnetic field line mapping are much larger than average.^33–35 QSLs are a natural extension of the 2D concept of a separatrix,⁴¹ and QSLs identify regions with large field line compression in one direction and large expansion in the orthogonal direction, having the geometry of a hyperbolic flux tube.⁴² The QSL is parametrized by the squashing factor, Q, defined by tracing field lines from an initial position (x, y, z) to a final position (X, Y, Z)

where a guide field in y has been assumed here. The basic concept of the squashing factor is closely related to the field line exponentiation factor that has been used in other work.⁴³ QSLs are identified as likely locations where magnetic reconnection can occur in 3D,⁴⁴ but temporal changes in the field line topology also require large values of Ξ. In the flux rope formation simulations of Ref. 14, the location of the QSL and peaks in the quasi-potential, Ξ, are well correlated.⁴⁵ Investigations of flux rope dynamics on LAPD first applied the concept of QSL to experimental data,^19,46 and several numerical studies have used the concept of QSL in order to interpret simulation results.^47,48

We model the flux ropes using VPIC (vector particle-in-cell)⁴⁹, which solves the relativistic Vlasov-Maxwell system of equations for a finite number of macro-particles for species s. Here, we consider only ions and electrons, s = i, e. More details of the VPIC implementation may be found in Refs. 49 and 50. In the collisionless limit, the electron fluid momentum equation can be rewritten to solve for the electric field

where n_e, U_e, and $P ¯ e$ represent moments of the electron distribution function corresponding to the density, the fluid velocity, and the pressure tensor, respectively. Here, following Ref. 31, we have decomposed the electron pressure tensor into a gyrotropic and non-gyrotropic piece, $P ¯ e = P ¯ e g + P ¯ e n g$⁠, where $P ¯ e g ≡ P e ∥ I ¯ + ( P e ∥ − P e ⊥ ) b ̂ b ̂ , P e ∥ ≡ b ̂ · P ¯ e · b ̂$⁠, $P e ⊥ ≡ [ Tr ( P ¯ e ) − P e ∥ ] / 2$⁠, and the magnetic field unit vector is $b ̂ ≡ B / | B |$⁠.

The left-hand side of Eq. (3) represents the ideal electric field. If the terms on the right-hand side are negligible (or if the dominant terms can be written as the gradient of a scalar), then the magnetic field is said to be frozen into the electron fluid. The first two terms on the right-hand side are from bulk electron inertia, and the last two terms represent contributions from the gyrotropic and non-gyrotropic parts of the electron pressure tensor. We identify the terms on the right-hand side that are responsible for supporting the reconnecting electric field in our simulations.

Our simulations begin from the Fadeev equilibrium,⁵¹ an extension of the well-known Harris sheet that includes a periodic chain of magnetic islands. The mutual attraction between parallel currents in the islands results in a coalescence instability,^52–56 and magnetic reconnection occurs in the absence of any external forcing.⁵⁷ The magnetic field takes the form $B ≡ B g y ̂ + y ̂ × ∇ ψ$ with

Here, B₀ is the asymptotic magnetic field strength, B_g is the guide field strength, L_s is the half-thickness of the initial current sheet, and ϵ is a dimensionless parameter controlling the island width. The current density is $J = B 0 μ 0 L s Ψ y ̂$⁠, where

We consider the restricted domain $x ∈ [ − L x / 2 , L x / 2 ]$ and $z ∈ [ − L z / 2 , L z / 2 ]$⁠, where L_x = 4πL_s and L_z = 2πL_s. There are two neighboring islands at x = ±πL_s, z = 0, and the island separation is Δ_island = 2πL_s. Our 3D simulations consider y ∈ [–L_y/2, L_y/2] where L_y = 4L_x.

With the prescribed magnetic field, it follows that the pressure is $p = p bg + B 0 2 2 μ 0 Ψ$⁠, where p_bg is a constant. There is a natural decomposition into a uniform background population and a flux rope population with spatial structure, represented by the second term on the right, where each population may be represented by a drifting Maxwellian distribution. We further assume that each species is in thermal equilibrium; that is, the background and rope ions have the same temperature, T_i, and the background and rope electrons have the same temperature, T_e. The ratio T_i/T_e remains a free parameter, but the ratio of rope ion drift velocity to rope electron drift velocity is chosen to match this to provide an exact Vlasov equilibrium.⁵⁸ The density profile then follows as n = n_bg + n₀Ψ.

We use normalized PIC units, where mass is normalized to m_e, charge is normalized to the elementary charge e, time is normalized to $ω p e − 1$⁠, and velocity is normalized to c. The ratio of ion to electron mass, m_i/m_e, is a free parameter which sets the ratio of d_i/d_e. The in-plane magnetic field strength is determined from the ratio of electron plasma frequency to electron cyclotron frequency, B₀ ∼ ω_ce/ω_pe, and it is a free parameter in our simulations. The guide field strength is set in relation to the asymptotic field, B_g/B₀. The island width, ϵ, and current sheet half-thickness, L_s/d_i, are additional free parameters. The ratio of background density to rope density, n_bg/n₀, is another free parameter of the model, as is the ion to electron temperature ratio, T_i/T_e. The initial equilibrium is completely characterized by the following seven free parameters: ω_pe/ω_ce, B_g/B₀, ϵ, L_s/d_i, T_i/T_e, n_bg/n₀, and m_i/m_e. All of the simulations presented here use ω_pe/ω_ce = 10, ϵ  =  0.4, T_i/T_e = 0.6, n_bg/n₀ = 0.2, and m_i/m_e = 100 for computational practicalities.

In this work, we choose certain key parameters to overlap with recent flux rope experiments on LAPD. There, the flux ropes are separated by a distance of Δ_rope ∼10.5 cm, and the densities vary between 10¹² cm⁻³ in the background plasma and 4 × 10¹² cm⁻³ in the ropes themselves. The ion skin depth in the ropes for a hydrogen plasma is d_i ∼11.4 cm. Matching the island separation in the Fadeev equilibrium normalized to the hydrogen ion skin depth, Δ_island/d_i = 2πL_s/d_i, to the rope separation, Δ_rope/d_i ∼0.92, we find L_s/d_i = 0.15. (Although helium plasmas are predominantly formed in the experiment, if the helium skin depth is used, the island separation becomes even smaller, L_s/d_i = 0.075.) In units of ion skin depth, our computational domain is then x ∈ [–0.94d_i, 0.94d_i] and z ∈ [–0.47d_i, 0.47d_i].

The axial magnetic field strength is roughly 300 Gauss, so that the ratio of ion sound gyroradius to ion skin depth is ρ_s/d_i ∼0.15 at the rope center (β ∼0.026). Here, we vary only the guide field strength B_g/B₀. For the Fadeev equilibrium, a guide field strength B_g/B₀ = 1 corresponds to ρ_s/d_i ∼1.5 at the rope center, roughly a factor of 10 larger than in experiment. As B_g/B₀ increases to 3, ρ_s/d_i decreases to ∼0.5 at the rope center, and for B_g/B₀ = 5, we find ρ_s/d_i ∼0.3, a much better match to experimental parameters. In our 2D simulations, B_g/B₀ ∈ [1,10], while in our 3D simulations, we consider both B_g/B₀ = 3 and B_g/B₀ = 5.

Periodic boundary conditions are applied at x = ±L_x/2, while perfectly conducting boundary conditions are applied at z = ±L_z/2. For the 3D simulations, two types of boundary conditions are considered in the y-direction: periodic and line-tied. In the laboratory experiments, the magnetic field is line-tied at the cathode end, but the anode end is more complicated, and sheath resistance is believed to be partly responsible for breaking the line-tying condition there.⁴⁶ Since the exact boundary conditions at the anode are not well known, in this work, we employ simple line-tied (perfectly conducting) boundary conditions at both ends. Although the VPIC code has the capability of treating Coulomb collisions,⁵⁹ in this manuscript, we only consider the collisionless limit.

For 3D simulations with line-tied boundaries, particles are absorbed and periodically injected at the two axial boundaries. The particle injection scheme is described in greater detail in Ref. 60, though here the injection acts to drive the system toward the initial Fadeev equilibrium rather than acting as a free boundary. The background and rope populations are tracked and injected separately, and we have verified that the particle injection preserves the initial equilibrium with either a large guide field B_g/B₀ or a short axial extent L_y in order to inhibit kinking of the flux ropes.

Our 2D simulations utilize a uniform grid of 512 × 256 points in the x and z directions, respectively, and the grid spacing relative to the Debye length is Δx/Λ_de ∼ 0.66. In 3D, we use L_y = 4L_x with n_y = 1024, so that Δy/Λ_de ∼1.33 in the guide field direction.

In order to better place our 3D simulations into context, it is first useful to consider the 2D reconnection dynamics for the same initial conditions. Here, we vary the guide field strength, B_g/B₀, while all other parameters are as specified in Sec. II C. The 2D system is initially perturbed with

where $B ̃ = 0.05 B 0$⁠. Similar 2D simulations in the Fadeev equilibrium have been performed, albeit at higher values of L/d_i.^26–28,61 To filter out the effects of statistical noise, these simulations have been time-averaged over a window of roughly Δt/τ_A0 ≈ 0.055 (where τ_A0 is defined below).

As the guide field strength increases, the maximum reconnection rate decreases and the time of peak reconnection is also delayed, as can be seen in Fig. 1. The reconnection rate is determined from the magnetic vector potential for in-plane components of magnetic field, A_y, at the O and X points, and the reconnection rate is normalized as

where B₀ is the asymptotic field strength, $V A 0 = B 0 / m i n 0 μ 0$⁠, and τ_A0 = 4πL_s/V_A0. This normalization differs slightly from that used in previous work,^28,61 which used a normalization based on the maximum in-plane field strength, B_m, between the two island O-points. For our simulations, B_m = ϵB₀, so that the rate defined here is related to the previous rate by a factor ϵ².

Similar 2D simulations presented in Ref. 61 find that the reconnection rate increases with B_g/B₀ up to B_g/B_m = 2.83, which corresponds to our case of B_g/B₀ = 1, and decreases for larger values of B_g. We also observe decreasing reconnection rates with increasing B_g/B₀ above unity, but for B_g/B₀ < 1, we find that the reconnection rate is essentially independent of guide field strength (not shown), in contrast to the previous results. We attribute this to the fact that when B_g/B₀ < 1, the ion sound gyroradius is significantly larger than the system size in our simulations, while Ref. 61 considered much larger systems with L_s/d_i ≥ 1.

The contributions to the reconnecting (out-of-plane) electric field, broken down according to Eq. (3), for our 2D simulations are shown in Fig. 2. The dashed lines indicate electron density contours in order to highlight the merging islands, and the solid contour in each panel shows where $( 〈 E 〉 + 〈 U e × B 〉 ) · y ̂$ is half of its minimum value (the out-of-plane component of the reconnecting electric field is negative). At low guide field, B_g/B₀ = 1, the reconnecting electric field is almost entirely supported by contributions from the electron pressure tensor. At the X-point, contributions from the $∇ · ( n e U e U e )$ term in the inertia and the gyrotropic electron pressure tensor vanish, as expected analytically. The non-gyrotropic piece supports reconnection near the X-point, while the gyrotropic portion of the pressure tensor provides modest contributions away from the X-point.

For guide field B_g/B₀ = 3, the reconnection region becomes narrower and more elongated. The reconnecting electric field is still supported by the non-gyrotropic electron pressure at the X-point. However, away from the X-point, there are significant contributions from both the gyrotropic electron pressure and the inertia. In particular, we note that the electron inertia acts primarily to counteract the electron pressure contributions (primarily gyrotropic) in the outflow region. A similar behavior is observed as the guide field strength increases further to B_g/B₀ = 5. These results are consistent with previous 2D simulations in both the Harris sheet geometry²⁹ and island coalescence in the Fadeev equilibrium,⁶¹ which find a transition to electron inertia supporting reconnection away from the X-point with an increasing guide field.

The electron agyrotropy, A∅_e, is a scalar measure of the deviation of $P ¯ e$ from cylindrical symmetry about the local magnetic field (defined in Ref. 62), and it is also used to identify regions where the electron pressure tensor is non-gyrotropic. In the B_g/B₀ = 1 case, we find that A∅_e is large across the domain, reflecting the fact that the electrons are not well magnetized for this low guide field. For the B_g/B₀ = 3 and 5 cases, however, A∅_e is strongly peaked in the reconnection region with a maximum value at the time of peak reconnection of ∼0.20 and ∼0.17 for the two cases, respectively.

One interesting dynamical feature of these 2D guide field simulations is that the islands undergo a slight rotation about the central X-point, and the rotation is directly related to the guide field strength. In simulations with B_g = 0 (not shown), the islands simply coalesce with no rotation about this central X-point, and switching the direction of the guide field causes the rotation to occur in the opposite direction. This rotation may be understood as follows. Reconnection between the islands results in a quadrupole out-of-plane magnetic field structure, $B ̃ y$⁠, created by in-plane currents, $J ̃ x$ and $J ̃ z$⁠. These in-plane currents couple to the out-of-plane guide field, B_y = B_g, and exert a force which causes the islands to rotate as they coalesce. This same basic mechanism is also active in our 3D simulations, but the growth of the kink instability results in significantly larger displacements leading to even stronger rotation.

In this section, we consider 3D simulations with the same initial conditions as used in Sec. III. Our focus in this work is to consider similar parameter regimes as the LAPD experiments, which involve the interplay between kink-unstable flux ropes. Consequently, we consider B_g/B₀ = 3 or B_g/B₀ = 5 in order to allow the ropes to kink within the axial length L_y = 4L_x. Results from the latter case are shown in many of the following figures, though comparable results are obtained for the lower guide field. To filter out the effects of statistical noise, these 3D simulations have been time-averaged over a window of roughly Δt/τ_A0 ≈ 0.064, which is slightly longer than in the comparable 2D simulations. The global evolution of these 3D simulations is first presented in Sec. IV A. The 3D reconnection rate for these simulations is computed and comparisons with the 2D rate are made in Sec. IV B. The contributions of individual terms to the reconnecting electric field are analyzed in Sec. IV C, and comparisons with experimental measurements on LAPD are presented in Sec. IV D.

Here, we consider three different sets of 3D simulations. The first two are periodic in the axial direction. In the first case, the perturbation is given exactly by Eq. (6) with no axial dependence. These are quasi-2D simulations, as the coalescence instability occurs before any significant 3D kink features have time to develop, and they allow us to make direct comparisons between the 2D and 3D measurements of the reconnection rate. In the second set of periodic simulations, we include an axial $cos ( π y L y )$ dependence in the initial perturbation, and the ropes both kink and coalesce as they reconnect. These are referred to as kinking periodic simulations. The axial dependence used here excites many individual kink modes rather than just a single one, but it vanishes at the axial endpoints, y = ±L_y/2. This is to allow for direct comparisons with the third set of simulations, which are line-tied at both axial ends of the system.

The essentially 2D nature of reconnection in the quasi-2D simulations is evident in Fig. 3, which shows a volume rendering of the magnitude of current density along with representative magnetic field lines beginning at the periodic boundary at several different times for the B_g/B₀ = 5 case (similar results are observed for B_g/B₀ = 3). The two flux ropes drift together as a result of the coalescence instability, and there is no evidence of any significant kinking behavior. As in the comparable 2D simulations, the ropes rotate slightly about the central axis (X-point in 2D) as they reconnect.

When an axial dependence is included in the initial perturbation, the kink instability is excited in addition to the coalescence instability. The kinking is clearly evident in Fig. 4, which shows a volume rendering of $| J |$ along with representative field lines for both the B_g/B₀ = 3 (top) and B_g/B₀ = 5 (bottom) kinking periodic simulations. There is considerably more kinking evident in the B_g/B₀ = 3 case, as the larger guide field tends to stabilize the kink mode. Nevertheless, both cases exhibit kinking behavior, and the flux ropes first begin to reconnect near the axial midplane, y = 0.

The extent to which the kink instability impacts the 3D system may be quantified by examining the axial decomposition of the magnetic energy, $∫ B 2 / 2 μ 0 d V$⁠, shown in Fig. 5. For the quasi-2D cases, shown in the top subplots of Fig. 5, we see little significant growth of non-axisymmetric components in these simulations. Only a few low n values are shown, though the spectrum is resolved up to n = 512 (for n_y = 1024), and higher modes have significantly lower energy. In the quasi-2D simulations, there is a small decrease in the axisymmetric (n = 0) contribution to the total magnetic energy during reconnection, but it is not visible here, as the large guide field dominates over the reconnecting field. For B_g/B₀ = 3, the initial magnetic energy in the n = 0 component decreases by 0.7% of its initial value during reconnection, while for B_g/B₀ = 5, only 0.3% gets converted.

In contrast, when the kink instability is excited initially, we see significant magnetic energy associated with the non-axisymmetric Fourier components, as can be seen in the bottom subplots of Fig. 5. Although the initial magnetic energy in the n = 1 component is roughly 5 orders of magnitude lower than in the n = 0 component, it grows rapidly as the n = 1 kink mode develops. The n = 1 kink is the fastest growing non-axisymmetric component during the early linear phase of our simulations, and it always has an order of magnitude more magnetic energy than the next fastest n = 2 component. Around t/τ_A0 = 0.4 in both the B_g/B₀ = 3 and B_g/B₀ = 5 cases, we see evidence of nonlinear coupling into the n = 2 component, but saturation occurs before it overtakes the n = 1 component. The nonlinear coupling to higher n components and their subsequent nonlinear saturation occurs at slightly later times for both cases.

In the line-tied simulations, the flux ropes kink and coalesce as in the previous periodic case as can be seen in the volume rendering of $| J |$ shown in Fig. 6. As before, the B_g/B₀ = 3 case displays more violent kinking motions, but in both cases, the ropes initially interact near the axial midplane (y = 0). The representative field lines in Fig. 6 have fixed footpoints, so it is clearly evident that reconnection is occurring: field lines that initially begin at one footpoint reconnect in the central region and end at the opposing footpoint of the other flux rope.

In this line-tied system, the squashing factor, Q, defined in Eq. (2), is a useful measure of the topology of the magnetic field. It has been postulated that reconnection in 3D occurs within quasi-separatrix layers, defined as regions of large Q.⁴⁴ Here, we visualize the QSL as the region where Q is much larger than the average value in the domain, which varies depending on the guide field strength. The shaded purple surfaces in Fig. 6 correspond to regions where $log 10 Q ≥ 4$ for B_g/B₀ = 3 and $log 10 Q ≥ 3.5$ for B_g/B₀ = 5 (although the precise values do not matter much for visualization purposes). Note that Q is considerably larger for lower B_g, as the field lines are able to wander farther from their initial positions while traversing the domain axially. The hyperbolic nature of the QSL is clearly evident as well. Significant changes in the magnetic topology of the system occur as the field lines reconnect, and this is reflected in the evolution of the QSL.

The 3D reconnection rate (normalized by L_yV_A0B₀) and the time of peak reconnection in our periodic, quasi-2D simulations are remarkably similar to the strictly 2D simulations as can be seen by comparing Fig. 7 with Fig. 1. The yellow circles correspond to the quasi-2D case with B_g/B₀ = 3, and green squares correspond to the same case with B_g/B₀ = 5. The reconnection rate in our 3D simulations is determined from the maximum value of the quasi-potential, $Ξ ≡ − ∫ E · d ℓ$⁠, where the path is taken along the magnetic field lines.³⁹ Here, we utilize a uniform grid of 512 × 202 seed points spanning $[ − L x / 2 , L x / 2 ] × [ − 0.42 L z , 0.42 L z ]$⁠, and we integrate from y = –L_y/2 + Δy/2 to y = L_y/2 – Δy/2.

For B_g/B₀ = 3, the 2D reconnection rate peaks at t/τ_A0 ≈ 0.875 with a peak $E ̂ R ≈ 0.168$⁠, while in the quasi-2D periodic case, the peak rate of $E ̂ R ≈ 0.163$ occurs at t/τ_A0 ≈ 0.902. There is a slightly larger difference in the B_g/B₀ = 5 cases: in 2D, the peak $E ̂ R ≈ 0.146$ occurs at t/τ_A0 ≈ 0.930, while in the quasi-2D periodic case, the peak $E ̂ R ≈ 0.158$ occurs at t/τ_A0 ≈ 0.967.

We find slightly lower peak reconnection rates for the periodic kinking simulations than their quasi-2D counterparts. The red triangles correspond to the periodic kinking simulation with B_g/B₀ = 3 and the purple inverted triangles correspond to the case with B_g/B₀ = 5. For B_g/B₀ = 3, the peak $E ̂ R ≈ 0.122$ is roughly 75% of the comparable quasi-2D case, and for B_g/B₀ = 5, the peak $E ̂ R$ is about 90% of its quasi-2D counterpart. The lower peak rate is attributed to the fact that the kinking of the ropes prevents reconnection from happening along the entire axial length simultaneously, and there is more kinking at the lower guide field. However, the time of peak reconnection is generally similar, although this is expected to vary depending on the relative time scales of the kink and coalescence instabilities.

In the line-tied simulations, we find that the reconnection rate nearly matches that of the kinking periodic simulations for the first ∼0.8 τ_A0, but the peak reconnection rate is considerably lower, as the rope footpoints are fixed. For B_g/B₀ = 3 and B_g/B₀ = 5, the peak $E ̂ R$ in the line-tied simulations is roughly 93% and 75%, respectively, of the comparable kinking periodic simulations. We find a larger peak reconnection rate for the lower guide field simulation, in contrast to the kinking periodic results. Here, some kinking is required to drive the ropes together and reconnect, which is better facilitated with a lower guide field. Lastly, we note that the duration of reconnection is considerably shorter in the line-tied simulations. The peak reconnection rates are summarized in Table I.

We now examine which terms in the electric field are responsible for supporting reconnection. Each term in Eq. (3) is integrated separately along the magnetic field in our 3D simulations with the resulting value being a property of the field line itself. A map of these field line integrals to three discrete axial planes for the kinking periodic simulation with B_g/B₀ = 5 is shown in Fig. 8. The dashed lines here show the electron density contours and indicate the extent of the flux ropes on the different axial slices, highlighting the kinking motion. The ropes first begin to merge around the axial midpoint, y = 0. The first panel shows the ideal electric field integrated along field lines, $∫ ( 〈 E 〉 + 〈 U e × B 〉 ) · d ℓ$⁠, and this is essentially equivalent to the quasi-potential, $Ξ = − ∫ E · d ℓ$⁠. Contributions from $∫ 〈 U e × B 〉 · d ℓ$⁠, which may be non-zero along field lines only as a result of correlated fluctuations, are verified to be small. The reconnection region may be delimited by the region of large Ξ, and the hyperbolic shape is evident. The region of peak Ξ closely mirrors the electron pressure tensor contributions. The non-gyrotropic contributions are significant across a narrow region near the central axis where the gyrotropic contributions are small. Away from the central axis, we find significant competition between the electron inertia and the electron pressure tensor.

There is remarkable similarity between the field line integrated quantities in the kinking periodic simulation and the line-tied simulation at B_g/B₀ = 5 as can be seen by comparing Fig. 9 with Fig. 8. In both cases, the ropes first begin to merge near the axial midplane (y = 0), and the structure of the contributions to the reconnecting electric field is remarkably similar. The dominant contribution in the line-tied case is again the non-gyrotropic piece of the electron pressure tensor, and the gyrotropic contributions are small throughout the peak Ξ region. In the line-tied case, the QSL coincides with the region of large Ξ as can be seen in the solid contour of Fig. 9, which shows where log₁₀Q = 2.8 and approximates the QSL.

It is instructive to compare these 3D results at the axial midplane (middle subplots) with the strictly 2D results. We note that the electron inertia and the non-gyrotropic electron pressure tensor contributions both have a quadrupolar structure around the central axis (X-point in 2D). The regions directly above and below this central point correspond to the outflow region, where electron inertia inhibits reconnection driven by non-gyrotropic electron pressure contributions, as in our purely 2D simulations. In addition, the electron inertia appears to contribute significantly to the reconnecting electric field away from the central axis/X-point in the inflow region, where the pressure contributions drop off.

The peak reconnection rate in our line-tied simulations is comparable to the peak rate measured in LAPD.²² The reconnection rate in the experiment is defined in a similar manner as R_r ≡ Ξ/(LB_θV_A), where L is the axial length, B_θ is the in-plane field strength, and V_A is the Alfvén speed computed with the in-plane magnetic field. They find a peak reconnection rate of roughly R_r ∼ 0.1, on par with the results from the line-tied simulations.

Measurements on LAPD infer that the reconnecting electric field is supported by contributions from the divergence of the electron pressure tensor, in agreement with our simulation results. They find that $∫ E · d ℓ$ is between –0.5 and –1.0 V, while $∫ ( − 1 / n e e ) ( ∇ · P ¯ e ) · d ℓ$ is between –0.5 and –1.2 V.²² The electron inertia contributions are estimated to be small, and the remainder of the reconnecting electric field is believed to be supported by resistivity, although this remains to be measured in experiment (and is not included in our collisionless model).

Lastly, we note that the flux ropes in LAPD are observed to rotate,^19,22 but the causal mechanism invoked to explain that rotation is distinct from the small rotation previously discussed in Sec. III. In the experiment, the rotation is associated with the real frequency of the ideal external kink mode with line-tied boundary conditions at one end and free motion at the other, where the ideal MHD system is not self-adjoint.²⁰ The observed rotation frequency closely matches the theoretical predictions. We expect that the two-fluid rotation observed in our simulations is also at work in the kinetic-scale flux ropes in LAPD, although this effect is likely subdominant to the ideal MHD rotation.

In this manuscript, the physics of reconnecting kinetic-scale flux ropes is examined in both 2D simulations (classic island coalescence problem) and 3D simulations in which the flux ropes are kink unstable and boundary conditions play a key role. The 3D reconnection rate is computed from the maximum of the quasi-potential, Ξ, and the topology is characterized using the squashing factor. As expected, for 3D periodic simulations in which the kink instability is suppressed, the reconnection rate and dissipation physics are both in excellent agreement with 2D results. When kinking is excited in the periodic simulations via a suitable initial perturbation, we find that the peak reconnection rate is reduced by 25% and 10% for B_g/B₀ = 3 and 5, respectively, relative to the quasi-2D simulations. This reduction is larger for the lower guide field case, B_g/B₀ = 3, since the amplitude of the n = 1 kink mode is greater and less of the flux ropes are in contact at any one time. With a larger guide field, B_g/B₀ = 5, there is less kinking, and larger portions of the two ropes are able to reconnect simultaneously.

In the line-tied simulations, the peak reconnection rate is further reduced by 30% and 35% for B_g/B₀ = 3 and 5, respectively, relative to the quasi-2D simulations. Reconnection also ends earlier since less of the flux ropes are free to move and reconnect. Here, the peak reconnection rate is lower with a larger guide field, in agreement with the 2D results, but in contrast to the 3D periodic kinking results. In the periodic computations, kinking results in less of the ropes being able to reconnect simultaneously and lowers the peak reconnection rate, but with line-tied boundary conditions, some kinking is required to drive the ropes together and reconnect. The peak reconnection rates in our line-tied simulations are in good agreement with recent experimental measurements of reconnecting flux ropes.²² 

In all cases, the non-gyrotropic portion of the electron pressure tensor is predominantly responsible for supporting the reconnecting electric field along the magnetic field line with maximum Ξ. Contributions from the gyrotropic portion of the electron pressure tensor in 3D are verified to be small for field lines in the large Ξ region. Recent experimental measurements on LAPD also found that the electron pressure tensor can account for most of the reconnecting electric field,²² though no attempt was made to decompose this into gyrotropic and non-gyrotropic contributions. In addition, the effects of collisions were not measured in the experiment. However, we find that the electron inertia becomes comparable to the non-gyrotropic pressure tensor contributions away from the peak Ξ field line. This conclusion is consistent with a wide range of previous 2D simulations which identify a transition to electron inertia supporting the reconnection outside of the X-point for increasing guide field.⁶¹ In the line-tied simulations, the quasi-separatrix layer, defined as the region of large Q, coincides with the region of large Ξ, in agreement with laboratory experiments,²² theoretical expectations,⁴⁴ and previous simulation results.⁴⁵ 

For kinetic-scale flux ropes of relevance to both laboratory experiments and magnetospheric plasmas, these results demonstrate that key aspects of the 3D reconnection physics are remarkably similar to 2D models. This includes both the reconnection time scale (to within a factor of ∼2) and the dissipation physics. However, the extension of these results to larger astrophysical systems remains an outstanding problem, due to the potential influence of both turbulence and field line chaos. The field line exponentiation factor⁴³ in the present work is σ = ln(Q)/2 ∼5, owing to the small system size considered here. However, in many astrophysical systems, it may be as large as σ ∼20, and these large values have been theorized to lead to very different behavior.⁶³ The scaling of 3D magnetic reconnection up to much larger system sizes remains an interesting open question.

The authors acknowledge helpful discussions with Dr. Stanier, Dr. Le, and Dr. Ohia, and are grateful for the suggestions from the reviewer of this manuscript. This work was funded by the Basic Plasma Science Program from the Department of Energy Office of Fusion Energy Sciences. 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. DE-AC52-06NA25396.