Guide field effects on the distribution of plasmoids in multiple scale reconnection

Impulsively fast magnetic reconnection is a ubiquitous phenomenon widely observed in magnetized space, solar, astrophysical, and laboratory fusion plasmas.^1–4 However, fast collisionless reconnection mechanisms based on non-MHD (Magnetohydrodynamic) effects such as two-fluid or kinetic effects (e.g., Ref. 5) are only applicable to scales much smaller than system sizes of these magnetized plasmas in space and astrophysics. In large systems, collisional MHD models such as the Sweet–Parker model are commonly used to describe reconnection processes. However, the predicted Sweet–Parker reconnection rate is much too slow to be consistent with many large scale phenomena like solar flares. At high Lundquist numbers, it has been realized that these current sheets can potentially fracture in a cascading process which results in a significantly increased reconnection rate.⁶ The instability responsible for this process is the plasmoid instability, which can grow rapidly leading to a reconnection rate that remains fast and is weakly dependent on the Lundquist number.^7–10 Thus, the plasmoid instability has been proposed as a promising mechanism to couple global system scale to local dissipation scales in reconnecting current sheets, although there exist alternative models based on MHD turbulence (e.g., Ref. 11).

In the MHD regime, this instability is found at sufficiently high Lundquist numbers,⁷ causing the breakup of long current sheets into chains of self-contained magnetic islands. These magnetic islands are highly dynamic, interacting with each other as they move through the length of the primary current sheet until they exit. Electrons can be efficiently accelerated to high energies by these dynamic magnetic islands^12,13 emitting observable radiation (e.g., Ref. 14) The acceleration process is based on the reflection of particles from the ends of each contracting magnetic island resulting in first order Fermi acceleration. Additionally, higher energy particles that are free to move throughout the primary current sheet across multiple islands can be mirrored and accelerated by an enhanced Fermi process.^15,16 A better understanding of the magnetic configuration of plasmoids and their statistical properties may help reveal how plasmoids can contract, the relative velocities with which they interact, and hence how the particles within and in between are accelerated.

Previous studies have shown how this instability can interrupt a Sweet–Parker current sheet which surpasses $Sc∼104$ depending on pre-existing noise.^9,17 A plasmoid unstable current sheet undergoes a cascading process which generates new plasmoids until their mean separation is low enough that the local Lundquist number of reconnecting current sheets between them has been reduced below S_c, thus stable to the plasmoid instability.⁶ These plasmoids, however, are highly dynamic through their interaction. Between the time when a plasmoid is born and the time when it advects out of the reconnection layer, it may absorb smaller plasmoids or be absorbed by larger plasmoids. Reconnection of their surrounding current sheets leads to the accretion of flux as well. Over time, the current sheet may produce very large plasmoids on the order of one tenth the length of the sheet or greater.¹⁸ 

Due to the importance of plasmoids in reconnection adjacent phenomena, statistical scalings have been sought both analytically and numerically to describe the distribution of plasmoids most commonly in terms of their reconnecting flux ψ, as $f(ψ)$⁠. In the Hall-MHD (magnetohydrodynamic) regime, $f∝ exp (−ψ)$ behavior has been predicted.¹⁹ This exponential dependence of reconnecting flux has been observed in the near-Earth space with limited in situ measurements by a few satellites,^20,21 by remote-sensing measurements of solar eruptions,²² and in laboratory experiments^23,24 but with limited resolutions. In the 3D MHD regime without a guide field, an entropy variational principle has been used to derive²⁵ a power law index of 3, while 2D relativistic particle-in-cell simulations have been paired with Monte Carlo methods to uncover a power law index in the range of 1–2. The often reported scaling is $f(ψ)∝ψ−2$⁠, seen in many 2D MHD simulations and justified in most cases by statistical approaches.^18,26–29 It has also been argued that the combination of a power law index of 1 with an exponential tail could appear as a power law index of 2 in numerical results, and analytical models explaining this behavior have been proposed.^22,27 Specifically, Huang et al.²⁷ developed a statistical approach that deftly demonstrates how the power law dependence of the distribution can be replicated by accounting for the essential behaviors of the unstable current sheet. Two models, one which allows for variation in the relative velocity to the mean outflow of these plasmoids and one which does not, produce power law indices of 1 and 2 of the reconnecting flux, respectively. All of these models do not explicitly take into account of the presence of a finite guide field in the reconnecting current direction. In many natural circumstances, however, a finite guide field is present during reconnection which may modify the distribution of plasmoids. Work by Ni et al.³⁰ found that the presence of a guide field does not have a major impact on the instability itself. We expand here the approach of Huang et al. to investigate effects of a finite guide field on the plasmoid distributions.

Both models of Huang et al.²⁷ are further developed here by adopting force-free field profiles internal to plasmoids as a result of Taylor relaxation³¹ to determine the distribution of plasmoids in the presence of a finite guide field. Without a guide field, the current sheet and the plasmoids within are essentially non-force-free, e.g., the incoming reconnecting field pressure is balanced by plasma pressure. With a finite guide field, however, plasmoids become magnetic flux ropes³² which are long known³³ to relax toward a force-free state possessing a strong core field in an approximately cylindrical shape.³⁴ Such force-free fields are also known in the laboratory pinch experiments as a result of Taylor relaxation,^31,35 during which magnetic energy is minimized while conserving magnetic helicity.³⁶ The force-free fields, $B→$⁠, are given by^35,37

where λ is the eigenvalue determined by boundary conditions. The lowest order solution to this equation was found to be $Bz=B0J0(λr)$ and $Bθ=B0J1(λr)$ in cylindrical coordinates. Here J₀ and J₁ are the zeroth and first order Bessel functions of the first kind, respectively. Taylor showed³⁵ that in toroidal pinch experiments with large Lundquist numbers and moderate toroidal fields, these profiles reasonably predicted the peaked toroidal (core) field at the center with a much-reduced magnitude or even reversed direction at the edge by extending beyond the first zero of J₀.

Magnetic reconnection with a guide field can be thought of as transporting magnetic helicity from the background into the plasmoids in the current sheet. If these plasmoids with a finite magnetic helicity are allowed to relax toward a minimum energy state while advecting toward the current sheet exit, their internal field structures will take a form of force-free profiles. When plasmoids merge, magnetic helicity can be dissipated in the secondary current sheets. However for high Lundquist numbers (⁠$≥Sc$⁠), this dissipation is low for moderate to weak guide fields³⁸ while magnetic energy is reduced significantly during merging. Therefore, the plasmoid merging process can be also regarded as a Taylor relaxation process to minimize magnetic energy while conserving the total magnetic helicity so that the resultant plasmoids also take the form of force-free profiles. Making use of these assumptions, the antiparallel reconnection model is modified to include the effects of Taylor relaxation with a finite but moderate guide field, following plasmoid mergers where adequate time for relaxation is available. These details are discussed in Secs. II B and IV. Hence, distribution of reconnecting and guide field plasmoid fluxes can be obtained. An alternative model is also provided which simply adds guide field fluxes upon plasmoid coalescence without Taylor relaxation.³⁹ This represents the strong guide field regime where Taylor relaxation is not allowed and the plasmoids are not force-free while non-negligible magnetic helicity is dissipated in the merging secondary current sheets. For both the relative velocity independent and velocity-dependent models, the non-force-free and force-free distributions are compared and their implications and differences are discussed.

In Sec. II A, we derive the relationships necessary to relate the reconnecting and guide field fluxes of a plasmoid to its magnetic helicity, providing the rate at which plasmoids gain magnetic helicity from the background reconnection. Some assumptions of the cylindrical plasmoid model are also discussed, along with their validity regimes. Section II B lays out further assumptions about magnetic helicity dissipation as relevant to the construction of a statistical equation that conserves magnetic helicity in plasmoid mergers. The non-force-free, guide field conserving model is also provided as the alternate regime of a strong guide field. Section III provides the numerical methods and solutions to the equations derived in Sec. II B. These distributions and their features are discussed in detail in Sec. IV. With solutions in hand, the assumptions made earlier can be verified or constrained. Section V summarizes the results, especially where relevant to particle acceleration. Possible future improvements to the statistical equations are also discussed.

Several simple but important relations between force-free plasmoids and the reconnecting current sheet where they are embedded are described in this section. Each isolated plasmoid either when they are born or as a result of merging of two plasmoids is assumed to take the force-free profiles, as predicted by Taylor through minimizing magnetic energy while conserving magnetic helicity. Unlike Taylor, however, both the guide field and reconnecting field fluxes are not conserved; rather they are determined by the embedded current sheet as follows. We assume that these plasmoids are simple axisymmetric cylinders as shown in Fig. 1(a) embedded in the reconnecting current sheet. The relations between these plasmoids and the background current sheet are derived in a straightforward manner by relating the plasmoid radius a to the eigenvalue λ. At the edge of the plasmoid r = a, we match the background reconnecting and guide fields (termed B_rec and B_g, respectively) to the azimuthal and axial magnetic components of the cylindrical plasmoids, $Bθ$ and B_z, respectively. We note that this is less valid for smaller, newly born plasmoids as their radius is not sufficiently large to reach the asymptotic values of B_rec and B_g. The impact of this assumption, however, is likely small and will be discussed in Sec. III. Eliminating B₀ from the profiles, we have that

There are infinitely many solutions to this equation. However, we assume that the solution with the lowest value of λa should be used. This is because solutions with higher λa involve a reversing guide field, yielding higher energy configurations. These configurations are unstable to kink modes past the critical $λa=3.176$⁠, leading to nonaxisymmetric dynamics which we do not seek to capture in this model.⁴⁰ This assumption is supported by observations that do not typically exhibit a reversed guide or core field within plasmoids. For the rest of this paper, the lowest solution with λa satisfying Eq. (2) will be used. Figure 1(b) shows two such examples of lowest λa: one for $Bg=0.75Brec$ and one for B_g = 0. For each choice of $Bg/Brec$⁠, the corresponding λa is determined as shown in Fig. 1(c). Furthermore, once λa is determined, B₀ may also be determined from both $Bg=B0J0(λa)$ and $Brec=B0J1(λa)$⁠.

The reconnecting flux per unit length (ψ) and the guide field flux (⁠$ϕ$⁠) contained within a plasmoid will be defined as follows:

For a toroidal system, ψ is the analog of the poloidal flux per unit axial length and $ϕ$ is the analog of the toroidal flux. In the construction of our statistical kinetic equation, it will become necessary to relate the reconnecting flux, guide field flux, and magnetic helicity through the background quantities by eliminating individual dependencies on λ and a. The eigenvalue equation can be integrated over an area with the $ẑ$ normal, which using Stoke's theorem and Eq. (3) gives

Performing the same integration with the $θ̂$ normal, we obtain

The loop integration was performed along the magnetic axis for a distance of l and returns on the cylindrical surface where the axial magnetic field component equals the background B_g. There is no contribution to the loop integration at the ends of the plasmoid due to the radial magnetic field being zero in the model. Dividing Eq. (4) by the square of Eq. (5) with rearrangements leads to

where the constant dimensionless term has been renamed as ζ for simplicity. The same conclusion can be arrived at from dimensional analysis when the goal is an expression relating ψ and $ϕ$ which does not expressly contain statistical variables such as the plasmoid radius a. Instead, only the combination of λa is involved as specified in Fig. 1(c) for a given $Bg/Brec$⁠. Using $∇×A→=B→$ the vector potential can be found for each plasmoid, namely,

where the gauge choice has been made such that $Aθ(λa)=0$ for each plasmoid to be isolated to avoid generating magnetic helicity through linking neighboring plasmoids and surroundings. (Otherwise, the concept of relative helicity^41,42 needs to be invoked.) The magnetic helicity for each plasmoid may then be calculated directly using $K=∫A→·B→d3r$⁠. The resulting relation can be manipulated to yield

where the constant coefficient has been renamed α. This gives a simple relationship for the magnetic helicity which depends only on ψ and the background parameters,

The growth rate of the plasmoid magnetic helicity due to reconnection if the plasmoid relaxes more rapidly than the rate by which the primary current sheet reconnects (an assumption to be justified in Sec. II B), is given by

The assignment $dψ/dt=γ$ has been made for brevity. To determine the final state of a plasmoid that has formed as the result of the merger of two plasmoids, we may then exploit either guide field flux conservation or magnetic helicity conservation depending on the background conditions.

A standard governing kinetic equation for plasmoid size distribution in terms of the reconnecting field flux includes contributions from the sources, growth, and sinks of plasmoids as they evolve in the dynamic current sheets. The source of plasmoids is always localized at zero flux while their growth is based on the overall reconnection rate. The latter is determined by the critical Lundquist number for a Sweet–Parker current sheet to become unstable:⁸ $Sc∼104$⁠. Therefore, the stable current sheets which lie in between plasmoids add flux to the plasmoids at a rate $dψ/dt=BVA/Sc=γ$⁠. Following Huang and Bhattacharjee,²⁷ with the assumption that mergers occur with a frequency independent of the plasmoid relative velocities, the statistical kinetic equation is given by

where

All of the solutions obtained are steady state with $∂t=0$⁠. On the right, the first term, $ξδ(ψ)$⁠, represents the creation of plasmoids where $δ(ψ)$ is Dirac δ-function. The second term represents the loss due to absorption by larger plasmoids with a rate proportional to their number, $n>(ψ)$⁠, assuming a typical relative speed on the order of V_A. The third is the sink term representing the advection of plasmoids out of the reconnection layer. The Alfvén time $τA=L/VA$ is that of the reconnecting component of the background current sheet. The solution for this equation was determined to be²⁷ 

where $C=1+2/N$ and N is the number of plasmoids. This solution possesses three regimes: a constant initial section where plasmoids grow due to reconnection not having had time to merge yet, followed by a power law dependence with the index of 2 where mergers dominate dynamics, then finally an exponential tail where advection loss dominates for the largest plasmoids. When accounting for plasmoids with varying velocities relative to the mean outflow, the statistical kinetic equation takes on a modified merger rate according to

where

The plasmoid creation term is now endowed with a velocity distribution, assumed to be $h(v)=(1/πVA) exp (−v2/VA2)$ (although as noted by Huang et al.²⁷ and confirmed by our own tests the solution is not very sensitive to the exact distribution chosen). The coalescence term corrects the rate at which mergers occur by the difference in velocity between the two merging plasmoids, approximating the resultant velocity as that of the larger plasmoid. An analytical solution to the steady state form of this equation is not known, it is instead solved numerically.

Equations (11) and (14) can be modified to account for the effects of Taylor relaxation with a guide field. As the reconnecting flux of a plasmoid grows by reconnection, its magnetic helicity also grows. This change in magnetic helicity calls for a continuous modification to the plasmoid internal profiles if there is sufficient time for the plasmoid to continuously relax. This condition is generally satisfied as the plasmoids grow on the tearing mode growth time while relaxation is mostly Alfvénic, evidenced for example from experimental observations.⁴³ Furthermore, the relaxation time is defined with respect to a given plasmoid's local Alfvén time. For current sheets with many small plasmoids, the vast bulk of the plasmoid population will have the local Alfvén times orders of magnitude faster than the global Alfvén time, thereby consistent with globally super-Alfvénic growth rates.^10,44,45 The relaxation time is also much shorter than the plasmoid advection time in the current sheet or equivalently the reconnection time on order of $ScτA=100τA$⁠. Therefore, it is well justified to assume that all plasmoids in the current sheet, other than those having very recently undergone coalescence, are always in a force-free state at any given time minimizing their magnetic energy while conserving magnetic helicity. For these plasmoids then, the reconnecting field flux and helicity are related according to Eq. (9), $α(ζ/B0)ψ3=K$⁠.

As hinted above, in many instances it will be assumed that there is ample time in between mergers for a given plasmoid to reach the force-free state. Small plasmoids rarely encounter plasmoids with a lower flux, and large plasmoids experience a small change in magnetic helicity when merging, thereby not deviating significantly from the force-free state. For intermediate plasmoids, however, there may exist a regime where plasmoids merge too frequently at comparable magnetic helicities to fully relax to a force-free state in between mergers. Therefore, we develop models in two opposing limits: one model for all plasmoids to be force-free after relaxation while conserving magnetic helicity even during plasmoid mergers, as described in Subsection II C. The other model on non-force-free plasmoids without relaxation while conserving guide field flux, as described in Subsection II D. The range of validity for the force-free assumption will require further discussion; however, it will be reserved for Sec. IV.

Regardless of whether the resultant plasmoid is force-free, magnetic helicity is conserved if the helicity change or dissipation in the secondary current sheet during coalescence is negligible. The change in magnetic helicity during plasmoid merging is due to the change in the linkage between the guide field flux contained in the secondary current sheet and the reconnecting field flux around the current sheet. An estimate of the rate of change in magnetic helicity per change in magnetic energy, W, is given by Ji,³⁸ 

where δ_s and L_s are thickness and length of the secondary current sheet, respectively. Since $δs/Ls$ is the steady state reconnection rate for the plasmoid merging, it should range from 0.1 to $1/Sc=0.01$⁠, depending on collisionality. Therefore, when the guide field is comparable to or weaker than the reconnecting field, no significant helicity is expected to be generated or dissipated by secondary current sheets. This permits our consideration of the plasmoid merger process as simultaneous Taylor relaxation. In a strong guide field, the dissipation of magnetic helicity would not be negligible and so the merging process would not be one of constant helicity. In fact, Taylor relaxation does not typically occur in this limit in the laboratory fusion plasmas and resultant profiles are generally not force-free.

With the assumptions above, the statistics of the new models can be described with modification of Eqs. (11) and (14). The nonvelocity-independent variable is chosen to be K, and the corresponding reconnection growth is given by Eq. (10). The loss term now includes advection and plasmoid mergers while a source is added to account for the resultant plasmoids which reached the K in question by merging. If we denote the sources and sinks as $ΣS(K)$ and note that the influx of probability density is given by $f(K)dK/dt$⁠, then the conservation of particle number for a dynamic distribution yields

Taking the steady state and equating the integrands, the distribution of plasmoids with a velocity independent merger frequency is therefore given by

where N is the total number of plasmoids. The growth rate is present inside the derivative in order to ensure that, if only reconnection growth were present, f is conserved. Here all mergers trigger a change in magnetic helicity, so the quantity $n>$ of Eq. (11) is instead integrated from 0 to $∞$⁠, becoming N. The merger term assumes the frequency of plasmoid mergers between $K−K′$ and $K′$ is $f(K′)dK′/τA$⁠, adding $f(K−K′)$ to f(K), in a similar manner to Fermo et al.¹⁹ where plasmoid area is conserved. It is integrated to $K/2$ in order to avoid double counting. Combined with the merger loss term, this kinetic model is constructed such that magnetic helicity is conserved by plasmoid mergers. Each merger transforms the resultant plasmoid to a new magnetic helicity which, ignoring reconnection and advection terms, keeps the total helicity in the distribution constant. Similar modifications can be made to the velocity-dependent model,

with

Note that the notation of F for velocity-dependent distributions and f for velocity-independent distributions will be used from here on. The same process is used as before to modify the merger rate to conserve magnetic helicity. Once again, H_t is simply H of Eq. (14) integrated from 0, and the resultant velocity of a merger is assumed to be that of the larger plasmoid. A proof of the helicity conserving property of the merger terms, Eq. (20) and the last line of Eq. (19), is presented in  Appendix by showing their first helicity moment is zero. It is easily extended to the velocity-independent case as well. Both of these force-free models are solved numerically in Sec. III.

As an opposing limit of the force-free plasmoid model, which is more applicable in the strong guide field regime, we offer an alternative model where plasmoid internal field profiles are not necessarily force-free but total guide field flux is conserved by mergers. In this case, the kinetic equation in ψ is unchanged from Huang et al. but the kinetic equations in $ϕ$ are similar to those of the relaxing model in terms of K,

and

for constant and nonconstant velocities, respectively. Here H_t is the same as Eq. (20) with the substitution $ψ→ϕ$⁠. Similar logic is employed for the source and loss terms, but the growth rate due to reconnection is modified. For a current sheet of thickness δ, while ψ is proportional to $Brecδ$ per unit length, $ϕ$ is proportional to $Bgδ2$⁠. The resultant growth rate is then $γg=Bgδ2/τA,local=Bg(l2/Sc)(VA/l)=BgVAL/NSc$⁠, where the local current sheet length $l=L/N$⁠. This solution is also pursued numerically.

The method of Jacobi relaxation was used to find the steady state solution of the force-free model, and forward Euler time stepping was used to find the quasi-steady state solution of the non-force-free model.⁴⁶ In both cases a logarithmically spaced grid was employed for the nonvelocity variable, so the magnetic helicity or guide field flux conserving source term needed to make use of interpolation. The value of $f(K−K′)$ (or $f(ϕ−ϕ′)$⁠) was found by using quadratic interpolation to $K−K′$ (⁠$ϕ−ϕ′$⁠), while $f(K′)$ (⁠$f(ϕ′)$⁠) was simply evaluated at the grid points. All K (⁠$ϕ$⁠) integrations were performed with Simpson's rule adapted for the logarithmic grid. The v integrations were performed using trapezoidal quadrature on a uniform grid. All integration source terms converge to at least second order. The convective gradient terms were calculated with Euler upwinding. To avoid roundoff error, the magnetic helicity solution was found as a function of $K1/3∼ψ$ which significantly lowered the number of decades spanned by the domain. When using Jacobi iteration the non-force-free solution suffered from unstable over-relaxation unless a very high precision was used, leading to the choice to iterate in time until a quasi-steady state was reached. In both cases, the lower bound was fixed to determine ξ and an upper bound was chosen equal to 0. The solutions are given assuming $Bg=Brec=L=vA=1$ (where B_rec is used for v_A), except for the strong guide field regime where $Bg=100Brec$ is used.

The distribution as a function of magnetic helicity for the velocity-independent model was calculated for approximately $S=106$⁠, which assumes $N∼S/Sc$⁠. It is important to stress that in practice when solving for the distributions, instead of choosing S we choose the left bound, which determines the number of plasmoids, which in turn is used to estimate S. A plot of the distribution can be found in Fig. 2 for N = 97.2. The helicity dependent solution can be converted to a solution as a function of the reconnecting field flux by using $ψ=(K/α(ζ/B0))1/3$ (equivalent to assuming all plasmoids are always force-free), and enforcing $∫fψ(ψ)dψ=∫fK(K)dK=N$ where $dK=3(αζ/B0)ψ2dψ$ (⁠$fψ$⁠, f_K, and $fϕ$ are simply f for each variable with subscripts to emphasize their respective normalizing scale factors). The similarities between the force-free and the non-force-free models at $S=106$ are apparent as shown in Fig. 3. The exponential tail has been incrementally stretched to higher ψ for the force-free model, and even more subtle is the lessening of the power law index from 2 to 7/4. These changes may in practice be hard to detect from an experimental or numerical standpoint. The flux distributions for multiple Lundquist numbers are shown in Fig. 4 for the force-free model. As expected, at higher Lundquist numbers, the power law region is extended similarly to the behavior of the non-force-free solution.²⁷ 

Further differences are revealed in the comparison between the force-free and non-force-free distributions as a function of $ϕ$⁠. Figure 5 shows both models for a similar number of plasmoids (⁠$Nff∼100,Nnff∼60$⁠) with significant differences between the assumption of magnetic helicity conservation and relaxation, vs guide field flux conservation. The force-free model ties the $ϕ$ distribution to the ψ distribution, so the power law index of 7/4 in ψ produces a power law index of 11/8 in $ϕ$ because of Eq. (6) (⁠$B0ϕ=ζψ2$⁠) and the requirement of consistent normalization (⁠$∫fψ(ψ)dψ=∫fϕ(ϕ)dϕ=N$⁠). Additionally, the nonconstant $ϕ$ growth rate of the force-free model results in a second low-$ϕ$ power law with index 1/2. In the nonforce-free case, a steeper power law index of 3/2 is obtained, with an approximately constant initial region. Both models begin to exhibit exponential decay around roughly the same value of $ϕ$⁠.

Using the numerical methods mentioned above, the solutions to the distribution of plasmoids with varying velocities were also found. For the Gaussian velocity distribution calculated in a domain of $ψ∈[0,0.1]$ and $v∈[−3,3]$ with $f(0,0)∼106$⁠, the solutions are shown in Figs. 6 and 7.

Figure 6 displays the solutions after translation to ψ using the force-free relations, as well as after integration in v. It is plotted alongside the integrated non-force-free distribution for comparison. The result of the inclusion of Taylor relaxation in this model differs from that of the velocity-independent model. The exponential tail is once again extended slightly and the transition to the power law region is largely unchanged, but the power law indices are both 1. Note that using $N∼S/Sc$⁠, the Sweet–Parker thickness δ can be found (⁠$δ∼l/Sc∼LSc/S$⁠) as $4×10−5$⁠. The average value of the magnetic field magnitude, B, for the Taylor profile is $0.575Brec$⁠, so plasmoid widths reach the order of magnitude of the current sheet thickness at $ψ=2.32×10−5$ when B scales linearly with ψ. This is relatively close to the constant premerger region, meaning the inability of small plasmoids to fill the current sheet may not play a significant role. Furthermore, plasmoids slightly below the current sheet width may already be sufficiently large so as to exhibit the essential properties sought for replicating the force-free behavior predicted here. Assuming the ψ at which a given plasmoid first merges with another is proportional to $l=L/N=LSc/S$⁠, the competing processes both scale as $1/S$ and hence this balance should not change for different S.

A two-dimensional distribution of the solution in ψ and v is shown in Fig. 7. The power law region begins where the approximately constant distribution begins to narrow, indicating that collisions occur frequently in this region. The narrowing in the velocity spread continues until a near delta-function distribution results. The comparison of guide field distributions in the force-free and non-force-free cases is shown in Fig. 8. Again the magnetic helicity conserving force-free distribution presents a less steep power law than the guide field flux conserving non-force-free distribution at higher $ϕ$⁠.

In all models shown as functions of ψ, one can immediately observe similarities between the numerical solutions obtained here and those of Huang et al.²⁷ The initial region appears approximately constant, which can be understood as most plasmoids growing due to reconnection prior to experiencing a single merger. If $dψ/dt=γ=0.01$ and we assume the average time between mergers is $τA/N$ which for 100 plasmoids is 0.01, then on average a plasmoid will not merge until they reach the size at $ψ≈10−4$⁠. Regarding the exponential decay, this occurs wherever the dominant loss mechanism becomes advection, i.e., the remaining number of plasmoids above a certain ψ, $n>$ (the integral of f from ψ to $∞$⁠), is on the order of one. This can be seen in any of the distributions by estimating $n>∼f(ψ)Δψ$ with $ψ∼Δψ$ near the start of the exponential decay and explains why steeper distributions transition to exponential decay at smaller ψ.

As can be seen in Fig. 3, the velocity-independent force-free distribution differs only slightly from the non-force-free distribution. Of these differences in the force-free distribution, the most visible is a weak reduction in slope, and an extension of the exponential tail. The change in slope can be understood as an effect of flux conversion from guide field to reconnecting as a result of Taylor relaxation conserving magnetic helicity. Guide field fluxes add during a merger while reconnecting field fluxes do not. To satisfy the relationship $ϕ=ζψ2/B0$ [Eq. (6)] for the force-free field after relaxation while conserving magnetic helicity, the excess guide field flux must be converted to reconnecting field flux. In practice, using these relationships one finds that at most 20.6% of the guide field flux will be converted after a merger between equivalent flux plasmoids. Given that smaller plasmoids are much more populous, mergers often result in little flux conversion, and hence the slope change in f is subtle. This can be understood more quantitatively by a comparison between the merger terms in the force-free and non-force-free statistical equations. These terms side by side can be rewritten for ψ as (dropping the τ_A)

where a scaling factor is present in the force-free expression as a result of the transition from K to ψ. Due to the monotonic decreasing nature of f, $f((ψ3−ψ′3)1/3)≥f(ψ)$ over the bounds of integration, increasing the integrand, and the overall value of the source in Eq. (18). This is further aided by the scaling factor, which within the bounds of the integral reaches a maximum of $22/3$⁠. The integration interval is slightly shortened by the upper bound of $ψ/21/3$⁠; however, due to the steepness of f and the scaling factor, the latter effect is outweighed by the former. In the case of the velocity-dependent model, the slope is far less steep and the collision rate is weighted by differential velocity (see Fig. 6). This more mild slope lowers the enhancement of $F((v,ψ3−ψ′3)1/3)$ in the integral, allowing it to be muted by an effect introduced with the inclusion of a velocity distribution. A comparison of the $ψ−v$ dependence of the non-force-free distribution with the force-free distribution in Fig. 9 shows that in the power law regime, the force-free distribution experienced an outward shift of plasmoids from the central low-v region to higher v regions. Overall, the average speed of plasmoids increased by 1.5%. While this change is small, this causes plasmoid mergers to occur more frequently in the power law region, increasing the steepness of F slightly. This may be an effect of flux conversion leading to more rapidly moving low flux plasmoids “leapfrogging” into the merger-dominated regime and shortly merging with other plasmoids. While this bolstered collisionality may increase the steepness, any merger resulting in flux conversion still boosts the ψ of one of the plasmoids. That is why in both the velocity independent and dependent force-free models, the exponential tail of the distribution is extended, albeit to a lesser extent in the velocity-dependent model.

From Fig. 5, the differences in the velocity-independent guide field distributions are much more pronounced than those of ψ. While one may expect flux conversion of guide field flux to reconnecting flux would cause the force-free model to have a steeper slope than the non-force-free model (at least without velocity), the results are quite the opposite. This is due primarily to the reconnection growth rates of $ϕ$ in the two models. In the non-force-free case the growth rate is the constant γ_g for all plasmoids, however, the relaxing model does not grow $ϕ$ at a constant rate. The nature of the guide field flux in this model is assumed to be passive, and a plasmoid's ψ grows in a manner that maintains a force-free state at the expense of the growth of guide field flux. Therefore, the plasmoid follows $ϕ=ζψ2/B0$ [Eq. (6)] so that $dϕ/dt=2γB0ϕ/ζ$⁠, growing faster the larger it gets. The force-free growth rate, although smaller than γ_g initially, surpasses its non-force-free counterpart near $ϕ≈10−6$⁠. This nonconstant growth rate leads to another distinguishing feature of the force-free model, a 1/2 power law index at low $ϕ$⁠. This two vs one power law difference between models is much more significant than that which occurs later in $ϕ$ and could significantly ease the difficulty in determining which model is most appropriate for a given data set. The same effects are seen in the velocity-dependent distribution comparison, although in this case, the non-force-free power law is the same as the velocity-independent version, even though the force-free power law has changed. This relates to the increase in steepness observed in the velocity-dependent force-free ψ distribution. Since $ϕ∼K2/3$ in the force-free model, mergers add $ϕ$ in a sort of two-thirds quadrature. However, directly adding $ϕ$ in the non-force-free model allows for a larger “jump” when merging, leading to a more significant difference in slopes.

In Sec. II B, the force-free model was chosen to be the limit where all plasmoids are in a force-free state at all times, even when merging. However, it is possible that a plasmoid that undergoes two rapid consecutive mergers may not reach a force-free state in between the mergers, even though magnetic helicity conservation is unaffected. Figure 10 shows the rate of change of a plasmoid's magnetic helicity, normalized to its helicity right before merging as a function of $ϕ$⁠, solely due to mergers for the force-free models. This reveals wherein the distribution a given plasmoid experiences the most frequent and significant mergers. The velocity-dependent model is shown evaluated at v = 0. As an exercise, if one considers relaxation to be possible when the helicity changes by no more than the given K in an Alfvén time, the v-independent model would exhibit force-free plasmoid configurations before $ϕ≈10−9$ and after $ϕ≈10−3$⁠, while slower plasmoids in the v-dependent model would be force-free after $ϕ≈10−4$ as well as before $ϕ≈10−9$⁠. The drop-off at high $ϕ$ indicates that the largest plasmoids in the exponential tail should be frequently found in a relaxed state. However, this is simply meant to illustrate wherein the distribution plasmoids undergo the most change. If relaxation can take place sufficiently rapidly, then one will observe characteristics of the force-free model in the region of rapid helicity change. Finally, we note that we assume that mergers of 3 or more plasmoids are rare enough so that their effects on the distribution are negligible.

These models do not necessarily represent the exact distributions of plasmoids in a current sheet that undergoes dynamic reconnection. Due to the disruption of relaxation by frequent mergers at some scales, it may be more appropriate to treat the non-force-free and force-free cases as asymptotic behaviors. From Fig. 10, we expect to observe force-free plasmoids at the smaller scales, as well as the largest. The suspicion that large ψ plasmoids may be able to relax is hinted at by the observation of Taylor-like relaxed flux ropes present in the solar wind as well as the earth's magnetotail.^47,48 While in the Hall-MHD regime, these results may still indicate that larger plasmoids are capable of achieving a force-free state, suggestive of some processes dominated by magnetic helicity conservation.

The results here indicate that guide field flux distributions may be noticeably distinct, especially due to the introduction of a second low-$ϕ$1/2 power law index in the force-free model. Experimentally, however, this can be a challenging quantity to measure. A more directly measurable characteristic would be the distribution of plasmoid sizes, or the quantity a in the force-free case. From the surface integration of the reconnecting field, we can easily demonstrate that $ψ∼a$ in the force-free case. Unfortunately, the distinction of this from the non-force-free model may also prove difficult. Using the conclusions of Uzdensky et al.,²⁶ in the incompressible limit of the antiparallel plasmoid instability the width of a plasmoid perpendicular to the current sheet would follow $wp∼ψ/Brec$⁠. This incompressible limit may also be thought of as the case of a strong guide field where the added magnetic pressure augments the pressure of the magnetized plasma within. Hence one may only be able to distinguish the physical size distributions, as well as the ψ distributions, can be distinguished between the force-free and non-force-free models. Instead, a more accessible indication that helicity conservation is, to some extent, playing a role may simply be the presence of an enhanced core field in line with that of the Taylor profiles internal to the plasmoids which form the backbone of the force-free model.^31,33,34,47 It should be noted that if the ψ distributions (or more importantly a or w_p distributions) are similar between the force-free and non-force-free models and the $ϕ$ distributions vary largely, then the internal plasmoid guide field pressures will also be significantly different. Here, guide field pressure would reach noticeably greater levels in the force-free model over the non-force-free model.

The models presented here may also prove relevant to the plasmoid distribution when the external guide field is negligible. Instabilities which arise in three dimensions such as kink modes add an out-of-plane magnetic field to an otherwise antiparallel configuration. The process of flux rope merging in 3D is a particularly important example of an event where kinking occurs in an otherwise straight cylindrical island in 2D. The condition of thin flux ropes being necessary should be mentioned to ensure that the characteristic length scale of a kink does not largely affect the assumption of approximately straight cylindrical Taylor states (as in high aspect ratio reverse field pinch devices). This negligible guide field regime is not forbidden in any way by our equations. Rather, it is the case where the quantity λa is the first zero of $J0(x)$⁠. This phenomenon of Taylor relaxation where the external guide field is zero has been observed experimentally.⁴⁹ 

The addition of relaxation has also allowed for an increase in the spread of plasmoid velocities. Although it is a small change of 1.5%, a Fermi-like process involving multiple reflections from plasmoids can add up to cause an increase in the highest achievable energy of charged particles undergoing this acceleration. More significantly, the acceleration which occurs during island contraction relies on how quickly the island is compressed. The enhancement of guide field pressures in the relaxing model and the increasing rate with which guide field flux is accreted onto a plasmoid can result in a greater effective mirror velocity in a first order process.⁵⁰ 

The force-free model's alterations to particle acceleration processes and the speed of plasmoids in the outflow are unlikely to affect the global rate of energy conversion in any considerable capacity. While energization may be enhanced for some high energy particles already capable of undergoing Fermi acceleration, they represent a minute fraction of the overall energy of the plasma and therefore are not likely to contribute to the global energy conversion rate in a significant manner. It is possible, however, that the dissipation of magnetic energy during the Taylor relaxation that accompanies plasmoid growth and mergers could result in a modified conversion rate. The assumption of perpetual relaxation as a plasmoid grows due to reconnection leaves the reconnecting flux growth rates the same, but the guide field flux growth rates differ significantly between models. In addition, it can be shown that a force-free plasmoid's total magnetic energy goes linearly with its guide field flux. For large plasmoids, the guide field flux growth rate is enhanced in the force-free model, but for small plasmoids, it is significantly hindered. Because of this, it is not clear whether the overall rate of energy conversion would increase or decrease between models. Therefore, we leave this calculation to future investigations.

In order to make a more robust model of the plasmoid distribution in a guide field, several more effects may be included in the future. Specifically, the assumption of merged plasmoids adopting the velocity of the larger allows for the narrowing of the velocity-dependent distribution to a near delta function in v. If instead flux weighted averaging was used, a greater spread in velocity at higher ψ would result. This broadening is also suggested by the spread in the velocity distribution found by Lingam and Comisso.²⁵ In parallel, a more thorough investigation of the effects of the choice of H(v) on the distribution of plasmoid velocities may be of interest. While this work has focused on the distribution in fluxes and helicity, more realistic merging rules may be accompanied by varied effects on the distribution of plasmoid velocities. It may also prove important to address the possibility of incomplete or particularly slow plasmoid mergers. Coalescence rates have been observed to stall in simulations due to a sloshing effect at high Lundquist numbers.⁵¹ This delayed or inhibited coalescence could give rise to more efficient production of high flux plasmoids in our statistical models. Additionally, there have been growth rates proposed for the plasmoid instability which would modify the Sweet–Parker reconnection rate γ used here. While they may not affect the power law slope, an increase or decrease in γ would affect where the transition from the constant region to the power law region occurs. Numerical or experimental investigations of these distributions may seek to determine which model is most appropriate for the problem of plasmoid unstable reconnection in a guide field, or whether they both remain contained within their expected regimes. Observation of power law behavior outside the limits of the force-free and non-force-free models may suggest that conservation of an alternative quantity holds over guide field flux or magnetic helicity, while something in between may, if appropriate, suggest an intermediate regime that possesses characteristics of both distributions.

This work is supported by the DOE through Contract No. DE-AC0209CH11466. The authors would like to thank Yi-Min Huang for his help with numerical methods, careful verification of results, and correction of mistakes, as well as William Fox for his helpful feedback and Amitava Bhattacharjee for his valuable conceptual discussions.

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

A required characteristic of our plasmoid kinetic equation is that the terms involving plasmoid mergers must conserve magnetic helicity in the distribution. In other words, the first moment of the following equation, from Eq. (19), must be zero:

In this section, we will prove that this is true. We begin by integrating the equation over the entire velocity space v. For readability, we will use the shorthand $∫v=∫−∞∞dv$ and $∫K=∫0∞dK$⁠, as well as $Δv̂=|v−v′|/VA$⁠,

Some reshuffling of the final term will prove useful. First, a substitution of $x=K−K′$ yields

Given that v, $v′$⁠, and x are integrated out, we can rename them as $v′$⁠, v, and $K′$⁠, respectively, without affecting the results

Therefore, the identical integrand on both side of the above equation permits the following relationship:

A property of our distributions is that $F(K,v)=0$ for K < 0 and equivalently $F(K−K′,v)=0$ for $K′>K$⁠. Therefore, we may extend the K integration out to infinity producing a convolution, and replace the result with our integrated merger terms

where we use the standard notation for a convolution $F(v′)*F(v)$ with explicit v dependence to differentiate between the primed and unprimed variables. The next step is to Laplace transform this equation from K to s, and take a derivative with respect to s

Taking the negative of Eq. (A7) and evaluating it at s = 0,

This proof can easily be generalized to the velocity-independent case, and is also of course applicable to the guide field conserving regime with K replaced with $ϕ$⁠.