Testing magnetic helicity conservation in ideal 3D MHD simulations

[Kalfeu]

Hi everyone,

I have a question regarding the integral conservation of magnetic helicity $H_m=∫A \cdot B \mathrm dV$ (where $A$ is the magnetic vector potential) in ideal MHD (without Ohmic resistivity) simulations using Athena++.

In Cartesian coordinates, the constrained transport algorithm preserves $\nabla \cdot B=0$ to machine precision. For smooth solutions and suitable boundary conditions (e.g., perfectly conducting walls), magnetic helicity is theoretically conserved in the ideal limit.

I am particularly interested in whether anyone has tested or benchmarked the conservation of $H_m$ in different 3D cases (2D being trivial), in particular cylindrical and spherical geometries. For example, using simple smooth or non-smooth field configurations, how well is magnetic helicity conserved compared with known analytic solutions? Any published references, code experiences, or known caveats would be greatly appreciated.

Two points that might be relevant for such tests:

1.Ambipolar diffusion and the Hall effect produce electric fields perpendicular to $B$; therefore, in the absence of Ohmic resistivity, they do not directly affect magnetic helicity conservation. If significant non‑conservation is observed, it should not be attributed to these effects but rather to numerical or boundary issues.

2.The spatial reconstruction method and the Riemann solver can indirectly influence magnetic helicity, because the evolution of $B$ is coupled to the velocity field through the induction equation. Discretization errors from these choices may accumulate and affect the integral.

Thank you.

[Kalfeu]

Let me add some further information to clarify the motivation behind this question.

1. What magnetic helicity really measures
Magnetic helicity $H_m = \int \mathbf{A}\cdot\mathbf{B} dV$ quantifies the topology of the magnetic field – how much the field lines are linked or knotted. For example, two linked flux rings give non‑zero helicity; unlinked ones give zero. In ideal MHD (zero resistivity) and with perfectly conducting walls ($\mathbf{B}\cdot\mathbf{n}=0,\mathbf{E}\times\mathbf{n}=0 $), $H_m$ is exactly conserved.

2. The helicity–energy inequality (Arnold 1974)
A very important relation is the Arnold inequality:

$\left| \int_{\Omega} \mathbf{A}\cdot\mathbf{B}dV \right|\le C_p \int_{\Omega} B^2dV$,

where $C_p$ is the Poincaré constant. This is a simple consequence of Cauchy–Schwarz and Poincaré inequalities, yet its physical meaning is profound:

• In ideal MHD, if the initial helicity is non‑zero, it stays constant in time.
• The magnetic energy cannot decay below the lower bound set by helicity; hence, even with dissipation, the system relaxes to a non‑zero energy steady state – typically a force‑free field. This is the foundation for understanding solar coronal magnetic fields and many astrophysical dynamos.

3. Numerical consequences: why the discretisation matters
From a numerical perspective, preserving helicity (or not) has dramatic effects on long‑time simulations.

• A helicity‑preserving (HP) scheme exactly conserves $H_m$ at the discrete level. Consequently, energy decays only down to the theoretical lower bound, giving a physically correct non‑zero steady state.
• A standard scheme that does not enforce helicity conservation suffers from numerical dissipation and rounding errors. The field topology is artificially broken (numerical reconnection), energy loses its topological lower bound, and eventually decays to an unphysical zero state.

This is not just a subtlety – numerical reconnection can spuriously release magnetic energy, heat the plasma, and trigger artificial instabilities, completely altering the solution.

4. A subtle example: magnetic braids with zero total helicity
Consider a magnetic braid composed of six strands: three positive and three negative interlinkings, so that the global helicity is zero.

• If we only enforce global helicity conservation, local braiding can reconnect and straighten out.
• But if we could enforce local helicity conservation (i.e., on a sub‑domain whose boundary satisfies $B\cdot n=0$), the braids remain tangled.

Which one corresponds to real physics? This leads to two important issues:

• Magnetic helicity is only a second‑order invariant. It does not capture all topological details. For the braid example, the topology is non‑trivial yet the total helicity vanishes. Higher‑order topological invariants are needed.
• Ideal MHD’s frozen‑in theorem implies infinitely many pointwise topological constraints – magnetic field lines cannot reconnect at all. In reality, however, a tiny resistivity or viscosity always exists. J. B. Taylor made a bold assumption: global helicity remains nearly conserved even while local reconnection occurs (Taylor relaxation). This assumption successfully explains the reversed‑field pinch phenomenon in fusion devices.

5. Why this matters for Athena++ (and any MHD code)
Numerical simulations are only trustworthy if we understand the intrinsic conservation properties of the discretisation. As the braid example shows, the choice between global‑only vs. local‑preserving (or high‑order) schemes can lead to completely different physical outcomes – false reconnection vs. realistic topological evolution.

I opened this discussion to understand how well Athena++ (especially in curvilinear coordinates) respects magnetic helicity conservation out‑of‑the‑box. If the code already does a good job, great – we can use it with confidence. If not, perhaps this thread can serve as a starting point for requesting new features or more advanced discretisations (e.g., structure‑preserving schemes) from the developers.

I’d love to hear your experiences, references, or ideas on:

• How to test helicity conservation in practice?
• Any known schemes that preserve local helicity or higher‑order invariants?
• Lessons from Taylor relaxation in numerical MHD.

Looking forward to your thoughts!