Energetics and 3D Structure of Elementary Events in Solar Coronal Heating

Loops in the solar corona are characterized by the steady and diffuse emission of X-rays and extreme-ultraviolet radiation. The coronal temperatures implied by this radiation indicate that in some way the corona is heated to temperatures significantly greater than the underlying chromosphere. Many different scenarios are at present under investigation to explain this heating. In all cases the ultimate energy source derives from convective motions in which there is more than enough energy to supply total coronal losses, the main question being how the energy is deposited in the corona. In the past few decades it has become more and more evident that the heating process involves the continuous formation and disruption of many small-scale current sheets. In each current sheet the energy is dissipated and transformed into internal energy with a consequent formation of temperature peaks. This local phenomenon can be termed the "elementary event" underlying coronal heating (Chiuderi 1993). In this paper we present numerical evidence relating to the extent, duration, and energetics of these elementary events.

In this introductory section we detail the history of the "elementary event" concept, beginning with the work of Parker (1972), who was the first to explore in detail the Gold & Hoyle (1960) idea that coronal heating could be the necessary outcome of an energy flux associated with the tangling of coronal field lines by photospheric flows. Parker (1983, 1988), Sturrock & Uchida (1981), Van Ballegooijen (1986), Gomez & Ferro Fontan (1988), and Berger (1991), among others, further explored the dynamics caused by such random shuffling of magnetic field lines, using weakly nonlinear analyses, turbulence phenomenologies, and/or extended dimensional arguments. The 3D simulations performed by Mikic et al. (1989), Strauss (1993), Longcope & Sudan (1994), Galsgaard & Nordlund (1996), and Hendrix & Van Hoven (1996) confirmed that this process causes the formation and exponential growth of local coronal currents.

To perform moderate- to high-resolution simulations of the system, allowing for much larger integration times, Einaudi et al. (1996) developed a 2D model, focusing on the dynamics in a given plane. The photospheric forcing function used to feed energy into the plasma has large-scale spatial structures (with amplitudes changing at the photospheric characteristic convection timescale), so that small scales are not directly excited. Therefore, in the 2D simulations, a magnetic forcing is imposed varying on the typical photospheric timescale τ_p , a first attempt at studying the Parker scenario with significant resolution, in an effort to detail the perpendicular dynamics without any assumptions on the properties of the resulting current sheets (Einaudi et al. 1996; Dmitruk et al. 1998; Georgoulis et al. 1998). The transverse dynamics for times longer than τ_p was shown to alter the magnetic field structure substantially from the spatial structure of the magnetic forcing, and the system evolution is turbulent and highly dynamical, exhibiting intermittency for both the spatial current distribution and the mean dissipation time series. Such simulations strongly suggest that the corona is heated by the dissipation processes occurring in a large number of small-scale current sheets. The corona is built up by a large number of small-scale current sheets where dissipative processes take place.

Based on these results, a series of 3D numerical simulations solving the simplified reduced magnetohydrodynamic (RMHD) equations, introduced by Strauss (1976), in Cartesian geometry were performed. The goal of these simulations was twofold: (1) to determine how a coronal loop responds to different photospheric velocity patterns, and (2) to investigate how the electric current sheets are formed, as well as the details of the heating that occurs through the dissipation of magnetic energy (Rappazzo et al. 2007, 2008, 2010; Rappazzo & Velli 2011; Rappazzo & Parker 2013; Rappazzo 2015). The results detailed in these papers confirmed the 2D results that the energy flux entering the corona because of photospheric motions causes a turbulent cascade of energy toward small scales, with electric current sheets continuously being created and destroyed throughout the coronal loop. The turbulence that develops is highly intermittent. The dissipation of energy occurs in current sheets, which are localized structures. Hence, the heating also occurs on small spatial scales. (Note that there is also particle acceleration, although that is beyond the scope of the present model.) When the loop is fully turbulent, it achieves a statistically steady state in which, on average, the Poynting flux induced at the boundaries by footpoint convection is balanced by the dissipation of energy in the electric current sheets. In this steady state, saturation occurs for the fluctuating magnetic energy, kinetic energy, mean square electric current, enstrophy, and other quantities that are then seen to fluctuate around their mean values. When the system becomes nonlinear, its behavior is highly dynamical and chaotic. This state occurs independently of the detailed form of the boundary velocity. A statistically steady state is achieved in which, although the magnetic field footpoints are convected continuously by the boundary flows, the resulting topology of the total magnetic field is not simply a mapping of those flows. By the same token, the turbulent dynamics that occur are due to the inherent nonlinear nature of the system, rather than being a consequence of the complexity of the magnetic field footpoint configuration.

For instance, a sheared boundary forcing constant in time (Rappazzo et al. 2010) or a boundary velocity made of distorted vortices constant in time (Rappazzo et al. 2008) generates nonlinear dynamics similar to that of distorted vortices changing in time on timescales larger than the coronal loop Alfvén crossing time (typical of X-ray-bright solar loops; Rappazzo et al. 2018). However, static large-scale boundary forcing does not change directly the mapping of the orthogonal magnetic field lines or create directly small-scale structures (and thus they do not force reconnection to occur); rather, once the orthogonal magnetic field intensity grows beyond a certain threshold, nonlinear dynamics develop because the coronal field is increasingly out of equilibrium in the vortices case (Parker 1994; Rappazzo & Parker 2013; Rappazzo 2015), or because an instability is triggered in the shear case (after which nonlinear dynamics is similar to the other cases, with the previous instability not occurring again).

In the fully developed stage, energy that is injected into the system at low spatial wavenumbers is redistributed by nonlinear interactions through MHD turbulent direct and inverse cascades. An inertial range develops in both the kinetic energy and magnetic energy. The spectra of these inertial ranges exhibit a power-law behavior. Note as well that the magnetic energy tends to be much larger than the kinetic energy. In configuration space the redistribution of energy has morphological consequences for both lower and higher spatial wavenumbers. The magnetic field comes to be organized in large-scale magnetic islands. These islands are separated by small-scale electric current sheets that extend along the axial magnetic field. Below a computationally very demanding threshold, but still small with respect to the coronal Reynolds numbers (viz., a magnetic Reynolds of about 800, corresponding to a linear grid resolution of about 8 km; in this manuscript we reach a resolution of about 15 km), average dissipation increases with increasing Reynolds number, as well as peak to average energy dissipation. Beyond this threshold average dissipation is independent of the Reynolds numbers (Rappazzo et al. 2007, 2008). The key point is that the energy input to the system is not solely due to the forcing at the boundary, i.e., the photospheric velocity, but also to the instantaneous state of the system, i.e., the configuration of the coronal magnetic field. Since the energy input is dependent on both the external forcing and the internal dynamics, the corona is a self-regulating forced system.

The recognition that the Reynolds number is very large in the corona, and the consequent continuous formation and disruption of current sheets where energy is dissipated, leads to a scenario that is qualitatively similar to the one proposed by Parker, namely, coronal heating is the observable signature of swarms of unobserved individual elementary events. Parker used order-of-magnitude arguments to estimate the magnitude and the timescale of production of the mean transverse component of the coronal magnetic field. In doing so, he neglected the internal dynamics of the system that are driven by the forcing at the boundaries. From a quantitative point of view, by spatially averaging over the entire system, as well as temporally averaging over a time that is long compared with the dynamical time, he derived the amount of energy liberated by magnetic dissipation in many "elementary events." The energy release that he computed was about 10¹⁷ J. He termed this averaged energy release event the "nanoflare." It is evident that the energy content of Parker "nanoflares" is an extreme upper limit of the energy content of an "elementary event."

While 2D MHD and RMHD do an excellent job of elucidating the dynamics of coronal heating, they cannot provide any information on the thermodynamical aspects. The energy due to ohmic and viscous dissipation is lost from the system since there is no accounting for thermal energy. Any effects due to thermal conduction and optically thin radiation are absent, and any feedback from the thermodynamics on the dynamics is neglected. The situation is reversed in 1D hydrodynamic models (Reale 2014). These models are often used to study the equilibrium state of heated loops, or the flows due to time-dependent evaporation or condensation flows. What is lacking in these 1D models is the important role of perpendicular dynamics. Hence, an ad hoc heating function is used to mimic the thermal energy deposition.

Simulation of the coronal heating problem requires a model that can provide a complete reproduction of the coronal heating energy cycle. Energy is injected in the corona through the convection of the magnetic field line footpoints. When the perpendicular magnetic field in the corona attains finite amplitude, intermittent turbulence develops in which nonlinear interactions lead to direct and inverse cascades. Energy dissipation is localized in electric current sheets in which magnetic energy is transformed into thermal energy (and fluctuating magnetic and kinetic energy as well). Thermal conduction efficiently transfers heat to the chromosphere, where it is lost by means of optically thin radiation, and may lead to chromospheric evaporation. To obtain this model of the full energy cycle, it is therefore crucial to augment the previously discussed numerical simulation models with an energy equation that accounts for thermal conduction and optically thin radiation. In addition, rather than an ad hoc "heating" term, the energy equation should include specific Joule heating and viscous heating terms. When gravitation effects are considered, it is also necessary to augment the momentum equation with pressure gradient and gravity terms, as well as adding in a density equation to account for the effects of stratification. This augmented system is implemented numerically in the HYPERION code, which solves the full set of compressible magnetohydrodynamic equations (Dahlburg et al. 2012, 2016a, 2018; see Section 3). HYPERION allows for the carrying out of simulations where footpoint shuffling heating causes heating due solely to the resistive and viscous dissipation induced in the corona. With the help of an atomic database (e.g., CHIANTI), temperatures and densities can then be used to synthesize observables, such as the emission of spectral lines, and construct differential emission measures, both useful to examine the radiative consequences of coronal heating.

Note that this added complexity comes with a cost—more variables and more equations imply that this system of equations is computationally more demanding. Hence, the range of Lundquist numbers that is achievable is somewhat less than for the 2D MHD and RMHD cases.

HYPERION is typically run with a chromospheric base of dimensions 4000 × 4000 km². The numerical simulations of our Cartesian coronal loop model neglect the geometry and topological complexity of the solar magnetic field but allow us to use higher numerical resolutions with respect to models that take into account entire active regions (e.g., Gudiksen & Nordlund 2002; Hansteen et al. 2010; Bingert & Peter 2011; Baumann et al. 2013; Hansteen et al. 2015; Guerreiro et al. 2015; Rempel 2017, and references therein), thus enabling a more detailed study of the energy dissipation processes at small scales and its thermodynamical and observational consequences that we describe in the following sections. The dynamics seen in the RMHD simulations is reproduced. At any given time a multitude of dissipative current sheets, which are elongated along the axial magnetic field, are observed to form (a visualization of these current sheets is found in Rappazzo et al. 2008). As for the temperature, the loop is found to be a multithermal system in which there is a spatial correlation between the electric current sheets and the temperature peaks. In general, high-temperature regions correspond to the locations of the electric current sheets. An important qualification is that the temperature, unlike the electric current, also depends on the density. Due to stratification, the loop is denser at its footpoints than at its apex. Hence, while the heating might be the same at both these locations, the temperature near the footpoints does not increase as much. Thus, we see that while the electric current sheets extend from one footpoint to the other, the high-temperature regions are more localized with respect to the loop apex. If the parameters of the loop, such as its length, magnetic field strength, or boundary velocity, are varied, then important features of the temperature distribution, such as the peak, extent, and duration, are found to vary.

The turbulent character of the heating has important implications for observations. The distribution of heating temperature in the loop is inhomogeneous. Furthermore, both spatial and temporal intermittency must be taken into account to interpret observations correctly; hence, significant heating occurs in only a small part of the loop volume at any given time. Thus, high-temperature plasma radiation originates from a very small fraction of the loop volume with a mass much lower than the total mass of the loop. These temperature structures cannot, at present, be resolved observationally owing to their small spatial extent. Using the computed number densities and temperatures, Dahlburg et al. (2016a) performed an emission measure analysis of their simulated loops. A loop length of 50,000 km and axial magnetic fields of 0.01, 0.02, and 0.04 T were examined. The computational emission measures derived from the simulated intensities were found to be in excellent agreement with emission-line intensities observed by the Extreme Ultraviolet Imaging Spectrometer (Culhane et al. 2007). Furthermore, in spite of the strong spatial and temporal intermittency of the loop system, the computed emission measures were found to retain the same form for the duration of the numerical simulations.

Dahlburg et al. (2018) studied the observational consequences for the temperature as the loop magnetic field strength is increased and/or the length of the loop decreases. This is equivalent to decreasing the Alfvén crossing (τ_A ) time of the loop. It was found that as τ_A decreases, both the peak temperature and the width of the temperature distribution of the loop exhibit increases. The emission measure broadening implies that the loop is more and more becoming a multitemperature system. Let L₀ denote the length of the coronal loop and B₀ denote the magnetic field magnitude of the loop. The predictions of how the effective temperature dependence varies with respect to the B₀/L₀ ratio are in agreement with the observational estimates of those quantities obtained by Xie et al. (2017), who examined 50 coronal loops in the <1–2 MK range. A larger observational temperature range is needed to determine the slope of the temperature dependence. Proportionality between temperature and emission width has been observed in coronal loops (Schmelz et al. 2014), but as yet there is no confirmation of the scaling with magnetic field strength.

The work of Rappazzo et al. (2018) elaborated the importance of the Alfvén crossing time (τ_A ) for the turbulent coronal heating process. They found that increasing the loop length, which decreases τ_A , leads to a decrease in the magnitude of the maximum electric current. This implies that the turbulent heating process is becoming less efficient. This conclusion was confirmed by Dahlburg et al. (2018), who found that when τ_p /τ_A is under a critical value of about 4 (recall that τ_p is the photospheric velocity time), the heating process will not produce enough radiation to permit, at present, any observation. This implies that long coronal loops with small magnetic fields will not be observed, unless some other heating process is present.

In this paper we present a more detailed investigation of the properties of elementary events by analyzing the data of numerical simulations of the same loop at two different resolutions. Not all elementary events have the same geometrical size, the same duration, and the same temperature peak, and we present their distribution in size, duration, and energy. In Section 2 we summarize the mechanism responsible for the onset of the cascade of energy to small scales and for the formation of the current sheets. We describe how to relate such a dynamics to the definition of an elementary event. In Section 3 we present the details of HYPERION and the initial and boundary conditions. In Section 4 Hyperion results in the case of B₀ = 0.01 T and L₀ = 50,000 km are analyzed. In Section 5 the statistical analysis performed on the data is used to derive the properties of the elementary events. We also present, for the sake of comparison, some results at higher resolution to show the trend that justifies the choices made in the analysis at lower resolution. Section 6 is devoted to discussion and conclusions.

The dynamics of the coronal loop are strongly influenced by the presence of a large-scale mean magnetic field (B₀). This field suppresses dynamics in the axial direction. Disturbances propagate into the system from the boundary along B₀ and cause the growth of perpendicular magnetic fields (b_⊥) and velocity fields (v_⊥). Dynamics in planes perpendicular to the axis then lead to the formation of electric current sheets. Note that the ratio of the rms b_⊥ to B₀ is much less than 1 and is also smaller than the aspect ratio. Therefore, an understanding of the energy dissipation can be obtained by studying a 2D section of the structure. The evolution of the system, in the limit of a large loop aspect ratio, can be modeled by using the RMHD equations (Strauss 1976). The RMHD equations are valid for plasmas with a small ratio of kinetic pressure to magnetic pressure, and in the limit of a small ratio of perpendicular magnetic field to axial magnetic field (b_⊥/B₀ ≤ l_⊥/L₀, where L₀/l_⊥ ≫ 1 is the loop aspect ratio). Consequently, the typical velocities are also sub-Alfvénic. In the RMHD model b_⊥ and v_⊥ depend on the axial coordinate (z). Their nonlinear interaction proceeds independently on different constant z surfaces along the loop. Long-wavelength (λ_A ≃ L) Alfvén waves provide communication across axial planes. These waves are the mechanism by which the energy originating from boundary motions propagates into the interior of the system.

The nonlinear dynamics do not depend on the pattern of the velocity forcing chosen to mimic photospheric motions, as appears evident by comparing the results obtained adopting very different photospheric velocity patterns. They have in common a photospheric velocity u_p = 1 km s⁻¹, the spatial scale l_p = 1000 km, and τ_A smaller than the photospheric convective timescale τ_p . Rappazzo et al. (2008) have adopted convection cell flow patterns constant in time, exciting all the wavenumbers k and l in the boundary planes included in the range $3\leqslant {({k}^{2}+{l}^{2})}^{1/2}\leqslant 4$. The results of the nonlinear dynamics induced by such forcing are very similar to those obtained by Rappazzo et al. (2010), where a photospheric velocity field in the form of a 1D shear flow pattern is present. This case is interesting because it is possible to follow step by step the physical mechanisms responsible for the evolution of the system. The current layers in the perpendicular planes formed during the linear phase are unstable to tearing modes (Furth et al. 1963), which are observed evolving during the linear stage. Once the system becomes fully nonlinear the dynamics is fundamentally different. The nonlinear terms do transport energy from the large to the small scales, and the magnetic field topology does not maintain any resemblance to the tearing-like instabilities.

Reduced MHD studies have shown that in the initial phase of development the perpendicular magnetic field grows linearly in time, obeying the relation b _⊥ = u _p t/τ_A . During this phase, b _⊥ is a simple map of the velocity forcing pattern at the photospheric boundary. The linear phase of growth continues until the perpendicular magnetic field attains a magnitude given by b_⊥ ∼ ℓ_⊥ B₀/L₀. This development occurs at a time given by t_ℓ ∼ τ_p > τ_A . At this time a cascade of energy toward small scales develops. The energy that is injected into the system at large spatial scales at the boundary is then balanced by the energy that is dissipated at small spatial scales in the interior of the system. The rms of b_⊥ and the Poynting flux also stop growing at this time, both fluctuating around a mean value in the subsequent nonlinear phase of evolution.

It is important to note that, for a given B₀, the Poynting flux across the photospheric boundaries depends only on b_⊥ and the photospheric velocity. The photospheric velocity is a boundary condition for the system, the average magnitude of which is fixed. The perpendicular magnetic field b_⊥, however, initially increases, leading to a similar increase in the Poynting flux. The Poynting flux then oscillates around a mean value in the subsequent nonlinear phase. Hence, the energy available to form the current sheets that heat the loop is determined by the maximum value attained by b_⊥. Since b_⊥ ∼ ℓ_⊥ B₀/L₀, this energy is ultimately dependent on the values of L₀ and B₀, the two significant parameters governing the perpendicular dynamics of the loop. As the loop parameters vary, different regimes of turbulence can develop. Strong turbulence is found for weak axial magnetic fields and long loops, leading to Kolmogorov-like spectra in the perpendicular direction, whereas weaker regimes (steeper spectral slopes of total energy) are found for strong axial magnetic fields and short loops (Rappazzo et al. 2007, 2008). Additionally, they found that the scaling of the heating rate with B₀ is a power law that depends on the spectral index of total energy, and therefore its index changes with the magnetic field intensity from ${B}_{0}^{3/2}$ for weak fields to ${B}_{0}^{2}$ for strong fields (for loops with same aspect ratio).

Parker recognized that the transverse component of the magnetic field b_⊥ plays a key role in the evolution of the system. In particular, he examined in great detail the inclination angle (θ) of the mean b_⊥ to the mean field direction due to photospheric motions. However, based on factors discussed earlier in this paper, no relationship exists between the inclination angle (θ) and the rate of magnetic energy dissipation (see, e.g., Galsgaard & Nordlund 1996).

The work of Dahlburg et al. (2016a) and Dahlburg et al. (2018) confirmed the dynamics described with the RMHD approach. In addition, this research was able to determine how the coronal plasma heats up and radiates energy as a consequence of photospheric convection of magnetic footpoints. As previously stated, the loop is found to be a multithermal system in which there is a spatial correlation between the electric current sheets and the temperature peaks. High-temperature regions are found to correspond to the locations of the electric current sheets and extend axially along the loop magnetic field. These isolated hot regions are surrounded by much larger regions of cooler plasma. The sites of enhanced temperature represent the elementary events defined above. They represent the building blocks of what is called "coronal heating," which is actually the superposition of the radiation emitted by all elementary events, averaged over timescales and spatial scales much longer than the typical dimensions and lifetime of a single elementary event. In the next sections we will discuss the extent, duration, and energetics of these elementary events.

The HYPERION code is a computational representation of the corona. In this representation the corona is modeled as a square cylinder (a figure that represents the computational box is discussed in Section 4) . Along the z-direction a DC magnetic field is applied. The boundaries in the z-direction represent the upper chromosphere. At each z boundary line-tied boundary conditions are enforced, and the footpoints of this magnetic field are convected by applied flows. Periodic boundary conditions are enforced in the x- and y-directions. Hence, the interior of the computational box represents the corona, while the z boundaries represent the upper chromosphere.

3.1. Governing Equations

As noted in the introduction, it is important to simulate all of the processes in the coronal energy cycle. Magnetic reconnection is a significant process, and hence the resistive terms must be included in the magnetic induction equation to ensure that it occurs correctly. The Joule heating term is needed in the energy equation to ensure that the thermal energy created by magnetic reconnection is deposited in the correct location and in the proper amount. Thermal conduction and optically thin radiation are needed to obtain the effects of the chromospheric response to heating, such as evaporation. In addition, compressibility must be factored in to allow for stratification effects due to gravitation. Hence, the equations used in the simulation model should allow for all of these effects. The HYPERION code evolves in time the 3D compressible magnetohydrodynamic equations. Viscous and resistive dissipation terms are included, as well as temperature-dependent thermal conduction and optically thin radiation. The governing equations are as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial n}{\partial t}=-{\rm{\nabla }}\cdot (n{\boldsymbol{v}}),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\frac{\partial n{\boldsymbol{v}}}{\partial t} & = & -{\rm{\nabla }}\cdot (n{\boldsymbol{vv}})-\beta {\rm{\nabla }}p+{\boldsymbol{J}}\times {\boldsymbol{B}}+\frac{1}{{S}_{v}}{\rm{\nabla }}\cdot {\boldsymbol{\zeta }}\\ & & +\ \frac{1}{{{Fr}}^{2}}\ n{\rm{\Gamma }}(z){\hat{e}}_{z},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\frac{\partial T}{\partial t} & = & -{\boldsymbol{v}}\cdot {\rm{\nabla }}T-(\gamma -1)({\rm{\nabla }}\cdot {\boldsymbol{v}})T\\ & & +\frac{1}{n}\left\{\frac{1}{{\Pr }\,{S}_{v}}\left[{\boldsymbol{B}}\cdot {\rm{\nabla }}\left({\kappa }_{\parallel }\ {T}^{5/2}\ \frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}T}{{B}^{2}}\right)\right]\right.,\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&+\displaystyle \frac{(\gamma -1)}{\beta }\left[\displaystyle \frac{1}{{S}_{v}}{\zeta }_{{ij}}\displaystyle \frac{\partial {v}_{i}}{\partial {x}_{j}}+\displaystyle \frac{1}{S}{\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}^{2}\right.\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\left.\left.-\displaystyle \frac{1}{{P}_{{rad}}{S}_{v}}{n}^{2}{\rm{\Lambda }}(T)+\displaystyle \frac{\beta }{(\gamma -1)}{{nC}}_{N}\right]\right\},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times {\boldsymbol{v}}\times {\boldsymbol{B}}-\displaystyle \frac{1}{S}{\rm{\nabla }}\times {\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{B}}=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&p={nT}.\end{eqnarray} \tag{ 8 }$$

As is typical in turbulence calculations, these equations are written in dimensionless form. The dimensional terms are lumped together to form dimensionless numbers. The dimensionless numbers are the viscous Lundquist number, S_v = n_* m_p V_A* L_*/μ (m_p = 1.673 × 10⁻²⁷ kg is the proton mass); the (resistive) Lundquist number, S = μ₀ V_A* L_*/η (μ₀ = 1.256 × 10⁻⁶ H m⁻¹ is the magnetic permeability); $\beta ={\mu }_{0}{p}_{* }/{B}_{* }^{2}\equiv$ boundary pressure ratio; the Prandtl number, ${\Pr }={C}_{v}\mu /{\kappa }_{\parallel }{T}_{* }^{5/2}$; and the radiative Prandtl number, P_rad , where $\mu /{\tau }_{A}^{2}{n}_{* }^{2}{\rm{\Lambda }}({T}_{* })$. C_v is the specific heat at constant volume, and the magnetohydrodynamic Froude number ${Fr}\,={V}_{A}/{({{gL}}_{* })}^{1/2}$, where the solar surface gravity is given by g = 274 m s⁻².

The nondimensional variables are defined in the following way: the number density is n( x , t), the flow velocity is v ( x , t) = (u, v, w), the thermal pressure is p( x , t), the magnetic induction field is B ( x , t) = (B_x , B_y , B_z ), the electric current density is J = ∇ × B , and the plasma temperature is T( x , t), Also note that the viscous stress tensor is ζ_ij = μ(∂_j v_i + ∂_i v_j ) − λ∇ · v δ_ij , the strain tensor is e_ij = (∂_j v_i + ∂_i v_j ), and the adiabatic ratio is γ.

The characteristic values of various quantities are chosen to represent values typical of the upper chromosphere. These characteristic values are then used to render the equations dimensionless. These values include the characteristic number density (n_*), the characteristic temperature (T_*), and the characteristic parallel Alfvén speed (V_A*). The perpendicular box width (L_*) is also used (here chosen to be 4 × 10⁶ m). Time (t) can then be measured in units of the dimensionless Alfvén time τ_A , which is given by τ_A = L_*/V_A*.

The shear viscosity and magnetic resistivity (μ and η) are spatially uniform and temporally constant. The bulk viscosity (λ) is set equal to (2/3)μ, i.e., the Stokes relationship is enforced. The term κ_∥ denotes the parallel thermal conductivity. κ_⊥, the perpendicular thermal conductivity, is neglected. The CHIANTI atomic database (Dere et al. 1997; Landi et al. 2012) is used to construct the radiation function. The base temperature of the loop (T_* = 20,000 K) is used to normalize the radiation function. An elliptical model is used to construct the loop gravity profile (Γ(z)). Near the bases of the loop Newton cooling (C_N ) is effective (Dorch & Nordlund 2001). Let T_i (z) denote the initial temperature profile; at the upper boundary we enforce ${C}_{N}=10\,[{T}_{i}(z)-T(z)]{e}^{-4(0.5{L}_{z}-z)}$, while we set ${C}_{N}=10\,[{T}_{i}(z)-T(z)]{e}^{-4(z+0.5{L}_{z})}$ at the lower boundary.

The numerical method used to solve employs a Fourier collocation finite-difference scheme in space. Time step splitting is employed in time. Super-time-stepping is used for thermal conduction (Meyer et al. 2012). All other terms are advanced in time with a low-storage fourth-order Runge–Kutta scheme (Carpenter & Kennedy 1994). More details on the numerics also can be found in our previous papers (Dahlburg et al. 2016a, 2016b, 2018).

3.2. Initial and Boundary Conditions

3.3. Initial Temperature, Number Density, and Gravity Specification

Heating is a consequence of magnetic reconnection in elementary events. Starting from an initial uniform potential magnetic field (B₀), it takes some time to build up a large-enough perpendicular magnetic field to permit this to happen. If thermal conduction and radiation are turned on from the beginning of the simulation, the initial temperature profile is significantly altered. In the absence of a heating function an energy balanced equilibrium cannot be specified. Hence, we begin with an ideal thermodynamic equilibrium and ramp up the thermal conduction while the magnetic stresses build up. Here the time to ramp up from zero to the full value is chosen to be 5 minutes. The initial temperature profile (T_i ) is set to 2 × 10⁴ K at the boundaries and 10⁶ K in the center (for the nondimensional values see Dahlburg et al. 2016a).

The loop gravity is given by an elliptical model that is described in Dahlburg et al. (2016a). With the initial temperature and gravity specified, the corresponding number density can be found numerically with a shooting method, as described in Dahlburg et al. (2016a). In this paper we use a dimensionless loop length consistent with a dimensional loop length of 50,000 km. We set a dimensionless value of 1 for the number density at the loop boundary. The shooting method then gives a dimensionless number density value of 2.298 × 10⁻⁴ at the loop apex. In dimensioned terms, the initial number density is 1.0 × 10¹⁸ m⁻³ at the loop boundary and 2.298 × 10¹⁴ m⁻³ at the loop apex. The profiles of the initial number density and temperature used in both numerical simulations are shown in Figure 1, in which they are restored to dimensional form.

Zoom In Zoom Out Reset image size

3.4. Selection of Parameters

In this section we describe the algorithms that we developed to identify elementary events in the HYPERION simulation runs. Given our experience studying the dynamics generated by footpoint-driven motions that braid the field, forming elongated current sheets that are followed by a sudden release of energy and increases in temperature, we expect heating events to be transient and localized changes in the MHD properties of the plasma.

In Section 2 we defined elementary events as the temperature enhancements associated with current sheet formation. A temperature increase is just a consequence of the heating, and therefore this definition does not allow us to call them elementary heating events from the outset. We will argue throughout this paper, however, that elementary events identified in temperature space can be a good proxy for certain properties of the elementary heating events, which otherwise are better described by terms such as the Joule heating.

There are several reasons to choose temperature over other plasma properties to identify elementary events in the corona: (1) the temperature increase is directly connected to the heating exchange that we are trying to investigate, (2) the temperature response of the plasma is also tied to the radiative response and therefore observational diagnostics, and (3) the temperature maps are smooth and easy to be interpreted by the algorithms we have developed. Density is important in the radiative response, but it is not as closely connected to heating because of its slower response (e.g., Reale 2014). For completeness we describe at the end of this section the results from identifying events as transients of the Joule heating term, and we discuss the similitudes and differences with the temperature-inferred results. While the Joule heating term, ηJ², directly measures the energy release, this term cannot be measured observationally. The relationship between the current and radiative emission is not straightforward. Furthermore, identifying events using this term presents challenges to the performance of the identification algorithms, which make it best to address it after the temperature case.

Identifying events in space and time within the simulation box involves locating neighboring pixels with a shared evolutionary history for a relevant property. Machine learning, and particularly clustering algorithms, have been a popular and efficient technique used to address this problem in a variety of contexts, including solar physics (Caballero & Aranda 2013), where automated data processing methods have blossomed in the past decade with the high data volumes of the Solar Dynamics Observatory (e.g., Aschwanden 2010; Martens et al. 2012; Bobra & Mason 2018). We describe below two clustering approaches that we have followed in assessing our ability to uniquely identify and characterize elementary events in space and time. We present both to show the robustness of the results. The first algorithm uses a well-known clustering technique called density-based clustering of applications with noise (DBSCAN). The second is our own developed brute-force cluster finder that relies on a distance search instead of the pixel search that characterizes our DBSCAN version. While the working philosophy is similar, definitions and search strategy make them distinct test beds for this experiment, and the comparison gives us some insight about the uncertainties of the process and the results.

4.1. Density-based Clustering (DBSCAN)

DBSCAN is a data clustering technique proposed by Ester et al. (1996). DBSCAN finds clusters of points within a sample with the assumption that a minimum number of neighbors are within a distance radius of a given point in the cluster, i.e., the density of points in the neighborhood is higher than a certain threshold.

We have developed an Interactive Data Language (IDL) algorithm that applies this technique to any n × m array, where n is the number of variables and m the number of observations, using as reference the definitions, assumptions, and algorithm description in Ester et al. (1996) and an existing Python implementation. ⁵ The search evaluates the number of points in the sample that are within a predefined neighborhood radius (in pixels) of a given point and makes decisions about whether that point is part of a new or existing cluster. If the number is equal to or larger than the required threshold (N_pts), then all the identified points are considered part of the cluster. Core points in the cluster are those that satisfy this condition about their neighborhood. Border points are those that have less than N_pts in their neighborhood but belong to the cluster because of the core point they are connected to. The algorithm evaluates the neighborhood of every point in the sample and stops when all have been classified as cluster or noise. Noise points do not have any core cluster neighbors in their . For the purpose of our study, noise points are discarded.

In our case n = 3, the three spatial dimensions of the simulation box in pixel units, and m is the number of points per time step that we aim to classify. As stated before, we make identifications based on the maximum temperature. We use a hard threshold (τ) to locate and extract those m locations and then apply DBSCAN to return the number of clusters within that set and the unique label for all points in each cluster. This is performed for all time steps. The results that we present here are made with = 1.5, N_pts = 2. We experimented with various temperature thresholds and settled with τ_{log T} = 6.10. That temperature captures temperature events that are visually discernible in a cross section at the loop apex, keeping the number of points and clusters in a manageable range to perform the temporal tracking and to compute their properties. This definition makes it more likely to detect events at the apex, because the apex of the loop generally has the highest (coronal) temperature. Additionally, the apex typically has the lowest density, so events occurring at this location will in general release more energy per particle and thus cause the largest temperature deviations. We should note that other temperature event definitions are possible, for instance, one that, instead of focusing on coronal events, looks for departures of temperature from a steady state. Figure 2 shows the cross section and reference thresholds for three time steps.

Zoom In Zoom Out Reset image size

Once we have identified all clusters of points, i.e., our definition of elementary events, in all time steps, we apply a second algorithm, which we call here "cluster tracker," that determines whether clusters from consecutive time steps are the same and should have the same label. The algorithm was also developed in IDL and tracks the events forward in time, by determining whether there is a spatial overlap between the clusters in consecutive time steps. Events that appear as a single cluster in a given time step can then split in multiple fragments as they evolve in time and retain a unique label, but events that form from coalescing clusters keep their independence before the merger. The cluster with the larger size dictates the label for the cluster going forward.

Figure 3 shows in the two leftmost panels the temperature isosurfaces at τ_{log T} = 6.10, together with the event identification. The algorithm correctly identifies four data clusters corresponding to four heating events. A movie version of the figure is available online. The movie goes with the paper.

Download figure:

Video Standard image High-resolution image

4.2. Physical Distance Clustering (PDC)

Being able to identify and track clusters of points in time allows us to determine the properties of the events as they evolve in the simulation run. As the HYPERION numerical solution provides the changes in temperature, density, velocity, current, etc., for every one of the pixels, we can compute the characteristic properties of the set of pixels that conform every event in our identification. In this section we present the results from this analysis. First, we present the properties of heating events identified in the temperature domain. Then, we discuss how the results compare when looking into the Joule heating space.

5.1. Temperature Events

Figure 4 shows the temperature changes in time along the Z-direction at the location of a temperature enhancement highlighted in Figure 2 with a blue circle. This is the first event at that location in this run, and the pre-event temperature profile is very similar to the profile at the start of the run. The heating resulting from the forcing at the footpoints produces a temperature spike of approximately 1.55 MK, offset toward the top boundary of the box (positive Z). This temperature profile is not the temperature distribution along an isolated field line, but it is representative of the general temperature changes observed along the box, i.e., along the loop, when a heating event takes place. The figure also shows the threshold level that sets the definition of a heating event in this analysis. Only voxels with temperatures above that threshold survive the identification process and temporal tracking depicted in Figure 3. During the evolution of the event, the density exhibits small departures from its initial profile.

Zoom In Zoom Out Reset image size

For every identified event we computed the total volume as the sum of all voxels, the extent in the three spatial dimensions, the duration, the centroid position, the maximum temperature and density, and several other derived quantities, such as aspect ratio, thickness, or internal energy. These quantities were calculated for every instant in the lifetime of an event. In Figure 5 we show the summary of these properties for the events identified by the two clustering algorithms. We present these results in CGS units, as opposed to the SI units of previous figures, to facilitate comparisons to observables that are usually reported in those units. In the case of volume, aspect ratio, centroid, thickness, and energies the values correspond to the time step of maximum temperature in each cluster, which marks the end of the heating. The purpose of showing the comparison of both algorithms is to illustrate the differences that the choice of algorithm has on the results. While DBSCAN and PDC identified a different number of events, 359 and 205, respectively, the overall distributions are very similar. The larger differences are in volume, as PDC sometimes clusters together events separated by DBSCAN, and internal energy, which is highly correlated to volume.

Zoom In Zoom Out Reset image size

The left column shows the distribution in the maximum density and temperature reached by any point within the cluster. The distribution is dominated by low-temperature and low-density events, close to the temperature threshold, suggesting that as we lower that threshold the number of events will increase and their peak temperature will go down. Figure 2 already hints that we are missing events with our choice of threshold. To verify this intuition, we run both algorithms with a threshold of τ_{log T} = 6.07, and the trend was confirmed as shown in Figure 6 for DBSCAN. The number of events increased to 555. We had difficulties working with lower thresholds because the increased number of points creates challenges to the "cluster tracking" part of the workflow. We are working on making the algorithms more efficient for future studies.

Zoom In Zoom Out Reset image size

One of the most fundamental quantities of the heating events that we want to know about is energy. In Figure 5 we show the maximum internal energy (U = κ_B N_e T_e ) reached by all clusters, computed as the sum of all voxels. This quantity, however, is only indicative of the energy in the system at this moment in time. It does not tell us what was gained with the most recent heating event. For that we calculate an internal energy increment from the difference in internal energy between the energy at peak temperature and the start of the event:

$$\begin{eqnarray}&&{\rm{\Delta }}U=\sum _{{ijz}}{\kappa }_{B}{N}_{e}{T}_{e}{\left|{}_{t({T}_{{\max }})}-\sum _{{ijz}}{\kappa }_{B}{N}_{e}{T}_{e}\right|}_{t=0},\end{eqnarray} \tag{ 9 }$$

where i, j, z are the spatial indices of the voxels that define the cluster at maximum temperature. This measure of the Joule heating is strictly less than the total energy release. This is because some of the energy is radiated away, conducted away, or advected away from the site of release. However, we argue that, despite this limitation, it is a useful proxy for the energy release because the bulk of the released energy goes into heating the plasma in the loop. In Section 5.2 we discuss further the validity of this assumption. In Section 5.3 we discuss the shape of the energy distributions.

Figure 5 shows the distribution of the cluster centroid positions along the loop (Z centroids), which tells us that the temperature events are distributed in both directions around the loop apex, but very close to it. This is not surprising given our choice of temperature as a proxy for heating. Coronal temperatures are more likely to be detected at the apex. As there can be differences between where the temperature peaks and where the heating is deposited, we examine identifying events with the current distribution in the next section. We also computed the thickness of the structures at the apex plane. The thickness is defined here as the minimum value of the median thicknesses in X and Y. This preserves the smaller dimension of the two in an elongated (e.g., elliptical) cross section. The duration of the events is shown in the rightmost column and reveals that most events are less than 100 s long, and about half that time for the time to reach the maximum temperature.

5.2. Joule Heating Events

Identifying the events in temperature space is relatively easy because the temperature distribution changes rather smoothly within the box and the clustering algorithm has no difficulty in identifying concentrations of points. The temperature changes, however, are not a direct diagnostic of heating, but a consequence of it. Temperature can change through direct heating, but also via thermal conduction from heating deposited elsewhere.

For this reason, we look here at the differences that arise when identifying events from the distribution of the Joule heating term. This quantity is directly connected to the energy released by the system and correlates with the Poynting flux changes (Dahlburg et al. 2016a).

For the identification we choose a threshold value of τ_Joule = 1.0 × 10⁻⁵ W m⁻³ that corresponds to 2.7 times the standard deviation σ of the electric current distribution (j_z = 0.159 A m⁻¹ and η = 5.59 × 10⁻⁵ ohm m), shown in the top panel of Figure 7 for one of the time steps. Figure 3 shows in the two rightmost panels the Joule heating isosurfaces at that threshold level and the DBSCAN event identification assuming = 10, N_pts = 10. The Joule heating term outlines the location of the current sheets where reconnection takes place. The currents extend all along the box, and therefore the clustering algorithm identifies events that are longer than the temperature events. Their configuration is also more complex. Being thinner, longer, and with peaks and troughs along their structure, it results in a fractionation of the long sheets that are sometimes identified as multiple events. The online movie version of Figure 3 shows the detections as a function of time.

Zoom In Zoom Out Reset image size

This discretization of the structures is in part due to the noise introduced by the high Reynolds number with respect to the resolution at 128 × 128 × 192. A better-resolved (see Servidio et al. 2010; Wan et al. 2010, for a more complete discussion on this topic) run with twice the resolution across (256x × 256 × 192) shows smooth coherent structures running from footpoint to footpoint (Figure 8). The clusters are better sampled and easier to track by our identification algorithms. In this particular case we looked at a threshold τ_Joule = 3.0 × 10⁻⁶ W m⁻³, also at the 2.7σ level (j_z = 0.123 A m⁻¹ and η = 2.79 × 10⁻⁵ ohm m). The corresponding current distribution is shown in the bottom panel of Figure 7. Unfortunately, the duration of this run did not provide sufficient statistics to be compared to the temperature or the Joule heating low-resolution case. This will be the subject of a future investigation. The observed trends reveal that the high-resolution case exhibits thinner current sheets and higher current peaks than the low-resolution case. This point was also noted by Baumann et al. (2013), who similarly found that current sheet widths depend on numerical resolution. The remainder of this section discusses the results of the low-resolution case.

Zoom In Zoom Out Reset image size

Figure 9 shows the Joule heating changes in time along the Z-direction at the same location as Figure 4 for the temperature. There are similarities in that the temperature enhancement takes place toward the top of the box (positive Z), which is also where the increase in Joule heating happens. The peak in the Joule heating term takes place nearly simultaneously (sometimes before in other events) to the peak in temperature, which is something that we had already seen in HYPERION for the electric current and maximum temperature (Dahlburg et al. 2016a).

Zoom In Zoom Out Reset image size

The properties for these events are shown in Figure 10. The most striking differences have already been described. The Joule heating events are longer, and therefore their aspect ratio numbers are higher. The events are thinner, as expected for the size of the current sheets, although that is dependent on the choice of threshold. The centroid positions along the Z-direction show that, in contrast to the temperature events, the Joule heating events are uniformly distributed along the length of the box. The figure also shows that there are no significant differences in estimating the duration of the events when they are identified in temperature and Joule heating space. The maximum density in the Joule heating events shows a peak at 10¹² cm⁻³. This is not surprising because the events are extending to the footpoints, where the density is higher. In the high-resolution case, where events are extended and do not break up into smaller fragments, we expect that distribution to skew almost completely to higher densities.

Zoom In Zoom Out Reset image size

In terms of the energy distribution, our results indicate that there are differences between the internal energies estimated for temperature events and those estimated for Joule heating events. While both cover the same energy range, the distributions are noticeably different. Section 5.3 discusses the shapes of the distributions in more detail. The differences can be explained in terms of volume. The internal energy is highly correlated with volume, and the volume distributions show similar discrepancies. The estimated volume depends on the number of voxels in the cluster, which in turn depends on the thresholds. While the temperature events seem to be a good proxy of where, across the flux tube, the heating is taking place, and for how long, the volume highlighted by the temperature increment does not correspond to the volume where the heating takes place.

We find, however, that the internal energy increment (ΔU) does not appear to be a bad proxy for the heating of the events, provided that we know where the heating is deposited. Figure 11 shows that for Joule events ΔU, already shown in Figure 10, compares well with the actual time-integrated Joule heating energy available in the clusters (E_j ). This distribution gives us the actual heating energies for the events in this run, which are capable of producing temperatures up to 4 MK, with the bulk under 2 MK. The distribution of ΔU for the temperature events shows a flatter distribution, particularly for a higher threshold in temperature (6.10), but reducing the threshold (6.07) and therefore reaching to lower-energy events skews the distribution in the direction of those low energies and brings it closer to the energy distribution in the Joule events.

Zoom In Zoom Out Reset image size

5.3. Kolmogorov–Smirnov Tests

We discuss here the shape of the energy distributions in some detail. The histograms visually appear consistent with lognormal distributions or perhaps power laws, so in order to determine their consistency with these distributions, we have performed Kolmogorov–Smirnov (K-S) tests. K-S tests work by comparing the cumulative distribution functions of the empirical distribution to a theoretical distribution and measuring their maximum separation D_max. If D_max is greater than a critical value, then we reject the null hypothesis and conclude that the empirical distribution is not consistent with the theoretical one.

To begin, we first perform maximum likelihood estimation (MLE) of the parameters of the distributions: mean μ and standard deviation σ for a lognormal distribution, and the minimum value x_min and slope b for a power law. In the former case, the MLE of the mean and standard deviation can be calculated directly. In the latter case, the slope b can be calculated directly for some assumed x_min, but the x_min value itself needs to be calculated through other means. To calculate x_min, we follow the method described by Clauset et al. (2007), which is a K-S minimization technique. For each value of x_i in the data set, assume that x_min = x_i , calculate the MLE of the slope, and then perform a K-S test on that fitted power law for all ${x}_{i}\geqslant {x}_{\min }$. The power law with the smallest D_max gives the best estimate of x_min. This method does not provide uncertainties for the power-law parameters, however.

The MLE parameters and K-S test results are shown in Table 1. The result depends on the variable used for clustering, the threshold, and the clustering method. When clustering based on temperature, DBSCAN tends to find distributions consistent with power laws, while PDC finds distributions consistent with lognormal. When clustering based on the Joule heating, the distributions appear to be neither consistent with power laws nor lognormal. Visually, these distributions do appear to have a noticeable skew to them. Because the results depend intimately on the choices in the clustering methods, we emphasize again that the derived values of heating are only estimates of the true energy release.

Table 1. MLE Parameters of the Empirical Distributions of Event Energies and K-S Tests to Determine Consistency with Lognormal Distributions (top) or Power-law Distributions (Bottom)

Case	μ	σ	Dmax	α1%	α5%	α10%	K-S Test Result
Lognormal, DBSCAN, ΔU ($\mathrm{log}T\gt 6.07$)	19.96 ± 0.08	1.33 ± 0.05	0.089	0.069	0.058	0.052	Reject H0
Lognormal, DBSCAN, ΔU ($\mathrm{log}T\gt 6.10$)	20.22 ± 0.08	1.33 ± 0.06	0.066	0.086	0.072	0.065	Accept H0 at 1% and 5% levels, reject 10%
Lognormal, PDC, ΔU ($\mathrm{log}T\gt 6.07$)	20.49 ± 0.14	1.55 ± 0.10	0.075	0.101	0.084	0.076	Accept H0
Lognormal, PDC, ΔU ($\mathrm{log}T\gt 6.10$)	20.69 ± 0.14	1.40 ± 0.11	0.071	0.114	0.095	0.085	Accept H0
Lognormal, DBSCAN, ΔU (J events)	18.89 ± 0.06	1.19 ± 0.04	0.081	0.058	0.048	0.043	Reject H0
Lognormal, DBSCAN, EJ (J events)	19.29 ± 0.04	0.87 ± 0.03	0.143	0.053	0.044	0.040	Reject H0
Case	xmin	b	Dmax	α1%	α5%	α10%	K-S Test Result
Power law, DBSCAN, ΔU ($\mathrm{log}T\gt 6.07$)	2.03 × 1018	1.24	0.427	0.072	0.060	0.054	Reject H0
Power law, DBSCAN, ΔU ($\mathrm{log}T\gt 6.10$)	1.47 × 1021	1.59	0.109	0.162	0.135	0.122	Accept H0
Power law, PDC, ΔU ($\mathrm{log}T\gt 6.07$)	8.64 × 1020	1.42	0.181	0.162	0.132	0.121	Accept H0 at 1% level, reject at 5 and 10%
Power law, PDC, ΔU ($\mathrm{log}T\gt 6.10$)	2.27 × 1020	1.35	0.196	0.150	0.125	0.113	Reject H0
Power law, DBSCAN, ΔU (J events)	8.95 × 1018	1.40	0.154	0.090	0.075	0.068	Reject H0
Power law, DBSCAN, EJ (J events)	8.81 × 1018	1.50	0.066	0.071	0.059	0.053	Accept H0 at 1% level, reject at 5 and 10%

Download table as:  ASCIITypeset image

In this paper we report the results of a detailed analysis of the data obtained running the HYPERION code on a specific case that represents a model of a typical coronal loop. We have assumed a loop base density n_base = 10¹⁷ m⁻³ (producing a loop apex density n_apex ≈ 3×10¹⁴ m⁻³), a loop base temperature T_base = 2 × 10⁴ K, perpendicular dimension L_⊥ = 4 × 10⁶ m, and a loop length of L₀ = 50,000 km. A basic loop magnetic field strength of B = 0.01 T is used. The dimensionless timescale of the forcing, t^*, is set to a value that corresponds to a 5-minute convection timescale, and the normalized driving velocity V_phot is 10³ m s⁻¹. Simulations were performed at two numerical resolutions. For the first simulation, the resolution is 128 × 128 × 192 (the 128 refers to resolution in the planes, and the 192 refers to the central difference grid points in the axial direction). At this resolution we were able to run the code for almost 1 hr of real time, a time sufficient to produce a large enough number of events to allow a statistical analysis. At this resolution the physical grid size is 31 km in the perpendicular direction and 260 km in the parallel one. We also performed a second simulation, at lower physical dissipation, with a resolution of 256 × 256 × 192. This simulation was run for a physical time of about 15 minutes, which was long enough to verify that the results show the same sort of trends as have been seen in the RMHD simulation as the resolution is increased.

At the beginning of the simulation the loop is in hydrostatic equilibrium, as specified in Section 3, and the thermodynamical equilibrium is ensured by ramping up conduction and radiation from zero, so that the loop temperature does not collapse before the reconnection-driven heating mechanism becomes efficient within the elementary events. Note that there is no "ad hoc" heating function—the only heating is due to the dissipation of the currents induced by the photospheric motions. The resulting loop is a multitemperature system with peaks over 4 MK within the elementary events and lows well under 1 MK elsewhere.

Adopting the machine learning algorithms described in the previous sections, we have been able to identify elementary events occurring continuously during the evolution of the system driven by the boundary motions. We have measured their spatial distribution and other important properties, namely:

• 1.  
  the size and duration of the events;
• 2.  
  the energy released, evaluated both as the variation of internal energy and as the energy due to current dissipation within the events;
• 3.  
  the variation of density with respect to the initial conditions.

The results detailed in this paper confirm the fundamental role of elementary events in the process of dissipating energy transferred from the photosphere into the corona. The simulations also demonstrate that the release of energy in the corona occurs through the dissipation of currents localized in current sheets resulting from the nonlinear coronal dynamics triggered by the large-scale photospheric forcing. Their existence appears clearly in Figures 3 and 8, which show (at both resolutions) that they extend along the entire loop and occupy a volume very small compared with the volume of the loop. From Figures 5, 6, and 10 it appears that their thickness is of a few pixels, which in dimensional units means less than 100 km. The thickness extending for a few pixels remains true also at higher resolution, less than 50 km. This fact can be easily understood considering that the current sheets are the manifestation of the energy cascade toward small spatial scales induced by the nonlinear effects dominating the perpendicular dynamics (Rappazzo et al. 2008). Typical of turbulent cascades, the spectral energy flow from large to small scales is stopped at the dissipative scale where dissipation dominates and provides a sink for the energy flow: here the dissipation rate equals the nonlinear transfer rate, leading to a dissipative scale (and therefore current sheet thickness) decreasing as a power of the Reynolds number (e.g., Biskamp 2003). The conclusion is that the thickness of the current sheets decreases with increasing resolution and in the solar conditions is far smaller than any possible observational resolution.

The current sheets form and disrupt, generating a localized increase of the plasma temperature (the elementary event), on timescales determined by the nonlinear dynamics induced by the growth of the component of the magnetic field (b_⊥) perpendicular to the main field. Figures 5, 6, and 10 show that the duration of each event is very rarely longer than 100 s, up to 200 s, and in most cases varies between a few and 50 s.

In the simulations the energy dissipated in the current sheets is transformed into internal energy through Joule heating. In this paper we measure the amount of energy involved in each event in two ways: (1) by determining the variation of internal energy during the occurrence of the event, and (2) by computing the Joule heating, ηj², created within the event by the current dissipation. We have identified the currents within the current sheets as those that form the non-Maxwellian tail. This is clearly seen in Figure 7, which shows the probability distribution function (pdf) of electric currents calculated at one time during the evolution. The variation of the internal energy method and the Joule heating method are in good agreement in the evaluation of the energies associated with the elementary events, exhibiting a distribution centered around 10¹⁸–10²¹ erg (see Figures 5, 6, and 10). Few events are found at higher energies up to 10²³ erg.

Looking at the pdf's of the electric current density, it is interesting to note that both the amount of currents in the tail of the distribution and their values increase considerably with resolution. At the higher resolution shown in this paper the currents are more concentrated within the currents sheets, and also the number of currents sheets increases and their thickness decreases. We expect this trend to continue when resolution and the corresponding Reynolds number are further increased (see Rappazzo et al. 2017), keeping the total energy associated with the events in the same range of 10¹⁸–10²¹ erg. This point must be verified by running much higher resolution simulations.

We also have found a good correlation between the energies, durations, and volumes of the events. The few events in the high-energy tail of the distribution are also in the long-duration range and big volume range.

As mentioned in Section 5.2, Joule heating events extend all the way to the bottom of the loop, where the density is the typical density above the chromosphere. They are much longer than the temperature events because of the density stratification; nevertheless, the much higher density at the bottom does not allow us to raise the temperature as significantly as at the loop apex with a similar local current dissipation rate. They also are thinner than the temperature events, a property related to the magnetic connectivity associated with current sheets. Recently, Rappazzo et al. (2017) have traced at a given time magnetic field lines from current sheets in a high-resolution reduced MHD simulation of a coronal loop with geometry and parameters similar to those shown in this paper. They have found that the field lines traced from such current sheets (elongated in the axial direction) fill a larger volume around the current sheets. As thermal conductivity is highly anisotropic and heat transport occurs essentially along magnetic field lines, it is clear that at each instant of time heating from current sheets is transported in these regions around current sheets. These are therefore the regions where heating occurs and our temperature events are found.

The density profile of the loop is practically unaltered (up to 5%) during the evolution. The fact that the maximum density within the temperature events has a distribution is mainly related to the fact that they extend down from the top of the loop, where the density increases owing to the stratification. The Joule heating events extending all the way to the bottom of the loop explains the appearance of a peak in the maximum density in their analysis. As a consequence of the small amount of energy involved in each event and the fact that the current sheets move across the loop during the current dissipation, along with the mentioned spread of conductive flux, very little thermal energy is conducting from the apex to the bottom of the loop. It follows that the evaporation from the chromosphere does not affect significantly the density of the loop.

As the numerical resolution is increased (and the physical dissipation decreased), we expect to find higher currents within the current sheets with smaller thickness and therefore smaller volumes. This is confirmed by the comparison of low- versus high-resolution simulations presented in this paper. We also expect an increase of the number of elementary events occurring during the same time. The energies involved should remain more or less the same because of the simultaneous increase of the current densities and decrease of the volumes. We will analyze Joule heating events in high-resolution RMHD simulations and compare them with the statistics presented here in an upcoming work.

It is important to emphasize that in our MHD framework all the dissipated energy goes into heating and the particle acceleration due to the well-known Fermi mechanisms (stochastic and systematic) is neglected. These mechanisms are at work at the very small scales the present paper indicates exist in the solar corona, contributing to the formation of a high-energy tail in the particle distribution (see Vlahos & Isliker 2019 and references therein).

The results of this paper indicate that it will still remain very challenging to observationally resolve elementary events because of their typical size and duration and the small amount of energy involved in each of them. At present only indirect observational tools will be able to verify the theory, as we have done in a previous paper through, for example, the use of emission measures. Moreover, a correct theoretical approach requires the use of 3D equations in order to properly take into account the perpendicular dynamics that are responsible for both the energy deposition and the continuous change of magnetic connectivity that makes the energy exchange between corona and chromosphere possible.

This material is based on work supported by the National Aeronautics and Space Administration under grant/contract/agreement No. NNH17AE96I issued through the Heliophysics Grand Challenge Research Program. R.B.D. was also supported by the Naval Research Laboratory 6.1 program. F.R. and M.V. acknowledge support of the NASA DRIVE SOLFer grant N.80NSSC20K0627. Computer simulations were performed on the LCP&FD Intel Core i7 cluster and the LCP&FD AMD Magny Cours cluster. We thank H. P. Warren for many helpful discussions and the referee for useful remarks.