Turbulent coronal heating mechanisms: coupling of dynamics and thermodynamics
^{1} Laboratory for Computational Physics and Fluid Dynamics, Naval Research Laboratory, Washington, DC 20375, USA
email: rdahlbur@lcp.nrl.navy.mil
^{2} Berkeley Research Associates, Inc., 6537 Mid Cities Avenue, Beltsville, MD 20705, USA
^{3} Bartol Research Institute, Department of Physics and Astronomy, University of Delaware, DE 19716, USA
^{4} Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
Received: 4 June 2012
Accepted: 23 July 2012
Context. Photospheric motions shuffle the footpoints of the strong axial magnetic field that threads coronal loops, which gives rise to turbulent nonlinear dynamics that are characterized by the continuous formation and dissipation of fieldaligned current sheets in which energy is deposited at smallscales and the heating occurs. Previous studies showed that the current sheet thickness is several orders of magnitude smaller than presentday stateoftheart observational resolution (~700 km).
Aims. To understand coronal heating and correctly interpret observations it is crucial to study the thermodynamics of such a system in which energy is deposited at unresolved smallscales.
Methods. Fully compressible threedimensional magnetohydrodynamic simulations were carried out to understand the thermodynamics of coronal heating in the magnetically confined solar corona.
Results. We show that temperature is highly structured at scales below observational resolution. It is also nonhomogeneously distributed so that only a fraction of the coronal mass and volume is heated at each time.
Conclusions. This is a multithermal system in which hotter and cooler plasma strands are also found next to each other at subresolution scales and exhibit a temporal dynamics.
Key words: Sun: corona / magnetohydrodynamics (MHD) / turbulence
© ESO, 2012
1. Introduction
Magnetohydrodynamic (MHD) turbulence has long been thought of as a key process in coronal heating (Einaudi et al. 1996) as well as fast and slow solar wind acceleration. Numerical simulations have been crucial in investigating the role of MHD turbulence.
To tackle the turbulence problem correctly with direct numerical simulations, it is necessary to resolve enough spatial scales to obtain an energycontaining range, an inertial range, and a dissipation range. Because of storage needs, the more the complexity of the problem is reduced, the more spatial scales are obtainable. Hence most previous research in this area has been restricted to dynamical effects, e.g., reduced MHD (RMHD; Strauss 1976) and cold plasma models. This considerably reduces the size of the problem. For example, for RMHD the dynamics of the system are determined by two variables: a stream function and a magnetic potential function (Rappazzo et al. 2007). For the cold plasma model it is customary to evolve three velocity field components and three magnetic field components, neglecting the pressure term (i.e., not correcting for the solenoidality of the velocity field) so that six variables are used (Dahlburg et al. 2009).
In all these approximations all thermodynamical aspects are lost, i.e., there is no information on the connection between the dynamics of the system, which is responsible for the heating mechanism, and the temperature and density behavior. In particular, energy lost through Ohmic and viscous diffusion is simply lost from the system and the thermodynamic physics involved in thermal conduction and optically thin radiation is neglected.
This lack of information prevents a complete reproduction of the energy cycle: photospheric motions shuffle magnetic fieldline footpoints, which injects energy into the corona and triggers nonlinear dynamics that form smallscales, which are organized in current sheets (Parker 1972), in which energy is then transformed into thermal (and kinetic and perturbed magnetic) energy by means of magnetic reconnection. Heat is then conducted from the hightemperature corona back toward the lowtemperature photosphere where it is lost via optically thin radiation, but at the same time causes chromospheric evaporation.
Thermodynamics are typically studied with onedimensional hydrodynamic models (Reale 2010) that study either the equilibrium condition of a heated loop, or the timedependent evaporative and possible condensation flows. However, they introduce an ad hoc heating function to mimic the energy deposition and neglect the crossfield dynamics (including fieldline reconnection). It is therefore crucial to use the computationally more demanding threedimensional compressible MHD (with eight variables), including an energy equation with thermal conduction and the energy sink provided by optically thin radiation.
The RMHD turbulence simulations (Rappazzo et al. 2007, 2008) have shown that in the planes orthogonal to the DC magnetic field the length of current sheets is on the order of the photospheric convective cells scale (i.e., ~1000 km), while their width is much smaller and limited only by numerical resolution. To understand the thermodynamics of this system, we carried out simulations that resolve the scales below the convection scale. The boundary conditions used do not allow evaporative flows to develop along the heated loops. Our focus here is on the relationship of turbulent current sheet dissipation and the formation of the coronal temperature structure.
Our new compressible code HYPERION has allowed us to make a start at examining the fully compressible, threedimensional Parker coronal heating model. HYPERION is a parallelized Fourier collocation – finite difference code with thirdorder RungeKutta time discretization that solves the compressible MHD equations with vertical thermal conduction and radiation included. The results provided by this new compressible code preserve the spatial intermittency of the kinetic and magnetic energies seen in earlier RMHD simulations, but also provide new information about the evolution of the internal energy, i.e., we can now determine how the coronal plasma heats up and radiates energy as a consequence of photospheric convection of magnetic footpoints (Dahlburg et al. 2010).
In Sect. 2 we describe the governing equations, the boundary and initial conditions, and the numerical method. In Sect. 3 we detail our numerical results. Section 4 contains a discussion of these results and some ideas about future directions of this research.
2. Formulation of the problem
2.1. Governing equations
We model the solar corona as a compressible, dissipative magnetofluid. The equations that govern such a system, written here in dimensionless form, are: with the solenoidality condition ∇ B = 0. The system is closed by the equation of state (5)The nondimensional variables are defined in the following way: n(x,t) is the numerical density, v(x,t) = (u,v,w) is the flow velocity, p(x,t) is the thermal pressure, B(x,t) = (B_{x},B_{y},B_{z}) is the magnetic induction field, J = ∇ × B is the electric current density, T(x,t) is the plasma temperature, ζ_{ij} = μ(∂_{j}v_{i} + ∂_{i}v_{j}) − λ∇·vδ_{ij} is the viscous stress tensor, e_{ij} = (∂_{j}v_{i} + ∂_{i}v_{j}) is the strain tensor, and γ is the adiabatic ratio. To render the equations dimensionless we used n_{0}, the photospheric numerical density, V_{A}, the vertical Alfvén speed at the photospheric boundaries, L_{0}, the orthogonal box length (= L_{x} = L_{y}) and T_{0}, the photospheric temperature. Therefore time (t) is measured in units of Alfvén cross time (τ_{A} = L_{0}/V_{A}).
The function Λ(T) that describes the temperature dependence of the radiation is taken following Hildner (1974), normalized with its value at the base temperature T_{0} = 10 000 K. The term C_{N} denotes a Newton thermal term, which is enforced at low temperatures to take care of the delicate chromospheric and transition region energy balance (Dorch & Nordlund 1999). We set C_{N} = 100 [T_{i}(z) − T(z)] e^{ − 4(z + 0.5 Lz)} at the lower wall and C_{N} = 100 [T_{i}(z) − T(z)] e^{ − 4(0.5 Lz − z)} at the upper wall, where T_{i}(z) is the initial temperature profile.
The important dimensionless numbers are: S_{v} = ρ_{0}V_{A}L_{0}/μ ≡ viscous Lundquist number, S = μ_{0}V_{A}L_{0}/η ≡ Lundquist number, pressure ratio at the wall, Prandtl number, and P_{rad}, the radiative Prandtl number , determines the strength of the radiation. C_{v} is the specific heat at constant volume. The term “κ” denotes the thermal conductivity. The magnetic resistivity (η) and shear viscosity (μ) are assumed to be constant and uniform, and Stokes relationship is assumed therefore the bulk viscosity λ = (2/3)μ.
2.2. Initial and boundary conditions
We solve the governing equations in a Cartesian box of dimensions 1 × 1 × L_{z}, where L_{z} is the loop aspect ratio (0 ≤ x,y ≤ 1, − L_{z}/2 ≤ z ≤ L_{z}/2), threaded by a DC magnetic field B_{0} in the zdirection. The system has periodic boundary conditions in x and y, and linetied boundary conditions at the top and bottom zplates where the vertical flows vanish (v_{z} = 0) and the forcing velocity in the xy plane is imposed as in previous studies (Hendrix et al. 1996; Einaudi et al. 1996; Einaudi & Velli 1999), evolving the stream function (6)where , and v = ∇ψ × ê_{z}. All wavenumbers with are excited, so that the typical lengthscale of the eddies is ~1/4. Because the typical convective cell size is ~1000 km and we use L_{z} = 5, this implies that our computational box in conventional units spans 4000^{2} × 20 000 km^{3}. Every t^{∗}, the coefficients and the are randomly changed alternatively for eddies 1 and 2. The magnetic field is expressed as B = B_{0}ê_{z} + b, where A is the vector potential associated with the fluctuating magnetic field b = ∇ × A. At the top and bottom zplates B_{z}, n and T are kept constant at their initial values B_{0}, n_{0} and T_{0}, while the magnetic vector potential is convected by the resulting flows.
Below we have assumed the normalizing quantities to be n_{0} = 10^{17} m^{3}, T_{0} = 10^{4} K, B_{0} = 10^{2} tesla, and L_{0} = 4 × 10^{6} m. It follows that the values of the various parameters appearing in the equations are β = 1.7 × 10^{4}, v_{A} = 6.9 × 10^{5} m/s, τ_{A} = 5.8 s, S = 2.7 × 10^{9}, S_{v} = 2.09 × 10^{9}, Pr = 1.52 × 10^{2}, P_{rad} = 3.35 × 10^{6}, , lnΛ = 10.
The normalized time scale of the forcing, t^{∗}, is set to 51.7 to represent the 5 min of a typical photospheric convection time scale. The normalized photospheric velocity is V_{0} = 1.45 × 10^{3}, corresponding to 10^{3} m/s.
We impose as initial conditions a temperature profile with a base temperature T_{0} = 10^{4} K and a top temperature 8 × 10^{5} K with the following zdependence for the nondimensional temperature T_{i}(z) = 1 + 79cos(πz/L_{z}), while the numerical density is given as n_{i}(z) = 1/T_{i}(z). This choice allows us to neglect the gravitational effects since for most of the loop the temperature remains sufficiently high for the gravitational lengthscale to exceed L_{z} at all times.
Fig. 1 Temperature contours at t = 300 in the midplane (z = 0). Maximum (130) and minimum (49) temperatures correspond to 1.3 million and 490 thousand degrees Kelvin in conventional units. Distances are normalized by L_{0} = 4 × 10^{6} m. 

Open with DEXTER 
3. Results
With previous definitions, Eq. (4) can be replaced by the magnetic vector potential equation: (7)We numerically solve Eqs. (1)–(3) and (7) together with Eq. (5). Space is discretized in x and y with a Fourier collocation scheme (Dahlburg & Picone 1989) with isotropic truncation dealiasing. Spatial derivatives are calculated in Fourier space, and nonlinear product terms are advanced in configuration space. A secondorder central difference technique (Dahlburg et al. 1986) is used for the discretization in z. A timestep splitting scheme is employed. The code was parallelized using MPI. A domain decomposition is employed in which the computational box is sliced into x–y planes along the z direction.
We present the results obtained running the code with S = S_{v} = 4 × 10^{4} with a resolution of 128^{3}. To keep the efficiency of the radiative and conductive terms in the energy equation the same as in the real corona, we rescaled Pr and P_{rad} accordingly with the choice of S_{v}, i.e., Pr = 793.39 and P_{rad} = 0.175. This choice was motivated by the result found in the RMHD model (Rappazzo et al. 2008), which is that turbulent dissipative processes are independent of viscosity and resistivity when an inertial range is well resolved.
The system begins in a ground state, threaded by a guide magnetic field: the interior of the channel has zero perturbed velocity and magnetic fields. The footpoints of the magnetic field are subjected to convection at the z boundaries with initial number density and temperature profiles assigned as described above. The behavior of the volumeaveraged quantities, such as kinetic and fluctuating magnetic energies and resistive and viscous dissipation, show a temporal behavior similar to the RMHD results previously obtained by Rappazzo et al. (2007, 2008, 2010).
The fluctuating kinetic and magnetic energies ( and ), as found in the RMHD model, and the total internal energy (U = ⟨ nT ⟩ β/(γ − 1)), not obtainable in the RMHD case, exhibit an intermittent behavior. After the system is driven by photospheric motions for more than 100 Alfvén times, a statistically stationary state is attained and all quantities fluctuate around their steady values.
Fig. 2 Maximum temperature (T_{max}) and maximum electric current  j  _{max} as functions of time. 

Open with DEXTER 
We verified that the dynamics of this problem are well represented by RMHD. We defined the volumeaveraged velocity θ_{v} and magnetic θ_{b} anisotropy angles as in Shebalin et al. (1983). For isotropic turbulence, θ_{v} and θ_{b} equal 45°. For fully anisotropic turbulence, θ_{v} and θ_{b} equal 90°, implying that their spectra would be normal to B_{0} ê_{z}. In the present computation, we find that, after t = 100, θ_{b} is always bigger than 87°, while the velocity field fluctuations are found to be slightly more anisotropic than those in the magnetic field, 88° < θ_{v} < 89°, indicating that the turbulence in the present system is highly anisotropic, i.e., the turbulence is strongly confined to planes orthogonal to the guide magnetic field. Moreover, we verified that the dynamics is almost incompressible.
The important new result of this letter concerns the thermodynamical behavior of the system, which results from the selfconsistent heating due to the magnetic fieldline tangling induced by photospheric motions, once the cooling effects of conduction and radiation are taken into account. The fact that the total internal energy, as mentioned above, fluctuates around a steady value means that the selfconsistent heating mechanism caused by photospheric motions is efficient enough to sustain a hot corona once the system has reacted to photospheric motion. Such a reaction leads to the formation of small scales in which energy can be efficiently dissipated. The small scales, corresponding to the presence of current sheets, are not uniformly distributed in the system, and as a result, the temperature is spatially intermittent, as shown in Fig. 1, where we show the midplane temperature contours at t = 300.
Figure 2 shows the time evolution of T_{max} and  j  _{max}, the maximum temperature and the maximum current density present in the loop at each time. Evidently the T_{max} temperature peaks oscillate in time and, after an initial transient, the T_{max} values are well above T = 10^{6} K, higher than the initial peak temperature of T = 8 × 10^{5} K degrees. An important fact to emphasize is that the peaks in T_{max} roughly coincide with bursts in the maximum electric current. Since currents are organized in sheets that are elongated along the mean magnetic field and the heating is caused by the dissipation of currents, we expect the heating not to be homogeneously distributed throughout the loop.
Fig. 3 Temperature is organized in quasi2D pancakelike structures aligned with the magnetic field. The two isosurfaces for T = 10^{6} K and T = 1.2 × 10^{6} K at time t = 300 are shown in yellow and gray (corresponding to T = 100,120 in nondimensional units), respectively. Sample magnetic fieldlines are drawn in red. The box has been rescaled in z to clarify the visualization. 

Open with DEXTER 
This expectation is confirmed by the threedimensional snapshot of the isosurfaces of two temperatures (T = 10^{6} and T = 1.2 × 10^{6} K) at t = 300 in Fig. 3. The spatial temperature morphology shows that the high temperatures are reached in selected quasi2D regions elongated along the mean magnetic field in which the system dynamics at each time places the current sheets in which energy is dissipated. The higher temperatures are found at the center of these structures (up to 1.35 × 10^{6} K at t = 300), while farther out the temperature decreases to a background value that in the central zregion averages about 6 × 10^{5} K. Snapshots taken at different times show similar behavior with different high temperatures locations.
4. Conclusions and discussion
We presented for the first time the thermodynamic behavior of a system, with the spatial scales below the convection scale resolved numerically, in which the conductive and radiative losses are balanced by the heating process due solely to the photospheric motions and the subsequent tangling of magnetic field lines, While previous compressible simulations (Gudiksen & Nordlund 2002; Bingert & Peter 2011) have considered entire active regions and resolved at most the convective scale, it is pivotal to resolve smaller scales to understand how the plasma is heated and the observational consequences of this.
Much work has been done in the past to show the reaction of a plasma that is subject to mass flow at its boundaries. It has been shown that the system reacts by developing a weak turbulent state, in which energy is dissipated in the current sheets that result from the dynamics induced by the boundary velocity forcing. The results of the present work show that this process leads to a loop in which the temperature is spatially highly structured.
A loop is therefore a multithermal system, in which the temperature peaks around current sheets and exhibits a distribution of temperatures. The characteristics of this distribution (peak, width, background) will vary depending on the parameters of the loop (length, boundary velocity, mean magnetic field).
Temperature and heating are highly nonhomogeneously distributed so that only a fraction of the loop volume exhibits a significant heating at each time. Therefore emission observed in a spectral band originates from the small subvolume at the corresponding temperature, whose mass is considerably lower than the total loop mass.
Observations currently have a spatial resolution of ~700 km and a temporal resolution of ~30 s. These temperature structures are therefore unresolved, and to interpret observations correctly the spatial and temporal intermittency must be taken into account. We intend to present in a following paper how our loop appears once the emission is averaged over the coarser observational resolutions. At those scales the loop might become multithermal or approximately isothermal, depending on the loop parameters.
It is very likely that the local thermodynamical equilibrium is lost and that the emission is caused by particles that are accelerated in the current sheets and not by direct transformation of magnetic energy into thermal energy. This is beyond the MHD model.
We have shown that the heating caused by photospheric motions is balanced by the loss mechanisms. Furthermore, the photospheric motions induce a multiscale dynamics that can account for the observed emission.
The density stratification in our simulation allows limited evaporation from the zboundary edges, where the density is higher but does not reach typical chromospheric values. In future work we plan to implement characteristic boundary conditions (e.g., Rappazzo et al. 2005) that allow higher density flows through the zboundaries, which can also give rise to condensation and prominences, and to study their impact on the thermodynamics.
Acknowledgments
This work was carried out in part at the Jet Propulsion Laboratory under a contract with NASA. R.B.D. thanks P. Byrne for helpful conversations.
References
 Bingert, S., & Peter, H. 2011, A&A, 530, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dahlburg, R. B., & Picone, J. M. 1989, Phys. Fluids B, 1, 2153 [NASA ADS] [CrossRef] [Google Scholar]
 Dahlburg, R. B., Montgomery, D., & Zang, T. A. 1986, J. Fluid Mech., 169, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Dahlburg, R. B., Liu, J.H., Klimchuk, J. A., & Nigro, G. 2009, ApJ, 704, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Dahlburg, R. B., Rappazzo, A. F., & Velli, M. 2010, AIP Conf. Proc., 1216, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Dorch, S.B.F., & Nordlund, Å. 2001, A&A, 365, 562 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Einaudi, G., & Velli, M. 1999, Phys. Plasmas, 6, 4146 [NASA ADS] [CrossRef] [Google Scholar]
 Einaudi, G., Velli, M., Politano, H., & Pouquet, A. 1996, ApJ, 457, L113 [NASA ADS] [CrossRef] [Google Scholar]
 Gudiksen, B. V., & Nordlund, Å. 2002, ApJ, 572, L113 [NASA ADS] [CrossRef] [Google Scholar]
 Hendrix, D. L., Van Hoven, G., Mikić, Z., & Schnack, D. D. 1996, ApJ, 470, 1192 [NASA ADS] [CrossRef] [Google Scholar]
 Hildner, E. 1974, Sol. Phys., 35, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1972, ApJ, 174, 499 [NASA ADS] [CrossRef] [Google Scholar]
 Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2005, ApJ, 633, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2007, ApJ, 657, L47 [NASA ADS] [CrossRef] [Google Scholar]
 Rappazzo, A. F., Velli, M., Einaudi, G., & Dahlburg, R. B. 2008, ApJ, 677, 1348 [NASA ADS] [CrossRef] [Google Scholar]
 Rappazzo, A. F., Velli, M., & Einaudi, G. 2010, ApJ, 722, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Reale, F. 2010, Liv. Rev. Sol. Phys., 7, 5 [Google Scholar]
 Shebalin, J. V., Matthaeus, W. H., & Montgomery, D. 1983, J. Plasma Phys., 29, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Strauss, H. R. 1976, Phys. Fluids, 19, 134 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Temperature contours at t = 300 in the midplane (z = 0). Maximum (130) and minimum (49) temperatures correspond to 1.3 million and 490 thousand degrees Kelvin in conventional units. Distances are normalized by L_{0} = 4 × 10^{6} m. 

Open with DEXTER  
In the text 
Fig. 2 Maximum temperature (T_{max}) and maximum electric current  j  _{max} as functions of time. 

Open with DEXTER  
In the text 
Fig. 3 Temperature is organized in quasi2D pancakelike structures aligned with the magnetic field. The two isosurfaces for T = 10^{6} K and T = 1.2 × 10^{6} K at time t = 300 are shown in yellow and gray (corresponding to T = 100,120 in nondimensional units), respectively. Sample magnetic fieldlines are drawn in red. The box has been rescaled in z to clarify the visualization. 

Open with DEXTER  
In the text 