# Scaling law for the heating of solar coronal loops

## 1 Introduction

Solar coronal loops are likely to be heated by ohmic dissipation of field-aligned electric currents. These currents are driven by photospheric motions which twist and shear the magnetic fieldlines at the loop footpoints. The generation of small scales in the spatial distribution of these currents has been observed and studied in recent simulations (Mikić et al. (1989), Longcope & Sudan (1994)). Small scale currents are required by most coronal heating theories, since they produce a major enhancement of the energy dissipation rate. The development of magnetohydrodynamic (MHD) turbulence in coronal loops has recently been proposed (Gómez & Ferro Fontán (1992), Heyvaerts & Priest (1992), Einaudi et al. (1996), Hendrix & Van Hoven (1996), Dmitruk & Gómez (1997) ) as a likely mechanism to provide such enhancement, since the energy being pumped by footpoint motions is naturally transferred to small scale structures by the associated energy cascade.

In the present paper we numerically test this scenario, describing the dynamics of a coronal loop through the reduced MHD (RMHD) approximation. In §2 we briefly describe the RMHD approximation and perform a dimensional analysis of the equations, which leads to an interesting prediction for the functional dependence of the heating rate with the physical parameters of the problem. The details of the numerical simulations are summarized in §3 and an overview of the results is given in §4. A quantitative scaling law for the heating rate is derived in §5, and the relevant results of this paper are listed in §6.

## 2 Dimensional analysis of the RMHD equations

When a coronal loop with an initially uniform magnetic field $`𝐁=B_0\widehat{𝐳}`$ is driven by photospheric motions at its footpoints, a rather complex dynamical evolution sets in, which is described by velocity and magnetic fields. For incompressible and elongated loops, i.e. loops of length $`L`$ and transverse section $`(2\pi l_0)\times (2\pi l_0)`$ such that $`l_0L`$, the velocity ($`𝐯=𝐯(x,y,z,t)`$) and magnetic $`𝐁=B_0\widehat{𝐳}+𝐛(x,y,z,t)`$) fields can be written in terms of a stream function $`\psi =\psi (x,y,z,t)`$ (i.e. $`𝐯=_{}\times (\psi \widehat{𝐳})`$) and a vector potential $`a=a(x,y,z,t)`$ (i.e. $`𝐛=_{}\times (a\widehat{𝐳})`$), respectively. The evolution of the fields $`a`$ and $`\psi `$ is determined by the RMHD equations (Strauss (1976)):

$$_ta=v_A_z\psi +[\psi ,a]+\eta _{}^2a$$
(1)   
$$_tw=v_A_zj+[\psi ,w][a,j]+\nu _{}^2w$$
(2)   
where $`v_A=B_0/\sqrt{4\pi \rho }`$ is the Alfvén speed, $`\nu `$ is the kinematic viscosity and $`\eta `$ is the plasma resistivity. The quantities $`w=_{}^2\psi `$ and $`j=_{}^2a`$ are the $`z`$-components of vorticity and electric current density, respectively. The non-linear terms are standard Poisson brackets, i.e. $`[u,v]=_xu_yv_yu_xv`$. We assume periodicity for the lateral boundary conditions, and specify the velocity fields at the photospheric boundaries,

$$\psi (z=0)=0,\psi (z=L)=\mathrm{\Psi }(x,y)$$
(3)   
where $`\mathrm{\Psi }(x,y)`$ is the stream function which describes stationary and incompressible footpoint motions. The strength of this external velocity field is therefore proportional to a typical photospheric velocity $`V_p`$.

To transform Eqs (1)-(2) into their dimensionless form, we choose $`l_0`$ and $`L`$ as the units for transverse and longitudinal distances. Since the dimensions of all physical quantities involved in these equations can be expressed as combinations of length and time, let us choose $`t_AL/v_A`$ as the time unit. The dimensionless RMHD equations are:

$$_ta=_z\psi +[\psi ,a]+\frac{1}{S}_{}^2a$$
(4)   
$$_tw=_zj+[\psi ,w][a,j]+\frac{1}{R}_{}^2w$$
(5)   
where $`S=\frac{l_0^2}{\eta t_A}`$ and $`R=\frac{l_0^2}{\nu t_A}`$ are respectively the magnetic and kinetic Reynolds numbers. Hereafter, we will consider the case $`S=R`$, and thus the (common) Reynolds number will be the only dimensionless parameter explicitly present in Eqs (4)-(5) . However, the boundary condition (Eqn (3)) brings an extra dimensionless parameter into play, which is the photospheric velocity $`V_p`$ divided by our velocity units, i.e. $`l_0/t_A`$. Since $`V_p=l_0/t_p`$ ($`t_p`$: photospheric turnover time), this velocity ratio is equivalent to the ratio between the Alfven time $`t_A`$ and the photospheric turnover time $`t_p`$. Therefore, from purely dimensional considerations, we can derive the following important result: for any physical quantity, its dimensionless version $`Q`$ should be an arbitrary function of the only two dimensionless parameters of the problem, i.e.

$$Q=(Q_1,Q_2),Q_1=\frac{t_A}{t_p},Q_2=S$$
(6)   
For instance, let us consider the important case of the heating rate per unit mass, i.e. $`ϵ/\rho `$. Its dimensionless version is,

$$Q=\frac{ϵ}{\rho }\frac{t_A^3}{l_0^2}=(\frac{t_A}{t_p},S)$$
(7)   
One of Kolmogorov’s hypothesis in his theory for stationary turbulent regimes at very large Reynolds numbers (Kolmogorov (1941)), states that the dissipation rate is independent of the Reynolds number (see also Frisch (1996)). In the next section we show that externally driven coronal loops eventually reach a stationary turbulent regime. Therefore, hereafter we assume that the dependence of the dissipation rate $`ϵ`$ with the Reynolds number $`S`$ in Eqn (7) can be neglected. For moderated values of the Reynolds numbers as the ones considered here, we observed only a mild dependence of the dissipation rate with S (see also Hendrix et al. (1996)), in accordance with Kolmogorov’s assumption. Therefore,

$$ϵ=\frac{\rho l_0^2}{t_A^3}(\frac{t_A}{t_p})$$
(8)   
We performed a sequence of numerical simulations for different values of the ratio $`\frac{t_A}{t_p}`$ to determine the function $``$.

## 3 Numerical simulations

For the numerical simulations of Eqs (4)-(5) , $`\psi `$ and $`a`$ are expanded in Fourier modes in each $`(x,y)`$ plane ($`0x,y2\pi `$ and $`0z1`$). The equations for the coefficients $`\psi _𝐤(z,t)`$ and $`a_𝐤(z,t)`$ are time-evolved using a semi-implicit scheme. Linear terms are treated in a fully implicit fashion, while nonlinear terms are evolved using a predictor-corrector scheme. Also, nonlinear terms are evaluated following a $`2/3`$ fully dealiased (see Canuto et al. (1988)) pseudo-spectral technique (see also Dmitruk, Gómez & DeLuca (1998)). To compute $`z`$-derivatives we use a standard method of finite differences in a staggered regular grid (see for instance Strauss (1976)).

We model the photospheric boundary motion in Eqn (3) as $`\mathrm{\Psi }_𝐤=\mathrm{\Psi }_0=`$ constant inside the ring $`3<kl_0<4`$, and $`\mathrm{\Psi }_𝐤=0`$ elsewhere. This choice is intended to simulate a stationary and isotropic pattern of photospheric granular motions of diameters between $`2\pi l_0/4`$ and $`2\pi l_0/3`$. We chose a narrowband and non-random forcing to make sure that the broadband energy spectra and the signatures of intermittency that we obtained (see below) are exclusively determined by the nonlinear nature of the MHD equations. The normalization factor $`\mathrm{\Psi }_0`$ is proportional to $`t_A/t_p`$, as mentioned above.

## 4 Development of MHD turbulence

We performed a sequence of numerical simulations of Eqs (4)-(5) with $`192\times 192\times 32`$ gridpoints, $`S=1500`$ and different values of the ratio $`t_A/t_p`$ in the range $`[0.01,0.15]`$. The typical behavior of the heating rate as a function of time, is shown in Figure 1 for the particular case of $`t_A/t_p=0.064`$. After an initial transient, the heating rate is seen to approach a stationary level. This level is fully consistent with the heating requirements for coronal active regions of $`10^7ergcm^2seg^1`$ (Withbroe & Noyes (1977)). The intermittent behavior of this time series, which is typical of turbulent systems, is also ubiquitous in all our simulations. Once the stationary regime is reached in each of the simulations, we determine the mean value and rms of the departure of the series from the mean.

A stronger indication of the presence of a turbulent regime, is the development of a broadband energy power spectrum, which behaves like $`E_kk^{3/2}`$ for both two and three dimensional MHD turbulence (Kraichnan spectrum). To compute the energy power spectrum, we performed a single simulation with better spatial resolution ($`384\times 384\times 32`$), to allow for a more extended inertial range. The parameter for this simulation is $`t_A/t_p=0.064`$ and $`S=5000`$. The spectrum shown in Fig. 2 was taken at $`t=5t_p`$, i.e. well in the stationary regime. This higher spatial resolution is still insufficient for the formation of a broad inertial range to determine the slope of these spectra with confidence. However, the observed slope is consistent with a Kraichnan spectrum, i.e. $`E_kk^{3/2}`$. Also, note that the kinetic and magnetic power spectra reach equipartition at large wavenumbers, even though the total kinetic energy is much smaller than the total magnetic energy.

Another important consequence of the energy power spectrum displayed in Fig. 2, is the fact that smaller scales dissipate more energy than the larger scales, since $`ϵ_kk^2E_kk^{1/2}`$. Three dimensional simulations of loops driven by random footpoint motions have shown that energy dissipation preferentially occurs in current sheets which form exponentially fast (Mikić et al. (1989), Longcope & Sudan (1994)). Within the framework of MHD turbulence, two dimensional simulations also show the ubiquitous presence of current sheets evolving in a rather dynamic fashion (Matthaeus & Lamkin (1986), Biskamp & Welter (1989), Dmitruk, Gómez & DeLuca (1998)). The presence of elongated current sheets as the most common dissipative structures is also apparent in these simulations (also Hendrix & Van Hoven (1996)). Figure 3 shows the spatial distribution of the electric current density $`j(x,y,z,t)`$ for our $`384\times 384\times 32`$ simulation ($`t_A/t_p=0.064`$ and $`S=5000`$). Intense currents going upward (downward) are shown in white (black) at various planes $`z=constant`$, and for $`t=5t_p`$. Although the reconnection in these turbulent regimes is expected to be fast (Priest & Forbes (1992), Hendrix & Van Hoven (1996), Milano et al. (1999)), the Reynolds numbers attained by current simulations cannot provide a definite answer to this question.

## 5 Determination of the heating rate scaling exponent

One of the main goals of this paper is to derive the functional dependence of the heating rate with the ratio between the two relevant timescales (see Eqn (8)). To this end, we performed various simulations for different values of the parameter $`t_A/t_p`$. Figure 4 shows that this functional dependence can be adequately fit by a power law, i.e.

$$ϵ=\frac{\rho l_0^2}{t_A^3}(\frac{t_A}{t_p})^s,s=1.51\pm 0.04$$
(9)   
Each diamond corresponds to the mean value of the heating rate for a particular simulation in its stationary regime. The error bars correspond to the rms value of the departure from the mean. The slope in this plot corresponds to the parameter $`s`$ quoted in Eqn (9), which was derived following a minimum squares procedure.

This result is fully consistent with the prediction arising from a two dimensional MHD model (Dmitruk & Gómez (1997)), $`s_{2D}=\frac{3}{2}`$, which is assumed to simulate the dynamics of a generic transverse slice of a loop. This coincidence between both scalings suggests that two dimensional simulations provide an adequate description of the dynamics of coronal loops (see also Hendrix & Van Hoven (1996)), provided that the external forcing in Eqs (1)-(2) is given by $`v_A_z\psi \frac{l_0V_p}{t_A}`$ and $`v_A_zj0`$ (Einaudi et al. (1996), Dmitruk & Gómez (1997)).

## 6 Discussion and conclusions

In the present paper we report preliminary results from a series of numerical simulations of the RMHD equations driven at their footpoints by stationary granule-size motions. The important results are summarized as follows:

* After a few photospheric turnover times, the system relaxes to a stationary turbulent regime, with dissipation rates consistent with the heating requirements of coronal active regions ($`10^610^7ergcm^2s^1`$).
* A dimensional analysis of the RMHD equations shows that the dimensionless heating rate does only depend on the ratio $`t_A/t_p`$ and the Reynolds number. However, for turbulent systems at very large Reynolds numbers, the heating rate is likely to be independent of the Reynolds number. From a series of numerical simulations, we find that the heating rate scales with $`t_A/t_p`$ as $`ϵ\frac{\rho l_0^2}{t_A^3}(\frac{t_A}{t_p})^{3/2}`$, i.e.

$$ϵ\frac{\rho ^{\frac{1}{4}}(B_0V_p)^{\frac{3}{2}}}{L}(\frac{l_0}{L})^{\frac{1}{2}}$$
(10)   
which might be useful for an observational test of this theoretical framework of coronal heating.
* The dissipative structures observed in our simulations are current sheets elongated along the axis of the loop. The spatial and temporal distribution of these structures is rather intermittent, as expected for turbulent regimes. We associate these dissipation events taking place inside coronal loops to the nanoflares envisioned by Parker (1988) , as a likely scenario for coronal heating.