Epstein Zeta Library

EpsteinLib

Authors: Andreas A. Buchheit, Jonathan K. Busse, Ruben Gutendorf, DevOps: Jan Schmitz

Contact: buchh.nosp@m.eit@.nosp@m.num.u.nosp@m.ni-s.nosp@m.b.de

EpsteinLib is a C library designed for the fast and efficient computation of the Epstein zeta function for arbitrary multidimensional lattices. In addition to the C library, we also offer a Python package, epsteinlib, which can be easily installed via pip. For more information on the properties of the Epstein zeta function and on the underlying algorithm, see our recent manuscript 10.1093/imanum/drag057.

Originally studied by Epstein [1,2], the Epstein zeta function forms the basis for computing general multidimensional lattice sums in classical and quantum physics applications [3]. Together with its regularization, it serves as the central ingredient in the singular Euler-Maclaurin (SEM) expansion, which generalizes the 300-year-old Euler summation formula to lattice sums in higher dimensions with physically relevant power-law interactions [4-5]. An efficiently computable representation of the Epstein zeta function is provided in [6,7,8]. In [8], we discuss in detail the analytical properties of the Epstein zeta function and present an algorithm for its computation, complete with error bounds.

The anisotropic Epstein zeta function that appears in vector derivatives of Epstein zeta functions is also included in this library.

Epstein zeta function

For a \(d\)-dimensional lattice \(\Lambda=A\mathbb Z^d\), with \(A\in \mathbb R^{d\times d}\) regular, \(\boldsymbol x,\boldsymbol y \in \mathbb R^d\), and \(\nu \in \mathbb C\), the Epstein zeta function is defined by the Dirichlet series

$$ Z_{\Lambda,\nu}(\boldsymbol x,\boldsymbol y) = \sum_{z \in \Lambda}{}^{'} \frac{e^{-2\pi i \boldsymbol y \cdot \boldsymbol z}}{\left| \boldsymbol x- \boldsymbol z\right|^\nu},\quad \mathrm{Re}(\nu)>d, $$

which can be meromorphically continued to \(\nu \in \mathbb C\). Here, the primed sum excludes the case \(\boldsymbol z = \boldsymbol x.\)

The Epstein zeta function is implemented in this library as

double complex epsteinZeta(double nu, unsigned int dim, const double *A, const double *x, const double *y);

In the Python package, it is implemented as

def epstein_zeta(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
) -> complex

In the Julia package, it is implemented as

epsteinzeta(ν::Float64, A::Matrix{Float64}, x::Vector{Float64}, y::Vector{Float64})::Complex{Float64}

and with optional keyword arguments as

epsteinzeta(ν; d, A, x, y)

where at least one of the arguments d, x, y, or A must be provided. By default, x and y are zero vectors of length d, and A is the d × d identity matrix.

In the Mathematica package, it is implemented as

EpsteinZeta[\[Nu],A,x,y]

and evaluates to full precision over the whole parameter range up to ten dimensions.

The Epstein zeta function admits singularities in the lattice \(\boldsymbol x\in\Lambda\) and in the reciprocal lattice \(\boldsymbol y\in\Lambda^*\). To ensure numerical stability when evaluating the Epstein zeta function, we implement the following cutoffs:

• If \((\boldsymbol x-\boldsymbol z)^2 < 10^{-64}\) for some \(\boldsymbol z \in \Lambda\), we numerically set \(\boldsymbol x = \boldsymbol z\).
• Similarly, if \((\boldsymbol y-\boldsymbol k)^2 < 10^{-64}\), for some \(\boldsymbol k \in\Lambda^*\), we numerically set \(\boldsymbol y = \boldsymbol k\).

When evaluating \(\boldsymbol x\in\Lambda\) or \(\boldsymbol y \in\Lambda^*\), users should set \(\boldsymbol x= \boldsymbol 0\) or \(\boldsymbol y = \boldsymbol 0\) and use the quasi-periodicity of the Epstein zeta function.

In addition, this library includes the regularized Epstein zeta function, which is analytic around \(\boldsymbol y= \boldsymbol 0\), and is defined via

$$ Z_{\Lambda,\nu}^{\mathrm{reg}}(\boldsymbol x,\boldsymbol y) = e^{2\pi i \boldsymbol x\cdot\boldsymbol y} Z_{\Lambda,\nu}(\boldsymbol x,\boldsymbol y ) -\frac{\hat{s}_{\nu}(\boldsymbol y)}{V_{\Lambda}}, $$

where \(V_{\Lambda}=|\det A|\) is the volume of the elementary lattice cell, and

$$ \hat{s}_\nu(\boldsymbol y) = \frac{\pi^{\nu/2}}{\Gamma(\nu/2)}\Gamma\big((d-\nu)/2\big) (\pi \boldsymbol y^2)^{(\nu - d)/2},\quad \nu \not\in (d+2\mathbb N_0) $$

is the distributional Fourier transform of \(\vert\boldsymbol z \vert^{-\nu}\), where \(\Gamma\) denotes the gamma function and we adopt the choice

$$ \hat s_{d+2k}(\boldsymbol y)= \frac{\pi^{k+d/2}}{\Gamma(k+d/2)}\frac{(-1)^{k+1}}{k!} ( \pi \boldsymbol y^2 )^{k} \log (\pi \boldsymbol y^{2}),\quad k\in \mathbb N_0. $$

In the c library, the regularized Epstein zeta function is included as

double complex epsteinZetaReg(double nu, unsigned int dim, const double *A, const double *x, const double *y);

in the Python package as

def epstein_zeta_reg(
    nu: Union[float, int],
    A: NDArray[Union[np.integer[Any], np.floating[Any]]],
    x: NDArray[Union[np.integer[Any], np.floating[Any]]],
    y: NDArray[Union[np.integer[Any], np.floating[Any]]],
) -> complex

in the Julia package as

epsteinzetareg(ν::Float64, A::Matrix{Float64}, x::Vector{Float64}, y::Vector{Float64})::Complex{Float64}

and with optional keyword arguments as

epsteinzetareg(ν; d, A, x, y)

with defaults for x, y, and A identical to those used in epsteinzeta, and in the Mathematica package as

EpsteinZetaReg[\[Nu],A,x,y]

To ensure numerical stability when evaluating the regularized Epstein zeta function as a function of \(\boldsymbol x\), we again implement the following cutoff:

• If \((\boldsymbol x-\boldsymbol z)^2 < 10^{-64}\) for some \(\boldsymbol z \in \Lambda\), we numerically set \(\boldsymbol x = \boldsymbol z\).

Anisotropic Epstein zeta function

Let \(\nu\in\mathbb C\) and signify by the multi-index \(\boldsymbol\alpha\in\mathbb N_0^d\) the anisotropy strength of

$$ V_{\nu,\boldsymbol \alpha}(\boldsymbol z) = \frac{\boldsymbol z^{\boldsymbol \alpha}}{\vert \boldsymbol z \vert^\nu} ,\qquad \boldsymbol z\in\mathbb R^d\setminus\{\boldsymbol 0