Conserving Local Magnetic Helicity in Numerical Simulations

Magnetic helicity is a conserved quantity in ideal MHD (Woltjer 1958) and one of the primary ways to quantify the complexity of a magnetic field. In a weakly nonideal turbulent system, it is still almost conserved, in the sense that the turbulent cascade is ineffective at transferring it to resistive scales, even while the magnetic and kinetic energies are dissipated efficiently. Consequently, Taylor (1974) offered the conjecture that a laboratory system would first relax to the minimum energy state with the same total magnetic helicity (see also Taylor 1986). For example, see Yamada (1999) for a review of the laboratory evidence accumulated in the following two decades. We note in particular the work of Ji et al. (1996), which shows that the dynamical transport of magnetic helicity is an important part of this balance. The conservation of magnetic helicity should be even more robust in astrophysical systems where the microscopic resistivity is usually negligible on large eddy scales. Minimizing energy for the same global magnetic helicity drives the creation of large-scale magnetic fields. We can see this in simple terms by noting that the magnetic helicity density has the dimensions of [length] × [energy density] = Wb². For large-length scales, the same magnetic helicity requires the least energy. This same argument shows why magnetic helicity is so resistant to dissipation. It cannot participate in the turbulent cascade, because the energy density corresponding to small-scale eddies decreases with decreasing eddy size. The tendency of robust magnetic helicity conservation to result in an inverse cascade motivates the proposal that the large-scale dynamo (LSD) is driven by magnetic helicity (Vishniac & Cho 2001) through a combination of turbulent transport and the inverse cascade. One of the main challenges in testing this proposal is to examine the flow of magnetic helicity spatially and between scales when numerical simulations tend to dissipate magnetic helicity at significant rates.

In addition to its importance in determining the relaxation of complex systems, magnetic helicity may play a role in the large-scale dynamo process, a point that was first raised by Pouquet et al. (1976). The evolution of the large-scale magnetic field B _l , where the subscript denotes low-pass filtering on the scale l, is given by

$$\begin{eqnarray}&&{\partial }_{t}{{\boldsymbol{B}}}_{{\bf{l}}}={\rm{\nabla }}\times ({{\boldsymbol{v}}}_{{\bf{l}}}\times {{\boldsymbol{B}}}_{{\bf{l}}}-\eta {{\boldsymbol{J}}}_{{\bf{l}}})+{\rm{\nabla }}\times {{\boldsymbol{e}}}_{{\bf{l}}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{e}}}_{{\bf{l}}}\equiv {\left({\boldsymbol{v}}\times {\boldsymbol{B}}\right)}_{l}-{{\boldsymbol{v}}}_{{\bf{l}}}\times {{\boldsymbol{B}}}_{{\bf{l}}}.\end{eqnarray} \tag{ 2 }$$

And Ohm's law is

$$\begin{eqnarray}&&{\boldsymbol{E}}\equiv -{{\boldsymbol{v}}}_{{\bf{l}}}\times {{\boldsymbol{B}}}_{{\bf{l}}}+\eta {{\boldsymbol{J}}}_{{\bf{l}}},\end{eqnarray} \tag{ 3 }$$

where η is the resistivity for a magnetic diffusivity. In a turbulent medium e _l , the electromotive force filtered on the scale l is approximately the large-scale average of the cross product of the turbulent velocity and magnetic fields. Pouquet et al. (1976) pointed out that e _l has contributions from the eddy-scale contributions to the averaged kinetic and current helicities multiplied by B _l , i.e., as a part of the α effect in dynamo theory. Neither helicity is conserved in ideal MHD, but in the Coulomb gauge, the eddy-scale current helicity, j · b, is approximately k² times the eddy-scale magnetic helicity, a · b, where k² is the mean square wavenumber associated with the large-scale eddies. While only the total magnetic helicity is conserved, the transfer of magnetic helicity between large and small scales depends on the parallel component of the electromotive force. This connection was exploited by Gruzinov & Diamond (1996) to show that a dynamo driven by the kinetic helicity saturates due to the local accumulation of magnetic helicity, and it turns off long before the large-scale field reaches equipartition with the turbulent energy density ("α suppression"). Apparently, the conservation of magnetic helicity can have profound implications for the large-scale dynamo.

Given the existence of large-scale magnetic fields, it seems safe to assume that there is some way to avoid the negative conclusion offered by Gruzinov & Diamond (1996). It is possible that the kinematic dynamo operates, perhaps at a reduced efficiency, while the eddy-scale magnetic helicity is ejected from the system (for a discussion of this idea, see Brandenburg et al. 2009). Alternatively, it may be the case that the kinetic helicity is irrelevant and a magnetic helicity flux, induced by rotation, creates local accumulations of eddy-scale magnetic helicity, which then drive a large-scale dynamo as a side effect of being transferred to large-scale fields (Vishniac & Cho 2001). Finally, it is possible that a large-scale dynamo (LSD) could be produced through a combination of shear and fluctuating kinetic helicity (Vishniac & Brandenburg 1997; Mitra & Brandenburg 2012), although the relation between this mechanism and α suppression has never been fully explored.

Regardless of which of these mechanisms dominate in real systems, we can see that accurate simulations of dynamos require careful accounting of the distribution of magnetic helicity in a system. To accurately model Taylor relaxation, one needs only conserve the total magnetic helicity in a simulation box. To model dynamos, including α suppression and the transport of magnetic helicity, one needs to extend this down to scales not much larger than individual eddy scales. Moreover, because the magnetic helicity is a quadratic quantity, the average magnetic helicity on domain scales can arise from correlations on eddy scales. The dissipation of this contribution to the magnetic helicity proceeds at the ohmic dissipation time at the large eddy scale (or a bit faster; see below). Numerical dissipation on this scale has to be slower than the large-scale dynamo growth rate, for an accurate large-scale dynamo simulation. In practice, while this is a trivial threshold for real systems, it represents a substantial challenge for simulations. This is a particularly difficult issue for simulations of magnetorotational instability, where the large-scale eddy size is proportional to the strength of the magnetic field. What is worse is that the MRI operates primarily in the $\hat{r}\hat{\phi }$ plane, and the vertical wavenumber in unstratified simulations tend to the dissipation scale, which guarantees that magnetic helicity will not be conserved for more than an eddy correlation time unless the magnetic Prandtl number is greater than one (see Brandenburg & Subramanian 2005).

This raises the question of whether it might be useful to enforce magnetic helicity conservation in numerical simulations. In many cases, it is possible to enforce the conservation of an important conserved quantity by choosing an appropriate numerical algorithm, but magnetic helicity is nonlocal and no local numerical algorithm conserves it exactly. Evaluating magnetic helicity in simulations of the MRI in a stratified accretion disk, for example, shows that, in a single orbit, the total magnetic helicity in the box can change by a number comparable to a density scale height times the magnetic energy density (see Davis et al. 2010). What price do we pay for adjusting the small-scale magnetic field in real time to compensate for numerical errors? As a general rule, requiring strict conservation of any quantity will prioritize it above other conserved quantities. In particular, adjusting the magnetic field to conserve magnetic helicity will affect the stress–energy tensor, including the energy density. However, in a turbulent simulation, the fluctuating magnetic energy is continuously lost to the turbulent cascade, so energy conservation at small length scales is of secondary importance unless the energy input (or loss) is comparable to the turbulent energy dissipation rate. Any scheme that adjusts local magnetic field strengths will violate flux freezing, but this is already violated by turbulent transport in the turbulent cascade (Eyink et al. 2013). Finally, if there is no turbulence, then there is almost no numerical loss of magnetic helicity and any scheme for correcting the field will have a negligible effect. In principle, there seems to be no fundamental obstacle to implementing a correction scheme to enforce magnetic helicity conservation. Whether or not a particular scheme improves the accuracy of numerical simulations at an acceptable cost depends on the details of the algorithm. In this paper, we will suggest an approach to enforcing magnetic helicity conservation and explore some of the constraints on its use.

In Section 2 of this paper, we will review the definition of magnetic helicity, as well as other relevant quantities, and their evolution in a turbulent medium. We will review how fluctuations in the magnetic helicity are generated and dissipated, including resistive losses. Finally, we will present an algorithm that allows us to impose strict conservation of magnetic helicity with minimal effect on other properties of the magnetic field. In Section 3, we will test the convergence of our algorithm. In Section 4, we discuss the applications and limitations of this tool.

In a highly conducting fluid, the magnetic field evolves as

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{B}}={\rm{\nabla }}\times \left({\boldsymbol{v}}\times {\rm{B}}-\eta {\rm{J}}\right),\end{eqnarray} \tag{ 4 }$$

where the resistive term η J is small. Numerical effects may or may not mimic physical resistivity. We will assume here that they do, but note places where the difference might be important. Consequently, the vector potential evolves as

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{A}}={\boldsymbol{v}}\times {\rm{B}}-\eta {\rm{J}}-{\rm{\nabla }}\phi ,\end{eqnarray} \tag{ 5 }$$

where the electrostatic potential, ϕ, satisfies

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\phi ={\rm{\nabla }}\cdot ({\boldsymbol{v}}\times {\boldsymbol{B}}-\eta {\boldsymbol{J}})\end{eqnarray} \tag{ 6 }$$

in the Coulomb gauge. This gauge choice has the advantage of enforcing a tight correlation between the small-scale contribution to the total current helicity, J · B, and the small-scale contribution to the magnetic helicity. The former is gauge invariant and plays an important role in dynamo theory. We note that other gauge choices have been used for computational convenience or because the magnetic helicity has some compelling interpretation in that gauge (for examples, see Candelaresi et al. 2011; Hubbard & Brandenburg 2011; Del Sordo et al. 2013).

In this gauge, the magnetic helicity flux is

$$\begin{eqnarray}&&{{\boldsymbol{J}}}_{H}\equiv {\boldsymbol{A}}\times ({\boldsymbol{v}}\times {\rm{B}}-\eta {\rm{J}}+{\rm{\nabla }}\phi ),\end{eqnarray} \tag{ 7 }$$

and the magnetic helicity, H, evolves as

$$\begin{eqnarray}&&{\partial }_{t}H=-{\rm{\nabla }}\,\cdot \,{{\boldsymbol{J}}}_{H}-2\eta {\boldsymbol{J}}\cdot {\rm{B}}.\end{eqnarray} \tag{ 8 }$$

When the resistivity is small, its contribution to J _H will be negligible, but its role in breaking magnetic helicity conservation can still be important. We can define the small-scale contribution to the magnetic helicity as h_l ≡ H_l − A _l · B _l so that the evolution equation is

$$\begin{eqnarray}&&{\partial }_{t}{h}_{l}=-{\rm{\nabla }}\cdot {{\boldsymbol{J}}}_{{Hl}}-2\eta ([{\left({\boldsymbol{J}}\cdot {\rm{B}}\right)}_{l}-{{\boldsymbol{J}}}_{l}\cdot {{\boldsymbol{B}}}_{l}]-2{{\boldsymbol{B}}}_{{\bf{l}}}\cdot {{\boldsymbol{e}}}_{{\bf{l}}},\end{eqnarray} \tag{ 9 }$$

where the eddy-scale magnetic helicity flux, J_Hl , is defined, as for the other eddy-scale contributions, by subtracting from the total flux the contribution from large-scale fields, including e_l . The last term on the RHS of Equation (9) describes the transfer of magnetic helicity between small and large scales, and it is critical for understanding how magnetic helicity is dissipated when the resistivity is not exactly zero.

For this purpose, Equation (9) can be generalized by dividing the turbulent cascade into a series of shells in phase space, each with a typical wavenumber amplitude k (Aluie 2017). Then, instead of large- and small-scale contributions to the total helicity on the scale l, we have contributions from each shell, h_k . The transfer of magnetic helicity between large and small scales is still B · e_l , but now e _l is the sum of multiple terms, each given by ( v × b)_kl . The contribution to e _l from the shell k and each term gives the transfer of averaged helicity, H_l , between large scales and the scale given by k⁻¹; see Figure 1.

Zoom In Zoom Out Reset image size

The rate at which magnetic helicity on some large-scale l is dissipated is η〈 J · B 〉_l , or in other words, proportional to the current helicity averaged over the same scale. The large-scale average of the current helicity is the sum of contributions from smaller scales. We can find the expected distribution of current helicity over scales by considering the transfer rate between scales, 2 B _l · ( v × b)_kl , where a positive sign represents the transfer of magnetic helicity to large-scale structure. In a stationary state, this goes to zero as the resistivity goes to zero. The usual expansion of ( v × b)_l to first order in the eddy correlation time (Pouquet et al. 1976) gives four terms: a contribution proportional to large-scale average magnetic helicity embedded in small-scale eddies, h_kl ; a similar contribution due to the kinetic helicity, but with the opposite sign; a term proportional to the gradient of the large-scale field—usually written as a turbulent diffusivity; and a drift term of the form V _eff × B_l . Obviously, the last magnetic helicity between scale terms orthogonal to B _l does not contribute to the transfer of s. The first term leads to the transfer of h_kl to larger scales at a rate comparable to ${\tau }_{k}^{-1}$, the small-scale eddy turnover rate (Vishniac & Cho 2001). The second term depends on kinetic helicity, which is not a conserved quantity. Naively, one would expect it to be embedded in small scales only to the extent that the symmetry breaking of the large-scale environment biases the small-scale turbulence. However, the current helicity embedded on small scales induces a kinetic helicity, the effect of which is to drive the second term to cancel a fraction of the first term (Vishniac & Shapovalov 2014). Finally, the third term contributes an effect that is approximately −〈(v_i v_j + b_i b_j 〉_{l ist} B_ls ∂_j B_lt τ_k . (Here, we have neglected a contribution due to pressure in the force equation, which somewhat reduces the contribution from the correlated small-scale magnetic fields). This term is of order −〈v²〉_L τ_k times the large-scale contribution to the current helicity. This is the term that has to balance the first two terms. Because the small-scale contribution to the current helicity is k² h_kl , the ratio of the small-scale current helicity to the large-scale current helicity is $\,{({{kv}}_{k}{\tau }_{k})}^{2}$, i.e., of order unity. The current helicity, and consequently the ohmic dissipation of magnetic helicity, is distributed across the full range of the inertial cascade with a one-dimensional spectral index of −1. This in turn implies that the dissipation of magnetic helicity will be enhanced by some logarithmic factor above the ohmic dissipation rate on the large eddy scale. Simulations that use artificial forms of resistivity that depend on higher powers of the wavenumber will damp magnetic helicity at slower rates because of the concentration of dissipation at the smallest eddy scales and because those eddies will be smaller. However, the reduced range of scales subject to dissipation will tend to accumulate extra power, and presumably magnetic helicity, due to the bottleneck effect, which will work in the opposite direction.

In a turbulent cascade, the energy dissipation rate is the energy density times the large-scale eddy turnover rate, i.e., the characteristic timescale is just L/v_L . We have shown that the corresponding time for the dissipation of magnetic helicity is L²/η (divided by some logarithmic measure of the dynamic range of the turbulent cascade). Dividing the former by the latter gives us a measure of how well magnetic helicity is conserved on dynamic timescales. This is inversely proportional to the magnetic Reynolds number V_L L/η, which in most cases means that magnetic helicity is conserved to a level that cannot be adequately reproduced while relying only on numerical precision.

There is one final point worth stressing about the dissipation of magnetic helicity at different scales within the turbulent cascade. We can derive a maximum dissipation rate without making an explicit calculation of the transfer of magnetic helicity between scales. On dimensional grounds, the maximum current helicity at a given wavenumber will be $\sim {{kb}}_{k}^{2}$. However, highly magnetized turbulence will be highly anisotropic, and the current helicity depends on symmetry breaking in all three directions. Consequently, the actual limit is ${k}_{\parallel }{b}_{k}^{2}$. If we multiply this by the Alfvén speed, we get the turbulent energy density times the Alfvén frequency for typical eddies with wavenumber k. Invoking the condition of critical balance, ${\tau }_{\mathrm{nonlinear}}^{-1}\sim {\omega }_{A}$, which is a generic feature of MHD turbulence models, we see that the maximum average current helicity embedded in eddies of wavenumber k is proportional to the energy cascade rate at those scales. From the conservation of energy and locality of nonlinear transfer of energy in phase space, this shows that the maximum dissipation rate of magnetic helicity is constant within the turbulent cascade. In other words, the rate we derived in the previous paragraph, using a detailed but approximate analysis of the transfer of magnetic helicity in phase space, is roughly the same as the maximum rate we could expect in any model of MHD turbulence.

The problem we wish to address in this paper is that the distribution of magnetic helicity will diverge from the distribution given by evolving Equations (7) and (8) with η = 0, due to small numerical errors induced by discreteness effects and by the necessity of using η > 0 (or its equivalent) to promote numerical stability. This leads to a magnetic field B _calc that differs from the ideal result on small scales and whose associated magnetic helicity is not conserved (Bodo et al. 2017). Of course, we can always redefine our goal as a simulation with significant resistivity, but this fails to accurately model the dynamics of realistic astrophysical systems. We can define our error, i.e., the amount by which we have failed to treat the magnetic helicity as a locally conserved quantity, by calculating ∇ · J _H at every time step and advancing H accordingly. Tracking magnetic helicity in a simulation requires calculating ∇ · J _H , which implies a knowledge of A and ϕ everywhere. If the underlying code is pseudospectral or is at least based on a regular grid, then this is straightforward (Squire & Bhattacharjee 2016). Otherwise, the code has to be supplemented with routines, based on Green's function, that will calculate these quantities (Yousef et al. 2008; Singh & Jingade 2015).

Here, we propose a rigorous procedure that will bring the magnetic helicity back into alignment with its expected value while minimizing the size of the correction. We could write this as a straightforward Gaussian minimization problem by defining an action dependent on the mean square correction to the field, subject to constraints that ensure that the corrected field will have the correct magnetic helicity distribution and will satisfy the usual constraints on the magnetic field. This is

$$\begin{eqnarray}&&S\equiv \int [(1/2){\left(\delta {\boldsymbol{B}}\right)}^{2}+\lambda ({\boldsymbol{x}})(H-{\boldsymbol{A}}\cdot {\rm{B}})+\mu ({\boldsymbol{x}}){\rm{\nabla }}\cdot {\boldsymbol{B}}]{dV},\end{eqnarray} \tag{ 10 }$$

where H is the magnetic helicity distribution that we expect, and A and B are the corrected values of the vector potential and magnetic field. The last term in the integral enforces the usual zero divergence condition on the magnetic field. Without this term, the corrected B would no longer satisfy Maxwell's equations. It is not necessarily important to obtain the correct distribution of H on all scales. Just correcting it on scales comparable to, and larger than, the large-scale eddy size would improve the accuracy of dynamo simulations. Random fluctuations in the magnetic helicity flux produce random fluctuations in the amplitude of magnetic helicity modes with scales comparable to the large-scale eddies or smaller. The large-wavelength tail of these fluctuations may be interesting, but small errors in the amplitude of eddy-scale helicity fluctuations are unlikely to have any dynamical significance. Consequently, we are primarily concerned with correcting the distribution of magnetic helicity at wavenumbers significantly smaller than the inverse of a large eddy scale.

Minimizing the change to the magnetic field is not the only choice. It is more convenient and more physically reasonable to minimize the change to the vector potential. The same change in the vector potential at small scales produces a larger change in B, so the net effect is that this choice puts more of the adjustment in the magnetic field at small scales. (This is not equivalent to correcting the magnetic helicity distribution at small scales. The magnetic field on small scales contributes to the large-scale distribution of magnetic helicity.) This is a better choice because we are trying to counteract the effects of inaccuracies at small scales. In addition, an algorithm that compensates for lost magnetic helicity by altering the large-scale magnetic field runs the risk of driving a large-scale dynamo by fiat. Conversely, fixing errors in the magnetic helicity distribution by making correlated changes on small scales implies that the bulk of the corrections will occur on scales within the inertial range of the turbulent cascade. On these scales, magnetic helicity is passed between scales at eddy turnover rates and any excess energy is quickly dissipated. If we weight the algorithm by ∣δ A ∣², then the action we choose to minimize is

$$\begin{eqnarray}&&S\equiv \int [{(1/2)(\delta {\boldsymbol{A}})}^{2}+\lambda ({\boldsymbol{x}})(H-{\boldsymbol{A}}\cdot B)-\mu ({\boldsymbol{x}}){\rm{\nabla }}\cdot {\boldsymbol{A}}]{dV}.\end{eqnarray} \tag{ 11 }$$

The second Lagrangian constraint now fixes the gauge of A , because the usual divergence-free condition for B is automatically satisfied. Integrating this by parts involves a surface term

$$\begin{eqnarray}&&{\rm{\Delta }}S=\int \left[\lambda {\boldsymbol{A}}\times \delta A-\mu \delta A\right]\cdot {dS}.\end{eqnarray} \tag{ 12 }$$

In a finite computational volume, this gives us a pair of boundary conditions, and λ and μ vanish nearly everywhere on the boundary, the only exception being places where A is normal to the surface (for λ). This guarantees that there is no nontrivial solution for λ unless the calculated distribution of magnetic helicity is nonzero somewhere. Here, we will assume a periodic box and ignore this constraint. We obtain

$$\begin{eqnarray}&&\delta {\boldsymbol{A}}=2\lambda {\boldsymbol{B}}+{\rm{\nabla }}(\lambda )\times {\boldsymbol{A}}-{\rm{\nabla }}\mu ,\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&H={\boldsymbol{A}}\cdot B.\end{eqnarray} \tag{ 14 }$$

The second equation is just a restatement of our goal, that the field quantities should give the correct distribution of magnetic helicity. The first equation is not as simple as it appears, because the fields on the RHS are the corrected fields, not the calculated ones. The last term on the RHS arises from the Coulomb gauge and leads to

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}\mu ={\rm{\nabla }}\lambda \cdot {\boldsymbol{B}}={\rm{\nabla }}\cdot (\lambda {\boldsymbol{B}}).\end{eqnarray} \tag{ 15 }$$

In order to make this problem more tractable, and to reduce the numerical cost of implementing a correction scheme, we will assume that λ is small and that we can therefore ignore the difference between the corrected and computed fields on the RHS of the first equation. This should be reasonable if the correction scheme is applied often enough. If we further assume that the correction to the magnetic helicity can be approximated to the same order, i.e., linear in λ, then we have a second-order differential equation for λ(x) in the computational box. With some work, this can be shown to be

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}H & \equiv & H-({{\boldsymbol{A}}}_{\mathrm{calc}}\cdot {{\boldsymbol{B}}}_{\mathrm{calc}})\\ & \approx & \delta {\boldsymbol{A}}\cdot {\boldsymbol{B}}+{\boldsymbol{A}}\cdot ({\rm{\nabla }}\times \delta {\boldsymbol{A}})\\ & = & 2\lambda {B}^{2}+{\boldsymbol{B}}\cdot ({\rm{\nabla }}\lambda \times {\boldsymbol{A}})-{\boldsymbol{B}}\cdot {\rm{\nabla }}\mu +2\lambda {\boldsymbol{J}}\cdot A\\ & & +\,2{\boldsymbol{A}}\cdot ({\rm{\nabla }}\lambda \times {\boldsymbol{B}})+{A}_{i}{A}_{j}{\partial }_{i}{\partial }_{j}\lambda \\ & & -\,\displaystyle \frac{1}{2}{\rm{\nabla }}\lambda \cdot {\rm{\nabla }}({A}^{2})-{A}^{2}{{\rm{\nabla }}}^{2}\lambda \\ & = & 2\lambda ({B}^{2}+{\boldsymbol{A}}\cdot J)-{\rm{\nabla }}\lambda \cdot ({\boldsymbol{A}}\times {\rm{\nabla }}\times A)\\ & & -\,\displaystyle \frac{1}{2}{\rm{\nabla }}\lambda \cdot {\rm{\nabla }}({A}^{2})\\ & & +\,({A}_{i}{A}_{j}-{A}^{2}{\delta }_{{ij}}){\partial }_{i}{\partial }_{j}\lambda \\ & = & 2({B}^{2}+{\boldsymbol{A}}\cdot J)\lambda +{\partial }_{i}\left(({A}_{i}{A}_{j}-{A}^{2}{\delta }_{{ij}}){\partial }_{j}\lambda \right)-{\boldsymbol{B}}\cdot {\rm{\nabla }}\mu ,\end{array}\end{eqnarray} \tag{ 16 }$$

where we employ implicit summation over repeated indices. The last term can be written as

$$\begin{eqnarray}&&{\boldsymbol{B}}\cdot {\rm{\nabla }}\mu ={\boldsymbol{B}}\cdot {\rm{\nabla }}{{\rm{\nabla }}}^{-2}{\rm{\nabla }}\cdot (\lambda {\boldsymbol{B}}).\end{eqnarray} \tag{ 17 }$$

The operator ∇∇⁻²∇ is just $\hat{k}\hat{k}$ in Fourier space. In real space, it is an integral operator acting on everything to its right and can be written

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}{{\rm{\nabla }}}^{-2}{\rm{\nabla }}\cdot ({\boldsymbol{F}}) & = & \displaystyle \frac{{\boldsymbol{F}}}{3}+\,\displaystyle \int \displaystyle \frac{{d}^{3}{\boldsymbol{r}}^{\prime} }{4\pi {r}^{{\prime} 5}}\left(3{\boldsymbol{r}}^{\prime} {\boldsymbol{r}}^{\prime} \cdot F(r+r^{\prime} )\right)\\ & & -\,\left.{r}^{{\prime} 2}{\boldsymbol{F}}({\boldsymbol{r}}+{\boldsymbol{r}}^{\prime} )\right),\end{array}\end{eqnarray} \tag{ 18 }$$

where F is any differentiable vector field. As before, we are assuming that boundary conditions can be ignored. To linear order, we obtain λ from

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}H & = & 2({B}^{2}+{\boldsymbol{A}}\cdot J)\lambda +{\partial }_{i}\left(({A}_{i}{A}_{j}-{A}^{2}{\delta }_{{ij}}){\partial }_{j}\lambda \right)\\ & & -\,{\boldsymbol{B}}\cdot {\rm{\nabla }}{{\rm{\nabla }}}^{-2}{\rm{\nabla }}(\lambda {\boldsymbol{B}}),\end{array}\end{eqnarray} \tag{ 19 }$$

and we use this in Equation (13).

We need a general procedure for solving Equation (19). We can do this through iteration. If we split two of the coefficients on the RHS into varying and constant pieces defined by

$$\begin{eqnarray}&&{C}_{0}\equiv \langle {B}^{2}+{\boldsymbol{J}}\cdot {\boldsymbol{A}}\rangle =2\langle {B}^{2}\rangle ,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{C}_{1}\equiv {B}^{2}+{\boldsymbol{J}}\cdot {\boldsymbol{A}}-{C}_{0},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{D}_{0{ij}}\equiv \langle {\delta }_{{ij}}{A}^{2}-{A}_{i}{A}_{j}\rangle ,\end{eqnarray} \tag{ 22 }$$

and

$$\begin{eqnarray}&&{D}_{1{ij}}\equiv {\delta }_{{ij}}{A}^{2}-{A}_{i}{A}_{j}-{D}_{0{ij}},\end{eqnarray} \tag{ 23 }$$

then we can rewrite Equation (19) as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}H & - & 2{C}_{1}\lambda +{\partial }_{i}{D}_{1{ij}}{\partial }_{j}\lambda +{\boldsymbol{B}}\cdot {\rm{\nabla }}\left({{\rm{\nabla }}}^{-2}{\rm{\nabla }}\cdot (\lambda {\boldsymbol{B}})\right)\\ & & =\,2{C}_{0}\lambda -{\partial }_{i}{D}_{0{ij}}{\partial }_{j}\lambda .\end{array}\end{eqnarray} \tag{ 24 }$$

The fourth term on the LHS of Equation (24) resists any simple reduction, because the operator B ∇⁻² B always depends on the wavenumber of λ( k ). In Fourier space, the RHS of this equation is a simple multiplication by the positive definite coefficient, which is always greater than 2C₀:

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{n+1}({\boldsymbol{k}}) & = & \displaystyle \frac{1}{2{C}_{0}+{k}_{i}{k}_{j}{D}_{0{ij}}}\left[{\rm{\Delta }}H({\boldsymbol{k}})+\left(-2{C}_{1}{\lambda }_{n}+{\partial }_{i}{D}_{1{ij}}{\partial }_{j}{\lambda }_{n}\right.\right.\\ & & +\,\left.\left.{\boldsymbol{B}}\cdot {\rm{\nabla }}\left({{\rm{\nabla }}}^{-2}{\rm{\nabla }}\cdot ({\lambda }_{n}{\boldsymbol{B}})\right)\right)({\boldsymbol{k}})\right].\end{array}\end{eqnarray} \tag{ 25 }$$

The choice of λ₀ is not dictated by this method. One possibility is to use 〈ΔH〉/2C₀ as λ₀ so that

$$\begin{eqnarray}&&{\lambda }_{1}=\displaystyle \frac{1}{2{C}_{0}+{k}_{i}{k}_{j}{D}_{0{ij}}}\left[{\rm{\Delta }}H({\boldsymbol{k}})-{C}_{1}\langle {\rm{\Delta }}H\rangle /{C}_{0}\right].\end{eqnarray} \tag{ 26 }$$

Alternatively, we could choose

$$\begin{eqnarray}&&{\lambda }_{0}=\displaystyle \frac{{\rm{\Delta }}H}{2({C}_{0}+{C}_{1})}.\end{eqnarray} \tag{ 27 }$$

The former captures the direct response to a helicity error at a wavenumber k and a small part of the indirect response through the interaction with the spatial variation in the coefficient of λ. The latter captures the nonlinear part a bit better but misses the dependence on gradients of λ.

We can get a sense of how important the various contributions to λ are by considering two simple cases. It is worth noting that, due to random fluctuations in the magnetic helicity flux, we expect random fluctuations in the amplitude of modes with scales comparable to the large-scale eddies or smaller. Small errors in the amplitude of these fluctuations are unlikely to have any dynamic significance. Consequently, we are primarily concerned with correcting the distribution of magnetic helicity at wavenumbers significantly smaller than the inverse of a large eddy scale. For simplicity, we will assume that ΔH( k ) is nonzero only for these modes larger than the large-scale eddy modes. We expect the magnetic energy density will be largest at the large eddy scale, and the rms vector potential will be either largest at that scale or the large-scale magnetic field domain scale. In the former case, we have ${C}_{0}\sim {B}_{T}^{2}\sim {D}_{0{ij}}{k}_{T}^{2}$, where k_T is the wavenumber of the large-scale eddies. The fluctuating coefficients C₁ and D_1ij will be functions of scale, but at the large eddy scale, both are comparable to their nonfluctuating pieces. We see from Equation (24) that the contribution to λ from the large eddy scale will be comparable to that from the scales where ΔH is nonzero. On smaller scales, l << l_T dominant term in C₁ will be A · J, which will be of order A_T A_l /l² ∝ l^−2/3 for the Goldreich–Sridhar model of turbulence. Consequently, the contribution to λ from scales in the inertial cascade will decrease as l^4/3.

The magnetic domain scale vector potential can be larger than A_T even in the early stages of the large-scale dynamo. In this case, the previous argument goes through in a slightly different form. Now the contribution to λ from the turbulent scale is reduced by a factor of A_T /A₀ relative to the scale where ΔH is nonzero. As before, the contribution to λ falls off as l^4/3 on smaller scales.

We conclude that solutions for λ need to be calculated to scales somewhat smaller than the large eddy scale in order to recover the desired distribution of magnetic helicity on magnetic domain scales.

After solving for λ to the required accuracy, we correct the magnetic potential:

$$\begin{eqnarray}&&\delta {\boldsymbol{A}}={\rm{\nabla }}\times (\lambda {\boldsymbol{A}})+\lambda {\boldsymbol{B}}-{\rm{\nabla }}{{\rm{\nabla }}}^{-2}{\rm{\nabla }}\cdot (\lambda {\boldsymbol{B}}).\end{eqnarray} \tag{ 28 }$$

In phase space, the last term does not have to be evaluated explicitly; rather, the sum of the first two terms can simply be projected onto the plane perpendicular to k . If the initial error in the magnetic helicity distribution is large, then it may be necessary to iterate to allow for the nonlinear term (∇ × δ A) · δ A.

If we neglect the gradient of λ and consider the correction to magnetic helicity, we find that

$$\begin{eqnarray}&&{\rm{\Delta }}H\approx {\boldsymbol{B}}\cdot (2\lambda {\rm{B}})+{\boldsymbol{A}}\cdot ({\rm{\nabla }}\times 2\lambda {\rm{B}})=2\lambda ({{\rm{B}}}^{2}+{\boldsymbol{A}}\cdot {\rm{J}}).\end{eqnarray} \tag{ 29 }$$

The reconstituted magnetic helicity is not deposited strictly according to the energy density in the turbulent cascade, but on average it works out to the same thing as a function of scale. This implies that the power spectrum of the deposited magnetic helicity is much shallower than we would expect to find in a turbulent cascade. This is a consequence of Equation (28). There is some minimum scale, or maximum wavenumber, where the correction to the power spectrum is of order unity, ${k}_{\max }\lambda \sim 1$. In order to avoid highly unphysical effects on the small-scale turbulent spectrum, we require that the dissipation scale wavenumber satisfy

$$\begin{eqnarray}&&1\gt {k}_{\mathrm{diss}}\lambda \sim \displaystyle \frac{{k}_{\mathrm{diss}}{\rm{\Delta }}H}{4\langle {B}^{2}\rangle }.\end{eqnarray} \tag{ 30 }$$

We can use this to set a limit on the cadence of corrections, because we need

$$\begin{eqnarray}&&{\rm{\Delta }}H\lt 4\langle {B}^{2}\rangle {l}_{\mathrm{diss}}.\end{eqnarray} \tag{ 31 }$$

It is not immediately obvious that the procedure described in Equation (24) will converge efficiently—or at all. We have performed a simple test using a time slice from a simulation of MHD turbulence in the turbulence database of The Johns Hopkins University (Li et al. 2008; Eyink et al. 2013). In this simulation, the turbulence is driven at wavenumbers of 2 and contains an inertial range extending roughly from wavenumbers of 10¹–10². Rather than tracking the magnetic helicity over time and using our algorithm to correct to a predicted magnetic helicity distribution, we deliberately distort the magnetic field and treat the new magnetic helicity distribution as an error that the algorithm will correct. To do this, we take each Fourier mode in our randomly chosen time slice and change the phase randomly, with an rms value of 0.5. For computational convenience, we degrade the resolution by a factor of 2–512³. This eliminates about half of the strongly damped regime but leaves the inertial cascade untouched. We can check the scale distribution of helicity dissipation by plotting the one-dimensional spectrum of the current helicity or J ( − k ) · B ( k )k². We see a bump at small wavenumbers and rough agreement with our predicted slope from wavenumbers ∣ k ∣ ∼ 20 to well into the damped regime. We note that this implies that some of the dissipation is occurring well within the damped regime where the magnetic spectrum is sharply declining.

One complication is that the gauge term given in Equation (17) is nonlocal and potentially expensive to evaluate. In Fourier space, the central operator has a simple form, but calculating nonlinear terms in real space means that we need to perform two fast Fourier transforms every time we evaluate this expression. Instead, we use a local estimate of this term, following the method described by Wagner et al. (2017).

Applying the algorithm to the altered map and using the original helicity distribution as the desired goal, we obtain convergence to the limit of numerical accuracy for λ at all scales in between 6 and 30 time steps. In Figure 2, we show a comparison between the initial distribution of helicity and the magnetic potential as a function of wavenumber and the final values after λ converged via the repeated application of Equation (25). These plots show amplitudes and not phase information, but because phase changes lead to amplitude changes in the helicity, this is an adequate proxy for the ability of the algorithm to recover the magnetic helicity. From the second panel, we see that we recover the original distribution of magnetic helicity without introducing significant errors in the vector potential.

Zoom In Zoom Out Reset image size

Finally, we performed a test with a smaller shift in the phase of the modes (an rms phase shift of 0.05). In this case, the change in the magnetic helicity—and subsequent corrections—was small at higher wavenumbers, although they were still substantial at large wavelengths. This is illustrated in Figure 3.

Zoom In Zoom Out Reset image size

Magnetic helicity is an important topological quantity that constrains large-scale relaxation of magnetized systems and plays an important role in large-scale dynamos. Because the rate at which it is dissipated in a turbulent cascade is smaller than the large-scale turnover rate (or the energy dissipation rate), requiring a simulation to treat it as a locally conserved quantity is roughly as important as requiring mass conservation. In this study, we have constructed an algorithm for correcting errors in the magnetic helicity distribution in numerical simulations. For the purposes of this work, we define an error as any deviation between the evolution of the simulated system and one in which magnetic helicity is a locally conserved quantity, as it typically is in the systems being modeled. Numerical simulations that include turbulence frequently show a loss in global magnetic helicity. There is no published information on the ability of current algorithms to track magnetic helicity correctly on smaller scales, but it is reasonable to suppose that they do significantly worse at this more difficult task. While there is no unique way to restore lost information, the algorithm explored here has some attractive features, including the ability to embed the lost helicity to roughly the same scales where numerical errors eliminated it. Adopting this algorithm has two associated costs. The most obvious is the computational costs involved with the algorithm itself. However, this method cannot be applied without tracking local magnetic helicity in a simulation, i.e., without the use of Equations (7) and (8) with η = 0. The latter implies a fractional increase in the computational costs of every time step. The former is hard to predict without knowing the necessary cadence of corrections, but it is probably less significant. We note also that, in the event that the resistivity is large enough that its role in the system evolution is not entirely swamped by numerical noise, this algorithm will still function properly as long the resistivity is included in Equations (7) and (8).

We have examined the cleaning algorithm for magnetic helicity in the JHU 1024³ MHD turbulence simulation. The algorithm restored the large-scale distribution of magnetic helicity after it was distorted by randomly altering the phases of Fourier components of the magnetic field. The iterative procedure for constructing the solution converged, despite a lack of formal proof of convergence. As expected, the algorithm does increase the amount of very small-scale power in the turbulent cascade, a result that needs to be considered in setting the cadence of correction in a simulation. We also note that, while magnetic helicity is robustly conserved with any gauge choice, our work is aimed specifically at enforcing conservation in the Coulomb gauge. We justify this choice by pointing to the close correlation between the magnetic helicity in this gauge and the current helicity, which appears in analytic theories of the large-scale dynamo (and is not gauge dependent).

One complication is that we have not formally demonstrated that our procedure for deriving λ will always converge. We have shown, as a practical matter, that it converges in our test case, a simulation of homogeneous isotropic magnetized turbulence. This is a sufficiently robust test that it is plausible that it will always converge. However, we note that the fluctuating parts of the coefficients in the second-order equation for λ, Equation (19), will be roughly as large as the constant parts. Moreover, the gauge constraint term will typically be of order λ B², i.e., the same order as the first RHS term. We have already seen that, in real space, this term is the sum of λ B²/3 plus a local quadrupole contribution that we can reasonably expect to be smaller. We can compare it to the first RHS term, whose average value is 4B² λ. We see that the gauge constraint term will be smaller by an order of magnitude than the other terms on the RHS.

The difficulty of enforcing strict conservation of magnetic helicity in MHD simulations is a major obstacle in performing realistic simulations of large-scale dynamo processes in astrophysical objects. This is particularly true for magnetorotational instability, where the large-scale eddy size can vary dramatically over time. More generally, magnetic helicity conservation may play a key role in the shear-current effect (Squire & Bhattacharjee 2016). Moreover, for small-scale MHD fluctuations, this effect may define the off-diagonal turbulent diffusivity η_xy (Singh & Jingade 2015, 2015; Squire & Bhattacharjee 2016; Dewar et al. 2020; Jingade & Singh 2021). As a result, this could be the crucial bound of for MHD instabilities. Our results show that a viable correction scheme exists and can be used to produce a realistic magnetic helicity distribution. In future work, we will include this cleaning algorithm in a stratified periodic shearing-box simulation (see, for example, Davis et al. 2010) to follow the flow of magnetic helicity and its effect on the magnetorotational instability-driven dynamo.

This work makes use of the Johns Hopkins Turbulence Database. The authors wish to recognize and acknowledge helpful discussions with Amir Jafari and Greg Eyink. E.T.V. thanks the AAS for supporting the research. Y.Z. thanks Mor Rozner for stimulating discussions. Y.Z. thanks the CHE Israeli Excellence Fellowship for International Postdoctoral Researchers for supporting the research. Both authors wish to thank the referee, whose numerous questions led us add detailed explanations of why magnetic helicity conservation is an important goal. We believe this has led to a paper that is much more accessible to readers.