The reversibility of magnetic reconnection

The mean free path associated with classical collisions significantly exceeds the typical system size in many heliospheric and astrophysical plasmas. This simple observation raises interesting questions related to the dissipation of energy, particularly for phenomena such as magnetic reconnection,^1–3 turbulence,⁴ and shocks⁵ in which kinetic scales play an important role. In the complete absence of collisions, Boltzmann's equation, describing the evolution of the distribution function f, reduces to the Vlasov equation, which is formally reversible in time. This reversibility implies, among other things, that in a closed system the entropy—and in fact the integral of any function of f—remains constant during the system's evolution. Despite this restrictive mathematical constraint, Landau⁶ demonstrated that one-dimensional electrostatic plasma waves in the nominally collisionless limit are accompanied by a cascade to finer and finer structures in phase space which, in actual plasmas, leads to dissipation, irreversibility, and an increase in entropy. Landau damping has been experimentally confirmed,⁷ and subsequent work has shown that the same basic mechanism occurs in association with many other plasma oscillations.^8,9 While clever manipulation can resurrect temporary echoes of the original oscillation,¹⁰ in the long-time limit the existence of any non-zero amplitude phase-space scattering inevitably produces dissipation.

The example of Landau damping demonstrates that mathematically reversible processes, when coupled with any finite level of dissipation, can lead to irreversibility in real systems. An interesting question is whether that conclusion extends to the phenomenon of magnetic reconnection. During reconnection, a change in magnetic topology triggers a transfer of energy from the magnetic field as slowly inflowing plasma crosses magnetic separatrices and is accelerated into Alfvénic jets flowing away from an X-point.¹¹ In a strictly collisionless Vlasov-governed system, the transfer of energy from field to particles is completely reversible, implying that it should be possible for reconnection to run backwards, with narrow Alfvénic plasma jets converging at an X-point to produce broad sub-Alfvénic outflows and the entire system eventually evolving into a quiescent and structured current sheet. To our knowledge, such a sequence of events has not been reported in a real system, suggesting that, in practice if not in theory, collisionless magnetic reconnection is irreversible.

Data from the Magnetospheric Multiscale (MMS) Mission have demonstrated close agreement between in situ observations of reconnection and the results of particle-in-cell (PIC) simulations.^12–14 Such simulations follow macroparticles as they move through their collective self-consistent electromagnetic fields interpolated to a numerical grid and, in principle, capture all of the important dynamics in reconnecting plasmas. Since, as with the (assumed to be collisionless) physical system, the equations evolved in PIC codes are time-reversible, such simulations serve as useful proxies. The implementation of PIC codes typically includes features that can inject irreversibility such as numerical (round-off) errors and truncation errors arising from the discretization in space and time. If these errors were to be made infinitesimally small, time-reversibility would be preserved. However, we will show below that reconnection behaves irreversibly in typical PIC simulations, a result consistent with earlier work on time reversibility in simulations of molecular dynamics,¹⁵ and discuss key questions about what factors affect the reversal process. In Sec. II, we discuss test-particle trajectories in a simple analytical model of reconnection. Section III describes self-consistent particle-in-cell simulations of reconnection run forward and backward in time, while Sec. IV discusses the results and their broader implications.

Particle trajectories in self-consistent PIC simulations will be discussed below, but we first consider a simpler model that captures the key features of fields near a reconnection X-point. Observations and simulations have demonstrated that reconnection in a fully three-dimensional domain is often accompanied by turbulence and complicated field line trajectories.^13,16 Despite the complexity, however, many important features of reconnection, including the large-scale topology and the reconnection rate, remain similar to what is observed in simulations with reduced domains—sometimes called 2.5D—in which variations perpendicular to the reconnection plane are suppressed.¹⁷ Since, as we will show, reconnection is irreversible in such restricted domains it seems clear that it will also be irreversible in fully three-dimensional systems.

With that motivation, we seek a simple model of a reconnection X-point in a system with an invariant direction, here parallel to $ŷ$⁠. Two cases are of interest: Anti-parallel reconnection, in which the reconnecting fields, which lie parallel to $x̂$⁠, have a shear angle of $180°$⁠; and guide-field reconnection in which a spatially constant out-of-plane component of the magnetic field reduces the shear angle.

We consider a dimensionless model, i.e., one in which every variable is associated with a suppressed scaling factor that carries information about the proper units. Let the only non-zero component of the vector potential have the form

where α_B and α_E are positive constants and the electrostatic potential $ϕ=0$⁠. The magnetic and electric fields are then

Since E_y, which is equivalent to the reconnection rate, is constant in space, Faraday's Law and the invariance with respect to y imply B is stationary in time. No guide field is present so the un-reconnected fields on either side of the current layer are anti-parallel and $E∥=(E·B)/B$ trivially vanishes everywhere. Ellipses centered at the origin mark contours of constant magnetic field strength. Magnetic separatrices along the lines $z=±αBx$ divide space into upstream (⁠$|z|>αB |x|$⁠) and downstream (⁠$|z|<αB |x|$⁠) regions. Similar configurations have previously been investigated as models of reconnection.^18,19

Generally speaking, the two instances of α in Eq. (1) can take different values with $αB$ related to the angle of the separatrices and $αE$ to the reconnection rate. (If, for example, $αB=1$ and $αE=0$⁠, then the field lines meet at right angles at the X-point, the associated current density vanishes, and reconnection does not occur.) However,²⁰ it has been argued that constraints at magnetohydrodynamic scales impose the condition $αE≈αB$ when $αB≪1$⁠. In this work, we choose $αB=αE=α=0.1$⁠.

The symmetry inherent in the reconnection of anti-parallel fields is broken by an out-of-plane (guide) component, which can in practice have significant effects even when relatively small in magnitude.²¹ Including a constant $By≠0$ while maintaining the condition $E∥=0$ implies the existence of non-zero in-plane electric fields. An electrostatic potential $ϕ$⁠, with $E⊥=−∇⊥ϕ$⁠, generating these components satisfies

the solution to which can be found from the method of characteristics to be

with g an arbitrary function that can, in general, be chosen to satisfy boundary conditions but will be neglected here. The resulting in-plane components of E are

In-plane contours of the $E×B$ flow and magnetic field lines for the case B_y = 1 are shown in Fig. 1. Interestingly, despite the simplicity of the model, there are striking similarities to particle-in-cell simulations of guide-field reconnection,²² including strong flows across the $z=αx$ separatrix and weak sheared flows straddling $z=−αx$⁠.

The existence of singularities in $ϕ$⁠, E_x, and E_z along the $z=−αx$ separatrix has deeper roots than the specific choices of functional forms assumed in this model. First note that it is not possible to choose a g in Eq. (4) that eliminates the singularities. Since g is only a function of $z2−αx2$⁠, it has the same magnitude on both $z=αx$ and $z=−αx$ so that any cancelation along the latter will create a new singularity on the former. More generally, any 2.5D model with a time-stationary magnetic field that includes a spatially constant out-of-plane component and satisfies $E·B=0$ must have singularities associated with the vanishing of the in-plane field. In such a model, $Bx=Bz=0$ at the X-point, making it impossible to satisfy the other requirements there while simultaneously maintaining $E∥=0$⁠. When the model includes a guide field in this work, we only consider particle trajectories that do not cross the $z=−αx$ separatrix.

In order to explore the reversibility of individual trajectories, we introduce test particles into the model's magnetic and electric fields. The field configuration remains fixed as the positions and velocities of the particles are integrated forward in time using the non-relativistic Newton–Lorentz equations: $dx/dt=v$ and $dv/dt=q(E+v×B)$⁠. These equations are reversible, and so particles followed forward and backward in time should return to their initial locations. (As in Sec. II A, the variables are assumed to include a suppressed factor representing the dimensional information. For instance, $q=±1$ for protons or electrons, respectively, although in what follows we choose q = 1 for simplicity.)

In practice, however, the accumulation of small integration errors has the potential to alter trajectories. These alterations will be of particular significance when particles pass near locations such as X-points or separatrices where the $κ2$ parameter, introduced in Ref. 23 and given by the ratio of the radius of curvature of the magnetic field to the particle Larmor radius, most closely approaches values associated with chaotic motion.²⁴ Because X-points and separatrices do not fill the domain, the trajectory of every particle that passes near one effectively threads a needle such that any small change can lead to significant deviations when the orbit is tracked backwards in time. Here, the deviations arise from numerical errors in solving the governing equations, but in real plasmas, analogous behavior can arise from, for instance, small-amplitude electric field fluctuations. To emphasize the analogy with Landau damping: during reconnection in either an actual plasma or a numerical simulation, there will always be random fluctuations that sufficiently perturb particle trajectories enough to destroy phase-space structures at some scale. When such fluctuations coincide with regions associated with chaotic particle motions, irreversible particle trajectories result.

We follow trajectories via a trapezoidal-leapfrog method with the velocity integration employing a version of the Boris algorithm,^25–27 specifically the Boris-B solver described in Ref. 28. Exact numerical convergence is impossible (achieving it would effectively imply infinite phase-space resolution), since for any time step, $Δt$ there will always exist particles whose trajectories cannot be exactly reversed. In this work, we have chosen $Δt=10−3$ so that trajectories that do not pass close to the separatrices and X-point are reversible, but emphasize that tests with different values of $Δt$ yield qualitatively similar results. In general, reducing $Δt$ will increase the length of time and number of particles for which trajectories are reversible but will not eliminate irreversibility altogether. Here, we show trajectories evolved forward in time for $2.5×105$ steps and then reversed for an equal period.

Figure 2 shows two particle trajectories in the anti-parallel case of the model described in Sec. II A. At t = 0, both have representative thermal velocities, with the perpendicular velocity augmented with the local $E×B$ drift. The particles have the same guiding centers but different phases (offset by $π/2$⁠) so that they initially occupy different locations on the same Larmor orbit centered at x = 1 and z = 4.5. Both particles begin upstream and are tracked forward in time (blue lines) as they pass near the X-point and head downstream until, at t = 250, the particles have the locations given by the colored circles [top: $(x,z)≈(20,−3)$⁠; bottom: $(x,z)≈(27,6)$]. At the beginning of their trajectories, the particles exhibit the expected helical motions about the field with a pitch given by $v∥/v⊥$⁠. Conservation of the adiabatically invariant magnetic moment implies that $v⊥$ varies in proportion to $B$ so that the pitch decreases (i.e., the helix winds more tightly) at the mirror points and then increases as the particle moves back toward the origin. Magnetic moment conservation begins to break down near the current layer at z = 0 due to the decrease in B. Despite having initial conditions on the same Larmor orbit, the blue trajectories in the two panels clearly differ near the X-point. Once the particles move away from the current layer they re-magnetize, mirror in the increasing field, and again pass through the region of small field. The overall behavior is quite similar to that seen in particle trajectories taken from full PIC simulations of reconnection.^29,30

At t = 250, the particle positions are marked by circles, and their trajectories are then integrated backwards in time (red lines) by replacing $Δt$ with $−Δt$ in the discretization of the Newton–Lorentz equations. (Mathematically, identical results are found by keeping $Δt$ unchanged but mapping $v→−v$ and either $B→−B$ or both $q→−q$ and $E→−E$⁠.) While both particles at first closely re-trace their inbound trajectories—the red lines at first completely overlap the blue ones—each eventually acquires a significant deviation after which the incoming and outgoing trajectories are more or less independent.

Figure 3 shows μ, the first adiabatic invariant, for each particle as a function of time. As in Fig. 2, μ during the forward trajectory is shown in blue, while that for the reverse is plotted in red. The large spikes correspond to times when the particles approach the X-point and B approaches zero. The inset in each panel displays the times where the forward and backward values notably begin to deviate and the diamonds correspond to those plotted in Fig. 2. In the first case, the diamonds are shown slightly after the deviation becomes noticeable at t = 110 when the particles are oscillating across the current sheet at $(x,z)≈(7,−1)$⁠.

Figure 4 shows a similar plot for a configuration that includes a uniform guide field B_y = 1 and the associated in-plane electric fields given in Eq. (5). As in the anti-parallel case, the two particles begin with the same velocities and guiding centers (centered at x = 1 and z = 4.4) but gyrophases separated in phase by $π/2$⁠. The Larmor orbits again exhibit variations in pitch as the particles encounter regions where B changes. Notably, however, the presence of the guide field means that the particles' magnetic moment is mostly conserved, the $κ2$ parameter of Ref. 23 remains in the non-chaotic regime, and the particles never fully de-magnetize. As a consequence, passing through the central current sheet does not lead to significant differences in the trajectories as it did in the anti-parallel case. Furthermore, as expected, the reversed-in-time trajectories more closely track the forward-in-time motion. In the bottom panel, the reversal is nearly exact, with the red reversed trajectory essentially covering the forward blue trajectory.

The dependence on guide field can be considered in light of the $κ2$ parameter defined by Ref. 23. For this configuration, $κ2≡Rc/ρL$ where R_c and ρ_L represent the radius of curvature and the Larmor radius, respectively, can be evaluated analytically. The limit $κ2≫1$ corresponds to adiabatic motion, and as $κ2$ approaches unity, particle trajectories are expected to display signs of deterministic chaos. In the anti-parallel case, see Fig. 2, $κ2=0$ at the X-point, but for any finite B_y the minimum value occurs downstream of the X-point, on the line z = 0, with $κmin2∝By2$⁠. The Larmor radius component of $κ2$ brings in a dependence on $v⊥$⁠, but for the parameters shown in Fig. 4 $κmin2≈1$⁠. As discussed in Ref. 23, the onset of deterministic chaos, which is closely related to the reversibility of the trajectories, arises due to the overlap of resonances between the Larmor motion and the bounce motion across the current layer at z = 0. Such behavior is clear even in a qualitative comparison of the current sheet crossings shown in Figs. 2 and 4.

It should be emphasized that while the plotted trajectories are intended to be representative, changes in multiple factors—numerical method, length of integration, time step, initial position, or velocity—will produce different trajectories, although even numerical experiments with solvers employing much smaller error tolerances (not shown) exhibit deviations in the forward and reversed trajectories. A more comprehensive approach comes from examining self-consistent simulations of reconnection.

Particle-in-cell simulations follow the motions of many particles in their self-consistent electric and magnetic fields. A key question is whether the results found in Sec. II for individual trajectories continue to apply. We perform such simulations with the code p3d.²⁷ In its normalization, a reference magnetic field strength B₀ and density n₀ define the velocity unit $vA0=B0/4πmin0$ where m_i is the ion mass. Times are normalized to the inverse ion cyclotron frequency $Ωi0−1=mic/eB0$⁠, lengths to the ion inertial length $di0=c/ωpi0$ (where $ωpi0=4πn0e2/mi$ is the ion plasma frequency), electric fields to $vA0B0/c$⁠, and temperatures to $mivA02$⁠.

The system parameters follow those of the Geospace Environmental Modeling (GEM) reconnection challenge.³⁴ The computational domain measures $Lx×Lz=25.6×12.8di$⁠. The ion-to-electron mass ratio is taken to be 25, which is sufficient to separate the electron and ion scales (the electron inertial length $de0=0.2di0$⁠). The electron thermal speed $vth,e≈2$ is much less than the normalized speed of light, c = 20; the latter further implies that $ωpe/Ωce=4$⁠. The spatial grid has resolution $Δ=0.05$ in normalized units while the Debye length, $≈0.04$⁠, is the smallest physical scale. The time step is $Δt=0.001$⁠. Each grid cell in the asymptotic plasma contains $≈1500$ macroparticles. The boundaries are periodic in the horizontal direction and perfectly reflective conducting walls at the top and bottom of the domain. A small perturbation is made to the center of the current sheet at t = 0 to begin reconnection. As with many PIC codes, p3d's basic algorithm does not naturally ensure the satisfaction of Gauss's law and so usually employs divergence cleaning to match $∇·E$ and ρ. However, since this correction breaks the invariance with respect to time of the Boris algorithm, it is not used in this work. Nevertheless, the observed violation of Gauss's law is small, and companion runs that included divergence cleaning (not shown) suggest that the effect is insignificant.

Figure 5 displays the out-of-plane electron current density J_ey at four times during the run. Panels (a) and (b) display the system at $t=20Ωci−1$ and $t=24Ωci−1$ during the forward time integration. The expected features are present: slowly inflowing plasma, growing islands of reconnected flux, and Alfvénic outflow jets (not shown). At $t=24Ωci−1$⁠, panel (b), the run is stopped and the backwards integration begun. Depending on the details of the numerical implementation, small errors can be introduced when starting the backwards integration. However, tests with multiple schemes taking more or less care to ensure exact reversibility during the first few time steps yielded essentially indistinguishable results. In the following, we simply set $Δt→−Δt$⁠. Panel (c) shows the system after it has been integrated back to $t=20Ωci−1$ and should be compared to panel (a). While some reversal has occurred—in particular, the island widths are similar in panels (a) and (c)—the structure of the central current layer is clearly different. Panel (d) shows the system after it has been evolved further backwards in time to $t=12Ωci−1$⁠. If the system were fully reversible, the islands would continue to shrink as flux un-reconnects and the system would eventually contain a nearly uniform current layer. Instead, reconnection has begun again, as demonstrated by the increase in the island size. The overall morphology closely resembles the state shown in panel (b) when the time reversal began.

Figure 6 displays the out-of-plane current density at four times for a system that includes an initial uniform guide field, B_y = 1, but is otherwise identical to that shown in Fig. 5. As before, panels (a) and (b) show $t=20Ωci−1$ and $t=24Ωci−1$ from the forward time integration. In this case, a plasmoid forms in the slightly narrower current layer near the X-point and moves downstream, nearly fully reconnecting with the flux in the island by the time of panel (b). As is typical in guide-field reconnection, one of the separatrices, in this case the one stretching from lower-left to upper-right, is stronger than the other. After the time shown in panel (b), the backwards integration begins. In panel (c) [which should be compared to panel (a)], the plasmoid has emerged from the downstream island and propagated back toward the X-point. The two panels are nearly identical although small differences can be seen, e.g., at $x≈11$ and $z≈6.4$⁠. However, as the reversed integration continues, reconnection eventually begins again as can be seen in panel (d).

A curious feature of the new phase of reconnection is that the strong and weak separatrices switch [compare panels (b) and (d)]. The morphology of the separatrices before the reversal is a well-understood by-product of guide-field reconnection and arises from the interaction of streaming electrons and $E∥$ along the separatrices.³⁷ As noted above, mapping $t→−t$ in the equations governing the system's evolution is mathematically equivalent to keeping t unchanged while mapping $B→−B$ and $v→−v$⁠. In other words, reversing the flow of time in the system is synonymous with the instantaneous reversal of the velocity of every particle and the direction of the magnetic field. When the system again begins to reconnect [panel (d)], this change effectively changes the sign of $E∥$ and hence the structure of the separatrices.

Figure 7 shows the reconnected flux, as measured by the difference in A_y between the X-point and O-point, vs time for three simulations: The two already discussed and an additional one with B_y = 2. The slope of the curve gives the reconnection rate, which is $≈0.1$ for all three simulations at the time of the reversal. The left panel shows the forward integration, with the vertical offset at t = 0 due to the small initial perturbation. Each simulation is run until t = 24, the time shown in panel (b) of Figs. 5 and 6. The right panel displays the result of the backward integration. In each case, the system retraces its trajectory (un-reconnects) immediately after the reversal, but eventually stops and then begins to reconnect new flux. As was the case with the test particles discussed in Sec. II B, the systems with guide fields show better reversibility. A quantitative measure of reversibility is the ratio $r=(ψ24−ψmin)/(ψ24−ψ0)$ where ψ is the reconnected flux, ψ₀ and ψ₂₄ are the values at t = 0 and t = 24, and ψ_min is the minimal value reached during the backward integration. When r = 0, the system exhibits no reversal, while for r = 1, it returns to its initial state. For the runs of Fig. 7, $r(By=0)=0.27, r(By=1)=0.59$⁠, and $r(By=2)=0.68$⁠. As expected from the particle trajectories discussed in Sec. II, increasing the guide field increases the magnetization (i.e., $κ2$⁠) and the degree to which the system is able to reverse.

The number of macroparticles—invariably small compared to the number of particles in a real plasma—is a potential source of noise in PIC simulations that could plausibly make an anomalously large contribution to the reversibility of a system. However, Fig. 8 shows results from additional runs of the B_y = 1 case where the number of macroparticles per cell (ppc) varied from 25 to 12800 (compared to $≈1500$ in Fig. 7). In those cases, $r={0.32,0.41,0.55,0.59,0.59,0.59}$ for $ppc={25,100,400,1600,6400,12 800}$⁠, respectively. For the smallest values of ppc, PIC noise appears to play a significant role, but above $ppc≈400$ the effect is minimal. For all of the runs, the rates at which flux un-reconnects (the slopes of the curves in the right panel of Fig. 7) maintain the same nominal dimensionless value of $≈0.1$ observed in the forwards integration and typically observed in PIC simulations.

As a first approximation, if the reversed reconnection rate does not strongly vary from $≈0.1$⁠, the time at which minimum flux occurs in the right panels of Figs. 7 and 8 should be proportional to r—a larger r implies a longer time-to-minimum and vice versa. Interestingly however, the observed time-to-minimum exhibits less variation than might be expected for increases in both the guide field and ppc. [For example, in Fig. 7 the difference between $r(By=0)=0.27$ and $r(By=1)=0.59$ might suggest a time-to-minimum that differed by $≈2$⁠, in contrast to the observed $≈1.3$⁠.] However, corrections in the time-to-minimum arise due to the shape of the flux minima. The source of these corrections is not clear, although the generation of plasmoids, such as in the one seen in Fig. 6 for the B_y = 1 case may play a role.

Unlike the test particles discussed in Sec. II B, particles in PIC simulations move on trajectories that evolve based on self-consistently determined electromagnetic fields. In a perfectly reversible simulation, the trajectories from the forward and backward integration would completely overlap but, as might be expected from the irreversibility shown in Fig. 7, they do not here. Figure 9 shows the average total displacement between the forward and backward integration of fifty electrons [panel (a)] and ions [panel (b)] as a function of time for reconnection with B_y = 0 (black curves), B_y = 1 (red) and B_y = 2 (blue). The particles initially occupied a small region just upstream of the reconnection current sheet. Most passed the separatrices during the simulation and were located in the reconnection outflow when the backwards integration began (this point in time denotes t = 0 on the horizontal axis). Because the system is periodic in the horizontal direction and the particles cannot move too far upstream in the y direction due to canonical momentum conservation, the displacement between particles cannot grow without bound. Instead, in all cases the particles initially track their doppelgangers closely before exponentially diverging and then eventually reaching an asymptotic steady-state. As expected, the deviations occur substantially earlier when B_y = 0 due to passages through the demagnetized current sheet. As the guide field increases the correlation between the forwards and backwards particles is maintained for longer periods, although eventually divergence occurs for all three cases. While the larger mass of ions [panel (b)] means larger Larmor radii (and hence less magnetization than electrons), it also increases the gyroperiod; the net result is that the timescale for ion deviation is longer than that for electrons.

Previous work⁴⁰ has demonstrated fully reversible magnetic reconnection in a gyrokinetic simulation of a collisionless plasma. The gyrokinetic equations effectively consider the case of very strong guide field, i.e., a fully magnetized plasma, and hence can be considered a limiting case of the work discussed here. The observed reversibility is consistent with the trend discussed above that increasing the strength of the guide field increases the system's reversibility. In addition, the continuum nature of the gyrokinetic model caused the noise levels in the simulations to be quite small (near machine precision), although artificially adding noise led to poorer reversibility. Observations of reconnection have shown it to be accompanied by small-scale electric-field fluctuations that have some similarities to this type of noise.^41,42 It would be of interest for future work to compare the characteristics (amplitude, spectrum, etc.) of this noise to that in different simulation models in order to assess its possible impact on reversibility.

Irreversibility is closely tied to the question of dissipation of energy. Circumstantial evidence for dissipation during reconnection comes in the form of the observed increase in the downstream plasma temperature.^43,44 However, definitively proving the existence of irreversible dissipation is not straightforward. A frequently used measure, the identification of regions where $J·E>0$ (where E is measured in the frame of the electrons⁴⁵) can be misleading since reversible processes can generate such signals. On the other hand, such signals can be tied to the development of complex structures in phase space. As the complexity increases, weaker and weaker non-ideal processes are sufficient to trigger irreversible heating.⁴⁶ 

Determining whether a process is reversible is equivalent to asking whether the entropy remains constant. The generation of entropy in PIC simulations has been a subject of a recent study.⁴⁷ The kinetic entropy was calculated for collisionless 2.5D anti-parallel reconnection in a closed system. For a simulation with excellent conservation of the total energy, the authors found a monotonic increase in the total entropy of the system, albeit at a low level (⁠$≈3%$ for the entire run). Since the system lacked physical collisions, the entropy increase was attributed to numerical effects. Later work⁴⁸ modeled the entropy increase as an effective numerical collisionality, with the collision frequency dependent on such numerical factors as the resolution, time step, and number of particles. When considered in the context of this work, these results suggest that reversibility in magnetic reconnection behaves in a manner analogous to Landau damping. In principle, every particle in the system follows a trajectory given by the collisionless, reversible equations of motion while the system itself follows a narrow trajectory through a high-dimensional phase space. However, during reconnection the orbits of particles passing near magnetic separatrices and magnetic nulls are sensitive to small perturbations. Deviations from collisionless trajectories (e.g., arising from numerical errors in a simulation or a non-zero effective collisionality in a real system) push the plasma off its narrow path in phase space onto an irreversible state. The size of the perturbation needed to effect this change varies, although we have shown that a larger guide field is linked to higher reversibility (and hence smaller entropy increase and energy dissipations).

Several open questions remain. To what degree does reversibility depend on other plasma parameters (e.g., three-dimensionality, β, the magnetization parameter σ that parameterizes relativistic reconnection,⁴⁹ or asymmetry across the reconnecting current sheet)? Is it possible for PIC simulations to run in the truly collisionless regime for general reconnection configurations and hence exhibit true reversible behavior? If not, does this have implications for the widely accepted notion that PIC simulations provide a nearly complete picture of magnetic reconnection? Finally, is there a definitive measure of irreversibility and entropy generation that can be evaluated with data obtainable by either spacecraft or laboratory experiments?

The authors acknowledge support from NSF Grant No. PHY1805829 and NASA Grant No. 80NSSC19K0396. This research used resources of the National Energy Research Scientific Computing Center, 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.

Raw data were generated at the National Energy Research Scientific Computing Center large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.