Zeta expansion for long-range interactions under periodic boundary conditions with applications to micromagnetics

Technical Summary: Zeta Expansion for Long-Range Interactions Under Periodic Boundary Conditions

Problem Statement
The paper addresses the computational challenge of evaluating power-law-based interaction potentials (specifically Riesz potentials of the form $|r|^{-\nu}$) for homogeneous $d$-dimensional bodies interacting with an infinite $n$-dimensional array of copies under periodic boundary conditions (PBC). This problem is central to micromagnetics, where accurate simulation of magnetic textures in wires ($n=1$) or films ($n=2$) requires handling infinite lattice sums of dipolar interactions ($\nu=1$) and their derivatives.

Current standard practices in popular micromagnetic software (e.g., MuMax3, OOMMF) rely on truncating these infinite lattice sums to a finite number of images. This approach introduces uncontrolled errors because the interaction kernels decay algebraically rather than exponentially. Consequently, the periodicity of local interactions is exact, but the non-local interactions are inaccurate. These discrepancies cannot be eliminated by mesh refinement. While Ewald summation exists for point charges, it is not directly applicable to the finite-volume domains and higher-order derivatives required in micromagnetics.

Methodology: The Zeta Expansion
The authors propose a "Zeta Expansion" method that decomposes the infinite lattice sum into two parts:

1. A Small Direct Sum: A truncated sum over a "near-field" lattice ($L_{near}$), similar to standard implementations.
2. A Correction Term: An analytically derived term that eliminates the truncation error.

The core theoretical insight is that the correction term can be expressed as a convergent Taylor series involving derivatives of a generalized set zeta function (a meromorphic continuation of lattice sums). For a compact domain $\Omega$, the potential $U^{(m)}(r)$ is given by:
$$U^{(m)}(r) = \sum_{z \in L_{near}} S^{(m)}(r+z) + D_\Omega Z^{(m)}_{L_{far}, \nu}(r)$$
where $S^{(m)}$ is the direct interaction integral, $Z$ is the set zeta function for the far-field lattice, and $D_\Omega$ is a shape-dependent differential operator defined by the moments of the domain $\Omega$. For axis-aligned cuboids, this operator takes an explicit form involving even-order derivatives.

To evaluate the zeta function and its derivatives efficiently, the authors utilize a Crandall representation. This representation splits the zeta function into:

• A real-space sum involving "upper Crandall functions" ($G_\nu$).
• A reciprocal-space sum involving "lower Crandall functions" ($g_\nu$) and incomplete Bessel functions ($G_\nu(k, r)$).

Both sums in the Crandall representation decay superexponentially, allowing for truncation with negligible error. The authors provide algorithms to compute these special functions (including incomplete Bessel functions) to full machine precision, handling singularities and derivatives analytically via recurrence relations and reflection formulas.

Key Contributions

1. Exact Representation: The paper derives an exact formula for infinite lattice sums of Riesz potentials and their derivatives by complementing a direct sum with a zeta-function-based correction term.
2. Superexponential Convergence: The proposed algorithm achieves machine precision with a computational cost comparable to standard truncated summation schemes, as the correction term converges superexponentially.
3. General Applicability: The framework supports arbitrary interaction exponents $\nu$, arbitrary dimensions $d$ and periodic dimensions $n$ ($n < d$), and arbitrary source/target geometries (including curved boundaries via moment integrals).
4. Special Function Evaluation: The authors provide a robust, superexponentially convergent algorithm for evaluating incomplete Bessel functions and their derivatives, which are critical for the reciprocal-space component of the method.
5. New Asymptotic Corrections: The method yields new corrections to the asymptotic limit of the demagnetization field for finite-size cuboids in periodic arrays, specifically identifying a correction term scaling as $N^{-7}$ for large $N$.

Results and Benchmarks

• Precision: The method consistently achieves full machine precision (relative errors $\approx 10^{-16}$) across various geometries and derivative orders.
• Comparison with Direct Summation: For large interaction exponents ($\nu \ge 10$), the zeta expansion matches direct summation over large grids. For the micromagnetic case ($\nu=1$), the authors demonstrate that standard truncation (e.g., $N_{cut}=10$) results in relative errors exceeding 15%, whereas the zeta expansion remains exact.
• Demagnetization Factors: The authors tabulate high-precision values for the demagnetization factor $D_z$ of 3D cuboids in 2D periodic arrays for various aspect ratios, serving as a reference for future solvers.
• Asymptotic Behavior: The method successfully reproduces known asymptotic limits (e.g., $D_z \to 1/3$ for point dipoles, $D_z = 1$ for continuous layers) and provides the next-order correction term.

Significance and Claims
The authors claim that this work solves the long-standing problem of uncontrolled truncation errors in periodic micromagnetic simulations. By leveraging generalized zeta functions, the method allows for the evaluation of interaction energies and fields to machine precision at a computational cost no greater than current truncated schemes.

The paper positions the method as a "pre-routine" for micromagnetic solvers: the interaction kernels can be precomputed to high precision and then used in standard iterative algorithms (e.g., FFT-based convolution). The authors emphasize that the technique is broadly applicable beyond micromagnetics, citing potential applications in:

• Ferroelectrics: Due to electromagnetic duality, the method applies to polarization textures.
• Atomistic Spin Dynamics: For modeling long-range interactions like RKKY in multilayers.
• Molecular Dynamics: For handling periodic boundary conditions in unbounded dimensions.

The work does not propose new physical phenomena but rather provides a rigorous numerical tool to eliminate numerical artifacts in existing physical models, thereby enabling more reliable studies of systems where minute energy differences (e.g., in frustrated systems) are critical.