PaperlayerView exact revision on arXiv

Source-linked AI summary

Computation of anisotropic singular sums from high-order derivatives of Epstein zeta functions

Andreas A. Buchheit, Jonathan K. Busse

arXiv:2609.28282v1math.NAcond-mat.str-elcs.NAmath.NT

TL;DR

Stable computation of anisotropic power-law lattice sums has lagged behind isotropic methods, despite its relevance to long-range interacting systems. The paper develops anisotropic Epstein zeta representations and an algorithm based on harmonic decomposition. Benchmarks show machine precision up to anisotropy order |α|=60, while the regularized functions also enter a lattice Euler–Maclaurin expansion.

  • Problem

    Existing anisotropic lattice-sum methods cover restricted exponents, low anisotropy orders, or vanishing shifts and wave vectors, while direct differentiation is numerically unstable at high orders.

    Source: §1 · p. 4

  • Method

    The paper derives stable anisotropic Epstein zeta representations using harmonic decomposition of the differential operator and backward-stable evaluation of harmonic polynomials.

    Source: §1 · p. 4

  • Results

    Benchmarks against independent reference values show machine precision up to |α|=60, with linear runtime growth for low- and mid-order anisotropy in one to three dimensions.

    Source: §6 · p. 44

  • Takeaways & Limitations

    Regularized anisotropic Epstein zeta functions provide coefficients for the singular Euler–Maclaurin expansion relating lattice sums to integrals for power-law-singular functions.

    Source: §1 · p. 3

    Near symmetry-related zeros, algorithm error grows with inverse distance to the zero and with anisotropy order.

    Source: §4.2 · p. 30

Abstract

from arXiv · show

The precise and efficient evaluation of large-scale lattice sums involving power-law kernels is a fundamental computational problem in the simulation of classical and quantum systems with long-range interactions. While methods for spatially isotropic kernels, some based on Epstein zeta functions, have advanced considerably in recent years, the anisotropic case has lagged behind, despite its broad relevance to both fundamental and effective interactions such as the dipole interaction in magnetic materials. In this work, we solve this issue by defining and analyzing anisotropic Epstein zeta functions for which we derive stably computable representations obtained from wave vector derivatives of lattice sums over isotropic interaction kernels. These functions find direct application in the analytical and numerical study of anisotropically interacting lattice systems. Going further, they provide the correction term in an exact equivalence between discrete lattices and their continuous analogs in a recent generalization of the classical Euler-Maclaurin summation formula to lattices and summands involving power-law kernels. Their connection to high-order derivatives of zeta functions can be used to improve convergence rates of numerical algorithms, for instance in micromagnetics, or to provide rapidly convergent expansions suitable for precomputations of generalized zeta functions. We derive a stably computable representation of anisotropic Epstein zeta functions, including the possibility for analytically removing Rayleigh--Wood singularities, and we develop a numerical algorithm for their stable evaluation for any lattice, power-law decay exponent and anisotropy order. We benchmark the algorithm against closed-form identities, direct summation and multi-precision results, obtaining machine precision across various lattices, power-law exponents, and anisotropy orders.

1. Introduction

Anisotropic power-law lattice sums matter in long-range interacting systems, but existing computational treatments have been limited and direct differentiation is unstable at high anisotropy. This work introduces anisotropic Epstein zeta functions and stable representations that support their computation and applications to lattice summation.

1. Introduction

  • Anisotropic Epstein zeta functions extend the isotropic Epstein zeta framework used in studies of magnetism, long-range Ising models, and theoretical chemistry.The paper connects these functions to shift and wave-vector derivatives, including derivatives in which the anisotropic monomial appears.Source: §1 · p. 3
  • Regularized anisotropic Epstein zeta functions supply coefficients in the singular Euler–Maclaurin expansion, representing the difference between a lattice sum and an integral through derivatives of the smooth factor.The expansion applies to integrable functions with a power-law singularity and band-limited smooth factor.Source: §1 · p. 3
  • The resulting representation enables efficient computation of singular sums at a numerical cost independent of the number of non-zero summands.The paper presents the singular Euler–Maclaurin expansion for a non-trivial three-dimensional geometry.Source: §1 · p. 4
  • Existing computational treatments covered only selected exponents, low-order anisotropy, or vanishing shift and wave vectors, while direct differentiation is unstable at high order.This left stable computation and a complete account of analytic properties inaccessible for general anisotropic kernels.Source: §1 · p. 4
  • The algorithm reaches strongly anisotropic regimes and obtains machine precision against independent references up to |α|=60, with linear runtime growth at low and mid orders in one to three dimensions.The representation uses a harmonic decomposition of the differential operator and a backward-stable harmonic-polynomial algorithm.Source: §1 · p. 4

2. Analytic properties

This section characterizes the analytic properties of anisotropic Epstein zeta functions, including their holomorphy, pole conditions, and symmetries, and connects them to wave-vector derivatives of set zeta functions. It develops Crandall-based representations and regularizations for continuation, singularity handling, and analytic extension.

2.1. Elementary properties.

  • For shifted lattices, anisotropic Epstein zeta functions agree with set zeta derivatives up to an oscillatory prefactor, while derivatives can be ill-defined at reciprocal-lattice points.The derivative correspondence has an exceptional case when ν = d + |α|, y lies in the reciprocal lattice, and all components of α are even.Source: §2.1 · p. 6
  • The anisotropic Epstein zeta function has inversion, translation, and scaling symmetries, and can vanish under specified inversion or mirror-symmetry conditions.The stated zeros include odd total derivative order with x in the lattice and 2y in its reciprocal lattice, as well as certain componentwise mirror conditions.Source: §2.1 · pp. 7–8

2.2. Derivative Crandall representation.

  • The derivative Crandall representation expresses the anisotropic zeta function through two superexponentially convergent series and gives its continuation, but floating-point precision degrades for large anisotropy order.The exact representation also permits singularity removal without cancellation error; the paper points to a harmonic-decomposition representation for stable evaluation at arbitrary anisotropy.Source: §2.2 · pp. 10, 12

2.3. Holomorphic extension.

  • The section establishes a holomorphic extension of the anisotropic Epstein zeta function while identifying possible singularities at ν = d + |α| and at lattice or reciprocal-lattice arguments.The extension is stated in Theorem 2.11; the supplied passage does not specify its full domain.Source: §2.3 · p. 12
  • The anisotropic zeta function extends holomorphically in its vector arguments on the stated lattice-complement domains; at lattice points, a simple pole at ν = d + |α| occurs only when χ(α) = 1.When both vectors lie in their respective lattices, it is holomorphic in ν if χ(α) = 0 and otherwise has that simple pole.Source: §2.3 · pp. 12–14

2.4. Singularity handling.

  • The original singularities near lattice shifts and reciprocal-lattice wave vectors arise from non-analyticity of the upper Crandall function at zero argument.Regularizing the power-law singularity separates it from a regular part, allowing integration without cancellation error.Source: §2.4 · p. 14
  • Replacing a diverging upper Crandall summand with the lower Crandall function yields a derivative Crandall representation for the regularized anisotropic zeta function.The regularization removes the wave-vector power-law singularity at the origin; the lower Crandall function is continuously extended there.Source: §2.4 · pp. 15–17
  • Explicit lower Crandall derivative representations extend the regularized theory to arbitrary anisotropy and complex vector arguments, generalizing the α = 0 result.The authors provide a new proof that extends to complex vector arguments.Source: §2.4 · p. 19
  • The regularized anisotropic zeta function extends holomorphically around the wave-vector origin, while the non-regularized derivative Crandall representation suffers catastrophic cancellation at high anisotropy.Derivatives of the regularized Crandall function can be stably computed for arbitrary anisotropy; the following section develops stable representations for the non-regularized case.Source: §2.4 · pp. 20–21

3. Stably computable representations

Direct floating-point derivatives of Crandall representations lose precision at high anisotropy; harmonic decomposition instead yields stable representations for anisotropic and regularized anisotropic Epstein zeta functions.

3. Stably computable representations

  • Directly computing upper Crandall derivatives in floating point introduces systematic precision loss for high-order anisotropy, motivating a different representation.The paper says stable representations cannot be obtained by applying partial derivatives directly to the derivative Crandall representation.Source: §3 · p. 21

3.1. Harmonic decomposition.

  • Harmonic decomposition splits partial derivatives into poly-Laplacians and homogeneous harmonic-polynomial factors that isolate anisotropy.The poly-Laplacians can be applied without differentiating, while the remaining operators have stably computable closed forms on upper Crandall functions.Source: §3 · p. 22, §3.1 · p. 22

3.2. Anisotropic Epstein zeta functions.

  • The resulting harmonic-decomposition Crandall representation provides stable evaluation of anisotropic Epstein zeta functions for arbitrary anisotropy.Its main result supports stable and efficient evaluation given a suitable implementation of the harmonic polynomials.Source: §3.2 · pp. 23, 25

3.3. Regularized anisotropic Epstein zeta functions.

  • The regularized anisotropic Epstein zeta function also has a Crandall-type representation based on harmonic decomposition, enabling efficient computation at high-order anisotropy.The representation is obtained by acting on the power-law singularity with the poly-Laplacian and the operator derived from harmonic polynomials.Source: §3.3 · pp. 25–26

4. Algorithm

The algorithm evaluates anisotropic Epstein zeta functions for arbitrary anisotropy using the harmonic-decomposition Crandall representation, with stable harmonic-polynomial evaluation as a key component.

4. Algorithm

  • The algorithm computes anisotropic Epstein zeta functions for arbitrary anisotropy from the harmonic-decomposition Crandall representation.The described implementation focuses on real parameters and sets the splitting parameter to λ = 1.Source: §4 · p. 27
  • The implementation addresses backward-stable harmonic-polynomial evaluation and cancellation errors that arise for large exponents.The coefficient sums can be computed exactly, and the algorithm also handles cancellation near shifted poles of the Gamma or upper Crandall function.Source: §4 · p. 27, §4.1 · p. 27

4.1. Harmonic polynomials.

  • Harmonic polynomials in the decomposition admit an explicit representation whose coefficients are dyadic rationals.The paper presents this as an efficient-computation representation for the polynomials.Source: §4.1 · p. 27

1. Rescale lattice matrix and vectors

The algorithm rescales and projects its lattice and vector inputs, then handles special zeros and sums in real and reciprocal space; several numerical cases require special care.

1. Rescale lattice matrix and vectors

  • The algorithm projects x and y to elementary lattice cells before its summation steps.It rescales the lattice matrix and vectors and forms projected coordinates from A^-1x and A^T y.Source: §1 · p. 28

4.2. Numerical stability.

  • For ν > 10 and |α| ≥ 2, nearby real-space summands can catastrophically cancel, so the algorithm extracts the singularity and changes summation order.Kahan summation is also used to reduce roundoff errors.Source: §4.2 · p. 29
  • Near symmetry-related zeros, error grows with inverse distance and anisotropy order; recentering shells may mitigate it, while symmetric z and −z addition is used near x = y = 0.The authors do not pursue recentering shells further.Source: §4.2 · p. 30

5. Numerical Experiments

The numerical experiments assess accuracy, stability, runtime, harmonic-polynomial evaluation, and application to a singular Euler–Maclaurin expansion. Reported errors reach machine precision, including for the order-six expansion, and the harmonic-decomposition algorithm remains stable at high anisotropy orders.

5.1. Machine precision for real parameters against independent reference values.

  • Independent-reference tests achieve E < 7 · 10^-14 for dimensions 1–3, and the all-equal-argument test across more than three million values has maximum error below 2 · 10^-13.The first bound is for a grid with anisotropy orders 2n from 2 to 12; the second uses second-order anisotropy across dimensions 1–6.Source: §5.1 · pp. 31–32, Table 1 · p. 32, Table 2 · p. 33
  • The eight-dimensional Ramanujan tau representation has maximum error E < 3 · 10^-14, while direct-summation benchmarks report E < 2 · 10^-15 for low-order and E < 3 · 10^-16 for high-order anisotropy.For generic arguments, comparison with a high-precision implementation gives typical errors E ≲ 10^-15 for |α| ≤ 10, E ≲ 10^-14 for |α| ≤ 20, and E ≲ 10^-12 for |α| ≤ 40.Source: §5.1 · pp. 33–34, Figure 2 · p. 34, Figure 3 · p. 34

5.2. Stability for high-order anisotropy.

  • The harmonic-decomposition Crandall algorithm stays flat with E < 2 · 10^-14 through anisotropy order 60, unlike the derivative representation whose error grows exponentially with order.Taylor reconstruction of the regularized Epstein zeta function converges exponentially, reaching Emax < 10^-15 for Taylor order n ≥ 54 in dimensions 1–3.Source: §5.2 · pp. 35, 37, Figure 5 · p. 36

5.3. Linear runtime for low-order anisotropy.

  • For low- to mid-order anisotropy, evaluation time increases roughly linearly with order; in dimensions 2–3, increasing the order from zero to ten increases runtime by about one order of magnitude.For orders up to 10, reported evaluation times are below 0.2 milliseconds in dimension 2 and below 0.9 milliseconds in dimension 3.Source: §5.3 · p. 37, Figure 7 · p. 38

5.4. Runtime and stability of harmonic polynomial evaluation.

  • Harmonic-polynomial evaluation has backward error below 2εM across anisotropy orders up to 30, while maximum evaluation time remains below 2 · 10^-7 after coefficient precomputation.The precomputation time grows by four orders of magnitude across orders 0–30, reaching tpre ≤ 2 · 10^-3 at α = (10, 20).Source: §5.4 · p. 40

5.5. Singular Euler-Maclaurin expansion.

  • The algorithm is applied to coefficients of the singular Euler–Maclaurin expansion for singular kernels sν(·) = |·|^-ν, and the reported correction-term error falls from E < 2 · 10^-2 at order ℓ = 0 to E < 2 · 10^-15 at ℓ = 6.The expansion is described as a generalization of classical Euler–Maclaurin summation, with applications to topological superconductors [13].Source: §5.5 · p. 40, §5.4 · p. 40, Figure 8 · p. 40
  • For tests with |α| = 24 across dimensions d = 1–4, the median condition numbers were below 29.The implementation precomputed and reused coefficients for the harmonic polynomials.Source: Table 3 · p. 41, §5.5 · p. 41
  • At α = (10, 20), precomputation took tpre ≤ 2·10^-3 s and maximum evaluation time was tmax eval < 2·10^-7 s; backward error stayed below 4.2·10^-16.These measurements cover the 2D high-performance harmonic-polynomial implementation based on Theorem 4.1.Source: Table 4 · p. 42, §5.5 · p. 42
  • The SEM analysis applies to singular Gaussian summands and uses anisotropic Epstein zeta representations for its correction terms.For ν ∉ d + 2N, the regularized anisotropic Epstein zeta function can be replaced by the non-regularized variant in the SEM operator representation.Source: §5.5 · p. 43
  • In the 3D fcc Gaussian test, the zero-order SEM error was less than 2%, while order six achieved machine precision, E < 2·10^-15.The test used σ = 10 and a truncated reference sum over z ∈ Afcc{−50, …, 50}^3.Source: §5.5 · p. 43

6. Outlook and Conclusions

The paper establishes stable representations for anisotropic Epstein zeta functions, including a harmonic-decomposition method for strongly anisotropic cases. These results support high-order computation and applications to lattice sums and precomputation.

6. Outlook and Conclusions

  • Two analytic representations provide meromorphic continuation, and the derivative Crandall representation establishes joint holomorphy with a conditional simple pole in ν.The pole occurs when y ∈Λ∗ and every component of α is even.Source: §6 · p. 44
  • Near y = 0, decomposition into a power-law singularity and a holomorphic part supports interpolation precomputation and avoids singularity-related cancellation errors.Source: §6 · p. 44
  • Machine-precision benchmarks reach anisotropy order |α| = 60, while runtime grows linearly for low- and mid-order anisotropy in one to three dimensions.The harmonic-decomposition representation relies on backward-stable evaluation of harmonic polynomials.Source: §6 · p. 44
  • The anisotropic Epstein zeta function gives the exact difference between lattice sums and integrals through the singular Euler-Maclaurin expansion, and supplies Taylor-series coefficients for precomputation and stable near-singularity calculations.Source: §6 · p. 44

Declarations

The declarations report funding acknowledgments and the authors’ use of generative AI for language editing and adversarial proof checking.

Declarations

  • J.B. acknowledges Quantum Fellowship Program of the German Aerospace Center funding, and A.B. acknowledges Klaus-Tschira Stiftung Grant No. 00.025.2025.Source: Funding · p. 44
  • The authors used Claude AI for grammar checking and minor language enhancements during manuscript preparation.Source: Funding · p. 44
  • Claude Opus 5 was also tasked adversarially with identifying errors and inconsistencies in proof drafts.Source: Funding · pp. 44–45

Appendix A. Absolute uniform convergence of set zeta derivatives

Appendix A establishes absolute and uniform convergence of differentiated set-zeta series under a decay condition on the exponent. Its proof reduces the tail to a convergent lattice sum using uniform discreteness.

Appendix A. Absolute uniform convergence of set zeta derivatives

  • For a uniformly discrete set L, the differentiated series is absolutely and uniformly convergent when Re(ν) > |α| + d.Source: Appendix A · p. 45
  • Uniform convergence of the differentiated series justifies iterated termwise differentiation.The proof applies the Weierstrass M-test to derivatives through order α.Source: Appendix A · p. 45
  • The tail is bounded by comparison with a scaled square-lattice sum, while the truncated region contains only finitely many points.Source: Appendix A · pp. 45–46

Appendix B. Derivatives of Crandall functions

Appendix B derives formulas for derivatives of Crandall functions by applying the chain rule to a squared argument. Homogeneity controls continuity at the origin, while the resulting derivative identities specify when meromorphic continuation has a simple pole.

Appendix B. Derivatives of Crandall functions

  • The Faà di Bruno formula organizes derivatives of a function composed with a squared argument as weighted sums of derivatives of the function.Source: Appendix B · p. 46
  • For a homogeneous polynomial h_m of degree m and Re(t) < m, h_m(u)G_t(u) extends continuously to u = 0, with a vanishing limit when m > 0.For m = 0, the limit is −2h_0(0)/t.Source: Appendix B · p. 46
  • Differentiating the upper Crandall representation yields a meromorphic continuation with a simple pole at ν = −|α| exactly when every component of α is even.Source: Appendix B · p. 47
  • The lower Crandall derivative identity extends holomorphically in its stated region and has a simple pole at ν = −|α| exactly when every component of α is even.Source: Appendix B · p. 48

Appendix C. Harmonic decomposition of the singularity

Appendix C develops the identities needed to apply harmonic-polynomial differential operators to singular power-law terms, including logarithmic cases. It uses these results to decompose the singularity and identify when algorithm branches vanish.

Appendix C. Harmonic decomposition of the singularity

  • Derivatives of the falling Pochhammer symbol simplify calculations involving derivatives of polynomials with logarithmic terms.The appendix derives an identity for these derivatives and uses it in subsequent harmonic-polynomial calculations.Source: Appendix C · pp. 48–49
  • For a homogeneous harmonic polynomial, its differential operator acting on a function of |y|^2 reduces to a multiple of the polynomial evaluated at y and a derivative of the function.The stated induction yields (a · ∇y)^n f(y2) = 2^n(a · y)^n f^(n)(y2) for the relevant harmonic components.Source: Appendix C · p. 50
  • The singularity decomposition treats generic exponents and integer cases, including terms with a logarithmic factor.For the logarithmic case, the proof distinguishes whether the integer exponent is below or above the combined operator order.Source: Appendix C · pp. 51–52
  • Algorithm branches can be skipped when the harmonic polynomial vanishes at zero, including when its degree is positive or χ(α) = 0.The lemma states these conditions imply hk,α(0) = 0.Source: Appendix C · p. 52

Appendix D. Derivation of the details on the algorithm

Appendix D justifies the harmonic-decomposition representation and details computations used by the algorithm. It establishes dyadic structure in its coefficients and convergence of a Taylor-series representation in anisotropic Epstein zeta functions.

Appendix D. Derivation of the details on the algorithm

  • The coefficients ck,α,γ are dyadic rationals, since each summand is an integer product divided by a power of two.The passage notes that the Pochhammer symbols, binomial and multinomial coefficients, and q coefficients involved are integers.Source: Appendix D · p. 53
  • Reordering the terms in the harmonic-decomposition Crandall representation is justified for the representation used to stabilize computations at large exponents.The proof begins by restricting the sums to k outside ν/2 + N, then establishes the stated identity.Source: Appendix D · pp. 53–54
  • The Taylor-series representation can be written in terms of anisotropic Epstein zeta functions and converges pointwise at the specified y.The proof establishes normal convergence in a polydisc within the region of holomorphy.Source: Appendix D · p. 55
  • The algorithm's early zero return is equivalent to the lattice conditions given in Lemma 2.6.The proof relates the reciprocal-lattice condition to the existence of a lattice vector that cancels the relevant component.Source: Appendix D · p. 55

Appendix E. Truncation

Appendix E derives a truncation remainder bound for the harmonic-decomposition Crandall representation and applies it to cutoff selection. It relates the bound for general lattices to their condition number and gives a square-lattice error target.

Appendix E. Truncation

  • The remainder estimate generalizes the Crandall-formula remainder estimate for the Epstein zeta function.The estimate applies to truncation of the harmonic-decomposition Crandall representation for regular lattices Λ = AZ^d with det(A) = 1.Source: Appendix E · pp. 58–59
  • After rescaling the truncation radius, the remainder bound depends on the lattice through a single global power of its condition number.The bound concerns general lattices Λ = AZ^d.Source: Appendix E · p. 60
  • Table 5 gives square-lattice truncation values targeting absolute error below 10^-18 for |ν| ≤ 10 and all |α| ≤ n.The bound uses ε = 1/50 and lists cutoff values as a function of dimension and anisotropy order.Source: Table 5 · p. 61, Appendix E · p. 61
  • For the tabulated cutoffs, the remainder is below machine precision for general lattices when κ(A)^(d+n+1) ≤ 10^2.The condition is stated for the error bound obtained using the Table 5 values r0.Source: Appendix E · p. 61

Appendix F. Identities for numerical experiments

Appendix F supplies identities used in the numerical experiments, including a relation involving a Ramanujan tau L-function and an identity for a three-dimensional singular Gaussian integral. It also establishes an equivalence for vanishing wave vectors and a holomorphic continuation for the Gaussian-integral identity.

Appendix F. Identities for numerical experiments

  • For vanishing wave vectors, the anisotropic Epstein zeta function corresponds to its regularized counterpart.The proof compares their derivative Crandall representations and treats the additional reciprocal-space summand.Source: Appendix F · p. 61
  • A theta-Mellin derivation yields the stated identity involving an anisotropic Epstein zeta value and the Ramanujan tau L-function.The argument first establishes the identity in a restricted half-plane and then uses holomorphicity to extend the equality.Source: Appendix F · pp. 62–63
  • The three-dimensional singular Gaussian integral identity extends holomorphically from 0 < Re(ν) < 3 to ν ∈ C \ (3 + 2N).The passage establishes the right-hand side as the holomorphic continuation of the integral on that domain.Source: Appendix F · pp. 63–64
Loading 2609.28282v1…