Do chaotic field lines cause fast reconnection in coronal loops?

Magnetic reconnection is a fundamental mechanism that changes the topology of magnetic field lines and converts the magnetic energy to the plasma thermal and non-thermal energy.^1–7 It is generally believed that this mechanism drives explosive phenomena in astrophysical, space, and laboratory plasmas, including solar flares, coronal mass ejections, geomagnetic substorms, and sawtooth crashes in fusion devices.

Magnetic reconnection can be classified into two-dimensional (2D) and three-dimensional (3D) reconnection. Real-world magnetic reconnection takes place in 3D, but 2D reconnection is commonly employed as an approximation by assuming that the whole process depends only on two spatial coordinates. Two-dimensional reconnection occurs at an X-point (or X-line when extending along the direction of symmetry) where separatrices, which separate topologically different magnetic field lines, intersect. When magnetic stress builds up at the X-point prior to reconnection, a thin current sheet forms. Then, during reconnection, the field-line velocity (i.e., a velocity field that carries field lines from one time to another) diverges at the X-point.⁸ In other words, the field line is cut and rejoined with another field line at the X-point.

Compared with 2D problems, magnetic reconnection in 3D remains a conceptual challenge, especially when topological structures, such as magnetic null points and closed magnetic field lines, are absent.⁹ In this situation, all magnetic field lines are topologically equivalent, and a continuous velocity field that preserves the magnetic field line connectivity can always be found.¹⁰ That raises a fundamental question: how and where does magnetic reconnection occur if all field lines are topologically identical? Several ideas have been proposed to address this question, including the general magnetic reconnection theory^11–13 that uses parallel voltage as a metric for 3D reconnection rate and the concept of quasi-separatrix-layers (QSLs) defined as regions with high squashing factors of the field-line mapping.^14–16 Even though 3D reconnection is a vast and ongoing research topic, most theories share one common aspect: just like in 2D reconnection, intense thin current sheets play a critical role in 3D reconnection.

Through a series of publications over the past decade, Boozer has advocated a paradigm shift regarding 3D reconnection.^17–24 The gist of Boozer's proposal can be described as follows. In a 3D magnetic field, neighboring magnetic field lines generically exponentiate away from each other. The field-line flow, which is Hamiltonian, becomes chaotic. If we follow two field lines initially separated by an infinitesimal distance $δr(0)$⁠, in most cases, the separation grows exponentially as $δr(ℓ)=eσ(ℓ)δr(0)$⁠, where $ℓ$ is the distance along the field line, and $σ(ℓ)$ is an overall (but in general, non-monotonically) increasing function over distance. Boozer has argued that under the condition of large field-line exponentiation, an exponentially small non-ideal effect will completely scramble the field-line mapping, leading to fast reconnection without the necessity of intense thin current sheets. Indeed, that fast reconnection can occur without intense thin current sheets is what sets Boozer's theory apart from the “traditional” theories, to use Boozer's terminology.²⁴ 

Recently, Boozer and Elder²⁵ have proposed a simple coronal loop model to test this new paradigm. In their model, the coronal loop is enclosed in a perfectly conducting cylinder of a radius a and a finite length $0≤z≤L$⁠. The initial magnetic field is uniform and pointing along the $ẑ$ direction. The magnetic field lines are line-tied to the top and the bottom boundaries, which represent the photosphere. On the top boundary at z = L, a time-dependent flow is imposed; all other boundaries are static. The imposed boundary flow mimics photospheric convection and gradually entangles (or “braids”) the field lines in the coronal loop. The braided magnetic field lines eventually reconnect and release the stored magnetic energy into plasma kinetic energy and heat.

On the surface, the Boozer–Elder model is similar to Parker's model^26,27 of coronal heating, but their predictions for the scenarios of magnetic reconnection and energy release are starkly different. Parker's scenario predicts that intense thin current sheets will develop throughout the coronal loop as a consequence of field line braiding. These thin current sheets cause numerous small-scale reconnection events, or “nanoflares” (as Parker called them), that heat the solar corona to millions of degrees. Thin current sheets play a crucial role in Parker's scenario of coronal heating. In fact, Parker argued that the thin current sheets will become singular (i.e., take the form of Dirac $δ$-functions) in the ideal-MHD limit when the resistivity vanishes.²⁸ Parker's prediction of ideal singular current sheets, often dubbed the “Parker problem,” has remained controversial for several decades, continuing to this day.^29–49 

Contrary to Parker's nanoflare scenario, intense thin current sheets play no significant role in the Boozer–Elder scenario. While Boozer and Elder also predict that the electric current distribution will form thin ribbons, the current density does not become very intense and increases only linearly in time. In their scenario, separation of neighboring field lines, which increases exponentially in time, plays a dominant role in triggering the onset of fast reconnection.

The Boozer–Elder model is a welcome new development of Boozer's theory because it provides concrete and testable predictions. The primary objective of this study is to test some of the predictions. Specifically, we will focus on two distinct predictions, one related to ideal evolution and the other to resistive evolution. For the ideal evolution, the model predicts that the current density will increase linearly in time, while the separation of neighboring field lines will grow exponentially in time. With a small but finite resistivity, because the exponential field-line separation amplifies the field-line velocity, thereby speeding up reconnection, the timescale for the onset of fast reconnection and, therefore, the current density will scale logarithmically with the Lundquist number.⁵⁰ This logarithmic scaling relation for the current density is perhaps the most striking difference of Boozer's theory compared to traditional reconnection theories.

This paper is organized as follows. Section II outlines the reduced magnetohydrodynamic (RMHD) model and the imposed footpoint motions. Section III lays out the governing equations for the ideal and resistive quasi-static evolution of a coronal loop driven by footpoint motions. In Sec. IV, we test Boozer and Elder's current density calculation with an ideal RMHD simulation. In Sec. V, we present resistive simulations to test the scaling of current density with the Lundquist number and investigate whether chaotic field-line separation causes the onset of reconnection. We conclude in Sec. VI.

We employ the standard reduced magnetohydrodynamics (RMHD) model^29,51,52 in this study. The RMHD model assumes a strong uniform guide field along the z direction and that spatial scales along the guide field are much longer than that in perpendicular directions. Under these assumptions, the MHD equations can be simplified to a set of two equations

Here, we normalize the strength of the guide field to unity. The operator $∇⊥≡x̂∂x+ŷ∂y$ denotes the gradient operator in the perpendicular directions of the guide field. The magnetic field is expressed in terms of the flux function A through the relation $B=ẑ+∇⊥A×ẑ$⁠. The plasma velocity u is expressed in terms of the stream function $ϕ$ as $u=∇⊥ϕ×ẑ$⁠. The vorticity and the electric current density along the z direction are given by $Ω≡−∇⊥2ϕ$ and $J≡−∇⊥2A$⁠, respectively. The Poisson bracket is defined as $[ f,g]=∂yf∂xg−∂xf∂yg$⁠. Dissipation is introduced by including the resistivity η, the viscosity ν, and a friction coefficient λ.

The RMHD model is widely used in analytic and numerical studies of Parker's coronal heating model^{29,33,53–62} and also in the study of Boozer and Elder.²⁵ From the definitions of B and u and the Poisson bracket, we have the following useful relations:

and

for an arbitrary variable f.

The RMHD equations are solved with the DEBSRX code,⁶³ which is a reduced version of the compressible MHD code DEBS.⁶⁴ The x–y plane is discretized with a Fourier pseudospectral method. The z direction is discretized with a finite-difference scheme where $ϕ$ and A reside on staggered grids. The timestepping scheme is a semi-implicit, predictor-corrector leapfrog method, where $ϕ$ and A are staggered in time. While the semi-implicit scheme allows time steps to be larger than the limit set by the Courant–Friedrichs–Lewy (CFL) condition,⁶⁵ too large a time step may compromise accuracy. For this reason, we have been conservative in setting the time step. We determine the time step dynamically by using the CFL condition of the Alfvén wave along the guide field and that of the perpendicular flow speeds, whichever gives a smaller time step. We have tested the accuracy of this choice by reducing the time steps by a factor of two for selected runs and have seen no significant difference.

We assume that the system is bounded in the z direction by two conducting plates at z = 0 and z = L. The x–y plane is a 1 × 1 box (⁠$0≤x,y≤1$⁠) with doubly periodic boundary conditions. We take L = 10 in this study. The bottom boundary at z = 0 is stationary for all time; i.e., we impose the boundary condition $ϕ=ϕb=0$⁠. At the top boundary of the simulation box (z = L), we impose a time-dependent flow given by the stream function

where c_i and ω_i are constants. This boundary flow is similar, but not identical, to the flow employed by Boozer and Elder.²⁵ Because the DEBSRX code is limited to doubly-periodic systems in the perpendicular directions whereas Boozer and Elder consider a cylindrical domain, we cannot impose the same boundary conditions as they do. The difference in the boundary conditions, however, is not germane to the issue we will discuss. Note that the boundary flow is localized to the middle and vanishes at the edges of the simulation box in the perpendicular directions. This feature is qualitatively similar to the flow employed by Boozer and Elder.

We start the simulations with only the guide field. The imposed boundary flow drags the footpoints and gradually entangles the field lines. From solar observations, it is known that the timescales of footpoint motions are much longer compared with the Alfvén transit time along the coronal loop. We set the constants to $c0=0, c1=c2=c3=0.000 25$⁠, $ω1=0.01252, ω2=0.01253$⁠, and $ω3=0.01255$⁠. For this set of parameters, the maximal footpoint speed is typically of the order of 0.01, corresponding to an advection timescale on the order of 100. In contrast, the Alfvén transit time along the coronal loop is 10. Therefore, the condition for the separation of timescales is met. We also make the ratios between frequencies $ωi/ωj$ irrational so that the flow pattern will not repeat itself, but this choice is not crucial for the purpose of this study.

Under the condition that the characteristic timescales of footpoint motions are much longer than the Alfvén transit time, the coronal loop will evolve quasi-statically, provided that the coronal loop is not close to an instability threshold. For a quasi-static coronal loop, we may assume that the plasma inertia as well as the viscous and frictional force are negligible; therefore, the magnetic force-free condition (within the RMHD approximation)

is satisfied for all time. The governing equations for the quasi-static evolution now are the force-free condition (6) together with the induction equation (2). In this set of equations, the time derivative appears only in the induction equation (2), whereas the force-free condition plays a role analogous to the incompressibility constraint, $∇·u=0$⁠.

The quasi-static evolution of coronal loops driven by slow footpoint motions can be determined as follows. First, taking the time derivative $∂t$ of Eq. (6) yields

We can obtain an equation for $∂tJ$ by applying the $−∇⊥2$ operator on Eq. (2), yielding

Now, we can use Eqs. (2) and (8) to eliminate the time derivatives in Eq. (7) and obtain the following equation for $ϕ$⁠:

where the linear operator $L$ is defined as

At a given instant, if the operator $L$ is invertible subject to the boundary conditions $ϕ|z=0=ϕb$ and $ϕ|z=L=ϕt$⁠, we can in principle obtain the stream function $ϕ$ (and the flow $u=∇⊥ϕ×ẑ$⁠) that will carry the system to another force-free equilibrium at the next instant. In the ideal limit $η→0$⁠, Eq. (9) is identical to the equation derived by van Ballegooijen in his study of the Parker problem.^29,33

The operator $L$ is the one that appears, unsurprisingly, also in the ideal linear stability problem of the instantaneous equilibrium, which takes the following form:

Here, a $ϕ∼eγt$ time-dependence and homogeneous boundary conditions $ϕ|z=0=0$ and $ϕ|z=L=0$ are assumed. The operators $L$ and $∇⊥2$ in Eq. (11) are self-adjoint; therefore, the eigenvalues $γ2$ are real, and the eigenfunctions form a complete set; functions $ϕm$ and $ϕn$ of different eigenvalues satisfy the orthogonal condition

while those with the same eigenvalues can be orthogonalized as well.⁶⁶ 

Suppose the complete set of eigenfunctions and eigenvalues ${ϕn,γn2}$ is known. We can formally solve Eq. (9) by making an eigenfunction expansion⁶⁷ 

Here, $ϕ̃$ is an arbitrary smooth function satisfying the inhomogeneous boundary conditions $ϕ|z=0=ϕb$ and $ϕ|z=L=ϕt$ [e.g., $ϕ̃=ϕb+(ϕt−ϕb)z/L$]. Using the orthogonality of eigenfunctions, we can determine the coefficients a_n as

From Eq. (14), we conclude that Eq. (9) is solvable provided that the system is not at marginal stability, i.e., none of the eigenvalues $γn2$ vanish.

The preceding discussion shows that to determine the quasi-static evolution of this model requires solving Eq. (9) in the whole domain in each time step. On the other hand, Boozer and Elder²⁵ take a different approach for the quasi-static ideal evolution that only involves solving an equation for the current density at the top boundary. Their approach is as follows. First, using Eqs. (3) and (4), and the definition of the Poisson bracket, the term $−∇⊥2(B·∇ϕ)$ in Eq. (8) can be written as

Boozer and Elder ignore the terms $−∂yB⊥·∂xu+∂xB⊥·∂yu$⁠, whereupon Eq. (8) becomes (here, we set η = 0 for the ideal evolution)

Because of the force-free condition $B·∇J=0$ and the ideal frozen-in condition, the left-hand side of Eq. (16) is constant along a field line. Consequently, the right-hand side $B·∇Ω$ is also constant along a field line. For the Boozer–Elder model, the vorticity at the bottom boundary vanishes, and the vorticity at the top boundary is prescribed. Therefore, along a field line, $B·∇Ω=Ω|z=L/L$ at any instant. Boozer and Elder then consider Eq. (16) at z = L, yielding

This equation (hereafter the BE equation) plays a critical role in the study of Boozer and Elder.⁶⁸ 

The BE equation takes the form of an advection equation for the current density J on the left-hand side, whereas the right-hand side serves as a source term. Because the flow velocity u (and, therefore, the vorticity Ω) is prescribed at the top boundary, the BE equation can be solved to determine the current density distribution at the top boundary without further knowledge of what occurs in the remaining part of the system. That is not the case if the missing terms are included, because $B·∇Ω$ will no longer be a constant along a field line. In other words, it is not possible to amend the BE equation simply by including those terms.

The advection term on the left-hand side of the BE equation only redistributes the current density without changing its magnitude. The peak current density can increase only through the source term on the right-hand side. Hence, the current density is bounded by $|J|≤Ωmaxt/L$⁠, where $Ωmax$ is an upper bound of $|Ω|$ at the top boundary. In other words, the current density increases at most linearly in time according to the BE equation.

Because some terms are neglected when deriving the BE equation, the question now is: are those terms negligible? We address this question by comparing the solutions of the BE equation and that of the RMHD equations. The BE equation is solved using a pseudospectral method implemented with the Dedalus framework;⁶⁹ the grid resolution is 1024². The RMHD equations are solved with the DEBSRX code with a grid resolution 1024³. We set η = 0 in the DEBSRX simulation, but a small viscosity $ν=10−6$ is applied for numerical stability. The numerical algorithms of the Dedalus implementation for the BE equation and the DEBSRX code for the RMHD equations are similar. Both use a Fourier pseudospectral method dealiased by the Orszag two-thirds rule⁷⁰ in the x–y plane, so the numerical errors should be similar. Both simulations have been compared with numerical solutions at lower resolutions (256² and 512² for BE; 256³ and 512³ for RMHD) to ensure that the results presented here are well-resolved and converged.

Even though the timescales of footpoint motions are longer than the Alfvén transit time by approximately one order of magnitude, the RMHD calculation is not exactly force-free. To ensure that the RMHD solution remains close to force-free, we restart from a few selected snapshots, set the plasma speed to zero at footpoints and across the entire domain, and then turn on the friction force and let the system relax to a force-free equilibrium. This experiment shows that the RMHD solution remains approximately force-free up to t = 200, thereby ensuring a fair comparison between the BE and the RMHD solutions.

We now compare the RMHD solution at the top boundary with the BE solution and summarize the results in Fig. 1. Panel (a) shows the time histories of the maximum current density obtained from both sets of solutions. The maximum current densities from both solutions agree until t = 150 and then depart significantly afterward. Furthermore, the current density in the RMHD solution increases significantly faster than predicted by the BE equation. Panel (b) shows snapshots from both sets of solutions at two representative times, at t = 140 and 200. Although the two solutions give essentially the same maximum current density at t = 140, we can already see the differences between the two solutions. The differences become quite pronounced at t = 200. By that time, thin current sheets have developed in the RMHD solution but are absent in the BE solution.

Our results show that the neglected terms in the derivation of the BE equation are not negligible. Therefore, to determine the current density at the top boundary, we must solve for it over the entire 3D domain. Importantly, the BE equation significantly under-predicts the current density of the RMHD solution. As we will see in Sec. V, these intensifying current sheets eventually lead to the onset of reconnection.

We continue with the resistive evolution of the model problem. Our objectives are to test the prediction that the peak current density will scale logarithmically with respect to the Lundquist number S, as well as to assess whether the “exponential” field-line separation causes the onset of reconnection.

To address these questions, we perform a series of simulations with the Lundquist number S varying from 10⁴ to 10⁶. Here, the Lundquist number $S≡aVA/η$ is defined through the box size a in the perpendicular direction, the Alfvén speed V_A of the guide field, and resistivity η. In our normalized units, the box size a = 1 and the Alfvén speed V_A = 1; therefore, the Lundquist number is simply $S=1/η$⁠. The viscosity is another free parameter. For simplicity, we set the viscosity $ν=η$ for all cases; i.e., the magnetic Prandtl number $Pm≡ν/η=1$⁠. (See the  Appendix for a discussion of the effect of viscosity.) The friction coefficient λ is set to zero. Table I lists the Lundquist numbers and grid resolutions for all the simulations we have performed. The simulation time is t = 1000 for all cases, corresponding to 100 Alfvén transit times and approximately ten footpoint advection times.

Table I also shows the total resistive dissipation $Wη$ and viscous dissipation $Wν$ during the whole period of each simulation. Over the range of Lundquist numbers S that spans two orders of magnitude, the total dissipation $Wη+Wν$ stays remarkably close to constant. The dissipation is predominantly due to resistivity, although the portion of viscous dissipation slightly increases as the Lundquist number increases.

To ensure that our simulations have sufficient numerical resolution, we perform three simulations (D1–D3) for the case $S=4×105$ with resolutions ranging from 512³ to 1024³ and do not find significant difference between them. Therefore, a grid resolution of 512³ appears to be adequate for $S=4×105$⁠. This finding gives us some confidence that the highest Lundquist number case, run E with $S=106$⁠, should be reasonably resolved by a grid resolution of 1024³. In addition, we may also assess the accuracy of our numerical simulations by how precisely the energy is conserved. The energy conservation requires that the total Poynting energy input through the top boundary W_P should be equal to the sum of the increase in the magnetic energy E_M, the increase in the kinetic energy E_K, and the total dissipation. For all the cases, the error in energy conservation is less than 1% of the total Poynting energy input over the entire simulation period.

The column $J1/2$ in Table I indicates that regions with $|J|>J1/2$ contribute half of the resistive dissipation over the entire simulation period, and $V1/2/V$ is the time averaged volumetric ratio between the regions with $|J|>J1/2$ and the entire domain. These numbers show that resistive dissipation is increasingly concentrated in small regions of high current density as the Lundquist number increases. For example, approximately 5% of the total volume contributes one half of the resistive dissipation at $S=104$⁠, whereas less than 0.2% of the total volume accounts for the same portion at $S=106$⁠.

We plot the probability distributions $P(|J|)$ of the current density over the four-dimensional spacetime during the period $0≤t≤1000$ for various Lundquist numbers S in Fig. 2. The probability distribution exhibits a strong dependence on the Lundquist number. Furthermore, the maximal current density roughly scales linearly with S. This $J∼S$ scaling relation is stronger than the Sweet–Parker^71,72 scaling relation $J∼S1/2$ and is on par with that of plasmoid-mediated reconnection.⁷³ Therefore, the prediction that the current density will depend logarithmically on S appears to be inconsistent with our numerical findings.

We now address the question whether the exponential separation of the neighboring field lines causes the onset of fast reconnection by taking a closer look at the highest-S simulation, run E.

During the early phase of the simulation when $t≤190$⁠, the system evolves quasi-statically, while the current density gradually increases as in the ideal evolution. After t = 190, which is approximately twice the footpoint advection time, an onset of activity leads to fast plasma flow and further intensifies the thin current sheets. The peak plasma flow speed after the onset is comparable to the typical in-plane Alfvén speed and is faster than the typical footpoint speed by an order of magnitude. In contrast, the plasma speed in the interior is typically slower than or comparable to the typical footpoint speed during the quasi-static phase.

The flow pattern and current density distribution after the onset clearly show signatures of magnetic reconnection, such as outflow jets within current sheets, as can be seen in the 2D slice at t = 260 shown in Fig. 3 (Multimedia view). Figure 4 shows a 3D visualization of the entire domain as well as a zoom-in view of a sub-domain. From the zoom-in view in panel (b), we can see that the primary current sheet with J < 0 (the one with dark colors) has become unstable to the plasmoid (or tearing) instability and developed flux-rope-like structures despite the stabilizing effects of the line-tied boundary conditions. The plasmoid instability has continued to be present in numerous reconnecting current sheets throughout this simulation. Similar to the systems without line-tied boundary conditions, the Lundquist number needs to exceed a threshold for the plasmoid instability to occur.^73–76 In all other simulations with lower Lundquist numbers, the plasmoid instability has not been observed.

Previous studies indicate that the stabilizing effect of the line-tied boundary condition for the tearing mode is negligible when the “geometric” width $δgeo∼(λtearingBz/LBup)δ$ is thinner than the inner layer width of the tearing mode, $δtearing$⁠.^77–80 Here, δ is the current sheet thickness, $Bup$ is the upstream in-plane magnetic field of the current sheet, and $λtearing$ is the wavelength of the tearing mode. Taking the current sheet at t = 250, immediately prior to the onset of the plasmoid instability, we estimate $Bup≃0.1$ and $δ≃0.005$⁠. For this set of parameters, the fastest growing mode wavelength $λtearing≃0.15$ and the corresponding inner layer width $δtearing≃0.001$⁠,⁷⁵ whereas the geometric width $δgeo≃0.000 75$⁠. Therefore, the geometric width is comparable to the inner layer width, and the stabilizing effect of line-tying should be marginal. This estimate is consistent with the fact that the current sheet becomes unstable. In contrast, for cases of lower Lundquist numbers, the corresponding current sheets at the same time do not satisfy the criteria $δgeo<δtearing$ and remain stable.

To disentangle the relationship between reconnection and the exponential separation of the neighboring field lines, we have implemented a suite of diagnostics together with field-line tracing. Along each field line, we calculate (a) the squashing factor Q,^14,15,80 (b) the parallel voltage, (c) the plasma flow velocity, and (d) the velocity of the field line relative to the plasma. The squashing factor Q quantifies the extent of the neighboring field-line separation; the parallel voltage is often employed as a metric for 3D reconnection rate;^11,13 finally, a large separation between the field-line velocity and the plasma velocity may also indicate magnetic reconnection.

We calculate the squashing factor Q by simultaneously integrating the equation for the magnetic field-line flow

and the equation for an infinitesimal separation $δx⊥$ between the neighboring field lines

Equation (19) can be written in a matrix form as follows:

Because Eq. (20) is linear in $δx$⁠, it is sufficient to integrate it with respect to two linearly independent initial condition. For that purpose, we integrate the equation

with the initial condition

Then, for any initial separation $δx|z=0$ we can obtain the separation at an arbitrary z as

The singular value decomposition (SVD)⁸¹ of N(z) is of the form

where both U(z) and V(z) are unitary matrices. The singular values $λmax$ and $λmin$ have a geometrical interpretation as follows. If we follow a infinitesimally thin flux tube starting with a circular cross section of radius δr at z = 0, the cross section becomes an ellipse at z > 0, with the semi-major axis $δrmax=λmaxδr$ and the semi-minor-axis $δrmin=λminδr$⁠. Because the field-line mapping in RMHD preserves the area, the singular values $λmax$ and $λmin$ satisfy the relation $λmaxλmin=1$⁠. The squashing factor Q is then defined as

evaluated at the top plate. In the limit $Q≫1$⁠, the squashing factor is approximately the ratio between the semi-major and the semi-minor axes; i.e., $Q≃λmax/λmin$⁠.⁸² 

Next, we calculate the field-line velocity as follows. The electric field in our RMHD model is given by the resistive Ohm's law

If we can express the electric field by

for some velocity field v and a scalar potential $Φ$⁠, then the evolution of B is formally governed by the ideal equation $∂tB=∇×(v×B)$ such that the magnetic field lines are frozen-in to the velocity field $v$⁠. The velocity field $v$ is not uniquely determined because we can add to $v$ an arbitrary component parallel to the magnetic field B without changing Eq. (27). Moreover, taking the inner product of Eq. (27) and B yields

which can be integrated along the magnetic field lines to determine the potential $Φ$ up to a free function $Φ(x⊥)|z=0$ defined on the bottom boundary. To uniquely specify $v$⁠, we impose the conditions $vz=0$ and $Φ(x⊥)|z=0=0$⁠. The field-line velocity $v⊥$ is then determined by the relation

and the relative velocity between the field line and the plasma is given by

The imposed boundary condition $Φ(x⊥)|z=0=0$ makes the plasma velocity and the field-line velocity coincide at the bottom boundary. If the two velocities deviate significantly from each other as we follow the field lines, it may indicate that reconnection is ongoing.

To calculate the field-line velocity together with field-line tracing, we adopt a method similar to the calculation of Q. Applying the $∇⊥$ operator on Eq. (28) yields

The left-hand side of Eq. (31) corresponds to the variation in $∇⊥Ψ$ along the field lines. Using Eq. (30) to replace $∇⊥Ψ$ by $w⊥$⁠, we can rewrite Eq. (31) in a matrix form

Here, the derivative d/dz on the left-hand side is a derivative along the magnetic field lines. We can now integrate Eq. (32) together with field-line tracing to obtain the field-line velocity.

Figures 5–7 show the results of the diagnostics at three representative times. In each figure, the upper-left panel shows the time histories of the maximum current density $Jmax$ and the maximum plasma speed $umax$⁠, with the vertical line indicating the time of the snapshot. We label the field lines by their footpoints at the bottom boundary. The lower-left panel shows the projected image of 5 × 5 squares at the top boundary to the bottom boundary following the field lines. The remaining four panels show the squashing factor Q, the parallel voltage $∫0LηJ dz$⁠, the maximum plasma speed $|u|$⁠, and the maximum field line relative speed $|w|$ along each field line. At each time slice, we trace 1000 × 1000 field lines uniformly distributed at the bottom boundary. Additionally, we provide the entire time history of these diagnostics for run E, together with that for runs A, B, C, and D3 as animations [see Fig. 7 (Multimedia view) and the supplementary material].

Figure 5 shows that the maximum Q is approaching 10⁵ at t = 190. The “slippage” speed between the field line and the plasma also become approximately one order of magnitude larger than the plasma speed for field lines with high Q values. The maximum slippage of footpoint mapping relative to that of the ideal run is approximately 30% of the simulation box size in the perpendicular direction. Although one might interpret that significant magnetic reconnection has already occurred based on this information, it may be more appropriate to attribute the footpoint slippage to diffusion rather than reconnection because the evolution remains close to quasi-static up to this time.

After t = 190, there is a rapid onset of activity with the current density and plasma speed both increasing substantially. Notably, the onset occurs in approximately two footpoint advection times, much sooner than in ten advection times predicted by Boozer and Elder. Subsequently, at t = 290 shown in Fig. 6, the plasma speed is approximately one order of magnitude higher than that during the quasi-static phase, the squashing factor Q goes above 10¹⁰, and the field line relative speed is above 10², which is more than three orders of magnitude higher than the plasma speed. The correlations between the squashing factor Q, the parallel voltage, and the plasma relative speed are also evident. The same general features are also present in Fig. 7 for t = 620. At this time, the squashing factor Q goes beyond 10¹⁵ at some locations, and the maximal field line relative speed is above 10⁴, more than five orders of magnitude higher than the plasma speed.

Based on the result that both the squashing factor and the field-line speed become extremely high, would one conclude that the chaotic field-line separation causes fast reconnection? Upon close examination, our simulation results do not appear to support this viewpoint. If the chaotic field-line separation is the cause of fast reconnection, we expect the squashing factor to reach a maximum when the coronal loop “breaks loose”; i.e., when it starts to deviate from the quasi-static evolution; subsequently, reconnection will simplify the field-line mapping, and the squashing factor will decrease. That is not what the simulation shows. As we can see from the time sequence shown in Figs. 5–7 and the associated animation, after the coronal loop breaks loose after t = 190, thin current sheets intensify, and the squashing factor increases tremendously around them. Unlike Boozer's claim that the chaotic field-line separation can speed up reconnection without intense current sheets, the simulation shows that current sheets must intensify to release the built-up magnetic stress; the intensified current sheets further enhance the chaotic field-line separation as reflected in an even higher squashing factor Q.⁸³ 

For a complicated 3D evolution with reconnection occurring at multiple locations and at different times, it is difficult, if not impossible, to quantify the speed of the reconnection process with a single reconnection rate. However, we may quantify the effectiveness of reconnection through its effect of converting the magnetic energy into plasma heating through dissipation. Because the dissipation rate is nearly independent of the Lundquist number S, it is not unreasonable to think that reconnection proceed approximately at the same rate for different S. This conclusion is consistent with the fact that the peak values of the parallel voltage, often employed as a metric for 3D reconnection rate, mostly fluctuate between $2×10−4$ and $6×10−4$⁠, regardless of the Lundquist number. In contrast, the field-line speed exhibits a strong dependence on the Lundquist number S. At $S=104$⁠, the field line relative speed stays below $10−2$⁠, whereas at $S=106$⁠, the plasma relative speed can go above 10⁴. Our findings suggest that the parallel voltage is a more reliable indicator of the reconnection rate than the field-line speed.

In conclusion, we have tested Boozer's reconnection theory using the Boozer–Elder model for a coronal loop driven by footpoint motions. Our simulation results significantly differ from their predictions in both ideal and resistive evolution.

For the ideal evolution, we show that Boozer and Elder significantly under-predict the intensity of the electric current in the coronal loop due to missing terms in their equations. For the resistive evolution, our simulations show that the maximal current density roughly scales linearly with the Lundquist number S, in stark contrast to the prediction of a logarithmic dependence on S. Because of the formation of intense thin current sheets, the onset of fast reconnection occurs much sooner than predicted by Boozer and Elder. Therefore, our simulation results do not appear to support Boozer's theory. Moreover, thin current sheets become unstable to the plasmoid instability when the Lundquist number is sufficiently high.

A precise definition of 3D reconnection remains an open question. Boozer's definition of reconnection relies entirely on the connections between the fluid elements, and he attributes any changes in the connections to reconnection. This definition, while precise, is overly general and blurs the distinction between reconnection and diffusion. As an illuminating example, let us first consider a stable line-tied screw pinch^84,85 undergoing resistive diffusion. Because the footpoint mapping between the boundaries changes as time progresses, this process is reconnection according to Boozer's definition, although it occurs on a slow, resistive timescale. Next, let us consider the same process in a stable coronal loop with chaotic field lines. The time evolution of the magnetic field remains slow, but the field-line velocity is fast due to the amplification caused by the chaotic field-line mapping. This enhanced field-line speed is fast reconnection according to Boozer's definition. Examples of slow resistive diffusion with rapid changes in the field-line connections have been demonstrated in Ref. 63. However, the same study also shows that when the slow resistive diffusion leads to an unstable configuration, the time evolution becomes dynamic, intense thin current sheets form, and reconnection outflow jets ensue. At that time, field-line connections change even faster.

These results of Ref. 63 and this study indicate that the field-line velocity is not a reliable metric for the reconnection rate. We should keep in mind that the field-line velocity is a mathematical construct rather than a physical velocity. Furthermore, the field-line speed has no upper limit and can even exceed the speed of light! Therefore, we should be cautious about using the field-line velocity to draw conclusions. In comparison, the parallel voltage, which is the integration of the parallel electric field, appears to be a better indicator of the reconnection rate. Because the parallel voltage $Φ∼ηJL$ and the maximal current density $Jmax∝S∝1/η$⁠, the maximal voltage $Φmax$ is approximately independent of S. Therefore, the reconnection rate in the coronal loop is approximately independent of S. This conclusion is consistent with the dissipation rate being insensitive to S.

Our simulation results show that intense thin current sheets are a natural outcome of a coronal loop forced by slow footpoint motions. This conclusion appears to be more consistent with Parker's nanoflare scenario, where thin current sheets play a crucial role, than with Boozer's theory. However, Boozer's theory does provide some new perspectives and poses new challenges to Parker's scenario. In particular, Boozer points out the fragility of the ideal MHD frozen-in constraint in the presence of the chaotic field lines. This important aspect warrants broader attention and further investigation. Traditionally, the Parker problem is usually posed as a question regarding whether singular current sheets can form in a coronal loop under the frozen-in constraint and line-tied boundary conditions. This viewpoint attributes the formation of the current sheets solely as an ideal MHD effect. However, because the footpoint motion at the photosphere naturally renders the field line chaotic, the ideal MHD frozen-in constraint may be overly restrictive for the Parker problem. In this study, we have observed that resistive simulations can develop more intensive current sheets than an ideal simulation at the same time when both are still in a quasi-static phase, indicating that non-ideal effects may play an active role in the formation of the current sheets.⁶³ A similar suggestion has been made by Bhattacharjee and Wang, who proposed that helicity-conserving reconnection processes might facilitate the formation of the current sheets without rendering the footpoint velocity discontinuous.⁸⁶ The roles of non-ideal effects in current sheet formation and the onset of reconnection in coronal loops will be a topic of future study.

See the supplementary material for diagnostic results for runs A, B, C, and D3, similar to those shown in Fig. 7, which are available as animations.

This research was supported by the U.S. Department of Energy, Grant No. DE-SC0021205, and the National Aeronautics and Space Administration, Grant Nos. 80NSSC18K1285 and 88NSSC21K1326. Computations were performed on facilities at National Energy Research Scientific Computing Center. We thank Professor Allen Boozer for numerous beneficial discussions about his view on 3D reconnection, as well as the anonymous referees for insightful comments. This paper is dedicated to the memory of Aad van Ballegooijen, who made incisive contributions to the problem of current sheets in coronal loops and whose untimely passing in 2021 has deprived the community of a gentleman and a scholar.

The authors have no conflicts to disclose.

Yi-Min Huang: Conceptualization (equal); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Resources (equal); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Amitava Bhattacharjee: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Resources (equal); Writing – review & editing (equal).

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

After the initial submission of this paper, Boozer suggested that the observed signatures of reconnection, in particular, the intense thin current sheets, may be attributed to damping of Alfvén waves after changes of the field-line connections release the stored magnetic energy as the plasma kinetic energy. He further suggested that by increasing the viscosity to damp the kinetic energy, intense thin current sheets may disappear.⁸⁷ This hypothesis is interesting and has practical relevance as well, because observational evidence suggests that the magnetic Prandtl number $Pm≡ν/η$ of coronal loops may be orders of magnitudes larger than unity.⁸⁸ 

We test this hypothesis by performing additional simulations with $Pm≫1$⁠. We keep the viscosity at a constant value $ν=10−3$⁠, which is as high as possible without significantly compromising the force-free approximation of the magnetic field as the coronal loop evolves quasi-statically. The values of η vary from $10−4$ to $10−6$⁠.

Table II summarizes the parameters and diagnostic results of the high-P_m runs. Remarkably, resistive dissipation still dominates over viscous dissipation even when $Pm≫1$⁠. The resistive dissipation remains close to constant as η varies. Moreover, compared with the results in Table I, the resistive dissipation for $Pm≫1$ is similar to that for P_m = 1.

Figure 8 shows the probability distribution of the current density for various Lundquist numbers S. Similar to Fig. 2 for P_m = 1, here, the probability distribution also exhibits a strong dependence on the Lundquist number S, although the maximal current density is lower than the corresponding maximal current density with the same S and P_m = 1. A 2D slice of run H with $S=106$ and $Pm=103$ is shown in Fig. 9. Compared with Fig. 3, the plasma velocity is significantly smoother and slower due to the higher viscosity, whereas the current density is only slightly smoother.

Even though increasing the magnetic Prandtl number does smooth the current sheets, the maximal current density still exhibits a strong dependence on the Lundquist number S. Therefore, intense thin current sheets appear to be a natural consequence of this coronal loop model with high Lundquist numbers, regardless of the magnetic Prandtl number.