Free Access
Issue
A&A
Volume 660, April 2022
Article Number A115
Number of page(s) 12
Section The Sun and the Heliosphere
DOI https://doi.org/10.1051/0004-6361/202142644
Published online 22 April 2022

© 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 Rem = 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, Prm = Rem/Re = ν/η, with Prm in the 10−7 − 10−4 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, Reeff and Rem, 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−6 − 10−2 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 Prm ≳ 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 Prm ≪ 1 (Schekochihin et al. 2007). In this context, it is still unclear if limPrm → 0Rem, c = ∞, with Rem, c the critical magnetic Reynolds number for dynamo action, or if instead Rem, c saturates to some finite value in the limit Prm → 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, Prm, 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 Prm, 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 Prm, 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, Prm, eff. Second, we investigate how the magnetic field resulting from small-scale solar dynamo simulations depends on Reeff and Rem, 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 CO5BOLD code (Freytag et al. 2012), with the aim of characterising the magnetic energy growth rate and the saturated magnetic field in terms of Reeff and Rem, 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 CO5BOLD simulations, and we discuss the results. Small-scale dynamo simulations and their characterisation in terms of Reeff and Rem, 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 {Oi(v, B)} with weights {wi}. 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 Bh, such that ϵh = ∥B − Bh∥=C ⋅ hp + 𝒪(hp + 1), where ϵh is the numerical error affecting Bh, C is a constant independent of h, ∥ ⋅ ∥ is a designed norm, and 𝒪(hp + 1) denotes a term that decreases to zero as hp + 1 in the asymptotic regime. The same notation is used for the fluid velocity, where the numerical solution is denoted as vh. We also assume that we have a post-processing scheme of order of accuracy q ≥ p, where the discretised counterparts of ∂t(⋅) and {Oi(⋅, ⋅)} 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 Bh and vh, .

The next step consists of writing

(4)

where rh({δwi}) is the residual error and the δwi are interpreted as changes of the weights from wi to wi + δwi due to the numerical error introduced by the simulation code. The δwi can be estimated by minimising the residual error over a set of N ≥ 2 space-time points x1, …, xN. 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 wi, that is, the δwi 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 Bh and vh. Moreover, assuming that a non-negligible part of the numerical error can be modelled with a diffusion term of the form ηeff2B, we can extend the list of operators to include a discretised counterpart of ∇2B, (∇2)pp, hBh 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 aj, 3 = (∇2)pp, hBh(xj) 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 δwi = 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 wi 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 {Oi(ρ, 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 Prm, 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 CO5BOLD 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 ei. Here, etot = ρei + ρv2/2 + B2/2 + ρΦ is the total energy, Φ the gravitational potential, and Frad 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 (vA) is high, vA 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 CS 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 CO5BOLD, 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−2. The thermal energy equation was activated for β < 0.1 and vA was limited at a maximum of 40 km s−1. We used test calculations setting the limit to 90 km s−1 to verify that limiting vA has only a negligible impact on the results discussed in Sects. 4 and 51. 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 ∂tvz = −cszvz and ∂zvh = 0, with cs being the local sound speed and vh and vz 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 CsChange = 0.1 and CpChange = 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, tchar = Δz/⟨vf + |vz|⟩h a characteristic timescale with Δz being the grid spacing in the vertical direction and vf the flux velocity as obtained from the HLL scheme, and ϵ a tiny computer number, which is of order 10−300 cm s−1 in double precision format, introduced to avoid a vanishing denominator. The notation ⟨ − ⟩h denotes a horizontal average, and Clin and Csqrt are two input parameters (in the following Clin = Csqrt = 0 if not specified otherwise). Moreover, we have set B = 0h = 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 CO5BOLD 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 CO5BOLD 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 104 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 Bz = 100 G and Bz = 1 mG.

The first set of simulations (initial vertical magnetic field Bz = 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 × 1010 cm2 s−1; four simulations with h = 8 km and ηH = 4, 8, 12, and 16 × 1010 cm2 s−1; one simulation with h = 6 km and νH = 4.1 × 1010 cm2 s−1; one simulation with h = 6 km and CS = 1.55; and one simulation with h = 6 km and ηH = 3.5 × 1010 cm2 s−1. A summary of the numerical parameters for the first set of simulations is presented in Table 1.

Table 1.

Overview of simulations with an initial uniform magnetic field Bz = 100 G.

The second set of simulations, with a seed field Bz = 1 mG, is used in Sect. 5 to investigate dynamo action at Prm, 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, CS = 1.55, and ηH = 2.2,  2.5, and 3.2 × 1010 cm2 s−1; and two simulations with h = 6 km and ηH = 1.2 and 1.4 × 1010 cm2 s−1. 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 Clin = 0.0025 and Csqrt = 0.002, which were further increased to Clin = 0.008 and Csqrt = 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.

Table 2.

Overview of simulations with an initial uniform magnetic field Bz = 1 mG.

4. Evaluation of ν and η

In Sect. 2, it is explained that, after having obtained ρh, vh, and Bh (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 CS are input parameters of two versions of the stress tensor implemented in CO5BOLD 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 η2B and ∇ ⋅ τ, respectively. To investigate the validity of this assumption for the induction equation implemented in CO5BOLD, we computed the Pearson correlation coefficient rB between (∇2)pp, hBh 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 rB > 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 ρ2v 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 rB 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.

Table 3.

Pearson correlation coefficient rB.

Table 4.

Pearson correlation coefficient .

Table 5.

Pearson correlation coefficient .

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)efft, 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 coefficients2. 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%.

thumbnail 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 CO5BOLD, 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 δwi 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 CO5BOLD, 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 CO5BOLD. 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.

thumbnail 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, λ = (ν3/ε)1/4, where ε is the rate of kinetic energy dissipation per mass unit (Kolmogorov 1941). Assuming , with the velocity spectrum, kh 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 Rem, 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 Reeff(t, z) = Lvrms(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 Rem, eff(t, z) = Lvrms(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 Reeff and Rem, eff. To investigate this hypothesis, we carried out the simulations h6νH and h6νS, adjusting νH and CS 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 ⟨Reefft and ⟨Rem, efft 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 Reeff and Rem, 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.

thumbnail 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, Prm, 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 ⟨Prm, efft 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 ⟨Prm, efft slightly increasing in the convection zone with finer grids, whereas the opposite occurs in the photosphere.

thumbnail 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 CO5BOLD. 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 Reeff(z, t), Rem, eff(z, t), and Prm, eff(z, t) = Rem, eff(z, t)/Reeff(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 ⟨Em⟩ 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 ⟨Em⟩ 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 Rem, eff, with γE ∼ Rem, 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.

thumbnail 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 Tmd 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 Tms (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 Tms peaks on scales, λ = 2π/kh, 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 Tms is mostly compensated by Tma (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)

thumbnail 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 Tkt, 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 Tms) and λkt (the peak of Tkt), 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 Rem, 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 Clin = Csqrt = 0, obtaining τE = 1200 s and ⟨|Bz(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 Reeff and Rem, eff profiles of simulation h12 when using CS = 1.55 and ηH = 3.5 × 1010 cm2 s−1 for the simulation with h = 6 km. Therefore, we now investigate whether simulations of similar Reeff and Rem, 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 CS = 1.55 and ηH = 2.2,  2.5, and 3.2 × 1010 cm2 s−1 (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 × 1010 cm2 s−1 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 × 1010 cm2 s−1 or less, corresponding to Rem, 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 −∂tEm + Tms + Tma, 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 Tms, Tma, and Tkt are rather similar, the main differences being at very small scales. On the other hand, while for the simulations dh12 and dh6 the term Tmd overestimates −∂tEm + Tms + Tma 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, Rem, c, and that the small-scale dynamo operates even at Prm, eff ≃ 0.65. Comparing these results to the critical magnetic Reynolds numbers of Rem, c ≃ 200−300 obtained by Schekochihin et al. (2005) for Re ≃ 1000−2000, we note that CO5BOLD simulations have about a factor four larger Rem, 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 CO5BOLD 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 CO5BOLD 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 (Prm, eff ≃ 0.65), although these values are still far from the Prm ≪ 1 regime characterising the solar atmosphere.


1

In particular for the dynamo runs in Sect. 5, this limit becomes active in the high photosphere close to the top boundary only.

2

In the vicinity of shock fronts, a shock capturing numerical scheme automatically reduces the order of accuracy to first order, which increases the effective diffusivity of the scheme.

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

  1. Batchelor, G. K., & Taylor, G. I. 1950, Proc. Royal Soc. London Ser. A Math. Phys. Sci., 201, 405 [NASA ADS] [Google Scholar]
  2. Bodenheimer, P., Laughlin, G., Rozyczka, M., et al. 2006, Numerical Methods in Astrophysics: An Introduction, Series in Astronomy and Astrophysics (Boca Raton: CRC Press) [CrossRef] [Google Scholar]
  3. Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
  4. Buehler, D., Lagg, A., & Solanki, S. K. 2013, A&A, 555, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Cartier-Michaud, T., Ghendrih, P., Sarazin, Y., et al. 2016, Phys. Plasmas, 23, 020702 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cartier-Michaud, T., Galassi, D., Ghendrih, P., et al. 2020, Phys. Plasmas, 27, 052507 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cattaneo, F. 1999, ApJ, 515, L39 [Google Scholar]
  8. Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
  9. Freytag, B. 2013, Mem. Soc. Astron. It., 24, 26 [NASA ADS] [Google Scholar]
  10. Freytag, B. 2017, Mem. Soc. Astron. It., 88, 12 [NASA ADS] [Google Scholar]
  11. Freytag, B., Steffen, M., Ludwig, H.-G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
  12. Fromang, S., & Stone, J. M. 2009, A&A, 507, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [Google Scholar]
  14. Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2017, A&A, 604, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Kitiashvili, I. N., Kosovichev, A. G., Mansour, N. N., & Wray, A. A. 2015, ApJ, 809, 84 [NASA ADS] [CrossRef] [Google Scholar]
  16. Kleint, L., Berdyugina, S. V., Shapiro, A. I., & Bianda, M. 2010, A&A, 524, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 [Google Scholar]
  18. Laney, C. B. 1998, Computational Gasdynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  19. Lites, B. W., Centeno, R., & McIntosh, S. W. 2014, PASJ, 66, S4 [Google Scholar]
  20. Margolin, L. G. 2019, Shock Waves, 29, 27 [NASA ADS] [CrossRef] [Google Scholar]
  21. Martínez Pillet, V. 2013, Space Sci. Rev., 178, 141 [CrossRef] [Google Scholar]
  22. Nordlund, A. 1982, A&A, 107, 1 [Google Scholar]
  23. Orozco Suárez, D., & Bellot Rubio, L. R. 2012, ApJ, 751, 2 [Google Scholar]
  24. Petrovay, K., & Szakaly, G. 1993, A&A, 274, 543 [Google Scholar]
  25. Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606 [Google Scholar]
  26. Ramelli, R., Bianda, M., Berdyugina, S., Belluzzi, L., & Kleint, L. 2019, in Solar Polarization 8, eds. L. Belluzzi, R. Casini, M. Romoli, & J. Trujillo Bueno, ASP Conf. Ser., 526, 283 [Google Scholar]
  27. Rembiasz, T., Obergaulinger, M., Cerdá-Durán, P., Aloy, M.-Á., & Müller, E. 2017, ApJS, 230, 18 [Google Scholar]
  28. Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
  29. Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337 [Google Scholar]
  30. Sánchez Almeida, J., & Martínez González, M. 2011, in Solar Polarization 6, eds. J. R. Kuhn, D. M. Harrington, H. Lin, et al., ASP Conf. Ser., 437, 451 [Google Scholar]
  31. Sánchez Almeida, J., Márquez, I., Bonet, J. A., Domínguez Cerdeña, I., & Muller, R. 2004, ApJ, 609, L91 [CrossRef] [Google Scholar]
  32. Schekochihin, A. A., Haugen, N. E. L., Brandenburg, A., et al. 2005, ApJ, 625, L115 [NASA ADS] [CrossRef] [Google Scholar]
  33. Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, New J. Phys., 9, 300 [Google Scholar]
  34. Smagorinsky, J. 1963, Monthly Weather Rev., 91, 99 [NASA ADS] [CrossRef] [Google Scholar]
  35. Thaler, I., & Spruit, H. C. 2015, A&A, 578, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Trujillo Bueno, J., Shchukina, N., & Asensio Ramos, A. 2004, Nature, 430, 326 [Google Scholar]
  37. Vögler, A., & Schüssler, M. 2007, A&A, 465, L43 [Google Scholar]

All Tables

Table 1.

Overview of simulations with an initial uniform magnetic field Bz = 100 G.

Table 2.

Overview of simulations with an initial uniform magnetic field Bz = 1 mG.

Table 3.

Pearson correlation coefficient rB.

Table 4.

Pearson correlation coefficient .

Table 5.

Pearson correlation coefficient .

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.