SHEAR PHOTOSPHERIC FORCING AND THE ORIGIN OF TURBULENCE IN CORONAL LOOPS

In recent papers (Rappazzo et al. 2007, 2008), we have described reduced magnetohydrodynamic (RMHD) simulations of the Parker problem (Parker 1972, 1988, 1994) for coronal loops in Cartesian geometry. We have shown that the system develops small scales, organized in current sheets elongated in the direction of the DC magnetic field, through an MHD turbulent cascade, and that a well-defined power-law spectrum is developed for total energy. The energy spectra develop steep slopes (seen also in the similar simulations of a line-tied RMHD system by Dmitruk et al. 2003) with spectral indices going from the classical −5/3 Kolmogorov spectrum up to almost −3. The energy spectral indices (slopes) have a bearing on the heating rate of the boundary-forced system. The heating rates are significantly increased for realistic values of the magnetic field intensity and loop length compared to the scaling laws with fixed indices (Dmitruk & Gómez 1999), the last being recovered when a Kolmogorov-like spectrum develops (i.e., for weak fields or long loops).

In the published simulations, we have always used photospheric velocity patterns made up of large spatial scale projected convection cells (large-scale eddies), mimicking disordered photospheric motions. As a consequence, the magnetic field developing in the corona was not in equilibrium, but dynamical evolution occurred from the outset with a timescale set by the interplay of forcing and nonlinearity. In this paper, in order to clarify the origin of the MHD turbulent dynamics found in our previous work, we explore the dynamics of coronal loops system when a shear velocity flow stirs the footpoints of the magnetic field lines. It is commonly thought that the topology of the photospheric driver should strongly influence the dynamics of a coronal loop, and that the magnetic field lines anchored to the photospheric planes should passively follow their footpoints motions. In this picture, the electric currents should develop along neighboring field lines whose footpoints have a relative shear motion. In particular, it might be argued that in our previous simulations the turbulent dynamics of the coronal medium originated from the "complexity" of the photospheric forcing patterns, with their large-scale eddies and stagnation points.

The simple unidimensional shear forcing velocity used here allows a very clear-cut numerical experiment: if the field lines were to passively follow the imposed footpoint motion only a sheared magnetic field would develop inside the volume. The dynamics of these current layers are subject to tearing instabilities (Furth et al. 1963). If this were the physical process at work the topology of the magnetic field would remain a mapping of the forcing velocity pattern, periodically disrupted by tearing-like instabilities when the shear grew beyond a threshold amount. Simulations show that initially the magnetic field is sheared and a tearing instability develops, but afterward turbulent dynamics similar to those with more complex boundary patterns develop, clarifying that turbulence does not stem from a direct influence of the photospheric velocity pattern but it is due to the inherent nonlinear properties of the system.

The Parker Scenario for coronal heating has been the subject of intense research. Both analytical (Sturrock & Uchida 1981; van Ballegooijen 1986; Antiochos 1987; Berger 1991; Gómez & Ferro-Fontan 1992; Cowley et al. 1997; Ng & Bhattacharjee 1998; Uzdensky 2007; Janse & Low 2009; Low & Janse 2009; Berger & Asgari-Targhi 2009; Aly & Amari 2010) and numerical (Mikic et al. 1989; Longcope & Sudan 1994; Hendrix & Van Hoven 1996; Galsgaard & Nordlund 1996; Huang & Zweibel 2009; Huang et al. 2010) investigations have been carried out. Given the complexity of such boundary-forced field-line tangling systems simplified models have been developed, including two-dimensional incompressible MHD models with magnetic forcing (Einaudi et al. 1996; Dmitruk & Gómez 1997; Georgoulis et al. 1998; Dmitruk et al. 1998; Einaudi & Velli 1999) and shell models (Nigro et al. 2004; Buchlin & Velli 2007). Insights gleaned from these models include the important result that dissipation in forced MHD turbulence occurs in the form of "bursty" events with well-defined power-law distributions in energy release, peak dissipation, and duration, in a way reminiscent of, and consistent with, the distribution of flares in the solar corona.

An analytical model of a forced system very similar to the simulation presented here was proposed by Heyvaerts & Priest (1992), and recently extended to the anisotropic turbulence regime by Bigot et al. (2008). They started with same MHD system threaded by a strong axial magnetic field in Cartesian geometry and apply at the top and bottom boundaries two one-dimensional velocity fields of opposite direction and assumed that the sheared structure that develops in the corona then dissipates via an effective "turbulent resistivity" provided by a cascade, so that a dissipative equilibrium is set up in which shearing is balanced by slippage provided by the turbulence. This amounts essentially to a one-point closure model of MHD turbulence (Biskamp 2003), where turbulence acts only on very small scales while the large scales remain laminar and indeed with the same large-scale magnetic structure. They then use the eddy-damped quasi-normal Markovian approximation (EDQNM; Pouquet et al. 1976) to estimate the effective cascade and dissipation for the given driving shear, and this allows them to develop a heating theory in which the only free parameter is the equivalent Kolmogorov constant. In contrast, our simulations will show that nonlinearity cannot be neglected even at the large scales, so that the turbulent self-consistent state has little resemblance to the imposed shear flow. As a result, Heyvaerts & Priest (1992) overestimate the heating rate as laminar dynamics would lead to a higher energy injection (Poynting flux).

More recently, Dahlburg et al. (2005, 2009) have proposed the so-called secondary instability as a leading mechanism operating in the Parker Scenario, responsible for the rapid release of energy. In their view, disruption on ideal timescale must arise after some time while slow quasi-steady reconnection allows magnetic energy to continue to accumulate in the system. In their view, the system evolution may be described by a sequence of equilibria destabilized by magnetic reconnection. The bulk of numerical simulations performed by Einaudi et al. (1996), Dmitruk & Gómez (1997), Georgoulis et al. (1998), Dmitruk & Gómez (1999), Einaudi & Velli (1999), Dmitruk et al. (2003), and Rappazzo et al. (2007, 2008) has proven that the system does not evolve through a sequence of equilibria, rather more complex dynamics develop.

The initial setup of the simulation presented in Dahlburg et al. (2009) is also very similar to the one implemented in this paper but, besides the lower resolution and therefore the higher influence of numerical diffusion, the time interval for which they advance the equations is too short compared with coronal loops and active region timescales. These lead them to claim as representative of the dynamics what is actually a transient event taking place only during the early stage of the dynamics, in our case the multiple tearing mode current sheet collapse. This evolution is not generic, but is only representative of a small class of very symmetric boundary velocity patterns, those which admit coronal equilibria at all times. We will return to this question and a more detailed discussion in the conclusion.

The paper is organized as follows. In Section 2, we describe the basic governing equations and boundary conditions, as well as the numerical code used to integrate them. In Section 3, we discuss the initial conditions for our simulations and briefly summarize the linear stage dynamics more extensively detailed in Rappazzo et al. (2008), while in Section 4 we outline the main points of Heyvaerts & Priest (1992) relevant to this work. The results of our numerical simulations are presented in Section 5, while the final section is devoted to our conclusions and discussion of the impact of this work on coronal physics.

We model a coronal loop as an axially elongated Cartesian box with an orthogonal cross section of size ℓ and an axial length L embedded in an homogeneous and uniform axial magnetic field $\mathbf {B_0} = B_0\, \mathbf {\hat{e}_z}$ aligned along the z-direction. Any curvature effect is neglected.

The top and bottom plates (z = 0 and L) represent the photospheric surfaces where we impose, as boundary conditions, velocity patterns mimicking photospheric motions. Along the x- and y-directions periodic boundary conditions are implemented.

At the top plate z = L we impose a sinusoidal shear flow with wavenumber 4

$$\begin{equation} \mathbf {u^L} \left(x, y \right) = \sin \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y}. \end{equation} \tag{ 1 }$$

At the bottom plate z = 0 we generally impose a vanishing velocity

$$\begin{equation} \mathbf {u^0} \left(x, y \right) = 0, \end{equation} \tag{ 2 }$$

except in one simulation (run F (Table 1)) where, in order to compare with previous simulations implementing a vortical velocity forcing applied at both plates (Rappazzo et al. 2007, 2008), the reversed pattern of the top plate (Equation (1)) is applied

$$\begin{equation} \mathbf {u^0} \left(x, y \right) = {-} {\rm sin} \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y}. \end{equation} \tag{ 3 }$$

Here, $\mathbf {\hat{e}_y}$ is the unitary vector directed along the y-direction, while the flow is sheared along x.

The dynamics are integrated, as in our previous works, using the equations of RMHD (Kadomtsev & Pogutse 1974; Strauss 1976; Montgomery 1982), which are well suited for a plasma embedded in a strong axial magnetic field. In dimensionless form, they are given by

$$\begin{eqnarray} \frac{\partial \mathbf {u_{_\perp }}}{\partial t} + (\mathbf {u_{_\perp }} \cdot \mathbf {\nabla _{_{\!\perp }}}) \mathbf {u_{_\perp }} &=& {-} \mathbf {\nabla _{_{\!\perp }}} \left(p + \frac{ \mathbf {b_{_\perp }}^2 }{2} \right) + (\mathbf {b_{_\perp }} \cdot \mathbf {\nabla _{_{\!\perp }}}) \mathbf {b_{_\perp }} \nonumber \\ &&+ c_{A}\, \frac{\partial \mathbf {b_{_\perp }}}{\partial z} + \frac{(-1)^{n+1}}{{\rm Re}_n} \mathbf {\nabla _{_{\!\perp }}^{2n}}\, \mathbf {u_{_\perp }}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \frac{\partial \mathbf {b_{_\perp }}}{\partial t} + (\mathbf {u_{_\perp }} \cdot \mathbf {\nabla _{_{\!\perp }}}) \mathbf {b_{_\perp }} &=& (\mathbf {b_{_\perp }} \cdot \mathbf {\nabla _{_{\!\perp }}}) \mathbf {u_{_\perp }} + c_{A}\, \frac{\partial \mathbf {u_{_\perp }}}{\partial z} \nonumber \\ && + \frac{(-1)^{n+1}}{{\rm Re}_n} \mathbf {\nabla _{_{\!\perp }}^{2n}}\, \mathbf {b_{_\perp }}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\mathbf {\nabla _{_{\!\perp }}} \cdot \mathbf {u_{_\perp }} = 0, \qquad \mathbf {\nabla _{_{\!\perp }}} \cdot \mathbf {b_{_\perp }} = 0, \end{eqnarray} \tag{ 6 }$$

where $\mathbf {u_{_\perp }}$ and $\mathbf {b_{_\perp }}$ are the velocity and magnetic fields components orthogonal to the axial field, p is the kinetic pressure. The gradient operator has components only in the perpendicular x–y planes

$$\begin{equation} \mathbf {\nabla _{_{\!\perp }}} = \mathbf {\hat{e}_x} \frac{\partial }{\partial x} + \mathbf {\hat{e}_y} \frac{\partial }{\partial y} \end{equation} \tag{ 7 }$$

while the linear term ∝∂_z couples the planes along the axial direction through a wave-like propagation at the Alfvén speed c_A. Incompressibility in RMHD equations follows from the large value of the axial magnetic fields (Strauss 1976) and they remain valid also for low β systems (Zank & Matthaeus 1992; Bhattacharjee et al. 1998) such as the corona.

To render the equations nondimensional, we have first expressed the magnetic field as an Alfvén velocity [$b \rightarrow b/\sqrt{4\pi \rho _0}$], where ρ₀ is the density supposed homogeneous and constant, and then all velocities have been normalized to the velocity u* = 1 km s⁻¹, the order of magnitude of photospheric convective motions.

Lengths and times are expressed in units of the perpendicular length of the computational box ℓ* = ℓ and its related crossing time t* = ℓ*/u*. As a result, the linear terms ∝∂_z are multiplied by the dimensionless Alfvén velocity c_A = v_A/u*, where $v_A = B_0/\sqrt{4\pi \rho _0}$ is the Alfvén velocity associated with the axial magnetic field.

The majority of the simulations performed (specifically runs A–E (see Table 1)) use a standard simplified diffusion model, in which both the magnetic resistivity η and viscosity ν are constant and uniform. The kinetic and magnetic Reynolds numbers are then given by

$$\begin{equation} {\rm Re} = \frac{\rho _0\, \ell ^{\ast } u^{\ast }}{\nu }, \qquad {\rm Re}_{_m} = \frac{4\pi \, \ell ^{\ast } u^{\ast }}{\eta c^2}, \end{equation} \tag{ 8 }$$

where c is the speed of light, and numerically they are given the same value Re = Re_m. In Equations (4) and (5), this case is realized for n = 1 with Re₁ = Re.

The index n is called dissipativity and for n > 1 the dissipative terms in Equations (4) and (5) correspond to the so-called hyperdiffusion (Biskamp 2003). We use hyperdiffusion, with n = 4, only in runs F and G (Table 1) dedicated to study the energy spectra. Hyperdiffusion is used because, even with a grid of 512 × 512 points in the x–y plane (the highest resolution grid we used for the plane), the timescales associated with ordinary diffusion are small enough to affect the large-scale dynamics and render difficult the resolution of an inertial range. The diffusive time τ_n at the scale λ associated with the dissipative terms used in Equations (4) and (5) is given by

$$\begin{equation} \tau _{_n} \sim {\rm Re}_{_n}\, \lambda ^{2n}. \end{equation} \tag{ 9 }$$

For n = 1, the diffusive time decreases relatively slowly toward smaller scales, while for n = 4 it decreases far more rapidly. As a result, for n = 4 we have longer diffusive timescales at large spatial scales and diffusive timescales similar to the case with n = 1 at the resolution scale. Numerically, we require the diffusion time at the resolution scale λ_min = 1/N, where N is the number of grid points, to be of the same order of magnitude for both normal and hyperdiffusion, i.e.,

$$\begin{equation} \frac{{\rm Re}_{_1}}{N^2} \sim \frac{{\rm Re}_{_n}}{N^{2n}} \quad \longrightarrow \quad {\rm Re}_{_n} \sim {\rm Re}_{_1}\, N^{2(n-1)}. \end{equation} \tag{ 10 }$$

Then for a numerical grid with N = 512 points that requires a Reynolds number Re₁ = 800 with ordinary diffusion we can implement Re₄ ∼ 10¹⁹ (Table 1), removing diffusive effects at the large scales and allowing (if present) the resolution of an inertial range.

We solve numerically Equations (4)–(6) written in terms of the potentials of the orthogonal velocity and magnetic fields (see Rappazzo et al. (2007, 2008) for a more detailed description of the numerical code) in Fourier space, i.e., we advance the Fourier components in the x- and y-directions of the scalar potentials. Along the z-direction, no Fourier transform is performed so that we can impose non-periodic boundary conditions (Section 3), and a central second-order finite-difference scheme is used. In the x–y plane, a Fourier pseudospectral method is implemented. Time is discretized with a third-order Runge–Kutta method.

At time t = 0, we start our simulations with a uniform and homogeneous magnetic field along the axial direction $\mathbf {B} = B_0\, \mathbf {\hat{e}_z}$. The orthogonal component of the velocity and magnetic fields are zero inside our computational box $\mathbf {u_{_\perp }}=\mathbf {b_{_\perp }}=0$, while at the top and bottom planes a large-scale velocity pattern is imposed (Equations (1) and (2) or Equations (1)–(3)) and kept constant in time.

We briefly summarize and extend to the shear forcing considered in this paper the linear stage analysis covered in more detail in Rappazzo et al. (2008). In general for an initial interval of time smaller than the nonlinear timescale t < τ_nl, nonlinear terms in Equations (4)–(6) can be neglected and the equations linearized. For simplicity, we will at first neglect also the diffusive terms and consider their effect in the second part of this section. The solution during the linear stage for generic boundary velocity forcings, u^L and u⁰, respectively, at the top and bottom planes z = L and 0, is given by

$$\begin{eqnarray} &&\mathbf {b_{_\perp }} (x,y,z,t) = [ \mathbf {u^L} (x,y) - \mathbf {u^0} (x,y) ] \frac{t}{\tau _A},\qquad\ \ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&\mathbf {u_{_\perp }} (x,y,z,t) = \mathbf {u^L} (x,y)\, \frac{z}{L} + \mathbf {u^0} (x,y) \left(1 - \frac{z}{L} \right), \end{eqnarray} \tag{ 12 }$$

where τ_A = L/v_A is the Alfvén crossing time along the axial direction z. The magnetic field grows linearly in time, while the velocity field is stationary and the order of magnitude of its rms is determined by the boundary velocity profiles. Both are linear combinations (mapping) of the boundary velocity fields.

Considering the boundary conditions (1) and (2) imposed in most of the following simulations, we obtain

$$\begin{eqnarray} &&\mathbf {b_{_\perp }} (x,y,z,t) = \frac{t}{\tau _A}\, \sin \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y}, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\mathbf {u_{_\perp }} (x,y,z,t) = \frac{z}{L}\, \sin \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y}. \end{eqnarray} \tag{ 14 }$$

Both fields are a clear mapping of the shear velocity at the boundary, with the magnetic field increasing its magnitude linearly in time.

When we shear the field lines from both photospheric plates (z = 0 and L) using the forcing (1)–(3) we obtain a very similar result, still a mapping but with different amplitudes:

$$\begin{eqnarray} &&\mathbf {b_{_\perp }} (x,y,z,t) = 2\, \frac{t}{\tau _A}\, \sin \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y},\quad\ \ \ \ \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\mathbf {u_{_\perp }} (x,y,z,t) = \left(2\frac{z}{L} -1 \right) \, \sin \left(4\, \frac{2\pi }{\ell }\, x +1 \right) \, \mathbf {\hat{e}_y}. \end{eqnarray} \tag{ 16 }$$

For a generic forcing, the solution (11) and (12) is valid only during the linear stage, while for t > τ_nl when the fields are big enough the nonlinear terms cannot be neglected.

Nevertheless, there is a singular subset of velocity forcing patterns for which the generated coronal fields (11) and (12) have a vanishing Lorentz force and the nonlinear terms vanish exactly. This subset of patterns is characterize by having the vorticity constant along the streamlines (Rappazzo et al. 2008). In this case, the solutions (11) and (12) are an exact solution at all times, not an approximation valid only for t < τ_nl.

As can be proved by direct substitution the sheared forcing (1) and (2) (and also (1)–(3)) is one of these degenerate patterns, and the solution (13) and (14) (or (15) and (16)) is exact at all times.

So far we have neglected the diffusive terms in the RMHD Equations (4)–(6). Magnetic reconnection develops even for small values of the resistivity, but in this section we are interested in the diffusive effects on the linear dynamics, i.e., when nonlinear terms are negligible or artificially suppressed (we will discuss such a case in Section 6). We now consider the effect of standard diffusion (case n = 1 in Equations (4) and (5)) on the solutions (11)–(16): these are the solutions of the linearized equations obtained from Equations (4) and (5) retaining also the diffusive terms.

In the linear regime, as the magnetic field grows in time (Equations (11), (13), and (15)), the diffusive term $[\mathbf {\nabla _{_{\!\perp }}^2}\, \mathbf {b_{{_\perp }}} \! \propto \! \mathbf {b_{_\perp }}/\ell ^2]$ becomes increasingly bigger until diffusion balances the magnetic field growth, and the system reaches a saturated equilibrium state. Considering diffusion, the magnetic field will evolve as

$$\begin{eqnarray} \mathbf {b_{_\perp }} (x,y,z,t) &=& [ \mathbf {u^L} (x,y) - \mathbf {u^0} (x,y) ] \nonumber \\ && \times \, \frac{\tau _R}{\tau _A} \, \left[ 1 - \exp \left(- \frac{t}{\tau _R} \right) \right]. \end{eqnarray} \tag{ 17 }$$

The diffusive timescale associated with the Reynolds number Re is τ_R = ℓ²_c Re/(2π)² where ℓ_c is the length scale of the forcing pattern, that for the pattern (1) and (2) is given by ℓ_c = ℓ/4 where ℓ is the orthogonal computational box length.

The total magnetic energy E_M and Ohmic dissipation rate J will then be given by

$$\begin{eqnarray} &&E_M = \frac{1}{2}\, \int _V\! {d}^{^3}\hspace{-.4em} x \ \mathbf {b_{_\perp }}^2 = E_M^{\rm sat} \, \left[ 1 - \exp \left(- \frac{t}{\tau _R} \right) \right]^2,\quad \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&\hspace{-12pt}J = \frac{1}{R}\, \int _V\! {d}^{^3}\hspace{-.4em} x \ \mathbf {j}^2 = J^{\rm sat} \, \left[ 1 - \exp \left(- \frac{t}{\tau _R} \right) \right]^2, \end{eqnarray} \tag{ 19 }$$

where E^sat_M and J^sat are the saturations value reached for t ≳ 2 τ_R, whose values are given by

$$\begin{equation} E_M^{\rm sat} = \frac{\ell ^6 c_A^2 u_{\rm ph}^2 {\rm Re}^2}{2L(8\pi)^4}, \qquad J^{\rm sat} = \frac{\ell ^4 c_A^2 u_{\rm ph}^2 {\rm Re}}{L(8\pi)^2}. \end{equation} \tag{ 20 }$$

Magnetic energy saturates to a value proportional to the square of both the Reynolds number and the Alfvén velocity, while the heating rate saturates to a value that is proportional to the Reynolds number and the square of the axial Alfvén velocity. In general, the equations of RMHD are valid as far as the orthogonal magnetic field $\mathbf {b_{_\perp }}$ is small compared to the dominant axial field $\mathbf {B_0} = B_0\, \mathbf {\hat{e}_z}$. In particular incompressibility holds as far as the perturbed magnetic pressure can be neglected compared to that of the strong field $\mathbf {b^2_{_\perp }} \ll \mathbf {c^2_A}$. Therefore, the solutions found in this section are valid as far as the saturated values of the magnetic field satisfy the previous condition. In all the simulations presented here this condition is satisfied.

Given the complexity of the Parker problem, simplified models have been derived. Heyvaerts & Priest (1992) have developed an effective diffusivity model that is in effect a one-point closure model. In order to discuss the impact of our work on their results, we briefly summarize the relevant one-point closure theory (Biskamp 2003).

In order to investigate basic properties of the Parker model and compare them with observational constraints, such as the global heating rate and required energy flux, the detailed dynamics of turbulent fluctuations (that are essential to determine how the individual field lines are heated and hence how radiation is emitted) contain more informations than actually needed.

It may then be attempted to split the velocity and magnetic fields into average and fluctuating parts:

$$\begin{equation} \mathbf {B} =\, \langle \mathbf {B} \rangle +\, \mathbf {\widetilde{b}}, \qquad \mathbf {u} =\, \langle \mathbf {u} \rangle +\, \mathbf {\widetilde{u}}. \end{equation} \tag{ 21 }$$

Incompressible MHD equations give for mean fields the following equations:

$$\begin{eqnarray} \partial _t \langle \mathbf {u} \rangle + \langle \mathbf {u} \rangle \cdot \mathbf {\nabla } \langle \mathbf {u} \rangle &=& {-} \mathbf {\nabla } \langle P\rangle + \langle \mathbf {B} \rangle \cdot \mathbf {\nabla } \langle \mathbf {B} \rangle \nonumber \\ && \, {-} \mathbf {\nabla } \cdot \langle \mathbf {\widetilde{u}} \mathbf {\widetilde{u}} - \mathbf {\widetilde{b}} \mathbf {\widetilde{b}} \rangle + \nu \, \mathbf {\nabla ^2} \langle \mathbf {u} \rangle, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \partial _t \langle \mathbf {B} \rangle + \langle \mathbf {u} \rangle \cdot \mathbf {\nabla } \langle \mathbf {B} \rangle &=& \langle \mathbf {B} \rangle \cdot \mathbf {\nabla } \langle \mathbf {u} \rangle + \mathbf {\nabla } \times \langle \mathbf {\widetilde{u}} \times \mathbf {\widetilde{b}} \rangle + \eta \, \mathbf {\nabla ^2}\, \langle \mathbf {{B}} \rangle,\nonumber\\ \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} &&\mathbf {\nabla } \cdot \langle \mathbf {u} \rangle =\mathbf {\nabla } \cdot \mathbf {\widetilde{u}} = \mathbf {\nabla } \cdot \langle \mathbf {B} \rangle =\mathbf {\nabla } \cdot \mathbf {\widetilde{b}} = 0, \end{eqnarray} \tag{ 24 }$$

where symbols have the usual meaning, and in particular ν and η are the microscopic viscosity and resistivity of the plasma. If it is possible to model the terms that contain the small-scale fluctuations then Equations (22)–(24) allow to advance the mean fields.

This is a two-scale approach, and the average and fluctuating parts can be represented as the large-scale and small-scale fields or, introducing a suitable cutoff wavenumber K (for instance, the grid resolution if this is applied to numerical simulations), as the low- and high-pass filtered "lesser" and "greater" functions

$$\begin{eqnarray} \langle \mathbf {B} \rangle & = & \mathbf {B}^{<}_{K} = \sum _{k \le K} \mathbf {\widehat{B}_{k}}\, e^{i \mathbf {k} \cdot \mathbf {x}}, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \mathbf {\widetilde{b}} & = & \mathbf {B}^{>}_{K} = \sum _{k > K} \mathbf {\widehat{B}_{k}}\, e^{i \mathbf {k} \cdot \mathbf {x}}. \end{eqnarray} \tag{ 26 }$$

The relevant quantities to be modeled in Equations (22)–(24) are the "turbulent" stress tensors

$$\begin{equation} R_{ij} = - \langle \widetilde{u}_i \widetilde{u}_j - \widetilde{b}_i \widetilde{b}_j \rangle, \qquad S_{ij} = - \langle \widetilde{u}_i \widetilde{b}_j - \widetilde{u}_j \widetilde{b}_i \rangle. \end{equation} \tag{ 27 }$$

Phenomenological modeling and numerical simulations (Yoshizawa 1990, 1991) have shown that these stress tensor can be approximated with

$$\begin{eqnarray} &&\hspace{-5pt}R_{ij} \sim \nu _t (\partial _i \langle u_j \rangle + \partial _j \langle u_i \rangle), \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&S_{ij} \sim \eta _t (\partial _i \langle B_j \rangle - \partial _j \langle B_i \rangle), \end{eqnarray} \tag{ 29 }$$

where the coefficients ν_t and η_t are called turbulent viscosity and resistivity, and for a plasma in coronal conditions have much higher values than ν and η. Inserting Equations (28) and (29) in the equations for the mean fields (Equations (22) and (23)), we obtain a set of equations for the mean fields that has the same structure of incompressible MHD equations except that ν → ν_t and η → η_t.

Then for a system embedded in a strong axial field (Montgomery 1982) Equations (4)–(6) that we use to model the system are also a one-point closure model if the dissipative coefficients are considered effective ones.

Heyvaerts & Priest (1992) use the results of an EDQNM (Pouquet et al. 1976) to express the effective diffusivity coefficients ν_t and η_t as a function of the energy flux flowing along the inertial range.

Additionally, they suppose that the one-dimensional boundary forcing velocity leads the large scales to evolve in a laminar (fields are one-dimensional and directed along y, the same direction of the forcing) steady state (∂_t = 0). Equations (4)–(6) and their boundary forcing (of which our forcing (1) and (2) is representative) are used to compute the flux of energy S entering the system due to the dragging of the field-line footpoints.

S is a function of the system parameters and the effective diffusivity coefficients ν_t and η_t that are the only unknown variables in both energy fluxes S and . The solution of the problem is achieved requiring that the flux S entering the system at the large scales is equal to the flux flowing from the large to the small scales.

As a matter of fact their laminar steady state coincides with the saturation linear state that we have computed in the previous section (Section 3). This was obtained neglecting the nonlinear terms in Equations (4) and (5) and retaining the diffusive terms. They obtain it in a similar way, by supposing that the induced fields will retain the one-dimensional symmetry of the forcing. Thus, the nonlinear terms vanish, as can be proved by direct substitutions of the generic fields $\mathbf {b_{_\perp }} = f(x)\, \mathbf {\hat{e}_y}$, $\mathbf {u_{_\perp }} = g(x)\, \mathbf {\hat{e}_y}$, with f and g generic functions.

Our simulations investigate the large-scale dynamics and as we show in Section 5 the large-scale flow is not laminar and it is steady only in a statistical sense. Turbulence cannot be confined only to the small scales. We discuss the implications for the findings and scaling laws of Heyvaerts & Priest (1992) in section Section 6 devoted to our conclusions.

In this section, we present a series of numerical simulations, summarized in Table 1. We will first describe the results of simulations A–E that model with different resolutions and associated different Reynolds numbers a coronal layer driven by the sheared velocity pattern (1) at the top plate z = L and a vanishing velocity (2) at the bottom plate z = 0. In all simulations, the computational box has an aspect ratio of 10 with ℓ = 1 and L = 10.

5.1. Shear Forcing: Run A

We present here the results of run A, a simulation performed with a numerical grid of 512 × 512 × 200 points, normal (n = 1) diffusion with a Reynolds number Re = 800. The Alfvén velocity is v_A = 200 km s⁻¹ corresponding to a ratio c_A = v_A/u_ph = 200. The total duration is 600 axial Alfvén crossing times τ_A = L/v_A.

We add to the system a perturbation (naturally present in the coronal environment) for the magnetic and velocity fields inside the computational box at time t = 0. If there were no perturbation the system would follow the linear saturation curves plotted with dashed lines in Figures 1 and 2, as discussed in Section 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have not performed a full linear instability analysis, but we have used as perturbations either an Orszag–Tang vortex (Orszag & Tang 1979), i.e.,

$$\begin{eqnarray} &&\delta \mathbf {u} = \epsilon \, [ - \sin \left(2\pi y \right) \mathbf {\hat{e}_x} + \sin \left(2\pi x \right) \mathbf {\hat{e}_y} ], \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} &&\delta \mathbf {b} = \epsilon \, [ - \sin \left(4\pi y \right) \mathbf {\hat{e}_x} + \sin \left(2\pi x \right) \mathbf {\hat{e}_y} ], \end{eqnarray} \tag{ 31 }$$

or a "white noise" (i.e., the value of the potentials associated with the fields are given in each grid point a random number included between 0 and 1), with different amplitudes . We have computed the resulting growth rates γ with our nonlinear code.

For both perturbations, the value of the growth rates is high and similar to each other. The white noise has a slightly higher growth rate at γτ_A ∼ 0.69. This timescale is of the order of the Alfvén crossing time τ_A, which corresponds to an ideal timescale. This implies that the Orszag–Tang vortex is close to the most unstable mode, which is selected from all the modes that are excited by the random perturbation. As we shall see in the following, the instability is a multiple tearing mode.

For each kind of perturbation, the smaller the amplitude the later the system becomes unstable. For the simulation presented in this section, we have used a "white noise" perturbation with an amplitude = 10⁻¹⁶, very small compared to the boundary imposed velocity Equation (1) which is $\mathcal {O}(1)$.

In Figures 1 and 2, we plot the total magnetic and kinetic energies

$$\begin{equation} E_M = \frac{1}{2} \int \! {d}V\, \mathbf {b_{_\perp }}^2, \qquad E_K = \frac{1}{2} \int \! {d}V\, \mathbf {u_{_\perp }}^2, \end{equation} \tag{ 32 }$$

and the total Ohmic and viscous dissipation rates

$$\begin{equation} J = \frac{1}{{\rm Re}} \int \! {d}V\, \mathbf {j}^2, \qquad \Omega = \frac{1}{{\rm Re}} \int \! {d}V\, \mathbf {\omega }^2, \end{equation} \tag{ 33 }$$

along with the saturation curves (Equations (18) and (19)). For smaller values of the perturbation amplitude , the system becomes unstable sooner but the ensuing dynamics in the fully nonlinear stage are similar with same average values for all relevant physical quantities.

As shown in Figures 1–3, the shear velocity at the top boundary (z = 10) initially induces velocity and magnetic fields inside the volume that follow the linear behavior (Equation (12) and (17)). The magnetic energy and the Ohmic dissipation rate also follow initially the linear diffusive saturation curves (Equation (18) and (19)). From Equation (20) with the values for this simulation we have that the magnetic energy and Ohmic dissipation would reach the saturation values

$$\begin{equation} E_M^{\rm sat} = 1604, \qquad J^{\rm sat} = 2533, \end{equation} \tag{ 34 }$$

with a diffusion time τ_R ∼ 25 τ_A, if the system were unperturbed.

Zoom In Zoom Out Reset image size

Figure 3 shows snapshots of the magnetic field lines (orthogonal component $\mathbf {b_{_\perp }}$) and electric current at selected times in the midplane z = 5. Until time t ∼ 75 τ_A, the velocity and magnetic fields in the volume are a mapping of the velocity at the boundary (Equation (12) and (17)), i.e., the sheared velocity at the boundary induces a sheared magnetic field in the volume.

This magnetic configuration is well known to be unstable to tearing instabilities (Furth et al. 1963), and in fact around time t ∼ 79 τ_A a multiple tearing instability develops when the system transitions from the linear to the nonlinear stage. Magnetic field lines reconnect in correspondence of X points while the characteristic magnetic islands are formed. The randomness of the perturbations is reflected in the lack of symmetry of the reconnecting field.

The magnetic energy (Figure 1) accumulated up to this point is then released in a big burst of Ohmic dissipation (see Figure 2). Note that correspondingly there is a peak in the viscous dissipation (i.e., the integrated vorticity), but it is smaller than the electric current peak as in the following dynamical evolution. This is in fact a magnetically dominated system also when considering only the magnetic fluctuations $\mathbf {b_{_\perp }}$ with respect to velocity fluctuations $\mathbf {u_{_\perp }}$, without taking into account the big axial field B₀.

Up to this point, the dynamics are not surprising, the magnetic field gets sheared and the shear increases in time. Such a configuration is very well known to be unstable to reconnection, and then it is expected that a tearing-like instability should develop.

What is surprising are the dynamics at later times. In fact as can be seen in Figures 1 and 2 for t > 90 τ_A the system reaches a magnetically dominated statistically steady state where integrated quantities fluctuate around an average value. In particular velocity fluctuations are smaller than magnetic fluctuations, which in turn are small compared to the axial magnetic field $(\langle \mathbf {b}_{_\perp }^2 \rangle ^{1/2}/B_0 \sim 0.027)$. Therefore, the orthogonal magnetic and velocity fields satisfy the RMHD ordering that requires them to be small compared to the axial field B₀.

The energy power entering the system at the boundaries as a result of the work done by photospheric motions on the footpoints of magnetic field lines is given by the integrated Pointing flux

$$\begin{equation} S = c_A \! \int \limits _{z=L} \! \! {d}a\, \mathbf {b_{_\perp }} \cdot \mathbf {u_{_\perp }^L}, \end{equation} \tag{ 35 }$$

where $\mathbf {u_{_\perp }^L}$ is the photospheric forcing velocity (1).

Dissipation rates and the Poynting flux also fluctuate around a mean value. In particular, as shown in the inset in Figure 2, the Poynting flux and the total dissipation (Ohmic plus viscous) balance each other on the average, although on shorter timescales a lag between the two signals is present. In this steady state, the energy that is injected into the system is then, on the average, completely dissipated.

The surprising feature is that the shearing in the magnetic field inside the volume is not recreated, as shown in Figure 3. It is natural to think that after the first big dissipative event around t ∼ 82 τ_A, the shear forcing velocity at the boundary would recreate over time a sheared magnetic field in the system that should then lead to another big dissipative event and so on.

The dynamics of this system are in fact commonly approximated as a sequence of equilibria, each destabilized by magnetic reconnection. This approximation is attained by neglecting the velocity and kinetic pressure in the MHD equations whose solution is then bound to be a static force-free equilibrium. Our simulations confirm that the system is magnetically dominated, and in particular magnetic energy is bigger than kinetic energy, as shown in Figure 1 where on the average E_M ∼ 61 E_k. But the self-consistent evolution of the kinetic pressure and velocity, although small compared with the dominant axial magnetic field B₀, does not bind the system to force-free equilibria and allows the possible development of alternative dynamics.

In the following sections, we will analyze further aspects of the dynamics and the spectral properties of the system. But first we illustrate the topology of the magnetic and velocity fields, to understand why a sheared magnetic field is not recreated.

5.1.1. Magnetic Field Topology and the Origin of Turbulence

As shown in Figure 3 at time t ∼ 79 τ_A reconnection starts to develop, enhancing the Ohmic dissipation, that reaches a peak around t ∼ 82 τ_A. This big dissipative event burns a large fraction of the magnetic energy previously accumulated by the system (Figure 1). But not all the magnetic energy is dissipated, which would otherwise bring the system to a configuration similar to the initial condition at t = 0.

This dissipative event is in fact due to magnetic reconnection, that during its evolution produces a component of the magnetic field along x, the cross-shear direction, forming magnetic islands (see Figure 3 at times t ∼ 79 and 82 τ_A). In fact around t ∼ 90 τ_A, at the end of the big dissipative event, the topology of the orthogonal component of the magnetic field is characterized by magnetic islands. Naturally, the Lorentz force does not vanish now and the vorticity is not constant along the streamlines. As typical of magnetic reconnection vorticity forms misaligned quadrupolar structures around current sheets (see Rappazzo et al. 2008).

Although the forcing velocity at the boundary is always a shear (Equations (1) and (2)), nonlinear terms do not vanish as they do during the linear stage for t < 79 τ_A. When they vanish magnetic energy can be stored without getting dissipated (see Section 3). But now nonlinearity can redistribute along the cross-shear (x) direction part of the energy associated with the shear-aligned (y-oriented) field along which the forcing injects energy, and continuously cascades to lower scales as described in the following sections.

The three-dimensional structures are shown in Figure 4. Although the magnetic energy dominates over the kinetic energy, the ratio of the rms of the orthogonal magnetic field over the axial dominant field B₀ is quite small. For c_A = 200 it is ∼3%, so that the average inclination of the magnetic field lines with respect to the axial direction is just ∼2°, it is only for lower value of c_A that this ratio increases and the angle increases accordingly (Rappazzo et al. 2008). The field lines of the total magnetic field at time 550 τ_A are shown in Figure 4 (top row). The computational box has been rescaled for an improved viewing, and to attain the original aspect ratio, the box should be stretched 10 times along the axial direction. The magnetic topology for the total field is quite simple, as the lines appear slightly bent.

Zoom In Zoom Out Reset image size

Figure 4 (bottom row) also shows a view from the side and the top of the three-dimensional current sheets at time 550 τ_A. The current sheets, elongated along the axial direction, look space filling when watched from the side of the computational box, but the view from the top shows that the filling factor is actually small, as they are almost two-dimensional structures.

So far we have analyzed the topology of the field lines only in the midplane z = 5. In Figure 5, we show the current density and magnetic field lines of $\mathbf {b}_{_\perp }$ in the midplane (z = 5) and in two other x–y planes close to the boundaries (z = 1 and 9, the axial length L = 10). The behavior is similar at different heights although in the plane z = 9, closer to the forced boundary z = 10, the topology of the field appears to be affected to some extent by the sheared velocity forcing (1) directed along the y-direction. The field lines in fact show a small alignment directed along y close to the boundary z = 10.

Zoom In Zoom Out Reset image size

The influence of the boundary forcing over the magnetic field can be expressed quantitatively through the correlation between the magnetic field $\mathbf {b_{_\perp }}$ in the plane z and the boundary forcing velocity u^L (Equation (1)):

$$\begin{equation} \mathrm{Cor} [ \mathbf {b}_{_{\!\perp }}, \mathbf {u^L} ] \left(z \right) = \frac{ \int \! \!\! \int \! {d}x {d}y \, \mathbf {b_{_\perp }} \! \cdot \mathbf {u^L} }{ \big [ \!\int \! \!\! \int \! {d}x {d}y \, \mathbf {b^2_{_\perp }} \cdot \int \! \!\! \int \! {d}x {d}y \, (\mathbf {u^L})^2 \big ]^{\frac{1}{2}} }. \end{equation} \tag{ 36 }$$

In Figure 6, we plot the correlation as a function of the axial coordinate z at selected times. In the linear stage, until time slightly bigger than t = 76 τ_A whereafter the system transitions to the nonlinear stage (Figures 1–3), the magnetic field is a mapping of the boundary velocity therefore as expected the correlation is equal to 1. In fact for the simulation presented in this section (run A), for which we have imposed the shear velocity profile (1) at the top plate z = 10 and a vanishing velocity at the bottom plate z = 0, the magnetic field in the linear stage is given by Equation (13) (or Equation (17) with u⁰ = 0 including diffusion) therefore the correlation is 1 as b_⊥ is proportional to the boundary velocity u^L.

Zoom In Zoom Out Reset image size

Next as the system transitions to the nonlinear stage releasing most of the accumulated magnetic energy the correlation between the magnetic field and the boundary forcing velocity decreases swiftly, at a faster pace the farther from the forced boundary z = 10, as shown by the curves at times 82.19 ⩽ t/τ_A ⩽ 88.28.

The correlation during the fully nonlinear stage is shown with color lines at 10 selected times separated by Δ t = 40 τ_A in the interval 200 τ_A ⩽ t ⩽ 600 τ_A. The correlation vanishes near the bottom boundary and then grows almost linearly with z up to ∼0.6 at the top boundary. As expected the correlation is bigger near the forced boundary and fades toward the interior of the computational box. The magnetic field is overall weakly correlated with the forcing velocity, and at most reaches a mild correlation close to the forcing boundary, but it is never close to a strong correlation Cor = 1.

In Figure 7, we show the two-dimensional averages in the x–y planes of the magnetic and velocity fields and of the Ohmic dissipation j²/R plotted as a function of z at different times. The behavior is very similar to our previous simulations with different (vortical) forcing patterns (Rappazzo et al. 2008). These macroscopic quantities are smooth and present almost no variation along the axial direction. The velocity must approach its boundary values at z = 0 and 10, and in the volume grows to values higher than the boundary velocity u_L. In fact also the velocity inside the volume is not a mapping of the boundary forcing but develops self-consistently, in particular the plasma jets at reconnection locations contribute too to its average. The reason of the overall smooth behavior of these quantities is that every disturbance or gradient along the axial direction is smoothed out by the fast propagation of Alfvén waves along this direction; their propagation time τ_A is in fact the fastest timescale present (in particular faster than the nonlinear timescale τ_A < τ_nl), and then the system tends to be homogeneous along this direction.

Zoom In Zoom Out Reset image size

5.2. Dissipation versus Reynolds Number and Transition to Turbulence

Simulations A–E differ only for the value of the Reynolds number (and corresponding grid resolution), all other parameters are the same including the amplitude of the perturbations. In Figure 8, we plot the total dissipation for the five simulations as a function of time.

Zoom In Zoom Out Reset image size

The simulation described in the previous paragraph had Re = 800, after the first big dissipative peak its signal displays a complex temporal structure that arises from the underlying turbulent dynamics. Decreasing the value of the Reynolds number to Re = 400 the structure of the signal has a simpler structure, with an almost sinusoidal form of different amplitude and period ∼8τ_A in some time intervals.

The dotted lines represent the linear diffusive behavior described in Section 3. These solutions (Equation (19)) are obtained when the nonlinear terms can be neglected. This can happen either because they vanish exactly due to the symmetry of the forcing velocity (as discussed in Section 3) or equivalently when nonlinear terms are depleted by diffusion, i.e., when diffusion dominates the dynamics.

At Re = 200 diffusion affects the system substantially, the instability takes longer to develop, and afterward the actual signal and linear saturation curve differ little. There is no first big dissipative peak, for a long time up to t ∼ 550 τ_A the signal follows the linear saturation curve (19). Afterward it departs slightly displaying a sinusoidal behavior of very small amplitude.

At lower Reynolds number (Re = 100 and 10 are shown) diffusion dominates the dynamics, and no small scales is formed, while Ohmic dissipation and energy follow the curve (Equations (18) and (19)) describing the diffusive equilibrium that is formed, as nonlinear terms are completely depleted.

It is very interesting to notice that total dissipation for simulations A, B, and C with Re = 800, 400, and 200 substantially overlap each other. As typical with spectral numerical codes, which use Fourier transforms to compute derivatives, dissipation due to numerical implicit diffusion is very small. We have checked that the energy conservation equation, which can be derived from Equations (4)–(6), is numerically verified within a very small error: it is around 2% for the lower resolution simulations and decreases below 0.5% at higher resolutions (see Rappazzo 2006, Figure 5.5). Therefore, in the simulations presented here the dissipation rates are substantially due to the explicit dissipative terms present in Equations (4) and (5), while the diffusion due to the numerical schemes is negligible.

We had already seen this behavior in our previous simulations with vortex forcing (Rappazzo et al. 2008), but in that case the boundary velocity was slightly different in each simulation. In fact, it was built by a linear combination of Fourier modes with random amplitudes normalized to have the same rms, but the spatial pattern was different for every simulation as determined by the random amplitudes. In the simulations presented here, the forcing is exactly the same for each simulation, as it is simply a single Fourier mode (Equation (1)).

As the forcing is the same for each simulation the overlap is more evident, and makes stronger the claim that total dissipation is independent of the Reynolds number beyond a threshold. We conjectured this hypothesis (Rappazzo et al. 2008) because the energy flux entering the system (Poynting flux) and the energy transport from large to small scales due to the development of a turbulent dynamics are both independent from the Reynolds number, where the turbulent transport exhibits this property only for a sufficiently high value of the Reynolds number.

This unfortunately does not imply that the thermodynamical and radiative outcome is independent of the Reynolds number. In fact how field lines are heated strongly depends on the dynamics and properties around the single current sheet. These become thinner and thinner at higher Reynolds numbers and the overall dynamics more chaotic. Hence, the thermodynamical properties of the field lines that cross the current sheets (elongated along the axial direction) and get so heated impulsively are all to be explored.

An open question is whether the dissipation is independent of the Reynolds number also in the single current sheets, or their number and properties change to attain independence for the total dissipation. Research in this area is active (Loureiro et al. 2007, 2009; Lapenta 2008; Servidio et al. 2009, 2010; Samtaney et al. 2009; Cassak et al. 2009; Bhattacharjee et al. 2009; Huang & Bhattacharjee 2010) and has already shown that at relatively high Reynolds numbers' reconnection departs the classic Sweet–Parker scalings in the MHD regime. Furthermore, for a plasma in coronal conditions kinetic effects cannot be excluded a priori, in fact they might play an important role in dissipating energy through particle acceleration. Also the high value of the magnetic Prandtl number could have a bearing (Schekochihin et al. 2004). Nonetheless the total dissipation, integrated over the whole volume, is likely to not depend on the detailed small-scale dissipative mechanism. In fact the energy injection rate (Equation (37)) depends only on the large-scale fields and the transfer of this energy toward the small scales appears to be local (Rappazzo & Velli 2010), i.e., it is determined by the fields at neighboring large scales, thus making also the transfer energy rate toward the small scales independent of the Reynolds number (provided it is beyond a minimal threshold to have scale separation).

In Figure 9, we show a close-up of total dissipation for an equal interval of time Δ t = 600 τ_A on the same scale in order to highlight their temporal structures. At higher Reynolds numbers smaller time frequencies are present, a clear indication of a transition to turbulence (Frisch 1995). For Reynolds numbers lower than 100, the signal is completely flat following exactly the linear saturation curve (Equation (19)).

Zoom In Zoom Out Reset image size

5.3. Spectral Properties

In order to study the spectral properties of the system and compare them with those of previous simulations, we have performed a new simulation where the sheared forcing is applied on both the top and bottom plates (run F) with reversed direction (Equations (1) and (3)). We compare these results with those obtained from a previous simulation (run G), that has all the same parameters, except that on the two boundary planes a large-scale "vortex-like" velocity pattern with the same rms (〈u²〉 = 1/2) is applied. Both simulations have been performed with a numerical resolution of 512 × 512 × 200 grid points, and hyperdiffusion with dissipativaty n = 4 and R₄ = 10¹⁹.

In Figure 10, we plot the energy spectra obtained from runs F and G. They are very similar to each other. Both the velocity and magnetic fields develop inertial ranges following power laws, and overlap each other. For both runs, the spectral index of the kinetic spectrum (∼−0.5) is much smaller than that of the magnetic energy (∼−2.1) that is steeper than Kolmogorov (−5/3). Also the kinetic energy has a lower value than the magnetic energy, as already noticed for the integrated quantities (Figure 1).

Zoom In Zoom Out Reset image size

The energy that is injected into the system for unit time is the integrated Poynting flux

$$\begin{equation} S = c_A \int _{z=L} {d}a \big(\mathbf {u_{_\perp }^L} \cdot \mathbf {b_{_\perp }}) - c_A \int _{z=0} {d}a (\mathbf {u_{_\perp }^0} \cdot \mathbf {b_{_\perp }}\big), \end{equation} \tag{ 37 }$$

where $\mathbf {u_{_\perp }^L}$ and $\mathbf {u_{_\perp }^0}$ are the imposed velocity patterns at the top and bottom planes.

For run G, we excite all wavenumbers $3 \le n_{_\perp } \le 4$, while for run F we excite only one Fourier component as $\mathbf {u_{_\perp }^L} = - \mathbf {u_{_\perp }^0} = \sin \left(8\pi x + 1 \right) \mathbf {\hat{e}_y}$, i.e., we are injecting energy in the system only at $\mathbf {n_{in}} = 4 \cdot 2\pi \, \mathbf {\hat{e}_x}$, the wavenumber 4 along x. This can be noticed also in Figure 10, where the kinetic spectrum for run G at n = 3 is higher than for run F, as part of the energy is injected also at n = 3 in the vortical case.

The lower level for the kinetic spectrum is due to the boundary conditions that roughly set the value or the velocity at the injection wavenumbers inside the volume. In the simple linear case, this is given by Equation (12). In the shear case, we would have E_K(4) = 1/2 · V · 〈u²〉 = 2.5 (V = 10 is the volume) in the linear regime, and from Figure 10 we notice that also in the nonlinear regime E_K(4) ∼ 2.5. On the other hand, the magnetic field grows linearly in time (Equation (11)) until a balance is reached between the energy flux that is injected at this scale and the flux of energy flowing toward smaller scales through a turbulent cascade.

The magnetic energy spectra of the two simulations are slightly different at the large scales with $n_{_\perp } \le 5$. The large-scale dynamics is in fact slightly different in the two cases.

In the vortical case (run G), energy is injected in all modes with wavenumbers $3 \le n_{_\perp } \le 4$ that then cascades toward smaller scales. In the shear case (run F), energy is injected only at one wavenumber $\mathbf {n_{_\perp }} = (4, 0)$. We have already noticed in Section 5.1.1 that although we continue shearing the footpoints of the field lines with our one-dimensional forcing (Equation (1)) in the nonlinear stage the orthogonal magnetic field is organized in magnetic islands (Figure 3), so that it is no longer a mapping of the boundary velocity. Although energy is injected only in the wavenumber 4 along x, energy is then redistributed by the nonlinear terms also to modes with wavenumbers along y at the large scales, and a small inverse cascade is present as in run G. This is the basic mechanism by which magnetic islands are sustained throughout the simulation in the nonlinear stage.

In this paper, we have investigated the dynamics of the Parker problem for the heating of coronal loops when the footpoints of the magnetic field lines are stirred by a one-dimensional shear velocity pattern at the photosphere-mimicking boundary, and compared these results with those previously obtained when a more complex "vortex-like" velocity pattern was imposed (Rappazzo et al. 2008). This very simple forcing is ideal to investigate the origin of turbulence in coronal loops and the influence of the boundary velocity forcing on the dynamics of the system.

We will also compare our results with those of Heyvaerts & Priest (1992) and of the more recent simulations of Dahlburg et al. (2005, 2009).

In summary, the main results presented in this paper are the following.

• 1.  
  Initially, the sheared velocity forcing induces a sheared perpendicular magnetic field inside the volume. The resulting current layers are known to be unstable to tearing modes (Furth et al. 1963). In fact when the system transitions from the linear to the nonlinear stage it is due to a multiple tearing instability, as shown in Figure 3. But once the system has become fully nonlinear the dynamics are fundamentally different. As the nonlinear terms no longer vanish they now do transport energy from the large to the small scales where in correspondence of the X points nonlinear magnetic reconnection takes place, without going through a series of equilibria disrupted by tearing-like instabilities. Similarly to the case with disordered vortical boundary forcing velocities (Rappazzo et al. 2007, 2008) in the fully nonlinear stage, the system is highly dynamical and chaotic (and increasingly so at higher Reynolds numbers). For this, we do not observe secondary tearing of the current sheets as in two-dimensional high-resolution simulations of decaying MHD turbulence (Biskamp & Welter 1989), as now at the small scales fast turbulent dynamics take place.
• 2.  
  The dynamics of the Parker model do not depend strongly on the pattern of the velocity forcing that mimics photospheric motions, as far as they are constant in time (we defer the investigation of time-dependent boundary forcing to a future work). The shear forcing (Equation (1)) applied only at the top plate is a very simple and ordered one-dimensional forcing. We have shown that the resulting dynamics are very similar to those developed when a more complex and disordered "vortex-type" forcing velocity is applied. We conclude that the turbulent properties of the system are not induced by the complexity of the path that the footpoints follow. It is rather the system itself to be intrinsically turbulent, and turbulence develops as we continuously inject energy at the scale of photospheric motions (∼1000 km).
• 3.  
  The system reaches a statistically steady state where, although the footpoints of the field lines are continuously dragged by the forcing shear, this does not induce a sheared magnetic field in the computational box. In fact, the topology of the magnetic field is not a mapping of the forcing velocity field. Nonlinear interactions are able to redistribute the energy that is injected only at the wavenumber $\mathbf {n_{_\perp }} = 4\, \mathbf {\hat{e}_x}$ also to perpendicular wavenumbers and to smaller wavenumbers through an MHD turbulent cascade. In physical space, this corresponds to the magnetic field being organized in magnetic islands, to small-scale formation (current sheets elongated along the axial direction) and to magnetic reconnection taking place.
• 4.  
  Kinetic and magnetic energies develop an inertial range where spectra exhibit a power-law behavior. Fluctuating magnetic energy dominates over kinetic energy. Spectra and integrated quantities, like energies and total dissipative rates, have values similar to those obtained with a vortex-type forcing velocity. In particular, the total dissipation rate of the same system simulated with different Reynolds number appears to overlap each other beyond Re = 200, suggesting that this is independent of the Reynolds number beyond a threshold.

As shown in Figures 1–3 initially the system until time t ∼ 79 τ_A follows the linear curves (12) and (17). Up to this point, the shear velocity at the top boundary induces a sheared magnetic field in the volume. As discussed in Section 5.1, we have introduced a perturbation mimicking those naturally present in the corona. With no perturbation the system would relax over the resistive diffusive timescale τ_R (∼25 τ_A for run A) in a saturated diffusive equilibrium as described in Equations (12) and (17). While the simulation presented here used a very small amplitude for the perturbation ( = 10⁻¹⁶), we have performed shortest simulations with different values for the amplitude. As expected for higher values of , the instability develops sooner and for smaller values later, always following the linear curves until the instability transitions to the nonlinear stage. The more complete and systematic analysis of Romeou et al. (2004, 2009) in two-dimensional confirms this behavior.

Dahlburg et al. (2009) have performed a similar simulation with a lower resolution and with a fixed value for the perturbation and for a time interval that covers only the initial stage of our simulations. They in fact stop right after the first big dissipative peak that in our Figures 1–3 corresponds at t ∼ 100 τ_A.

Continuing the simulation, and using a higher numerical resolution, the system reaches a statistically steady state where magnetic energy consistently fluctuates around a mean value and the shear is not recreated in the topology of the orthogonal magnetic field. Their analysis is then limited to a transient event taking place only during the early stages of the dynamics, and that afterward does not repeat.

As shown in Rappazzo et al. (2008) during the linear stage the system is able to accumulate energy well beyond the average value maintained in the nonlinear stage only if the boundary forcing velocity satisfies the condition that its vorticity is constant along the streamlines. The sheared profiles used in this paper satisfy this condition as well the profile used by Dahlburg et al. (2009) (a linear combination of six sheared profiles).

These profiles are a very small subset of all the possible forcing profiles, and while they are very useful to get insight into the origin of turbulence in coronal loops they are not representative of the disordered photospheric motions, for which the strong stress buildup required for secondary instability to develop does not take place. The significance of their conclusions is then strongly diminished.

Furthermore, the Parker angle for this system cannot be defined as the relative angle between magnetic field lines at which the system becomes unstable. This is not a well-posed definition. In fact for given initial conditions the angle or equivalently the time (as the linear Equations (12) and (17) imply) at which the instability develops depends on the value of the amplitude of the perturbation that we add to the system. Depending on the value of the perturbation the Parker angle so defined is not unique.

On the other hand, in the fully nonlinear stage the average magnetic field line magnitude fluctuates around a mean value. It is then possible to give a unique value for the Parker angle, defined now as the average inclination of the magnetic field lines respect to the axial direction as done in Rappazzo et al. (2007, 2008), and as originally introduced by Parker (1988).

As summarized in Section 4, the one-point closure model developed by Heyvaerts & Priest (1992) splits the domain into large and small scales. They conjecture that the large-scale fields evolve into a stationary laminar regime, the field magnitudes determined by the effective diffusion coefficients. These laminar regimes correspond to our linear saturated diffusive regimes computed in Section 3. In Figure 8, the dotted lines show such diffusive curves for different values of the Reynolds numbers.

In their model, the large-scale fields computed in this way are used to obtain S, the energy flowing into the system for unit time at the boundary (the power) due do the work done by photospheric motions on the magnetic field-line footpoints. They also calculate, through an EDQNM approximation, the value of the spectral energy flux flowing along the inertial range at the small scales. Both S and are functions of the effective diffusion coefficients, and the solution of the problem results requiring balance between the two powers S = (S and have both the dimension of a power, i.e., energy over time, as S is the Poynting flux integrated over the boundary surface and is integrated over the whole volume as in Rappazzo et al. (2007)).

As shown in our simulations the large-scale fields are not laminar, and they are stationary only statistically. Nevertheless, it is useful to use Heyvaerts & Priest (1992) model in order to understand why it is not applicable. From Figure 8, we can estimate that the effective Reynolds number for which the diffusive regime dissipation matches the dissipation of the simulated system is R_eff = 150. Unfortunately, this values is too low for their model to work. In fact for R = 150 the dynamics are so diffusive that only a few modes of the order of the injection scale (∼1000 km) are not suppressed but only reduced in magnitude. Therefore, there is no flux of energy at the small scales = 0.

At a more fundamental level, the idea to split the domain into large and small scales does not work because nonlinearity cannot be confined only at the small scales. As shown by our simulations nonlinearity is at work at all scales, and unfortunately this fundamental aspect cannot be circumvented.

Finally, the use of RMHD equations is valid as far as the magnetic field fluctuations $\mathbf {b_{_\perp }}$ are small compared to the axial magnetic field B₀. This seems particularly apt to describe the dynamics of long-lived slender loops that apparently show no dynamics while shining bright at the resolution scale (∼800 km) of current state-of-the-art X-ray and EUV imagers on board Hinode and STEREO. Clearly, these results do not apply to highly dynamical active regions where dynamics cannot be modeled as fluctuations about an equilibrium configuration.

The series of simulations that we have performed proves that dragging the footpoints of magnetic field lines in the Parker problem quickly triggers nonlinear dynamics for small values of the orthogonal magnetic fields, and that these small magnetic field fluctuations are able to transport a considerable amount of energy toward the small scale with the overall energy flux ∼1.6 × 10⁶ erg cm⁻² s⁻¹ (Rappazzo et al. 2008) in the lower range of the observed constraint ∼10⁷ erg cm⁻² s⁻¹.

This prevents the orthogonal magnetic fluctuations to grow to an arbitrarily high value, self-consisting limiting the dynamics of the Parker problem to small fluctuations if the initial conditions are given by a uniform strong axial magnetic field.

We thank Russ Dahlburg (NRL) for useful discussions. A.F.R. was supported by the NASA Postdoctoral Program. This research was supported in part by the Jet Propulsion Laboratory, California Institute of Technology under contract with NASA, and in part by ASI contract n. I/015/07/0 Exploration of the Solar System. Financial support by the European Commission through the SOLAIRE Network (MRTN-CT-2006-035484) and by the Spanish Ministry of Research and Innovation through projects AYA2007-66502 and CSD2007-00050 is gratefully acknowledged. Simulations have been performed through the NASA Advanced Supercomputing SMD award 09-1112 and at CINECA (Italy). A.F.R. thanks the Leverhulme Trust International Network for Magnetized Plasma Turbulence for travel support.