An Adjoint Projection Formulation for Enforcing the divergence-free Constraint in Smoothed Particle Magnetohydrodynamics

Technical Summary: An Adjoint Projection Formulation for Enforcing the $\nabla \cdot \mathbf{B} = 0$ Constraint in Smoothed Particle Magnetohydrodynamics

Problem Statement
In magnetohydrodynamics (MHD), the solenoidal constraint $\nabla \cdot \mathbf{B} = 0$ is a fundamental physical requirement. While mesh-based codes often employ constrained transport (CT) to preserve this constraint discretely, meshless methods like Smoothed Particle Magnetohydrodynamics (SPMHD) lack a natural face–edge topology to apply CT directly. Consequently, the mainstream approach in SPMHD has been "divergence cleaning," which evolves an auxiliary scalar field to propagate and damp divergence errors. However, cleaning methods introduce physical parameters (propagation speed and damping timescale) that can affect simulation outcomes if not chosen consistently, particularly in systems with large dynamical ranges. Previous projection attempts in SPMHD, such as those by Price and Monaghan, utilized gravity-like summations for corrections but suffered from inconsistencies between the divergence operator used to measure errors and the gradient operator used to compute corrections, leading to poor performance on short-wavelength noise.

Methodology
The paper proposes a projection method that enforces the discrete divergence-free constraint by solving an elliptic problem after each MHD update. The core innovation is the construction of a correction operator that is strictly the adjoint of the discrete divergence operator used to measure the error, defined within a volume-weighted metric.

1. Formulation: The method seeks a corrected magnetic field $\mathbf{B}^{n+1}$ that minimizes the change in the discrete magnetic energy norm (defined by particle volumes $V_i = m_i/\rho_i$) subject to the constraint $D\mathbf{B}^{n+1} = 0$, where $D$ is the discrete SPH divergence operator.
2. Adjoint Gradient: By introducing a Lagrange multiplier field $\pi$, the correction is derived as $\mathbf{B}^{n+1} = \mathbf{B}^* - \mathcal{G}\pi$. The operator $\mathcal{G}$ is explicitly constructed as the volume-metric adjoint of $D$ ($\mathcal{G} = M_V^{-1} D^T$). This ensures that the correction space is orthogonal to the divergence-free subspace with respect to the volume inner product.
3. Discretization: The paper derives explicit SPH formulas for both the divergence operator $D$ and its adjoint gradient $\mathcal{G}$. Crucially, the adjoint gradient accounts for the asymmetry in SPH derivative coefficients (arising from varying smoothing lengths and densities) by using the specific neighbor-list symmetry and volume weighting.
4. Numerical Solution: The resulting projection equation $L\pi = r$ (where $L = D\mathcal{G}$) forms a symmetric positive semidefinite linear system. This allows the system to be solved efficiently using the Preconditioned Conjugate-Gradient (PCG) method without explicitly assembling the global matrix. The implementation distinguishes between active particles (where $\pi$ is solved) and fixed boundary particles.

Key Contributions

• Consistent Adjoint Formulation: The paper provides the first SPMHD projection formulation where the gradient correction is rigorously constructed as the adjoint of the chosen discrete divergence operator under a volume-weighted metric. This resolves the discretization-level inconsistencies found in prior projection attempts.
• Energy Minimizing Property: The volume metric is shown to be the natural metric for magnetic energy. Consequently, the projection yields an energy-minimizing correction that does not increase the discrete magnetic energy, a critical property for physical fidelity.
• Efficient Solver: The formulation leads to a symmetric positive semidefinite system solvable by PCG, avoiding the $O(N^2)$ cost of direct summation methods used in earlier approaches.
• Practical Stopping Criteria: The authors introduce practical stopping criteria based on the accumulation of divergence error between projection steps, rather than requiring convergence to machine precision for every step, balancing accuracy with computational cost.

Results
The method was tested using two-dimensional Dedner-type divergence tests and three-dimensional magnetized gravitational collapse simulations.

• Divergence Tests: In 2D tests with both perturbed lattices and random particle distributions, the method successfully reduced the discrete divergence residual to the floating-point roundoff level ($\sim 10^{-15}$) with sufficient PCG iterations. Notably, it corrected errors at wavelengths comparable to the smoothing length ($h$), which previous methods struggled with.
• Collapse Simulations: In 3D collapse runs, the projection method significantly outperformed standard hyperbolic divergence cleaning.
  • Error Reduction: With a practical stopping criterion (reducing accumulated error to a fraction of the MHD update error), the projection method suppressed normalized divergence errors well below those obtained with divergence cleaning.
  • Robustness: Varying the projection interval ($N_{interval} = 1, 10, 100$) produced consistent density and plasma-$\beta$ structures. In contrast, the divergence-cleaning run showed quantitative differences in the outflow structure and plasma $\beta$.
  • Cost: The computational cost of the projection was modest, ranging from roughly 1% to 10% of the SPMHD update cost for practical intervals ($N_{interval} \approx 10$). Even in extreme stress tests, the scheme could drive errors to roundoff levels.

Significance and Claims
The paper claims that this adjoint projection formulation is a robust and attractive alternative to divergence cleaning for controlling $\nabla \cdot \mathbf{B}$ errors in SPMHD and related particle or meshless MHD schemes.

Key advantages highlighted include:

1. Parameter Independence: Unlike divergence cleaning, which depends on physical parameters (wave speed, damping scale) that can alter astrophysical results, the projection method's parameters (stopping criteria, interval) are dimensionless and shown to have minimal impact on the solution within the tested ranges.
2. Physical Consistency: The method guarantees that the correction does not inject artificial magnetic energy.
3. Scalability: The approach is not limited to SPMHD; the authors suggest the operator-based projection idea is applicable to other meshless or moving-mesh schemes (e.g., GIZMO, AREPO) with appropriate modifications.

The authors conclude that while the method requires choosing free parameters for termination and frequency, its ability to produce consistent, high-fidelity results without the sensitivity issues of cleaning methods makes it a superior choice for controlling numerical divergence in particle-based MHD.