REALISTIC MODELING OF LOCAL DYNAMO PROCESSES ON THE SUN

The origin of magnetic field generation is a key problem for understanding solar variability across a wide range of scales. Modern high-resolution observations of global magnetic fields by the Helioseismic and Magnetic Imager (HMI) on NASA's Solar Dynamics Observatory (Scherrer et al. 2012), as well as investigations of small-scale magnetic fields in selected areas using the Solar Optical Telescope on Hinode (Tsuneta et al. 2008) and the balloon-borne SUNRISE observatory (Solanki et al. 2010), have revealed complicated dynamics of magnetic fields on the solar surface and their tendency to self-organize into compact magnetic structures.

Traditionally, the solar dynamo problem is divided according to two timescales: the global dynamo, operating on the scale of the 22 year solar cycle and controlling the global toroidal and poloidal fields (e.g., Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Brandenburg et al. 2012), and the local dynamo (e.g., Nordlund et al. 1992; Cattaneo 1999; Cattaneo et al. 2003; Vögler & Schüssler 2007; Kitiashvili et al. 2013a; Rempel 2014), which operates on the timescales of granulation and super-granulation from a few minutes to a few hours and is believed to be responsible for the Sun's "magnetic carpet" (Schrijver et al. 1998; Schrijver & Title 2002). However, the energetic and dynamic connections between the two dynamos are not clear.

In the global dynamo, the kinetic helicity is produced by the action of the Coriolis force on convective turbulence (for recent formulations of the global dynamo theory, see, e.g., Brandenburg et al. 2005; Pipin 2008; Pipin & Kosovichev 2014). However, in the local dynamo problem, without global rotation, the mean helicity integrated over a volume is zero. This has caused a debate about the existence of a "fluctuation dynamo" driven by non-helical forced turbulence. Recent numerical simulations convincingly demonstrate the existence of such a dynamo (Cattaneo 1999; Schekochihin et al. 2004b; Iskakov et al. 2007). The exact physical mechanism of the fluctuating dynamo is not yet understood, and current theoretical interpretations include analytical Kazantsev (1968)–type models, models with turbulent anisotropy induced by a mean velocity shear ("shear-current" effect; Rogachevskii & Kleeorin 2003), and other ideas (for a review, see Schekochihin et al. 2007).

It is important to emphasize that the local dynamo mechanism in realistic solar conditions is different from both these cases. While the mean kinetic helicity on local dynamo scales is very small (zero in our model), strong stratification and highly anisotropic near-surface turbulence lead to strong local helical and shearing flows, particularly in downdrafts located in the intergranular lanes. Our simulations show that these flows are the primary source of the local dynamo in quiet-Sun regions.

Observational techniques are now able to resolve flows and magnetic fields associated with magnetic flux emergence events on very small subgranular scales (e.g., Centeno et al. 2007; Orozco Suárez et al. 2008; Martínez González & Bellot Rubio 2009; Thornton & Parnell 2011; Vargas Domínguez et al. 2014), which probably reflect the local dynamo at work in shallow subsurface layers. In addition, resolving magnetic features though observations on smaller and smaller scales has renewed interest in the small-scale magnetic field generation problem and raised a question about the existence of unresolved (or "hidden") magnetic flux in the quiet Sun (e.g., Stenflo 1982, 2012; Trujillo Bueno et al. 2004; Shchukina & Trujillo Bueno 2011).

Theories of small-scale (or turbulent) dynamo processes have been intensively studied in the past. In particular, the similarity between the MHD induction and vorticity equations (Batchelor 1950) initiated discussions about magnetic field amplification due to turbulent flows. By analogy with vortices when stretching makes them stronger, twisting and stretching of magnetic field lines can amplify the field strength. Thus, repeating stretching, twisting, and folding of magnetic field lines can provide a local increase in the magnetic flux and its recirculation in a dynamo process (Vainshtein & Ruzmaikin 1972).

Previous studies were mostly done with direct numerical simulations (DNSs) performed for artificially forced turbulent flows (e.g., Meneguzzi & Pouquet 1989; Schekochihin et al. 2004b, 2005) and also for convectively driven flows for a wide range of parameters (e.g., Nordlund et al. 1992; Cattaneo 1999; Brandenburg et al. 2012). In addition, recent "realistic"-type radiative MHD simulations have reproduced the solar surface conditions with a high degree of realism and demonstrated that magnetic fields can be quickly amplified by local dynamos in the upper convection zone from a very weak (∼10⁻² G) seed field, leading to more than 1 kG magnetic elements (Vögler & Schüssler 2007; Kitiashvili et al. 2013a; Rempel 2014).

In this paper, we present new results of 3D radiative MHD simulations of the local solar dynamo for various seed-field strengths, from 10⁻⁶ to 10⁻² G (in five different simulation runs), and we investigate the development and properties of the dynamo process in the realistic regime of solar magnetoconvection, when the turbulent magnetic Prandtl number is not prescribed a priori as a constant. Our goal is to investigate the primary physical mechanism of the magnetic field generation by local dynamo in the realistic solar conditions. A detailed statistical study and comparison are planned for a future paper.

Numerical simulations of local dynamos on the Sun have two specific features that must be taken into account: (1) the solar turbulence is driven by near-surface convective flows in a highly stratified medium; (2) the hydrodynamic and magnetic Reynolds numbers are so high that DNSs with fully resolved turbulent scales are not possible.

The first feature requires accurate modeling of solar conditions including radiative transfer with realistic opacities and a realistic equation of state. This means that DNSs of flows with artificially forced turbulence have limited applicability to the solar dynamo problem. For our simulations, we use a radiative MHD code ("SolarBox") specifically developed for modeling solar conditions with a high degree of realism. The basic formulation is presented in Section 2.1.

The second feature of solar simulations requires implementation of subgrid scale (SGS) turbulence models. Previous solar magnetoconvection simulations (e.g., Stein et al. 2003; Vögler & Schüssler 2007; Rempel 2014) used either artificial or numerical viscosity to model, explicitly or implicitly, dissipation due to turbulence at unresolved scales. In this work, we use a large-eddy simulation (LES) approach implemented in the SolarBox code. This approach is widely used in engineering applications and tested with laboratory experiments and is still a topic of extensive study. The SGS models adopted in our solar simulations are developed and tested by using DNS MHD simulations including a fluctuating dynamo driven by forced turbulence. This approach provides a verification of the LES models. Some initial results of this project are described by Balarac et al. (2010), and the LES models that currently provide the best performance (G. Balarac 2014, private communication) are described in Section 2.1.

Perhaps the most significant implication of the LES approach for dynamo modeling is that in our simulations the characteristic turbulent dissipation parameters are not prescribed as constants, but are modeled according to the local flow and magnetic field properties. This means that the critical non-dimensional parameters, such as the turbulent magnetic Prandtl number, are not prescribed but change with local conditions.

In Section 3, for the first time, we describe the physical mechanism of individual magnetic field generation events in terms of local plasma dynamics. We find a strong correlation between such events and local helical flows and vortex tubes. It is commonly believed that the flow helicity induced by rotation is responsible for the generation of the solar global (poloidal) magnetic field. No rotation is included in our simulations, but, because of the strong stratification of the turbulent surface convection, vortex tubes with helical motion are generated, mostly in high-speed (near-sonic) downdrafts in the intergranular lanes. The mechanisms of vortex tube generation, dynamics, and associated atmospheric effects were described in our previous papers (Kitiashvili et al. 2012a, 2012b, 2013a, 2013b). We show here that the local turbulent flow helicity in the solar magnetoconvection regime plays a key role in the small-scale dynamo process.

In Section 4, we present our results on the statistical relationship between dynamo properties and the local turbulent magnetic Prandtl number. Such investigations have not been performed before. Also, previous investigations were limited to descriptions of stationary dynamo states.

In Section 5, we investigate the effects of dynamo-generated magnetic fields on the turbulent properties of solar magnetoconvection (back reaction of magnetic fields). A surprising result is that the resulting kinetic energy spectrum is closer to the Kolmogorov law than the initial spectrum without magnetic field. We discuss the evolution of the kinetic and magnetic energy spectra from seed-field initiation to the dynamo stationary state and also the variation of the spectra with depth. This allows us to investigate the process of energy transfer in the dynamo mechanism.

In Section 6, we describe the development of local bipolar magnetic structures that probably play an important role in the observed "salt-and-pepper" structure of the quiet-Sun magnetograms. Another significant new result is that the local dynamo can produce relatively long-lived compact magnetic structures. Their formation is related to the amplification of shearing and twisting flows caused by growing electric currents, and thus magnetic forces, around magnetic elements.

In Section 7, we discuss the topological properties of dynamo-generated magnetic fields in the solar atmosphere and show that they form a canopy-like structure consisting of small-scale magnetic loops. Such a topology results in a rapid change with height in the solar atmosphere of the ratio between the vertical and horizontal field components and can explain the discrepancies in the observed ratios (e.g., Danilovic et al. 2010). The mean strength of the dynamo-generated magnetic fields corresponds to the observed values from Hinode (Lites et al. 2008). In Section 8, we summarize the principal results and discuss their implications for our understanding of solar magnetism, as well as future plans.

2.1. 3D Radiative MHD "SolarBox" Code

To perform the local dynamo simulations we use the 3D radiative MHD code "SolarBox" (Wray et al. 2015). The code is based on an LES formulation for compressible flow and includes a fully coupled radiation solver in which local thermodynamic equilibrium is assumed. Radiative transfer between fluid elements is calculated using a 3D multi-spectral-bin method with long characteristics. For the initial conditions, we use a standard solar model of the interior structure and lower atmosphere. The code has been carefully tested and was previously used for studying the excitation of solar acoustic oscillations by turbulent convection in the upper convection zone (Jacoutot et al. 2008b, 2008a) and several other problems of solar dynamics and MHD (e.g., Kitiashvili et al. 2010, 2011, 2013b).

We solve the grid-cell-averaged equations of the conservation of mass (1), momentum (2), energy (3), and magnetic flux (4):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{(\rho {u}_{i})}_{,i}=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {u}_{i}}{\partial t}+{(\rho {u}_{i}{u}_{j}+({P}_{{ij}}+\rho {\tau }_{{ij}}))}_{,j}=-\rho {\varphi }_{,i},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+[{{Eu}}_{i}+({P}_{{ij}}+\rho {\tau }_{{ij}}){u}_{j}-(\kappa +{\kappa }^{t}){T}_{,i}\\ &&\quad +{{(\displaystyle \frac{c}{4\pi })}^{2}\displaystyle \frac{1}{\sigma +{\sigma }^{t}}({B}_{i,j}-{B}_{j,i}){B}_{j}+{F}_{i}^{\mathrm{rad}}]}_{,i}=0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{i}}{\partial t}+{[{u}_{j}{B}_{i}-{u}_{i}{B}_{j}-\displaystyle \frac{{c}^{2}}{4\pi (\sigma +{\sigma }^{t})}({B}_{i,j}-{B}_{j,i})]}_{,j}=0,\end{eqnarray} \tag{ 4 }$$

where ρ is the averaged mass density, u_i is the Favre-averaged (density-weighted) velocity, B_i is the magnetic field, and E is the averaged total energy density $E=\frac{1}{2}\rho {u}_{i}{u}_{i}+\rho e+\rho \varphi +\frac{1}{8\pi }{B}_{i}{B}_{i}$, where φ is the gravitational potential and e is the Favre-averaged internal energy density per unit mass. F_i^rad is the radiative flux, which is calculated by solving the radiative transfer equation, and P_ij is the averaged stress tensor ${P}_{{ij}}=(p+\displaystyle \frac{2}{3}\mu {u}_{k,k}+\frac{1}{8\pi }{B}_{k}{B}_{k}){\delta }_{{ij}}-\mu ({u}_{i,j}+{u}_{j,i})-\frac{1}{4\pi }{B}_{i}{B}_{j}$, where μ is the dynamic viscosity. Gas pressure p is a function of e and ρ calculated through a tabulated equation of state (Rogers et al. 1996); τ_ij is the Reynolds stress, κ is the molecular thermal conductivity; κ^t is the turbulent thermal conductivity, σ is the molecular electrical conductivity, and ${\sigma }^{t}={c}^{2}/4\pi {\eta }_{m}^{t}$ is the turbulent electrical conductivity in Gaussian units. The magnetic diffusivity, ${\eta }_{m}^{t}$, is calculated using the subgrid-scale (SGS) model (Theobald et al. 1994; Balarac et al. 2010): ${\eta }_{m}^{t}=a{\triangle }^{2}| {\rm{\nabla }}\times \bar{B}| /\sqrt{\rho }$, where a = 2, △ is the grid spacing, and overbars denote the resolved field. The turbulent heat conductivity κ^t is given in terms of the turbulent thermal Prandtl number, Pr^t, as ${\kappa }^{t}=\rho {c}_{p}{\nu }^{t}/{{Pr}}^{t}$, where ${\nu }^{t}={\mu }^{t}/\rho$. We assume that Pr^t is unity.

Currently, it is impossible to achieve a realistic Reynolds number in the DNSs of solar MHD phenomena; therefore, modeling of the dynamical properties of solar convection is performed through implementation of SGS LES turbulence models. These models create a realistic turbulent Reynolds number at which the simulations can be accurate, and so they provide a better representation of small-scale motions (e.g., Zhou et al. 2014). Here, we used a Smagorinsky eddy–viscosity model (Smagorinsky 1963), in which the compressible Reynolds stresses were calculated in the form (Moin et al. 1991) ${\tau }_{{ij}}=-2{C}_{{\rm{S}}}{\triangle }^{2}| S| ({S}_{i,j}-{u}_{k,k}{\delta }_{{ij}}/3)$ + $2{C}_{C}{\triangle }^{2}| S{| }^{2}{\delta }_{{ij}}/3$, where the Smagorinsky coefficients are C_S = C_C = 0.001; S_ij is the large-scale stress tensor.

The simulation results presented here were obtained for a computational domain of 6.4 × 6.4 × 6.2 Mm, which includes 1 Mm of the lower atmosphere. The grid size is 12.5 km in the horizontal direction. Above the solar surface, the vertical grid size is constant (12 km) but stretched out in the deeper layers. The lateral boundary conditions are periodic. The top boundary is open for mass, momentum, and energy fluxes, and also for the radiation flux. The top boundary is handled through characteristic boundary conditions; hence, fluxes of all quantities can occur through this boundary in either direction. The bottom boundary is handled as follows: (1) the vertical velocity component is set to zero; (2) horizontally uniform numerical fluxes are created at the bottom boundary to compensate for any net inward or outward flux through the top boundary, for mass, momentum, energy, and B-field, except that the radiative energy flux through the top is not counted in this process; and (3) an energy flux equal to the stellar luminosity is introduced at the bottom boundary, also uniformly in the horizontal dimensions. As a result, all conserved quantities retain constant integrals over the computational volume, except that the total energy is allowed to vary due to radiative loss and the steady bottom-boundary flux that mimics energy flow from the core. This seems to us to be a physically natural treatment of the energy, and it also serves as a check on the radiative transfer calculation. We have verified that over a sufficiently long time interval, the constant power input at the bottom boundary comes to be compensated by the radiative flux through the top boundary, resulting in a total energy that oscillates randomly and with very small amplitude around a constant value. In summary, (1) mass, momentum, energy, and B-field can advect and diffuse through the top boundary; (2) these losses or gains are compensated at the bottom boundary; and (3) power to satisfy the stellar luminosity is introduced at the bottom, and radiation through the top quickly comes to balance this input within a few percent and eventually balances it exactly in time average.

Compared to early magnetoconvection simulations by Stein et al. (2003) and later work by Vögler & Schüssler (2007) and Rempel (2014), our simulation domain is significantly deeper: 5.2 Mm versus 2.5 Mm of Stein et al. (2003), 0.9 Mm of Vögler & Schüssler (2007), and 2.3 Mm of Rempel (2014) in his high-resolution runs. The deeper domain reduces the effects of the bottom-boundary conditions. The simulation results presented in Section 3 show that the primary region of magnetic field generation and energy exchange is located in the upper 1 Mm layer, but slowly extends into deeper layers. Therefore, it is important that the computational domain has a significant depth.

2.2. Initial Conditions

Magnetic field amplification in turbulent solar convection can roughly be divided into three basic mechanisms: (1) magnetic field concentration due to converging flows (e.g., Nordlund 1983), (2) convective collapse (Parker 1978; Spruit & Zweibel 1979), and (3) dynamo processes driven by helical or shearing motions (e.g., Batchelor 1950; Vainshtein & Ruzmaikin 1972; Brandenburg 1995). The turbulent nature of the photospheric layers, where the intense radiative cooling drives downward convective motions, makes it impossible to separate different sources of the locally growing magnetic energy (Spruit 1984), and therefore the magnetic field amplification can be discussed only in terms of a dominant mechanism.

To avoid confusion in separating dynamo and non-dynamo field amplification, no additional magnetic flux was introduced into the computational domain after the "seed"-field initialization. This provides an important test for the dynamo action: if there were no dynamo process, then local magnetic patches would only be formed due to converging flows in the intergranular lanes, and then convective downdrafts would transport these patches into deeper layers, from which a part of the flux would be recycled by turbulent motions. In contrast, with small-scale dynamo action, magnetic field patches are continuously formed, resulting in a rapid growth of magnetic energy.

We note that during very early times of the "dynamo" runs (shorter than one overturning time) the evolution of the "seed"-field magnetic elements behaves similarly to corks in distributed turbulent flows: they tend to collect in the intergranular lanes where the convective flows converge. Figure 1 shows the distribution of the vertical velocity (panel (a)) and magnetic field (panel (b)) in the photosphere shortly after introducing the "seed" magnetic field as a checkerboard structure (case B; see Table 1). Because the weak field mostly follows the turbulent flow, the regular field distribution becomes deformed. During this stage, local magnetic field amplification is primarily caused through compression by converging flows in the intergranular lanes. Shortly after this initial phase, we see the start of the dynamo process: small-scale helical motions in the intergranular lanes drag magnetic field lines and stretch them and can even reverse the initial local polarity of the field. The dynamo process becomes stronger with time as turbulent flows on larger scales become involved. The mean magnetic energy density becomes saturated after 7–8 solar hours. Snapshots of the vertical velocity and magnetic field in the developed state at z = 0 are shown in Figures 1(c) and (d). The maximum magnetic field magnitude reaches 2 kG, and the photospheric mean magnitude reaches ∼20 G (Figure 2(c)).

Zoom In Zoom Out Reset image size

The dynamics of the small-scale dynamo is illustrated in the animation of Figure 1, in which the semi-transparent color scale corresponds to the magnetic field strength. The color scale is saturated at 500 G, but the field strength can reach more than 2000 G at the photosphere. After the seed-field initialization, the first signatures of the amplified field (light green-blue diffuse structures) appear above the photosphere due to the short dynamical timescales there. Shortly thereafter, the first magnetic patches appear at the photosphere. The appearance of opposite-polarity patches reveals a "salt-and-pepper" distribution of magnetic fields (Figure 1(d)). The patches evolve rather chaotically, but nevertheless the maximum magnetic field strength grows continuously until it finally saturates. These magnetic elements are highly twisted and show a very dynamic and non-stationary behavior, and none of them survive for long.

Because the local dynamo is driven by turbulent motions that cause magnetic field amplification due to the stretching and twisting of field lines, it is natural to consider the relationship between turbulent plasma motions and the local evolution of the magnetic field. For instance, a relationship between the dynamo process and vortex motions produced by the Coriolis force in a rotating stratified medium was demonstrated numerically by Nordlund et al. (1992). Our simulations reveal small-scale vortex tubes in the intergranular lanes; however, no external rotation has been introduced. These vortex tubes are a result of overturning convection and Kelvin–Helmholtz type instabilities (Kitiashvili et al. 2012b). They play a significant role in the small-scale dynamo mechanism.

To quantify the apparent close relation between the appearance of magnetic patches and local twisting motions observed in our simulations, we calculate the kinetic helicity density, $H={\boldsymbol{v}}\cdot ({\rm{\nabla }}\times {\boldsymbol{v}})$, and investigate its correlation with the time derivative of the magnetic energy density. A similar relationship between the kinetic helicity density and local amplification of magnetic flux (analogous to the shear-twist-fold dynamo of Vainshtein & Zel'dovich 1972) was previously noticed by Brandenburg et al. (1996). In their simulations, the vorticity was induced by rotating the simulation box. In our case without global rotation, the mechanism of vorticity generation due to shearing instabilities and compression has been discussed in our previous paper (Kitiashvili et al. 2012a).

Vertical profiles of the rms velocity and the kinetic helicity density, $H={\boldsymbol{v}}\cdot ({\rm{\nabla }}\times {\boldsymbol{v}})$, displayed in Figure 3, show that the strongest turbulent flows occupy a relatively thin 1 Mm deep layer just below the solar surface. Despite the similarity of the vertical profiles for velocity and helicity, the kinetic helicity density is more concentrated in the upper 1 Mm subsurface layer. Our previous studies of helical motions (vortex tubes) in the near-surface layers (Kitiashvili et al. 2012b) showed that the lifetime of the vertically oriented vortex tubes can be much longer (up to 1 hr) than the granulation lifetime (∼10 minutes). The vortex tubes are characterized by strong flows reaching sonic speeds, by a sharp decrease of the gas pressure, and by intense radiative cooling. The strong horizontal and vertical vortex motions are able to attract small magnetic patches, amplify magnetic field by stretching the field lines, and, via downdrafts, transport the field into deeper layers, where the field can be further amplified by compression (Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To investigate the statistical relationship between twisting motions and magnetic field generation, we calculate cross-correlations of the squared magnetic field strength, B², with the kinetic helicity density and with the magnitude of the shear stress, $| S| ={(2{S}_{{ij}}{S}_{{ji}})}^{1/2}$ (Moin et al. 1991; Figure 5). The result is that the distributed magnetic field patches have better correlation with the kinetic helicity density (Figure 5(a)) than with the magnitude of shear stress (Figure 5(b)). These correlation functions decrease with depth more slowly for the kinetic helicity density due to transport of the magnetic field by helical downflows (Figure 4). Thus, we conclude that the small-scale dynamo under solar conditions primarily develops due to local fluctuations of the kinetic helicity density and that shearing fluctuations might also contribute to the small-scale dynamo.

Zoom In Zoom Out Reset image size

Because magnetic field generation and dissipation are highly inhomogeneous and depth dependent, we analyze areas where the magnetic energy grows and dissipates (Figure 6). We find that the time derivative of the magnetic energy density is small above the photosphere, while at the photosphere the growth rate of the magnetic energy is faster than dissipation. In the subsurface layers, the rates of growth and dissipation of magnetic energy are very similar. In layers deeper than 1 Mm, the rates of energy variation slowly decrease and are characterized by increasing amplitude of the local energy fluctuations. We plan to investigate the atmospheric effects in a separate paper.

Zoom In Zoom Out Reset image size

Magnetic field amplification due to swirling and shearing flows has been investigated by many authors (see, e.g., Batchelor 1950; Vainshtein & Ruzmaikin 1972; Nordlund et al. 1992; Petrovay & Szakaly 1993; Childress & Gilbert 1995; Brandenburg et al. 1996, 2012). Stretching, twisting, and folding of magnetic field lines induce growth of the magnetic flux under certain plasma conditions, such as a sufficiently high magnetic Reynolds number. The magnetic Prandtl number, which is the ratio of the diffusion coefficients of momentum and magnetic field, also plays a critical role in dynamo conditions (Nordlund et al. 1992; Vögler & Schüssler 2007). Previously, it was shown that magnetic field can be easily amplified for large magnetic Prandtl numbers (≥1; e.g., see Pao 1963; Meneguzzi & Pouquet 1989; Brandenburg et al. 1996; Schekochihin et al. 2004b). However, estimates of the magnetic Prandtl number for the Sun, based on molecular diffusivities, vary with depth from ∼10⁻¹ to 10⁻² near the base of the convective zone to ∼10⁻⁵ at the surface (Rieutord & Rincon 2010). This suggests that the dynamo cannot operate efficiently since, under these conditions, magnetic diffusion is much higher than momentum diffusion. Recently, however, it was found that dynamos can work for small magnetic Prandtl numbers (e.g., Ponty et al. 2004; Schekochihin et al. 2005; Iskakov et al. 2007) and that the transition between the regimes where magnetic field is either amplified or diffused by turbulent plasma motions is described by a critical magnetic Prandtl number, Pr^c_m, below which the dynamo processes do not occur. For theoretical discussions of this problem in the case of incompressible artificially forced turbulence, see papers by Rogachevskii & Kleeorin (1997), Schekochihin et al. (2004a, 2007), and Tobias et al. (2011).

In our numerical simulations, the molecular scales are unresolved, and LES turbulence models are used to estimate the plasma dynamics and diffusivities due to motion on subgrid scales. Therefore, in our simulations, the magnetic Prandtl number reflects the properties of plasma and magnetic fields only on turbulent scales and is defined as the ratio of the local turbulent plasma and magnetic diffusion coefficients (${{Pr}}_{m}^{t}={\nu }^{t}/{\eta }_{m}^{t}$). Figure 7 illustrates typical snapshots of the vertical velocity, magnetic field strength, electric current density, the time derivative of magnetic energy density, the turbulent magnetic Prandtl number, and the turbulent magnetic diffusivity in a fraction of the simulation domain at the photosphere (z = 0). The simulations reveal concentration of the magnetic field in the intergranular lanes and highly inhomogeneous time evolution of the magnetic energy on scales smaller than the width of the lanes. In particular, we find that the magnetic energy grows mostly at the edge of granules and at the periphery of turbulent vortex tubes, where shear motions are strongest, whereas the decrease of magnetic energy mostly occurs in areas where strong downflows transport magnetic field into the deeper layers. This inhomogeneity is also reflected in the distributions of turbulent magnetic Prandtl number ${{Pr}}_{m}^{t}$ and magnetic diffusivity ${\eta }_{m}^{t}$. Small values of ${{Pr}}_{m}^{t}$ are mostly associated with areas where shearing flows are present and the magnetic diffusivity has its largest values.

Zoom In Zoom Out Reset image size

A statistical distribution of the turbulent magnetic Prandtl number (${{Pr}}_{m}^{t}$) for a simulation of 1 hr of solar time is shown in Figures 2 and 8. According to our simulation results, ${{Pr}}_{m}^{t}$ can vary by several orders of magnitude, and most values are in the range of 0.1–1 (Figure 2(b)). The local growth and decay of the magnetic energy density correlate with turbulent Prandtl values ${{Pr}}_{m}^{t}\sim {10}^{-2}-1$ (Figure 8(a)). The distribution of the magnetic energy density rate of change for different ranges of the magnetic Prandtl number (Figure 8) shows a distinct asymmetry between the decay and amplification of magnetic field, in particular, for smaller ${{Pr}}_{m}^{t}$. This is clear evidence that the local dynamo processes are particularly efficient in areas characterized by a low magnetic Prandtl number.

Zoom In Zoom Out Reset image size

The quiet-Sun regions are characterized by relatively weak mean magnetic flux on the solar surface, but this flux does affect the turbulent properties of convection on small scales, where magnetic field patches can partially suppress plasma motions. For instance, comparison of the turbulent spectra of the photosphere, calculated from simulations with and without magnetic field, shows that the slope of the inertial range changes from −11/5 in the purely hydrodynamic case to −5/3 in the presence of magnetic fields (Kitiashvili et al. 2013c). Thus, the dynamo-generated magnetic fields locally act on surrounding turbulent flows (the so-called "back reaction"). In the current simulations, which account for turbulent dynamics on subgrid scales, we are able to capture the complicated interaction and energy exchange between the small-scale fields and flows.

The dynamo process produces magnetic patches with a wide range of scales of magnetic field distribution (Figure 2(c)) and involves the redistribution of energy among different scales. Therefore, we consider the development of dynamo as an evolution of the turbulent kinetic energy density spectra, which show a redistribution of energy through different scales during the transition from the "seed" field to a saturated magnetoconvective regime (Figure 9(a)). In particular, we find that the initial development of magnetic patches on the smallest scales increases small-scale electric currents and kinetic energy (red curve on panel (a)). Then the energy shift from small to larger scales changes the slope of the spectra so that they approach k^−5/3 (green curve). At the same time, the continuing amplification of magnetic fields at large wavenumbers (small scales) causes suppression of the kinetic energy on these scales. Thus, the evolution of the kinetic energy spectra during the development of the dynamo process can be interpreted as partly energy transport from smaller to larger scales and partly as a transformation of turbulent kinetic energy into magnetic energy. Note that because of the relatively small size of the computational domain and occasional appearance of intermittent "bursts" of the local magnetic field generation (similar to previously observed bursts in DNS simulations by Pratt et al. 2013), variations of the energy spectra can be due to transient energy redistributions across turbulent scales (e.g., see the spectrum shown by the blue curve in Figure 9(a)). Curiously, in the saturation phase (Figure 9(b)), the kinetic energy density spectrum at the photosphere is of the Kolmogorov type (k^−5/3), whereas the results by Rempel (2014) do not reproduce this slope. Also, the kinetic energy density in Rempel's simulations at k = 10 Mm⁻¹ is an order of magnitude higher and at k = 100 Mm⁻¹ is almost two orders of magnitude higher compared to our simulations. We do not know the reason for this difference, but perhaps the higher kinetic energy density is a result of the shallower simulation domain (2.4 Mm) in Rempel's simulations. The characteristic scale of the solar granulation at the photosphere is in agreement with both observations and Rempel's simulations (note the definition of the horizontal scale f = k_h/2π in his paper). In deeper layers, the slope of the spectrum in the inertial range decreases in magnitude: the power law is closer to k^−6/5 at a depth of 1 Mm, and the slope is even smaller in deeper layers.

Zoom In Zoom Out Reset image size

Figure 10 shows time and depth variations of the power spectra of the magnetic energy density from 2 Mm below the solar surface up to the photosphere. We found significant variations of the magnetic energy spectrum during the first half-hour after the "seed" field initialization, i.e., during the exponential growth phase. Later in time, the shape of the spectrum does not change, and the spectral power shows a steady energy increase at all scales. The slope of the magnetic energy spectra for small wavenumbers changes after field initialization. It reaches a value of 1/3 after ∼3–4 hr in the subsurface layers and a bit earlier in the photosphere. Our simulation results show the existence of a k⁻¹ power law in the saturation phase only for depths of about of 2 Mm and below, where the distribution of the magnetic field is more homogeneous. We do not reproduce the theoretical k^−11/3 law (Golitsyn 1960; Moffatt 1961). This possibly reflects deviations from the theoretical assumptions of magnetic field weakness, and thus the absence of magnetic field back reaction, in previous models. Note that Rempel (2014) did not see a clear indication of a power law for the magnetic energy in his simulations. Figure 11 shows the vertical distribution of the Poynting flux averaged over one hour and the horizontal magnetic flux density transported by upflows and downflows. The magnitude of the Poynting flux is similar to the values found in Rempel's simulations for the open boundary/zero field boundary conditions. However, in our simulations, the Poynting flux is decreasing at depths below 3 Mm (Figure 11(a)), while in Rempel's simulations it continues to increase. The reason for this difference is unclear and deserves further detailed investigation. Most of the horizontal flux is transported downflows into the deeper convection zone (Figure 11(b)). The mean flux recirculating by upflows is several times smaller than the flux in downflows.

Zoom In Zoom Out Reset image size

Understanding the local small-scale magnetic field generation process is critical for studying complex solar MHD problems, such as the interaction of magnetic fields and flows through different scales. In addition to the statistical properties described in previous sections, it is important to consider this process in terms of local physical properties. Here, we present two characteristic cases of local magnetic field amplification.

The first case, demonstrated in Figure 12, shows the evolution of the vertical velocity, the vertical magnetic field (color background), the horizontal velocity field (arrows), the electric current density (color map), and enstrophy (contour lines) in a selected 400 km × 400 km region of the photosphere with a 15 s cadence.

Zoom In Zoom Out Reset image size

The development of a local bipolar magnetic structure is associated with strong (∼5 km s⁻¹) swirling motions, on scales from ∼12–25 to ∼300–400 km, stretching and twisting the magnetic field, which is also compressed by converging flows in the downdrafts. Such small-scale swirling motions in the intergranular lanes are typically associated with strong downflows, ∼6–8 km s⁻¹ (Kitiashvili et al. 2010, 2011). The appearance of a bipolar magnetic structure (with a prominent negative (blue) polarity in Figure 12(b)) is also a result of swirling flows driven by a vortex tube oriented along the solar surface (indicated by the white arrow in Figure 12(a)), similar to one previously described by Steiner et al. (2010). The horizontally oriented vortex tube captures magnetic field lines and drags them into the subsurface layers. The dynamics of the positive-polarity patch is mostly related to the vertically oriented helical motions similar to the process described by Kitiashvili et al. (2010, 2011), who showed that the process of spontaneous formation of stable magnetic structures consists of two basic steps: (1) formation of small-scale filamentary magnetic structures associated with concentrations of vorticity and whirlpool-type motions, and (2) merging of these structures due to vortex attraction, caused by converging downdrafts around magnetic concentration below the surface. Such complicated flow dynamics is caused by interaction of differently oriented vortex tubes, which, in fact, have orientations varying in space and only locally can be regarded as having "vertical" or "horizontal" orientation (Kitiashvili et al. 2012b). The generated electric current (Figure 12(c)) shows a clear correlation with the swirling dynamics of convective flows. The strongest electric current density corresponds to areas on the periphery of helical motions or between the vortices, where stretching of the magnetic field lines is the strongest.

Another example illustrates the development of a bipolar magnetic structure due to helical flows and shearing flows along an intergranular lane. Figure 13 shows a time sequence of the vertical velocity, the vertical magnetic field, the magnitude of the electric current density, the time derivative of the vertical magnetic field with overlaid contour lines of enstrophy, and the vertical component of the electric current density (dashed curves correspond to negative values). An interesting feature of this example is a different scenario for the development of magnetic elements of opposite polarity (determined by the sign of the vertical field component): the positive-polarity magnetic patch evolves following the "classical" scenario of field amplification due to helical flows described in the first example, whereas the negative-polarity patch starts forming in the intergranular lane mostly due to converging flows (magnetic collapse). The evolution of the magnetic elements of this bipolar structure is accompanied by locally growing electric current (Figures 13(c) and (d)). Comparison of the velocity, pattern, and time derivative of the vertical magnetic field shows that some swirling motions induce magnetic field dissipation, probably due to scattering of the field lines, and often the field amplification takes place on the periphery of swirls or regions with shearing flows.

Zoom In Zoom Out Reset image size

The distribution of vertical electric current density (Figure 14) shows a strong correlation with the horizontal velocity of swirling flows in the subsurface layers. In deeper layers, the magnitude of the electric current increases, but its distribution becomes diffuse and does not show a clear association with the near-surface dynamics. The kinetic helicity density patterns show that the scale of the swirling motions increases with depth, from ∼50 to 60 km at the photospheric layer (Figure 15) up to 120 km at a depth of 300 km below the photosphere (panel (d)). In the deeper layers, the scale of helical motions continues to increase (to larger than 150 km), but the distribution of kinetic helicity density becomes complicated and consists of opposite-sign helical flows, which continue swirling together (Figure 15(e)), and disappear at a depth of about 500 km below the solar surface (panel (f)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From previous numerical studies, it is known that important dynamical and energetic links between subsurface turbulent convective flows and the low atmosphere are established through small-scale vortex tubes. In the presence of a magnetic field, the vortex tubes represent channels of energy exchange between the convective layers and the chromosphere and can result in the heating of chromospheric layers. This may be a source of small-scale spicule-like eruptions and Alfvén waves (Kitiashvili et al. 2012b, 2013b).

In addition to these effects, our simulations of the local dynamo show that the anisotropy between the vertical and horizontal magnetic field components takes place only in the photosphere and above. The nature of this anisotropy has been the subject of recent debate (e.g., Orozco Suárez et al. 2007; Lites et al. 2008; Danilovic et al. 2010; Ishikawa & Tsuneta 2010; Steiner & Rezaei 2012; Stenflo 2013). The anisotropy changes with height, but there is no dependence on the "seed" field. The vertical distribution of the mean unsigned magnetic field components (Figure 16(a)) shows a slow increase of their strength with depth that is a result of compression. In the convection zone, the distribution of the vertical and horizontal magnetic fields is statistically similar. This difference between the vertical and horizontal field strengths is small with a weak dominance of the horizontal fields below the solar surface. Indeed, the ratio between the transverse field, defined as ${B}_{h}=\sqrt{{B}_{x}^{2}+{B}_{y}^{2}}$, and the unsigned vertical field is only slightly above $\sqrt{2}$ (Figure 16(b)), meaning that the strength of all the field components is similar. Near the surface layers, the vertical magnetic flux becomes dominant. Above the photosphere, the mean vertical unsigned flux slowly decreases, whereas the horizontal field increases and becomes dominant (Figure 16). Note that our results are different from the previous simulations by Schüssler & Vögler (2008) and Rempel (2014), who found that the mean horizontal component of magnetic field is larger than the vertical component in the whole range, from the deep photosphere to the chromosphere. Our results indicate that the horizontal component is dominant only in the upper photosphere and low chromosphere, from about 300 to 700 km.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This dominance of the vertical or horizontal fields in the different layers reflects the topological properties of the dynamo-generated magnetic fields, which are characterized by field lines organized in small-scale magnetic loops above the photosphere (Figure 17). Such a topological structure resembles the magnetic canopy suggested from observations (Giovanelli 1980; Jones & Giovanelli 1982; Schrijver & Title 2002). According to our simulation results, in the intergranular lanes, the magnetic field is mostly vertical and becomes horizontal above the photosphere. The height of these loops is greater when the magnetic field is stronger at the loop footpoints. Such a topology of the magnetic field lines was suggested both by observations (e.g., Lites et al. 2008) and simulations (e.g., Schüssler & Vögler 2008; Steiner et al. 2008). The horizontal magnetic fields are stronger at the granule edges and in the regions where strong turbulent motions are present. This tendency of the distribution of transverse magnetic fields was previously found in observations (Lites et al. 2008). For a detailed comparison of the simulation results with the observational data, it is necessary to perform spectropolarimetric analyses for the simulated data, taking into account the instrumental characteristics. We plan this work for a future paper.

Zoom In Zoom Out Reset image size

In addition, it is interesting to note that the closest opposite-polarity patches may not even be connected by magnetic field lines above the solar surface as simple magnetic loops, but instead interact through electric currents above and below the photosphere, forming very complicated structures. Figure 18 illustrates the topological structure of the electric current density streamlines above and below the photosphere (shown as a horizontal semi-transparent plane). Each streamline is tracked from a point in the region of positive-polarity (orange streamlines) and negative-polarity (dark blue) patches. The topological structure of the electric currents above the photosphere is often characterized by spirals, arcs, and large swirls. Below the solar surface, such a topology can represent highly turbulent flows, as in the case of the positive patch (orange streamlines in Figure 18), or as a very regular spiral structure, as in the case of the negative patch (blue lines). For instance, in Figure 18(a), the current streamline originating in the positive-polarity patch (orange) is strongly twisted around the negative patch.

Zoom In Zoom Out Reset image size

In this paper, we addressed the problem of small-scale (local) dynamos that are probably responsible for a significant fraction of the background magnetic fields on the solar surface. To investigate this problem, we used a 3D radiative MHD code and performed several simulation runs for different strengths (10⁻⁶–10⁻² G) and spatial distributions of the initial seed field (Table 1). After injection of the seed field, no magnetic flux was added to or removed from the computational domain.We find that, in all simulation cases, magnetic field generation is qualitatively different only during the first few minutes after initialization. This transient reflects differences in the initial seed fields. However, after about 4 hr of solar time, the statistical distributions of the generated fields are similar. Therefore, in this paper, we mostly focused on an analysis of case E, in which the "seed" field was distributed as white noise with amplitude ±10⁻⁶ G. Our simulations show that, due to turbulent dynamo action, the magnetic field can be locally magnified above the equipartition strength (∼600 G), reaching more than 2000 G in the photosphere.

Our simulation results show that magnetic field amplification is driven by converging flows into the intergranular lanes, shearing flows, and helical motions. All these mechanisms of magnetic field amplification are present in our numerical model and are linked to each other. Thus, the local dynamo process represents a complicated interplay of multiple mechanisms that are difficult to separate from each other and that can contribute differently in individual magnetic field amplification events.

To model the local dynamo process, we use the LES approach, which allows us to resolve the essential scales of the turbulent solar plasma and model subgrid scales. Unlike in many other dynamo simulations, our numerical setup does not restrict the magnetic turbulent Prandtl number ${{Pr}}_{m}^{t}$ to a constant value. For instance, Figure 7 illustrates a local distribution of various properties on the solar surface and shows a correlation between the areas of magnetic energy growth and small values of the magnetic Prandtl number. The statistical distribution of the turbulent magnetic Prandtl number shows that local variations can be of several orders of magnitude, but most of the distribution is in the range ${{Pr}}_{m}^{t}\sim {10}^{-1}-1$ (Figure 2(b)).

The kinetic and magnetic energy density spectra during the development of the dynamo process show an intense energy exchange between different scales and between the kinetic and magnetic energies (Figures 9 and 10). The energy redistribution causes changes in the slope of the power law, which in the developed regime reaches the Kolmogorov law (k^−5/3).

The development of dynamo action can be seen in the animation association with Figures 1 and 13. Shortly after initialization of the "seed" field, magnetic patches appear in the photosphere. They are highly twisted and interact with each other. The magnetic energy in these patches increases with time. During the evolution of these magnetic elements, strong helical downflows transport the magnetic field from the subsurface into deeper layers. Some of these patches diffuse; however, new magnetic patches are continuously formed. Our simulations show that the magnetic patches can have magnetic fields stronger than 2 kG (Figure 2(c)).

Our analysis also shows that the local dynamo works most efficiently in the 1 Mm deep subsurface layer, where the turbulent flows are strongest (Figures 3 and 5). We confirm the role of vortical motions of plasma in magnetic field amplification due to twisting and stretching (Nordlund et al. 1992; Brandenburg et al. 1995, 1996). The magnetic patches generated due to shearing and vortical motions are transported into deeper layers by downdrafts. In the deeper layers, the magnetic field can be further amplified by compression.

We presented two characteristic examples in detailed studies: a bipolar magnetic structure generated by the interaction of vertical and horizontal vortex tubes, and magnetic field generation by a combined action of a vertical vortex tube and shearing flows. These examples illustrate three important properties of the process of magnetic field amplification. The first is a strong coupling of various magnetic field amplification mechanisms, which leads to increased local magnetic energy. The second key property of local dynamo action is the multi-scale nature of this process, which evolves in turbulent flows from the smallest resolved scales to granular scales (Figures 14 and 15). The third property is the complex topological and dynamical structure of the process, illustrated in Figure 18, and which is reflected in the interaction of individual magnetic patches and surrounding magnetic fields and in self-organization into a magnetic network. In this paper, we described just some of these properties of magnetic field evolution; further investigation is required. The primary topology of the dynamo-generated magnetic field is represented by compact magnetic loops appearing as bipolar structures in the intergranular lanes (Figure 17); these loops reach higher levels for stronger magnetic field concentrations at the footpoints of the loops.

Formation of magnetic loops in the solar atmosphere reflects the observed height-dependent anisotropy of vertical and transverse small-scale magnetic fields (Figure 16). We found an equipartition of the field component in the layers deeper than 2 Mm below the solar surface, but above the solar surface the transverse magnetic fields are dominant because of compact loop topology. This variation of the vertical and horizontal field anisotropy can explain discrepancies among different observations (e.g Orozco Suárez et al. 2007; Lites et al. 2008; Danilovic et al. 2010; Ishikawa & Tsuneta 2010; Stenflo 2013) and support previous numerical analysis (Schüssler & Vögler 2008; Steiner et al. 2008).

Investigation of the small-scale dynamo action is a key element for understanding dynamical and eruptive processes on the Sun and forms one of the building blocks of the global dynamics. In this paper, we suggest a new view of the problem through analysis of realistic-type radiative MHD simulations of the solar subsurface and chromosphere. We plan to continue analysis of these data, make detailed comparisons with observations, and clarify the physics of the small-scale dynamo action and organization of the Sun's magnetic network.

The simulation results were obtained on NASA's Pleiades supercomputer at NASA Ames Research Center. This work was partially supported by NASA grants NNX10AC55G and NNH11ZDA001N-LWSCSW, and Oak Ridge Associated Universities.