Draft version July 22, 2020

Preprint typeset using L^{A}TEX style emulateapj v. 08/22/09

NUMERICAL AND PHYSICAL MAGNETIC HELICITY PRODUCTION IN A KINEMATIC DYNAMO

Axel Brandenburg^{1,2,3} & Evan O’Connor^{2}

1Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-10691 Stockholm, Sweden 2Department of Astronomy, AlbaNova University Center, Stockholm University, SE-10691 Stockholm, Sweden 3McWilliams Center for Cosmology & Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

(Revision: 1.28 )

Draft version July 22, 2020

ABSTRACT

... Subject headings: dynamo — magnetic fields — MHD — turbulence — methods: numerical

1. INTRODUCTION

Magnetic helicity is conserved by the magnetohydrodynamic (MHD) equations in the limit of large magnetic Reynolds numbers. This is a remarkable property of MHD, because something analogous does not exist for the kinetic helicity, which does change—even when the Reynolds number is very large. When the flow is forced with a helical forcing function, it leads to an efficient dynamo process. Since magnetic helicity is conserved, the magnetic field tends to become bihelical with opposite signs of magnetic helicity at large and small length scales, whose contributions cancel.

Numerical codes tend to obey magnetic helicity conservation reasonably well when magnetic energy and helicity dissipation is accomplished through Ohm’s law by microphysical resistivity (Brandenburg & Scannapieco 2020). Magnetic energy dissipation through purely numerical means, on the other hand, yields a magnetic helicity evolution that strongly depends on the details of the numerical scheme. With the eight-wave solver of the FLASH code, for example, spurious magnetic helicity is generated both over long and short time scales (Brandenburg & Scannapieco 2020). By contrast, the SPARK solver in FLASH produces much less spurious magnetic helicity, as will be demonstrated here in detail.

At finite magnetic Reynolds numbers, some level of magnetic helicity is always generated. There has been significant progress in characterizing magnetic helicity production during the late evolution of the dynamo process in its saturation phase (Brandenburg 2001; Field & Blackman 2002; Blackman & Brandenburg 2002), but not much is known about the early evolution before nonlinear effects have played a role. An exception is the work of Brandenburg et al. (2002), who found that for the ABC flow, the normalized magnetic helicity decreases with magnetic Reynolds number Re_{M} like Re_{M}^{−1/2}; see their Figure 5. Assuming Re_{M} to be proportional to the number of mesh points N, this suggests a relatively slow convergence with increasing resolution. It turns out, however, that the actual magnetic helicity produced decreases faster—approximately like N^{−3/2}. Investigating and understanding this in detail is the purpose of the present work.

2. DEFINITION OF THE MODEL

The evolution of the magnetic field B is governed by the induction equation,

\frac{dB}{dt} = ∇ × (U × B − ηJ) ,  (1)

where U(x, t) is the velocity field, which is in general time-dependent, J = ∇ × B/μ_{0} is the current density, with μ_{0} being the vacuum permeability, and η is the magnetic diffusivity. The mean magnetic energy density is given by ⟨B^{2}⟩/2μ_{0}, where angle brackets denote volume averaging. A more convenient quantity is the mean squared magnetic field,

M = ⟨B^{2}⟩,  (2)

as well as the magnetic and normalized current helicities

H = ⟨A · B⟩ and C = ν_{0}⟨J · B⟩.  (3)

We consider a periodic domain of size L^{3}, so the smallest wavenumber is k_{1} = 2π/L. Because of periodicity, there are no magnetic helicity fluxes and the mean magnetic helicity density ⟨A · B⟩ is governed by the equation

\frac{d}{dt}⟨A · B⟩ = −2η⟨J · B⟩  (4)

During the kinematic dynamo phase of interest, the rms magnetic field increases exponentially with the growth rate λ and the mean squared magnetic field ⟨B^{2}⟩ grows at a rate 2λ. It is then useful to introduce normalized quantities

ϵ_{M} = k_{1}⟨A · B⟩/⟨B^{2}⟩, k_{C} = ⟨J · B⟩/⟨B^{2}⟩, (5)

where ϵ_{M} is non-dimensional and characterizes the net magnetic helicity in the system, and k_{C} has the units of a wavenumber that characterizes the magnetic energy spectrum. We then have

2λϵ_{M} = −2ηk_{C}k_{1} (kinematic phase).  (6)

Thus, we see that ϵ_{M} < 0, when the magnetic field grows (λ > 0). It is conceivable that k_{C} increases with Re_{M} to some power n, so we write k_{C} = k_{1}Re_{M}^{n}, where n is an empirically determined parameter. For Kolmogorov turbulence, for example, one expects n = 1/4 (Blackman & Brandenburg 2002). In general, however, we have

ϵ_{M} = −ηk_{C}k_{1}^{2}/λ = −λ̃^{−1}k̃_{f}^{−2}Re_{M}^{−(1−n)},  (7)

Fig. 1.— FLASH .

where λ̃ = λ/u_{rms}k_{f} is the nondimensional growth rate and k̃_{f} = k_{f}/k_{1} is the nondimensional forcing wavenumber. The findings of Brandenburg et al. (2002) are reproduced when n = 1/2. In an ideal code, we expect there to be an effective magnetic Reynolds number proportional to N, so we define

R_{eff} = R_{0}N,  (8)

where R_{0} is the mesh Reynolds number. We have dropped here the specification to the magnetic Reynolds number, because for ideal codes one expects the fluid and magnetic Reynolds numbers to be the same. The empirical finding ϵ_{M} ∝ N^{−3/2} then implies then for n a negative value, n = −1/2, and thus

ϵ_{M} = −λ̃^{−1}R_{eff}^{−3/2}.  (9)

This allows us to compute the effective grid Reynolds number as

R_{0} = −N^{−1} (λ̃ k̃_{f}^{2}ϵ_{M})^{−2/3}.  (10)

Figure 1 shows magnetic helicity production with the SPARK solver in FLASH. The early exponential growth phase of the dynamo corresponds to the time interval from 20 to 140 (in code units). In the following, we average the instantaneous growth rate, λ(t) = d ln B_{rms}/dt and various other normalized quantities such as ϵ_{M} and k_{C}, as well as u_{rms}, over this time interval. Figure 2 shows that the magnetic helicity production decreases with increasing resolution N like N^{−1.5}. (Not clear why.)

In Figure 3 we show the magnetic helicity evolution in DNS for different values of Re_{M}; see also Table 1 for a summary of quantities averaged over the linear growth phase. The data of the DNS correspond to an inhomogeneous dataset because both Pr_{M} and Re_{M} vary. We also note that in the DNS, u_{rms} is only about half the value used in the FLASH code results. Figure 4 shows that the magnetic helicity production decreases with Re_{M} like Re_{M}^{−0.7}, at least for the range of values shown here.

In Figure 5 we combine the FLASH and Pencil Code results. For each value of Re_{M}, we show here separate lines with different resolutions.

Using Equation (10), we compute R_{0}, which is shown in Figure 7 as a function of N.

Fig. 2.— Results with FLASH .

Fig. 3.— Results with Pencil Code.

Fig. 4.— Results with Pencil Code.

Looking at Table 1, we see that k_{C}/k_{f} converges with N more slowly than ϵ_{M}.

3. TWO-SCALE APPROACH

In a closed or periodic domain, magnetic helicity can only be produced by the ⟨J ·B⟩ term, even if ⟨A·B⟩ = 0 initially. To understand this, we now compute ⟨J · B⟩

TABLE 1 FLASH (first three rows) and Pencil Code results.

Run  N k_{f} u_{rms}  ν  η  Pr_{M}  Re  Re_{M} −ϵ_{M}  k_{C}/k_{1}  λ
O  16 16.5 0.405  0  0  –  –  – 2.4e-02 -2.5e-01 7.8e-02
O  32 16.5 0.427  0  0  –  –  – 8.7e-03 4.8e-02 1.2e-01
O  64 16.5 0.408  0  0  –  –  – 2.8e-03 5.4e-02 1.3e-01
I  128 16.5 0.211 5.0e-05 5.0e-05  1 255.3 255.3 3.5e-04 9.4e-01 2.0e-01
I  256 16.5 0.212 5.0e-05 5.0e-05  1 256.5 256.5 3.1e-04 1.1e+00 1.9e-01
I  512 16.5 0.211 5.0e-05 5.0e-05  1 255.8 255.8 3.0e-04 1.1e+00 1.9e-01
P  256 16.5 0.206 2.0e-04 5.0e-06  40 62.5 2501.7 5.8e-05 6.2e-01 3.5e-01
P  512 16.5 0.206 2.0e-04 5.0e-06  40 62.3 2493.6 6.3e-05 4.2e+00 3.5e-01
X  64 28.3 0.085 2.0e-03 5.0e-05  40  1.5  59.8 4.3e-03 1.2e+01 1.5e-01
X  128 28.3 0.085 2.0e-03 5.0e-05  40  1.5  60.1 4.1e-03 1.3e+01 1.6e-01
X  256 28.3 0.084 2.0e-03 5.0e-05  40  1.5  59.6 4.3e-03 1.3e+01 1.5e-01
Y  128 28.3 0.178 4.0e-04 1.0e-05  40 15.7 629.0 1.3e-04 -2.7e+00 4.0e-01
Y  256 28.3 0.180 4.0e-04 1.0e-05  40 15.9 635.6 2.8e-04 9.5e+00 3.8e-01
Y  512 28.3 0.180 4.0e-04 1.0e-05  40 15.8 634.0 2.7e-04 1.0e+01 3.8e-01
Z  128 28.3 0.194 2.0e-04 5.0e-06  40 34.3 1372.4 8.1e-05 -7.7e+00 1.4e+00
Z  256 28.3 0.199 2.0e-04 5.0e-06  40 35.2 1407.5 8.0e-05 5.1e-01 4.9e-01
Z  512 28.3 0.199 2.0e-04 5.0e-06  40 35.2 1407.7 9.3e-05 8.6e+00 4.9e-01
Z 1024 28.3 0.196 2.0e-04 5.0e-06  40 34.5 1380.9 9.5e-05 9.3e+00 4.9e-01

Fig. 5.— Combined presentation of the results both with FLASH and the Pencil Code.

Fig. 6.— Combined presentation of the results both with FLASH and the Pencil Code versus Re_{M}. For the FLASH results, we have assumed Re_{M} = N^{7/3}/2000.

under the two-scale approximation, i.e.,

⟨J · B⟩ = ⟨\overline{J} · \overline{B}⟩ + ⟨j · b⟩  (11)

⟨J · B⟩ = k_{1}^{2}⟨\overline{A} · \overline{B}⟩ + k_{f}^{2}⟨a · b⟩  (12)

Fig. 7.— Effective grid Reynolds number. (Dubious?)

Fig. 8.— Dependence of k_{C} on N.

Making use of H_{M} = ⟨\overline{A} · \overline{B}⟩ + ⟨a · b⟩, we have ⟨\overline{A} · \overline{B}⟩ = H_{M} − ⟨a · b⟩, so

⟨J · B⟩ = k_{1}^{2}(H_{M} − ⟨a · b⟩) + k_{f}^{2}⟨a · b⟩  (13)
⟨J · B⟩ = k_{1}^{2}H_{M} + (k_{f}^{2} − k_{1}^{2})⟨a · b⟩  (14)

\frac{d}{dt}⟨A · B⟩ = −2η⟨J · B⟩  (15)

Fig. 9.— Dependence of k_{C} on Re_{M}.

\frac{dH_{M}}{dt} = −2ηH_{M} − 2η(k_{f}^{2} − k_{1}^{2})⟨a · b⟩  (16)

Using k_{f}^{2}⟨a · b⟩ ≈ ⟨j · b⟩ ≈ μ_{0}ρ_{0}⟨ω · u⟩, we have

\frac{dH_{M}}{dt} = −2ηH_{M} − 2η(1 − k_{1}^{2}/k_{f}^{2})μ_{0}ρ_{0}⟨ω · u⟩ (17)

4. CONCLUSIONS SO FAR

The magnetic helicity produced during the early kinematic phase of a dynamo is being reduced with increasing resolution. Both SPARK and the Pencil Code reduce the magnetic helicity production approximately like N^{−1.5}.

REFERENCES

Blackman, E. G., & Brandenburg, A. 2002, ApJ, 579, 359 Brandenburg, A. 2001, ApJ, 550, 824 Brandenburg, A., & Scannapieco, E. 2020, ApJ, 889, 55 Brandenburg, A., Dobler, W., & Subramanian, K. 2002, AN, 323, Field, G. B., & Blackman, E. G. 2002, ApJ, 572, 685