On the Efficient Computation of Large Scale Singular Sums with Applications to Long-Range Forces in Crystal Lattices

We develop a new expansion for representing singular sums in terms of integrals and vice versa. This method provides a powerful tool for the efficient computation of large singular sums that appear in long-range interacting systems in condensed matter and quantum physics. It also offers a generalised trapezoidal rule for the precise computation of singular integrals. In both cases, the difference between sum and integral is approximated by derivatives of the non-singular factor of the summand function, where the coefficients in turn depend on the singularity. We show that for a physically meaningful set of functions, the error decays exponentially with the expansion order. For a fixed expansion order, the error decays algebraically both with the grid size, if the method is used for quadrature, or the characteristic length scale of the summand function in case the sum over a fixed grid is approximated by an integral. In absence of a singularity, the method reduces to the Euler–Maclaurin summation formula. We demonstrate the numerical performance of our new expansion by applying it to the computation of the full nonlinear long-range forces inside a domain wall in a macroscopic one-dimensional crystal with 2×1010\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\times 10^{10}$$\end{document} particles. The code of our implementation in Mathematica is provided online. For particles that interact via the Coulomb repulsion, we demonstrate that finite size effects remain relevant even in the thermodynamic limit of macroscopic particle numbers. Our results show that widely-used continuum limits in condensed matter physics are not applicable for quantitative predictions in this case.


Introduction
Large sums appear everywhere in nature; our macroscopic world is composed of microscopic particles whose interaction forces determine the properties of the world we live in. Sums with singularities describe discrete long-range interacting systems in condensed matter and quantum physics [1], with examples ranging from the computation of forces and energies in atomic crystals [2] to the study of charge transfer in DNA strings [3]. Making predictions about these sums is however in general a difficult task. In many cases, the evaluation of an integral is easier than the evaluation of a sum, either because there are more tools available on the analytical side or because on the numerical side, efficient quadrature schemes are known. The question arises how sums and integrals are related and how we can approximate one by the other.
A part of the answer to this question was given independently by Leonard Euler in 1736 and by Colin Maclaurin in 1742 [4]. Consider Fig. 1, which provides an illustration for the approximation of a sum (rectangles) by an integral (blue region). In the red parts, the sum dominates the integral, whereas in the green parts, the integral is larger than the sum. The Euler-Maclaurin (EM) expansion describes this difference between sum and integral of a sufficiently differentiable function in terms of derivatives evaluated at the limits of integration plus a remainder integral.
Before we state the EM expansion, we introduce the following standard notation: The set of integers is denoted by Z,  defined analogously, where only the existence of the corresponding one-sided limit is needed at the end points. Regulated functions are continuous up to a countable number of points.
For a, b ∈ Z, a < b, δ ∈ (0, 1], ∈ N, and for a function f ∈ C +1 [a + δ, b + δ], the EM expansion reads [4,5] where y is the smallest integer larger than or equal to y. The functions B are the Bernoulli polynomials, which are uniquely defined by the recurrence relation For a derivation of the EM expansion as well as a brief introduction to its history we refer to [4]. The applicability of the EM expansion for the approximation of a particular sum is determined by the scaling of the remainder integral on the right hand side of (1.1) with the expansion order . The remainder integral cannot be evaluated in a numerically feasible way, as the integrand is a piecewise defined function. In the numerical application, an order is chosen, the remainder integral is discarded and the error made by discarding the remainder integral is estimated. If the addend is based on an entire function whose derivatives in addition satisfy certain bounds (more details in the next section), the remainder integral decreases exponentially with . This is the ideal case for the EM expansion. For most functions however, even for smooth functions, the error begins to diverge after a certain threshold value of is reached, thus limiting the precision that the EM expansion can offer. There have been a number of valuable extensions of the classic works of Euler and Maclaurin, which make the expansion applicable to a larger set of functions. Important progress has been made by Navot [6], where the EM expansion is generalised to functions with an algebraic singularity at the boundaries of the integration interval. Further generalisations have been developed by Monegato and Lyness [5] and Sidi [7,8], where divergent integrals are regularised by the use of Hadamard finite part integrals. Furthermore, a higher dimensional generalisation of the EM for simple lattice polytopes has been found [9]. We also mention here a recent alternative approach to the EM expansion where the difference between sum and integral is written in terms of integrals only [10].
One particular set of functions, for which the EM expansion fails to converge and which are extremely important in practice are functions that involve an asymptotically smooth singularity, see Definition 2.1. Unfortunately, these functions are of strong practical interest as all physical interactions belong to this kind. In this work, we present the singular Euler-Maclaurin expansion (SEM), which makes the classic expansion applicable to functions that involve an asymptotically smooth singularity. This paper targets a wide range of different audiences. Therefore it is organised as follows. Section 2 concludes the main results regarding the SEM, offering all the tools needed for application. Section 3 covers details to apply the SEM as a numerical tool. We discuss its numerical performance and apply it to a macroscopic long-range interacting crystal as a physically relevant example. Note that a basic implementation of the SEM in Mathematica is provided along with the article. 1 Section 4 provides the main derivation of the SEM. It includes the proofs of the most important theorems and propositions. Finally, in Sect. 5, we make our concluding remarks.

Main Result and Notation
In this chapter, we formulate the SEM expansion for intervals [a + δ, b + δ], with a, b ∈ Z, δ ∈ (0, 1], which is in particular applicable, if the function splits into two factors where x ∈ Z, s ∈ C ∞ (R\{0}) has a singularity at 0 and g : [a + δ, b + δ] → C is a sufficiently differentiable function. We apply the following general strategy: singularities at x or other smooth functions whose derivatives increase quickly with the derivative order are included in s. The function s limits the applicability of the standard EM expansion to f and therefore requires a special treatment. In practice, the function s often represents a pairwise long-ranged interaction potential or the one-dimensional forces generated by such a potential. We refer to s in the following as the interaction. The remaining factor g includes well-behaved functions, whose derivatives increase sufficiently slowly with the derivative order. The slower the derivatives increase, the better are the convergence rates. We briefly outline our presentation of the SEM expansion. We first discuss the properties of the function s that are required for the expansion. The interaction is subsequently made integrable by introducing an exponential regularisation. We then make use of the integrability of the regularised interaction and define the Bernoulli-A functions, a generalisation of the periodised Bernoulli polynomials in which we encode s. These functions then form the coefficients of the differential operator of the SEM, replacing the Bernoulli polynomials in the standard EM expansion (1.1). The differential operator acts on g only, avoiding the fast increase in the derivatives of s that causes the breakdown of the standard EM expansion. A finite-order approximation of this differential operator leads to the finite-order SEM expansion. For smooth g, whose derivatives fulfil certain bounds, we take the order of the expansion to infinity, leading to the infinite order SEM.
Before moving on to the formulation of the SEM expansion, we need to specify the admissible set of functions for the interaction: the function s has to be asymptotically smooth [11,Sec. 3.2]. for all y ∈ R\{0} and ∈ N. We denote the vector space of all asymptotically smooth functions by S. It is the union of the linear subspaces S α , α ∈ R, where s ∈ S α if and only if for all ξ > 0 there is c ξ > 0 such that This set includes entire functions like polynomials, but also a broad set of functions with singularities. Typical examples for asymptotically smooth functions are for ν ∈ R which are singular for ν > 0. From 2.1, we find that the th derivative of the interaction may scale with the factorial of . It is this fast increase in the derivatives that causes the breakdown of the EM expansion, and therefore, taking derivatives of s has to be avoided. It turns out that we can integrate s instead. However, s is in general not integrable on [1, ∞); take for instance ν = 1 in (2.3). This challenge is overcome by using an exponentially decaying regularisation of the interaction.

Notation 1 Let s ∈ S. The exponentially weighted interaction reads
For β 0, the weighting is gradually removed and the interaction regains its original range. It is crucial here that the reduction of the interaction range is introduced not as a sharp cut-off, but in a smooth way, such that s β remains asymptotically smooth.
We now use regularised differences between sums and integrals in order to define a replacement for the Bernoulli polynomials (1.1). We call this replacement the Bernoulli-A functions, in which we encode all information about the interaction s. In the following, calligraphic symbols indicate an explicit dependence on the interaction.

Definition 2.2 (Bernoulli-A functions) Suppose s ∈ S. The Bernoulli-A functions are defined as
The shifted monomials in sum and integral guarantee that the Bernoulli-A functions are subsequent antiderivatives, A +1 = A . A special role is assigned to the function A 0 that exhibits jump discontinuities at the positive integers. Both the size of the jump discontinuities and the left sided derivative of A 0 are equal to value of the interaction, The properties of the Bernoulli-A functions are discussed in detail in Proposition 4.2.
A different recursive definition of generalised Bernoulli functions on bounded intervals [a, b] starting from a bounded function by iterated integration on [a, b] is presented in [12]. Due to the chosen normalisation, requesting that the integrals of these functions over [a, b] to vanish, the generalisation to nonintegrable singularities and unbounded domains is not possible. Our definition is more general and, in particular, allows for arbitrary algebraic singularities at the origin.
The Bernoulli-A functions replace the periodic extension of the Bernoulli polynomials in the differential operator part and the remainder of the EM expansion. We display them for s(y) = |y| −1 in Fig. 2. For a constant interaction s = 1, we recover a rescaled version of the periodised Bernoulli polynomials, a proof of which is given in "Appendix 1". We now define the SEM differential operator as follows:

Definition 2.3 (SEM operator)
For s ∈ S and y ∈ R + , we define the differential operator of infinite order where D is the derivative operator. We call D y the SEM operator. For ∈ N the finite order approximations D ( ) y are given by Finally, we move on to the formulation of the new expansion.
The SEM operator only acts on g, not on f , and therefore differentiation of s is avoided. In order to perform the limit → ∞ in Theorem 2.1, the function g has to belong to the set of functions of exponential type. This restricts the growth of its derivatives.

Definition 2.4 (Band-limited functions)
Let g ∈ C ∞ (R). We say that g is band-limited with bandwidth σ and write g ∈ E σ , if for every bounded interval I ⊂ R, there exists M > 0 such that The admissible range for σ depends on s, in particular on γ from Definition 2.1.
Theorem 2.2 (SEM error scaling) Assume the conditions of Theorem 2.1 and let g ∈ E σ , σ > 0. Then the remainder in the SEM expansion, for τ = 2π/(γ + 1) and a polynomial P whose degree only depends on s.
For σ < τ, the approximation error of a sum by means of the SEM expansion falls off exponentially in the expansion order and polynomially in the parameter σ . In particular, the limit → ∞ is well-defined. Hence our new expansion succeeds in providing a reliable approximation where the standard EM expansion, which does not converge for singular functions, fails.

Remark 2.1
For simplicity, we have formulated Theorem 2.1 such that x is positioned to the left of the interval [a + 1, b]. Of course, the theorem can also be applied in case that x is positioned to the right of this interval by simply reflecting the interval about x. Consider and thus x ≤ã <b. Using the reflected interval, we can transform the sum as follows with the functionss ∈ S andg such that Thus the SEM becomes applicable to the right hand side of (2.6).

Model Description
We demonstrate the numerical performance of the SEM expansion by applying it to the calculation of long-range forces in a macroscopic one-dimensional crystal lattice of charged particles. The SEM expansion naturally provides the answer to the question how to correctly include the discreteness of the lattice within a continuum formulation that avoids discrete sums and is therefore efficiently computable for all asymptotically smooth interaction forces.
We consider a particularly difficult scenario, the case where the interaction potential is the Coulomb repulsion, which decays algebraically with an exponent equal to the system dimension. Then the discreteness of the crystal has an observable effect on the forces at all scales, which makes a continuum approximation challenging [2]. We consider a one-dimensional crystal of 2N +1 particles, N ∈ N, and denote the particle positions as x j ∈ R, j = −N , . . . , N . The particles are displaced from an equidistant grid with lattice constant h > 0, through a smooth displacement function u.
If the interaction energy V ∈ S between two particles decays algebraically with their distance x, the force acting on the particle with reference position x reads where the function f factors into f (y) = s(y − x)g(y) with s ∈ S and g smooth such that . (3.5) We now have determined a suitable factorisation for the summand function. In order to subsequently apply the SEM expansion to above sum, we only need to determine the corresponding SEM operator, which amounts to computing its coefficient functions A .

Bernoulli-A Functions for Physical Interactions
We move on to the computation of the Bernoulli-A functions for the algebraic interaction | · | −ν , which is of high relevance in physical systems. In case of ν = 1, it either describes the Coulomb interaction between charged or the gravitational interaction between massive particles. For ν = 3, it describes the dipole interaction between spins. The Bernoulli-A functions determine the SEM differential operator and hence provide the tool needed for the efficient computation of large singular sums. It turns out that the coefficient functions can be compactly written in terms of the Hurwitz zeta function. A proof is given in Appendix A.
Theorem 3.1 For s = | · | −ν and ν ∈ R\N + , the Bernoulli-A functions read with ζ the Hurwitz zeta function, For ν ∈ N + , the coefficients are well-defined, where γ e is the Euler-Mascheroni constant and H k denotes the kth harmonic number, After having determined analytic formulas for the Bernoulli-A functions for a large set of physically relevant interaction functions, we proceed by demonstrating that the SEM expansion is able to describe long-range interactings systems in a very precise way, where, in contrast, standard integral approximations fail.

Numerical Results
For our numerical study, we investigate a chain of particles interacting via the Coulomb repulsion. The displacement function is chosen as the integral of a normalised Lorentzian which describes the simplified profile of a kink, a domain wall in the crystal, with the asymptotic behaviour lim y→∞ u(y) = 1 and lim y→−∞ u(y) = 0. Kinks typically arise when an additional nonlinear potential is applied to the crystal. The parameter λ > 0 controls the width of the kink. When adjusting the Definition 2.4 of band-limited functions such that only up to + 1 derivatives need to fulfil the estimate then the error estimates for the SEM expansion in Theorem 2.2 continue to hold. The function u is thus effectively band-limited with a band-width σ proportional to 1/λ. For further details regarding kinks in condensed matter physics, see [13]. A schematic depiction of a kink is displayed in Fig. 3. In the following, we compute the Coulomb force that is exerted on an arbitrary particle (red circle) by all other particles (black circles). We compute the Coulomb forces in (3.4) by exact summation and subsequently compare the result to the integral approximation, which is standard in condensed matter physics, choosing the integral offset δ = 1. For a mesoscopic chain with N = 10 3 and a kink width λ = 10, the resulting forces are displayed in the left panel of Fig. 4, where the blue dots show the exact results and the black curve shows the integral approximation. The integral approximation in the mesoscopic system reproduces the correct qualitative behaviour, showing that the particles are drawn towards the kink, as it constitutes a delocalised particle hole. However, it strongly underestimates the absolute value of the forces and thus cannot be used for reliable quantitative predictions. We can now improve significantly upon the precision of the integral approximation by using the SEM expansion of order = 1. The result for the SEM expansion is shown in red, which is visually indistinguishable from the result obtained by exact summation.The SEM reproduces the exact forces very precisely both at the chain edges as well as inside the kink in the centre with an absolute error smaller than 3 × 10 −7 for all particles with 4 digits of precision. Moreover, the runtime for the adaptive numerical integration and for the computation of the derivatives is essentially independent of N and λ for a single force evaluation. The SEM expansion thus combines the precision of exact summation with the O(1) scaling of the integral approximation. 2 In the next step, we can now apply the SEM expansion to the computation of the full nonlinear long-range Coulomb forces in a crystal of macroscopic size, where the computation of the exact forces has become infeasible. We choose λ = 10 5 and N = 10 10 , which for a typical lattice constant h ≈ 10 −10 m, corresponds to a crystal with a total length of two metres. The first order SEM (red line) is compared to the integral approximation (black line). Even at the macro scale, the SEM shows a visible difference to the integral approximation. This difference can be highly relevant if we want to compute material properties ab-inito. The finite-size effects in the crystal remain relevant even in for the macroscopic system. We have therefore demonstated that the often made claim, that a lattice sum may be replaced by an integral if the underlying charge distribution is sufficiently broad, see e.g. the discussion around Eq. (5.3) in [14], is not correct in case of ν = 1 in one dimension.
We now analyse the scaling of the absolute error in the maximum norm for different SEM orders and kink widths λ for N = 200. The results are displayed in Fig. 5. The smaller particle number is chosen, such that the exact forces can still be computed efficiently for all particles. We find that in case of the integral approximation (black dots), the maximum absolute error does not scale with λ. The maximum error occurs at the chain edges, the error in the centre scales as λ −2 . Note that the inclusion of the zero order SEM contribution already compensates the error at the edges, with the maximum error now appearing close to the kink centre, in the region where the derivatives of u are large. For = 1 (red dots), the first order SEM, the error scales approximately as λ −4 . For odd orders we find that the error scaling coefficient is approximately + 3. The exact scaling coefficients calculated from a linear fit of the last 5 data point is given in Fig. 5. For = 7 (purple dots) and λ = 25, the SEM offers an absolute error smaller than 10 −17 which corresponds to at least 13 digits of precision for all forces. In conclusion, we find that the SEM expansion yields an error that falls off exponentially in the expansion order and polynomially in the width λ, reproducing the theoretically predicted error scaling.

Derivation of the Singular Euler-Maclaurin Expansion
We begin the derivation of the new expansion by showing the well-definedness of the Bernoulli-A functions, on which it is based.

Proposition 4.1 For ∈ N, the function A is well-defined and obeys the estimate
with τ = 2π/(1 + γ ) and P a polynomial that depends only on s and ξ .
For the proof of above proposition, we need two lemmata that provide estimates on the derivatives of functions that involve the interaction function s.

Lemma 4.1 Let s ∈ S with constants c
for all z > 0, ∈ N and β > 0.
Proof For ∈ N and y > 0, we compute The behaviour of the derivatives of the integrand in Definition 2.2 of the Bernoulli-A functions is investigated in the following.
Proof For y > 0 we compute In case of z = y, only the term with j = contributes and we can insert the estimate from 4.1. For z > y and k ≥ , we can set the upper limit of the sum to and find with C = c c ξ with ξ = y and the lemma follows after noticing that |y/z − 1|γ < γ and that k!/ ! ≤ k k− for k ≥ .
With above lemma, we now prove the well-definedness of the Bernoulli-A functions and provide a bound.
Proof (Proposition 4.1) For β > 0, we first introduce the auxiliary function A ,β (y) = ∞ n= y s y, ,β (n) − ∞ y s y, ,β (z) dz, y > 0, and subsequently prove the existence of the limit β 0 using the standard EM expansion, which takes the simplest form for the integral offset δ = 1 + y − y , namely with the remainder y, ,β (z) dz.
Clearly, the limit β 0 is well defined for the sum over k. Moreover, setting the expansion order m = + max{0, α + 1 }, garantees, by Lemma 4.2, that a β-independent and integrable majorant for the integrand in the remainder can be established. Hence, the limit and thus A is well-defined by the dominated convergence theorem.
We now use above representation in order to provide a bound to A that falls off exponentially in . To this end, we first recall the following bound of the Bernoulli polynomials, see e.g. Eq. (19) and following discussion in [15], Inserting this bound as well as the estimate from Lemma 4.2 into (4.1), we find with m α = max{0, α + 1} and C = c c ξ . Thus, with τ = 2π/(1 + γ ), and a constant C > 0 and a polynomial P that depend only on ξ > 0 and s.
Using above estimates, we are now in the position to prove the properties of the Bernoulli-A functions.

Proposition 4.2 (Properties of
Then A is smooth on R + \N + and A ∈ C −1 (R + \N + ). A +1 is an antiderivative of A . Finally, the function A 0 obeys the jump relation lim ε 0 A 0 (n) − A 0 (n + ε) = s(n), n ∈ N + , and the derivative relation Proof The binomial theorem yields the following representation of the Bernoulli-A functions, where we have used that the limit in β on the right hand side exists due to asymptotic smoothness of y k s β (y). Using that the β-limit is well-defined, we already find that the βlimit is uniform in y on any compact subset of R + by choosing m ∈ N + with m > y and writing with I β the y-independent remaining regularised difference between sum and integral. We now show that A ∈ C 0 (R + ) for all ≥ 1. Above representation implies smoothness of A on R + \N + , hence we only need to investigate y ∈ N + . As · is left continuous, it is sufficient to study the finite difference from the right, finding which vanishes for ≥ 1. For = 0, we obtain the jump condition We now consider differentiability of A . For y ∈ R + \N + , the function A is smooth as the composition of smooth functions. With the same argument, it is left differentiable on R + . We now compute the left-handed difference quotient of A . We find where δ 0, denotes the Kronecker delta. For ≥ 1, the first term on the right hand side equals A −1 and the second term vanishes, showing that A is an antiderivative of A −1 . Thus, . For = 0 on the other hand, the first term vanishes and the second term yields the derivative relation.
which extends to the undirected limit ε → 0 on y ∈ R + \N + as A 0 is smooth outside of integers.
Using Proposition 4.2, we can now formulate the zero order SEM expansion.
Subsequently (4.3) is divided into two separate sums. On the right hand side, an index shift is performed in the sum that includes the terms A 0 (n − x − ε) resulting in the expression (4.4) where the second term results from an adjustment of the differing summation intervals. Going back to the integral on the left hand side of (4.2), we use the derivative property in Proposition 4.2 and find b+δ a+δ

A recombination of the two sums yields
Integration by parts on all integrals in (4.5) in order to remove the derivatives of A 0 from the expression yields (4.5) where the integrals have been combined to a single one by taking the limit ε → 0. Subtracting (4.5) from (4.4), we obtain As A 0 is left continuous ( · is left continuous), the proposition follows after performing the limit ε → 0.
We proceed with the proof of the main theorem.
The case = 0 is proved in Proposition 4.3. The case of higher expansion orders ≥ 1 readily follows via iterated integration by parts, where successive antiderivatives of A 0 are given by (A k ) k∈N , see Proposition 4.2.
Finally, the bound on the remainder formulated in Theorem 2.2 is a direct consequence of the bound on A from Proposition 4.1 and the defining property of the band-limited function g, see Definition 2.4.

Outlook
In this work, we have developed the singular Euler-Maclaurin (SEM) expansion that, first, allows for the precise and efficient evaluation of large sums with singularities and, second, provides precise quadrature rules for singular integrals. With the new method, we can quantify finite-size effects in condensed matter systems with nonlinear long-range interactions, improving upon standard approximations in the field. Using the forces in a domain wall in a crystal of charged particles as a physically relevant example, we have shown that the method provides an approximation to the force sum in a runtime that does not depend on the particle number and with an approximation error that falls of exponentially with the expansion order and algebraically in the characteristic lengthscale of the function. Hence we are able to provide reliable predictions about the properties of solids even for macroscopic particle numbers. As a first important result, we have shown that the often-used integral approximation tends to underestimate the absolute value of forces inside a domain wall if the particles interact via the Coulomb repulsion. We expect that the SEM expansion will serve as an important tool in condensed matter and quantum physics, helping us to replace qualitative estimates by reliable predictions. It would for instance be interesting to determine precise time evolutions and stationary states through the solution of the integro-differential equations that arise from the application of the SEM.
Acknowledgements We would first like to thank Daniel Seibel and Darya Apushkinskaya for inspiring discussions. Our colleagues Peter Schuhmacher and Daniel Seibel have taken the time to proof-read earlier version of this manuscript and have offered valuable suggestions which improved this work; we acknowledge your support. We would like to express our gratitude towards our supervisor Prof. Sergej Rjasanow, whose support and guidance made this work possible. Finally, we are grateful for the suggestions and comments of the anonymous referee, which improved this work.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data Availability Statement
The code needed to reproduce the results of this publication is available in our github repository https://github.com/andreasbuchheit/singular_euler_maclaurin

Conflict of interest
The authors declare that they have no conflict of interest.

Availability of Data and Materials Available on github.
Code Available on github.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Inserting this in Definition 2.2 shows
for ∈ N and ν ∈ R\N. For ν ∈ N + , above expression has a removable singularity. Given k ∈ N, we write By Eq. (9)