THE SPATIAL AND TEMPORAL DEPENDENCE OF CORONAL HEATING BY ALFVÉN WAVE TURBULENCE

The solar corona is two orders of magnitude hotter than the underlying photosphere. This has been puzzling scientists over many decades as they try to identify the mechanism(s) behind these extreme conditions. Many different ideas on coronal heating have been proposed (see reviews by Aschwanden 2005; Klimchuk 2006), but solar observations have not yet advanced to the point where they allow us to clearly identify the physical mechanisms responsible for coronal heating. One plausible idea is that the corona is heated by dissipation of Alfvén waves (Hollweg 1981, 1986; Heyvaerts & Priest 1983; Kudoh & Shibata 1999; Moriyasu et al. 2004; Matsumoto & Shibata 2010; Antolin & Shibata 2010) or by turbulence (Rappazzo et al. 2008; Dmitruk & Gomez 1997; Cranmer et al. 2007). Alfvén waves have indeed been observed at various heights in the solar atmosphere (Ulrich 1996; Fujimura & Tsuneta 2009; De Pontieu et al. 2007b; Tomczyk & McIntosh 2009; McIntosh et al. 2011), but the role of such waves in coronal heating has not yet been clearly demonstrated.

Recently, the present authors developed a three-dimensional (3D) magnetohydrodynamics (MHD) model describing the propagation and dissipation of Alfvén waves in a coronal loop, including the lower atmospheres at the two ends of the loop (van Ballegooijen et al. 2011, hereafter Paper I). According to this model, the waves are generated by the interactions of granule-scale convective flows within kilogauss flux tubes in the photosphere. The flows produce non-axisymmetric kink-type waves on a transverse length scale less than 100 km in the photosphere, below the resolution of present-day telescopes. The waves propagate upward along the tubes and strongly reflect in the chromosphere and transition region (TR), producing counter-propagating waves that are subject to strong nonlinear wave–wave interactions (e.g., Shebalin et al. 1983; Higdon 1984; Oughton & Matthaeus 1995; Goldreich & Sridhar 1995, 1997; Bhattacharjee & Ng 2001; Cho et al. 2002; Oughton et al. 2004). These nonlinear interactions result in turbulent transfer of wave energy to small spatial scales and heating of the chromospheric plasma. A fraction of the wave energy is transmitted through the TR into the corona, and produces turbulence and heating there. It was found that this model can quantitatively explain both the chromospheric and coronal heating rates observed in active regions, and is consistent with observational constraints on the degree of braiding of the coronal loops.

A key feature of the Alfvén wave turbulence model is that it naturally predicts the spatial variation of heating rate Q(s) with position s along a coronal loop, so the heating profile is not arbitrarily prescribed. The thermal stability of a coronal loop depends strongly on the heating profile (e.g., Antiochos & Klimchuk 1991; Müller et al. 2003; Testa et al. 2005; Karpen et al. 2006; Mok et al. 2008; Klimchuk et al. 2010). If most of the heat is deposited near the loop footpoints, much of the energy is conducted downward into the upper chromosphere, producing strong chromospheric evaporation and large coronal densities. The radiative losses at the loop top may then become too large to be balanced by local heating. This results in the formation of cool condensations that move along the loop under the influence of gravity and pressure gradients. Modeling predicts that the loop undergoes periodic convulsions as it searches for a nonexistent equilibrium state. Klimchuk et al. (2010) examine the possibility that the observed warm (1–1.5 MK) loops in active regions can be explained by such thermal nonequilibrium.

Asgari-Targhi & van Ballegooijen (2012, hereafter Paper II) applied the Alfvén wave turbulence model to an active region observed with the Solar Dynamics Observatory (SDO). The dynamics of Alfvén waves in selected coronal loops was computed by the model of Paper I, based on magnetic field strengths taken from a 3D magnetic model of the active region. Wave heating rates were derived from these simulations, and average coronal temperatures and densities were computed self-consistently with the modeled heating rates. It was found that the loops in the core of the active region were thermally stable with peak temperatures of 2–3 MK. The expected temperature fluctuations of such loops are quite small, ΔT ∼ 0.1 MK, implying that the loops never cool off into the 1–1.5 MK range. However, other loops at the periphery of the active region were predicted to be thermally unstable, and presumably exist in a state of thermal nonequilibrium with large temperature and density fluctuations. The thermal evolution of such unstable loops was not addressed in Paper II.

To further test the Alfvén wave turbulence model, it is necessary to construct detailed numerical models of wave-heated active regions and then compare such models with solar observations. Such models must describe not only the magnetic field but also the temperature and density of the plasma. Both thermally stable and unstable coronal loops must be considered. To follow the evolution of the unstable loops, the time evolution of temperature and density must be simulated with high spatial resolution as neighboring loops evolve almost independently. Such simulations will be very challenging and are beyond the scope of the present paper. However, the simulations will require detailed knowledge of how the local wave heating rate depends on coronal field strength, loop length, and coronal density. It is also important to know how the heating rate varies with position along a coronal loop because this is an important factor in determining whether a loop is stable or unstable.

The main objective of the present paper is to derive an expression for the wave heating rate that can be used in future 3D MHD simulations of active region structure and evolution. A second objective is to study the internal motions of a wave-heated loop and discuss the effects of such motions on plasma temperature. These results are important for understanding the observable consequences of the Alfvén wave turbulence model.

Our previous work indicates that the wave heating rate Q depends strongly on magnetic structure, and to lesser extent on coronal density (see papers I and II). Therefore, it is important that the Alfvén wave simulations are done using a realistic model for the background magnetic field. In the present paper we use data from the Helioseismic and Magnetic Imager (HMI) on board SDO to construct a 3D potential field model of an observed active region, and we trace 22 field lines through this model. For each field line we construct a detailed numerical model of the wave dynamics, which yields the time-averaged heating rate Q(s) as a function of position along that field line. Finally, we combine the results from different field lines into a single expression for the heating rate and its dependence on magnetic field strength, density, and the loop length.

The paper is organized as follows. Section 2 presents the SDO observations of an active region, and the construction of the 3D magnetic model for this region. Section 3 describes the Alfvén wave turbulence modeling for individual flux tubes. Section 4 describes the simulation results. Section 5 presents the combined formula for the wave heating rate. Section 6 discusses the internal wave dynamics. The modeling results are further discussed in Section 7.

The active region used in the present study is NOAA 11428, which crossed the central meridian on 2012 March 7 and has a relatively simple, bipolar magnetic structure. We use SDO data taken around 18:30 UT on that day. Figure 1 shows the magnetic flux distribution in the photosphere (red and green contours) superposed on the AIA 171 Å image from 18:30:01 UT. The flux distribution is derived from HMI line-of-sight magnetograms taken over a 15 minute period. The method for extracting the flux distribution from a magnetogram is described by Su et al. (2009). Sunspots are present in both the leader and follower parts of the region. The AIA image shows a set of coronal loops and fans emanating from the two spots and surrounding fluxes.

Zoom In Zoom Out Reset image size

This particular active region does not exhibit much "sigmoidal twist," and therefore it is modeled based on Potential Field Source Surface (PFSS) model. We construct a PFSS model of the observed region, based on the flux distribution derived from the HMI magnetogram and the SOLIS synoptic map for Carrington Rotation 2121. The PFSS model was computed using the Coronal Modeling System software (Su et al. 2011). The model includes both a high resolution grid covering the target active region (grid spacing 10⁻³R_☉) and a low resolution grid covering the entire Sun (1° spacing). The source surface is assumed to be located at radial distance r_s = 2.05R_☉ from Sun center. We trace 22 randomly selected field lines through the PFSS model. Figure 1 shows the selected field lines superposed on the AIA 171 Å image. For each field line we determine the magnetic field strength B(s) and the height z(s) above the photosphere as functions of position s along the field line, which are needed for the wave modeling described in the next section.

The Alfvén waves are simulated using the 3D MHD model developed in Papers I and II. In this model a thin flux tube surrounding a traced magnetic field line is considered. The tube is assumed to have a circular cross-section with radius R(s) that varies with position s along the tube, which extends from the photosphere at one end of the coronal loop to the photosphere at the other end. Even though our PFSS model uses high spatial resolution, the field strengths in the lower atmosphere predicted from the magnetogram are not as high as those known to exist on the Sun (Stenflo 1973). Therefore, we artificially increase the field strengths at the ends of the flux tube to produce kilogauss field strengths in the photosphere: B_Phot = 1400 [G]. This is done on each computed field line, but only at heights below the TR. The corrected field strength is

$$\begin{equation} B(s) = \sqrt{B_{\rm Phot} ^2 \tilde{f}(s) + B_0^2 (s) [1 - \tilde{f}(s) ] }, \end{equation} \tag{ 1 }$$

where $\tilde{f}(s)$ is a weighting factor derived from the gas pressure:

$$\begin{equation} \tilde{f}(s) = \frac{p(s) - p_{\rm TR}}{p_{\rm Phot} - p_{\rm TR}}, \end{equation} \tag{ 2 }$$

and B₀(s) is the field strength as derived from the PFSS model. The height dependent gas pressure p(s) decreases monotonically from its photospheric value (p_Phot) to its value at the TR (p_TR). The background field strength B(s) and plasma density ρ(s) are assumed to be constant over the tube's cross-section, so the Alfvén speed v_A(s) ($\equiv B / \sqrt{4 \pi \rho }$) is also constant over the cross-section. The plasma temperature T(s) is fixed in time, and mass flows along the flux tube are neglected. Therefore, in its present form the model is not suitable for simulating thermally unstable loops, which have time-variable temperature, density, and Alfvén speed.

At the two ends of the flux tube the field lines are subjected to random footpoint motions that simulate the effects of convective flows on kilogauss flux elements in the photosphere. The imposed footpoint motions are assumed to have an rms velocity of 1.5 km s⁻¹ and correlation time of 60 s, similar to the observed random motions of G-band bright points in the photosphere (e.g., Chitta et al. 2012). The footpoint motions produce Alfvén waves that travel upward inside the flux tube and reflect at various heights due to spatial variations in Alfvén speed v_A(s). Such reflections result in counter-propagating Alfvén waves that interact nonlinearly with each other and produce wave turbulence (Kraichnan 1965; Shebalin et al. 1983; Oughton & Matthaeus 1995; Cho et al. 2002; Cranmer 2010).

The wave dynamics are simulated numerically using the "reduced MHD" approximation (e.g., Strauss 1976; Paper I), which means that slow and fast MHD modes are filtered out and only the Alfvén mode is retained. Furthermore, only the magnetic- and velocity fluctuations of the waves are simulated, not their effects on temperature and density. The magnetic fluctuations B₁ are assumed to be small compared to the background field B, and can be approximated as B₁ = ∇_⊥h × B, where h(r, t) is the magnetic flux function and t is the time. Similarly, the velocity fluctuations are approximated as ${\bf v}_1 = \nabla _\perp f \times \hat{\bf B}$, where f(r, t) is the velocity stream function and $\hat{\bf B} (x,y,s)$ is the unit vector along the background field. Flows along the background field are neglected. In the present version of the reduced MHD model the functions f(r, t) and h(r, t) satisfy the following coupled equations:

$$\begin{equation} \frac{\partial \omega }{\partial t} + \hat{\bf B} \cdot (\nabla _\perp \omega \times \nabla _\perp f) = v_A^2 [ \hat{\bf B} \cdot \nabla \alpha + \hat{\bf B} \cdot (\nabla _\perp \alpha \times \nabla _\perp h) ] + D_v, \end{equation} \tag{ 3 }$$

$$\begin{equation} \frac{\partial h}{\partial t} = \hat{\bf B} \cdot \nabla f + \frac{f}{H_B} + \hat{\bf B} \cdot (\nabla _\perp f \times \nabla _\perp h) + D_m, \end{equation} \tag{ 4 }$$

where $\omega \equiv - \nabla _\perp ^2 f$ is the parallel component of vorticity, $\alpha \equiv -\nabla _\perp ^2 h$ is the magnetic torsion parameter, and H_B(s) ≡ B/(dB/ds) is a magnetic scale length. The terms D_v and D_m describe the effects of viscosity and resistivity on the high wavenumber modes. Derivations of the above equations and descriptions of their numerical implementation are given in Paper I.

The spatial dependence of the waves on the transverse coordinates x and y is described in terms of the eigenfunction of the $\nabla _\perp ^2$ operator on the domain x² + y² < R². These eigenfunctions have perpendicular wavenumbers k_⊥ = a_k/R(s), where the a_k are given by the zeros of Bessel functions, and k is the mode index (k = 1, ..., N). In the present work we use a_max = 30, which corresponds to a spatial resolution of about one-tenth of the tube radius and requires N = 209 modes. Also, the viscous and resistive damping rates are assumed to be equal, and vary with the fourth power of the perpendicular wavenumber ("hyperdiffusion"):

$$\begin{equation} \nu _k (s,t) = \nu _{\rm max} (s,t) \left(\frac{a_k}{a_{\rm max}} \right)^4. \end{equation} \tag{ 5 }$$

The maximum damping rate is given by $\nu _{\rm max} = 70 \overline{v_{\rm rms}}(s,t) / R(s)$, where v_rms(s, t) is the rms velocity of the waves, and the bar represent a running time average over a time interval of 300 s.

In the TRs at the two ends of the coronal loop the temperature and density vary rapidly with position, causing strong variation in Alfvén speed and therefore strong wave reflection. To avoid having to use very small time steps, we treat the TRs as discontinuities where the waves can reflect. In the corona the temperature is assumed to be a function of height only:

$$\begin{equation} T(s) = T_{\rm max} \left[ 0.01 + 0.99 \left(\frac{z(s) - z_{\rm TR}}{z_{\rm max} - z_{\rm TR}} \right) \right]^{2/7}, \end{equation} \tag{ 6 }$$

where z(s) is the height as derived from the PFSS model, z_max is the maximum height, and z_TR is the height of the two TRs. The peak temperature T_max is estimated from the coronal loop length L and coronal pressure p, using the RTV scaling law (Rosner et al. 1978). The coronal pressure p is the only free parameter of the model. The density ρ(s) is computed from the temperature by solving the equation of hydrostatic equilibrium along the field line.

The way we have set up the model is such that the temperature and density are fixed in time, so the radiative and conductive losses are also fixed. However, these radiative and conductive loss rates are not necessarily equal to the wave heating rate Q(s), which is determined by averaging Q(s, t) over the duration of the simulation. Therefore, the loops simulated here are generally not in thermal equilibrium. When the heating rate is higher (lower) than the radiative and conductive losses, the coronal temperature is expected to increase (decrease) with time, but this variation with time is not taken into account in our simulations.

After Equations (3) and (4) have been solved, the volumetric heating rate can be computed as a sum of viscous- and resistive contributions: Q = Q_kin + Q_mag. The averages of these quantities over the loop cross-section can be written as sums over perpendicular modes:

$$\begin{eqnarray} \hspace{-9pt}Q_{\rm kin} (s,t) & = & \frac{\rho }{R^2} \sum _{k=1}^N \nu _k a_k^2 f_k^2, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} Q_{\rm mag} (s,t) & = & \frac{B^2}{4 \pi R^2} \sum _{k=1}^N \nu _k a_k^2 h_k^2, \end{eqnarray} \tag{ 8 }$$

where f_k(s, t) and h_k(s, t) are mode amplitudes for the velocity stream function and magnetic flux function, respectively (see Appendix C in Paper I). In this paper we also compute the local heating rate Q(x, y, s, t) at any point within the loop. This quantity is the sum of the following non-negative expressions:

$$\begin{eqnarray} Q_{\rm kin} (\xi,\varphi,s,t) & = & \frac{\rho }{R^2} \nu _{\rm max} \left[ \sum _{k=1}^N \left(\frac{a_k}{a_{\rm max}} \right)^2 a_k f_k F_k (\xi,\varphi) \right]^2, \nonumber \\ \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} Q_{\rm mag} (\xi,\varphi,s,t) & = & \frac{B^2}{4 \pi R^2} \nu _{\rm max} \left[ \sum _{k=1}^N \left(\frac{a_k}{a_{\rm max}} \right)^2 a_k h_k F_k (\xi,\varphi) \right]^2, \nonumber \\ \end{eqnarray} \tag{ 10 }$$

where ξ ≡ r/R is a dimensionless radial coordinate, φ is the azimuth angle, and F_k(ξ, φ) is a set of orthogonal basis functions (see Appendix B in Paper I). Integrating these expressions over the cross-section and using the orthogonality of the eigenmodes, Equations (7) and (8) are recovered.

The above turbulence modeling was applied to 22 field lines traced through the PFSS model for NOAA region 11428, where each field line represents a particular coronal loop. The PFSS model predicts that most loops expand significantly with height; the loop expansion factor is defined by Γ ≡ B_TR/B_min, where B_TR is the magnetic field strength at the TR (average of the two TRs) and B_min is the minimum field strength along the loop. For each loop we construct several models for the background atmosphere with different values of the coronal pressure p, leading to a total of 77 different models shown in Table 1. The models use a broad range of plasma pressures because the results will be used to construct a heating rate formula (see Section 5) to be used in future dynamical loop modeling. We anticipate that for thermally unstable loops the temperature and density will vary strongly with time, and our formula should capture the effects of density on the wave heating rate. For thermally stable loops the highest value of pressure was chosen to lie above the value for thermal equilibrium (i.e., wave heating balanced by radiative and conductive losses).

For each model we simulate the 3D dynamics of the Alfvén waves in the neighborhood of the selected field line. The waves are described with modest spatial resolution in the direction perpendicular to the mean magnetic field (ℓ_⊥ ∼ 100 km in the corona). The waves are launched from the photosphere, and are simulated for a period of 3000 s. The time resolution is typically about 0.1 s.

To illustrate the simulation results, Figure 2 shows various parameters as function of position along the flux tube for model f19r2 (see Table 1). The model name f19r2 contains the field line number and the run number refers to different values for coronal pressure. Positions are expressed in terms of the Alfvén travel time τ along the loop:

$$\begin{equation} \tau (s) \equiv \int _{s_0}^s \frac{ds^\prime }{v_A (s^\prime)}, \end{equation} \tag{ 11 }$$

where s₀ denotes the position of the positive polarity footpoint. Figure 2(a) shows the position s (solid curve) and height z (dashed curve) plotted as functions of τ. Figure 2(b) shows the temperature T(s), which varies from about 6000 K in the photosphere and 8000 K in the chromosphere to about 1 MK in the corona. The magnetic field strength B(s) varies from 1400 G at the base of the photosphere to 8.7 G in the corona (see Figure 2(c)). The values of τ on the lower axes of these plots indicate that the time for an Alfvén wave to propagate from the base of the photosphere to the base of the corona is about 145 s at the positive polarity end of the loop (left side of each figure) and 136 s at the negative polarity end (right side). This is comparable to the time of 173 s for the wave to travel the much larger distance through the corona (L = 178 Mm). Therefore, to understand Alfvén wave dynamics in coronal loops it is important to include the lower atmospheres at the two ends of a loop. Figure 2(d) shows that the Alfvén speed varies from about 15 km s⁻¹ in the photosphere to more than 1000 km s⁻¹ in the low corona. Therefore, waves traveling to the corona suffer strong wave reflection (Hollweg 1981). Nevertheless, we find that the waves in the corona build up to significant amplitudes and have strong nonlinear interactions.

Zoom In Zoom Out Reset image size

Figure 2(e) shows the time-averaged, rms velocity of the waves, v_rms(s), as predicted by the reduced MHD model. The velocity is averaged over the cross-section of the flux tube and over the period 200 < t < 3000 s. Note that v_rms increases from about 1.5 km s⁻¹ in the photosphere (the imposed footpoint velocity) to about 35 km s⁻¹ in the corona. This amplification is mainly due to the stratification of density, which drops by six orders of magnitude over these layers. The great advantage of the reduced MHD model is that it can easily deal with such large density differences. The time-averaged heating rate Q(s) is shown in Figure 2(f), together with the magnetic and kinetic contributions to the heating. Note that the heating rate in the coronal part of the loop is about 10⁻³erg cm⁻³ s⁻¹, much smaller than the heating rates in the lower atmospheres at the two ends of the loop. In this paper we focus on the heating in the coronal part of the loop. The coronal heating rate can be fit to a power law:

$$\begin{equation} Q(s) = Q_{\rm TR} \left[ \frac{B(s)}{B_{\rm TR}} \right]^m, \end{equation} \tag{ 12 }$$

where Q_TR is the average heating rate at the two TRs. The exponent m and heating rate Q_TR are determined by the fit.

Figure 3 shows the Alfvén wave energy flux F_A(s) for model f19r2, as computed with Equation (C12) of Paper I. The non-radiative energy flux entering the corona is F_A ∼ 10⁷ erg cm⁻² s⁻¹ consistent with Withbroe & Noyes (1977). This energy enters a coronal loop at both ends and is dissipated within the loop. For most models only 5–7 % of the available energy is dissipated in the corona; the remainder is dissipated in the photosphere and chromosphere at the two ends of the coronal loop.

Zoom In Zoom Out Reset image size

Results similar to those shown in Figure 2 are obtained for all 77 models. This modeling provides information on the turbulent heating produced by the waves along the 22 selected field lines for different values of coronal pressure p. The model parameters are listed in Table 1. The parameters include the coronal loop length L, transition region height z_TR, coronal loop expansion factor Γ, minimum field strength B_min in the corona, average field strength B_TR at the two TRs, and coronal pressure p used in setting up the background atmosphere, τ_AC is the coronal Alfvén travel time between the two TRs, τ_AC = τ(s_TR2) − τ(s_TR1) (see Equation (11)), exponent m from Equation (12), the average heating rate Q_TR, and the heating ratio η_Q (see Equation (13)). In most cases the exponent m > 0.5, so the heating rate Q(s) decreases significantly from the legs of the coronal loop where the magnetic field is strongest to the top where it is weakest.

The purpose of this section is to develop a formula describing the spatial distribution of the time-averaged heating rate due to Alfvén wave turbulence in the observed active region. The formula is intended to be used in future studies of both thermally stable and unstable coronal loops, and therefore must describe not only how the heating depends on position along each field line, but also how it varies from one field line to another. The formula will be derived from the numerical results described in the previous section (see Table 1); therefore, it is based on a realistic model for the 3D magnetic structure of the observed active region. The simulations provide heating rates Q_n(s) for all models, where n is the model index (n = 1, ..., N), and N is the total number of models (N = 77). The analysis proceeds in two steps. First, we consider the dependence of the heating rate on position along a loop. Then, we consider the variations between loops, i.e., the dependence of the heating on loop parameters such as coronal loop length L_n and coronal pressure p_n.

5.1. Variations along Field Lines

Our simulation results indicate that the heating rate Q_n(s) depends on position along the loop. Figure 4(a) shows the heating rates at the TR, Q_{TR, n} and at the loop top, Q_{min, n}, plotted as function of loop length L_n for all 77 models. Note that Q_TR is much larger than Q_min, and that both quantities decrease with loop length, but Q_min drops faster than Q_TR. The dependence of Q(s) on position is usually well described by a power law in terms of the magnetic field strength B_n(s), Equation (12). The magnitude of this variation can be parameterized by

$$\begin{equation} \eta _{Q,n} \equiv Q_{{\rm TR},n} / Q_{{\rm min},n}, \end{equation} \tag{ 13 }$$

where Q_{TR, n} is computed from a fit as described in Section 4, and Q_{min, n} is the minimum heating rate obtained from the same fit with B(s) = B_min. The exponent m in the relationship between Q_n(s) and B_n(s) is given by

$$\begin{equation} m_n = \frac{\log \eta _{Q,n}}{\log \Gamma _n}, \end{equation} \tag{ 14 }$$

where Γ_n ≡ B_{TR, n}/B_{min, n} is the expansion factor for model n. The ratio η_{Q, n} has been measured for each of the 77 models listed in Table 1. Figure 4(b) shows this quantity plotted as a function of coronal loop length L_n. Note that there is a strong trend with longer loops having larger η_{Q, n} ratios, i.e., stronger variations of heating rate Q_n(s) along the loop. We also plotted η_{Q, n} against other loop parameters, but found that the loop length L_n provides the strongest correlation.

Zoom In Zoom Out Reset image size

The simulation results have been fitted with the following formula:

$$\begin{equation} \eta _Q (L_n) = a_0 L_{5,n}^{a_1} + a_2 L_{5,n}^{a_3}, \end{equation} \tag{ 15 }$$

where L_{5, n} is the loop length in units of 10⁵ km for model n. We also assume that short loops (not simulated) must have nearly uniform heating, so we impose the further constraint that η_Q ≈ 1 for L = 10⁴ km and Equation (15) only applies to loops with L ⩾ 10⁴ km. The resulting fit has the following values for the coefficients:

$$\begin{equation} a_0 = 3.021, a_1 = 0.479, a_2 = 2.014, a_3 = 3.442, \end{equation} \tag{ 16 }$$

and is shown by the solid curve in Figure 4(b). Together with the PFSS model, Equations (15) and (16) can be used to predict the spatial variation of the heating (η_Q) for any field line in the observed active region.

5.2. Variations Across Field Lines

We now consider all spatial variations of the heating rate Q in the corona, including variations across magnetic field lines. Different formulae for Q with different dependencies on the coronal loop parameters were tested. As an example, let us consider the case where Q depends on the coronal loop length L and coronal pressure p. In this case the heating rate could be written in the following logarithmic form:

$$\begin{eqnarray} &&\log Q(s) = C_0 + C_1 \log L_{5} + C_2 \log p + m(L, \Gamma) \log \left[ \frac{B(s)}{B_{\rm TR}} \right], \nonumber \\ \end{eqnarray} \tag{ 17 }$$

where L₅ is the coronal loop length (in units of 10⁵ km), p is the coronal pressure (in dyne cm⁻²), and m(L, Γ) is computed in accordance with Equations (14) and (15):

$$\begin{equation} m(L, \Gamma) = \frac{\log \big(a_0 L_{5}^{a_1} + a_2 L_{5}^{a_3}\big)}{\log \Gamma }. \end{equation} \tag{ 18 }$$

The last term in Equation (17) describes the variations of the heating rate along the field lines, and is consistent with Equation (12).

We now describe how the constants C₀, C₁, and C₂ in Equation (17) are determined. Bringing the last term on the left-hand side of the equation, and applying the formula to model n, we can define

$$\begin{equation} Y_n (s) \equiv \log Q_n (s) - m(L_n, \Gamma _n) \log \left[ \frac{B_n (s)}{B_{{\rm TR},n}} \right], \end{equation} \tag{ 19 }$$

where Q_n(s) and B_n(s) are the heating rate and magnetic field strength as functions of position along the loop, B_{TR, n} is the field strength at the TR, and L_n is the loop length. By construction Y_n(s) is nearly constant along the coronal part of the loop; therefore, we average it over the coronal grid points to obtain a single value $\overline{Y}_n$ for each model. Then we fit the following formula to the measured values of $\overline{Y}_n$:

$$\begin{equation} \overline{Y}_{{\rm fit},n} = C_0 + C_1 \log L_{5,n} + C_2 \log p_n, \end{equation} \tag{ 20 }$$

where L_{5, n} is the loop length (in units of 10⁵ km), and p_n is the coronal pressure (in dyne cm⁻²). The coefficients C₀, C₁, and C₂ are computed by minimizing the following quantity:

$$\begin{equation} \chi ^2 = \frac{1}{N} \sum _{n=1}^{N} \left(\overline{Y}_{n} - \overline{Y}_{{\rm fit},n} \right)^2, \end{equation} \tag{ 21 }$$

which is a standard regression analysis. In the present case we find:

$$\begin{equation} C_0 = -2.43, C_1 = -0.95, C_2 = 0.19. \end{equation} \tag{ 22 }$$

The error of the fit for a particular model is $\Delta \overline{Y}_n = \overline{Y}_{{\rm fit},n} - \overline{Y}_{n}$, and averaged over models the error equals the minimum value of χ. The variances of the fit are defined by

$$\begin{equation} \sigma _{1}^2 \equiv C_{1}^2 \frac{1}{N} \sum _{n=1}^{N} \left(X_{1n} - \overline{X_1} \right)^2 = 0.05, \end{equation} \tag{ 23 }$$

and

$$\begin{equation} \sigma _{2}^2 \equiv C_{2}^2 \frac{1}{N} \sum _{n=1}^{N} \left(X_{2n} - \overline{X_2} \right)^2 = 0.0087, \end{equation} \tag{ 24 }$$

where X_1n = log L_{5, n}, X_2n = log p_n, and $\overline{X_1}$ and $\overline{X_2}$ are the averages of these quantities. Figure 5 shows the relationship between $\overline{Y}_n$ and $\overline{Y}_{{\rm fit},n}$. Note that the average error of the fit is less then 0.1 in the base 10 logarithm. The scatter in Figure 5 is due to several effects. First, the above Equation (17) is only an approximation and does not fit the numerical results perfectly. Second, there are real differences between loops that have the same length and pressure. For example, f9r3 and f19r2 have the same loop length and coronal pressure, but the heating for f9r3 is a factor of 2.75 larger than that for f19r2. This is due to the fact that field line F9 is rooted in strong magnetic field (sunspot), while F19 is low-lying and rooted in weaker fields.

Zoom In Zoom Out Reset image size

In the above discussion we used the coronal loop length L and pressure p as an example, but the same approach can be used with any combination of loop parameters.

We examined the dependency of the heating rate on different parameters of the field line. The result is illustrated in Table 2. The first four rows show the cases with two parameters from the selection of L, p, ρ_min, and τ_AC. The last four rows are cases with three model parameters including the averaged transition region magnetic field strength B_TR. C₀ is the constant in the heating equation (Equation (17)). The exponents in the equation are C₁, C₂, and C₃. The variances are σ₁², σ₂², and σ₃². The fit error χ is represented in the last column.

Table 2. Testing the Parameter Dependence of the Coronal Heating Equation (17)

Case	Par1	Par2	Par3	C0	C1	C2	C3	$\sigma _{1}^2$	$\sigma _{2}^2$	$\sigma _{3}^2$	χ
1	L	p	...	−2.43	−0.95 ± 0.24	0.19 ± 0.12	...	0.05	0.0087	...	0.092
2	L	ρmin	...	−2.28	−0.81 ± 0.26	0.24 ± 0.15	...	0.04	0.01	...	0.093
3	τAC	p	...	−2.60	−0.38 ± 0.09	0.32 ± 0.12	...	0.05	0.02	...	0.10
4	τAC	ρmin	...	−2.31	−0.30 ± 0.09	0.39 ± 0.14	...	0.03	0.03	...	0.089
5	L	p	BTR	−2.44	−0.95 ± 0.25	0.19 ± 0.12	0.01 ± 0.17	0.05	0.009	0.00002	0.091
6	L	ρmin	BTR	−2.30	−0.80 ± 0.27	0.24 ± 0.15	0.04 ± 0.17	0.03	0.01	0.0002	0.092
7	τAC	p	BTR	−2.51	−0.43 ± 0.11	0.35 ± 0.12	−0.19 ± 0.18	0.06	0.03	0.005	0.082
8	τAC	ρmin	BTR	−2.27	−0.31 ± 0.10	0.39 ± 0.14	−0.09 ± 0.18	0.04	0.03	0.001	0.085

Download table as:  ASCIITypeset image

From Table 2, case 1 is described in detail in this section. In case 2, we replace the gas pressure p with minimum density ρ_min. The exponent C₁ does not change very much. In fact C₁ is always negative for all cases, pointing to the fact that the heating rate varies inversely with the length of the loop. The exponent C₂ is positive, confirming the direct correlation between the heating rate and gas pressure. In cases 3 and 4, we substitute the length of the loop with τ_AC, the exponents C₁ and C₂ have similar values compared to the cases 1 and 2. The values of χ² indicate that models 1–4 are equally good.

For cases 5–8, we include B_TR. The interesting point is that the heating structure does not change significantly due to the effect of the inclusion of B_TR. The exponent C₃ corresponding to B_TR is very small and the error in the calculation of this term is larger than the exponent itself. This shows that the equation of the heating is better described by cases 1–4. Apparently, the magnetic field strength dependency of the heating rate equation is well described by the term B(s)/B_TR, and including the term B_TR, does not improve the fit.

Papers I and II demonstrated that our model of Alfvén wave turbulence predicts that the heating rate Q varies strongly with space and time within a given loop. Here we explore "macroscopic" variations in s and t (averaged over x and y) in Section 6.1, and then smaller "microscopic" variations resolved over the loop cross section in Section 6.2.

6.1. Variations over Time and Loop Length

The Alfvén wave turbulence model predicts that the heating rate Q varies strongly in space and time. We study the spatial and temporal intermittency of heating events along a representative coronal loop by computing the total energy released by heating in each discrete zone of the model. In other words, for the same model f19r2 shown in Figure 2, we compute the energy per zone

$$\begin{equation} E(s,t) =Q(s,t) \, \Delta s\, \pi R^{2} \, \Delta t \end{equation} \tag{ 25 }$$

where Δs is the distance between neighboring loop length grid zones (which varies as a function of s), πR² is the cross sectional area of the loop, and Δt = 1.08 s is the interval between output time steps. Figure 6 (top panel) shows the energy per zone as a function of position along the flux tube, averaged over time in the same way as was done in Figure 2. Figure 6 (middle panel) illustrates the dependence of E on position s (plotted as Alfvén wave travel time on the vertical axis) and time t (on horizontal axis). The grayscale image shows E(s, t) on a linear scale in units of 10²¹ erg. Because the vertical axis is given in terms of Alfvén wave travel time, the trajectories of individual wave "packets" are straight lines, and it is clear that the strongest heating events occur where upward and downward packets collide with one another.

Zoom In Zoom Out Reset image size

Figure 6 (left side of middle panel) also shows representative "boxes" that are defined to extend for 19.4 s in simulation time and 19.4 s in Alfvén wave travel time. This box size was chosen because it appears to encompass the duration and spatial extent of most individual discrete heating events that appear in the coronal part of the simulations. In order to study the intermittency of nanoflare-like energy release events, we binned the coronal part of the simulation into a grid of boxes of this size. There were 8 boxes along the loop (the first starting at travel time 153.6 s; the last ending at travel time 308.5 s), and there were 141 boxes in the time domain (the first starting at t = 248.4 s; the last ending at t = 2988.1 s). Thus, there were a total of 1128 boxes, each containing a given amount of energy dissipation. The results discussed here do not depend on our choice of starting or ending locations for the definition of the box edges.

The distribution of summed energy values per event (i.e., per box) range over two orders of magnitude, with a minimum value of 1.3 × 10²² erg, a maximum value of 1.4 × 10²⁴ erg, and a mean value of 2.2 × 10²³. These values are highly reminiscent of what is expected for classical nanoflares (Hannah et al. 2011). Figure 6 (bottom panel) shows the statistical distribution of energy values collected into bins of equal logarithmic spacing. The roughly Gaussian appearance of the distribution suggests a lognormal distribution of energy release events along this particular loop. This kind of distribution has been suggested both for statistics of solar irradiance variability (Pauluhn & Solanki 2007) and for magnetic field variations in the heliosphere (Burlaga 2001). However, we believe that the low energy cutoff of the distribution in this case is due to the limited dynamic range of the present simulations (i.e., a spatial resolution of about 10% of the tube radius). Thus, we anticipate that when all relevant scales are included, the distribution of energy per event should become broader and possibly approach a power-law shape as has been found to occur in many other forms of intermittent energy release.

6.2. Variations over Loop Cross Section

To understand how the plasma inside the loop is heated, it is important to consider the local heating rate Q(x, y, s, t), which is also a function of the transverse coordinates x and y. This rate is computed using expressions (9) and (10). Figure 7 shows results for model f19r2 for a particular instant at the end of the simulation (see also the animated version included in the online journal). Figure 7(a) shows the normalized heating rate integrated over the x coordinate:

$$\begin{equation} I(s,y) = \frac{\int Q(x,y,s,t) dx}{R(s) Q(s)}, \end{equation} \tag{ 26 }$$

where s is the horizontal coordinate in the figure, y is the vertical coordinate in the figure, and Q(s) is the time-averaged heating rate. In this figure the transverse and longitudinal variations of the heating are shown on the same spatial scale, but the loop curvature is neglected. Only the coronal portion of the loop is shown; the TRs are located at the left and right edges of the diagram. Figure 7(b) shows the same quantity on an expanded scale with s replaced by the Alfvén travel time τ(s) and the y coordinate normalized to the tube radius R(s), and Figure 7(c) shows a similar plot for the heating integrated over the y coordinate. In these figures the spatial variations of the heating within the loop can be seen more clearly. Note that there is a strong heating pulse at position s ≈ 0.4L from the left TR, and weaker wave-like disturbances at other places along the loop.

Zoom In Zoom Out Reset image size

Figure 7(d) shows the normalized heating rate in six different cross-sections along the loop, I(x, y) = Q(x, y, s, t)/Q(s). The positions of these cross-sections are indicated by thin lines between panels (c) and (d). Note that the heating occurs in thin ridges and rings that are located near current sheets and shear layers (Oughton et al. 2001; Mininni & Pouquet 2009; Servidio et al. 2011), i.e., sites of strong magnetic and velocity shear. The shapes of these rings are different in the different cross-sections, and vary rapidly with time. In our simulations the ridges cover a significant fraction of the loop cross-section, so the heating is not highly localized within the loop. This may be an artifact of our limited spatial resolution of the waves, together with the fact that we use hyperdiffusion to describe the wave damping. We have simulated the Alfvén wave turbulence with the use of hyperdiffusion as in Papers I and II.

An animated version of Figure 7 is included in the online journal. The animation covers a period of about 432 s near the end of the 3000 s simulation, and shows the spatio-temporal variations of the heating Q(x, y, s, t) as seen from the side (panels (a), (b), and (c)) and in the six cross-sections (panel (d)). The animation shows wave-like disturbances emerging from the two TRs and propagating into the corona. The injection of waves into the corona is a highly intermittent process, which is a consequence of the random nature of the footpoint motions and the turbulence in the chromosphere. The strongest brightenings occur when oppositely directed wave packets "collide" with each other. This is a result of nonlinear wave–wave interactions, which produce a rapid cascade to higher wavenumbers and dissipation of the wave energy.

In our modeling, we simulated the dissipation and propagation of Alfvén waves in coronal loops. In order to fully understand Alfvén wave dynamics in coronal loops we included the lower atmospheres at the two ends of a loop as in Papers I and II. The waves are launched from the photosphere, and are simulated for a period of 3000 s. Along the field lines, the Alfvén speed varies from about 10 to 15 km s⁻¹ in the photosphere to more than 1000 km s⁻¹ in the low corona. Therefore, waves traveling to the corona encounter strong wave reflection. We find that the waves in the corona build up to significant amplitudes and have strong nonlinear interactions. This results in the dissipation of Alfvén waves. As shown in Paper II, the dissipation rate of Alfvén wave turbulence is sufficient to produce the observed rates of chromospheric and coronal heating in active regions.

In this work as well as identifying the mechanism responsible for the energy input in the corona and the subsequent heating of the coronal loops, we predicted the heating rate based on a photospheric magnetogram and derived an equation describing the dependence of the heating rate on loop parameters. Specifically, we simulated the dynamics of Alfvén waves for 22 field lines in an active region observed on 2012 March 7. We looked at the variations of the wave heating rate for these 22 field lines with field strength, loop length and coronal density. The results for different loops were combined into a formula describing the average heating rate Q(s) as function of position s along the loop. As we show in the Alfvén wave modeling of the coronal heating, the wave dissipation varies strongly with position along the loop. In the heating Equation (17), the dependence on the position is represented by the term B(s) where magnetic field strength is a function of position along the loop and the realistic models of coronal loops all point to this variation. As Equation (17) indicates, the heating also depends on the loop parameters such as coronal loop length L and coronal pressure p. Results for other combinations of parameters are shown in Table 2. The spatial variation of the heating rate is not arbitrarily prescribed but is a natural consequence of the Alfvén wave turbulence model. It is important to note that the loops heating profile plays an important role in determining whether or not the loops are thermally stable.

Equation (17) and similar expressions listed in Table 2 are intended to be used for constructing 3D time-dependent models of coronal plasma in an observed active region, and could also be used for other regions. Such modeling requires solving the energy transport equation for coronal loops, i.e., the balance between non-thermal heating, radiative losses, and thermal conduction. In general both thermally stable and unstable coronal loops will be present, and for the latter the plasma temperature and density will vary strongly with time. Therefore, in the present modeling the coronal pressure p was treated as a free parameter, and the loops were not constrained to be in thermal equilibrium. Computations were done for a wide range of coronal pressures, so that the heating-rate formulae would accurately describe the effects of the coronal density on the propagation and dissipation of the Alfvén waves for a wide range of coronal conditions.

Mandrini et al. (2000) investigate the coronal heating dependence on the magnetic field strength and coronal loop length. Schrijver et al. (2004) compare various models of the heating with TRACE observations of the corona, and found that the heating flux density F_H varies as B/L, in agreement with one of the models studied by Mandrini et al. (2000) and Démoulin et al. (2003). This corresponds to a heating rate per unite volume that depends strongly on loop length: Q∝L⁻². However, the above mentioned authors do not include the variation along the loop. In our studies, the heating rate formula has a weaker dependency on the loop length than that found by Schrijver et al. (2004) and the magnetic field dependency of our heating expression is well represented by the term B(s)/B_TR in Equation (17).

In our present work we do not address the question how the coronal plasma responds to the heating events. A comprehensive overview of the coronal heating problem is presented by Klimchuk (2006), who argues that loops are heated by "storms" of nanoflares. Warren et al. (2002) explored the possibility that active region loops are a collection of small-scale filaments that have been impulsively heated. They found that the density variations lag the temperature variations by several hundred seconds, and during the cooling phase the filaments are significantly denser than steadily heated loops with the same temperature. This is important for understanding the brightness of observed loops at EUV wavelengths. Reale et al. (2005) were the first to have simulated the response of a loop to nanoflare heating associated with MHD turbulence, and they simulated the emission in EUV passbands. The modeled coronal temperatures and pressures are somewhat lower than those observed in active regions; we suggest this may be due to the fact that this model does not include wave amplification in the lower atmosphere. Future modeling of Alfvén wave turbulence should take into account the dynamic response of the plasma to the heating events, including the effects of spatial variations of temperature and density over the loop cross-section.

In the present model the TR is treated as a discontinuity, and energy exchange between the chromosphere and corona is not included. However, thermal conduction and mass flows along the field lines are likely to play an important role in the structure of the TR. The wave turbulence model predicts that the wave heating rate in the chromosphere is orders of magnitude larger than that in the corona, e.g., see Figure 2(f). In reality some of this energy may be injected into the low corona by spicules or similar dynamic phenomena (De Pontieu et al. 2007a; McIntosh et al. 2011). This may alter the way Alfvén waves reflect at the TR and may increase the fraction of wave energy dissipated in the corona.

In our future work we are intending to construct models of the thermal structure of coronal loops heated by Alfvén wave turbulence where we include the effect of conduction and radiative losses and also make use of Doppler shift measurements as an important constraints on our coronal loop modeling.

We thank James Klimchuk and Peter Cargill for discussions related to this work. We are most grateful to Alex Voss from the School of Computer Science at the University of St. Andrews for his support with the computational work, which was funded by the UK's Engineering and Physical Sciences Research Council (EP/I034327/1). The HMI and AIA data have been used courtesy of NASA/SDO and the AIA, EVE, and HMI science teams.

This project is supported under contract NNM07AB07C from NASA to the Smithsonian Astrophysical Observatory (SAO) and SP02H1701R from LMSAL to SAO.