© 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 energy-containing 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 field-line footpoints, which injects energy into the corona and triggers non-linear dynamics that form small-scales, 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 high-temperature corona back toward the low-temperature photosphere where it is lost via optically thin radiation, but at the same time causes chromospheric evaporation.

Thermodynamics are typically studied with one-dimensional hydrodynamic models (Reale 2010) that study either the equilibrium condition of a heated loop, or the time-dependent evaporative and possible condensation flows. However, they introduce an ad hoc heating function to mimic the energy deposition and neglect the cross-field dynamics (including field-line reconnection). It is therefore crucial to use the computationally more demanding three-dimensional 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, three-dimensional Parker coronal heating model. HYPERION is a parallelized Fourier collocation – finite difference code with third-order Runge-Kutta 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: $\begin{array}{ccc}\frac{\mathit{\partial n}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{-}\mathrm{\nabla }\hspace{0.17em}\mathrm{\left(}\mathit{n}{v}\mathrm{\right)}\mathit{,}\\ \frac{\mathit{\partial n}{v}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{-}\mathrm{\nabla }\hspace{0.17em}\mathrm{\left(}\mathit{n}{vv}\mathrm{\right)}\mathrm{-}\mathit{\beta }\mathrm{\nabla }\mathit{p}\mathrm{+}{J}\mathrm{\times }{B}\mathrm{+}\frac{\mathrm{1}}{{\mathit{S}}_{\mathit{v}}}\mathrm{\nabla }\hspace{0.17em}{\zeta }\mathit{,}\\ \frac{\mathit{\partial T}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{-}{v}\hspace{0.17em}\mathrm{\nabla }\mathit{T}\mathrm{-}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\left(}\mathrm{\nabla }\hspace{0.17em}{v}\mathrm{\right)}\mathit{T}\mathrm{+}\frac{\mathrm{1}}{\mathit{n}}\{\frac{\mathrm{1}}{\mathit{Pr}\hspace{0.17em}{\mathit{S}}_{\mathit{v}}}\frac{\mathit{\partial }}{\mathit{\partial z}}\left({\mathit{T}}^{\frac{\mathrm{5}}{\mathrm{2}}}\frac{\mathit{\partial T}}{\mathit{\partial z}}\right)\\ & & \mathrm{+}\frac{\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{\beta }}[\frac{\mathrm{1}}{{\mathit{S}}_{\mathit{v}}}{\mathit{\zeta }}_{\mathit{ij}}{\mathit{e}}_{\mathit{ij}}\mathrm{+}\frac{\mathrm{1}}{\mathit{S}}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{B}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{P}}_{\mathrm{rad}}{\mathit{S}}_{\mathit{v}}}{\mathit{n}}^{\mathrm{2}}\mathrm{\Lambda }\mathrm{(}\mathit{T}\mathrm{)}\\ & & \mathrm{+}\frac{\mathit{\beta }}{\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{n}{\mathit{C}}_{\mathrm{N}}]\}\mathit{,}\\ \frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& \end{array}$with the solenoidality condition ∇   B = 0. The system is closed by the equation of state $\mathit{p}\mathrm{=}\mathit{nT}\mathit{.}$(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 = μ(∂_jv_i + ∂_iv_j) − λ∇·vδ_ij is the viscous stress tensor, e_ij = (∂_jv_i + ∂_iv_j) is the strain tensor, and γ is the adiabatic ratio. To render the equations dimensionless we used n₀, the photospheric numerical density, V_A, the vertical Alfvén speed at the photospheric boundaries, L₀, the orthogonal box length (= L_x = L_y) and T₀, the photospheric temperature. Therefore time (t) is measured in units of Alfvén cross time (τ_A = L₀/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₀ = 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   L_z)} at the lower wall and C_N = 100    [T_i(z) − T(z)] e ^{− 4(0.5   L_z − z)} at the upper wall, where T_i(z) is the initial temperature profile.

The important dimensionless numbers are: S_v = ρ₀V_AL₀/μ ≡  viscous Lundquist number, S = μ₀V_AL₀/η ≡  Lundquist number, $\mathit{\beta }\mathrm{=}{\mathit{\mu }}_{\mathrm{0}}{\mathit{p}}_{\mathrm{0}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\equiv }$ pressure ratio at the wall, $\mathit{Pr}\mathrm{=}{\mathit{C}}_{\mathit{v}}\mathit{\mu }\mathit{/}\mathit{\kappa }{\mathit{T}}_{\mathrm{0}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{\equiv }$ Prandtl number, and P_rad, the radiative Prandtl number $\mathit{\mu }\mathit{/}{\mathit{\tau }}_{\mathrm{A}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{\Lambda }\mathrm{\left(}{\mathit{T}}_{\mathrm{0}}\mathrm{\right)}$, 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₀ in the z-direction. The system has periodic boundary conditions in x and y, and line-tied boundary conditions at the top and bottom z-plates where the vertical flows vanish (v_z = 0) and the forcing velocity in the x-y plane is imposed as in previous studies (Hendrix et al. 1996; Einaudi et al. 1996; Einaudi & Velli 1999), evolving the stream function $\mathit{\psi }\mathrm{\left(}\mathit{x,y,t}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{1}}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\pi t}}{\mathrm{2}{\mathit{t}}^{\mathrm{\ast }}}\right)\mathrm{+}{\mathit{f}}_{\mathrm{2}}{\sin }^{\mathrm{2}}\left(\frac{\mathit{\pi t}}{\mathrm{2}{\mathit{t}}^{\mathrm{\ast }}}\mathrm{+}\frac{\mathit{\pi }}{\mathrm{2}}\right)\mathit{,}$(6)where ${\mathit{f}}_{\mathit{i}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}{\sum }_{\mathit{n,m}}{\mathit{a}}_{\mathit{nm}}^{\mathit{i}}\sin \mathrm{(}{\mathit{k}}_{\mathit{n}}\mathit{x}\mathrm{+}{\mathit{k}}_{\mathit{m}}\mathit{y}\mathrm{+}{\mathit{\zeta }}_{\mathit{nm}}^{\mathit{i}}\mathrm{)}$, and v = ∇ψ × ê_z. All wave-numbers with $\mathrm{3}\mathrm{\le }\mathrm{(}{\mathit{k}}_{\mathit{n}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{m}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\le }\mathrm{4}$ are excited, so that the typical length-scale 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²  ×  20   000 km³. Every t^∗, the coefficients ${\mathit{a}}_{\mathit{nm}}^{\mathit{i}}$ and the ${\mathit{\zeta }}_{\mathit{nm}}^{\mathit{i}}$ are randomly changed alternatively for eddies 1 and 2. The magnetic field is expressed as B = B₀ê_z + b, where A is the vector potential associated with the fluctuating magnetic field b = ∇ × A. At the top and bottom z-plates B_z, n and T are kept constant at their initial values B₀, n₀ and T₀, while the magnetic vector potential is convected by the resulting flows.

Below we have assumed the normalizing quantities to be n₀ = 10¹⁷   m^-3, T₀ = 10⁴   K, B₀ = 10^-2 tesla, and L₀ = 4 × 10⁶   m. It follows that the values of the various parameters appearing in the equations are β = 1.7 × 10^-4, v_A = 6.9 × 10⁵   m/s, τ_A = 5.8   s, S = 2.7 × 10⁹, S_v = 2.09 × 10⁹, Pr = 1.52 × 10^-2, P_rad = 3.35 × 10^-6, $\mathit{\kappa }\mathrm{\times }{\mathit{T}}_{\mathrm{0}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{0.18}\hspace{0.17em}\mathrm{W}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$, 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₀ = 1.45  ×  10^-3, corresponding to 10³   m/s.

We impose as initial conditions a temperature profile with a base temperature T₀ = 10⁴   K and a top temperature 8 × 10⁵   K with the following z-dependence 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 length-scale 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 L0 = 4 × 106   m.

3. Results

With previous definitions, Eq. (4) can be replaced by the magnetic vector potential equation: $\frac{\mathit{\partial }{A}}{\mathit{\partial t}}\mathrm{=}{v}\mathrm{\times }\mathrm{(}{\mathit{B}}_{\mathrm{0}}\hspace{0.17em}{ê}\mathit{z}\mathrm{+}\mathrm{\nabla }\mathrm{\times }{A}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathit{S}} \mathrm{\nabla }\mathrm{\times }\mathrm{\nabla }\mathrm{\times }{A}\mathit{.}$(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 de-aliasing. Spatial derivatives are calculated in Fourier space, and nonlinear product terms are advanced in configuration space. A second-order central difference technique (Dahlburg et al. 1986) is used for the discretization in z. A time-step 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⁴ with a resolution of 128³. 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 volume-averaged 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 (${\mathit{e}}_{\mathit{v}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}\mathit{n}{{v}}^{\mathrm{2}}\mathrm{⟩}$ and ${\mathit{e}}_{\mathit{b}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}{{b}}^{\mathrm{2}}\mathrm{⟩}$), 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 (Tmax) and maximum electric current  | j | max as functions of time.

We verified that the dynamics of this problem are well represented by RMHD. We defined the volume-averaged 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₀   ê_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 self-consistent heating due to the magnetic field-line 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 self-consistent 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⁶   K, higher than the initial peak temperature of T = 8 × 10⁵   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 quasi-2D pancake-like structures aligned with the magnetic field. The two isosurfaces for T = 106   K and T = 1.2 × 106   K at time t = 300 are shown in yellow and gray (corresponding to T = 100,120 in nondimensional units), respectively. Sample magnetic field-lines are drawn in red. The box has been rescaled in z to clarify the visualization.

This expectation is confirmed by the three-dimensional snapshot of the isosurfaces of two temperatures (T = 10⁶ and T = 1.2 × 10⁶   K) at t = 300 in Fig. 3. The spatial temperature morphology shows that the high temperatures are reached in selected quasi-2D 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⁶   K at t = 300), while farther out the temperature decreases to a background value that in the central z-region averages about 6 × 10⁵   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 multi-thermal 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 sub-volume 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 multi-scale dynamics that can account for the observed emission.

The density stratification in our simulation allows limited evaporation from the z-boundary 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 z-boundaries, 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

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 L0 = 4 × 106   m.
In the text

	Fig. 2 Maximum temperature (Tmax) and maximum electric current  | j | max as functions of time.
In the text

	Fig. 3 Temperature is organized in quasi-2D pancake-like structures aligned with the magnetic field. The two isosurfaces for T = 106   K and T = 1.2 × 106   K at time t = 300 are shown in yellow and gray (corresponding to T = 100,120 in nondimensional units), respectively. Sample magnetic field-lines are drawn in red. The box has been rescaled in z to clarify the visualization.
In the text