© ESO, 2010

1. Introduction

The magnetic field at the visible surface of the Sun is known to take the form of bipolar regions. Above these magnetic concentrations the field continues in an arch-like fashion. These formations appear usually as twisted loop-like structures. These loops can be thought of as a continuation of more concentrated flux ropes in the bulk of the solar convection zone. However, this interpretation is problematic because we cannot be certain that the magnetic field in the Sun is generated in the form of flux ropes. Indeed, simulations of successful large-scale dynamos suggest that concentrated tube-like structures are more typical of the early kinematic stage, but in the final nonlinear stage the field becomes more space-filling (Brandenburg 2005; Käpylä et al. 2008). The idea that the dynamics of such tubes is governed by magnetic buoyancy is problematic too, because the solar convection zone is strongly stratified with concentrated downdrafts and broader upwellings. This leads to efficient downward pumping of magnetic field toward the bottom of the convection zone (Nordlund et al. 1992; Tobias et al. 1998). This downward pumping is generally found to dominate over magnetic buoyancy. The question then emerges whether magnetic buoyancy can still be invoked as the main mechanism for causing magnetic flux emergence at the solar surface. Another possible mechanism for the emergence of magnetic field at the solar surface might simply be the relaxation of a strongly twisted magnetic field in the bulk of the convection zone. Twisted magnetic fields are produced by a large-scale dynamo mechanism that is generally believed to be the motor of solar activity (Parker 1979). One such dynamo mechanism is the α effect that produces a large-scale poloidal magnetic field from a toroidal one. However, this mechanism is known to produce magnetic fields of opposite helicity (Seehafer 1996; Ji 1999). This magnetic helicity of opposite sign is an unwanted by-product, because it quenches the dynamo effect (Pouquet et al. 1976). A commonly discussed remedy is therefore to allow the helicity of small-scale field to leave the domain, possibly in the form of coronal mass ejections (Blackman & Brandenburg 2003).

In order to study the emergence of helical magnetic fields from a dynamo, we consider a model that combines a direct simulation of a turbulent large-scale dynamo with a simple treatment of the evolution of nearly force-free magnetic fields above the surface of the dynamo. An additional benefit of such a study is that it alleviates the need for adopting a boundary condition for the magnetic field at the top of the dynamo region. This is important, because it is known that the properties of the generated large-scale magnetic field strongly depend on boundary conditions. A common choice for the outer boundary condition is to assume that the magnetic field can be matched smoothly to a potential field. Such a condition is relatively easily implemented in calculations employing spherical harmonic functions (see, e.g., Krause & Rädler 1980). A more realistic boundary condition might be an extrapolation to a force-free magnetic field where the Lorentz force vanishes in the exterior. This means that the current density is proportional to the local magnetic field, but the constant of proportionality depends generally on the magnetic field itself. This renders the magnetic boundary condition nonlinear, which is therefore not easy to implement. Moreover, a perfectly force-free magnetic field is not completely realistic either. Instead, we know that above the solar surface, magnetic fields drive flares and coronal mass ejections through the Lorentz force. A more comprehensive approach would be to include in the calculations the exterior regions above the solar or stellar surface. This can be computationally prohibitive and a realistic treatment of that region may not even be necessary. It may therefore make sense to look for simplifying alternatives. One possibility is therefore to attempt an iteration toward a nearly force-free magnetic field such that the field can deviate from a force-free state locally in regions where the field cannot easily be made force-free. This could be done by solving the induction equation with an additional ambipolar diffusion term, which implies the presence of an effective velocity correction proportional to the local Lorentz force. This is sometimes called the magneto-frictional method and has been introduced by Yang et al. (1986) and Klimchuk & Sturrock (1992). In this approach the electromotive force attains not only a term in the direction of B, but also a term in the direction of J (Brandenburg & Zweibel 1994). The latter corresponds to a diffusion term, which explains the diffusive aspects of this effect. However, the resulting ambipolar diffusivity coefficient is proportional to and can locally become so large that the computational time step becomes significantly reduced. This is a typical problem of parabolic equations. A better method is therefore to turn the problem into a hyperbolic one and to solve an additional evolution equation for the velocity correction where the driving force is the Lorentz force. This approach is sometimes called the “force-free model” (FFM), even though the field in this model is never exactly force-free anywhere (Mikić et al. 1988; Ortolani & Schnack 1993). In the context of force-free magnetic field extrapolations this method is also known as the stress-and-relax method (Valori et al. 2005).

2. Equations for the force-free model

The equation for the velocity correction in the force-free model (FFM) is similar to the usual momentum equation, except that there is no pressure, gravity, or other driving forces on the right-hand side. Thus, we just have (1)where J × B is the Lorentz force, J = ∇ × B / μ₀ is the current density, μ₀ is the vacuum permeability, is the viscous force, and ρ is here treated as a constant corresponding to a prescribed density. Equation (1) is solved together with the induction equation. In order to preserve ∇·B = 0, we write B = ∇ × A in terms of the vector potential A and solve the induction equation in the form (2)where we have adopted the so-called resistive gauge in which the electrostatic potential is equal to  − η∇·A and the magnetic diffusivity η is assumed constant. The value of η will be given in terms of the magnetic Reynolds number, whose value will be specified below. No continuity equation for ρ is solved in this part of the domain, because there is no pressure gradient in the momentum equation.

In the following we couple such a model for the magnetic field above the solar photosphere to a simulation of a large-scale dynamo. In order to keep matters simple, we restrict ourselves to the case of an isothermal equation of state with constant sound speed c_s. Our goal is then to analyze the appearance of the resulting magnetic field above the surface of the dynamo and to study also the effects on the dynamo itself.

3. The model

The idea is to combine the evolution equations for the dynamo interior with those of the region above by simply turning off those terms that are not to be included in the upper part of the domain. We do this with error function profiles of the form (3)where w is the width of the transition. Thus, the momentum equation is assumed to take the form (4)where is the viscous force, is the traceless rate-of-strain tensor, commas denote partial differentiation, is the specific pseudo-enthalpy, c_s = const is the isothermal sound speed, and f is a forcing function that drives turbulence in the interior. The pseudo-enthalpy term emerges from the fact that for an isothermal equation of state the pressure is given by , so the pressure gradient force is given by . The continuity equation can either be written in terms of h(5)where we have inserted the θ_w(z) factor to terminate the evolution of h in the exterior, or in terms of ρ, (6)which serves the same purpose, but also preserves total mass. Most of the runs presented below are carried out using Eq. (5), but comparisons using Eq. (6) resulted in rather similar behavior.

The forcing function consists of random plane helical transversal waves with wavenumbers that lie in a band around an average forcing wavenumber k_f. These waves are maximally helical with ∇ × f ≈ k_ff, so the helicity is positive. This type of forcing was also adopted in Brandenburg (2001) and many other recent papers. The profile function θ_w(z) in front of the forcing term restricts the occurrence of turbulence mostly to the dynamo region, z < 0. The forcing amplitude is chosen such that the rms velocity in this region, u_rms, is about 4% of the sound speed.

We adopt non-dimensional units by measuring density in units of the initially constant density ρ₀, velocity in units of c_s, and length in units of , where k₁ = 2π / L_x is the minimal wavenumber in the x direction with an extent L_x. In most of the cases reported below, the vertical extent is L_z1 ≤ z ≤ L_z2 with L_z1 =  − L_x / 3 and L_z2 = 2L_x / 3. In a few cases we shall consider larger domains where the domain is two or four times larger in the z direction than in the horizontal directions. The extent of the domain in the y direction is L_y = L_x. We adopt periodic boundary conditions in the x and y directions. For the velocity we employ stress-free boundary conditions at top and bottom, i.e. (7)for the magnetic field we adopt perfect conductor boundary conditions at the bottom, which corresponds to (8)and vertical-field or pseudo-vacuum conditions at the top, i.e., (9)Note that no mass is allowed to escape at the top. Although this restriction does not seem to affect the results of our simulations significantly, it might be useful to adopt in future applications outflow boundary conditions instead.

Our model is characterized by several dimensionless parameters. Of particular importance is the magnetic Reynolds number, (10)where k_f is the wavenumber of the energy-carrying eddies. The ratio of viscosity to magnetic diffusivity is the magnetic Prandtl number, Pr_M = ν / η. In most of our simulations we use Pr_M = 1. The typical forcing wavenumber, expressed in units of the box wavenumber, k_f / k₁, is another important input parameter. In our simulations this value is 10. For the profile functions we take a transition width w with k₁w = 0.1 in most of the runs, but in some cases it is 0.2. The magnetic field is expressed in terms of the equipartition value, B_eq, where , and the average is taken over the turbulence zone. We measure time in non-dimensional units τ = tu_rmsk_f, which is the time normalized to the eddy turnover time of the turbulence. As initial condition we choose a hydrostatic state, U = 0, with constant density ρ = ρ₀, and the components of the magnetic vector potential are random white noise in space with Gaussian statistics and low amplitude (10^-4 below equipartition).

In this paper we present both direct numerical simulations and mean-field calculations. In both cases we use the Pencil Code¹, which is a modular high-order code (sixth order in space and third-order in time) for solving a large range of different partial differential equations.

4. Results

We begin by considering first hydrodynamic and hydromagnetic properties of the model. In the turbulence zone the velocity reaches quickly a statistically steady value, while in the exterior it takes about 1000 turnover times before a statistically steady state is reached. This is seen in Fig. 1, where we show u_rms(z) at different times. In the following we discuss the properties of the magnetic field that is generated by the turbulence.

4.1. Dynamo saturation

Dynamo action is possible when Re_M reaches a critical value that is about 0.5 in our case. The situation is only slightly modified compared with dynamo saturation in a periodic domain. For not too large values of Re_M the dynamo saturates relatively swiftly, but for larger values of Re_M the magnetic field strength may decline with increasing value of Re_M (Brandenburg & Subramanian 2005). An example of the saturation behavior is shown in Fig. 2, where Re_M ≈ 3.4. The initial saturation phase (100 ≤ τ ≤ 500) is suggestive of the resistively slow saturation found for periodic domains (Brandenburg 2001), but then the field declines somewhat. Such a decline is not normally seen in periodic domains, but is typical of dynamo action in domains with open boundaries or an external halo (see Fig. 5 of Hubbard & Brandenburg 2010a). The field strength is about 78% of the equipartition field strength, B_eq.

Fig. 1 Vertical dependence of the rms velocity at different times. Note the drop of urms(z) from the turbulence zone to the exterior by a factor of about 3 − 5. The inset shows urms(z) in the exterior at different times.

Fig. 2 Initial exponential growth and subsequent saturation behavior of the magnetic field in the interior for forced turbulence with dynamo action.

Fig. 3 Magnetic field structure in the dynamo exterior at τ = 1601. Field lines are shown in red and the modulus of the current density is shown in pink with semi-transparent opacity. Note the formation of a vertical current sheet above the arcade.

Fig. 4 Structure of magnetic field lines in the exterior, together with a representation of the normal component of the field at the interface at z = 0 between turbulence zone and exterior at τ = 1601. Green represents a positive and blue a negative value of Bz.

In all cases the magnetic field is strongest in the turbulence zone, but it always shows a systematic variation in one of the two horizontal directions. It is a matter of chance whether this variation is in the x or in the y direction. Comparison of different runs shows that both directions are about equally possible (see below). Also, the magnetic field pattern shows sometimes a slow horizontal migration, but this too seems to be a matter of chance, as will be discussed below.

4.2. Arcade formation

After some time the magnetic field extends well into the exterior regions where it tends to produce an arcade-like structure, as seen in Figs. 3 and 4. The arcade opens up in the middle above the line where the vertical field component vanishes at the surface. This leads to the formation of anti-aligned field lines with a current sheet in the middle; see Figs. 3 and 4. The dynamical evolution is seen clearly in a sequence of field line images in Fig. 5 where anti-aligned vertical field lines reconnect above the neutral line and form a closed arch with plasmoid ejection above. This arch then changes its connectivity at the foot points in the sideways direction (here the y direction), making the field lines bulge upward to produce a new reconnection site with anti-aligned field lines some distance above the surface. Note that this sideways motion takes the form of a slowly propagating wave. However, it is a matter of chance whether this wave propagates in the positive or negative coordinate direction, as will be shown below in Sect. 4.4.

Field line reconnection is best seen for two-dimensional magnetic fields, because it is then possible to compute a flux function whose contours correspond to field lines in the corresponding plane. In the present case the magnetic field varies only little in the x direction, so it makes sense to visualize the field averaged in the x direction. Since the averaging commutes with the curl operator, we can also average the x component of the magnetic vector potential, i.e. we compute  ⟨ A_x ⟩ _x, where the second subscript indicates averaging along the x direction. This function corresponds then to the flux function of the magnetic field in the yz plane and averaged along the x direction. In Fig. 6 we plot contours of  ⟨ A_x ⟩ _x, which correspond to poloidal field lines of  ⟨ B ⟩ _x in the yz plane. This figure shows clearly the recurrent reconnection events with subsequent plasmoid ejection. We also compare with a color/grey scale representation of the x component of the x-averaged magnetic field,  ⟨ B_x ⟩ _x. Note that in the exterior the contours of  ⟨ B_x ⟩ _x trace nearly perfectly those of  ⟨ A_x ⟩ _x. This suggests that the x-averaged magnetic field has nearly maximal magnetic helicity in the exterior. This is also in agreement with other indicators that will be considered below.

4.3. Averaged field properties

The magnetic field is largely confined to the turbulence zone where it shows a periodic, nearly sinusoidal variation in the y direction. Away from the turbulence zone the field falls off, as can be seen from the vertical slice shown in Fig. 7. Near the top boundary, some components of the field become stronger again, but this is probably an artifact of the vertical field condition employed in this particular case.

In order to describe the vertical variation of the magnetic field in an effective manner, it is appropriate to Fourier-decompose the field in the two horizontal directions and to define a complex-valued mean field as (11)where superscripts l and m indicate a suitable Fourier mode. In Fig. 8 we plot absolute values of the three components of as a function of z for a time representing the final saturated state. This figure shows quite clearly a relatively rapid decline of with z, while and level off at values that are still about 40% of that in the turbulence zone. This suggests that our model is suitable to describe the evolution of magnetic fields in the dynamo exterior. Earlier simulations of coronal loops and coronal heating (Gudiksen & Nordlund 2002, 2005; Peter et al. 2004) demonstrate that the dynamics of such fields is controlled by the velocity properties at their footpoints, which is here the interface between the turbulence zone and the dynamo exterior.

Fig. 5 Time series of arcade formation and decay. Field lines are colored by their local field strength which increases from pink to green. The plane shows Bz increasing from red (positive) to pink (negative). The normalized time τ is giving in each panel.

Fig. 6 Time series of the formation of a plasmoid ejection. Contours of  ⟨ Ax ⟩ x are shown together with a color-scale representation of  ⟨ Bx ⟩ x; dark blue stands for negative and red for positive values. The contours of  ⟨ Ax ⟩ x correspond to field lines of  ⟨ B ⟩ x in the yz plane. The dotted horizontal lines show the location of the surface at z = 0.

4.4. Force-free versus current free

Already the straightforward inspection of magnetic field lines viewed from the top suggests that the magnetic field is twisted and forms a left-handed spiral; see Fig. 9. This is indeed the orientation expected for turbulence with positive kinetic helicity, producing a negative α effect and hence magnetic spirals with negative helicity at the scale of the domain.

We expect the magnetic field in the dynamo exterior to be nearly force free, i.e., we expect to be small compared with . Here,  ⟨ · ⟩ _H denotes an xy average. In order to characterize the degree to which this is the case, we define the quantities (12)and note that (13)with . In Fig. 10 we show and as functions of z. Given that has values below 0.1, it is evident that the field is indeed nearly force free in the exterior. In the turbulence zone, on the other hand, the Lorentz force is quite significant.

To prove the existence of a force-free structure, it is convenient to calculate the angle χ between J and B. We expect χ to be 0 or π for an ideal force-free environment. In Fig. 11 we show the distribution of values of χ plotted over the height z. One sees that, in the exterior, the angle χ is close to π, but it drops to π / 2 at the very upper part, where the normal field condition enforces that J and B are at right angles to each other.

Fig. 7 Slice of Bz through an arbitrarily chosen cross-section x = const. Note the periodicity with nearly sinusoidal variation in the y direction, and the more nearly monotonous fall-off in the z direction.

Fig. 8 Vertical dependence of the moduli of the components of , as defined in Eq. (11).

Fig. 9 Three-dimensional visualization of the magnetic field viewed from above. The vertical magnetic field component is color-coded (yellow pointing upward and blue pointing downward). Note that the field lines form a left-handed spiral over the scale of the domain, as expected for turbulence with positive helicity at small scales.

Fig. 10 Vertical dependence of and . Note the decline of the normalized Lorentz force from more than 60% in the turbulence zone to less than 10% in the exterior.

Fig. 11 Two-dimensional histogram of the distribution of angles, p(χ), where χ = arccos(J·B) is the angle between J and B at different heights. p(χ) is normalized such that ∫p(χ)dχ = 1. The dashed line gives the location of the maximum position of the distribution.

Fig. 12 Dependence of versus time τ and height z for Lz = 6.4 with ReM = 3.4 (Run A).

Fig. 13 Similar to Fig. 12, but for Lz = 8π and ReM = 6.7 (Run B1).

In order to demonstrate that plasmoid ejection is a recurrent phenomenon, it is convenient to look at the evolution of the ratio versus τ and z. This is done in Fig. 12 for L_z = 6.4 and Re_M = 3.4 (Run A) and in Fig. 13 for L_z = 8π and Re_M = 6.7 (Run B1). It turns out that in both cases the typical speed of plasmoid ejecta is about 1/2 of the rms velocity of the turbulence in the interior region. However, the time interval δτ between plasmoid ejections increases from  ≈ 250 to  ≈ 570 turnover times as we increase the kinetic and magnetic Reynolds numbers. At higher magnetic Reynolds numbers, the length of the interval increases to  ≈ 800 turnover times. A summary of all runs is given in Table 1. Here we also give the l and m values of the leading mode of the mean field in Eq. (11) and indicate explicitly whether the large-scale magnetic field varies in the x or the y direction. Indeed, both directions are possible, confirming that it is a matter of chance. In most cases the magnetic field pattern shows a slow horizontal migration, whose direction appears to be random. The sign in the table indicates whether the wave migrates in the positive or negative coordinate direction.

4.5. Interpretation in terms of a mean-field model

The magnetic field found in the simulations displays a clear large-scale structure. One may have expected that the magnetic field varies mainly in the z direction, because this is also the direction in which the various profile functions vary. However, this is not the case. Instead, the main variation is in one of the horizontal directions (see Fig. 7, where the field varies mainly in the y direction). The magnetic field does of course also vary in the z direction, but this happens without sign change in . Above the surface at z = 0, the field gradually decays and retains only rather smooth variations.

In order to compare with dynamo theory, we now solve the usual set of mean-field equations for the mean magnetic vector potential , where is the mean magnetic field and is the mean current density, (14)We recall that η = const is the microscopic magnetic diffusivity, which is not negligible, even though it is usually much smaller than η_t. We consider first the kinematic regime where α and η_t are independent of . In order to account for the fact that there is no turbulence above the turbulence zone, we adopt the profile θ_w(z) for α and η_t, i.e., we write (15)where α₀ and η_t0 are constants, and for w we take the same value as for the other profile functions used in the direct simulations. The excitation condition for the dynamo can be quantified in terms of a dynamo number that we define here as (16)where η_T0 = η_t0 + η is the total magnetic diffusivity and k₁ was defined in Sect. 3 as the smallest horizontal wavenumber that fits into the domain. If the turbulence zone were homogeneous and periodic in the z direction, the critical value of C_α is unity, but now the domain is open in the z direction, so one expects the dynamo to be harder to excite. In the models presented below we therefore adopt the value C_α = 2.5, which is also compatible with estimates of the critical value from the simulations, if we write α ≈ u_rms/3 and η_t ≈ u_rms / 3k_f.

Next, we consider the nonlinear regime by employing the dynamical quenching model (Kleeorin & Ruzmaikin 1982; Blackman & Brandenburg 2002). We assume that α = α_K + α_M, where now α_K(z) = α₀θ_w(z) is the kinetic α effect profile used earlier in the solution to the kinematic equations, and α_M is the solution to the dynamical quenching equation, (17)where is the mean electromotive force. Note that we have here ignored the possibility of magnetic helicity fluxes that must become important at larger values of Re_M.

Fig. 14 Comparing the average in the x direction of the magnetic field with the 2D mean-field model. The structures are very similar.

In Fig. 14 we compare the resulting magnetic field geometry with that of the direct simulations. In both cases the horizontal variation of the field is similar. However, in the direct simulations the field extends more freely into the exterior. This is probably caused by a vertical upward flow that appears to be driven by the magnetic field. In Fig. 15 we see the vertical dependence of an upward flow, which soon reaches a statistically steady state.

Fig. 15 Horizontally averaged rms velocity field as a function of height z for 5 different times for a run with Lz = 8π. Note the development of a statistically steady state after about 1000 turnover times.

5. Conclusions

Our first results are promising in that the dynamics of the magnetic field in the exterior is indeed found to mimic open boundary conditions at the interface between the turbulence zone and the exterior at z = 0. In particular, it turns out that a twisted magnetic field generated by a helical dynamo beneath the surface is able to produce flux emergence in ways that are reminiscent of that found in the Sun.

Some of the important questions that still remain open include the presence and magnitude of magnetic helicity fluxes. In the present model we expect there to be diffusive magnetic helicity fluxes associated with the vertical gradient of magnetic helicity density. A related question concerns the dependence on the magnetic Reynolds number. One expects magnetic helicity fluxes to become more important at large values of Re_M.

One of the future extensions of this model includes the addition of shear. In that case one might expect there to be strong magnetic helicity fluxes associated with the Vishniac & Cho (2001) mechanism that may transports magnetic helicity along the lines of constant shear, although more recent considerations now cast doubt on this possibility (Hubbard & Brandenburg 2010b). One should also keep in mind that the magnetic field cannot really be expected to be fully helical, as was assumed here in order to promote large-scale dynamo action under relatively simple conditions. Reducing the degree of helicity makes the dynamo harder to excite. On the other hand, shear helps to lower the excitation conditions, making it again feasible to obtain large-scale dynamo action even at low relative helicity of the driving. Another promising extension would be to move to a more global geometry, including the effects of curvature and gravity. This would allow for the emergence of a Parker-like wind that turns into a supersonic flow at sufficiently large radii. This would also facilitate the removal of magnetic field through the sonic point.

Acknowledgments

We thank Jaime de la Cruz Rodriguez, Gustavo Guerrero, and Göran Scharmer for discussions at the early stages of this work. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and the National Supercomputer Centers in Linköping. This work was supported in part by the European Research Council under the AstroDyn Research Project 227952 and the Swedish Research Council grant 621-2007-4064.

References

All Tables

All Figures

	Fig. 1 Vertical dependence of the rms velocity at different times. Note the drop of urms(z) from the turbulence zone to the exterior by a factor of about 3 − 5. The inset shows urms(z) in the exterior at different times.
In the text

	Fig. 2 Initial exponential growth and subsequent saturation behavior of the magnetic field in the interior for forced turbulence with dynamo action.
In the text

	Fig. 3 Magnetic field structure in the dynamo exterior at τ = 1601. Field lines are shown in red and the modulus of the current density is shown in pink with semi-transparent opacity. Note the formation of a vertical current sheet above the arcade.
In the text

	Fig. 4 Structure of magnetic field lines in the exterior, together with a representation of the normal component of the field at the interface at z = 0 between turbulence zone and exterior at τ = 1601. Green represents a positive and blue a negative value of Bz.
In the text

	Fig. 5 Time series of arcade formation and decay. Field lines are colored by their local field strength which increases from pink to green. The plane shows Bz increasing from red (positive) to pink (negative). The normalized time τ is giving in each panel.
In the text

	Fig. 6 Time series of the formation of a plasmoid ejection. Contours of  ⟨ Ax ⟩ x are shown together with a color-scale representation of  ⟨ Bx ⟩ x; dark blue stands for negative and red for positive values. The contours of  ⟨ Ax ⟩ x correspond to field lines of  ⟨ B ⟩ x in the yz plane. The dotted horizontal lines show the location of the surface at z = 0.
In the text

	Fig. 7 Slice of Bz through an arbitrarily chosen cross-section x = const. Note the periodicity with nearly sinusoidal variation in the y direction, and the more nearly monotonous fall-off in the z direction.
In the text

	Fig. 8 Vertical dependence of the moduli of the components of , as defined in Eq. (11).
In the text

	Fig. 9 Three-dimensional visualization of the magnetic field viewed from above. The vertical magnetic field component is color-coded (yellow pointing upward and blue pointing downward). Note that the field lines form a left-handed spiral over the scale of the domain, as expected for turbulence with positive helicity at small scales.
In the text

	Fig. 10 Vertical dependence of and . Note the decline of the normalized Lorentz force from more than 60% in the turbulence zone to less than 10% in the exterior.
In the text

	Fig. 11 Two-dimensional histogram of the distribution of angles, p(χ), where χ = arccos(J·B) is the angle between J and B at different heights. p(χ) is normalized such that ∫p(χ)dχ = 1. The dashed line gives the location of the maximum position of the distribution.
In the text

	Fig. 12 Dependence of versus time τ and height z for Lz = 6.4 with ReM = 3.4 (Run A).
In the text

	Fig. 13 Similar to Fig. 12, but for Lz = 8π and ReM = 6.7 (Run B1).
In the text

	Fig. 14 Comparing the average in the x direction of the magnetic field with the 2D mean-field model. The structures are very similar.
In the text

	Fig. 15 Horizontally averaged rms velocity field as a function of height z for 5 different times for a run with Lz = 8π. Note the development of a statistically steady state after about 1000 turnover times.
In the text