{'fraction_non_alphanumeric': 0.025714605150374537, 'fraction_numerical': 0.009838631535795476, 'mean_word_length': 4.339037007560685, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Surface appearance of dynamo-generated large-scale fields

Twisted magnetic fields are frequently seen to emerge above the visible surface of the Sun. This emergence is usually associated with the rise of buoyant magnetic flux structures. Here we ask how magnetic fields from a turbulent large-scale dynamo appear above the surface if there is no magnetic buoyancy. The computational domain is split into two parts. In the lower part, which we refer to as the turbulence zone, the flow is driven by an assumed helical forcing function leading to dynamo action. Above this region, which we refer to as the exterior, a nearly force-free magnetic field is computed at each time step using the stress-and-relax method. Twisted arcade-like field structures are found to emerge in the exterior above the turbulence zone. Strong current sheets tend to form above the neutral line, where the vertical field component vanishes. Time series of the magnetic field structure show recurrent plasmoid ejections. The degree to which the exterior field is force free is estimated as the ratio of the dot product of current density and magnetic field strength to their respective rms values. This ratio reaches values of up to 95% in the exterior. A weak outward flow is driven by the residual Lorentz force.


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 mag-netic fields of opposite helicity (Seehafer 1996;Ji 1999). This magnetic helicity of opposite sign is an unwanted byproduct, 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 forcefree 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 B 2 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).

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 DU where J × B is the Lorentz force, J = ∇ × B/µ 0 is the current density, µ 0 is the vacuum permeability, F visc 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 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 largescale 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.

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 where w is the width of the transition. Thus, the momentum equation is assumed to take the form where is the traceless rate-of-strain tensor, commas denote partial differentiation, h = c 2 s ln ρ 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 p = c 2 s ρ, so the pressure gradient force is given by ρ −1 ∇p = c 2 s ∇ ln ρ = ∇h. The continuity equation can either be written in terms of h where we have inserted the θ w (z) factor to terminate the evolution of h in the exterior, or in terms of ρ, which serves the same purpose, but also preserves total mass. Most of the runs presented below are carried out using Equation (5), but comparisons using Equation (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 f f , 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 ρ 0 , velocity in units of c s , and length in units of k −1 1 , where k 1 = 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.
for the magnetic field we adopt perfect conductor boundary conditions at the bottom, which corresponds to and vertical-field or pseudo-vacuum conditions at the top, i.e., 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, 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 all our simulations we use Pr M = 1. The typical forcing wavenumber, expressed in units of the box wavenumber, k f /k 1 , is another important input parameter. In our simulations this value is 10.
For the profile functions we take a transition width w with k 1 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 B 2 eq = µ 0 ρu 2 , and the average is taken over the turbulence zone. We measure time in nondimensional units τ = tu rms k 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 ρ = ρ 0 , 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 1 , 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.

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 Figure 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.   Figure 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 . 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.

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 Figures 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 Figures 3 and 4. The dynamical evolution is seen clearly in a sequence of field line images in Figure 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 sur-  face. 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 Section 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 Figure 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.

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 Figure 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 where superscripts l and m indicate a suitable Fourier mode. In Figure 8 we plot absolute values of the three components of B 01 (z, t) as a function of z for a time representing the final saturated state. This figure shows quite clearly a relatively rapid decline of |B y | with z, while |B x | and |B z | 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, 2005Peter 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. 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.

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 Figure 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 (J × B) 2 H to be small compared with B 2 H J 2 H . Here, · H denotes an xy average. In order to characterize the degree to which this is the case, we define the quantities and note that with k 2 JB = µ 2 0 J 2 H / B 2 H . In Figure 10 we show k 2 J×B and k 2 J·B as functions of z. Given that k 2 J×B /k 2 JB 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.  J×B /k 2 JB and k 2 J·B /k 2 JB . Note the decline of the normalized Lorentz force from more than 60% in the turbulence zone to less than 10% in the exterior.
In Figure 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, 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.
where the normal field condition enforces that J and B are at right angles to each other.
In order to demonstrate that plasmoid ejection is a recurrent phenomenon, it is convenient to look at the evolution of the ratio J ·B H / B 2 H versus t and z. This is done in Figure 12 for L z = 6.4 and Re M = 3.4 (Run A) and in  Table 1. Summary of runs. Run A is for a cubic domain while Runs B1 to B7 are for taller domains at different values of Re M . Note that the nondimensional interval length δτ as well as the ejection speeds cannot be determined accurately for most of the runs, but the values suggests that there is no systematic dependence on the value of Re M . The column 'lm' gives the l and m values of the leading mode of the mean field in Equation (11) and the columns 'direction' gives the direction of propagation of this mean field, confirming that it is a matter of chance.  Figure 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 Equation (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.

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 Figure 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 B z . 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 A, where B = ∇×A is the mean magnetic field and J = ∇ × B/µ 0 is the mean current density, 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 B. 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 where α 0 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 where η T0 = η t0 + η is the total magnetic diffusivity and k 1 was defined in Section 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) = α 0 θ 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, where E = αB − η t µ 0 J 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 . In Figure 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 Figure 15 we see the vertical dependence of an upward flow, which soon reaches a statistically steady state.

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.