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/00046361/202142644  
Published online  22 April 2022 
Methodology for estimating the magnetic Prandtl number and application to solar surface smallscale dynamo simulations
^{1}
Istituto Ricerche Solari Locarno (IRSOL), Università della Svizzera italiana (USI), 6605 LocarnoMonti, Switzerland
email: fabio.riva@irsol.usi.ch
^{2}
LeibnizInstitut für Sonnenphysik (KIS), Schöneckstrasse 6, 79104 Freiburg i.Br., Germany
Received:
11
November
2021
Accepted:
6
February
2022
Context. A crucial step in the numerical investigation of smallscale dynamos in the solar atmosphere consists of an accurate determination of the magnetic Prandtl number, Pr_{m}, stemming from radiative magnetohydrodynamic (MHD) simulations.
Aims. The aims are to provide a reliable methodology for estimating the effective Reynolds and magnetic Reynolds numbers, Re and Re_{m}, and their ratio Pr_{m} = Re_{m}/Re (the magnetic Prandlt number), that characterise MHD simulations and to categorise smallscale dynamo simulations in terms of these dimensionless parameters.
Methods. The methodology proposed for computing Re and Re_{m} is based on the method of projection on proper elements and it relies on a postprocessing step carried out using higher order accurate numerical operators than the ones in the simulation code. A number of radiative MHD simulations with different effective viscosities and plasma resistivities were carried out with the CO^{5}BOLD code, and the resulting growth rate of the magnetic energy and saturated magnetic field strengths were characterised in terms of Re and Re_{m}.
Results. Overall, the proposed methodology provides a solid estimate of the dissipation coefficients affecting the momentum and induction equations of MHD simulation codes, and consequently also a reliable evaluation of the magnetic Prandtl number characterising the numerical results. Additionally, it is found that smallscale dynamos are active and can amplify a small seed magnetic field up to significant values in CO^{5}BOLD simulations with a grid spacing smaller than h = 12 km, even at Pr_{m} ≃ 0.65. However, it is also evident that it is difficult to categorise dynamo simulations in terms of Pr_{m} alone, because it is not only important to estimate the amplitude of the dissipation coefficients, but also at which scales energy dissipation takes place.
Key words: Sun: magnetic fields / magnetohydrodynamics (MHD) / dynamo / Sun: photosphere / methods: numerical / turbulence
© ESO 2022
1. Introduction
Over the past two decades, leveraging on the work done by Nordlund in the 1980s (Nordlund 1982), radiative magnetohydrodynamic (MHD) simulations took on an increasingly important role in improving our understanding of the nonlinear 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^{−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, Re_{eff} and Re_{m, eff}, are generally unknown.
One of the major open questions in solar physics concerns the origin of the ubiquitous smallscale magnetic field observed in the solar photosphere (see, e.g. Sánchez Almeida & Martínez González 2011). Since polarimetric signals measured in the internetwork 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 smallscale magnetic field of the quiet Sun mainly originates from a turbulent cascade acting on fields produced by a globalscale dynamo, or if it is mainly generated locally by a smallscale 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 smallscale dynamos in the solar atmosphere were performed with realistic solarlike numerical simulations by Vögler & Schüssler (2007), showing that a considerable amount of magnetic energy in the solar photosphere could be selfsustained 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 smallscale 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 smallscale 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 smallscale 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 smallscale 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 twodimensional problems for which we know the exact solution of the resistiveviscous 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 smallscale 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 CartierMichaud et al. (2016) and later extended to the independent PoPe (CartierMichaud 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 postprocess 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 leastsquares fit. The proposed methodology is applied to a number of smallscale dynamo simulations carried out with the CO^{5}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^{5}BOLD simulations, and we discuss the results. Smallscale 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 CartierMichaud 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,
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:
where
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 postprocessing 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 postprocessing 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
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 spacetime points x_{1}, …, x_{N}. In matrix form, this linear leastsquares regression writes as
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 postprocessing 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 nonnegligible part of the numerical error can be modelled with a diffusion term of the form η_{eff}∇^{2}B, we can extend the list of operators to include a discretised counterpart of ∇^{2}B, (∇^{2})^{pp, h}B^{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
to obtain , and η_{eff}. This leads to adding a column to the matrix A with elements a_{j, 3} = (∇^{2})^{pp, h}B^{h}(x_{j}) and a row to the vector with element and solving Eq. (5) for , and η_{eff} over a set of N ≥ 3 spacetime grid points. The second method, (ii), assumes that δw_{i} = 0, and requires minimising the distance
over a set of N ≥ 1 spacetime 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 CartierMichaud 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
respectively, and i = 1, 2, 3, 4, with ρ being the mass density, P the plasma pressure, and g the gravity field. A postprocessing 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 spacetime 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^{5}BOLD code, which is designed to simulate solar and stellar atmospheres, as well as their interior, by solving the timedependent ideal MHD equations:
in combination with a nongrey 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}/2 + B^{2}/2 + ρΦ is the total energy, Φ the gravitational potential, and F_{rad} the frequencyintegrated 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 fluxconserving Riemann type solvers, with a Roe solver (Roe 1986) for nonmagnetic simulations, or an HLL solver (Harten et al. 1983) that can be used both for the magnetic and the nonmagnetic 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 divergencefree magnetic field, the constrainedtransport 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 subgridscale 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 longcharacteristic 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. Pretabulated 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 longcharacteristic transfer modules were extended to a hybrid, Message Passing Interface (MPI) and OpenMP, parallelisation. For additional details on CO^{5}BOLD, we refer the reader to Freytag et al. (2012).
For the present paper, we considered only local simulations (boxinastar 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 longcharacteristic 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 v_{A} 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 v_{A} has only a negligible impact on the results discussed in Sects. 4 and 5^{1}. 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 ∂_{t}v_{z} = −c_{s}∂_{z}v_{z} and ∂_{z}v_{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 ∂_{z}v = 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:
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^{−300} cm s^{−1} 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 smallscale 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^{5}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 nonmagnetic CO^{5}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 nonmagnetic run for at least 10^{4} s. From the final state of the five relaxed nonmagnetic 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^{10} cm^{2} s^{−1}; four simulations with h = 8 km and η_{H} = 4, 8, 12, and 16 × 10^{10} cm^{2} s^{−1}; one simulation with h = 6 km and ν_{H} = 4.1 × 10^{10} cm^{2} s^{−1}; one simulation with h = 6 km and C_{S} = 1.55; and one simulation with h = 6 km and η_{H} = 3.5 × 10^{10} cm^{2} s^{−1}. A summary of the numerical parameters for the first set of simulations is presented in Table 1.
Overview of simulations with an initial uniform magnetic field B_{z} = 100 G.
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^{10} cm^{2} s^{−1}; and two simulations with h = 6 km and η_{H} = 1.2 and 1.4 × 10^{10} cm^{2} 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 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.
Overview of simulations with an initial uniform magnetic field B_{z} = 1 mG.
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 postprocessing scheme to discretise the operators , and ∇^{pp, h} and apply the methodology we propose. To ensure that the choice of the postprocessing 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:
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 postprocessing 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^{5}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 η∇^{2}B and ∇ ⋅ τ, respectively. To investigate the validity of this assumption for the induction equation implemented in CO^{5}BOLD, we computed the Pearson correlation coefficient r_{B} between (∇^{2})^{pp, h}B^{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 postprocessing scheme FD2FD6 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 ρ∇^{2}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.
Pearson correlation coefficient r_{B}.
Pearson correlation coefficient .
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)_{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^{2}. We also see that the results do not have a strong dependence on the postprocessing 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 postprocessing 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. Timeaveraged 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 postprocessing scheme FD2, FD4, and FD6, respectively, and purple, green, and light blue for method (ii) with postprocessing scheme FD2, FD4, and FD6, respectively. 
Given these results, in the following we consider only the postprocessing 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^{5}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^{5}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^{5}BOLD. We remind the reader that this result is independent of the order of accuracy of the postprocessing 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ν_{H}4, h8ν_{H}8, h8ν_{H}12, and h8ν_{H}16 (Fig. 2, centre), and for the simulations h8η_{H}4, h8η_{H}8, h8η_{H}12, and, h8η_{H}16 (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 frozenin 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ν_{H}4, h8ν_{H}8, h8ν_{H}12, and h8ν_{H}16 as function of ν_{H} (centre panel); and for the simulations h8, h8η_{H}4, h8η_{H}8, h8η_{H}12, and h8η_{H}16 as function of η_{H} (right panel). The dashed and dotted black lines in the left panel are proportional to h^{2} 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 nonlinear 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, 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 rootmeansquare 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 highresolution 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: timeaveraged profiles of Re_{eff}(t, z) as function of geometrical height z for the simulations h12, h6, h6ν_{H}, and h6ν_{S}. Centre: timeaveraged profiles of Re_{m, 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 k_{h}. The spectra of are multiplied by 10^{4} 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 timeaveraged 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. Timeaveraged profiles of the effective magnetic Prandtl number Pr_{m, eff} as a function of geometrical height z for the simulations h24, h16, h12, h8, and h6. 
5. Smallscale 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 smallscale dynamo action with CO^{5}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 nonlinear (saturation) phase occurring once the impact of the Lorentz force on the fluid motion starts to be significant. The magnetic energy efolding 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 smallscale 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 efolding 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}η_{H}25, and dh6η_{H}12. 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 smallscale dynamo action. However, without further analysis it cannot be excluded that the observed amplification of smallscale magnetic fields originates from mechanisms other than the action of a smallscale dynamo, such as a turbulent tangling (or shredding) of the largescale magnetic fields that could be produced by a meanflow (1 Mm scale) dynamo or by the Alfvénic response of largescale fields to smallscale velocity fluctuations (turbulent induction). To rule out these possibilities and disentangle the different sources of smallscale 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
where
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 secondorder, 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 dashdotted 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,
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}η_{H}25. 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 −∂_{t}E_{m} + T_{ms} + T_{ma}. Black vertical lines denote the position of λ_{ms} (dashed), λ_{kt} (dotted), and λ_{md} (dashdotted). 
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 smallscale velocity fluctuations interacting with a largescale magnetic field to produce smallscale magnetic energy, we conclude that the magnetic energy is indeed amplified by smallscale 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 magnetictokinetic 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 timeaveraged vertical Re_{eff} and Re_{m, eff} profiles of simulation h12 when using C_{S} = 1.55 and η_{H} = 3.5 × 10^{10} cm^{2} s^{−1} for the simulation with h = 6 km. Therefore, we now investigate whether simulations of similar Re_{eff} and Re_{m, eff} produce similar efolding 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^{10} cm^{2} s^{−1} (simulations dh6ν_{S}η_{H}22, dh6ν_{S}η_{H}25, and dh6ν_{S}η_{H}32, respectively). From a study of the time evolution of the total magnetic energy in the boxes, it results that with η_{H} = 3.2 × 10^{10} cm^{2} s^{−1} or larger, it is not possible to observe dynamo action (in Fig. 5, left, we show only the timeevolution profile of simulation dh6ν_{S}η_{H}25 for clarity). On the other hand, the magnetic energy grows exponentially when η_{H} = 2.5 × 10^{10} cm^{2} s^{−1} 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 −∂_{t}E_{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}η_{H}25 (right panel). Comparing the results of the simulations dh12 and dh6ν_{S}η_{H}25, 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 −∂_{t}E_{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}η_{H}25. 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η_{H}12 and dh6η_{H}14 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 efolding times are reported in Table 2, while the time evolution of the magnetic energy resulting from the simulation dh6η_{H}12 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 smallscale 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^{5}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 CartierMichaud 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 postprocessing 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^{5}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 postprocessing 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 smallscale dynamo action with CO^{5}BOLD was investigated. It is found that, for resolutions higher than or equal to h = 12 km, smallscale 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 smallscale 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.
In particular for the dynamo runs in Sect. 5, this limit becomes active in the high photosphere close to the top boundary only.
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
 Batchelor, G. K., & Taylor, G. I. 1950, Proc. Royal Soc. London Ser. A Math. Phys. Sci., 201, 405 [NASA ADS] [Google Scholar]
 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]
 Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Buehler, D., Lagg, A., & Solanki, S. K. 2013, A&A, 555, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 CartierMichaud, T., Ghendrih, P., Sarazin, Y., et al. 2016, Phys. Plasmas, 23, 020702 [NASA ADS] [CrossRef] [Google Scholar]
 CartierMichaud, T., Galassi, D., Ghendrih, P., et al. 2020, Phys. Plasmas, 27, 052507 [NASA ADS] [CrossRef] [Google Scholar]
 Cattaneo, F. 1999, ApJ, 515, L39 [Google Scholar]
 Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 [Google Scholar]
 Freytag, B. 2013, Mem. Soc. Astron. It., 24, 26 [NASA ADS] [Google Scholar]
 Freytag, B. 2017, Mem. Soc. Astron. It., 88, 12 [NASA ADS] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H.G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Fromang, S., & Stone, J. M. 2009, A&A, 507, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harten, A., Lax, P. D., & van Leer, B. 1983, SIAM Rev., 25, 35 [Google Scholar]
 Khomenko, E., Vitas, N., Collados, M., & de Vicente, A. 2017, A&A, 604, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kitiashvili, I. N., Kosovichev, A. G., Mansour, N. N., & Wray, A. A. 2015, ApJ, 809, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Kleint, L., Berdyugina, S. V., Shapiro, A. I., & Bianda, M. 2010, A&A, 524, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 [Google Scholar]
 Laney, C. B. 1998, Computational Gasdynamics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Lites, B. W., Centeno, R., & McIntosh, S. W. 2014, PASJ, 66, S4 [Google Scholar]
 Margolin, L. G. 2019, Shock Waves, 29, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Martínez Pillet, V. 2013, Space Sci. Rev., 178, 141 [CrossRef] [Google Scholar]
 Nordlund, A. 1982, A&A, 107, 1 [Google Scholar]
 Orozco Suárez, D., & Bellot Rubio, L. R. 2012, ApJ, 751, 2 [Google Scholar]
 Petrovay, K., & Szakaly, G. 1993, A&A, 274, 543 [Google Scholar]
 Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606 [Google Scholar]
 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]
 Rembiasz, T., Obergaulinger, M., CerdáDurán, P., Aloy, M.Á., & Müller, E. 2017, ApJS, 230, 18 [Google Scholar]
 Rempel, M. 2014, ApJ, 789, 132 [Google Scholar]
 Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337 [Google Scholar]
 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]
 Sánchez Almeida, J., Márquez, I., Bonet, J. A., Domínguez Cerdeña, I., & Muller, R. 2004, ApJ, 609, L91 [CrossRef] [Google Scholar]
 Schekochihin, A. A., Haugen, N. E. L., Brandenburg, A., et al. 2005, ApJ, 625, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Schekochihin, A. A., Iskakov, A. B., Cowley, S. C., et al. 2007, New J. Phys., 9, 300 [Google Scholar]
 Smagorinsky, J. 1963, Monthly Weather Rev., 91, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Thaler, I., & Spruit, H. C. 2015, A&A, 578, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trujillo Bueno, J., Shchukina, N., & Asensio Ramos, A. 2004, Nature, 430, 326 [Google Scholar]
 Vögler, A., & Schüssler, M. 2007, A&A, 465, L43 [Google Scholar]
All Tables
All Figures
Fig. 1. Effective viscosities and plasma resistivity as a function of height. Timeaveraged 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 postprocessing scheme FD2, FD4, and FD6, respectively, and purple, green, and light blue for method (ii) with postprocessing 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ν_{H}4, h8ν_{H}8, h8ν_{H}12, and h8ν_{H}16 as function of ν_{H} (centre panel); and for the simulations h8, h8η_{H}4, h8η_{H}8, h8η_{H}12, and h8η_{H}16 as function of η_{H} (right panel). The dashed and dotted black lines in the left panel are proportional to h^{2} 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: timeaveraged profiles of Re_{eff}(t, z) as function of geometrical height z for the simulations h12, h6, h6ν_{H}, and h6ν_{S}. Centre: timeaveraged profiles of Re_{m, 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 k_{h}. The spectra of are multiplied by 10^{4} 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. Timeaveraged profiles of the effective magnetic Prandtl number Pr_{m, 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}η_{H}25, and dh6η_{H}12. 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}η_{H}25. 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 −∂_{t}E_{m} + T_{ms} + T_{ma}. Black vertical lines denote the position of λ_{ms} (dashed), λ_{kt} (dotted), and λ_{md} (dashdotted). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.