© ESO 2022

1. Introduction

Over the past two decades, leveraging on the work done by Nordlund in the 1980s (Nordlund 1982), radiative magneto-hydrodynamic (MHD) simulations took on an increasingly important role in improving our understanding of the non-linear processes observed on the Sun. Unfortunately, the discretisation of an MHD model and its implementation in a simulation code always introduce numerical errors in the system, which in turn add numerical viscosity and magnetic diffusivity to the model equations (see, e.g. Laney 1998; Bodenheimer et al. 2006; Margolin 2019). This poses two critical problems. First, it is currently impossible to simulate the solar atmosphere at realistic Reynolds and magnetic Reynolds numbers, Re = LU/ν and Re_m = LU/η, with L and U the characteristic length scale and velocity of the plasma flow and ν and η the molecular viscosity and magnetic diffusivity (the latter also known as plasma resistivity). The reason is that the plasma flows in the solar atmosphere are characterised by large Reynolds and magnetic Reynolds numbers and by very small magnetic Prandtl numbers, Pr_m = Re_m/Re = ν/η, with Pr_m in the 10⁻⁷ − 10⁻⁴ range depending on the depth (see, e.g. Brandenburg & Subramanian 2005; Schekochihin et al. 2005). Therefore, it is necessary to resolve the wide range of length scales L ≫ l_η ≫ l_ν to simulate the solar atmosphere reliably (Schekochihin et al. 2007). Here, l_ν and l_η are the viscous and resistive length scales, respectively. The second issue is related to the fact that the intrinsic diffusivities of a simulation code complicate the interpretation of the results obtained with it, since the effective Reynolds and magnetic Reynolds numbers stemming from a simulation, Re_eff and Re_m, eff, are generally unknown.

One of the major open questions in solar physics concerns the origin of the ubiquitous small-scale magnetic field observed in the solar photosphere (see, e.g. Sánchez Almeida & Martínez González 2011). Since polarimetric signals measured in the inter-network of the quiet Sun are at the limit of present telescope capabilities and different observations of different resolutions and measurement techniques provide different results, it is unclear what the exact mean strength and spatial and angular distribution of the magnetic field in the solar atmosphere are (see, e.g. Martínez Pillet 2013; Khomenko et al. 2017, and references therein). Given these uncertainties, it is difficult to understand if the small-scale magnetic field of the quiet Sun mainly originates from a turbulent cascade acting on fields produced by a global-scale dynamo, or if it is mainly generated locally by a small-scale turbulent dynamo, as first suggested by Petrovay and Szakaly in 1993. The fact that properties of the quiet Sun’s magnetic field seem independent of the solar cycle (Trujillo Bueno et al. 2004; Sánchez Almeida et al. 2004; Buehler et al. 2013; Lites et al. 2014; Ramelli et al. 2019) would favour the latter explanation. However, hints of contradictory results were obtained by Kleint et al. (2010).

To shed light on this controversy, a number of numerical studies have been carried out in the past two decades. Initially, the possibility of driving dynamo action in the quiet photosphere through turbulent convection was demonstrated numerically by Cattaneo (1999). Later studies of small-scale dynamos in the solar atmosphere were performed with realistic solar-like numerical simulations by Vögler & Schüssler (2007), showing that a considerable amount of magnetic energy in the solar photosphere could be self-sustained through dynamo action. More recently, the possibility of amplifying a small seeded magnetic field (∼10⁻⁶ − 10⁻² G) in the solar photosphere through dynamo action was investigated by means of radiative MHD simulations in a number of additional publications (see, e.g. Pietarila Graham et al. 2010; Rempel 2014; Kitiashvili et al. 2015; Thaler & Spruit 2015; Khomenko et al. 2017).

Despite these numerical studies providing extremely useful information on many aspects of small-scale dynamo action in turbulent plasma flows, it is still unclear how to extrapolate these results to the Sun. This is for two reasons. First, the classical picture discussed in Batchelor & Taylor (1950) and generally used to explain field amplification at Pr_m ≳ 1, the regime typically accessible by current small-scale solar dynamo simulations, does not apply to the solar atmosphere. This is because the fundamental assumption that the scale of the fluid motion that does the stretching is larger than the scale of the field that is stretched, is not valid for Pr_m ≪ 1 (Schekochihin et al. 2007). In this context, it is still unclear if lim_{Prm → 0}Re_m, c = ∞, with Re_m, c the critical magnetic Reynolds number for dynamo action, or if instead Re_m, c saturates to some finite value in the limit Pr_m → 0. Second, the saturated averaged magnetic field resulting from small-scale dynamo simulations appears to strongly depend on the details of the simulation setup, and in particular on the effective magnetic Prandtl number, Pr_m, eff, which is generally unknown.

A crucial step in the study of small-scale dynamos in the solar atmosphere is therefore a reliable evaluation of Pr_m, eff and, consequently, of the effective dissipation coefficients stemming from a radiative MHD simulation, ν_eff and η_eff. In this respect, a common approach consists of adding explicit, artificial diffusion terms to the model equations and assuming that these terms are much larger than the intrinsic diffusivities of the simulation code. Several methodologies have also been proposed in the recent past for estimating the effective diffusivities stemming from an MHD simulation, or at least the resulting magnetic Prandtl number. These include: simulating simple one- or two-dimensional problems for which we know the exact solution of the resistive-viscous MHD equations and estimating the value of ν_eff and η_eff by fitting the results (Rembiasz et al. 2017); adding forcing terms to the model equations and associating the resulting saturated magnetic field to a turbulent magnetic diffusivity (Fromang & Stone 2009); and relating Pr_m, eff to the ratio of the velocity and magnetic energy Taylor microscales by means of a homogeneous turbulence simulation (Pietarila Graham et al. 2010). However, while these procedures provide reliable results under some particular conditions, they are difficult to generalise to different flow regimes and simulation setups. This motivates the present work, with a twofold objective. First, we propose a general and rigorous methodology for estimating ν_eff and η_eff, and, consequently, Pr_m, eff. Second, we investigate how the magnetic field resulting from small-scale solar dynamo simulations depends on Re_eff and Re_m, eff.

Our investigation is based on the method of projection on proper elements (PoPe), first introduced in Cartier-Michaud et al. (2016) and later extended to the independent PoPe (Cartier-Michaud et al. 2020), which has been used in the plasma physics community to quantify the error affecting numerical results and to verify plasma turbulence simulation codes. The estimate of the numerical error relies on using different, higher order accuracy numerical operators than the ones in the simulation code to post-process the simulation results. Then, we assume that the resulting numerical error can be, at least partially, modelled with a diffusion operator, and we compute the diffusivities through a least-squares fit. The proposed methodology is applied to a number of small-scale dynamo simulations carried out with the CO⁵BOLD code (Freytag et al. 2012), with the aim of characterising the magnetic energy growth rate and the saturated magnetic field in terms of Re_eff and Re_m, eff.

The remainder of this paper is structured as follows. In Sect. 2, we review the method of PoPe and we discuss how to extend it to compute the effective dissipation coefficients stemming from an MHD simulation and the resulting magnetic Prandtl number. In Sect. 3, we illustrate the numerical setup used for the simulations discussed in the present paper. In Sect. 4, we apply the methodology proposed in Sect. 2 to a number of CO⁵BOLD simulations, and we discuss the results. Small-scale dynamo simulations and their characterisation in terms of Re_eff and Re_m, eff are the subjects of Sect. 5. Finally, in Sect. 6, we report our conclusions.

2. Overview of the methodology

The methodology we propose for estimating the effective viscosity and magnetic diffusivity stemming from a radiative MHD simulation is based on the PoPe method recently detailed in Cartier-Michaud et al. (2016), which is briefly summarised here for completeness. To illustrate the main steps of the procedure, we consider the induction equation for an ideal plasma,

(1)

as an example, with v and B being the fluid velocity and the magnetic field, respectively. We note that v and B are functions of time and space, but this dependence is omitted in Eq. (1) to lighten the notation.

We rewrite Eq. (1) in the following more abstract form:

(2)

where

(3)

are the sequences of weights and operators, respectively, and i = 1, 2. This can be considered as the projection of the time derivative of B on the basis of elements {O_i(v, B)} with weights {w_i}. Moreover, we denote the magnetic field resulting from a simulation code implementing the discretised version of Eq. (1) of order of accuracy p and obtained on a mesh with degree of refinement h as B^h, such that ϵ^h = ∥B − B^h∥=C ⋅ h^p + 𝒪(h^p + 1), where ϵ^h is the numerical error affecting B^h, C is a constant independent of h, ∥ ⋅ ∥ is a designed norm, and 𝒪(h^p + 1) denotes a term that decreases to zero as h^p + 1 in the asymptotic regime. The same notation is used for the fluid velocity, where the numerical solution is denoted as v^h. We also assume that we have a post-processing scheme of order of accuracy q ≥ p, where the discretised counterparts of ∂_t(⋅) and {O_i(⋅, ⋅)} are denoted as and , respectively. We note that, in general, if the post-processing scheme does not correspond exactly to the numerical algorithm implemented in the simulation code used to obtain B^h and v^h, .

The next step consists of writing

(4)

where r^h({δw_i}) is the residual error and the δw_i are interpreted as changes of the weights from w_i to w_i + δw_i due to the numerical error introduced by the simulation code. The δw_i can be estimated by minimising the residual error over a set of N ≥ 2 space-time points x₁, …, x_N. In matrix form, this linear least-squares regression writes as

(5)

where A is a matrix with N × 2 elements , is a column vector with two elements , and b is a column vector with N elements . The PoPe methodology states that, if a code is bug free, the expected value of converges to w_i, that is, the δw_i vanish at rate p in the limit h → 0.

To proceed further, we note that if q > p and h is sufficiently small, the numerical error introduced by the post-processing scheme is negligible with respect to the numerical error introduced by the simulation code used to obtain B^h and v^h. Moreover, assuming that a non-negligible part of the numerical error can be modelled with a diffusion term of the form η_eff∇²B, we can extend the list of operators to include a discretised counterpart of ∇²B, (∇²)^pp, hB^h and project the residual error on it to compute the effective magnetic diffusivity η_eff. In the following, we consider two different methods for performing this projection. The first method, denoted in the following as method (i), consists of minimising

(6)

to obtain , and η_eff. This leads to adding a column to the matrix A with elements a_j, 3 = (∇²)^pp, hB^h(x_j) and a row to the vector with element and solving Eq. (5) for , and η_eff over a set of N ≥ 3 space-time grid points. The second method, (ii), assumes that δw_i = 0, and requires minimising the distance

(7)

over a set of N ≥ 1 space-time grid points to obtain η_eff. In other words, the first approach consists of splitting the numerical error in two sources, one that modifies the weights from w_i to and another that introduces numerical dissipation into the induction equation, whereas the second method assumes that the weights are unchanged and the numerical errors just increase the diffusion of the scheme. We note that the second approach can be understood as a more conservative estimate of the intrinsic magnetic diffusion and shares some similarities with the approach discussed in Khomenko et al. (2017) for the evaluation of η, and it was used in Cartier-Michaud et al. (2020) for estimating the effective diffusion affecting the continuity equation in the TOKAM3X code. We also note that, in general, we can solve Eq. (5) multiple times, each time considering a few different grid points, and therefore obtaining a statistical distribution of η_eff, or, we can solve Eq. (5) considering a large number of points at once, thus reducing the statistical error in the estimate of η_eff.

A similar procedure can be used to estimate the intrinsic kinematic viscosity affecting the momentum equation. More precisely, the momentum equation is written in a form similar to Eq. (2) as , where the lists of weights and operators are now

(8)

respectively, and i = 1, 2, 3, 4, with ρ being the mass density, P the plasma pressure, and g the gravity field. A post-processing scheme is then used to compute the discretised counterparts of these operators and build a matrix with N × 4 elements, with the rows corresponding to different space-time points and the columns to different operators. The same steps done to obtain η_eff are then repeated to obtain the effective kinematic viscosity, ν_eff; the only difference is the form of the diffusion operator being added to the list {O_i(ρ, v, B)}, which is denoted in the following as ∇ ⋅ τ, where τ is the stress tensor, which depends on ν. Finally, the effective magnetic Prandtl number stemming from an MHD simulation is computed as Pr_m, eff = ν_eff/η_eff. More details on the explicit form of τ, as well as a validation of the present methodology, are discussed in Sect. 4.

3. Numerical setup

The simulations discussed in the present paper were obtained with the CO⁵BOLD code, which is designed to simulate solar and stellar atmospheres, as well as their interior, by solving the time-dependent ideal MHD equations:

(9)

(10)

(11)

(12)

in combination with a non-grey radiative transfer equation and an equation of state (eos) that connects the mass density ρ and the gas pressure P to the internal energy e_i. Here, e_tot = ρe_i + ρv²/2 + B²/2 + ρΦ is the total energy, Φ the gravitational potential, and F_rad the frequency-integrated radiative energy flux vector. I is the identity matrix. The code can be used to carry out both global simulations of an entire star from the central core to the visible stellar surface, or at least a shell encompassing the convection zone, and local simulations of a small partial volume of the star under investigation, typically encompassing the surface layers where energy transport changes from predominantly convective to purely radiative.

The (magneto-)hydrodynamic equations are evolved using modern flux-conserving Riemann type solvers, with a Roe solver (Roe 1986) for non-magnetic simulations, or an HLL solver (Harten et al. 1983) that can be used both for the magnetic and the non-magnetic cases. In regions of small plasma-β, with β the ratio of the thermal to the magnetic pressure of the plasma, the equation of the thermal energy is used instead of the equation of the total energy to avoid negative gas pressure, this at the expense of strict energy conservation. Additionally, to avoid very small time steps when the Alfvén speed (v_A) is high, v_A can be limited by artificially reducing the strength of the Lorentz force. To ensure a divergence-free magnetic field, the constrained-transport method of Evans & Hawley (1988) is used. Furthermore, explicit artificial viscous and resistive terms can be added to the momentum and induction equations to increase the dissipation of the numerical scheme. In particular, in the following we consider two different forms of artificial viscosity and one of artificial magnetic diffusivity: a homogeneous viscous coefficient, ν_H; a turbulent sub-grid-scale viscosity as given by the Smagorinsky model (Smagorinsky 1963), where we denote as C_S the Smagorinsky coefficient; and a homogeneous magnetic diffusivity η_H.

Concerning the radiative transfer equation, this can be solved using either a short- or a long-characteristic scheme. Additionally, to reduce the computational cost of a simulation, the diffusion approximation can be used to model the radiative transport in the deep, optically thick layers of the numerical domain. Pre-tabulated values as functions of density and internal energy, accounting for the ionisation balance of hydrogen, helium, and a representative metal, are used to solve the eos and compute the plasma pressure and temperature.

The code is written in a modular form using Fortran 90 (and some Fortran 77) language, and it is parallelised using OpenMP directives. To carry out the simulations discussed in the present paper, the parallelisations of the MHD and of the long-characteristic transfer modules were extended to a hybrid, Message Passing Interface (MPI) and OpenMP, parallelisation. For additional details on CO⁵BOLD, we refer the reader to Freytag et al. (2012).

For the present paper, we considered only local simulations (box-in-a-star setup) with a Cartesian computational domain encompassing a volume of 6.0 × 6.0 × 2.4 Mm, with approximately 1.6 Mm below and 0.8 Mm above the mean Rosseland optical depth τ_R = 1. The transfer equation was solved with the long-characteristic radiation transport module, using Rosseland mean opacities (grey opacity), and employing the diffusion approximation for approximately the first 1 Mm from the bottom boundary of the domain. To compute the fluxes for the HLL solver, we used a FRweno reconstruction scheme (Freytag 2013). Moreover, the gravitational field, g, was vertical and uniform with a constant value of g = 27 500 cm s⁻². The thermal energy equation was activated for β < 0.1 and v_A was limited at a maximum of 40 km s⁻¹. We used test calculations setting the limit to 90 km s⁻¹ to verify that limiting v_A has only a negligible impact on the results discussed in Sects. 4 and 5¹. The equations were then evolved in time with an Hancock predictor step.

The side boundaries were periodic in both horizontal directions. The top boundary was open for fluid flows and outward radiation, with the density decreasing exponentially in the boundary (ghost) cells outside the domain, whereas magnetic fields were forced to be vertical. Moreover, to avoid reflection of acoustic waves at the top boundary, we enforced ∂_tv_z = −c_s∂_zv_z and ∂_zv_h = 0, with c_s being the local sound speed and v_h and v_z the horizontal and vertical velocities, respectively. At the bottom we have set an open boundary condition for fluid variables, enforcing ∂_zv = 0, vanishing horizontally averaged vertical mass flux, and a control on the specific entropy of the inflowing material that ensures that the model has a resulting effective temperature close to 5770 K. The two parameters C_sChange = 0.1 and C_pChange = 0.3 were used to reduce deviations of the entropy and the pressure from the horizontal means, as detailed in Freytag et al. (2012). Similarly, a damping of the vertical velocity was introduced at the bottom boundary in the following form:

(13)

where Δt is the internal time step, t_char = Δz/⟨v_f + |v_z|⟩_h a characteristic timescale with Δz being the grid spacing in the vertical direction and v_f the flux velocity as obtained from the HLL scheme, and ϵ a tiny computer number, which is of order 10⁻³⁰⁰ cm s⁻¹ in double precision format, introduced to avoid a vanishing denominator. The notation ⟨ − ⟩_h denotes a horizontal average, and C_lin and C_sqrt are two input parameters (in the following C_lin = C_sqrt = 0 if not specified otherwise). Moreover, we have set B = 0_h = 0 in inflow regions at the bottom boundary, ensuring that no magnetic energy enters the domain at these locations. With regard to the small-scale dynamo simulations of Sect. 5, this is a very conservative condition, because it allows for the removal of magnetic energy without replenishing it from outside of the computational domain. More details on the boundary conditions used in CO⁵BOLD are given in Freytag (2017).

In this work, we considered the five grid resolutions h = 24, 16, 12, 8, and 6 km, where h denotes a uniform grid spacing in the three spatial coordinates. We started from a previous non-magnetic CO⁵BOLD simulation of the solar atmosphere, adapted it to a computational domain encompassing a volume of 6.0 × 6.0 × 2.4 Mm, interpolated the results to a grid resolution of h = 24 km, and evolved the resulting plasma fields for more than 10 000 s to damp out fluctuations introduced by this interpolation. Then, for each finer resolution, we started from a relaxed, coarser simulation, we interpolated the numerical results to the new grid, and we carried out a non-magnetic run for at least 10⁴ s. From the final state of the five relaxed non-magnetic simulations, we started two sets of magnetic simulations, with an initial uniform vertical magnetic field B_z = 100 G and B_z = 1 mG.

The first set of simulations (initial vertical magnetic field B_z = 100 G) is used in Sect. 4 to investigate the methodology presented in Sect. 2. To this aim, we carried out a total of 16 simulations, each for 1000 s: five reference simulations with h = 24, 16, 12, 8, and 6 km; four simulations with h = 8 km and ν_H = 4, 8, 12, and 16 × 10¹⁰ cm² s⁻¹; four simulations with h = 8 km and η_H = 4, 8, 12, and 16 × 10¹⁰ cm² s⁻¹; one simulation with h = 6 km and ν_H = 4.1 × 10¹⁰ cm² s⁻¹; one simulation with h = 6 km and C_S = 1.55; and one simulation with h = 6 km and η_H = 3.5 × 10¹⁰ cm² s⁻¹. A summary of the numerical parameters for the first set of simulations is presented in Table 1.

The second set of simulations, with a seed field B_z = 1 mG, is used in Sect. 5 to investigate dynamo action at Pr_m, eff ≲ 1. To this aim, we carried out a total of ten simulations: five reference simulations with h = 24, 16, 12, 8, and 6 km; three simulations with h = 6 km, C_S = 1.55, and η_H = 2.2,  2.5, and 3.2 × 10¹⁰ cm² s⁻¹; and two simulations with h = 6 km and η_H = 1.2 and 1.4 × 10¹⁰ cm² s⁻¹. Because of some numerical instabilities emerging at low grid spacings, for all simulations with h < 16 km we introduced artificial numerical dissipation at the bottom boundary by setting C_lin = 0.0025 and C_sqrt = 0.002, which were further increased to C_lin = 0.008 and C_sqrt = 0.004 for the last 2 h of simulation dh6. An investigation of the impact of these two parameters on dynamo action is discussed in Sect. 5. Table 2 summarises the numerical setups of the dynamo runs.

4. Evaluation of ν and η

In Sect. 2, it is explained that, after having obtained ρ^h, v^h, and B^h (the numerical results from the simulations described in Sect. 3), we need a post-processing scheme to discretise the operators , and ∇^pp, h and apply the methodology we propose. To ensure that the choice of the post-processing scheme does not influence our conclusions, in the following we consider three finite difference schemes of different order of accuracy: second, fourth, and sixth order, denoted as FD2, FD4, and FD6, respectively. Moreover, we need to choose how to express the term ∇ ⋅ τ. To investigate if the particular form of ∇ ⋅ τ has an impact on the resulting effective viscosity, in the following we consider the two expressions:

(14)

(15)

The first expression corresponds to the typical form used to express the divergence of the stress tensor in the context of viscous flows, whereas the second expression represents a simple diffusion term in the velocity equation with diffusion coefficient ν^b. We note that ν^a and ν^b are free parameters of the post-processing methodology that are used to estimate the effective viscosity affecting the numerical results, whereas ν_H and C_S are input parameters of two versions of the stress tensor implemented in CO⁵BOLD and used to add explicit viscosity to the momentum equation.

In Sect. 2 it is assumed that part of the numerical errors introduced by the discretisation of the induction and momentum equations can be expressed as η∇²B and ∇ ⋅ τ, respectively. To investigate the validity of this assumption for the induction equation implemented in CO⁵BOLD, we computed the Pearson correlation coefficient r_B between (∇²)^pp, hB^h and (with weights and operators in Eq. (3)) for the simulations h24, h16, h12, h8, and h6 and at t = 500 s and t = 1000 s.

The results are displayed in Table 3. We find that r_B > 0.76 independently of the order of accuracy of the post-processing scheme FD2-FD6 used for the computation. We repeated the same analysis for the momentum equation, for which we computed the Pearson correlation coefficients and between ∇ ⋅ {ρ[∇v + (∇v)^T − 2∇⋅v/3]} and , and between ρ∇²v and , respectively (with weights and operators in Eq. (8)), again for the simulations h24, h16, h12, h8, and h6 and at t = 500 s and t = 1000 s. The results are displayed in Tables 4 and 5. It is found that . Although the Pearson correlation coefficients for the momentum equation are smaller than for the induction equation, , , and r_B are noticeably larger than 0, implying that a considerable linear correlation exists between the dissipation terms and the residuals of the corresponding equations. This is consistent with our assumption that the residuals can be, to a large degree, expressed in the form of diffusion terms.

We can now apply the methodology of Sect. 2 and estimate ν_eff and η_eff. To this aim, in the following, we minimise the residual with the procedures (i) and (ii), Eqs. (6) and (7), considering all the spatial grid points on a horizontal plane at once and plane by plane independently, thus obtaining effective dissipation coefficients that depend on time and height, ν_eff(t, z) and η_eff(t, z). In Fig. 1 we display ⟨ν_eff(t, z)⟩_t and ⟨η(t, z)_eff⟩_t, where the time averages are performed over the time instants t = 500 s and t = 1000 s, computed with both procedures (i) and (ii) and considering both forms of ∇ ⋅ τ in Eqs. (14) and (15), for the simulations h24, h12, and h6. As we could expect, the effective viscosity and magnetic diffusivity decrease by refining the grid. Moreover, it appears that the profiles vary less in the convection zone (z < 0 km) than in the photosphere (z > 0 km), the main exception being ν_eff near the bottom boundary. We think that the steep increase in the upper photosphere, above z = 300 km, is due to the presence of shock waves in this region, which lead to larger gradients of the plasma velocity and of the magnetic field, and therefore to larger dissipation coefficients². We also see that the results do not have a strong dependence on the post-processing scheme employed to compute the coefficients, except for some small local deviations and for the evaluation of η_eff with the procedure (ii). Also, ν_eff is almost independent of the form chosen to express ∇ ⋅ τ, that is, the left and middle panels of Fig. 1 are almost identical. Concerning the differences resulting from employing Eq. (6) or Eq. (7), procedures (i) or (ii), we observe that the second approach generally produces slightly larger coefficients, as we could expect according to our remark in Sect. 2. Nevertheless, the differences are rather small, in particular in the convection zone where they are typically smaller than 15%. We note that very similar trends are obtained for the simulations h16 and h8 (not displayed in Fig. 1 for clarity). We also tested a different approach for minimising Eq. (6) or Eq. (7): we considered all the spatial grid points in the numerical domain at once, instead of independently plane by plane, thus obtaining global dissipation coefficients that only depend on time (not shown here for conciseness). They confirm the trends discussed above, with marginal dependence on the post-processing scheme and resulting coefficients that are slightly larger for procedure (ii) than for (i), with differences that are generally smaller than 10%.

Fig. 1. Effective viscosities and plasma resistivity as a function of height. Time-averaged profiles of (left panel), (centre panel), and ηeff(t, z) (right panel) for the simulations h24 (solid lines), h12 (dotted lines), and h6 (dashed lines). Colour code: blue, red, and yellow for method (i) with post-processing scheme FD2, FD4, and FD6, respectively, and purple, green, and light blue for method (ii) with post-processing scheme FD2, FD4, and FD6, respectively.

Given these results, in the following we consider only the post-processing scheme FD2, the methodology (ii), and we express ∇ ⋅ τ according to Eq. (15) (we omit the superscript b to lighten the notation). We note that, in order to solve Eq. (7) for ν_eff(t, z) and η_eff(t, z), we implemented the necessary routines directly in CO⁵BOLD, such that we can compute the dissipation coefficients at run time with much higher output frequency, without the need for storing the full snapshots at these high output frequencies. In the remainder of this section, we consider an output frequency of the numerical results of 5 s, with time averages over the interval t ∈ [500, 1000] s. The standard deviations of the time signals of the dissipation coefficients ν_eff(t, z) and η_eff(t, z) during this interval are rather small. They are generally below 5%, except at certain locations, in particular near the top and bottom boundaries, where they can increase to 10–20%.

The PoPe methodology states that if a code is bug free, then δw_i should vanish at rate p in the asymptotic limit h → 0. This should also hold true for ν_eff and η_eff. To verify if it is the case for CO⁵BOLD, Fig. 2 (left) displays the effective viscosity and magnetic diffusivity for the simulations h24, h16, h12, h8, and h6, averaged over the vertical coordinate z, as a function of the grid spacing h. The dissipation coefficients are clearly decreasing as we refine the grid, at a rate between 1 (linear) and 2 (quadratic). This is consistent with the order of accuracy of the FRweno scheme implemented in CO⁵BOLD. We remind the reader that this result is independent of the order of accuracy of the post-processing scheme used to compute ν_eff and η_eff; thus, it indicates that we are not measuring an extra dissipation introduced by our methodology. We repeated the same analysis for the simulations h8ν_H4, h8ν_H8, h8ν_H12, and h8ν_H16 (Fig. 2, centre), and for the simulations h8η_H4, h8η_H8, h8η_H12, and, h8η_H16 (Fig. 2, right), displaying the results as function of ν_H and η_H, respectively. We obtain η_eff ≃ η_H for large enough artificial magnetic diffusivities, whereas ν_eff and η_eff converge to the intrinsic numerical diffusivity of the scheme for ν_H, η_H → 0. The effective numerical viscosity ν_eff is slightly larger than ν_H for finite artificial viscosities. In Fig. 2, we also observe that η_eff decreases with increasing ν_H, and similarly ν_eff decreases with increasing η_H. We speculate that this is related to the frozen-in condition exhibited by plasmas under almost ideal conditions, which entails the smooth behaviour of the magnetic field with increasing viscosity and the smooth behaviour of the flow with increasing magnetic diffusivity.

Fig. 2. Time and spatial averages of νeff (red squares) and ηeff (blue circles) for the simulations h24, h16, h12, h8, and h6 as function of grid spacing h (left panel); for the simulations h8, h8νH4, h8νH8, h8νH12, and h8νH16 as function of νH (centre panel); and for the simulations h8, h8ηH4, h8ηH8, h8ηH12, and h8ηH16 as function of ηH (right panel). The dashed and dotted black lines in the left panel are proportional to h2 and h, respectively, whereas the dashed black lines in the centre and right panels are proportional to νH and ηH, respectively. The error bars denote one standard deviation computed over time and vertical direction z.

To gain a deeper insight into the methodology of Sect. 2, it is useful to remember that, according to the classical picture of turbulent energy cascade from Kolmogorov’s theory, the kinetic energy that resides mostly at large scales is transferred across scales by non-linear processes and is eventually dissipated in the form of heat at about the Kolmogorov dissipation scale, λ = (ν³/ε)^1/4, where ε is the rate of kinetic energy dissipation per mass unit (Kolmogorov 1941). Assuming , with the velocity spectrum, k_h the horizontal wave number, and and the Fourier transform of v and its complex conjugate, respectively, we can compute the Kolmogorov dissipation scale for the simulations h24, h16, h12, h8, and h6. We obtain ⟨λ⟩_t, z = 23.6, 16.5, 12.6, 8.8, and 6.9 km, that is, ⟨λ⟩_t, z ≈ h. This is consistent with an intrinsic viscosity acting approximately at the grid scale.

In the context of Kolmogorov’s theory it is also often assumed that, for large enough Re and Re_m, the statistics of velocity fluctuations at small scales is uniquely determined by ν and ε, and that η affects the distribution of magnetic fields only at the smallest scales (Thaler & Spruit 2015). Accordingly, two simulations with similar Re_eff(t, z) = Lv_rms(t, z)/ν_eff(t, z), where L = 1 Mm is the typical length scale and is the root-mean-square value of the fluid velocity, should exhibit similar velocity spectra; analogously, two simulations with similar Re_m, eff(t, z) = Lv_rms(t, z)/η_eff(t, z) should have similar magnetic energy spectra. Consequently, a high-resolution simulation with explicit viscosity and magnetic diffusivity is expected to behave like a lower resolution simulation of similar Re_eff and Re_m, eff. To investigate this hypothesis, we carried out the simulations h6ν_H and h6ν_S, adjusting ν_H and C_S to approximately replicate the effective viscosity of the simulation h12 (similarly, we adjusted η_H in h6η_H to replicate the effective magnetic diffusivity of simulation h12). The resulting vertical profiles of ⟨Re_eff⟩_t and ⟨Re_m, eff⟩_t are displayed in Fig. 3, left and centre, respectively. While the profiles exhibit some local differences, which can be expected due to the different nature of the dissipation terms introduced in the different simulations, the overall agreement is rather good, in particular if compared to the results of simulation h6, displayed in red in the aforementioned panels of Fig. 3. Figure 3 (right) displays the velocity and magnetic energy spectra, normalised to their integral area, and , averaged over time and between z = −200 km and z = 0 km. Concerning the velocity spectra, there is indeed good agreement among the simulations h12, h6ν_H and h6ν_S, although some local differences are present. These spectra are clearly offset from the velocity spectrum of simulation h6. Concerning the magnetic spectra, there is a better agreement between h6η_H and h12 than between h6 and h12, although noticeable differences appear at small scales. These results reveal that it is not enough to replicate the values of Re_eff and Re_m, eff, or equivalently ν_eff and η_eff, to obtain the same spectra; the nature of the dissipation terms plays a role in setting the profiles of and . Evaluating the average velocity and magnetic spectra in different depths in the convection zone (not shown here) yielded the same results.

Fig. 3. Reynolds and magnetic Reynolds numbers as a function of height and kinetic and magnetic energy spectra. Left: time-averaged profiles of Reeff(t, z) as function of geometrical height z for the simulations h12, h6, h6νH, and h6νS. Centre: time-averaged profiles of Rem, eff(t, z) for the simulations h12, h6, and h6ηH. Right: Spectra of (solid lines) and (dashed lines) averaged over time and in space between z = −200 km and z = 0 km as a function of horizontal wave number kh. The spectra of are multiplied by 104 for clarity. The shaded grey areas of the left and centre panels denote the region between z = −200 km and z = 0 km, used for the spatial average of the energy spectra.

In light of the results discussed in the present section, we are now confident that our methodology provides a reliable estimate of ν_eff and η_eff. We can therefore compute the magnetic Prandtl number, Pr_m, eff = ν_eff/η_eff, as discussed in Sect. 2. The time-averaged results for the simulations h24, h16, h12, h8, and h6 are displayed in Fig. 4. The profiles of ⟨Pr_m, eff⟩_t are quite flat in the convection zone, with a value slightly above one, the only exception being near the bottom boundary, where the magnetic Prandtl number drops below one. There is a trend related to the grid resolution, with ⟨Pr_m, eff⟩_t slightly increasing in the convection zone with finer grids, whereas the opposite occurs in the photosphere.

Fig. 4. Time-averaged profiles of the effective magnetic Prandtl number Prm, eff as a function of geometrical height z for the simulations h24, h16, h12, h8, and h6.

5. Small-scale dynamo simulations

Dynamo runs with increasing spatial resolution, and consequently increasing Reynolds and magnetic Reynolds numbers, were carried out in order to investigate the possibility of simulating small-scale dynamo action with CO⁵BOLD. The corresponding runs and numerical parameters are summarised in Table 2. There, the Reynolds and magnetic Reynolds numbers, as well as the magnetic Prandtl numbers, are estimated as the time averages over the kinematic phase and the spatial averages for z < 0 km of Re_eff(z, t), Re_m, eff(z, t), and Pr_m, eff(z, t) = Re_m, eff(z, t)/Re_eff(z, t), respectively.

The resulting time evolution of the magnetic energy density, averaged over the full domain, for the five simulations dh24, dh16, dh12, dh8, and dh6, is displayed in Fig. 5, left. The simulations start with a short transient phase of a few minutes where magnetic flux is expelled from the granules and transported to the intergranular space. After that, the temporal profiles of ⟨E_m⟩ are rather flat for h = 16 km and h = 24 km (dashed lines), and the overall magnetic energy density remains extremely low. On the other hand, for grids with h < 16 km, we clearly recognise a linear (kinematic) phase, during which magnetic energy is amplified exponentially, followed by a non-linear (saturation) phase occurring once the impact of the Lorentz force on the fluid motion starts to be significant. The magnetic energy e-folding times reported in Table 2, τ_E, are then obtained by fitting the ⟨E_m⟩ temporal profiles during the kinematic phase with an exponential function. As one would expect for a small-scale dynamo, the growth rate, γ_E = 1/τ_E, is an increasing function of Re_m, eff, with γ_E ∼ Re_m, eff, which is consistent with what was found in Pietarila Graham et al. (2010) and Rempel (2014). A detailed investigation of the time evolution of the magnetic energy density in planes of different heights in the convection zone (not shown here for conciseness) yields energy e-folding times very similar to those presented in Table 2, which refer to the full domain. On the other hand, the magnetic energy is amplified much slower in the upper photosphere. For this reason, in the following we focus our attention to the convection zone alone, where we expect a more effective dynamo.

Fig. 5. Growth of magnetic energy with time and kinetic and magnetic spectra. Left: temporal profiles of the magnetic energy density averaged over the full domain as function of time for simulations dh24, dh16, dh12, dh8, dh6, dh6νsηH25, and dh6ηH12. Solid lines denote simulations showing dynamo action, whereas dashed lines denote simulations without dynamo action. Right: spectra of the normalised kinetic energy density (solid lines) and magnetic energy density (dashed lines) averaged over the time span of the kinematic phase and in space between z = −900 km and z = −700 km.

Figure 5 (right) shows the kinetic and magnetic energy spectra and , normalised to their integral area and averaged over time during the kinematic phase and in space between z = −900 km and z = −700 km, for the simulations dh12, dh8, and dh6. Clearly, the magnetic energy spectra peak at scales much smaller than the kinetic energy spectra, which peak at a scale of about 1.5 Mm, and their peaks shift to smaller scales when increasing the resolution. This is consistent with small-scale dynamo action. However, without further analysis it cannot be excluded that the observed amplification of small-scale magnetic fields originates from mechanisms other than the action of a small-scale dynamo, such as a turbulent tangling (or shredding) of the large-scale magnetic fields that could be produced by a mean-flow (1 Mm scale) dynamo or by the Alfvénic response of large-scale fields to small-scale velocity fluctuations (turbulent induction). To rule out these possibilities and disentangle the different sources of small-scale magnetic energy, we follow the approach proposed in Pietarila Graham et al. (2010) and Rempel (2014) and carry out a spectral analysis of the energy exchanges between kinetic and magnetic energy reservoirs. More precisely, as explained in the appendices of Pietarila Graham et al. (2010) and Rempel (2014), starting from the induction equation, we can decompose the time evolution of the spectral magnetic energy, , in three contributions converting kinetic energy into magnetic energy and vice versa as

(16)

where

(17)

(18)

(19)

are the energy transfers to and from the magnetic energy due to the misalignment between velocity shear and magnetic field lines (stretching), to advection and compression, and to magnetic diffusivity, respectively, with denoting the Fourier transform, ^* the complex conjugate, and c.c. the respective complex conjugate expression. The approximation symbol in Eq. (16) highlights the fact that this is not exactly equal, because T_md is an estimate of the magnetic energy losses due to numerical dissipation obtained under the assumption that the residual of the induction equation can be modelled with a Laplacian magnetic diffusivity. In the following, all the terms related to the energy transfer rates are computed using second-order, centred, finite difference schemes.

Figure 6, left and centre, displays the spectral energy transfer rates, normalised to the magnetic energy density and averaged over time during the kinematic phase and in space between z = −900 km and z = −700 km, for the two simulations dh12 and dh6, respectively. From T_ms (blue solid curve), we see that sources of magnetic energy due to the stretching of magnetic field lines are present at all scales of the simulations and that T_ms peaks on scales, λ = 2π/k_h, of about λ_ms ≃ 200 km and λ_ms ≃ 140 km for simulations dh12 and dh6, respectively; that is, the dominant scale for magnetic energy production by stretching is approximately λ_ms ≃ 17 h to λ_ms ≃ 23 h (black vertical dashed lines in Fig. 6). On the other hand, energy transfer by advection and compression results in a removal of magnetic energy from large scales (yellow dashed curves), breaking it down into smaller scale magnetic fields (solid red curves). At scales larger than 28 h, the term T_ms is mostly compensated by T_ma (solid blue versus dashed yellow curves). The contributions from numerical dissipation, shown as dashed purple curves in Fig. 6, dominate at small scales, approximately at λ_md ≃ 10−14 h (black vertical dash-dotted lines in Fig. 6). Additionally, from the momentum equation, we compute the energy transfer term representing kinetic energy losses due to the work of the Lorenz force via magnetic tension,

(20)

Fig. 6. Normalised energy transfer functions averaged over the time span of the kinematic phase and in space between z = −900 km and z = −700 km. Left: for the simulation dh12; centre: for the simulation dh6; and right: for the simulation dh6νSηH25. Solid colour curves represent sources of magnetic energy due to stretching (blue curves) and advection and compression (red curves) of magnetic fields, dashed colour curves represent magnetic energy losses due to advection and compression (yellow curves) and a Laplacian diffusion (purple curves) of magnetic fields. The green dashed curves represent the kinetic energy losses due to stretching of magnetic fields, whereas the light blue dotted lines represent the energy transfer imbalance −∂tEm + Tms + Tma. Black vertical lines denote the position of λms (dashed), λkt (dotted), and λmd (dash-dotted).

This is displayed in Fig. 6 in green. From T_kt, we see that the stretching of magnetic field lines is responsible for losses of kinetic energy, transferred to the magnetic energy reservoir, at scales of approximately λ_kt ≃ 24 h to λ_kt ≃ 36 h (black vertical dotted lines in Fig. 6).

From this analysis, it is clear that a range of scales exists between λ_ms (the peak of T_ms) and λ_kt (the peak of T_kt), where magnetic energy injection is mostly due to turbulent stretching of magnetic field lines, with the injected energy coming from the kinetic energy reservoir. Moreover, the magnetic energy spectra in Fig. 5, right, peak on scales of approximately λ_Em ≃ 27 h. Since these three scales, λ_ms, λ_Em, and λ_kt, all lie in the inertial range (see Fig. 5, right) and are incompatible with small-scale velocity fluctuations interacting with a large-scale magnetic field to produce small-scale magnetic energy, we conclude that the magnetic energy is indeed amplified by small-scale dynamo action (Pietarila Graham et al. 2010). We note that similar results are obtained for the simulation dh8, although they are not displayed here for conciseness.

In Table 2, we also report the time- and horizontally averaged unsigned vertical magnetic field and magnetic-to-kinetic energy density ratio, obtained during the saturation phase at z = 0. For the case of h = 12 km, for which the magnetic Reynolds number probably does not exceed Re_m, c enough, the magnetic energy and the unsigned vertical magnetic field settle to rather low values. On the other hand, the temporally and horizontally averaged unsigned vertical magnetic field reaches 49 G for h = 6 km, with a magnetic to kinetic energy ratio of about 1.4%. We remind the reader that, while these values are about a factor of 1.5–3 smaller than the magnetic field strengths observed on the quiet Sun (see, e.g. Trujillo Bueno et al. 2004; Orozco Suárez & Bellot Rubio 2012), the overall efficiency of our dynamo is strongly affected by the details of the numerical setup used in these simulations. As a matter of fact, we repeated the simulation dh8 with C_lin = C_sqrt = 0, obtaining τ_E = 1200 s and ⟨|B_z(z = 0)|⟩ = 30 G. More importantly, the bottom boundary condition used in these simulations is very conservative, since we assumed a zero horizontal magnetic field in inflow regions. As also discussed in Rempel (2014), this means that all magnetic field that is continuously pumped out of the computational domain across the bottom boundary is lost for the dynamo, and it is not replenished by magnetic field carrying upflows.

In Sect. 4, we see that it is possible to approximately replicate the time-averaged vertical Re_eff and Re_m, eff profiles of simulation h12 when using C_S = 1.55 and η_H = 3.5 × 10¹⁰ cm² s⁻¹ for the simulation with h = 6 km. Therefore, we now investigate whether simulations of similar Re_eff and Re_m, eff produce similar e-folding times. To this aim, we carried out three additional simulations, always starting from a seed magnetic field of 1 mG, with C_S = 1.55 and η_H = 2.2,  2.5, and 3.2 × 10¹⁰ cm² s⁻¹ (simulations dh6ν_Sη_H22, dh6ν_Sη_H25, and dh6ν_Sη_H32, respectively). From a study of the time evolution of the total magnetic energy in the boxes, it results that with η_H = 3.2 × 10¹⁰ cm² s⁻¹ or larger, it is not possible to observe dynamo action (in Fig. 5, left, we show only the time-evolution profile of simulation dh6ν_Sη_H25 for clarity). On the other hand, the magnetic energy grows exponentially when η_H = 2.5 × 10¹⁰ cm² s⁻¹ or less, corresponding to Re_m, eff = 860 or larger, although we could not afford to run until saturation, because of the extremely large computational cost of these simulations. Comparing the results in Table 2 to simulation dh12, we find that, when using explicit viscosities and magnetic diffusivities, we need larger effective magnetic Reynolds numbers with respect to the case of using only implicit numerical dissipations to obtain dynamo action. To further analyse this kind of hysteresis, Fig. 6 also displays the imbalance −∂_tE_m + T_ms + T_ma, with , indicated with light blue dotted curves, as an estimate of the numerical dissipation, as well as the energy transfer rates for simulation dh6ν_Sη_H25 (right panel). Comparing the results of the simulations dh12 and dh6ν_Sη_H25, we see that the profiles of T_ms, T_ma, and T_kt are rather similar, the main differences being at very small scales. On the other hand, while for the simulations dh12 and dh6 the term T_md overestimates −∂_tE_m + T_ms + T_ma at large scales, but much better agrees at smaller scales (purple versus light blue curves), we see a much better agreement between the purple dashed and light blue dotted curves for the simulation dh6ν_Sη_H25. This implies that, to obtain an equal growth rate, it is not enough to have the same effective Reynolds and magnetic Reynolds numbers, it is also important that the dissipation terms operate at the same length scales.

Finally, the simulations dh6η_H12 and dh6η_H14 were carried out to investigate the possibility of amplifying magnetic energy through dynamo action at smaller magnetic Prandtl numbers than 0.9. The corresponding magnetic energy e-folding times are reported in Table 2, while the time evolution of the magnetic energy resulting from the simulation dh6η_H12 is displayed in Fig. 5. The results suggest that the magnetic Reynolds numbers of these two simulations are very close to, but larger than, Re_m, c, and that the small-scale dynamo operates even at Pr_m, eff ≃ 0.65. Comparing these results to the critical magnetic Reynolds numbers of Re_m, c ≃ 200−300 obtained by Schekochihin et al. (2005) for Re ≃ 1000−2000, we note that CO⁵BOLD simulations have about a factor four larger Re_m, c. This means that we need a correspondingly smaller magnetic diffusivity to obtain dynamo action in these simulations.

6. Conclusions

In the present paper, we propose a methodology, based on the method of PoPe of Cartier-Michaud et al. (2016), for estimating the effective magnetic Prandtl number stemming from radiative MHD simulations. The methodology is simple, general, and can be applied as a post-processing step after having obtained the simulation results, without the need for modifying the simulation code or knowing the exact details of the employed numerical scheme. It just requires introducing a different numerical scheme than the one used in the simulation code (in particular it should be of higher order of accuracy), recomputing the operators of the momentum and induction equations with the new scheme, and then minimising a residual through classical residual methods.

The application of the proposed methodology to a number of CO⁵BOLD simulations allowed us to investigate the advantages and drawbacks of the procedure. In particular, it is found that the results are almost insensitive to the order of accuracy of the post-processing scheme we tested. Also, the results are insensitive to the exact form used to express the stress tensor. Moreover, it is shown that it is also possible to obtain dissipation coefficients that depend on height, that the resulting effective dissipations are in rather good agreement with explicit diffusion coefficients if these are large enough, and that our results are consistent with an intrinsic viscosity acting approximately at the grid scale. Overall, we are now confident that the proposed methodology provides a solid estimate of the dissipation coefficients affecting the momentum and induction equations of MHD simulation codes and, consequently, a reliable evaluation of the magnetic Prandtl number characterising the numerical results. It also emerged that having the same effective dissipation coefficients does not ensure we obtain the same magnetic and kinetic energy spectra, as implicit and explicit dissipations might operate at different scales.

Finally, the possibility of simulating small-scale dynamo action with CO⁵BOLD was investigated. It is found that, for resolutions higher than or equal to h = 12 km, small-scale dynamos are active and can amplify a small seed magnetic field up to significant values. However, as also found by other authors in the past, the small-scale dynamo strongly depends on the details of the numerical setup, and in particular of the boundary conditions used at the bottom of the numerical boxes and the exact nature of the magnetic diffusivity. It was also possible to identify dynamo action at rather low magnetic Prandtl numbers (Pr_m, eff ≃ 0.65), although these values are still far from the Pr_m ≪ 1 regime characterising the solar atmosphere.

Acknowledgments

This work was supported by the Swiss National Science Foundation under grant ID 200020_182094. Part of the numerical simulations were carried out on Piz Daint at CSCS under project IDs sm51, s1059, and u14, and the rest was carried out on the HPC ICS cluster at USI. We sincerely thank the anonymous referee for the careful reading of the manuscript and the many comments, which greatly helped us to improve the quality of the paper.

References

All Tables

All Figures

Fig. 1. Effective viscosities and plasma resistivity as a function of height. Time-averaged profiles of (left panel), (centre panel), and ηeff(t, z) (right panel) for the simulations h24 (solid lines), h12 (dotted lines), and h6 (dashed lines). Colour code: blue, red, and yellow for method (i) with post-processing scheme FD2, FD4, and FD6, respectively, and purple, green, and light blue for method (ii) with post-processing scheme FD2, FD4, and FD6, respectively.

In the text

Fig. 2. Time and spatial averages of νeff (red squares) and ηeff (blue circles) for the simulations h24, h16, h12, h8, and h6 as function of grid spacing h (left panel); for the simulations h8, h8νH4, h8νH8, h8νH12, and h8νH16 as function of νH (centre panel); and for the simulations h8, h8ηH4, h8ηH8, h8ηH12, and h8ηH16 as function of ηH (right panel). The dashed and dotted black lines in the left panel are proportional to h2 and h, respectively, whereas the dashed black lines in the centre and right panels are proportional to νH and ηH, respectively. The error bars denote one standard deviation computed over time and vertical direction z.

In the text

Fig. 3. Reynolds and magnetic Reynolds numbers as a function of height and kinetic and magnetic energy spectra. Left: time-averaged profiles of Reeff(t, z) as function of geometrical height z for the simulations h12, h6, h6νH, and h6νS. Centre: time-averaged profiles of Rem, eff(t, z) for the simulations h12, h6, and h6ηH. Right: Spectra of (solid lines) and (dashed lines) averaged over time and in space between z = −200 km and z = 0 km as a function of horizontal wave number kh. The spectra of are multiplied by 104 for clarity. The shaded grey areas of the left and centre panels denote the region between z = −200 km and z = 0 km, used for the spatial average of the energy spectra.

In the text

	Fig. 4. Time-averaged profiles of the effective magnetic Prandtl number Prm, eff as a function of geometrical height z for the simulations h24, h16, h12, h8, and h6.
In the text

Fig. 5. Growth of magnetic energy with time and kinetic and magnetic spectra. Left: temporal profiles of the magnetic energy density averaged over the full domain as function of time for simulations dh24, dh16, dh12, dh8, dh6, dh6νsηH25, and dh6ηH12. Solid lines denote simulations showing dynamo action, whereas dashed lines denote simulations without dynamo action. Right: spectra of the normalised kinetic energy density (solid lines) and magnetic energy density (dashed lines) averaged over the time span of the kinematic phase and in space between z = −900 km and z = −700 km.

In the text

Fig. 6. Normalised energy transfer functions averaged over the time span of the kinematic phase and in space between z = −900 km and z = −700 km. Left: for the simulation dh12; centre: for the simulation dh6; and right: for the simulation dh6νSηH25. Solid colour curves represent sources of magnetic energy due to stretching (blue curves) and advection and compression (red curves) of magnetic fields, dashed colour curves represent magnetic energy losses due to advection and compression (yellow curves) and a Laplacian diffusion (purple curves) of magnetic fields. The green dashed curves represent the kinetic energy losses due to stretching of magnetic fields, whereas the light blue dotted lines represent the energy transfer imbalance −∂tEm + Tms + Tma. Black vertical lines denote the position of λms (dashed), λkt (dotted), and λmd (dash-dotted).

In the text