Singular Euler-Maclaurin expansion

We present the singular Euler--Maclaurin expansion, a new method for the efficient computation of large singular sums that appear in long-range interacting systems in condensed matter and quantum physics. In contrast to the traditional Euler--Maclaurin summation formula, the new method is applicable also to the product of a differentiable function and a singularity. For suitable non-singular functions, we show that the approximation error decays exponentially in the expansion order and polynomially in the characteristic length scale of the non-singular function, where precise error estimates are provided. The sum is approximated by an integral plus a differential operator acting on the non-singular function factor only. The singularity furthermore is included in a generalisation of the Bernoulli polynomials that form the coefficients of the differential operator. We demonstrate the numerical performance of the singular Euler--Maclaurin expansion by applying it to the computation of the full non-linear long-range forces inside a macroscopic one-dimensional crystal with $2\times 10^{10}$ particles. A reference implementation in Mathematica is provided online.


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, N = {0, 1, 2, . . . } are the nonnegative integers, and N + = N \ {0} are the positive integers. Similarly, R are the real numbers, denote the one-sided limits from left and right. The vector space C −1 (Ī) is 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 (2) B 0 (y) = 1, B (y) = B −1 (y), 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) 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 in [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, where divergent integrals are regularised by the use of Hadamard finite part integrals [5]. Furthermore, a higher dimensional generalisation of the EM for simple lattice polytopes has been found [7]. 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 [8].
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 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. The details of this derivation, including the proofs of rather technical lemmas, are shown in Section 5. Finally, in Section 6, 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 where x ∈ Z, s ∈ C ∞ (R * ) 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). 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 [9,Sec. 3.2].
Definition 1 (Asymptotically smooth functions). A function s ∈ C ∞ (R * ) is called asymptotically smooth if there exist c > 0 and γ ≥ 1 such that for all y ∈ R * and ∈ N. We denote the vector space of all asymptotically smooth functions by S.
This set includes entire functions like polynomials, but also a broad set of functions with singularities. Typical examples for asymptotically smooth functions are (5) s(y) = |y| −ν , y ∈ R * , for ν ∈ R which are singular for ν > 0.
We can further classify asymptotically smooth functions by their growth rate at infinity.
Remark 2. From Grönwall's lemma follows immediately that See Section 5 for a proof.
From Definition 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 (5). This challenge is overcome by using an exponentially decaying regularisation of the interaction. Notation 1. Let s ∈ S. The exponentially weighted interaction reads (8) s β (y) = s(y)e −β|y| , y ∈ R * , with β ≥ 0.
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.
The basic object from which the SEM expansion is deduced is the function C.
Definition 3. Let s ∈ S. We define C : R + × R + → C as In the following, calligraphic symbols indicate an explicit dependence on the interaction s. The function C quantifies the difference between sum and integral of the regularised interaction. It can be used to generate a replacement for the periodised Bernoulli polynomials in (1), in which we encode all information about the interaction s.

Definition 4 (Bernoulli-A functions). Let s ∈ S. The Bernoulli-A functions,
are defined as the coefficients in the power series We say that (A (ξ)) ∈N is exponentially generated by and refer to G ξ as the generating function.
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.
Remark 3. For s = 1, we recover the periodised Bernoulli polynomials, Proof. For y > 0 and β > 0, The first term in brackets is the exponential generating function for the Bernoulli polynomials evaluated at 1 + y − y [10, Sec. 1.13, Eq. (2)]. Since B 0 = 1, we have We now define the SEM differential operator as follows: Definition 5 (SEM operator). For s ∈ S and ξ ∈ R + , we define the differential operator of infinite order with D is the derivative operator. We call D ξ the SEM operator. It formally reads and for ∈ N the finite order approximations D ( ) ξ are given by We present the first theorem of this work, the finite order SEM 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 1, the function g has to belong to the set of functions of exponential type. This sets a growth condition on its derivatives. Definition 6. Let g be entire. If there exists σ > 0 such that for every > 0, there is M > 0 with we say that g is of exponential type σ. By E σ we denote the vector space of all functions of exponential type σ.
The classical paper [11] gives an exhaustive review of functions of exponential type. The admissible range for σ depends on s, in particular on γ from Definition 1.
where s ∈ S and g ∈ E σ with σ < 2π/(1 + γ). Then, Remark 4. For simplicity, we have formulated the Theorems 1 and 2 such that x is positioned to the left of the interval [a + 1, b]. The Theorems can however also be applied in case that x is positioned to the right of this interval by simply reflecting it about x. Consider a < b < x. Then the reflected interval is [ã + 1,b], wherẽ 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 (4).

Remark 5. For the prototypical asymptotically smooth interactions
with ν ∈ R, the upper bound on σ can be improved to 2π. This is a consequence of the next example: the radius of convergence for the generating function G y is 2π. The functions (A ) ∈N are given by with ζ(·, ·) the Hurwitz zeta function, and analytically continued to the complex plane for z = 1 [10, Sec. 1.10]. For an integral ν, the coefficients are well-defined in the limit, where γ e is the Euler-Mascheroni constant and H k denotes the kth harmonic number,

Numerical Application
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. The SEM expansion naturally provides the answer to the question how to correctly include the discreteness of a 1D crystal within a continuum formulation that avoids discrete lattice sums and is therefore numerically feasible for all asymptotically smooth interaction forces. We consider a particularly difficult scenario, the case where the interaction potential is the 3D Coulomb repulsion, which decays algebraically with an exponent equal to the system dimension. Then the discreteness Figure 3. Forces F as a function of distance x in centre of a kink for different choices of the kink width λ and the particle number N . The red line shows the first order approximation of the singular Euler-Maclaurin expansion, the blue dots display the exact forces. In panel (a), the approximation error in the maximum norm over the whole chain is smaller than 3 × 10 −7 and the relative error is smaller than 8 × 10 −5 . The black line shows the approximation of the discrete sum by an integral only.
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, (14) x j = jh + u(jh), j ∈ −N, . . . , N , 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 All physical dimensions are from now on removed, where we write positions in units of h and forces in units of V (h)h. Then the forces follow as where the function f factors into with s ∈ S and g smooth such that .
For our numerical study, we make the parameter choice ν = 1, corresponding to the 3D Coulomb interaction restricted to 1D. The displacement function is chosen as the integral of a normalised Lorentzian which describes the simplified profile of a kink, an extended defect in the crystal, where Kinks typically arise when an additional nonlinear potential is applied to the crystal. The parameter λ > 0 controls the width of the kink. For further details regarding kinks in condensed matter physics, see [12]. We compute the forces in (19) by using the SEM expansion in Theorem 1 up to order , choosing δ = 1. The differential operator is first evaluated symbolically and is then subsequently applied to the function g. The runtime for the adaptive numerical integration and for the computation of the derivatives is then essentially independent of N for a single force evaluation. The SEM is implemented in Mathematica and a working example is provided along with this article.
We apply the SEM expansion to the computation of the full nonlinear long-range Coulomb forces in a crystal of macroscopic size. In Fig. 3 we first display the forces in the centre of the kink in case of a microscopic crystal with λ = 10 and N = 1000 in (a) and then for a macroscopic crystal with λ = 10 5 and N = 10 10 in (b). For a typical lattice constant h ≈ 10 −10 m, the crystal in (b) exhibits a total length of two metres. The first order SEM (red line) is compared to the approximation of the sum by an integral only (black line). The blue dots in (a) show the exact discrete particle forces. The calculation of the exact forces in (b) is not numerically possible anymore due to the large number of particles. We find that the SEM reproduces the exact forces in (a) very precisely both at the chain edges as well as inside the kink in the centre. The absolute error is less than 3 × 10 −7 for all particles with 4 digits of precision. The integral approximation however shows a significant error. In (b) both the particle number and the size of the kink are increased to macroscopic scales. Even at the macro scale, the SEM shows a visible and important difference to the integral approximation. The maximum forces at the kink centre scale as log(λ)/λ 2 while the first order SEM contribution in the centre scales as λ −2 . Both are independent of N as N → ∞. As the logarithm increases too slowly with λ in order to dominate the first order SEM contribution, the discrete particle nature of the crystal remains very much relevant even in the thermodynamic limit. We have therefore shown 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. [13,Eq. (5.3)]), is incorrect 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. 4. 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. 4. 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.
The analysis of the error scaling shows that an inclusion of the SEM correction is important for the correct prediction of long-range forces, even more so if finite chains with edges are considered, where the integral approximation suffers a complete break down independent of the choice of λ. A very regular scaling of the error with λ is observed, with the first order SEM already offering a λ −4 scaling.
The SEM expansion is applicable to all asymptotically smooth interactions. In particular, this includes all standard interaction forces and energies that appear in nature.

Derivation of the singular Euler-Maclaurin expansion
In this section, we lay out the proofs of the main results whilst skipping technical details. To this end, we moved the proofs of most lemmas in a separate section. Before we present the basis for the theorems, namely the zero order SEM expansion, we collect several properties of the function C and its derivatives. With above lemmas, especially the jump relations, we can formulate the zero order SEM expansion. Proof. First note that by Lemma 2, the function C(·, β) exhibits discontinuities at n ∈ N + but is smooth for y ∈ R + \ N with the properties where the weighting of s is removed by the limit β 0. Subsequently (25) is divided into two separate sums. An index shift is performed in the sum that includes the terms C(n − x − , β) resulting in the expression

Recombination of the two sums yields
where the second term results from an adjustment of the differing summation intervals. Going back to the integral on the left hand side of (22), we use property (24) and rewrite it as where the integrals have been combined to a single one by taking the limit → 0. Substracting Before continuing the transformation of the right hand side (22) into a differential operator, more properties of the function C are needed. We start with an investigation of the derivatives with respect to the weighting β, which are subsequently put into connection with antiderivatives of C with respect to y. Note that the definition below is well-defined as by Lemma 1 the function C(·, β) decays exponentially for β > 0.
The iterated antiderivatives of C can be expressed explicitly by derivatives with respect to the regularisation parameter. This form is very useful for deriving bounds on (C ) ∈N . Moreover, we are able to extend their definition to β = 0. Lemma 3. Let s ∈ S and ∈ N. The function C admits the explicit form which is is also valid in the limit β 0. In addition, C (·, β) ∈ C −1 (0, ∞) for all β ≥ 0. Relation (29) can be compactly written as Above lemma is the basis for bounds on (C ) ∈N . They are needed in the proof of Theorem 2 to perform the limit → ∞. Lemma 4. Let s ∈ S α for α ∈ R with constants c, c 0 , γ ≥ 1. Set α = max{0, α }+ 1. Then for all ∈ N, and y, β > 0, with c s > 0 only depending on s and The estimate holds in particular in the limit β 0.
As a direct consequence of Lemmas 3 and 4, we get Lemma 5. Let s ∈ S α , α ∈ R, and ∈ N. The functioñ lies in in C −1 (0, ∞) and is estimated by for y > 0. Here, the constants are the same as in Lemma 4.
We now show that (Ã ) ∈N are the Bernoulli-A functions from Definition 4.

Lemma 6.
We haveÃ Proof. We show that (Ã ) ∈N and (A ) ∈N have the same generating function. For ξ > 0 and ∈ N, we computẽ which proves We proceed with the proofs of the two main theorems. Proof of Theorem 1: Let x, a, b ∈ Z, x ≤ a < b and δ ∈ (0, 1]. Furthermore, we pick f : [a + δ, b + δ] → C which factors into for all ∈ N. For = 0, this is proved in Proposition 1. The case ≥ 1 readily follows via iterated integration by parts, where successive antiderivatives of C are given by (C k ) k∈N . In order to take the limit β 0, we note that C k (·, β) is uniformly bounded in β > 0 for all k ∈ N by Lemma 4. Since g ( +1) is continuous and therefore bounded on [a + δ, b + δ], the integrand is uniformly bounded in β and we conclude by the dominated convergence theorem that where we have used Lemma 6,

Proof of Theorem 2:
We now show that the limit → ∞ for (1) in Theorem 1 exists, which implies Theorem 2. First we analyse the behaviour of the remainder integral. We set and furthermore by definition of g, for all ε > 0. This implies As an immediate consequence of above estimates, the series converges uniformly on [a+δ, b+δ] for σ < τ and thus the limit → ∞ is well-defined.

Proof of Example 1:
For β > 0, y > 0 and s(y) = |y| −ν with ν ∈ C, where Γ(·, ·) denotes the incomplete gamma function and is subsequently continued to a meromorphic function on C × C [14, Chap. IX]. For ν ∈ R \ N, we have the following expansion for the series, which is valid for all β ∈ (0, 2π) [10, Sec. 1.11, Eq. (8)]. The incomplete gamma function admits the power series, valid for all β > 0, [10, Sec. 9.2, Eq. (5)], Substracting both terms, the singularities cancel and we get By above analysis, the radius of convergence of the series is 2π. To obtain the generating function of the coefficients (A ) ∈N , we compute e βy C(y, β) by means of the Cauchy product, This proves the form of the coefficients for ν ∈ R\N. For integral ν, above expression has an removable singularity. Given k ∈ N, we write By Eq. (9) from [10, Sec. 1.10], the first difference tends to where ψ is the digamma function [10, Sec. 1.7]. The last term is the differential quotient of the function ν → y −ν−k−1 , evaluated at k + 1. Therefore, the limit is equal to − log y. In total, we have The last term can be expressed as [10, Sec. 1.7.1, Eq. (9)], where γ e is the Euler-Mascheroni constant and H k is the kth harmonic number,

Technical Lemmas
Before we prove the remaining lemmas from Section 4, we give the proofs for the two remarks in Section 2.
For ∈ N, the th derivative of s is The product in above formula has the alternative form where ∼ means that the quotient of both sides tends to 1 for z → ∞. Applying this in case of → ∞, we get This shows that for ν ≤ 1, the quotient of the factorial and the prefactor is bounded. In case of ν > 1, it diverges algebraically which implies for all ε > 0. To summarise, we have shown with c > 0 only depending on ν and ε > 0. If ν ≤ 1, the estimate also holds for ε = 0. Albeit the definition of C is a priori not well-defined for β = 0, we can often prove β independent bounds and thereby extend the results in the limit β 0. A key role plays the following simple estimate.
Replacing y by −y in above estimate shows the result for all y ∈ R * .
We now investigate the behaviour of C both as a function in y and β. is differentiable, it suffices to show that it is differentiable on every compact subinterval [β 1 , β 2 ] of (0, ∞). The integrand is a smooth function in both variables with partial derivative Since s is asymptotically smooth, it admits at most polynomial growth, see Remark 2. Thus, there is c > 0 with |zs(z)|e −β1z ≤ c 1 + z 2 , z ≥ y. This means Because the majorant is independent of β ∈ [β 1 , β 2 ] and summable, respectively integrable on [y, ∞), we conclude that C(y, ·) is differentiable with derivative is asymptotically smooth for all ∈ N. Sinces e −β· coincides with ∂ β h(·, β) on (0, ∞) for all β > 0, we conclude inductively that C is infinitely differentiable with respect to β. Furthermore, it suffices to prove the rest of the lemma for = 0. Fix β > 0. Since s ∈ S, there are α ∈ R and c 0 > 0 with and c α,β > 0 such that Combining above estimates, we have To prove the existence of the limit β 0, we apply the EM expansion up to order ∈ N, The remainder reads For β 0, the sum of first two terms in (38) converges to so we only have to investigate R . From Remark 2 we know that there is α ∈ R with s ∈ S α . Combining this with estimate (4), we see that s ( +1) is integrable for all ∈ N larger than or equal to α = max{0, α } + 1.
With Lemma 7, the modulus of the integrand in R α is readily estimated by for all z ≥ 1. This upper bound is integrable on [1, ∞) as a product of a bounded and an integrable function. We can therefore employ the dominated convergence theorem which yields Lemma 2. Let s ∈ S and ∈ N, β > 0. From Lemma 1, we know that Using Eq. 38, we see where α ∈ R such that s ∈ S α and α = max{0, α } + + 1 as in the proof of Lemma 1. As shown in 39, above equation remains valid in the limit β 0. On (n, n + 1), n ∈ N, Eq.40 shows that ∂ β C is an antiderivative of a smooth function, where c n is a constant only depending on s and n. This shows that ∂ β C is a smooth function in y on R + \ N.
To prove the jump relation, fix n ∈ N + . For ε 1 , ε 2 ∈ (0, 1), it holds Taking the limit ε 1 , ε 2 0 yields for the right hand side In the last step, we applied the EM expansion (1) with a = n − 1, b = n and δ = 1 up to order α .
Above jump relations are the key for the proof of the alternative representation of (C ) ∈N . Lemma 3. Let s ∈ S. The base case = 0 holds by definition. Let ∈ N. We assume that the form (29) holds for , and now prove that it holds for + 1. Once we established this relation, we can extend the definition of C to β = 0 since the right hand side is well-defined for β = 0, see Lemma 1. We recall the jump relation for all n ∈ N + , k ∈ N and the formula for mixed derivatives, ∂ y ∂ k β C(y, β) = (−y) k s β (y) for all y ∈ R + \ N + . Note that these formulas hold for β ≥ 0.
The derivative of To show thatC +1 and C +1 differ at most by a constant, we need to show that both functions are continuous. Since C (·, β) is C −1 , its antiderivative C +1 (·, β) is C , and therefore continuous. We already know thatC +1 (·, β) is smooth on R + \ N + , so to show continuity, we have to investigate its behaviour at the positive integers. By the jump relation, we compute lim y nC for all n ∈ N + . Hence,C +1 (·, β) is continuous. Consequently, it differs from C +1 (·, β) only by constant which is zero because both functions tend to zero at infinity, see Lemma 1. The jump relations are also valid for β = 0, so combined with the induction hypothesis above argument shows lim β 0 C +1 (·, β) ∈ C (0, ∞).

Lemma 8.
Let s ∈ S and ∈ N. The function C (·, ·) takes the form Proof. Let s ∈ S. From Lemma 1, we know for all y > 0, β > 0. Combining this with from Lemma 3, the claim follows from direct application of the binomial theorem.
Lemma 9. Let s ∈ S with constants c > 0, γ ≥ 1. For y > 0, ξ ≥ 0 and y > ξ, we have for y ∈ R * , k, ∈ N with k > , and β > 0. For ξ = y, we have Proof. For y > 0 we compute In case of ξ = y, we obtain (45) from Lemma 7. Otherwise Then by Definition 1, we obtain for s ∈ S s (k) Lemma 4. Let s ∈ S α with α ∈ R with constants c > 0 and γ ≥ 1. In the following, c s > 0 denotes a generic constant that only depends on s and may change between different equations. For ∈ N, ξ > 0 and β ≥ 0, the function R * → C, y → s ξ, ,β (y) = 1 ! (y − ξ) s β (y), is asymptotically smooth and belongs to S α+ , which holds in particular for β = 0. We take the explicit form of the antiderivatives C (·, β) from Lemma 8, We apply the EM expansion to the right hand side up to order k α = α + with α = max{0, α } + 1, which yields y, ,β (z) dz.
We now derive uniform bounds β for C (y, β). We first consider y ∈ N + and then extend the result to y ∈ R + by means of a Taylor expansion. We have y , ,β ( y ) y , ,β (z) dz.

Outlook
We have tried to make this article as accessible as possible to a large audience without sacrificing mathematical rigour. It is our sincere hope that you, dear reader, will find the tools and results provided in this publication useful. We encourage you to develop the results further or use them to build something of interest. On the application side, the SEM expansion allows for a precise and fast evaluation of macroscopic lattice sums in solid state physics; here spin lattices come to mind. On the theory side, various extensions of the SEM are possible. Our next publication will be devoted to an extension of the SEM expansion to higher spatial dimensions. As a separate project, we aim at the development of fast and accurate algorithms for time evolutions of macroscopic crystals. Another point of interest is the determination of stationary states through the solution of the integro-differential equations that arise from the application of the SEM.