Coronal Heating Topology: the Interplay of Current Sheets and Magnetic Field Lines

TL;DR

This study explores the magnetic topology and field line behavior in a nanoflare-heated solar corona, revealing anisotropic diffusion along current sheets and implications for particle acceleration and heat transport.

Contribution

It introduces the concept of Current Sheet Connected (CSC) regions and analyzes their influence on magnetic field line diffusion in the corona.

Findings

CSC regions form coherent structures around current sheets.

Field line diffusion is strongly anisotropic along current sheets.

Enhanced magnetic field intensities increase displacements in CSC regions.

Abstract

The magnetic topology and field line random walk properties of a nanoflare-heated and magnetically confined corona are investigated in the reduced magnetohydrodynamic regime. Field lines originating from current sheets form coherent structures, called Current Sheet Connected (CSC) regions, extended around them. CSC field line random walk is strongly anisotropic, with preferential diffusion along the current sheets' in-plane length. CSC field line random walk properties remain similar to those of the entire ensemble but exhibit enhanced mean square displacements and separations due to the stronger magnetic field intensities in CSC regions. The implications for particle acceleration and heat transport in the solar corona and wind, and for solar moss formation are discussed.

Equations

$${\partial }_t\psi =\left[\phi ,\psi \right]+B_0{\partial }_z\phi +\eta {\nabla }^2\psi ,$$

$${\partial }_t\omega =\left[j,\psi \right]-\left[\omega ,\phi \right]+B_0{\partial }_{z}j+\nu {\nabla }^2\omega ,$$

$$\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}z}\left(z\right)=\frac{1}{B_0}\mathbf{b}\left(\mathbf{x}(z),z\right)\to \frac{\mathrm{d}}{\mathrm{d}z}{\mathbf{x}}^2=\frac{2}{B_0}\mathbf{x}\cdot \mathbf{b},$$

$$\frac{\mathrm{d}}{\mathrm{d}z}\left\langle\mathbf{x}^{2}(z)\right\rangle=\frac{2}{B_{0}^{2}}\int\limits_{0}^{z}\mathrm{d}z^{\prime}\ \big{\langle}\mathbf{b}\left(\mathbf{x}(z^{\prime}),z^{\prime}\right)\cdot\mathbf{b}\left(\mathbf{x}(z),z\right)\big{\rangle}.$$

$$C({\mathbf{x}}_1,z_1,{\mathbf{x}}_2,z_2)=⟨\mathbf{b}({\mathbf{x}}_1,z_1)\cdot \mathbf{b}({\mathbf{x}}_2,z_2)⟩.$$

$$C({\ell }_{\perp },|\Delta z|)=⟨\mathbf{b}(0,0)\cdot \mathbf{b}({\mathbf{\ell }}_{\perp },\Delta z)⟩,$$

$$\tilde{C}\left(z,\Delta z\right)=⟨\mathbf{b}\left(\mathbf{x}(z-\Delta z),z-\Delta z\right)\cdot \mathbf{b}\left(\mathbf{x}(z),z\right)⟩,$$

$$\frac{\mathrm{d}}{\mathrm{d}z}⟨{\mathbf{x}}^2(z)⟩=\frac{2}{B_0^2}\underset{0}{\overset{z}{\int }}\mathrm{d}\Delta z\tilde{C}\left(z,\Delta z\right).$$

$$\tilde{C}\left(\Delta z\right)\sim C(⟨{\mathbf{X}}^2(\Delta z){⟩}^{1/2},\Delta z).$$

$$E_{{\ell }_{\perp }}\propto {\ell }_{\perp }^{\alpha }⟶\delta b_{{\ell }_{\perp }}^2\propto {\ell }_{\perp }^{\alpha -1},$$

$$\frac{C({\ell }_{\perp },0)}{⟨b^2⟩}\sim \{\begin{array}{ll}1-{\left(\frac{{\ell }_{\perp }}{{\lambda }_{\perp }}\right)}^{\alpha -1}& \text{for}{\ell }_{\perp }\le {\lambda }_{\perp },\\ 0& \text{for}{\ell }_{\perp }\ge {\lambda }_{\perp },\end{array}$$

$$\frac{\tilde{C}(\Delta z)}{⟨b^2⟩}=\{\begin{array}{ll}1-{\left(\frac{\Delta z}{z_D}\right)}^{\alpha -1}& \text{for }\Delta z\le z_D,\\ 0& \text{for }\Delta z\ge z_D.\end{array}$$

$$D\sim \frac{\alpha -1}{2\alpha }{\lambda }_{\perp }\frac{⟨b^2{⟩}^{1/2}}{B_0},$$

$$\frac{\mathrm{d}\boldsymbol{\xi}}{\mathrm{d}z}(z)=\frac{1}{B_{0}}\Big{[}\mathbf{b}\left(\mathbf{x}_{2}(z),z\right)-\mathbf{b}\left(\mathbf{x}_{1}(z),z\right)\Big{]}.$$

$$\frac{\mathrm{d}⟨{\mathbf{\xi }}^2(z)⟩}{\mathrm{d}z}=\frac{2}{B_0^2}\underset{0}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }$$

$$\cdot$$

$$\frac{\mathrm{d}}{\mathrm{d}z}⟨{\mathbf{\xi }}^2(z)⟩=\frac{4}{B_0^2}\left(I_{11}-I_{12}\right),\mathrm{where}:$$

$$\displaystyle I_{11}=\int\limits_{0}^{z}\mathrm{d}z^{\prime}\ \big{\langle}\mathbf{b}\left(\mathbf{x}_{1}(z^{\prime}),z^{\prime}\right)\cdot\mathbf{b}\left(\mathbf{x}_{1}(z),z\right)\big{\rangle},$$

$$\displaystyle I_{12}=\int\limits_{0}^{z}\mathrm{d}z^{\prime}\ \big{\langle}\mathbf{b}\left(\mathbf{x}_{1}(z^{\prime}),z^{\prime}\right)\cdot\mathbf{b}\left(\mathbf{x}_{2}(z),z\right)\big{\rangle}.$$

$$\frac{\mathrm{d}}{\mathrm{d}z}⟨{\mathbf{\xi }}^2(z)⟩=2z\frac{⟨\delta b_{{\xi }_0}^2⟩}{B_0^2},\mathrm{hence}$$

$$⟨{\mathbf{\xi }}^2(z)⟩-{\mathbf{\xi }}_0^2=z^2\frac{⟨\delta b_{{\xi }_0}^2⟩}{B_0^2},$$

$$⟨{\mathbf{\xi }}^2(z)⟩-{\mathbf{\xi }}_0^2\approx \{\begin{array}{ll}\frac{2}{{\lambda }_{\perp }^{\alpha -1}}\frac{⟨b^2⟩}{B_0^2}{z}^2{\xi }_0^{\alpha -1},& \mathrm{for}{\xi }_0\le {\lambda }_{\perp }\\ 2\frac{⟨b^2⟩}{B_0^2}{z}^2,& \mathrm{for}{\xi }_0\ge {\lambda }_{\perp }\end{array}$$

Full text

Sun: activity — Sun: corona — Sun: magnetic fields — solar wind — magnetohydrodynamics (MHD) — turbulence

1 Introduction

The stochastic properties of magnetic fluctuations in turbulent plasmas are reflected in the stochastic character of magnetic field lines, giving rise to field line random walk (FLRW, Jokipii & Parker, 1968; Jokipii, 1973; Matthaeus et al., 1995) that strongly affects the propagation and cross-field transport of energetic particles. Additionally the intense electric fields associated to turbulent coherent structures, such as current sheets (and the related in- and out-flows), strongly contribute to particle acceleration (Swann, 1933; Drake et al., 2005; Dalena et al., 2014).

Because in current sheets particles are energized and plasma heated, the topology of field lines that originate from them determines how these accelerated particles and heat are transported. Furthermore in strong magnetic fields, where particle diffusion perpendicular to field lines is small and thermal conduction highly anisotropic (essentially parallel), heat and particles are in first approximation transported along field lines. It is therefore key to understand the interplay between current sheets and magnetic field lines.

It has become increasingly clear that the effective heating and particles acceleration occur at scales of the order of the ion (proton) inertial length $d_{i}$ (Sonnerup, 1979; Mozer et al., 2002; Runov et al., 2006; Phan et al., 2007; Xiao et al., 2007; Wu et al., 2011; Liu et al., 2014), that in the solar corona, for an ion density $n_{i}\sim 10^{8}$ cm*-3*, is $d_{i}=c/\omega_{pi}\sim 23$ m ($\omega_{pi}=\sqrt{4\pi n_{i}e^{2}/m_{i}}$ is the proton plasma frequency, $c$ the speed of light, $e$ the electron charge, and $m_{i}$ the proton mass). For typical hot coronal loops with temperatures $T\sim 10^{6}$ K and magnetic field intensities $B\sim 50$ G the ion gyroradius is much smaller than $d_{i}$ (reaching $d_{i}$ only in the higher $\beta$ regions typical of the solar wind).

In situ measurements in Earth’s magnetotail (Runov et al., 2005, 2006) and magnetosheath (Retinò et al., 2007), and laboratory experiments (Matthaeus et al., 2005; Yamada et al., 2006) show that current sheet thickness is generally somewhat larger than the ion inertial length, with activity increasing for thinner current sheets as their width approaches $d_{i}$.

Additionally PIC simulations of fully developed turbulence have shown that clustering of current sheet thickness occurs at scales of $\sim d_{i}$, with substructures down to the electron inertial length $d_{e}$ (Karimabadi et al., 2013). These are the natural scales at which kinetic effects will convert the energy coming from large scales into the different species thermal and non-thermal energies.

Note that these findings are also consistent with the recent understanding that thin current sheets are strongly unstable under the plasmoid instability (Bulanov et al., 1978; Biskamp, 1986; Loureiro et al., 2007; Lapenta, 2008; Bhattacharjee et al., 2009), with growth rates reaching fast “ideal” Alfvén values ($\gamma\tau_{A}\sim 1$) for sufficiently small thicknesses (Pucci & Velli, 2014; Tenerani et al., 2015; Landi et al., 2015). Including the Hall effect the instability becomes explosive as the current sheet thickness approaches $d_{i}$ (Pucci et al., 2017).

Although the aforementioned studies include at most a weak guide magnetic field, the formation of current sheets with the exponentially thinning widths have been observed in fully nonlinear 2D and 3D MHD simulations (Sulem et al., 1983; Frisch et al., 2003; Krstulovic et al., 2011; Brachet et al., 2013), and line-tied simulations with a strong guide field and vanishing initial velocity (Rappazzo & Parker, 2013). Although kinetic simulations with a strong guide field are still computationally challenging, we expect that the overall phenomenology and current sheet structure is not substantially modified in the strong guide field case of interest to the solar corona and inner heliosphere. We then consider the gyroradii of bulk ions and electrons to be generally smaller than the current sheet thickness, and the initial stage of their acceleration is thus strongly affected by the field line topology.

FLRW in turbulent fields is a topic of intense research (Jokipii & Kota, 1989; Zimbardo et al., 2000; Snodin et al., 2013; Beresnyak, 2013), but little attention has been dedicated to the effects of spatial intermittency and coherent structures (Pucci et al., 2016). However, the plasma thermodynamical properties are strongly affected by the topology of field lines originating in current sheets, both in the corona and solar wind. For instance the thermodynamics and high-energy radiative emission of coronal loops are determined by the temporal and spatial properties of energy dissipation along the field lines (Klimchuk, 2006; Reale, 2014). Additionally, energetic particles and heat transport toward the transition region at coronal loops footpoints give rise to a reticulated spongy pattern in X-rays and EUV, so-called moss, that could be explained by the complex trajectories of energetic particles in a stochastic magnetic field (Kittinaradorn et al., 2009).

Here we investigate the magnetic field lines random walk in a nanoflare-heated and magnetically confined corona to advance our understanding of transport of heat and energetic particles in stochastic magnetic fields, its relationship with coherent structures, and discuss its impact on coronal and solar wind dynamics.

2 Model

Our model coronal loops are “straightened-out” in a Cartesian elongated box with axial length L (along the z-direction) and orthogonal square cross section of size $\ell$ (x-y planes), with aspect ratio $L/\ell$=10. The system, with uniform density $\rho_{0}$, is threaded by a strong axial magnetic field $\mathbf{B}_{0}=B_{0}\mathbf{\hat{e}}_{z}$, and its dynamics are well described by the reduced magnetohydrodynamic (RMHD) equations (Kadomtsev & Pogutse, 1974; Strauss, 1976; Montgomery, 1982; Zank & Matthaeus, 1992), valid in the limit of a large loop aspect ratio ($\epsilon=\ell/L\ll 1$) and of a small ratio of orthogonal to axial magnetic field ($b/B_{0}\leq\epsilon$). The velocity ($\mathbf{u}$) and fluctuating magnetic field ($\mathbf{b}$) have only components perpendicular to the axial direction $z$, and indicating their potentials with $\varphi$ and $\psi$, they can be written as $\mathbf{u}=\nabla\varphi\times\mathbf{\hat{e}}_{z}$ and $\mathbf{b}=\nabla\psi\times\mathbf{\hat{e}}_{z}$, with the current density $j=-\nabla^{2}\psi$, and vorticity $\omega=-\nabla^{2}\varphi$. In non-dimensional form they are given by:

[TABLE]

where the magnetic field has been expressed as an Alfvén velocity (i.e., $b\rightarrow b/\sqrt{4\pi\rho_{0}}$), and then all velocities normalized to $u^{\ast}$ = 1 km s*-1* (the photospheric granulation velocity rms). The Poisson bracket is defined as, e.g., $[\varphi,\psi]=\partial_{x}\varphi\partial_{y}\psi-\partial_{y}\varphi\partial_{x}\psi=-\mathbf{u}\cdot\nabla\psi$, and the Laplacian operator has only orthogonal components. Lengths are normalized to the orthogonal box length, thus $\ell$=1 and $L$=10. Normalized resistivity and viscosity coefficients are set equal with $\eta$=$\nu$=1/$R$, where the Reynolds number $R$=$800$, numerical resolutions is 2048${}^{2}\times$512, and the guide field intensity $B_{0}$=$10^{3}$. As in previous simulations field lines are line-tied to the top and bottom plates $z=$0,10 where they are shuffled by a prescribed photospheric granulation-mimicking velocity constant in time with length-scale $\sim$1/4 (Rappazzo et al., 2008), while in x-y planes periodic boundary conditions are implemented. In the x-y planes a pseudo-spectral scheme with periodic boundary conditions is implemented, time is advanced with a third-order Runge-Kutta and an adaptive time-step. More details on the model and numerical code can be found in Rappazzo et al. (2007, 2008).

The numerical integration of Eqs. (1)-(2) cannot implement enough grid points to attain a realistic description of the internal structure of current sheets (that would additionally require the inclusion of the Hall term, or the integration of a kinetic model, that in turn could not describe properly the large scale dynamics). Nevertheless, while a good representation of the small-scale structure is very important for the acceleration of particles (and we have developed a hierarchical multi-scale model to this effect in Dalena et al., 2014), since FLRW is affected mostly by the large scale components of the magnetic field, its properties should not depend critically on the small-scale structure of the current sheets, and an MHD model represents a good starting point. For these reasons we adopt an empirical approach, selecting the value of resistivity for essentially numerical reasons, but subsequently associating the resulting dissipative scale with the ion inertial length.

3 Results

Our simulations start with the guide field $B_{0}$ (directed along z) and no magnetic or velocity fluctuations in the computational box. The imposed large-scale velocity at the boundaries z=0, and L twists the field lines and, once the twist exceeds a small critical threshold, the orthogonal magnetic field line tension is no longer balanced. Thus, as proposed by Parker (1972, 1994), the magnetic field $\mathbf{b}$ transitions to non-equilibrium (Rappazzo & Parker, 2013; Rappazzo, 2015), bringing about turbulent dynamics that transfers energy towards the small scales where it is dissipated in nanoflares (Parker, 1988). Line-tying keeps the velocity field in the computational box smaller than the magnetic field (far from equipartition). Nevertheless a cascade with preferential energy transfer in the x-y planes orthogonal to $B_{0}$ and a broad-band power-law magnetic energy spectrum develop, but the enhanced field lines stiffness introduced by line-tying gives rise to steeper magnetic energy spectra $E_{M}(k_{\perp})\propto k_{\perp}^{-\alpha}$ with $\alpha\in[5/3,3]$, with the steepest spectra corresponding to stronger guide fields $B_{0}$ (Einaudi et al., 1996; Dmitruk & Gómez, 1997; Dmitruk et al., 1998; Dmitruk & Gómez, 1999; Dmitruk et al., 2003; Rappazzo et al., 2007, 2008; Rappazzo & Velli, 2011). In the simulations considered here $\alpha\sim$ 3.

The magnetic field structure is characterized by approximately field-aligned current sheets. Although the overall physical conditions are markedly different between the line-tied and fully periodic reduced MHD, we hypothesize that the FLRW properties of field lines traced from current sheets are qualitatively the same, since in both cases current sheets are field aligned (e.g., Dmitruk et al., 2005) and their axial extension must be linked to the parallel correlation length. Clearly the current sheet extension in the axial direction can change depending on the particular type of forcing, but our conclusions on field line random walk and diffusion properties may be tentatively extended to the fully periodic case, as discussed in our Conclusions, section 4.

3.1 Field line diffusion from a single point

Since in reduced MHD the $z$-component ($B_{0}$) of the magnetic field is constant, the magnetic field line equation can be written as

[TABLE]

where $\mathbf{x}=(x,y)$ indicates the orthogonal coordinates. The mean square displacement $\langle\mathbf{X}^{2}(z)\rangle=\langle[\mathbf{x}(z)-\mathbf{x}_{0}]^{2}\rangle=\langle\mathbf{x}^{2}(z)\rangle-\mathbf{x}_{0}^{2}$, where $\langle\ldots\rangle$ indicates ensemble average, is shown in Fig. 1 (top panel). After an initial ballistic stage with $\langle\mathbf{X}^{2}\rangle\propto z^{2}$, it subsequently exhibits diffusion with $\langle\mathbf{X}^{2}\rangle=4Dz$. To understand this behavior, from Eq. (3) we can write

[TABLE]

Although the position vectors $\mathbf{x}(z)$ are random functions determined by the trajectory, the integrand in Eq. (4) is linked to the magnetic field two-point correlation function

[TABLE]

For homogeneity and isotropy this depends only on the relative parallel and orthogonal distances of the two points, i.e., indicating with $\boldsymbol{\ell}_{\perp}=\mathbf{x}_{2}-\mathbf{x}_{1}$, and $\Delta z=z_{2}-z_{1}$:

[TABLE]

independent of the origination point (as long as both points are within the z-span). As shown in Fig. 1 (center), the correlation decreases at larger $\ell_{\perp}$ and $z$. But while it vanishes in the perpendicular direction at the correlation length $\lambda_{\perp}\sim 0.11$, it does not vanish in the parallel direction (for $\ell_{\perp}=0$). Namely the parallel correlation length $\lambda_{\parallel}$ is larger than the box size $L$, i.e., the turbulent field has a strong 2D component. Clearly this is due to the low frequency of photospheric motions. Indeed, for typical hot loops, the field line footpoints are shuffled slowly compared to the fast Alfvén crossing timescale at which the induced magnetic field twist propagates along the loop axis.

The correlation in Eq. (4) is Lagrangian, i.e., it is computed along the field lines:

[TABLE]

with $0\leq\Delta z\leq z$, and it is shown in Fig. 1 (bottom panel). Introducing the change of variable $\Delta z=z-z^{\prime}$ we can then write

[TABLE]

Since the mean square displacement between two points along a field line at a parallel distance $\Delta z$ is to a good approximation given by $\langle\mathbf{X}^{2}(\Delta z)\rangle$, the two correlations are then approximately linked by

[TABLE]

Additionally the correlation function is connected to the second-order structure function as $C(0,0)-C(\ell_{\perp},0)=\langle\delta b_{\ell_{\perp}}^{2}\rangle/2$ (e.g., see Biskamp, 2003), that in turn, for values of $\ell_{\perp}$ in the inertial range is linked to the magnetic energy spectrum by $E_{\ell_{\perp}}\propto\ell_{\perp}\delta b_{\ell_{\perp}}^{2}$, consequently

[TABLE]

with $\alpha\in[5/3,3]$ for our boundary forced coronal loop model (Dmitruk et al., 2003; Rappazzo et al., 2007, 2008). Therefore extending the power-law behavior beyond the inertial range for all $\ell_{\perp}\leq\lambda_{\perp}$, and taking into account that the correlation vanishes at $\lambda_{\perp}$, we can approximate the magnetic correlation function with

[TABLE]

with the exponent ranging from 2/3 for $\alpha=5/3$ up to 2 for $\alpha=3$. This function is plotted in Fig. 1 (middle) for $\alpha=3$. Here we neglect the parallel variation of $C$ when used in Eq. (9) because $\langle\mathbf{X}^{2}(\Delta z)\rangle$ increases monotonically with $\Delta z$, and as shown in Fig. 1 (middle) the curves then tend to overlap quickly becoming approximately independent of $\Delta z$.

The behavior of the mean square displacement $\langle\mathbf{X}^{2}(z)\rangle$ can then be readily understood from the correlation function $C(\ell_{\perp},\Delta z)$ (Fig. 1). For small values of $z$ the integral in Eq. (8) can be Taylor-expanded, and since $C(0,0)=\langle b^{2}\rangle$, to the first order we obtain $\langle\mathbf{X}^{2}(z)\rangle\approx z^{2}\langle b^{2}\rangle/B_{0}^{2}$. On the other hand as soon as $\langle\mathbf{X}^{2}(z)\rangle^{1/2}\gtrsim\lambda_{\perp}$ exceeds the orthogonal correlation length $\lambda_{\perp}\sim 0.11$, the integral in Eq. (8) remains approximately constant, because the largest contribution comes from $\ell_{\perp}<\lambda_{\perp}$, hence $\langle\mathbf{X}^{2}(z)\rangle\sim 4Dz$ diffuses linearly. The transition from the ballistic to the diffusive stage occurs for $\langle\mathbf{X}^{2}(z)\rangle\sim z^{2}\langle b^{2}\rangle/B_{0}^{2}\sim\lambda_{\perp}^{2}$, i.e., for $z_{D}\sim\lambda_{\perp}B_{0}/\langle b^{2}\rangle^{1/2}$. In our case, since $\lambda_{\perp}\sim 0.11$, $B_{0}=10^{3}$, and $\langle b^{2}\rangle^{1/2}\sim 20$, the transition occurs at $z_{D}\approx 5.5$, as confirmed in Fig. 1 (top).

We can estimate the diffusion coefficient $D$ by using in Eq. (9) the approximation for the mean square displacement outlined in the previous paragraph (i.e., $\langle\mathbf{X}^{2}(\Delta z)\rangle\approx\Delta z^{2}\langle b^{2}\rangle/B_{0}^{2}$ for $z\leq z_{D}$). The Lagrangian correlation along the field lines can then be approximated from Eqs. (9) and (11) with

[TABLE]

Substituting in Eq. (8), and integrating it to obtain $\langle\mathbf{X}^{2}(z)\rangle$, the diffusion coefficient $D$ is then given by

[TABLE]

a functional form characteristic of Bohm diffusion (Ghilea et al., 2011), with the coefficient $(\alpha-1)/2\alpha$ ranging in the narrow interval $[1/5,1/3]$ as $\alpha\in[5/3,3]$. Since $\langle b^{2}\rangle^{1/2}\sim 20$, $B_{0}=10^{3}$, $\lambda_{\perp}\sim 0.11$ and $\alpha\sim 3$, we obtain $D\sim 7.3\times 10^{-4}$, corresponding to $D\sim 2.9$ km in conventional units, compatible with $D\sim 3.1$ km computed from our simulation (Fig. 1, top panel). The approximated Lagrangian correlation function $\widetilde{C}(\Delta z)$ (Eq. (12)) with $\alpha=3$ and $z_{D}\sim 5.5$ is shown in Fig. 1 (bottom panel).

3.2 Current sheet connected regions

The current density $j$ is intermittently distributed in space, as typical of turbulent systems. Its Probability density function (pdf), shown in Fig. 2, is not gaussian and exhibits typical large tails where the current is strong, corresponding to current sheets in physical space. The noticeable skewness in Fig. 2 results from the use of a single snapshot, and it is a fluctuation that vanishes when averaging over several snapshots, i.e., the time-averaged distribution is symmetric.

We define as current sheets all those spatial regions where current is larger than two standard deviations $|j|\geq 2\sigma$, with $\sigma=\langle j^{2}\rangle^{1/2}$, shown in yellow for the representative planes z=0, L/2 and L in Fig. 3.

The relationship between magnetic field topology and current sheets is then investigated by tracing field lines originating in current sheets. Specifically field lines are traced from all the grid points where $|j|\geq 2\sigma$ in 9 equispaced x-y planes (from z=0 up to z=L, separated by L/8). For grid points that are not at the boundaries z=0 or z=L the respective field lines are traced both forward and backward with respect to the z-direction. We trace a total of 1,107,242 field lines, extending from the bottom to the top plate, and their intersection with the selected plane is shown as a black dot in Fig. 3. Clearly field lines are present also in the yellow regions, both those traced from there plus others originating from current sheets in different planes.

Current sheets are elongated in the guide field direction $z$, and the field lines traced from them form similarly shaped coherent structures, that we indicate as Current Sheet Connected (CSC) regions. Although current sheets in reduced MHD have a complex structure with a cross-shear magnetic field component and mostly external X-points (Zhdankin et al., 2013; Wan et al., 2014), noticeably the presence of a strong magnetic shear in correspondence of current sheets makes the field line random walk strongly anisotropic, with field lines diffusing preferably along the in-plane sheet length and very little across it (Fig. 3). Since diffusion increases with distance its effects are most apparent in planes z=0 and L.

The coherence and strong anisotropy of the CSC regions are in stark contrast with the homogeneity of the stochastic properties typically associated with FLRW, and their well-know tendency to fragment flux tubes (Matthaeus et al., 1995; Servidio et al., 2014; Rappazzo et al., 2012). Nevertheless the mean square displacement $\langle\mathbf{X}^{2}(z)\rangle$ of the CSC field lines has properties similar to those of the entire ensemble, as shown in Fig. 1 (top). The higher value of the diffusion coefficient $D_{\rm CSC}/D_{\rm All}\sim 1.2$ is due to the higher magnetic field intensity in the CSC region as $\langle b^{2}\rangle_{CSC}^{1/2}/\langle b^{2}\rangle_{All}^{1/2}\sim 1.2$ in agreement with Eq. (13).

3.3 Pair separation

To further understand the magnetic topology we consider the separation of field line pairs. From Eq. (3) their orthogonal separation in the $x$–$y$ plane $\boldsymbol{\xi}(z)=\mathbf{x}_{2}(z)-\mathbf{x}_{1}(z)$, given the initial separation $\boldsymbol{\xi}(0)=\boldsymbol{\xi}_{0}$, is determined as a function of $z$ by

[TABLE]

Similarly to the single field line case we obtain

[TABLE]

that following Ruffolo et al. (2004) can be written as

[TABLE]

$I_{11}$ is same as the integral in Eq. (4) because it refers to single field lines, therefore we already understand its behavior. But $I_{12}$ differs as it considers a pair.

As indicated by Batchelor (1950) for the hydrodynamic case, for $z$ sufficiently small the mean-square separation will grow quadratically with z (see inset in Fig. 4, left panel) and it will be proportional to the second-order structure function. Indeed from Eq. (16), Taylor-expanding Eqs. (17)-(18) in z, and since $C(0,0)-C(\mathbf{\xi}_{0},0)=\langle\delta b_{\xi_{0}}^{2}\rangle/2$ we obtain

[TABLE]

where as usual $\xi_{0}=|\boldsymbol{\xi}_{0}|$. Additionally we can now approximate the second order structure function as in Eq. (11) with $\langle\delta b^{2}_{\xi_{0}}\rangle\sim 2\langle b^{2}\rangle(\xi_{0}/\lambda_{\perp})^{\alpha-1}$ for $\xi_{0}\leq\lambda_{\perp}$ and $\langle\delta b^{2}_{\xi_{0}}\rangle\sim 2\langle b^{2}\rangle$ for $\xi_{0}\geq\lambda_{\perp}$, with the scaling relation more accurate for values of $\xi_{0}$ in the inertial range. We can then write

[TABLE]

with $\alpha\in[5/3,3]$.

We plot $(\langle\boldsymbol{\xi}^{2}(z)\rangle-\boldsymbol{\xi}^{2}_{0})/z^{2}$ in Fig. 4 (right panel) as a function of $\xi_{0}$ for 11 values of $z\in[0,2]$ separated by 0.2, almost perfectly overlapping and showing that in the inertial range ($\xi_{0}\lesssim 0.1$) it scales approximately as $\xi_{0}^{2}$, compatible with $\alpha\sim 3$ in our simulations, and saturates correctly to $\sim 2\langle b^{2}\rangle/B_{0}^{2}\sim 8\times 10^{-4}$ with our parameters ($\langle b^{2}\rangle\sim 400$, $B_{0}=10^{3}$, $\lambda_{\perp}\sim 0.11$). The small departure from $\xi_{0}^{2}$ in the inertial range shown in Fig 4 occurs because the second order structure function $\langle\delta b_{\xi_{0}}^{2}/2\rangle=C(0,0)-C(\xi_{0},0)$ considers zero separation in z and is calculated in the bottom boundary plane $z=0$ (since the coefficients in the Taylor expansion in z of Eqs. (17)-(18) are calculated for z=0). Indeed the 2D magnetic energy spectrum averaged over the whole box $E_{\ell_{\perp}}\propto\ell_{\perp}\delta b_{\ell_{\perp}}\sim\ell_{\perp}^{3}$, corresponding to $\delta b_{\ell_{\perp}}\sim\ell_{\perp}^{2}$, and the same behavior is observed for the spectra in all z-planes, except those in proximity of the boundaries z=0 and L where line-tying boundary conditions are applied. At these boundaries the velocity field is prescribed and therefore the dynamics does not follow the same equations as in the interior, with the effect of slightly modifying the magnetic energy spectrum for the planes in their close proximity. Nevertheless the scaling for the separation remains very close to a $\propto\xi_{0}^{2}$ scaling, departing strongly in the inertial range from $\xi_{0}^{2/3}$ expected for a standard Kolmogorov spectrum with $\alpha=5/3$ (Batchelor, 1950; Bourgoin et al., 2006). Therefore the ratio of mean square separations for field lines with relative larger initial separations $\xi^{\prime}_{0}>\xi_{0}$ are increasingly bigger for steeper spectral indices, indeed from Eq. (21) the ratio of their separations in the ballistic range ($\propto z^{2}$) grows like $(\xi^{\prime}_{0}/\xi_{0})^{\alpha-1}$.

In general, as $\xi_{0}\rightarrow 0$ the two field lines tend to the same field line, i.e., $\mathbf{x}_{2}(z)\rightarrow\mathbf{x}_{1}(z)$, consequently $I_{12}\rightarrow I_{11}$, and the mean square separation vanishes in this limit. On the other hand the mean square separation in $I_{12}$, i.e. $\xi_{12}(z^{\prime},z)=\langle[\mathbf{x}_{2}(z^{\prime})-\mathbf{x}_{1}(z)]^{2}\rangle^{1/2}$ is always larger than $\xi_{0}$. For initial separations larger than $\lambda_{\perp}$ the correlation is small $\forall z\in[0,L]$ so that $I_{12}\approx 0$. In this case the diffusion coefficient for pair separation is double that of single field line diffusion, i.e., $\langle\boldsymbol{\xi}^{2}(z)\rangle-\boldsymbol{\xi}^{2}_{0}=4D_{\mathrm{p}}z$ with $D_{\mathrm{p}}=2D$, as shown in Fig. 4 (left panel) for $0.1<\xi_{0}<1$ by the continuous black line (we average these curves since they overlap). For any $\xi_{0}<\lambda_{\perp}$ there is always a critical height $z_{D}(\xi_{0})$ above which the mean separation between the two field lines in $I_{12}$ is larger than $\lambda_{\perp}$. Hence the increasingly larger negative contribution of the $I_{12}$ term to mean square pair separation will display diffusion at progressively lower heights, with a smaller total diffusion coefficients $D\leq D_{p}\leq 2D$, as shown in Fig. 4 by the red lines, that consider 5 bins with $\xi_{0}\in[0.015,0.095]$ with $\Delta\xi_{0}=0.02$. Clearly, for sufficiently small initial separations, field lines will not be able to display diffusion because our system is bounded in the axial direction z to maximum length L and separation cannot grow up to the perpendicular correlation length, as shown well by the blue continuous lines in Fig. 4.

A ballistic stage $\propto z^{2}$ is always present initially for $z\lesssim 3$, as shown in Fig. 4 (insets), but field lines with initial separation within the dissipative range ($\xi_{0}\lesssim 10$ dx, where dx=1/nx=1/2048 is the numerical grid step size) exhibit subsequently a Richardson-like superdiffusive stage (Beresnyak, 2013; Eyink et al., 2013; Servidio et al., 2016) with mean square separation up to $\propto z^{4}$ (blue continuous lines in Fig. 4, with respectively $\xi_{0}=2.5,$ 6.5, and 10.5 dx), while for larger separations they transition to the diffusive regime.

Even though FLRW in CSC regions is strongly anisotropic, CSC field line pair separation (shown in Fig. 4 with dashed lines) exhibit similar properties to the whole ensemble (shown with continuous lines, color code is the same for both line types). Their statistics are degraded for larger $\xi_{0}$ because in CSC regions the number of field line couples diminishes at larger separations so that the averages do not saturate yet to their ensemble value. This can be seen in the left inset, where the red dashed curves become negative, because the relation $\langle\boldsymbol{\xi}^{2}(z)\rangle-\boldsymbol{\xi}_{0}^{2}=\langle(\boldsymbol{\xi}(z)-\boldsymbol{\xi}_{0})^{2}\rangle\geq 0$ is valid only for a sufficiently high number of field lines, when $\langle\boldsymbol{\xi}(z)\rangle=\boldsymbol{\xi}_{0}$.

The main difference between CSC pair separations and the whole ensemble is that for same initial separation $\xi_{0}$ the CSC field lines exhibit higher separations. Similarly to mean square displacement (Fig. 1, top panel) the higher values are due to the greater magnetic field intensity in CSC regions ($\langle b^{2}\rangle_{CSC}/\langle b^{2}\rangle\sim 1.5$), and indeed the separation in the ballistic range is proportional to $\langle b^{2}\rangle$ (Eq. (21)). Notice that the coefficient in Eq. (21) includes also the perpendicular correlation length $\lambda_{\perp}$ that is not readily computable in the non-Cartesian CSC region, but from the data we can estimate that the separation for field lines with same initial separation is about 10 times larger in CSC regions than for the whole ensemble. As mentioned previously the statistics are degraded for CSC field lines with larger initial separations. Therefore while the larger separations for CSC field lines is very well demonstrated for smaller initial separations (blue dashed lines, and first red dashed line), we cannot yet fully verify this conclusion for larger initial separations, as the averages have not yet saturated to their ensemble values (this point will be further investigated in upcoming work).

4 Conclusions and Discussion

To gain insight into particle acceleration and heat transport in coronal loops, solar wind, and more in general for plasmas in the reduced MHD regime, we have investigated the magnetic topology of field lines originating from current sheets. We have found that they form coherent structures, dubbed Current Sheet Connected (CSC) regions, similarly to the current sheets they originate from. Field lines in these regions perform highly anisotropic FLRW, with diffusion occurring preferentially along the current sheet in-plane length. Nevertheless FLRW and diffusion coefficients have similar properties for CSC field lines and the whole ensemble, with larger displacement and separations occurring in CSC regions where the magnetic field intensity is higher.

This emerging picture has strong implications for particle acceleration and heat transport, particularly in the low corona, where all protons and electrons with temperatures below $10^{6}$ K have gyroradii smaller than the current sheet thickness. It indeed implies that in coronal loops particle and heat are transported almost exclusively within the CSC region, a small volume of plasma around current sheets with a small filling factor, while most of the volume is topologically disconnected from current sheets and the associated flow of particles and heat.

This picture is fully consistent with observations (Schmelz et al., 2010; Warren et al., 2008; Ugarte-Urra et al., 2009; Klimchuk, 2009) and recent thermodynamical 3D simulations (Dahlburg et al., 2012, 2016) strongly suggesting that coronal loops cannot be modeled with single isothermal flux tubes, as their radiative properties can only be explained by the presence of both hot and cold plasmas at observational sub-resolution scales (multi-temperature loops). Also, the complex topology in CSC regions, with enhanced magnetic field line displacements and separations, points to a complex stochastic nature for the heating function along the field lines.

Additionally the structure of the CSC regions in the top and bottom plates z=0 and L is consistent with that of so-called moss, the spongy reticulated pattern in X-rays and EUV formed at the coronal base of hot loops, confirming that the FLRW can play a strong role in the formation of these structures as recently proposed (Kittinaradorn et al., 2009).

The reduced MHD FLRW properties and topology strongly support the results of recent test-particle simulations (Dalena et al., 2014) with initial gyroradii smaller than current sheet widths propagated in similar magnetic fields to those discussed here. Particles are at first strongly accelerated along the z-direction in current sheets by the strong electric field associated with the current, until they pitch-angle scatter thus increasing their gyroradii. Subsequently, as long as the gyroradius is smaller than the orthogonal correlation length, those that remain close to the CSC region are accelerated by a (non-magnetic moment conserving) betatron-like mechanism due to the inhomogeneous $\mathbf{u}\times B_{0}\mathbf{\hat{z}}$ electric field associated to outflows in current sheets.

For parameters typical of hot solar coronal loops we have found that the parallel correlation length is longer than the loop length, and our results consider this specific case. Nevertheless they may also apply to the unbounded (e.g., periodic) case, which we plan to investigate thoroughly in upcoming work. In general we expect the current sheet length along z to be strongly correlated with the magnetic field parallel correlation length $\lambda_{\parallel}$. Therefore we expect that CSC regions connected to any such current sheet of length $\sim\lambda_{\parallel}$ to have a similar structure to those found here around the current sheet. But as field lines are traced further away at distances larger than $\lambda_{\parallel}$ we expect the CSC region to fragment and the associated FLRW to lose anisotropy and acquire more homogeneous properties, i.e., the CSC field lines will at that point connect and diffuse isotropically throughout the volume (hence mostly in regions with low current).

Therefore in unbounded systems we expect heat and accelerated particles to be initially confined to CSC regions around current sheets (of length $\sim\lambda_{\parallel}$), but further away heat and particles would distribute more uniformly throughout the plasma. This picture is strongly consistent with recent analyses of solar wind data (Osman et al., 2012), where temperature is found to peak in regions with high magnetic field gradients, while it rapidly descends to approach the ambient solar wind temperature as distance from those regions increases.

Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This work was partially supported by NASA through Grand Challenge Research grant NNX15AB88G, Living with a Star grant NNX14AI63G, as well as the Solar Probe Plus Observatory Scientist grant, and by the Thailand Research Fund (grant RTA5980003).