Turbulent coronal heating and the distribution of nanoflares

We perform direct numerical simulations of an externally driven two-dimensional magnetohydrodynamic system over extended periods of time to simulate the dynamics of a transverse section of a solar coronal loop. A stationary and large-scale magnetic forcing was imposed, to model the photospheric motions at the magnetic loop footpoints. A turbulent stationary regime is reached, which corresponds to energy dissipation rates consistent with the heating requirements of coronal loops. The temporal behavior of quantities such as the energy dissipation rate show clear indications of intermittency, which are exclusively due to the strong nonlinearity of the system. We tentatively associate these impulsive events of magnetic energy dissipation to the so-called nanoflares. A statistical analysis of these events yields a power law distribution as a function of their energies with a negative slope of 1.5, which is consistent with those obtained for flare energy distributions reported from X-ray observations.


Introduction
Coronal loops in active regions are likely to be heated by Joule dissipation of highly structured electric currents. The formation of small scales in the spatial distribution of electric currents is necessary to enhance magnetic energy dissipation and therefore provide sufficient heating to the plasma confined in these loops. Various scenarios of how these fine scale structures might be generated have been proposed, such as the spontaneous formation of tangential discontinuities (Parker 1972, Parker 1983, the development of an energy cascade driven by random footpoint motions on a force-free configuration (van Ballegooijen 1986), or the direct energy cascade associated to a fully turbulent magnetohydrodynamic (MHD) regime (Heyvaerts & Priest 1992, Gómez & Ferro Fontán 1992. These rather different models share in common the dominant role of nonlinearities in generating fine spatial structure. In this paper we assume that the dynamics of a coronal loop driven by footpoint motions is described by the MHD equations. Since the kinetic (R) and magnetic (S) Reynolds numbers in coronal active regions are extremely large (R ∼ S ∼ 10 10−12 ), we expect footpoint motions to drive the loop into a strongly turbulent MHD regime.
Footpoint motions whose lengthscales are much smaller than the loop length cause the coronal plasma to move in planes perpendicular to the axial magnetic field, generating a small transverse magnetic field component. In §2 we model this coupling to simulate the driving action of footpoint motions on a generic transverse section of the loop. The numerical technique used for the integration of the two-dimensional MHD equations is described in §3. In §4 we report the energy dissipation rate that we obtain, and a statistical analysis of dissipation events is presented in §5. Finally, the relevant results of this paper are summarized in §6.

Forced two-dimensional magnetohydrodynamics
The dynamics of a coronal loop with a uniform magnetic field B = B 0 z, length L and transverse section (2πl) × (2πl), can be modeled by the RMHD equations (Strauss 1976): (1) where v A = B 0 / √ 4πρ is the Alfven speed, ν is the kinematic viscosity, η is the plasma resistivity, ψ is the stream function, a is the vector potential, w = −∇ 2 ψ is the fluid vorticity, j = −∇ 2 a is the electric current density and [u, v] = z · ∇u × ∇v. For given photospheric motions applied at the footpoints (plates z = 0 and z = L) horizontal velocity and magnetic field components develop in the interior of the loop, given by v = ∇ × (zψ) and b = ∇ × (za).
The RMHD equations can be regarded as describing a set of two-dimensional MHD systems stacked along the loop axis and interacting among themselves through the v A ∂ z terms. For simplicity, hereafter we study the evolution of a generic two-dimensional slice of a loop to which end we model the v A ∂ z terms as external forces (see Einaudi et al. 1996 for a similar approach). We assume the vector potential to be independent of z and the stream function to interpolate linearly between ψ(z = 0) = 0 and ψ(z = L) = Ψ, where Ψ(x, y, t) is the stream function for the photospheric velocity field. These assumptions yield (1)) and v A ∂ z j = 0 (in Eqn (2)) and correspond to an idealized scenario where the magnetic stress distributes uniformly throughout the loop. The 2D equations for a generic transverse slice of a loop become,

Numerical procedures
We performed numerical simulations of Eqs (3)-(4) on a square box of sides 2π, with periodic boundary conditions. The magnetic vector potential and the stream function are expanded in Fourier series. To be able to perform long-time integrations, we worked with a moderate resolution version of the code (96 × 96 grid points). The code is of the pseudo-spectral type, with 2/3 de-aliasing (Canuto et al. 1988). The temporal integration scheme is a fifth order predictor-corrector, to achieve almost exact energy balance over our extended time simulations.

Energy dissipation rate
To restore the dimensions to our numerical results, we used typical values for the solar corona: L ∼ 5 × 10 9 cm, l ∼ 10 8 cm, v p ∼ 10 5 cm/s, B 0 ∼ 100 G, n ∼ 5 × 10 9 cm −3 and The energy dissipation rate is also a strongly intermittent quantity as shown in Fig. 2. For turbulent systems at large Reynolds numbers, the dissipation rate in the stationary regime is expected to be independent of the Reynolds number R (Kolmogorov 1941). For the rather moderate Reynolds number simulations reported here, a weak (monotonically increasing) dependence of the dissipation rate with R still remains (see also Einaudi et al. 1996).
EDITOR: PLACE FIGURE 2 HERE.
The total dissipation rate for the typical numbers listed above is ǫ ≃ 1.6 × 10 24 erg.s −1 .
We can transform the heating rate into an energy influx from the photosphere, by simply dividing by twice the transverse area (because we have two boundaries), i.e. F = ǫ/(2(2πl) 2 ). In Eqn (5) This energy flux compares quite favorably with the heating requirements for active regions, which span the range F = 3 × 10 5 − 10 7 erg cm −2 s −1 (Withbroe & Noyes 1977).

Distribution of nanoflares
In this section, we associate the peaks of energy dissipation displayed in Fig. 2 to the so-called nanoflares (Parker 1988). We estimate the occurrence rate for these nano-events, i.e. the number of events per unit energy and time P (E) = dN/dE so that, is the total number of events per unit time and is the total heating rate (in erg s −1 ) contributed by all events in the energy range [E min ; E max ]. A simple inspection to the ǫ(t) time series shown in Fig. 2 indicates that these events are in a highly concentrated or piled-up regime, i.e. their event rate R multiplied by their typical duration is much larger than unity. At any given time, many events are going on simultaneously.
This piled-up scenario is a serious drawback against performing any kind of statistical analysis. To overcome this difficulty we define an event in the following fashion: (1) we set a threshold heating rate ǫ 0 on the time series displayed in Fig. 2, of the order of its time average, (2) events are excesses of dissipation which start when ǫ(t) surpasses ǫ 0 and finish when ǫ(t) returns below ǫ 0 . Once a particular threshold is set, we perform a statistical analysis of the events, keeping track of their peak values, durations and total energy content.
The implicit assumption behind our working definition, is that the small fraction of events that emerge over the threshold are statistically representative of the whole set. We do not make any attempt to prove this assertion, which therefore remains as a working hypothesis.
Among the interesting results of this statistical analysis, we find a significant correlation between the energy and duration of events like E ≃ τ 2 (Dmitruk & Gómez 1997), which is consistent with the correlation reported by Lee et al. 1993 from hard X-ray observations and by Lu et al. 1993 for an avalanche model of flare occurrences. Perhaps the most important result is that the occurrence rate as a function of energy displays a power law behavior in the energy range spanned from E min ≃ 5 × 10 24 erg to E max ≃ 10 26 erg.
We define the constant A in Eqn (8) so that the heating rate computed from Eqn (7) matches the total heating rate from our simulation (see Eqn (5)). Figure 3 shows our occurrence rate, displaying the power law behavior indicated in Eqn (8). For comparison, we also plotted the occurrence rate for transient brightenings derived by Shimizu 1995 (slope between 1.5-1.6) from Yohkoh soft X-ray observations and by Crosby et al. 1993 (slope 1.53) from SMM hard X-ray data.
EDITOR: PLACE FIGURE 3 HERE.
Also, we obtained the distribution of events as a function of peak fluxes which is a power law with slope 1.7 ± 0.3. This slope is consistent with the one derived by Crosby et al. 1993 (1.68) from X-ray events and somewhat flatter than those reported by Hudson 1991 (1.8).
Note that the slopes derived from X-ray observations assume that the luminosity in X-rays is proportional to the dissipated energy. Porter et al. 1995 criticize this assumption after comparing the UV and X-ray emission for several microflares and find that slopes derived from X-rays become slightly steeper when corrected for this effect.
The remarkable correspondence between the rates plotted in Fig. 3, is indicative of the presence of a common physical process behind the dissipation events ranging from 10 24 erg or less, to 10 33 erg for the largest flares. Since in all cases the index of the power law remains smaller than two, Eqn (7) implies that the contribution to energy dissipation in a given energy range ([E min ; E max ]) is dominated by the most energetic events (i.e. E ≃ E max ).
According to this result, the relatively infrequent large energy events contribute more to the heating rate than the much more numerous small energy events. This assertion might only change if indications of a turn-up of the slope (to an index larger than two) are found at low energies (Hudson 1991).

Discussion and conclusions
In the present paper we simulate the dynamics of a transverse section of a solar coronal loop through an externally driven two-dimensional MHD code. The relevant results of this study are summarized as follows: (1) for an external forcing which is narrowband in wavenumber, we find that the system becomes strongly turbulent and after about twenty photospheric turnover times a stationary turbulent regime is reached, (2) the energy dissipation rate obtained for typical footpoint velocities is consistent with the power necessary to heat active region loops (F ≃ 2 × 10 6 erg cm −2 s −1 ), (3) the energy dissipation rate displays a highly intermittent behavior, which is a ubiquitous characteristic of turbulent systems, (4) we associate the elementary events that compose this intermittent heating rate to the so-called nanoflare events described by Parker 1988. A statistical analysis of the events performed on a long term numerical simulation (about 200 turnover times) shows a power law event rate, going like dn/dE ∼ E −1.5 , which is remarkably consistent with the statistics of flare occurrence derived from observations (Hudson 1991, Crosby et al. 1993, Lee et al. 1993, Shimizu 1995. We acknowledge financial support by the University of Buenos Aires (grant EX247) and by Fundación Antorchas.  Crosby et al. 1993