Singular Euler-Maclaurin expansion on multidimensional lattices
SSINGULAR EULER–MACLAURIN EXPANSION ONMULTIDIMENSIONAL LATTICES
ANDREAS A. BUCHHEIT AND TORSTEN KE(cid:223)LER
Abstract.
We extend the classical Euler–Maclaurin expansion to sums overmultidimensional lattices that involve functions with algebraic singularities.This offers a tool for the precise quantification of the effect of microscopicdiscreteness on macroscopic properties of a system. First, the Euler–Maclaurinsummation formula is generalised to lattices in higher dimensions, assuming asufficiently regular summand function. We then develop this new expansionfurther and construct the singular Euler–Maclaurin (SEM) expansion in higherdimensions, an extension of our previous work in one dimension, which remainsapplicable and useful even if the summand function includes a singular functionfactor. We connect our method to analytical number theory and show thatall operator coefficients can be efficiently computed from derivatives of theEpstein zeta function. Finally we demonstrate the numerical performance ofthe expansion and efficiently compute singular lattice sums in infinite two-dimensional lattices, which are of high relevance in solid state and quantumphysics. An implementation in Mathematica is provided online along with thisarticle. Introduction
Discrete particle systems with long-range interactions appear abundantly in na-ture. For instance, the microscopic electromagnetic interactions between atomsand molecules give rise to the macroscopic properties of a solid. The universe itselfis composed of subatomic particles and their interactions determine the evolutionof the whole. Some theories extend this notion of granularity to space-time it-self. In lattice quantum field theories, the discreteness of the space-time regularisesdivergent path-integrals and simplifies numerical predictions [30]. Other theoriesemploy discreteness at the Planck-scale in order to reconcile quantum mechanicswith general relativity, e.g. loop quantum gravity [7, 29]. Recently, it has beenconjectured that a granularity of space-time could explain dark energy, the causefor the accelerating expanse of the universe [28].Finding an efficient way to compute sums with a large number of addends, asthey appear for instance in the simulation of discrete particle systems in condensedmatter or in the evaluation of partition functions in statistical physics, is in generala challenging task. For instance, the number of summands required in the computa-tion of the energy of a lattice with long-range interactions scales quadratically withthe particle number and becomes unfeasible for solid state systems of macroscopicsize, where the particle number is in the range of N ≈ .A continuum is, in many cases, more accessible than a discrete system. Arisingintegrals can often be computed analytically or, at least, decent quadrature rulesfor numeric approximations exist. Moreover, our mathematical understanding of Mathematics Subject Classification.
Primary 65B1, 540H05; Secondary 46F10, 35B65,11E45.
Key words and phrases.
Euler–Maclaurin expansion, long-range interactions, quadrature, mul-tidimensional lattice sums, partial differential equations, elliptic regularity, analytic numbertheory. a r X i v : . [ m a t h . NA ] F e b ANDREAS A. BUCHHEIT AND TORSTEN KE(cid:223)LER continuous systems vastly exceeds the understanding of discrete systems, a notableexception being the tools developed in number theory. It is thus natural to ask, inhow far we can describe a discrete system by a related continuum, or in how far alattice sum can be approximated by an integral.The main tool for describing the difference between a sum and an integral in onedimension is the Euler–Maclaurin (EM) expansion (see [2] for a historic overview).For ‘ ∈ N , a, b ∈ Z , δ , δ ∈ (0 ,
1] and f ∈ C ‘ +1 ([ a + δ , b + δ ] , C ), we have [2, 5, 25] b X y = a +1 f ( y ) − b + δ Z a + δ f ( y ) d y = − ‘ X k =0 ( − k k ! B k +1 (1 + y − d y e ) k + 1 f ( k ) ( y ) (cid:12)(cid:12)(cid:12)(cid:12) y = b + δ y = a + δ + ( − ‘ ‘ ! b + δ Z a + δ B ‘ +1 (1 + y − d y e ) ‘ + 1 f ( ‘ +1) ( y ) d y, with d y e the ceiling of y and B ‘ the Bernoulli polynomials, which are defined bythe recurrence relation B ( y ) = 1 , B ‘ ( y ) = ‘B ‘ − ( y ) , Z B ‘ ( y ) d y = 0 , ‘ ≥ . The EM expansion has been extended to higher dimensions, where the main ap-proach is based on a tensorisation of the 1D expansion. For a higher-dimensionalequivalent of the EM expansion applicable to polynomials on simple polytopes, seee.g. [22]. Higher dimensional extensions of the EM expansion that do not rely onrepeated application of the 1D result are rare, a notable exception being the workby Müller in [26], where a two dimensional generalisation of the EM expansion wasderived that is valid for a larger set of functions and integration regions. A secondexample is the work by Freeden [13] in which a generalisation of the Müller resultto higher dimensions has been carried out. The above generalisations of the EMexpansion exhibits different advantages and disadvantages. While the results basedon tensorisation (or Todd operators) are easy to apply in practice, they are veryrestrictive in the set of functions and integration regions that are allowed. Also,they usually do not give error estimates. The non-tensorised results apply to a moregeneral set of functions and regions, however, the results are mainly theoretical anderror estimates are missing.All of the above generalisation share one critical disadvantage: They are notapplicable to summand functions, whose derivatives increase quickly with the de-rivative order, e.g. due to an algebraic singularity. Unfortunately, these functionsare of high importance in physical applications, as for instance all relevant interpar-ticle interaction, e.g. the Coulomb interaction between charged particles, belong tothis set of functions.In our previous work [5], we have developed the singular Euler–Maclaurin (SEM) expansion that makes the EM expansion applicable to physically relevant summandfunctions in one-dimension, including functions that exhibit an algebraic singularity.In this paper, we extend our previous work and generalise the SEM expansion fromone dimension to lattices in an arbitrary number of space dimensions. We avoida simple yet restrictive tensorisation of the 1D result, and offer a generalisation ofthe SEM expansion that can be applied to physically relevant interaction functionsand lattice structures. We show that our expansion can be used as a powerfulnumerical tool for the fast evaluation of forces and energies in lattices of macroscopicsize. Furthermore, we are able to precisely determine the effect of a microscopic
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 3 granularity on a macroscopic system, even if the ratio between the macro and themicro scale is astronomically large.This work is structured as follows. In Section 2, we provide a short overviewon distribution theory and elliptic regularity, which serve as important tools in thefollowing parts. We then derive a generalisation of the EM expansion to multidi-mensional lattices in Section 3, which is subsequently used as a stepping stone forarriving at the SEM expansion in higher dimensions in Section 4. Here we addressfirst the case that the singularity is positioned at a lattice point outside of theintegration region and then move on to the challenging but highly relevant casethat it is found inside of it. We then remove the only free parameter of the SEMfor interior lattice points in Section 5 by means of hypersingular integrals, whichleads to the hypersingular Euler–Maclaurin expansion. During this procedure, weunravel a deep connection of our theory to analytic number theory, which opensa way for us to efficiently compute all necessary operator coefficients. We thendemonstrate the numerical performance of the expansion in Section 6 and analysethe error. Finally, we draw our conclusions in Section 7.2.
Preliminaries
In the presentation of the material we mostly follow [20]. An extensively studythe Fourier transform and explicit expressions for many distributions can be foundin [15]. Finally, the first chapter of [33] presents an accessible introduction to ellipticregularity.2.1.
Distributions.
For an open set Ω ⊆ R d and k ∈ N = { , , , . . . } or k = ∞ , C k (Ω) denotes the set of k -times continuously differentiable functions f : Ω → C .A sequence ( u n ) n ∈ N in C k (Ω) is said to converge to u ∈ C k (Ω) if all derivatives upto order k converge compactly on Ω, i.e. for all K ⊆ Ω compact and α ∈ N d with | α | ≤ k , lim n →∞ sup x ∈ K (cid:12)(cid:12) D α u n ( x ) − D α u ( x ) (cid:12)(cid:12) = 0 . With C k (Ω) we denote the subspace of functions in u ∈ C k (Ω) whose supportsupp u = { x ∈ Ω : u ( x ) = 0 } is contained in a compact subset of Ω. Endowed with a stronger topology than theone inherited from C ∞ (Ω), C ∞ (Ω) is called the space of test functions. Its dualspace, denoted by D (Ω), is called the space of distributions on Ω. A prototypicalexample of a distribution is the Dirac distribution at x ∈ Ω which sends a testfunction to its point evaluation in x , h δ x , ψ i = ψ ( x ) , ψ ∈ C ∞ (Ω) . For p ∈ [1 , ∞ ], we set L p (Ω) as the Banach space of all measurable functions v : Ω → C with finite p -norm, k v k pp, Ω = Z Ω | v ( x ) | p d x < ∞ , for p < ∞ and k v k ∞ , Ω = ess sup x ∈ Ω | v ( x ) | < ∞ , in case of p = ∞ . Furthermore, the space of locally p -integrable functions is definedas L p loc (Ω) = (cid:8) v : Ω → C measurable : v | K ∈ L p ( K ) for all K ⊆ Ω compact (cid:9) . ANDREAS A. BUCHHEIT AND TORSTEN KE(cid:223)LER
Every function v ∈ L (Ω) defines a distribution by virtue of C ∞ (Ω) → C , ψ Z Ω v ( x ) ψ ( x ) d x . Of particular interest for us are the functions s ν : R d \ { } → C ,(2.1) s ν ( x ) = 1 | x | ν , for ν ∈ C . Clearly, s ν gives rise to a distribution on R d \ { } . Additionally, we have Theorem 2.1.
The function s ν has an extension to a distribution on R d . In casethat ν = d + 2 k , k ∈ N this extension is unique. Otherwise, for ν = d + 2 k , thereare infinitely many extensions and two of them differ by a linear combination ofderivatives of order k of δ . We denote the convolution of u ∈ D ( R d ) and ϕ ∈ C ∞ ( R d ), where one of themis assumed to have compact support, by u ∗ ϕ . It holds u ∗ ϕ ∈ C ∞ ( R d ) and( u ∗ ϕ )( x ) = u ( ϕ ( x − · )) , x ∈ R d . For v ∈ L ( R d ), we define the Fourier transform ˆ v = F v asˆ v ( ξ ) = F v ( ξ ) = Z R d e − πi h ξ , x i v ( x ) d x , ξ ∈ R d . The Fourier transform is an isomorphism on the Schwartz space S ( R d ) of rapidlydecaying smooth functions, S ( R d ) = n u ∈ C ∞ ( R d ) : sup x ∈ R d | x β D α u ( x ) | < ∞ ∀ α , β ∈ N d o . By duality, the definition of F extends to S ( R d ), the dual space of S ( R d ), calledthe space of tempered distributions. Theorem 2.2.
Let ν ∈ C . Any extension of s ν to a distribution on R d is atempered distribution. Its Fourier transform is a C ∞ -function on R d \ { } . For ν ( d + 2 N ) it holds ˆ s ν ( ξ ) = π ν − d Γ (cid:0) d − ν (cid:1) Γ (cid:0) ν (cid:1) | ξ | − d + ν , ξ ∈ R d \ { } , where Γ denotes the Gamma function. Elliptic regularity. A d -variate polynomial P of degree m ∈ N with complexcoefficients, P ( ξ ) = X α ∈ N d a α ξ α , ξ ∈ R d , is called elliptic if P m ( ξ ) = X | α | = m a α ξ α , ξ ∈ R d , does not vanish on R d \ { } . Moreover, a differential operator with constant coef-ficients is called elliptic if the associated polynomial is elliptic. Theorem 2.3.
An elliptic differential operator L with constant coefficients is hy-poelliptic, which means that if L u ∈ C ∞ (Ω) for u ∈ D (Ω) , then already u ∈ C ∞ (Ω) . Theorem 2.3 remains valid if smoothness is replaced by analyticity.
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 5
Theorem 2.4.
An elliptic differential operator L with constant coefficients isanalytic-hypoelliptic: If L u is analytic in Ω for u ∈ D (Ω) , then u is an analyticfunction on Ω . The following theorem is a generalisation of [20, Theorem 4.4.2].
Theorem 2.5.
Let ( u j ) j ∈ N be a sequence in D (Ω) that converges to u ∈ D (Ω) , lim j →∞ u j ( ψ ) = u ( ψ ) , ψ ∈ C ∞ (Ω) . Furthermore, let L be an elliptic differential operator with constant coefficients. Ifwe have L u j ∈ C ∞ (Ω) for all j ∈ N and the sequence ( L u j ) j ∈ N converges in C ∞ (Ω) to v ∈ C ∞ (Ω) , then ( u j ) j ∈ N and u are C ∞ -functions and ( u j ) j ∈ N converges to u in C ∞ (Ω) , that is compactly on Ω in all derivatives.Proof. Since L u j ∈ C ∞ (Ω) for all j ∈ N and L u = lim j →∞ L u j = v in D (Ω) , Theorem 2.3 holds u j ∈ C ∞ (Ω) and u ∈ C ∞ (Ω). We now show that alreadylim j →∞ u j = u in C ∞ (Ω) . We first choose an open neighbourhood Y ⊆ Ω of K such that there is χ ∈ C ∞ (Ω)with χ = 1 on Y . Let E ∈ D (Ω) be the fundamental solution of L which exists byTheorem 7.3.10 in [20]. If we write f j = L u j and f = L u , then u j − u = (cid:0) u j − E ∗ ( χf j ) (cid:1) − (cid:0) u − E ∗ ( χf ) (cid:1) + E ∗ (cid:0) χ · ( f j − f ) (cid:1) . On Y , we have L (cid:0) u j − E ∗ ( χf j ) (cid:1) = f j − χf j = 0 . Furthermore, since the convolution E ∗ · : C ∞ (Ω) → C ∞ (Ω)is continuous [32, Theorem 27.3], it holdslim j →∞ (cid:0) u j − E ∗ ( χf j ) (cid:1) = u − E ∗ ( χf ) in D ( Y ) . By [20, Theorem 4.4.2], above limit also holds in C ∞ ( Y ). Again, by using thecontinuity of the convolution, we conclude thatlim j →∞ E ∗ (cid:0) χ · ( f j − f ) (cid:1) = 0 in C ∞ ( Y ) . Therefore, for all α ∈ N d , D α ( u j − u ) = D α (cid:0) u j − E ∗ ( χf j ) (cid:1) − D α (cid:0) u − E ∗ ( χf ) (cid:1) + E ∗ D α (cid:0) χ · ( f j − f ) (cid:1) → K for j → ∞ . (cid:3) Band-limited functions.Definition 2.6 (Band-limited functions) . A function f : R d → C is said to beband-limited with bandwidth σ > h ∈ C (cid:0) B σ (cid:1) , f = F h. Here, B r denotes the Euclidean ball of radius r >
0. The vector space of allband-limited functions with bandwidth σ is denoted by E σ .The next lemma follows readily from the usual calculation rules for the Fouriertransform. ANDREAS A. BUCHHEIT AND TORSTEN KE(cid:223)LER
Lemma 2.7.
Let f ∈ E σ with σ > and f = ˆ h and let Ω ⊆ R d be open andbounded. Then k ∆ ‘ f k , Ω ≤ (2 πσ ) ‘ vol(Ω) k f k . Here, vol(Ω) denotes the Lebesgue measure of Ω . Euler–Maclaurin expansion in higher dimensions
In this section, we extend the EM expansion to lattices in d ∈ N + dimensions.In the following sections, we then move on to an higher dimensional extensionof the SEM expansion, allowing the summand function to also include singularfactors. Before deriving the new EM expansion, we discuss necessary definitionsand notation. We first introduce multidimensional lattices. Definition 3.1 (Lattices and related properties) . We call Λ ⊆ R d a lattice if thereexists M Λ ∈ R d × d with det( M Λ ) = 0 such thatΛ = M Λ Z d . We denote the set of all lattices in R d as L ( R d ). The elementary lattice cell E Λ isdefined as E Λ = M Λ [ − / , / d . The volume of the elementary cell is equal to the covolume V Λ of the lattice, V Λ = vol( E Λ ) = (cid:12)(cid:12) det( M Λ ) (cid:12)(cid:12) . We furthermore set a Λ > a Λ = min x ∈ Λ \{ } | x | . The dual lattice Λ ∗ , also called the reciprocal lattice, is defined asΛ ∗ = M Λ ∗ Z d , where M Λ ∗ = M −> Λ , with M −> Λ = (cid:0) M − (cid:1) > . It then holds that h y , x i ∈ Z ∀ x ∈ Λ , y ∈ Λ ∗ . Finally, we denote by n Λ ∈ N + the number of elements of Λ ∗ with norm a Λ ∗ .The dual lattice Λ ∗ is connected to Λ via the Fourier transform. Lemma 3.2 (Poisson summation formula) . Let Λ ∈ L ( R d ) and f ∈ L ( R d ) . Ifthere exist C, ε > such that (cid:12)(cid:12) f ( z ) (cid:12)(cid:12) + (cid:12)(cid:12) ˆ f ( z ) (cid:12)(cid:12) ≤ C (cid:0) | z | (cid:1) − ( d + ε ) , z ∈ R d , then V Λ X z ∈ Λ f ( z ) e − πi h z , y i = X z ∈ Λ ∗ ˆ f ( z + y ) , y ∈ R d . Proof.
The identity is well known for Λ = Z d , see Corollary 2.6 in [31, Chapter VII]or [20, Section 7.2]. For the case of a general lattice Λ, observe that for f ∈ L ( R d ), F (cid:0) f ◦ M Λ (cid:1) = 1 | det M Λ | ˆ f ◦ M −> Λ = 1 V Λ ˆ f ◦ M Λ ∗ , so the general case follows from the Poisson summation formula for Z d . (cid:3) INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 7
We subsequently introduce a new mathematical operator, the sum-integral XZ that quantifies the difference between a multidimensional lattice sum and a relatedintegral. Notation . Let Λ ∈ L ( R d ) and Ω ⊆ R d measurable. For f ∈ L (Ω) summable on Ω ∩ Λ we denote the difference between the sum of f over alllattice points in Ω and the integral of f over Ω per lattice covolume as XZ Ω , Λ f = X y ∈ Ω ∩ Λ f ( y ) − V Λ Z Ω f ( y ) d y . In case that the sum-integral is applied to longer expressions, we explicitly specifythe variable over which summation and integration take place for better readability,namely XZ y ∈ Ω , Λ f ( y ) = XZ Ω , Λ f. For ∂ Ω and f sufficiently regular, we aim at expressing the difference betweenthe lattice sum and the integral XZ Ω , Λ f as a surface integral over derivatives of f plus a remainder. Too this end, we definea higher-dimensional analogue of the periodised Bernoulli functions that appearin the one-dimensional EM expansion. The sums that appear in their definitionusually do not converge a priori, so we need to include a regularisation by meansof smooth cutoff functions as discussed below. Definition 3.4 (Mollifiers and smooth cutoff functions) . Let χ ∈ C ∞ ( B ) berotationally invariant with χ ≥ R d . For β >
0, set χ β = β − d χ ( · /β )with supp χ β ⊆ ¯ B β . We call χ β a mollifier and its Fourier transform ˆ χ β a smoothcutoff function.The basic approximation result for convolution with a mollifier is the following Lemma 3.5.
Let Ω ⊆ R d open and u ∈ C k (Ω) , k ∈ N ∪ {∞} . The convolution of u with the mollifier χ β results in the smooth function u β , u β : Ω β → C , x χ β ∗ u ( x ) = Z B χ ( y ) u ( x − β y ) d y , with Ω β = { x ∈ Ω : dist( x , ∂ Ω) > β } . Then for every β > we have that u β → u as β → in C k (Ω β ) . The following lemma is a direct consequence of ˆ χ β = ˆ χ ( β · ). Lemma 3.6.
Let ˆ χ β , β > , be a family of smooth cutoff functions. Then ˆ χ β → ˆ χ β ( ) = 1 , β → , in C ∞ ( R d ) . Lemma 3.7.
Let ˆ χ β , β ∈ (0 , , be a family of smooth cutoff functions. For all α ∈ N d there exists a constant C > with sup β ∈ (0 , (cid:12)(cid:12) D α ˆ χ β ( ξ ) (cid:12)(cid:12) ≤ C (cid:0) | ξ | (cid:1) −| α | , ξ ∈ R d . ANDREAS A. BUCHHEIT AND TORSTEN KE(cid:223)LER
Proof.
First note that as χ β = β − d χ ( · /β ), D α ˆ χ β ( ξ ) = F (cid:16) ( − πi · ) α χ β (cid:17) ( ξ ) = β | α | F (cid:16) ( − πi · ) α χ (cid:17) ( β ξ ) , ξ ∈ R d . We then set h α ( x ) = ( − πi x ) α χ ( x ), x ∈ R d , and find for γ ∈ N d and ξ ∈ R d that(2 πiβ ξ ) γ D α ˆ χ β ( ξ ) = β | α | F (cid:0) D γ h α (cid:1) ( β ξ ) . Because β ∈ (0 , (cid:12)(cid:12) ξ γ (cid:12)(cid:12) · (cid:12)(cid:12) D α ˆ χ β ( ξ ) (cid:12)(cid:12) ≤ | π | −| γ | · (cid:13)(cid:13) D γ h α (cid:13)(cid:13) , the right hand side being independent of ξ and β . The desired estimate now followsfrom (cid:0) | ξ | (cid:1) | α | (cid:12)(cid:12) D α ˆ χ β ( ξ ) (cid:12)(cid:12) ≤ d X k =1 | ξ k | ! | α | (cid:12)(cid:12) D α ˆ χ β ( ξ ) (cid:12)(cid:12) , ξ ∈ R d , noting that, after expanding the polynomial on the right hand side by the binomialtheorem, all terms are uniformly bounded in ξ and β . (cid:3) By means of the smooth cutoff functions, we can define lattice sums over thefunction s ν defined in (2.1) and study their behaviour in the limit of vanishingregularisation, β →
0. With this technique, we can now present the fundamentaltheorem of this section, from which the multidimensional EM expansion is derived.
Theorem 3.8.
Let Λ ∈ L ( R d ) and ν ∈ C . We define Z Λ ,ν : R d \ Λ → C , y V Λ ∗ lim β → X z ∈ Λ ∗ ˆ χ β ( z ) e − πi h z , y i | z | ν , where the primed sum excludes z = 0 . The function Z Λ ,ν is well-defined, i.e. thelimit exists for all y ∈ R d \ Λ and is independent of the chosen regularisation. Thefunction can be extended to a tempered distribution on R d by virtue of hZ Λ ,ν , ψ i = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) | z | ν , ψ ∈ S ( R d ) . Furthermore, the function Z Λ ,ν is analytic and the limit β → is compact in allderivatives. We split the proof into several parts. First, we show that Z Λ ,ν defines a tempereddistribution. Lemma 3.9. Z Λ ,ν as in Theorem 3.8 defines a tempered distribution.Proof. For a finite regularisation parameter β >
0, we define the auxiliary function Z Λ ,ν,β : R d → C , Z Λ ,ν,β ( y ) = V Λ ∗ X z ∈ Λ ∗ ˆ χ β ( z ) e − πi h z , y i | z | ν . The series is well-defined since ˆ χ β is a Schwartz function and thus the terms insidethe sum decay superpolynomially as | z | → ∞ . Since Z Λ ,ν,β is a bounded function,its action as a distribution is given by hZ Λ ,ν,β , ψ i = Z R d Z Λ ,ν,β ( y ) ψ ( y ) d y = V Λ ∗ X z ∈ Λ ∗ ˆ χ β ( z ) ˆ ψ ( z ) | z | ν INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 9 for a Schwartz function ψ ∈ S ( R d ). Due to | ˆ χ β | ≤ χ β → β →
0, we haveby the dominated convergence theoremlim β → hZ Λ ,ν,β , ψ i = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) | z | ν = hZ Λ ,ν , ψ i . Hence Z Λ ,ν defines a tempered distribution. (cid:3) The next lemma quantifies the convergence of terms that will arise in the proofof Theorem 3.8 after Poisson summation of the auxiliary functions.
Lemma 3.10.
Let Λ ∈ L ( R d ) and ν ∈ C with Re( ν ) > d . For a family of mollifiers χ β , β > , the functions h β : R d \ Λ β → C , y X z ∈ Λ χ β ∗ s ν ( z + y ) , with Λ β = Λ + ¯ B β , reside in C ∞ ( R d \ Λ β ) and converge to h : R d \ Λ → C , y X z ∈ Λ | z + y | − ν . in C ∞ ( R d \ Λ β ) as β → for any β > . All statements remain true if Λ isreplaced by a subset Λ of the lattice.Proof. First note that s ν has an extension to a holomorphic function ˜ s ν defined ona conic complex neighbourhood U of R d \ { } . As Re( ν ) > d , the series X z ∈ Λ ˜ s ν ( z + y )converges compactly in y on U \ Λ by the Weierstraß M-test since | ˜ s ν ( z + y ) | ≤ | z / | − Re( ν ) for sufficiently large z ∈ Λ. Therefore, h is analytic on R d \ Λ as the compact limitof analytic functions. For β > h β as h β ( y ) = χ β ∗ h ( y ) , y ∈ R d \ Λ β . Lemma 3.5 shows that h β → h in C ∞ ( R d \ Λ β ) as β → β >
0. Finally,observe that all above arguments remain valid if Λ is replaced by a subset Λ of thelattice. (cid:3) Using the previous two lemmas, we now prove the main result of this section.
Proof of Theorem 3.8.
We have previously shown in Lemma 3.9 that Z Λ ,ν definesa tempered distribution. Now, we prove that this distribution can be identified asan analytic function on R d \ Λ. We start with the auxiliary functions Z Λ ,ν,β , β > Z Λ ,ν,β ( y ) = V Λ ∗ X z ∈ Λ ∗ f β ( z ) e − πi h y , z i , y ∈ R d , with f β = ˆ χ β s ν , and transform the Dirichlet series over the reciprocal lattice intoa sum over Λ by means of Poisson summation. To this end, we at first add therestriction Re( ν ) < − ( d + 1) . Then, f β can extended to a function in C d +1 ( R d ) with f β ( ) = 0. Thus, z = cannow be included in the definition of Z Λ ,ν,β . The conditions for Poisson summation are then fulfilled as, firstly, f β inherits the superpolynomially decay of the smoothcutoff function ˆ χ β , and, secondly, |F f β ( z ) | ≤ C (cid:0) | z | (cid:1) − ( d +1) , z ∈ R d , due to f β ∈ C d +1 ( R d ). Hence by Poisson summation Z Λ ,ν,β ( y ) = X z ∈ Λ ˆ f β ( z + y ) , y ∈ R d . For y ∈ R d \ Λ β with Λ β = Λ + ¯ B β , we can express above formula in terms of theconvolution that appears in Lemma 3.10, Z Λ ,ν,β ( y ) = X z ∈ Λ χ β ∗ ˆ s ν ( z + y ) = c ν,d X z ∈ Λ χ β ∗ s d − ν ( z + y ) , where the prefactor c ν,d ∈ C is given in Theorem 2.2. The restriction on ν nowyields Re( d − ν ) > d , such that with Lemma 3.10,lim β → Z Λ ,ν,β = Z Λ ,ν in C ∞ ( R d \ Λ) , noting that for any compact set K ⊆ R d \ Λ there exists β > K ⊆ R d \ Λ β for all β < β . The lemma furthermore shows that Z Λ ,ν is analytic.As a final step, we extend the result to all ν ∈ C through elliptic regularity. For ‘ ∈ N we have(3.1) ∆ ‘ Z Λ ,ν,β = (2 πi ) ‘ Z Λ ,ν − ‘,β . Hence, we can choose ‘ large enough such thatRe( ν − ‘ ) < − ( d + 1)and find from our previous considerations that the right hand side of (3.1) convergesin C ∞ ( R d \ Λ) as β →
0. Theorem 2.5 then yields that already Z Λ ,ν,β , whichapriori converges only weakly to a distribution for β →
0, converges compactly inall derivatives to a smooth function on R d \ Λ. Finally, Theorem 2.4 shows that Z Λ ,ν is analytic. (cid:3) We then introduce the Bernoulli functions for multidimensional lattices.
Definition 3.11 (Bernoulli functions) . Let Λ ∈ L ( R d ) and ‘ ∈ N . We define theBernoulli functions B ( ‘ )Λ : R d \ Λ → R as follows B ( ‘ )Λ ( y ) = Z Λ , ‘ +1) ( y )(2 πi ) ‘ +1) . In analogy to Z Λ ,ν , they define tempered distributions via hB ( ‘ )Λ , ψ i = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) (cid:0) πi | z | (cid:1) ‘ +1) , for ψ ∈ S ( R d ). Remark . Clearly, B ( ‘ )Λ is Λ-periodic, B ( ‘ )Λ ( · + x ) = B ( ‘ )Λ , x ∈ Λ . We now introduce the central distributional property, on which the EM expan-sion in higher dimensions is based.
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 11
Proposition 3.13 (Sum-integral property of B ( ‘ )Λ ) . Let Λ ∈ L ( R d ) and ‘ ∈ N .Then for ψ ∈ S ( R d ) , (cid:10) ∆ ‘ +1 B ( ‘ )Λ , ψ (cid:11) = (cid:10) X Λ − V − , ψ (cid:11) = XZ R d , Λ ψ, where X Λ is the Dirac comb for the lattice Λ , X Λ = X z ∈ Λ δ z , and where δ z is the Dirac delta distribution.Proof. For ‘ ∈ N , the action of B ( ‘ )Λ on ψ ∈ S ( R d ) reads hB ( ‘ )Λ , ψ i = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z )(2 πi | z | ) ‘ +1) . We now compute the distributional poly-Laplacian ∆ ‘ +1 B ( ‘ )Λ , h ∆ ‘ +1 B ( ‘ )Λ , ψ i = hB ( ‘ )Λ , ∆ ‘ +1 ψ i = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) = V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) − V Λ ∗ ˆ ψ ( ) . From the Poisson summation formula follows V Λ ∗ X z ∈ Λ ∗ ˆ ψ ( z ) − V Λ ∗ ˆ ψ ( ) = X z ∈ Λ ψ ( z ) − V Λ ∗ Z R d ψ ( z ) d z . Then, as V Λ ∗ = V − , (cid:10) ∆ ‘ +1 B ( ‘ )Λ , ψ (cid:11) = (cid:10) X Λ − V − , ψ (cid:11) = XZ R d , Λ ψ. (cid:3) In the next step, we determine the maximum norm of the Bernoulli functions ofsufficiently high order ‘ , which plays an important role in the error scaling of theEM expansion in higher dimensions. Corollary 3.14 (Maximum norm of B ( ‘ )Λ ) . Let Λ ∈ L ( R d ) and ‘ ∈ N . For ‘ +1) >d , the functions B ( ‘ )Λ can be continuously extended to R d , with maximum norm kB ( ‘ )Λ k ∞ = 1 V Λ X z ∈ Λ ∗ | π z | ‘ +1) , and the scaling as ‘ → ∞ is determined by a Λ ∗ , lim ‘ →∞ (2 πa Λ ∗ ) ‘ +1) kB ( ‘ )Λ k ∞ = n Λ V Λ , with n Λ the number of elements of Λ ∗ with norm a Λ ∗ .Proof. Let k = 2( ‘ + 1) − d >
0. Then the Dirichlet series in Definition 3.11converges absolutely without regularisation on R d . Now, (cid:13)(cid:13) B ( ‘ )Λ (cid:13)(cid:13) ∞ ≤ V Λ X z ∈ Λ ∗ | π z | ‘ +1) , and inequality can be replaced by equality as the upper bound is attained on Λ.Concerning the scaling as ‘ → ∞ , observe that by the monotone convergencetheorem lim ‘ →∞ (2 πa Λ ∗ ) ‘ +1) kB ( ‘ )Λ k ∞ = 1 V Λ lim ‘ →∞ X z ∈ Λ ∗ (cid:18) a Λ ∗ | z | (cid:19) ‘ +1) = n Λ V Λ . Figure 1.
Multidimensional Bernoulli function B (0)Λ for d = 2 andΛ = Z . (cid:3) The Bernoulli function B (0)Λ for a square lattice in d = 2 dimensions is displayedin Fig. 1. It exhibits logarithmic singularities at at all lattice points, which originatefrom the fundamental solution of the Laplace equation in two dimensions.The Bernoulli functions form the coefficients of the EM differential operator. Definition 3.15 (EM operator) . For Λ ∈ L ( R d ), ‘ ∈ N , and y ∈ R d \ Λ, we definethe ‘ th order EM operator D ( ‘ )Λ , , y as D ( ‘ )Λ , , y = ‘ X k =0 (cid:16) ∇ ∆ ‘ − k B ( ‘ )Λ ( y ) − ∆ ‘ − k B ( ‘ )Λ ( y ) ∇ (cid:17) ∆ k . With D Λ , , y , we denote the infinite order EM operator obtained by setting ‘ = ∞ in above equation.We will show later that the infinite order operator is well-defined for band-limitedfunctions with bandwidth σ < a Λ ∗ .In the expansion, we will compute surface integrals that involve the SEM oper-ator over the boundary of a domain Ω. Notation . In the following a domain Ω ⊆ R d shall always denote anon-empty and connected open set with Lipschitz boundary ∂ Ω.The Bernoulli function of order ‘ can be interpreted as an infinite linear combi-nation of parametrices for the poly-Laplace operator. Here, a parametrix for ∆ ‘ +1 is a distribution E ∈ D (cid:0) R d (cid:1) with∆ ‘ +1 E = δ − ψ for a smooth function ψ [20, Definition 7.1.21]. Green’s third identity, or repre-sentation formula, also holds with a small modification if the fundamental solutionis replaced with a parametrix [21, p. 235, Eq. (20.1.6)]. We only state the resultfor ‘ = 0 where the case of general ‘ follows from successive application of Green’ssecond identity. Lemma 3.17.
Suppose E is a parametrix for the Laplacian, ∆ E = δ − ψ , with ψ ∈ C ∞ (cid:0) R d (cid:1) . Then E ∈ C ∞ (cid:0) R d \ { } (cid:1) and the following representation formula INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 13 holds for f ∈ C (cid:0) ¯Ω (cid:1) , a domain Ω ⊆ R d , and x ∈ Ω : f ( x ) − Z Ω ψ ( x − y ) f ( y ) d y = Z ∂ Ω (cid:16) ∂ n y E ( x − y ) f ( y ) − E ( x − y ) ∂ n y f ( y ) (cid:17) d S y + Z Ω E ( x − y )∆ f ( y ) d y , where ∂ n y = h n y , ∇ y i denotes the normal derivative and n y is the outward normalvector to Ω at y ∈ ∂ Ω . Furthermore, f is assumed to have compact support on ¯Ω if Ω is unbounded. We now present the EM expansion on multidimensional lattices.
Theorem 3.18 (Multidimensional EM expansion) . Let Λ ∈ L ( R d ) and Ω ⊆ R d adomain such that ∂ Ω ∩ Λ = ∅ . If f ∈ C ‘ +1) ( ¯Ω) , ‘ ∈ N , with compact support in ¯Ω in case of an unbounded domain, then the sum-integral of f over (Ω , Λ) has therepresentation XZ Ω , Λ f = Z ∂ Ω D D ( ‘ )Λ , , y f ( y ) , n y E d S y + Z Ω B ( ‘ )Λ ( y )∆ ‘ +1 f ( y ) d y . If Ω is bounded and f ∈ E σ with σ < a Λ ∗ , then XZ Ω , Λ f = Z ∂ Ω h D Λ , , y f ( y ) , n y i d S y . Proof.
From Corollary 3.14, we know that the poly-Laplacian of B ( ‘ )Λ describes thefollowing tempered distribution∆ ‘ +1 B ( ‘ )Λ = X Λ − V − = ∆ (cid:0) ∆ ‘ B ( ‘ )Λ (cid:1) . By assumption on Ω and f , the sum X z ∈ Ω ∩ Λ f ( z )has only a finite number of nonzero summands. Thus we can apply Lemma 3.17 tothe sum-integral and obtain XZ Ω , Λ f = Z ∂ Ω (cid:16) ∂ n y ∆ ‘ B ( ‘ )Λ ( y ) − ∆ ‘ B ( ‘ )Λ ( y ) ∂ n y (cid:17) f ( y ) d S y + Z Ω ∆ ‘ B ( ‘ )Λ ( y )∆ f ( y ) d y , where we have used that ∆ ‘ − k B ( ‘ )Λ is Λ-periodic and symmetric, i.e.∆ ‘ − k B ( ‘ )Λ ( z − y ) = ∆ ‘ − k B ( ‘ )Λ ( y ) , z ∈ Λ , y ∈ R d \ Λ . We apply Green’s second identity ‘ times to the right hand side and find XZ Ω , Λ f = Z ∂ Ω ‘ X k =0 (cid:16) ∂ n y ∆ ‘ − k B ( ‘ )Λ ( y ) − ∆ ‘ − k B ( ‘ )Λ ( y ) ∂ n y (cid:17) ∆ k f ( y ) d S y + Z Ω B ( ‘ )Λ ( y )∆ ‘ +1 f ( y ) d y . Finally, given a bounded domain and f ∈ E σ , we derive an estimate for the re-mainder integral R ( ‘ )Λ = Z Ω B ( ‘ )Λ ( y )∆ ‘ +1 f ( y ) d y , for 2( ‘ + 1) > d . By Corollary 3.14, B ( ‘ )Λ is continuous and bounded, thus (cid:12)(cid:12) R ( ‘ )Λ (cid:12)(cid:12) ≤ kB ( ‘ )Λ k ∞ k ∆ ‘ +1 f k , Ω . As f ∈ E σ , with f = ˆ h , we have by Lemma 2.7 k ∆ ‘ +1 f k , Ω ≤ (2 πσ ) ‘ +1) vol(Ω) k h k . After inserting the asymptotic scaling of the Bernoulli functions from Corollary 3.14,we obtain lim ‘ →∞ (cid:0) σ/a Λ ∗ (cid:1) − ‘ +1) (cid:12)(cid:12) R ( ‘ )Λ (cid:12)(cid:12) ≤ n Λ vol(Ω) V Λ k h k . Hence, (cid:12)(cid:12) R ( ‘ )Λ (cid:12)(cid:12) ∼ (cid:0) σ/a Λ ∗ (cid:1) ‘ +1) , and the remainder vanishes as ‘ → ∞ for σ < a Λ ∗ . (cid:3) The requirement of compact support on unbounded domains can be relaxed toa sufficiently fast decay at infinity.
Corollary 3.19 (EM expansion on unbounded domains) . The EM expansion oforder ‘ ∈ N extends to unbounded domains Ω ⊆ R d with ∂ Ω ∩ Λ = ∅ and functions f ∈ C ‘ +1) ( ¯Ω) for which there exist C, ε > such that |h t , ∇i k f ( y ) | ≤ C (cid:0) | y | (cid:1) − ( d + ε ) , y ∈ ¯Ω , for all t ∈ ∂B and k ≤ ‘ + 1) .Proof. Let η ∈ C ∞ ( R d ) with η ( ) = 1. For n ∈ N , we set η n = η (cid:0) · / ( n + 1) (cid:1) . Since f n = η n f has compact support, we can expand the sum-integral as XZ Ω , Λ f n = Z ∂ Ω D D ( ‘ )Λ , , y f n ( y ) , n y E d S y + Z Ω B ( ‘ )Λ ( y )∆ ‘ +1 f n ( y ) d y , for all n ∈ N . Now, the derivatives of η n are bounded independently of n , k D α η n k ∞ ≤ n + 1) | α | k D α η k ∞ ≤ k D α η k ∞ , and the bounds on the derivative of f on ¯Ω provide us with an integrable (andsummable) majorant. The EM expansion for f then follows from the dominatedconvergence theorem. (cid:3) The error of the EM expansion is measured by the remainder R ( ‘ )Λ = Z Ω B ( ‘ )Λ ( y )∆ ‘ +1 f ( y ) d y . If we interpret the EM expansion from Theorem 3.18 as a quadrature rule, we thenfind for the error under refinement, Λ h = h Λ for h >
0, of the integral approximation V Λ h |R ( ‘ )Λ h | ≤ h ‘ +1) V Λ kB ( ‘ )Λ k ∞ k ∆ ‘ +1 f k , Ω , for 2( ‘ + 1) > d , where the bound for B ( ‘ )Λ on the scaled lattice follows from Corol-lary 3.14. Similarly, if we dilate the argument of f , f λ ( x ) = f ( x /λ ), for λ > |R ( ‘ )Λ | ≤ λ − ‘ +1) kB ( ‘ )Λ k ∞ k ∆ ‘ +1 f k , Ω . INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 15
This estimate is of importance in the opposite case when the lattice sum is approx-imated by an integral, see Section 6 for an example. Special care has to be taken ifwe are interested in the convergence of R ( ‘ ) for ‘ → ∞ . The proof of above theoremshows that for band-limited functions f = F h with bandwidth E σ , the error decaysexponentially if σ < a Λ ∗ , |R ( ‘ )Λ | ≤ C Λ , Ω k h k (cid:18) σa Λ ∗ (cid:19) ‘ +1) , for a constant C Λ , Ω > f includes an algebraic singularity. In this case, a more advancedversion of the multidimensional EM expansion is required, which we develop in thenext section.4. Multidimensional singular Euler–Maclaurin expansion
Let us consider a lattice Λ ∈ L ( R d ), a bounded domain Ω ⊆ R d , and a latticepoint x ∈ Λ outside of ¯Ω. For a function f x : ¯Ω → C that factors into f x ( y ) = s ( y − x ) g ( y ) , with s ∈ C ∞ ( R d \{ } ) exhibiting an algebraic singularity at , and g ∈ C ‘ +1) ( ¯Ω),the EM expansion in Theorem 3.18 of the sum-integral XZ Ω , Λ f x does not converge as we extend the order of the expansion to infinity, even if g ∈ E σ with σ < a Λ ∗ . In the following, we call s the interaction, and g the interpolatingfunction. This lack of convergence has its origin in the algebraic singularity of s , dueto which the derivatives of s increase quickly with the derivative order. Moreover,the error of the expansion is uncontrolled and typically significant, even for loworders, as the Fourier transform of the summand function is not constrained toa ball of radius a Λ ∗ . In the following, we focus on the physically most relevantinteraction s ν = | · | − ν , with ν ∈ C .The lack of convergence of the one-dimensional EM expansion for functions withsingularities is a well known problem [1], which can be overcome by means of the1D SEM expansion [5]. In the same way as the derivation of the 1D SEM expansionrelies on the standard EM expansion in d = 1, we will make use of the EM expansionin higher dimensions in order to show the existence of mathematical objects, be itfunctions or distributions, that appear in the multidimensional SEM expansion.As we have demonstrated in 1D, the key to making the EM expansion applicableto functions that include a singular interaction lies in the inclusion of the interactionin a generalisation of the Bernoulli functions. With these singular Bernoulli func-tions, the differential operator of the expansion will act only on the well-behavedfunction g and not on s , which avoids the divergence of the remainder integral andthus ultimately leads to an expansion that is useful in practice.We structure the derivation of the SEM expansion in the following way. Inthe first step, in Section 4.1, we introduce rotationally symmetric fundamentalsolutions to the poly-Laplace operator. In Section 4.2, we proceed by suitablycombining the interaction with the fundamental solution. We call the resultingfunction Bernoulli symbol, as it shares properties with symbols that appear in thestudy of pseudo-differential operators. We then construct the singular Bernoullifunctions as regularised sum-integrals of the Bernoulli symbols in Section 4.3 and derive their properties. We define the coefficients of the SEM operator by means ofthe singular Bernoulli functions and formulate the SEM expansion in Section 4.4 forthe case that the singularity lies outside of the region Ω. Finally, in Section 4.5, weshow that singularities inside Ω, a highly relevant case in practice, can be describedby an additional local SEM operator.4.1. Fundamental solutions to the poly-Laplace operator.
The constructionof the SEM expansion relies on fundamental solutions to the poly-Laplace operator.
Notation ‘ +1 ) . Let ‘ ∈ N .We denote by φ ‘ ∈ S ( R d ) a rotationally symmetric fundamental solution to ∆ ‘ +1 ,∆ ‘ +1 φ ‘ = δ . Remark . The choice for φ ‘ is unique up to a rotationally symmetric polynomialof order 2 ‘ , which follows from the representation of the poly-Laplace operator inspherical coordinates. Lemma 4.3.
Every fundamental solution φ ‘ can be identified as a C ∞ -functionon R d \ { } . One possible choice is given by φ ‘ ( x ) = C ‘,d | x | ‘ +1) − d , ‘ ∈ N and d odd , ( φ ‘ ( x ) = C ‘,d | x | ‘ +1) − d , ‘ = 1 , . . . , m − ,φ ‘ ( x ) = C (1) ‘,d | x | ‘ +1) − d − C (2) ‘,d | x | ‘ +1) − d log | x | , ‘ ≥ m for d = 2 m, for constants C ‘,d , C (1) ‘,d , C (2) ‘,d ∈ R . Details on the constants are given in [3, Chapter I.2].
Lemma 4.4.
For φ ‘ as in Lemma 4.3 and α ∈ N d there exists a constant C > such that for all x ∈ R d \ { }| D α φ ‘ ( x ) | ≤ C | x | ‘ +1) − d −| α | (cid:16)(cid:12)(cid:12) log | x | (cid:12)(cid:12) + 1 + log( ‘ + 1) (cid:17) . Lemma 4.5.
Let Ω ⊆ R d be a bounded domain. For ‘ ∈ N and g ∈ C ‘ +1) ( ¯Ω) ,we can express g ( x ) , x ∈ Ω , by the representation formula, g ( x ) = ‘ X k =0 Z ∂ Ω (cid:16) ∂ n y ∆ ‘ − k φ ‘ ( x − y ) − ∆ ‘ − k φ ‘ ( x − y ) ∂ n y (cid:17) ∆ k g ( y ) d S y + Z Ω φ ‘ ( x − y )∆ ‘ +1 g ( y ) d y . The following lemma is a direct consequence of the estimates for the derivativesof the fundamental solutions, see [17, Lemma 4.1] for a proof in case of the Laplaceoperator, ‘ = 0. Lemma 4.6.
Let Ω ⊆ R d open and bounded. For g ∈ C ( ¯Ω) and ‘ ∈ N , the Newtonpotential f ( x ) = Z Ω φ ‘ ( x − y ) g ( y ) d y , x ∈ Ω , defines a C ‘ +1 -function on Ω . The derivatives up to order ‘ + 1 are given by D α f ( x ) = Z Ω D α φ ‘ ( x − y ) g ( y ) d y , x ∈ Ω , for α ∈ N d with | α | ≤ ‘ + 1 . INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 17
The following representation of the poly-Laplace operator in terms of higherorder directional derivatives is known as Pizetti’s formula [15, p. 74].
Lemma 4.7 (Integral representation of the poly-Laplace operator) . Let y ∈ R d .Assume g ∈ C ‘ ( U ) on some open neighbourhood U of y for ‘ ∈ N . Then ∆ ‘ g ( y ) = p ‘,d ω d Z ∂B h z , ∇i ‘ g ( y ) d S z , with ω d the surface area of the unit sphere and where the prefactor is given by p ‘,d = ( d/ ‘ (1 / ‘ . Here, ( x ) ‘ = x ( x + 1) · · · ( x + ‘ − denotes the Pochhammer symbol.Proof. Clearly, the assertion is equivalent to the identity1 ω d Z ∂B h z , ξ i ‘ d S z = 1 p ‘,d = (1 / ‘ ( d/ ‘ , ξ ∈ ∂B . By the Funk–Hecke theorem [18, Theorem 3.4.1] the integral over the sphere reducesto a one-dimensional integral,1 ω d Z ∂B h z , ξ i ‘ d S z = ω d − ω d Z − t ‘ (cid:0) − t (cid:1) ( d − / d t, which can be evaluated in terms of Gamma functions, ω d − ω d Z − t ‘ (cid:0) − t (cid:1) ( d − / d t = 1 √ π Γ( d/ (cid:0) ( d − / (cid:1) Γ (cid:0) ( d − / (cid:1) Γ( ‘ + 1 / d/ ‘ ) . Since Γ(1 /
2) = √ π , we have1 √ π Γ( d/ d/ ‘ ) Γ( ‘ + 1 /
2) = (1 / ‘ ( d/ ‘ . (cid:3) Bernoulli symbols.
We define the Bernoulli symbol as the product of theinteraction with a function that includes both the fundamental solution as well asa precisely chosen number of its higher order directional derivatives.
Definition 4.8 (Bernoulli symbol) . For ν ∈ C and ‘ ∈ N , we set a ( ‘ ) ν : (cid:0) R d \ { } (cid:1) × R d \ n ( x , x ) : x ∈ R d o , with a ( ‘ ) ν ( y , z ) = 1 | z | ν (cid:16) φ ‘ ( y − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y ) (cid:17) , and call a ( ‘ ) ν Bernoulli symbol of order ‘ for the interaction exponent ν .As will become evident at a later point, it is important to make the correct choicefor the number of higher order derivatives of the fundamental solution. If we taketoo many derivatives, integrability in the first argument is lost. If, on the otherhand, we take too few derivatives, the symbol will depend on the particular choicefor the fundamental solution φ ‘ . Lemma 4.9.
The Bernoulli symbol does not depend on the choice of the funda-mental solution φ ‘ . Proof.
We consider two choices for the rotationally symmetric fundamental solu-tion, which we denote by φ ‘, and φ ‘, . Then by Remark 4.2, the difference betweenthe two is given by a polynomial P of order 2 ‘ , φ ‘, − φ ‘, = P. We now denote by a ( ‘ ) ν, and a ( ‘ ) ν, the Bernoulli symbols for the respective fundamen-tal solutions and show that they are equal. We have a ( ‘ ) ν, ( y , z ) − a ( ‘ ) ν, ( y , z ) = 1 | z | ν (cid:16) P ( y − z ) − ‘ +1 X k =0 k ! h− z , ∇i k P ( y ) (cid:17) = 0 , as P is a polynomial of order 2 ‘ and thus equals its Taylor series of order 2 ‘ . (cid:3) We now show that the spherical surface integral of the Bernoulli symbol withrespect to its second argument vanishes, if the radius of the sphere is chosen suffi-ciently small.
Lemma 4.10.
Let ν ∈ C , ‘ ∈ N . Then for y ∈ R d \ { } Z ∂B r a ( ‘ ) ν ( y , z ) d S z = 0 , r < | y | . Proof.
Let | z | < | y | . We can then expand the first term in the Bernoulli symbol ina Taylor series in z , which leads to a ( ‘ ) ν ( y , z ) = 1 | z | ν ∞ X k =2( ‘ +1) k ! h− z , ∇i k φ ‘ ( y ) , with uniform convergence in y on B r for 0 < r < | y | . We then integrate theBernoulli symbol over a sphere with radius r , Z ∂B r a ( ‘ ) ν ( y , z ) d S z = 1 r ν ∞ X k = ‘ +1 k )! Z ∂B r h z , ∇i k φ ‘ ( y ) d S z , where we have used that terms with odd powers of z vanish in the surface integral.Using the integral representation of the poly-Laplace operator from Lemma 4.7, wefind that Z ∂B r h z , ∇i k φ ‘ ( y ) d S z = ω d p ‘,d r ( d − k ) ∆ k φ ‘ ( y ) , but now as k ≥ ‘ + 1, we have that∆ k φ ‘ ( y ) = 0 , y ∈ R d \ { } . Hence all terms in above sum vanish. (cid:3)
We now discuss the scaling of the derivatives of the Bernoulli symbol with respectto its second argument.
Lemma 4.11.
Let ‘ ∈ N , ν ∈ C , α ∈ N d , and K ⊆ R d \ { } compact. Then thereexist R > and C > such that (cid:12)(cid:12) D αz a ( ‘ ) ν ( y , z ) (cid:12)(cid:12) ≤ C | z | ‘ +1) − Re( ν ) −| α | , | z | > R, y ∈ K, where C only depends on ‘ , ν , α , and K .Proof. We first recall the definition of the Bernoulli symbol, a ( ‘ ) ν ( y , z ) = 1 | z | ν (cid:16) φ ‘ ( y − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y ) (cid:17) . INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 19
Figure 2.
Singular Bernoulli function A (0)Λ ,ν for d = 2, Λ = Z ,and ν = 2 + 10 − .In the following, we use C > ‘ , ν , α , K and whose value may change during the proof. First, we consider derivatives of theterms in brackets. Then by Lemma 4.4 (cid:12)(cid:12) D αz φ ‘ ( y − z ) (cid:12)(cid:12) ≤ C | y − z | ‘ +1)+1 − d −| α | , | y − z | > , where the exponent is increased by 1 in order to bound the logarithmic terms. For R > n , sup y ∈ K | y | o , we have that for all γ ∈ R | y − z | γ ≤ | γ | | z | γ , | z | > R. Thus (cid:12)(cid:12) D αz φ ‘ ( y − z ) (cid:12)(cid:12) ≤ C | z | ‘ +1)+ ε − d −| α | ≤ C | z | ‘ +1) −| α | , | z | > R. Moreover, (cid:12)(cid:12)(cid:12)(cid:12) D αz ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C | z | ‘ +1) −| α | . Concerning derivatives of the singular prefactor, we find (cid:12)(cid:12)(cid:12) D αz | z | − ν (cid:12)(cid:12)(cid:12) ≤ C | z | − Re( ν ) −| α | , z ∈ R d \ { } . The estimate now follows from the Leibniz rule. (cid:3)
Singular Bernoulli functions.
We now construct the singular Bernoullifunctions by applying the regularised sum-integral to the Bernoulli symbol. Thefundamental solutions then provide the singularities that lead to Dirac distribu-tions at lattice points if the poly-Laplacian is applied. Here, the regularisation ofthe sum-integral with a smooth cutoff function allows us to add and integrate overan infinite number of fundamental solutions while keeping the difference betweensum and integral well-defined.
Definition 4.12 (Singular Bernoulli functions) . Let Λ ∈ L ( R d ), ν ∈ C , and ‘ ∈ N .We define A ( ‘ )Λ ,ν : R d \ Λ → C as A ( ‘ )Λ ,ν ( y ) = lim β → XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) a ( ‘ ) ν ( y , z ) , for a family of smooth cutoff functions ˆ χ β , β >
0, and an arbitrary δ ∈ (0 , a Λ ) suchthat δ < | y | .The zero order Bernoulli function A (0)Λ ,ν is displayed in Fig. 4.3 for the two-dimensional grid Λ = Z and an interaction coefficient ν = 2 + 10 − . The Bernoullifunction exhibits the characteristic logarithmic singularities of the fundamentalsolution at all lattice points. In addition, it has a singularity at the origin y = 0that depends on the order ‘ and on the interaction coefficient ν . Finally, theasymptotic behaviour of the function is determined by the interaction. We nowshow the existence of the singular Bernoulli functions and outline their centralproperties in the following fundamental theorem. Its proof relies on the Euler–Maclaurin expansion in higher dimensions. Theorem 4.13 (Fundamental theorem of the SEM expansion) . For Λ ∈ L ( R d ) , ‘ ∈ N , and ν ∈ C , the function A ( ‘ )Λ ,ν is well-defined, and independent of the choicesfor φ ‘ , δ , and χ . Furthermore, A ( ‘ )Λ ,ν is analytic and the limit β → in the definitionof A ( ‘ )Λ ,ν is compact in all derivatives. We split the proof into several propositions and lemmas, for all of which theconditions on Λ, ‘ , ν , and ˆ χ β from Definition 4.12 shall hold. The first propositionis concerned with the well-definedness of the sum-integral in the definition of thesingular Bernoulli functions for finite β > Proposition 4.14.
For β > , the auxiliary function A ( ‘ )Λ ,ν,β : R d \ Λ → C with A ( ‘ )Λ ,ν,β ( y ) = XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) a ( ‘ ) ν ( y , z ) , | y | > δ, for < δ < a Λ is analytic and independent of the choices for δ and φ ‘ .Proof. The sum-integral is well defined due to the superpolynomial decay of ˆ χ β ( z )as | z | → ∞ . This decay also permits to interchange differentiation with the sum-integral, so that A ( ‘ )Λ ,ν,β inherits the analyticity of the Bernoulli symbol.We now show that the auxiliary function does not depend on the particularchoice for δ . To that end, let y ∈ R d \ Λ and pick δ , δ from (0 , a Λ ) smaller than | y | . Without loss of generality, we assume δ > δ . We first note that the sum overΛ is independent of the choice for δ since we require that it is smaller than a Λ , theminimal distance of two points in Λ. The difference of the sum-integral for the twochoices δ , δ is then proportional to(4.1) Z B δ \ B δ ˆ χ β ( z ) a ( ‘ ) ν ( y , z ) d z . Now we know from Lemma 4.10 that Z ∂B r a ( ‘ ) ν ( y , z ) d S z = 0 , r < | y | , so the integral in (4.1) vanishes since ˆ χ β is rotationally symmetric. This provesthat the auxiliary function does not depend of δ . Furthermore, by Lemma 4.9, the INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 21
Bernoulli symbol and thus the auxiliary function do not depend on the choice ofthe fundamental solution. (cid:3)
Proposition 4.15.
The auxiliary functions A ( ‘ )Λ ,ν,β , β > , are locally integrableon R d \ { } with Z K A ( ‘ )Λ ,ν,β ( y ) d y = XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) Z K a ( ‘ ) ν ( y , z ) d y , β > , for K ⊆ R d \ { } compact and δ > such that δ < min (cid:0) a Λ , dist( , K ) (cid:1) . Furthermore, the auxiliary functions converge in L loc ( R d \ { } ) for β → to thelocally integrable function A ( ‘ )Λ ,ν , which is independent of the choice of χ . In partic-ular, Z K A ( ‘ )Λ ,ν ( y ) d y = lim β → XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) Z K a ( ‘ ) ν ( y , z ) d y . Proof. As a ( ‘ ) ν ( · , z ) ∈ L loc ( R d ) for z ∈ R d \{ } and due to superpolynomially decayof ˆ χ β , the integral of the auxiliary function over the compact set K ⊆ R d \ { } exists and can be written as Z K A ( ‘ )Λ ,ν,β ( y ) d y = XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) Z K a ( ‘ ) ν ( y , z ) d y . In the next step, we choose
R > K ⊆ B R , dist( K, ∂B R ) >δ and such that the estimates in Lemma 4.11 hold. Additionally, we request that ∂B R ∩ Λ = ∅ . We then split the auxiliary function into two parts, A ( ‘ )Λ ,ν,β ( y ) = XZ B R \ B δ , Λ ˆ χ β a ( ‘ ) ν ( y , · ) + XZ R d \ B R , Λ ˆ χ β a ( ‘ ) ν ( y , · ) , y ∈ K, and expand the sum-integral over the unbounded domain by the EM expansion inCorollary 3.19, XZ R d \ B R , Λ ˆ χ β a ( ‘ ) ν ( y , · ) = − Z ∂B R D D ( m )Λ , , z (cid:16) ˆ χ β a ( ‘ ) ν ( y , · ) (cid:17) ( z ) , n z E d S z + Z R d \ B R B ( m )Λ ( z )∆ m +1 (cid:0) ˆ χ β a ( ‘ ) ν ( y , · ) (cid:1) ( z ) d z , with a yet to be specificied order m ∈ N . As dist( ∂B R , Λ) >
0, the integrand in thesurface integral is smooth in a neighbourhood of ∂B R . Since ˆ χ β → β → C ∞ ( R d ), the sum-integral over B R \ B δ and the surface integral over ∂B R convergein L ( K ) to y XZ B R \ B δ , Λ a ( ‘ ) ν ( y , · ) − Z ∂B R D D ( m )Λ , , z a ( ‘ ) ν ( y , · )( z ) , n z E d S z by virtue of the dominated convergence theorem. We now consider the convergenceof the remainder, R ( m ) β ( y ) = Z R d \ B R B ( m )Λ ( z )∆ m +1 (cid:0) ˆ χ β a ( ‘ ) ν ( y , · ) (cid:1) ( z ) d z , y ∈ K. Pizetti’s formula for the poly-Laplacian in Lemma 4.7 yields∆ m +1 (cid:0) ˆ χ β a ( ‘ ) ν ( y , · ) (cid:1) = p m +1 ,d ω d m +1) X k =0 (cid:18) m + 1) k (cid:19) Z ∂B h t , ∇i m +1) − k ˆ χ β h t , ∇i k a ( ‘ ) ν ( y , · ) d S t . With the uniform estimates from Lemma 3.7 for ˆ χ β and Lemma 4.11 for a ( ‘ ) ν , wethen find | ∆ m +1 (cid:0) ˆ χ β a ( ‘ ) ν ( y , · ) (cid:1) ( z ) | ≤ C m +1) X k =0 (cid:18) m + 1) k (cid:19) | z | − m +1)+ k | z | ‘ +1) − Re( ν ) − k =2 m +1) C | z | ‘ +1) − Re( ν ) − m +1) for all z ∈ R d \ B R and where C > K , m , ‘ and ν . Choosing m sufficiently large such that2( m + 1) > max (cid:8) ‘ + 1) − Re( ν ) + d, ‘ + 1) } ensures both that our upper bound for (cid:12)(cid:12) ∆ m +1 (cid:0) ˆ χ β a ( ‘ ) ν (cid:1)(cid:12)(cid:12) is integrable on K × R d \ B R and that the Bernoulli function B ( m )Λ is bounded by virtue of Corollary 3.14. Thedominated convergence theorem now shows that R ( m ) β converges in L ( K ) to y Z R d \ B R B ( m )Λ ( z )∆ m +1 z a ( ‘ ) ν ( y , z ) d z . Therefore, A ( ‘ )Λ ,ν is independent of χ and locally integrable on R d \ { } as the L -limit of locally integrable functions. (cid:3) We now show that the poly-Laplacian of the auxiliary function yields a sum-integral that includes the regularised interaction. Outside of lattice points, weidentify the distribution as a smooth function whose limit β → C ∞ ( R d \ Λ).
Proposition 4.16 (Sum-integral property) . The distributional poly-Laplacian of A ( ‘ )Λ ,ν,β , β > , on R d \ { } reads ∆ ‘ +1 A ( ‘ )Λ ,ν,β = ˆ χ β X Λ − V − | · | ν . Furthermore, ∆ ‘ +1 A Λ ,ν,β → − V Λ | · | ν , β → , in C ∞ ( R d \ Λ) and lim β → (cid:10) ∆ ‘ +1 A ( ‘ )Λ ,ν,β , ψ i = (cid:28) X Λ − V − | · | ν , ψ (cid:29) = XZ R d , Λ ψ | · | ν for all ψ ∈ C ∞ ( R d \ { } ) .Proof. For ψ ∈ C ∞ ( R d \ { } , we have (cid:10) ∆ ‘ +1 A ( ‘ )Λ ,ν,β , ψ (cid:11) = XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) D ∆ ‘ +1 ξ a ( ‘ ) ν ( ξ , z ) , ψ ( ξ ) E , INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 23 where in this context ξ denotes a placeholder with respect to which the action ofthe distribution is applied. We now insert the Definition of a ( ‘ ) ν , a ( ‘ ) ν ( y , z ) = 1 | z | ν (cid:16) φ ‘ ( y − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y ) (cid:17) , and find due to the properties of the fundamental solution of the poly-Laplacianthat D ∆ ‘ +1 ξ a ( ‘ ) ν ( ξ , z ) , ψ ( ξ ) E = ψ ( z ) | z | ν , as ψ is not supported at . Then (cid:10) ∆ ‘ +1 A ( ‘ )Λ ,ν,β , ψ (cid:11) = XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) ψ ( z ) | z | ν = XZ z ∈ R d , Λ ˆ χ β ( z ) ψ ( z ) | z | ν , as ψ has no support in B δ . The dominated convergence theorem then yields(4.2) lim β → (cid:10) ∆ ‘ +1 A ( ‘ )Λ ,ν,β , ψ i = (cid:28) X Λ − V − | · | ν , ψ (cid:29) . On R d \ Λ, the distribution ∆ ‘ +1 A ( ‘ )Λ ,ν,β can be identified as a smooth function,∆ ‘ +1 A ( ‘ )Λ ,ν,β = − ˆ χ β V Λ | · | ν . Since ˆ χ β → β →
0, in C ∞ ( R d ) by Lemma 3.6, the limit in (4.2) does not onlyhold weakly but also in C ∞ ( R d \ Λ). (cid:3)
With Propositions 4.14, 4.15, and 4.16, we now prove the fundamental theoremof the SEM expansion using elliptic regularity.
Proof of Theorem 4.13.
By Proposition 4.15 the auxiliary functions A ( ‘ )Λ ,ν,β , β > L ( R d \ Λ), and thus weakly as distributions, to A ( ‘ )Λ ,ν . Proposition 4.16shows that ∆ ‘ +1 A ( ‘ )Λ ,ν,β resides in C ∞ ( R d \ Λ) for all β > C ∞ ( R d \ Λ) to the analytic function ∆ ‘ +1 A ( ‘ )Λ ,ν , so, by virtue of Theorem 2.5, A ( ‘ )Λ ,ν,β alreadyconverges in C ∞ ( R d \ Λ) to A ( ‘ )Λ ,ν . The analyticity of the latter now follows fromTheorem 2.4. Finally, the limit function A ( ‘ )Λ ,ν is independent of δ and φ ‘ as theauxiliary functions are independent of this choice by Proposition 4.14 and it isindependent of χ due to Proposition 4.15. (cid:3) From the fundamental theorem of the SEM expansion follows the central distri-butional property of the singular Bernoulli function, on which the SEM expansionis based.
Corollary 4.17.
Let Λ ∈ L ( R d ) , ‘ ∈ N , and ν ∈ C . Let ψ ∈ C ∞ ( R d \ { } ) , then (cid:10) ∆ ‘ +1 A ( ‘ )Λ ,ν , ψ (cid:11) = XZ R d , Λ ψ | · | ν . Singular Euler–Maclaurin expansion for exterior lattice points.
Thesingular Bernoulli functions form the coefficients of the SEM differential operator.
Definition 4.18 (SEM operator) . We define the ‘ th order SEM operator D ( ‘ )Λ ,ν, y as D ( ‘ )Λ ,ν, y = ‘ X k =0 (cid:16) ∇ ∆ ‘ − k A ( ‘ )Λ ,ν ( y ) − ∆ ‘ − k A ( ‘ )Λ ,ν ( y ) ∇ (cid:17) ∆ k . We define the infinite order SEM operator D Λ ,ν, y by setting ‘ = ∞ in the aboveequation.Finally, after introducing all necessary functions and operators, we present theSEM expansion on bounded sets. Theorem 4.19 (SEM expansion) . Let Λ ∈ L ( R d ) , Ω ⊆ R d a bounded domain suchthat ∂ Ω ∩ Λ = ∅ , and x ∈ Λ \ Ω . If f x : ¯Ω → C factors into f x ( y ) = g ( y ) | y − x | ν , with ν ∈ C and g ∈ C ‘ +1) ( ¯Ω) , ‘ ∈ N , then the sum-integral of f x over (Ω , Λ) hasthe representation XZ Ω , Λ f x = Z ∂ Ω D D ( ‘ )Λ ,ν, y − x g ( y ) , n y E d S y + Z Ω A ( ‘ )Λ ,ν ( y − x )∆ ‘ +1 g ( y ) d y . Proof.
We find from the sum-integral property in Corollary 4.17 and Lemma 3.17that XZ Ω , Λ f x = Z ∂ Ω (cid:16) ∂ n y ∆ ‘ A ‘ ( y − x ) − ∆ ‘ A ‘ ( y − x ) ∂ n y (cid:17) g ( y ) d S y + Z Ω ∆ ‘ A ‘ ( y − x )∆ g ( y ) d y . The theorem follows after applying Green’s second theorem ‘ times to the righthand side. (cid:3) In contrast to case of the EM expansion, in the SEM expansion, the differentialoperator acts only on g and not on the interaction. Thus the convergence problemsof the EM expansion are avoided.4.5. Singular Euler–Maclaurin expansion for interior lattice points.
Inthe previous section, we have derived the SEM expansion for a singularity that liesoutside of the set Ω. We now move on the case that the singularity is positioned at alattice point x inside Ω. From an application point of view, this is the most relevantscenario, as it describes singular interactions inside a lattice. We will show thatthese singular interactions can be described by local derivatives of the interpolatingfunction at the position of the singularity, where the corresponding differentialoperator does not depend on Ω. The arising contribution remains relevant even inthe case of an infinite lattice without boundaries. Furthermore, as the resultingterm does not rely on oscillating surface integrals, its computation can be readilyimplemented numerically. Theorem 4.20 (SEM expansion for interior lattice points) . Assume the conditionsof Theorem 4.19, however with x ∈ Λ ∩ Ω . Let in addition ε > with ε < a Λ smallenough such that ¯ B ε ( x ) ⊆ Ω . Then XZ Ω \ B ε ( x ) , Λ f x = = D ( ‘ )Λ ,ν,ε g ( x ) + S ( ‘ )Λ ,ν g ( x ) + R ( ‘ )Λ ,ν,ε g ( x ) , with the local SEM operator = D ( ‘ )Λ ,ν,ε , = D ( ‘ )Λ ,ν,ε g ( x ) = ‘ X k =0 k )! lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) h z , ∇i k | z | ν g ( x ) , INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 25 a surface integral over derivatives of g of up to order ‘ + 1 , S ( ‘ )Λ ,ν g ( x ) = Z ∂ Ω D D ( ‘ )Λ ,ν, y − x g ( y ) , n y E d S y , and the remainder R ( ‘ )Λ ,ν,ε g ( x ) = XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) Z Ω a ( ‘ ) ν ( y − x , z )∆ ‘ +1 g ( y ) d y . To simplify the terms that appear in the proof of the SEM expansion we needthe following lemma that readily follows from the representation formula for thepoly-Laplacian and the regularity of the associated Newton potential.
Lemma 4.21.
Let ‘ ∈ N , x ∈ R d , and g ∈ C ‘ +1) (cid:0) B δ ( x ) (cid:1) for δ > . Then forany linear differential operator P of order smaller or equal ‘ +1 , we have for ε < δ P g ( x ) = Z ∂B ε ‘ X m =0 (cid:16) P ∂ n y ∆ ‘ − m φ ‘ ( y ) − P ∆ ‘ − m φ ‘ ( y ) ∂ n y (cid:17) ∆ m g ( y + x ) d S y + Z B ε P φ ‘ ( y )∆ ‘ +1 g ( x + y ) d y . The volume term vanishes in the limit ε → . We proceed with the proof of the SEM expansion for interior lattice points.
Proof of Theorem 4.20. As x ∈ Ω and Ω is open, we can choose ε > B ε ( x ) ∈ Ω. We then apply the SEM expansion in Theorem 4.19 to the sum-integralof f x over (Ω \ B ε ( x ) , Λ). After dividing the SEM surface integral into the integralover ∂ Ω and ∂B ε ( x ), we separate the terms that do not depend on ε from thosewhich do, and obtain XZ Ω \ B ε ( x ) , Λ f x = Z ∂ Ω D D ( ‘ )Λ ,ν, y − x g ( y ) , n y E d S y − S ε + R (1) ε , with S ε = Z ∂B ε ( x ) D D ( ‘ )Λ ,ν, y − x g ( y ) , n y E d S y , and with the remainder R (1) ε = Z Ω \ B ε ( x ) A ( ‘ )Λ ,ν ( y − x )∆ ‘ +1 g ( y ) d y . After inserting the definition of the SEM operator in S ε , S ε = Z ∂B ε ‘ X m =0 (cid:16) ∂ n y ∆ ‘ − m y A ( ‘ )Λ ,ν ( y ) − ∆ ‘ − m y A ( ‘ )Λ ,ν ( y ) ∂ n y (cid:17) ∆ m g ( x + y ) d S y , it becomes evident that the value of the surface integral is determined by thebehaviour of A ( ‘ )Λ ,ν around . For 0 < δ < min { ε, a Λ } , we have A ( ‘ )Λ ,ν ( y ) = lim β → XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) a ( ‘ ) ν ( y , z ) . As a ( ‘ ) ν ( y , · ) is locally integrable on R d \ { } , and as the definition of A ( ‘ )Λ ,ν isindependent of the choice of δ for δ < min {| y | , a Λ } , we can replace R d \ B δ by R d \ ¯ B ε in the sum-integral if | y | < a Λ . Now remember from the fundamental theorem of the SEM expansion 4.13 that the limit in β is compact in all derivativeswith respect to y on R d \ Λ. Therefore, we can exchange the surface integral overthe ε -sphere, including all derivatives, with the limit in β and the sum-integral. Weinsert the definition of the Bernoulli symbol, a ( ‘ ) ν ( y , z ) = 1 | z | ν (cid:16) φ ‘ ( y − x − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y − x ) (cid:17) , and note that, due to symmetry of lattice and interaction, the odd derivatives onthe right hand side do not contribute in the sum-integral. The surface integral thenreads S ε = − lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) | z | ν (cid:16) T (1) ε, z + T (2) ε, z (cid:17) , where T (1) ε, z = ‘ X k =0 k )! Z ∂B ε ‘ X m =0 (cid:16) h∇ , z i k ∂ n y ∆ ‘ − m φ ‘ ( y ) − h∇ , z i k ∆ ‘ − m y φ ‘ ( y ) ∂ n y (cid:17) × ∆ m g ( y + x ) d S y , and T (2) ε, z = − Z ∂B ε ‘ X m =0 (cid:16) ∂ n y ∆ ‘ − m φ ‘ ( y − z ) − ∆ ‘ − m y φ ‘ ( y − z ) ∂ n y (cid:17) ∆ m g ( y + x ) d S y . We first concern ourselves with T (1) ε, z . By Lemma 4.21 we can rewrite it as T (1) ε, z = ‘ X k =0 k )! h∇ , z i k g ( x ) − Z B ε h∇ , z i k φ ‘ ( y )∆ ‘ +1 g ( y + x ) d y ! . Moving on to T (2) ε, z , we find from Green’s second identity T (2) ε, z = Z B ε φ ‘ ( y − z )∆ ‘ +1 g ( y + x ) d y , | z | > ε. Thus, we obtain for the surface integral over the ε -sphere S ε = − (cid:0) = D ( ‘ )Λ ,ν,ε g ( x ) + R (2) ε (cid:1) , with R (2) ε = lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) Z B ε ( x ) a ( ‘ ) ν ( y − x , z )∆ ‘ +1 g ( y ) d y . As, by Proposition 4.15, the limit β → L ( R d \ { } ) , we can exchange the integral in R (1) ε withthe limit in β . Then R (1) ε and R (2) ε can be merged into a single remainder term R ( ‘ )Λ ,ν,ε g ( x ) = R (1) ε + R (2) ε = XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) Z Ω a ( ‘ ) ν ( y − x , z )∆ ‘ +1 g ( y ) d y . In total, the ε dependent terms can be rewritten as − S ε + R (1) ε = = D ( ‘ )Λ ,ν,ε g ( x ) + R ( ‘ )Λ ,ν,ε g ( x ) . (cid:3) INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 27 Hypersingular expansion and connection to analytic numbertheory
The local differential operator in the SEM expansion for interior lattice pointsexhibits a single free parameter ε . In many cases, it proves advantageous to re-move free parameters, as it simplifies the resulting theory and leads to a deeperunderstanding of the underlying mechanisms that make it work. In this particularcase, the removal of ε reveals a connection between the SEM expansion and ana-lytic number theory. It is this very connection that will allow us to compute thecoefficients of the local SEM operator in an efficient way.5.1. Derivation of the hypersingular Euler–Maclaurin expansion.
We aimat eliminating the free parameter ε by taking the limit ε →
0. The algebraicallydecaying interaction, however, is in general not locally integrable. Therefore, we usethe Hadamard finite-part integral in order to give meaning to otherwise divergentintegrals [23, Chapter 5].
Definition 5.1 (Hadamard finite-part integral) . Let Ω ⊆ R d be a bounded domainand x ∈ Ω. Consider a function f x : ¯Ω \ { x } → C that factors into f x ( y ) = g ( y ) | x − y | ν with ν ∈ C and g ∈ C ∞ (Ω). The Hadamard finite-part integral is then defined asthe action of a homogeneous distribution,= Z Ω f x ( y ) d y = D | · | − ν , g ( x + · ) E . If ν ∈ C \ ( N + d ), the Hadamard integral can be uniquely extended to functions g ∈ C ‘ ( ¯Ω), ‘ ∈ N , with ‘ ≥ ‘ ν,d , ‘ ν,d = b Re( ν ) − d c , and b t c the nearest integer smaller than or equal to t . The extension reads= Z Ω f x ( y ) d y = lim ε → Z Ω \ B ε ( x ) f x ( y ) d y − (cid:0) H ν,ε g (cid:1) ( x ) ! with H ν,ε = ‘ ν,d X k =0 k ! Z R d \ B ε h y , ∇i k | y | ν d y , ν ∈ C \ ( N + d ) . For ν ∈ ( N + d ), the Hadamard integral is uniquely defined up to derivatives of g of order ν − d . One possible choice is H ν,ε = ‘ ν,d − X k =0 k ! Z R d \ B ε h y , ∇i k | y | ν d y + 1 ‘ ν,d ! Z B \ B ε h y , ∇i ‘ ν,d | y | ν d y . Note that due to the spherical symmetry of | · | − ν , the non-unique term vanishes if ν is odd. Other choices for the Hadamard integral are obtained by replacing B inabove equation by an arbitrary, bounded and open neighbourhood of the origin.Hadamard integrals are also referred to as hypersingular integrals. As the result-ing expansion relies on hypersingular integrals, we refer to it as the hypersingularEuler–Maclaurin expansion (HSEM).The SEM expansion in the previous section relies on a regularisation of sumintegral by means of smooth cutoff functions. This procedure gives meaning to the difference between lattice sum and integral, even if the function is not integrable,e.g. because it increases at a polynomial rate. The HSEM expansion now requires ageneralisation of the sum-integral that makes it applicable to functions that exhibitnon-integrable algebraic singularities inside the region Ω. This is achieved in anatural way by making use of the Hadamard integral. Definition 5.2 (Hadamard sum-integral) . Let Λ ∈ L ( R d ) and Ω ⊆ R d a boundeddomain such that Ω ∩ Λ = ∅ . Let f x : Ω → C factor into f x ( y ) = g ( y ) | x − y | ν , with ν ∈ C and g ∈ C ‘ (Ω) with ‘ ≥ ‘ ν,d . We then define the Hadamard sum-integralas X = Z Ω , Λ f x = X y ∈ Ω ∩ Λ f x ( y ) − V Λ = Z Ω f x ( y ) d y , where the primed sum excludes y = x .After the introduction of the Hadamard sum-integral, we present the HSEMdifferential operator. Definition 5.3 (Hypersingular Euler–Maclaurin operator) . Let Λ ∈ L ( R d ), ν ∈ C , and ‘ ∈ N . We define the ‘ th order hypersingular Euler–Maclaurin (HSEM)operator = D ( ‘ )Λ ,ν as = D ( ‘ )Λ ,ν = ‘ X k =0 k )! lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) h z , ∇i k | z | ν . The infinite order operator = D Λ ,ν is defined by setting ‘ = ∞ in the above equation. Theorem 5.4 (Hypersingular Euler–Maclaurin expansion) . Let Λ ∈ L ( R d ) and Ω ⊆ R d a bounded domain such that ∂ Ω ∩ Λ = ∅ . For x ∈ Λ ∩ Ω , let f x : ¯Ω → C factor into f x ( y ) = g ( y ) | x − y | ν , with ν ∈ C and g ∈ C m +3 ( ¯Ω) , m ∈ N , such that m + 1) ≥ ‘ ν,d = b Re( ν ) − d c . Then for ‘ ∈ N with ‘ ≤ m , X = Z Ω , Λ f x = = D ( ‘ )Λ ,ν g ( x ) + S ( ‘ )Λ ,ν g ( x ) + R ( ‘ )Λ ,ν g ( x ) . The expansion of the sum-integral consists of the local HSEM operator = D ( ‘ )Λ ,ν , asurface integral over derivatives of g of up to order ‘ + 1 , S ( ‘ )Λ ,ν g ( x ) = Z ∂ Ω D D ( ‘ )Λ ,ν, y − x g ( y ) , n y E d S y , and the remainder R ( ‘ )Λ ,ν g ( x ) = lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) Z Ω a ( ‘ ) ν ( y − x , z )∆ ‘ +1 g ( y ) d y . INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 29
Proof.
We begin with the restructured sum-integral in Lemma 4.20, XZ Ω \ B ε ( x ) , Λ f x = = D ( ‘ )Λ ,ν,ε g ( x ) + R ( ‘ )Λ ,ν,ε g ( x ) + S ( ‘ )Λ ,ν g ( x ) , and add the rescaled Hadamard regularisation H ε,ν g ( x ) V Λ = 1 V Λ ‘ ν,d X k =0 k ! Z R d \ B ε h y , ∇i k | y | ν g ( x ) d y , to both sides. If the case 2 k = ν arises, we replace R d \ B ε by B \ B ε in theassociated integral. We now show that the HSEM expansion follows as ε →
0. Theleft hand side is the definition of the Hadamard integral,lim ε → XZ Ω \ B ε ( x ) , Λ f x + H ε,ν g ( x ) V Λ ! = X = Z Ω , Λ f x . On the right hand side, we divide the Hadamard regularisation into two parts, H ε,ν g ( x ) = H (1) ε,ν g ( x ) + H (2) ε,ν g ( x ) , where H (1) ε,ν g ( x ) includes the directional derivatives of g of order smaller or equal2 ‘ + 1 and H (2) ε,ν g ( x ) includes any remaining derivatives of higher order (note thatodd derivatives do not contribute due to symmetry of lattice and interaction). Weshow that the first contribution is absorbed in the HSEM operator while the secondregularises the remainder. The first limit yieldslim ε → (cid:18) = D ( ‘ )Λ ,ν,ε + H (1) ε,ν g ( x ) V Λ (cid:19) = lim ε → ‘ X k =0 k )! lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) h∇ , z i k | z | ν g ( x ) + H (1) ε,ν g ( x ) V Λ ! = ‘ X k =0 k )! lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) h∇ , z i k | z | ν g ( x )= = D ( ‘ )Λ ,ν g ( x ) , where we have used that derivatives of ˆ χ β do not contribute in the Hadamardintegral as ˆ χ β → β → C ∞ ( R d ). We now rewrite the remainder as R ( ‘ )Λ ,ν,ε ( x ) = lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) | z | ν h x ( z ) , with the auxiliary function h x ( z ) = Z Ω (cid:16) φ ‘ ( y − x − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y − x ) (cid:17) ∆ ‘ +1 g ( y ) d y , and show that the appropriate Hadamard regularisation for the sum-integral in theremainder is already provided by the second Hadamard regularisation, namely H ε,ν h x ( ) = H (2) ε,ν g ( x ) . Notice that h z , ∇i k h z ( ) = 0 , k = 0 , . . . , ‘ + 1 , as we have substracted a truncated Taylor expansion of order 2 ‘ + 1 from thefundamental solution. Hence these orders are already regularised. Then by Pizetti’s formula in Lemma 4.7, we can rewrite the second Hadamard regularisation in termsof poly-Laplace operators H (2) ε,ν g ( x ) = b ‘ ν,d / c X k = ‘ +1 k )! Z R d \ B ε | y | k − ν p k,d d y ∆ k g ( x ) , for ν ∈ C \ ( N + d ), where we again replace R d \ B ε by B \ B ε in case that 2 k = ν .Now as k ≥ ‘ +1 we find by the representation formula for the poly-Laplace operatorthat∆ k h x ( ) = ∆ k z Z Ω (cid:16) φ ‘ ( y − x − z ) − ‘ +1 X k =0 k ! h− z , ∇i k φ ‘ ( y − x ) (cid:17) ∆ ‘ +1 g ( y ) d y (cid:12)(cid:12)(cid:12) z = = ∆ k − ( ‘ +1) z ∆ ‘ +1 z Z Ω φ ‘ ( y − x − z )∆ ‘ +1 g ( y ) d y (cid:12)(cid:12)(cid:12) z = = ∆ k − ( ‘ +1) ∆ ‘ +1 g ( x + z ) (cid:12)(cid:12)(cid:12) z = = ∆ k g ( x ) , as g ∈ C m +3 , with 2 m + 3 > k + 1). After evaluating the remaining ε -limit,lim ε → (cid:18) R ( ‘ )Λ ,ν,ε g ( x ) + H (2) ν,ε g ( x ) V Λ (cid:19) = lim ε → lim β → XZ z ∈ R d \ ¯ B ε , Λ ˆ χ β ( z ) | z | ν h x ( z ) + H ν,ε h x ( ) V Λ ! = lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) | z | ν h x ( z ) , we have recovered the HSEM expansion. (cid:3) Connection to analytic number theory.
We now establish a connectionbetween the coefficients of the HSEM differential operator and analytic numbertheory. To this end, we first provide an alternative representation of = D Λ ,ν . Theorem 5.5.
Let Λ ∈ L ( R d ) and ν ∈ C . We set Z (0)Λ ∗ ,ν : R d \ Λ ∗ → C , Z (0)Λ ∗ ,ν ( y ) = lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) e − πi h z , y i | z | ν , where the function depends on the choice of the Hadamard regularisation in casethat ν ∈ (2 N + d ) . For all choices, Z (0)Λ ∗ ,ν can be extended to an analytic functionon R d \ Λ ∗ ∪ { } and the infinite order HSEM operator admits the representation = D Λ ,ν = Z (0)Λ ∗ ,ν (cid:18) − ∇ πi (cid:19) , in the sense of a Taylor expansion of the function Z (0)Λ ∗ ,ν around zero. The finiteorder operators are found by truncating the Taylor expansion at the correspondingorder.Proof. The proof follows in close analogy to the one of Theorem 3.8. We firstshow that Z (0)Λ ∗ ,ν defines a distribution. The sum is discussed in the aforementionedproof, so we only need to investigate the Hadamard integral. For β >
0, we set
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 31 u β ∈ D ( R d ) with h u β , ψ i = Z R d ψ ( y ) = Z R d ˆ χ β ( z ) e − πi h z , y i | z | ν d z d y , ψ ∈ C ∞ ( R d ) . The superpolynomial decay of ˆ χ β allows the exchange of the Hadamard integralwith the integration over R d , from which we recover the Fourier transform h u β , ψ i = = Z R d ˆ χ β ( z ) ˆ ψ ( z ) | z | ν d z . Since the Hadamard integral defines an extension of the function s ν to a tempereddistribution ¯ s ν , we can rewrite the definition of u β as h u β , ψ i = (cid:10) ˆ χ β ¯ s ν , ˆ ψ (cid:11) . Now, as ˆ χ β → β → C ∞ ( R d ), it holds, by continuity of the multiplicationof a distribution with a smooth function,lim β → h u β , ψ i = (cid:10) ¯ s ν , ˆ ψ (cid:11) , ψ ∈ C ∞ ( R d ) . Hence, Z (0)Λ ∗ ,ν defines a distribution by virtue of hZ (0)Λ ∗ ,ν , ψ i = X z ∈ Λ ˆ ψ ( z ) | z | ν − V Λ (cid:10) ¯ s ν , ˆ ψ (cid:11) , ψ ∈ C ∞ ( R d ) . We subsequently show that Z (0)Λ ∗ ,ν can be identified as an analytic function on R d \ Λ ∗ ∪ { } . Under the familiar restriction Re( ν ) < − ( d + 1), Poisson summationyields X = Z z ∈ R d , Λ ˆ χ β ( z ) e − πi h z , y i | z | ν = V Λ ∗ X z ∈ Λ ∗ χ β ∗ ˆ s ν ( z + y ) , as the Hadamard integral coincides with the usual integration over R d for thischoice of ν . We already studied this type of series in the proof of Theorem 3.8,the only relevant difference being that the origin is excluded in above sum. AsLemma 3.10 also holds for Λ ∗ \ { } , the rest of the proof coincides with that ofTheorem 3.8. A Taylor expansion of Z (0)Λ ∗ ,ν at then readily yields the operatorcoefficients of = D Λ ,ν . (cid:3) We proceed by showing that the coefficients of the HSEM operator can be iden-tified as meromorphic continuations of lattice sums.
Theorem 5.6 (Meromorphic continuation of Dirichlet series) . Consider a polyno-mial P : R d → C of order m . Then for ν ∈ C with Re( ν ) > d + m , the sum-integralequals its associated Dirichlet series, lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) P ( z ) | z | ν = X z ∈ Λ P ( z ) | z | ν . The left hand side forms, as a function of ν , the meromorphic continuation of theDirichlet series for ν ∈ C , whose simple poles lie at ν = d +2 k , k ∈ N , with residues ω d V Λ (1 / k (2 k )!( d/ k ∆ k P ( ) . For the proof we need two lemmas that investigate holomorphy of the sum-integral and of the Hadamard integral with respect to ν . Lemma 5.7.
Let P : R d → C be a polynomial and fix δ > with δ < a Λ . Then ν lim β → XZ z ∈ R d \ B δ , Λ ˆ χ β ( z ) P ( z ) | z | ν is an entire function.Proof. Following the proof of Proposition 4.15 we can express the limit β → | · | − ν P by means of the EM expansion for unbounded domains of sufficientlyhigh order. The integrands are then entire functions of ν and so are the resultingintegrals. (cid:3) Lemma 5.8.
Let P : R d → C be a polynomial of degree m ∈ N and let δ > .Then for ν ∈ C \ ( d + 2 N ) , = Z B δ P ( z ) | z | ν d z = − b m/ c X k =0 ω d (1 / k (2 k )!( d/ k δ − ν +(2 k + d ) ν − (2 k + d ) ∆ k P ( ) , so the left hand side defines a meromorphic function with simple poles at ν = d +2 k , k ∈ N with residue − ω d (1 / k (2 k )!( d/ k ∆ k P ( ) . Proof.
The polynomial P is equal to its Taylor expansion of order m around , P ( z ) = m X k =0 k ! h z , ∇i k P ( ) . Inserting this into the definition of the Hadamard integral yields= Z B δ P ( z ) | z | ν d z = lim ε → Z B δ \ B ε m X k =0 k ! h z , ∇i k | z | ν P ( ) d z − Z R d \ B ε ‘ ν,d X k =0 k ! h z , ∇i k | z | ν P ( ) d z ! = − ‘ ν,d X k =0 k ! Z R d \ B δ h z , ∇i k | z | ν P ( ) d z + m X k =max { ,‘ ν,d +1 } k ! Z B δ h z , ∇i k | z | ν P ( ) d z , with ‘ ν,d = b Re( ν ) − d c . We now apply Pizetti’s formula Lemma 4.7 and obtain= Z B δ P ( z ) | z | ν d z = − b ‘ ν,d / c X k =0 k )! Z R d \ B δ | z | k − ν p k,d d z ∆ k P ( )+ b m/ c X k =max { , b ( ‘ ν,d +1) / c} k )! Z B δ | z | k − ν p k,d d z ∆ k P ( ) , where the derivatives of odd order vanish due to the rotational symmetry of R d \ B δ and B δ . The integrals on the right hand side are readily computed, yielding= Z B δ P ( z ) | z | ν d z = − b m/ c X k =0 k )! ω d p k,d δ − ν +(2 k + d ) ν − (2 k + d ) ∆ k P ( ) . The final expression for the residues follows after inserting the definition of p k,d from Lemma 4.7. (cid:3) With the previous two lemmas, we are now in the position to prove Theorem 5.6.
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 33
Proof of Theorem 5.6.
We first show that the Hadamard sum-integral is meromor-phic in ν . Define the auxiliary function P β = ˆ χ β P . Then for δ > δ < a Λ , the sum-integral can be separated as followslim β → X = Z z ∈ R d , Λ P β ( z ) | z | ν = lim β → XZ z ∈ R d \ B δ , Λ P β ( z ) | z | ν − V Λ lim β → = Z B δ P β ( z ) | z | ν d z = lim β → XZ z ∈ R d \ B δ , Λ P β ( z ) | z | ν − V Λ = Z B δ P ( z ) | z | ν d z , due to locality of the Hadamard integral and as ˆ χ β → C ∞ ( R d ). By Lemma 5.7,the first term on the right hand side defines an entire function in ν . Moreover, byLemma 5.8 the second term defines a meromorphic function with simple poles at ν ∈ ( d + 2 N ). Therefore, lim β → X = Z z ∈ R d , Λ P β ( z ) | z | ν is meromorphic in ν and the poles with associated residues are determined byLemma 5.8.We now show that the Hadamard sum-integral coincides with the Dirichlet series,in case that Re( ν ) > d + m . Under this restriction, the sum-integral convergesabsolutely without regularisation,lim β → X = Z z ∈ R d , Λ P β ( z ) | z | ν = X z ∈ Λ P ( z ) | z | ν − V Λ = Z R d P ( z ) | z | ν d z . As ‘ ν,d > m , the Hadamard integral on the right hand side vanishes,= Z R d P ( z ) | z | ν d z = lim ε → Z R d \ B ε P ( z ) | z | ν − ‘ ν,d X k =0 k ! h z , ∇i k | z | ν P ( ) ! d z = 0 , thus proving the equality of the sum-integral and the associated Dirichlet series. (cid:3) The coefficients of the HSEM operator are connected to the Epstein zeta function,which has been introduced by Epstein in [11, 12], and which can be efficientlycomputed in any number of space dimensions [8].
Definition 5.9 (Epstein zeta function) . For x , y ∈ R d , A ∈ R d × d symmetricpositive definite (s.p.d.) and ν ∈ C with Re( ν ) > d , the Epstein zeta function Z isdefined by the Dirichlet series Z (cid:12)(cid:12)(cid:12)(cid:12) xy (cid:12)(cid:12)(cid:12)(cid:12) ( A ; ν ) = X z ∈ Z d e − πi h z , y i | z + x | Aν , with | z | A = √ z > A z , and where the primed sum excludes z = − x . The Epstein zeta function can beanalytically continued to a holomorphic function in ν if not both x ∈ Z d and y ∈ Z d . If on the other hand x ∈ Z d and y ∈ Z d , then it can be continued to anmeromorphic function in ν ∈ C \ { d } with a simple pole at ν = d . We furthermoredefine the simple Epstein zeta function Z such that Z ( A ; ν ) = Z (cid:12)(cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12)(cid:12) ( A ; ν ) . The Epstein zeta function has been utilised by Emersleben [9, 10] to preciselyevaluate the electrostatic potential of ionic lattices. Series representations of theEpstein zeta function with exponentially fast convergence have been developed,among others the Chowla-Selberg formula, thus it can be computed efficiently, see[8] and references therein.
Example 5.10.
The Hadamard sum-integral offers a representation for the simpleEpstein zeta function via lim β → X = Z R d , Λ ˆ χ β | · | ν = Z (cid:0) M > Λ M Λ ; ν (cid:1) , which is a meromorphic function in ν for ν ∈ C with a simple pole at ν = d withresidue ω d V Λ = ω d q det( M > Λ M Λ ) . For the special case d = 1 and Λ = Z , the usual Riemann zeta function ζ isrecovered, lim β → X = Z R , Z ˆ χ β | · | ν = 2 ζ ( ν ) . Vice versa, any Epstein zeta function Z can be represented by a Hadamard sum-integral. To see this, assume A ∈ R d × d symmetric and positive definite. Thematrix A now admits the unique Cholesky factorisation A = L > L with a regularlower triangular matrix L ∈ R d × d . We set the lattice Λ L = L > Z d and obtain Z ( A ; ν ) = Z ( L > L ; ν ) = lim β → X = Z R d , Λ L ˆ χ β | · | ν for all ν ∈ C \ { d } .In case that we want to avoid the use of Hadamard finite-part integrals in nu-merical applications, we can apply the SEM for interior lattice points instead ofthe HSEM. In this case, the following corollary provides us with the necessary localSEM operator = D Λ ,ν,ε . Corollary 5.11.
Let ‘ ∈ N , ν ∈ C \ ( d + 2 N ) , and ε > with ε < a Λ . Then = D Λ ,ν,ε = = D Λ ,ν − ω d V Λ ‘ X k =0 (1 / k (2 k )!( d/ k ε − ν +(2 k + d ) ν − (2 k + d ) ∆ k . Proof.
The corollary is a direct consequence of Theorem 5.6 and of the representa-tion of the Hadamard integral in Lemma 5.8. (cid:3)
The observation that the Hadamard sum-integral generates meromorphic con-tinuations of multidimensional Dirichlet series provides a connection of our work toanalytic number theory. This connection is fruitful in numerical practice, as allowsus to utilise the vast theory of multi-dimensional lattice sums (see [4] and refer-ences therein) in order to compute the coefficients of the HSEM operator efficientlyand precisely. In the following section, we first lay out how the coefficients of theHSEM operator can be calculated and then analyse the numerical performance ofthe expansion.
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 35 Numerical application
Model description.
In the following, we implement the HSEM expansionand analyse the approximation error. Instead of focussing on one particular physicalapplication, we investigate a prototypical multidimensional sum that appears in awide range of physically relevant systems, and demonstrate that it can be preciselyand efficiently reproduced. We choose to approximate X y ∈ Z f x ( y )for x ∈ Z and f x ( y ) = g ( y ) | y − x | ν . and ν ∈ C . We furthermore take a Gaussian function with width λ > g ,(6.1) g ( y ) = e −| y | /λ . Sums of this kind can be found in many different areas of condensed matterand quantum physics. Among others, they appear in the computation of forcesand energies in lattices of long-range interacting atoms, which critically determinethe properties of the resulting materials (see [14] for a recent review on long-rangeinteracting nanoscale systems). They can arise in the simulations of spin-waves,also called magnons, which can be used as information carriers [16]. Finally, theycan be found in the evaluation of certain partition functions in quantum mechanicsand solid state physics, from which the thermodynamical properties of the systemsin question can be extracted [6].We have made the choice for a square lattice in two dimensions as to allow fora basic proof-of-principle implementation that can be easily verified and modified.A reference implementation of the HSEM expansion in Mathematica is providedonline alongside with this article . The method can readily be applied to higherdimensions after the technical task of implementing the Epstein zeta function andits derivatives in higher dimensions, following [8].6.2. Efficient computation of HSEM operator coefficients.
We briefly dis-cuss how to determine the coefficients of the HSEM operator in multidimensionallattices. The general approach in d dimensions is to implement the function Z (0)Λ ∗ ,ν from Proposition 5.5 and compute the HSEM operator coefficients from its Taylorseries at , = D Λ ,ν = Z (0)Λ ∗ ,ν (cid:18) − ∇ πi (cid:19) = ∞ X k =0 k )! lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) h z , ∇i k | z | ν , using the exponentially convergent summation formulas in [8] and their derivatives.From the the meromorphic continuation theorem 5.6, we know that if the Hadamardsum-integral converges without β -regularisation it is equal to its associated Dirichletseries, lim β → X = Z z ∈ R d , Λ ˆ χ β ( z ) h z , ∇i k | z | ν = X z ∈ Λ h z , ∇i k | z | ν . Otherwise the regularised Hadamard sum-integral generates the meromorphic con-tinuation of the Dirichlet series. The case k = 0 yields the zero order HSEM https://github.com/andreasbuchheit/hsem Figure 3.
Absolute error of HSEM expansion in the maximumnorm for an interaction function s ( y ) = | y | − and an interpolatingfunction g as in (6.1) as a function of the scaling parameter λ fordifferent expansion orders ‘ .operator coefficient, which is in many applications the dominant contribution. Itcorresponds to a simplified Epstein zeta function,lim β → X = Z R d , Λ ˆ χ β | · | ν = Z (cid:0) M > Λ M Λ ; ν (cid:1) , with Λ = M Λ Z d , which can be efficiently evaluated in arbitrary dimensions (see[8]), and has been analytically determined for cubic lattices in some dimensions.6.3. Results and discussion.
We now approximate the two-dimensional sum of f x over Z \ { x } by the respective Hadamard integral plus the HSEM operator oforder ‘ , X y ∈ Z f x ( y ) ≈ = D ( ‘ )Λ ,ν g ( x ) + = Z R d g ( y ) | x − y | ν d y . and analyse the error for different widths λ of the interpolating function. Here,the implementation of the HSEM operator in two dimensions is discussed in Ap-pendix A.1. For details on the computation of the Hadamard integral, see Appen-dix A.2.In this formulation an additional feature of the HSEM becomes apparent. Whereasthe approximate computation of the sum on the left hand side requires summationover a large number of particles N with a computational time that increases linearlywith the number of terms, the computation of the right hand side is practically inde-pent of N or at most scales sublinearly with N . The first term is a local differentialoperator whose coefficients depend on the lattice and can be reduced to evaluationsof derivatives of the Epstein zeta function. For a large class of interpolating func-tions, the Hadamard integral over the whole space may be computed analytically.If this is not possible then we can exploit the convolution structure of the integraland the fact that in the case of a band-limited interpolating function, quadraturerules converge exponentially in the number of quadrature nodes as these functionsare entire and band-limited. INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 37
As the scaling coefficient of the interaction, we choose the inverse square interac-tion with ν = 2 = d , which appears in many different quantum mechanical models(see Ref. [19] and references therein), and which forms the most numerically chal-lenging model, as both short and long-range contributions to the sum in generalremain relevant and neither can be discarded. The absolute error in the maximumnorm over Z is displayed in Fig. 3 for different orders ‘ of the HSEM operator asa function of the width λ of the interpolating function g . The same error scalingis found for the Coulomb interaction with ν = 1 and the dipole interaction with ν = 3. The corresponding plots are provided online in our github repository. Theconvergence graphs for other ν ∈ C can be efficiently generated by adapting therespective parameter for ν in the Mathematica code.We find that, to good approximation, the error scales as E ‘ ( λ ) ∼ λ − ‘ +1) , where the exact scaling coefficients obtained from a linear fit are given in Fig. 3.This scaling law is identical to the one that we predict for the EM expansion appliedto a band-limited function of sufficiently small bandwidth. Hence, the singularityof f x has been well absorbed in the coefficients of the HSEM operator, restoring theconvergence properties of the expansion. Thus, as λ → ∞ , the easily computablezero order HSEM contribution already yields a good approximation to the sum. Onthe other hand, approximations that simply replace the sum by an integral yield anerror that is independent of λ , and are therefore unreliable. For increasing HSEMorders ‘ , the EM scaling law for band-limited functions is well obeyed. Alreadyfor ‘ = 6 and λ = 10, an absolute error of less than 10 − is reached. The goodconvergence properties can be explained by the fact that the Fourier transformof ∆ ‘ +1 g has its mass concentrated inside the unit ball, with only a small highfrequency correction. Thus, the function essentially behaves as if it belonged to E σ with σ < a Λ ∗ = 1 and only for very large orders ‘ , the convergence breaks down.The HSEM works well, as long as the interpolating function g , which in many casesis an interpolation of discrete particle positions, does not exhibit oscillations in spacewhose inverse wavelengths exceed a Λ ∗ . These oscillations would exhibit wavelengthsthat are smaller than the unit lattice cell and would therefore be unphysical. Thus,the HSEM converges for physically meaningful interpolation functions g , allowing usto approximate singular sums in arbitrary dimensions independently of the particlenumber. 7. Conclusions and outlook
In this work, we have extended the Euler–Maclaurin (EM) summation formula tolattices in higher dimensions. We have then subsequently used the EM expansion asa tool in the proof of the multidimensional singular Euler–Maclaurin (SEM) expan-sion that allows for the inclusion of singular factors in the summand function andthus makes the formula applicable to interaction functions that appear in physicalapplications. This is for instance the case in a solid state system with Coulomb ordipolar particle interactions. In order to avoid the evaluation of oscillatory surfaceintegrals and make the expansion easy to use in practice, we have gone one step fur-ther and have introduced the hypersingular Euler–Maclaurin (HSEM) expansion,where the differential operator is local with coefficients that can be evaluated usingstandard techniques from number theory. We have designed the SEM and HSEMexpansion as mathematical tools that can immediately be applied to open questionsin physics. It is our hope that the SEM and HSEM expansion will find use in theprecise evaluation of long-range forces and energies in crystal and spin lattices, in We include a numerical offset of 10 − in order to facilitate the implementation. the evaluation of high-dimensional partition functions in statistical physics, and inthe quantification of discreteness effects in fundamental physics.We expect that the summation formulas developed in this paper can be furthergeneralised to quasi-crystals, like the Penrose lattice, and to statistical distributionsof particles. We also consider it to be worthwhile to investigate solution techniquesfor the integro-differential equations that follow from the HSEM expansion, forinstance in the context of spin waves. Acknowledgements
We would like to thank our colleagues Daniel Seibel, Christian Michel, and PeterSchuhmacher for proof-reading the manuscript and for helpful suggestions. Wethank Prof. Sergej Rjasanow for insightful discussions and for his support.
Appendix A. HSEM expansion in two dimensions
A.1.
HSEM operator coefficients.
For a square lattice in d = 2 dimensions,a simple approach for generating the HSEM operator coefficients is available thatavoids derivatives of Epstein zeta functions and uses efficient summation formulasthat have been found in the analysis of the Riemann hypothesis in higher dimensions[24]. For a two-dimensional square lattice, it has been shown that [34, Eq. (9)] X z ∈ Z | z | ν = 4 ζ ( ν/ β D ( ν/ , where β D is the Dirichlet beta function and where the Dirichlet series can be ex-tended to ν ∈ C \ { } . Furthermore, the meromorphic continuation of the latticesum C n ( ν ) = X z ∈ Z z n | z | ν +2 n , ν ∈ C \ { } , has been shown to be computable for n ∈ N + via the formula [24, Eq. (2.3)] C n ( ν ) = 2 √ π Γ( ν/ n − / ζ ( ν − ν/ n )+ 8 π ν/ Γ( ν/ n ) ∞ X z =1 ∞ X z =1 (cid:18) z z (cid:19) ( ν − / ( z z π ) n K ( ν − / n (2 πz z ) , with K ν ( x ) the modified Bessel function of the second kind. As the double sumconverges exponentially fast in both variables, the lattice sum can be efficientlyapproximated. We now show that, using the above two lattice sums, we can gener-ate the whole HSEM operator in d = 2 dimensions by using an expansion in solidharmonics.First note that for k >
0, only even k lead to a nonzero contribution due tosymmetry of the interaction. Setting k = 2 n , with n ∈ N , we then findlim β → X = Z z ∈ R , Z ˆ χ β ( z ) h z , ∇i n | z | ν = n X m =0 a (2 n )2 m lim β → X = Z z ∈ R , Z ˆ χ β ( z ) | z | n − m ) A m ( z ) | z | ν A n − m ) ( ∇ )∆ n − m , with the solid harmonic A k : R → R , A k ( y ) = Re (cid:16) ( y + iy ) k (cid:17) , INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 39 and a ( k )0 = 12 π π Z cos( φ ) k d φ, a ( k ) n = 1 π π Z cos( nφ ) cos( φ ) k d φ. Finally,(A.1) X z ∈ Z A m ( z ) | z | ν +2 m = Z ( I , ν ) + m X k =1 k )! T (2 k )2 m (0) C k ( ν ) , where I ∈ R × denotes the identity matrix and T m is the Chebyshev polynomialof the first kind of order m . To prove this, we use that A m ( z ) = | z | m cos(2 mφ )for the polar angle φ , z = | z | (cos φ, sin φ ). Now, cos(2 mφ ) can be expanded intopowers of cos φ by means of the Chebyshev polynomial T m ,cos(2 mφ ) = T m (cos φ ) = m X k =0 T ( k )2 m cos( φ ) k . Inserting this into the right hand side of (A.1), observing that odd orders vanishdue to the symmetry of the lattice and furthermore | z | k (cos φ ) k = z k , yields the desired equality.A.2. Evaluation of the Hadamard integral.
We briefly discuss how the Hadamardintegral is evaluated. For the special choice of g in (6.1) and d = 2, we can determinethe Hadamard integral analytically. As g ∈ S ( R d ), we have that= Z R d g ( y ) | x − y | ν d y = F (cid:16) ( F| · | − ν )( F g ) (cid:17) ( x ) , where we have applied the convolution theorem for distributions. We then find= Z R d g ( y ) | x − y | ν d y = π Γ(1 − ν/ λ ν − M (cid:0) ν/ , , −| x /λ | (cid:1) , with M the Kummer confluent hypergeometric function [27, Eq. (13.2.2)]. References [1] T. M. Apostol,
Introduction to Analytic Number Theory , Undergraduate Textsin Mathematics, Springer New York, 1998.[2] ,
An Elementary View of Euler’s Summation Formula , The AmericanMathematical Monthly (1999), no. 5, 409–418.[3] N. Aronszajn, T.M. Creese, and L.J. Lipkin,
Polyharmonic Functions , Oxfordmathematical monographs, Clarendon Press, 1983.[4] J. Borwein, M. Glasser, R. McPhedran, J. Wan, and I. Zucker,
Lattice SumsThen and Now , Encyclopedia of Mathematics and its Applications, CambridgeUniversity Press, 2013.[5] Andreas A. Buchheit and Torsten Keßler,
Singular Euler-Maclaurin expansion ,(2020).[6] A. Campa, T. Dauxois, D. Fanelli, and S. Ruffo,
Physics of long-range inter-acting systems , OUP Oxford, 2014.[7] M. Dupuis, J.P. Ryan, S. Speziale, et al.,
Discrete gravity models and loopquantum gravity: a short review , SIGMA. Symmetry, Integrability and Geom-etry: Methods and Applications (2012), 052. [8] E. Elizalde, Zeta functions: formulas and applications , Journal of Computa-tional and Applied Mathematics (2000), no. 1-2, 125–142.[9] O. Emersleben,
Zetafunktionen und elektrostatische Gitterpotentiale. I , Phys.Z (1923), 73–80.[10] , Zetafunktionen und elektrostatische Gitterpotentiale. II , Phys. Z (1923), 97–104.[11] P. Epstein, Zur Theorie allgemeiner Zetafunktionen. I , Math. Ann. (1903),615–644.[12] , Zur Theorie allgemeiner Zetafunktionen. II , Math. Ann. (1906),205–216.[13] W. Freeden, Metaharmonic lattice point theory , CRC Press, 2011.[14] R.H. French et al.,
Long range interactions in nanoscale science , Reviews ofModern Physics (2010), no. 2, 1887.[15] I. M. Gel’fand and G. E. Shilov, Generalized functions , Academic Press, 1964.[16] M. Gibertini, M. Koperski, A.F. Morpurgo, and K.S. Novoselov,
Magnetic2D materials and heterostructures , Nature nanotechnology (2019), no. 5,408–419.[17] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Sec-ond Order , Springer-Verlag, 1998.[18] H. Groemer,
Geometric Applications of Fourier Series and Spherical Harmon-ics , Cambridge University Press, 1996.[19] K.S. Gupta,
Quantum inverse square interaction , Modern Physics Letters A (2003), 2355–2362.[20] L. Hörmander, The Analysis of Linear Partial Differential Operators I: Distri-bution Theory and Fourier Analysis , Classics in Mathematics, Springer, 2003.[21] ,
The Analysis of Linear Partial Differential Operators III: Pseudo–Differential Operators , Classics in Mathematics, Springer, 2007.[22] Y. Karshon, S. Sternberg, and J. Weitsman,
Exact Euler–Maclaurin formulasfor simple lattice polytopes , Advances in Applied Mathematics (2007), no. 1,1–50.[23] W. McLean, Strongly elliptic systems and boundary integral equations , Cam-bridge University Press, 2000.[24] R.C. McPhedran, I.J. Zucker, L.C. Botten, and N.P. Nicorovici,
On the Rie-mann property of angular lattice sums and the one-dimensional limit of two-dimensional lattice sums , Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences (2008), no. 2100, 3327–3352.[25] G. Monegato and J. N. Lyness,
The Euler-Maclaurin expansion and finite-partintegrals , Numerische Mathematik (1998), no. 2, 273–291.[26] C. Müller, Eine Verallgemeinerung der Eulerschen Summenformel und ihreAnwendung auf Fragen der analytischen Zahlentheorie , Abhandlungen aus demMathematischen Seminar der Universität Hamburg, vol. 19, Springer, 1954,pp. 41–62.[27] F.W. Olver, D.W. Lozier, R.F. Boisvert, and C.W. Clark,
Nist handbook ofmathematical functions , Cambridge University Press, 2010.[28] A. Perez and D. Sudarsky,
Dark energy from quantum gravity discreteness ,Physical Review Letters (2019), no. 22, 221302.[29] C. Rovelli and S. Speziale,
Reconcile Planck-scale discreteness and the Lorentz-Fitzgerald contraction , Physical Review D (2003), no. 6, 064019.[30] J. Smit, Introduction to quantum fields on a lattice , Cambridge UniversityPress, 2002.[31] E. M. Stein and G. Weiss,
Introduction to Fourier Analysis on EuclideanSpaces , Princeton University Press, 1972.
INGULAR EULER–MACLAURIN EXPANSION ON MULTIDIMENSIONAL LATTICES 41 [32] F. Trèves,
Topological vector spaces, distributions and kernels , Academic Press,1967, Reprinted by Dover, 2006.[33] ,
Basic linear partial differential equations , Academic Press, 1975,Reprinted by Dover, 2006.[34] I.J. Zucker,
The Exact Evaluation of Some New Lattice Sums , Symmetry (2017), no. 12, 314. Department of Mathematics, Saarland University, PO 15 11 50, D-66041 Saarbrücken
Email address : [email protected] Department of Mathematics, Saarland University, PO 15 11 50, D-66041 Saarbrücken
Email address ::