Martin Costabel
email: martin.costabel@univ-rennes1.fr

Monique Dauge
email: monique.dauge@univ-rennes1.fr

Khadijeh Nedaiasl
email: knedaiasl85@gmail.com

Stability Analysis of a Simple Discretization Method for a Class of Strongly Singular Integral Equations

Keywords: Mathematics Subject Classification. 65R20, 45E10, 47A12, 65B10, 78M99 volume integral equation, strongly singular kernel, delta-delta discretization, discrete dipole approximation, numerical stability

HAL is

STABILITY ANALYSIS OF A SIMPLE DISCRETIZATION METHOD FOR A CLASS OF STRONGLY SINGULAR INTEGRAL EQUATIONS

MARTIN COSTABEL, MONIQUE DAUGE, AND KHADIJEH NEDAIASL ABSTRACT. Motivated by the discrete dipole approximation (DDA) for the scattering of electromagnetic waves by a dielectric obstacle that can be considered as a simple discretization of a Lippmann-Schwinger style volume integral equation for time-harmonic Maxwell equations, we analyze an analogous discretization of convolution operators with strongly singular kernels.

For a class of kernel functions that includes the finite Hilbert transformation in 1D and the principal part of the Maxwell volume integral operator used for DDA in dimensions 2 and 3, we show that the method, which does not fit into known frameworks of projection methods, can nevertheless be considered as a finite section method for an infinite block Toeplitz matrix. The symbol of this matrix is given by a Fourier series that does not converge absolutely. We use Ewald's method to obtain an exponentially fast convergent series representation of this symbol and show that it is a bounded function, thereby allowing to describe the spectrum and the numerical range of the matrix.

It turns out that this numerical range includes the numerical range of the integral operator, but that it is in some cases strictly larger. In these cases the discretization method does not provide a spectrally correct approximation, and while it is stable for a large range of the spectral parameter λ, there are values of λ for which the singular integral equation is well posed, but the discretization method is unstable. Introduced almost 50 years ago by Purcell and Pennypacker [START_REF] Purcell | Scattering and adsorption of light by nonspherical dielectric grains[END_REF], the Discrete Dipole Approximation (DDA) is a classical numerical method in computational electromagnetics that is the subject of a vast and still rapidly growing literature (see the surveys [START_REF] Yurkin | The discrete dipole approximation: An overview and recent developments[END_REF][START_REF] Chaumet | The discrete dipole approximation: A review[END_REF]), but is virtually unknown in the mathematical community. It can be considered as a numerical approximation scheme for a strongly singular volume integral equation that, however, is too simple to fit into any known framework for standard approximation schemes for such equations (Galerkin, collocation or Nyström methods etc). In particular, to the authors' knowledge, there does not exist any error estimate or convergence proof for this method.

In the paper [START_REF] Yurkin | Convergence of the discrete dipole approximation. I. Theoretical analysis[END_REF], estimates for a consistency error are derived, and it is observed that to complete the convergence analysis, a uniform estimate for the inverse of the matrix of the linear system (stability estimate) would be needed. In the present paper, we prove first results on the way to such stability estimates for the DDA and related numerical schemes. The class of singular integral equations considered here includes the quasi-static case (i.e. zero frequency limit) of the Maxwell volume integral equation that describes the scattering of electromagnetic waves by a penetrable dielectric body in the case of constant electric permittivity. Further stability results for the nonzero frequency case will be the subject of a forthcoming paper.

Because of the simplicity of the class of operators considered here (convolution operators with kernels positively homogeneous of degree -d on a bounded domain Ω ⊂ R d ), we are able to obtain rather sharp results on the region of stability, by estimating the numerical range of the discretized operator in comparison with the numerical range of the integral operator. It turns out that for some operators, including the quasi-static Maxwell case in dimension d ≥ 2, the stability region is smaller than what one would naïvely expect. This corresponds to the fact that the eigenvalues of the system matrix, as the mesh-width of the discretization tends to zero, accumulate on a set that is strictly larger than the convex hull of the essential spectrum of the integral operator.

In the paper [START_REF] Rahola | On the eigenvalues of the volume integral operator of electromagnetic scattering[END_REF], motivated by the convergence analysis of iterative solutions of the resulting large linear systems, the essential spectrum of the Maxwell volume integral operator was studied for the case of scattering by a dielectric ball in R 3 . This is a subset of the segment in the complex plane that corresponds to the essential numerical range of the integral operator. It is now known (see [START_REF] Costabel | The essential spectrum of the volume integral operator in electromagnetic scattering by a homogeneous body[END_REF][START_REF] Costabel | Volume integral equations for electromagnetic scattering in two dimensions[END_REF]) that the same form of the essential spectrum is valid for more general bounded Lipschitz domains. In [START_REF] Rahola | On the eigenvalues of the volume integral operator of electromagnetic scattering[END_REF], results of some numerical experiments are then shown that seem to indicate that the eigenvalues of the system matrices accumulate either at isolated points, corresponding to eigenvalues of the integral operator and hence to eigenvalues or resonances of the scattering problem, or at the points of the segment that is spanned by the essential spectrum of the integral operator. Looking closer at Figures 4.2-4.4 of [START_REF] Rahola | On the eigenvalues of the volume integral operator of electromagnetic scattering[END_REF], one can detect an "overshoot", namely that the observed segment of accumulation points is actually larger than the span of the essential spectrum.

In the paper [START_REF] Yurkin | Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived[END_REF], there is a discussion of the spectrum of the system matrices of the DDA scheme for the quasi-static Maxwell equations, motivated by the numerical modeling of the scattering of light by dust particles whose size is small with respect to the wavelength of the light ("Rayleigh particles"). Based on extensive experience with numerical computations using the DDA code ADDA, the authors are convinced that the DDA provides a faithful approximation of the solution of the volume integral equation in the sense that, among other things, the spectral measure of the DDA system matrices converges to the spectral measure of the volume integral operator. They provide plots of the spectral density of these matrices, including a zoom on a neighborhood of the lower end of the spectrum, see graph (a) in [START_REF] Yurkin | Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived[END_REF]FIG. 8]. There one can clearly see that there is an overshoot, namely a part of the spectrum below zero, and that its negative minimum does not disappear as the number of dipoles grows, but rather seems to converge to some number around -0.09. In a subsequent paper [START_REF] Smunev | Rectangular dipoles in the discrete dipole approximation[END_REF], the authors detect this "spill-out" of the spectrum of the DDA system matrices and relate it to an explosion of the needed iterations in an iterative solution method that they observed for large refractive indices. They study the behavior of this overshoot for anisotropic meshes, where it becomes larger, and for some recently introduced improvements of the DDA, where it seems to disappear.

For the quasi-static Maxwell case we prove below (see Proposition 3.15 and (3.51)) that for the classical DDA on a cubic grid such an overshoot indeed exists and that it amounts to an almost 20% increase of the length of the segment spanned by the essential spectrum. This somewhat unexpected result implies that the simple discretization scheme of the DDA does not provide a spectrally correct approximation of the strongly singular volume integral operator. The additional observation, supported by numerical experience, that this concerns only a small neighborhood of the essential spectrum or perhaps even only of the endpoints of this spectrum, whereas discrete eigenvalues and large parts of the spectral density nevertheless are correctly approximated, still awaits a precise description and proof.

It also implies that the DDA scheme is actually unstable in high-contrast situations, namely if the relative permittivity is very small (smaller than ∼ 0.093 ) or very large (larger than ∼ 11.8). We prove this here for the zero-frequency limit, but expect that it is also true for non-zero frequencies.

The Discrete Dipole Approximation.

As its name indicates, the DDA (sometimes called Coupled Dipole Approximation) can be considered as an approximation of a dielectric continuum described by Maxwell's equations by a different physical system consisting of a finite number of dipoles that are characterized by their polarizability, interacting via electromagnetic fields. The same mathematical system can be obtained by a procedure more amenable to arguments of numerical analysis, namely by transforming the Maxwell equations for the original dielectric continuum into an equivalent Lippmann-Schwinger style volume integral equation and then discretizing this integral equation by a simple delta-delta approximation on a regular grid

{x n | n ∈ Z d } of meshwidth h > 0.
Thus a linear integral equation on a bounded domain

Ω ⊂ R d (1.1) λu(x) - Ω K(x, y)u(y)dy = f (x) (x ∈ Ω)
will be approximated by the finite dimensional linear system

(1.2) λu m - xn∈Ω,n =m h d K(x m , x n )u n = f (x m ) (x m ∈ Ω) .
We omit the diagonal term m = n, because we shall have to do with singular kernels. Apart from this, (1.2) looks like a Galerkin method with Dirac deltas as trial and test functions.

Let us briefly describe the construction of the volume integral equation. A more detailed derivation can be found in [START_REF] Kirsch | An integral equation approach and the interior transmission problem for Maxwell's equations[END_REF] and, with special emphasis on the two-dimensional situation, in [START_REF] Costabel | Volume integral equations for electromagnetic scattering in two dimensions[END_REF]. We write the time harmonic Maxwell equations with normalized frequency κ ∈ C as a second order system for the electric field u.

(1.3) curl curl u -κ 2 u = iκJ.
Here it is assumed that the magnetic permeability is constant (normalized to 1) in the whole space. If one further assumes that the permittivity is equal to 1 outside of a bounded domain and the source current J has compact support, one can write this as a perturbation of the free-space situation

(1.4) curl curl u -κ 2 u = -κ 2 (1 -)u + iκJ.
Here the right hand side has compact support, and therefore convolution with the outgoing fundamental solution g κ of the Helmholtz equation and application of the operator ∇ div +κ 2 leads to the volume integral equation in distributional form

(1.5) u = -(∇ div +κ 2 )g κ (1 -)u + u inc .
Here the incoming field u inc combines the field generated by the current density with possible sourceless full space solutions of Maxwell's equations (plane waves etc.)

Equation (1.5) can be considered in any dimension d ≥ 2, but only d = 2 and d = 3 are relevant for electrodynamics. The equation can be written in the form of a second kind strongly singular integral equation with the d × d matrix valued kernel

(1.6) K(x, y) = -(D 2 + κ 2 )g κ (x -y) .
The integral operator thus defined involves second order distributional derivatives of the weakly singular kernel g κ (x -y). Instead of this form of an integro-differential operator, one can write the strongly singular integral operator also in the form of a Cauchy principal value integral, using the well-known relation (for more details, see section 3.5.1 below)

(1.7) D 2 R d g κ (x -y)u(y)dy = p.v. R d D 2 g κ (x -y)u(y)dy - 1 d u(x) .
If we further assume that the permittivity equals a constant r ∈ C \ {1} in Ω, we can divide by 1r and arrive at the final form (1.1) with the integral understood in the principal value sense, the kernel given by (1.6), and the spectral parameter λ defined by the relation

(1.8) λ = 1 1 -r - 1 d = d -1 + r d(1 -r ) . For d = 3, this relation λ = 2+ r 3(1-r )
is known in the DDA literature as Clausius-Mossotti polarizability, referring to the fact that 1 λ corresponds to the polarizability of the dipoles and to the Clausius-Mossotti equation between the molecular polarizability and the electric permittivity in a dielectric material, see for example [START_REF] Jackson | Classical Electrodynamics[END_REF]Section 4.5].

The principal part of the volume integral operator is obtained by taking the limit κ → 0, and we will refer to this situation as the quasi-static Maxwell case. The resulting kernel is homogeneous of degree -d, and this property allows to analyze the corresponding linear system (1.2) using Fourier analysis of Toeplitz matrices. For this reason we study in this paper a class of strongly singular kernels that includes the quasi-static Maxwell kernel.

Outline of the paper.

In Section 1.5 we define a class of strongly singular kernels that are homogeneous of degree -d and translation invariant, and we evoke the relation between the numerical range of the corresponding singular integral operator in L 2 and values of its symbol. The notion of numerical range allows to use a Lax-Milgram type argument to get a resolvent estimate for the restriction of the convolution operator to a bounded domain Ω.

After introducing in Section 1.6 the delta-delta discretization, we state in Theorem 1.6 the main stability result valid for our class of operators.

In Section 2 we study tools for proving stability results, namely infinite Toeplitz matrices and their symbols defined by Fourier series. Here a main difficulty is that one needs precise bounds for the values of a function (numerical symbol) defined by a Fourier series that is not absolutely convergent. We find that one can use Ewald's method for this purpose. The result is that the symbol of the Toeplitz matrix is a bounded function, and that its range is always a superset of the range of the symbol of the integral operator, but that it might be strictly larger. If this is the case, then stability of the delta-delta scheme implies well-posedness of the integral equation, but not vice versa: The numerical scheme does then not provide a spectrally correct approximation, and it might be unstable for values of the spectral parameter λ for which any Galerkin scheme of the integral equation would be stable.

In Section 3 we study in detail five representative examples.

Example 1 concerns the one-dimensional singular integral equation defined by the finite Hilbert transformation. Here the numerical symbol has a simple explicit expression, and this can be used to get estimates for the resolvent of the discretized operator by the resolvent of the integral operator, with constant equal to 1. This gives Theorem 3.4, which is the ideal stability result that subsequent results are measured against.

Examples 2 and 3 exhibit different behavior of the delta-delta scheme for two strongly singular integral operators in two dimensions. Whereas the two integral operators are equivalent, related by a simple rotation of the coordinate system, the two discrete systems show opposite behavior: We prove that in Example 2 the ranges of the symbol of the integral operator and of the symbol of the infinite Toeplitz matrix are identical, whereas in Example 3 there is an overshoot; the region of instability of the approximation scheme is strictly larger than the numerical range of the integral operator.

In Example 4, we graphically illustrate the relations, proved in Sections 1.5 and 2, between the spectrum and numerical range of the system matrices and the numerical range of the singular integral operator by considering a non-selfadjoint case. The kernel is a complex-valued function whose real and imaginary parts are given by the kernels of Examples 3 and 2, respectively. The kernels studied in Examples 2 and 3 are also the off-diagonal and diagonal terms, respectively, in the matrix-valued kernel of the quasi-static Maxwell volume integral operator, which is the subject of Example 5. We study this for dimensions d ≥ 2 and give more precise results for d = 2 and d = 3. In two dimensions we find the same overshoot of the numerical range of the numerical symbol versus the symbol of the integral operator as Example 3. In three dimensions this overshoot is even larger, and it can be verified numerically either by computing the numerical symbol using Ewald's method or by studying the asymptotic behavior of the smallest and largest eigenvalues of the matrix of the linear system as the mesh width tends to zero. 1.4. Notation for Fourier transforms and Fourier series.

We use the following convention for the Fourier transformation in R d .

(1.9)

f (ξ) = F f (ξ) = R d f (x)e iξ•x dx .
Inverse:

(1.10)

f (x) = F -1 f (x) = (2π) -d R d f (ξ)e -ix•ξ dξ .
For Fourier series, we use the following notation. For a sequence a : Z d → C, its Fourier series is defined as

(1.11) ã(τ ) = m∈Z d a(m)e im•τ , τ ∈ Q = [-π, π] d .
Inverse:

(1.12)

a(m) = (2π) -d Q ã(τ )e -im•τ dτ.
The definitions are extended in the usual way from convergent sums and integrals to suitable spaces of functions and distributions. In particular, we have Parseval's theorem

(1.13) f → (2π) -d 2 f : L 2 (R d ) → L 2 (R d ) and a → (2π) -d 2 ã : 2 (Z d ) → L 2 (Q)
are unitary (i.e. isometric Hilbert space isomorphisms).

Combining the Parseval formula and the convolution theorem gives

(1.14) R 2d u(x) k(x -y) v(y) dy dx = (2π) -d R d u(ξ) k(ξ) v(ξ) dξ , (1.15)
m.n∈Z d a(m) c(m -n) b(n) = (2π) -d Q ã(τ ) c(τ ) b(τ ) dτ .
From these formulas follows immediately that the operators of convolution with k in L 2 (R d ) and of discrete convolution with c in 2 (Z d ) are bounded if and only if the "symbols" k and c are bounded functions belonging to L ∞ (R d ) and L ∞ (Q), respectively.

Sufficient conditions for this are that k ∈ L 1 (R d ) and c ∈ 1 (Z d ). But these conditions are not necessary, and it is precisely the situation where they are not satisfied that will be relevant in the following.

We will use the Poisson summation formula in the form (1.16)

m∈Z d f (m)e im•τ = n∈Z d f (τ + 2πn) .
A sufficient (but in no way necessary) condition for (1.16) to hold for all τ is that

f Z d ∈ 1 (Z d ) and f ∈ L 1 (R d ).
If we do not assume f Z d ∈ 1 , but only f ∈ L 1 , then f is bounded, the left hand side of (1.16) converges in the distributional sense and the right hand side converges in L 1 (Q). Then (1.16) is true in a weaker sense, the distributional left hand side being equal to the L 1 (Q) right hand side.

Example: Gaussian with parameter s > 0.

(1.17)

f (x) = e -|x| 2 s ⇐⇒ f (ξ) = ( π s ) d 2 e -|ξ| 2 4s .
For this example the Poisson summation formula takes the form (for

τ ∈ R d ) (1.18) m∈Z d e -|m| 2 s e im•τ = n∈Z d ( π s ) d 2 e -|τ +2πn| 2 4s .
1.5. Kernels and their symbols.

1.5.1. Homogeneous kernels. Later on, we will consider a rather restricted class of strongly singular integral operators on R d that are convolutions with kernel functions of the form

(1.19) K(x) = p(x) |x| -d-2
where p is a homogeneous polynomial of degree 2 .

But first we recall some well-known general properties of homogeneous functions and distributions that can be found, for example, in [8, Chap. III].

Let K be a function on R d , positively homogeneous of degree -d and smooth outside of the origin. For a given > 0, one can define a distribution K ∈ S (R d ) that coincides with K on R d \ {0} by its action on a test function φ as

(1.20) K , φ = |x|< K(x)(φ(x) -φ(0))dx + |x|> K(x)φ(x) dx .
This is independent of if and only if K satisfies the cancellation condition on the unit sphere S d-1

(1.21)

S d-1 K ds = 0 .
In this case, we denote the distribution simply by K, and we can take the limit → 0, thus we get the Cauchy principal value.

(

1.22) K, φ = p.v. K(x)φ(x) dx = lim →0 |x|> K(x)φ(x) dx .
Another consequence of the cancellation condition (1.21) is that the Fourier transform K of the homogeneous distribution K is a bounded function homogeneous of degree 0, smooth outside of the origin and also satisfying the cancellation condition. The operator A of convolution with K is therefore bounded in L 2 (R d ). Note that in the absence of condition (1.21), K would have a logarithmic singularity at 0.

The operator A is diagonalized by Fourier transformation:

(1.23)

F Au = K u .
Therefore in L 2 (R d ), we can obtain information about the spectrum Sp(A) and about the numerical range W (A) from the corresponding easily checked information about the operator of multiplication by the symbol K.

We recall that the numerical range of A is defined by

W (A) = {(u, Au) | u = 1} ,
where (•, •) denotes the Hilbert space inner product. It is convex by the Toeplitz-Hausdorff theorem and it contains the spectrum of A. Denote by im( K) = K(R d ) the image (range) of K. This is a compact set. We note a first result implied by the unitary equivalence (1.23) with the multiplication operator.

Lemma 1.1. The spectrum Sp(A) is the image im( K), and the closure W (A) of the numerical range of A is the convex hull of im( K).

It is well known (and easy to prove) that the numerical range allows estimates for the operator norm of the resolvent: For any

λ ∈ C \ W (A), (λI -A) -1 ≤ dist(λ, W (A)) -1 .
It is also monotone with respect to inclusions of subspaces, a property not shared by the spectrum. Given an open set Ω ⊂ R d , we denote by A Ω the restriction of the convolution operator A to L 2 (Ω) and consider the strongly singular integral equation (λI

-A Ω )u = f , or in detail (1.24) λu(x) -p.v. Ω K(x -y)u(y) dy = f (x) (x ∈ Ω).
From the definition of the numerical range follows immediately the inclusion

W (A Ω ) ⊂ W (A).
We can summarize this discussion:

Proposition 1.2. Let C ⊂ C be a closed convex set such that K(ξ) ∈ C for all ξ ∈ S d-1
. Then for all λ ∈ C and any f ∈ L 2 (Ω), the integral equation (1.24) has a unique solution u ∈ L 2 (Ω), and there is a resolvent estimate in the L 2 (Ω) norm

(1.25) (λI -A Ω ) -1 ≤ dist(λ, C) -1 .
Remark 1.3. The same argument implies stability for any Galerkin method: Let X h be any closed subspace of L 2 (Ω), and let A h : X h → X h be the operator defined by restricting the sesquilinear form (u, Av) to X h × X h . Then the statement of Proposition 1.2 remains true if we replace A Ω by A h .

Remark 1.4. Whereas there is, in general, no simple relation between the spectra Sp(A Ω ) and Sp(A), for the numerical ranges of our convolution operators with homogeneous kernels we not only have the inclusion W (A Ω ) ⊂ W (A), but also the converse. Namely there holds

(1.26) W (A Ω ) = W (A) for any open subset Ω ⊂ C.
Proof. The set of Rayleigh quotients (u,Au) (u,u) , where u ∈ L 2 (R d ) \ {0} has compact support, is a dense subset of W (A). We show that it is a subset of W (A Ω ): Indeed, let u be such a function and let ρ > 0 and a ∈ R d be chosen such that the support of the function u ρ,a defined by u ρ,a (x

) = u(ρx + a) is contained in Ω. Then (u, Au) (u, u) = (u ρ,a , Au ρ,a ) (u ρ,a , u ρ,a ) ∈ W (A Ω ).
1.5.2. Special kernels. For d = 1, there is essentially only one non-trivial kernel homogeneous of degree -d, namely K(x) = 1 x . In R d for d ≥ 2, while some of the following analysis would be possible for more general homogeneous kernels, we focus now on the situation (1.19). This means that from now on, we fix a strongly singular kernel K and a homogeneous polynomial p of degree 2 with K(x) = p(x)|x| -d-2 , satisfying (1.21), considered as a distribution on R d according to (1.22), and we denote by K its Fourier transform. Lemma 1.5. Let K have the form (1.19) and satisfy (1.21). Then

(1.27) K(ξ) = -ν d p(ξ) |ξ| 2 , where ν d = 2π d 2 d Γ( d 2 )
is the volume of the unit ball in R d .

Proof. We first compute the Fourier transform of p(x)e -|x| 2 s , using (1.17)

F x →ξ [p(x)e -|x| 2 s ] = ( π s ) d 2 p(-i∂ ξ )e -|ξ| 2 4s .
The evaluation of these derivatives leads to the following simple result, as we will show:

(1.28)

F x →ξ [p(x)e -|x| 2 s ] = -( π s ) d 2 1 4s 2 p(ξ)e -|ξ| 2 4s .
For j, k ∈ {1, . . . , d} with j = k, let

(1.29) a jk (x) = x 2 j -x 2 k , b jk (x) = x j x k .
Any homogeneous polynomial of degree 2 satisfying the cancellation condition S d-1 p = 0 is a linear combination of the a jk and b jk , so we need to verify (1.28) only for these.

Note that ∂ ξ j e -|ξ| 2 4s = -1 2s ξ j e -|ξ| 2 4s and ∂ 2 ξ j e -|ξ| 2 4s = -1 2s + 1 4s 2 ξ 2 j e -|ξ| 2 4s .
Then for p = a jk , we see

(∂ 2 ξ k -∂ 2 ξ j )e -|ξ| 2 4s = 1 4s 2 (ξ 2 k -ξ 2 j )e -|ξ| 2 4s ,
and for p = b jk , we see

∂ ξ j ∂ ξ k e -|ξ| 2 4s = 1 4s 2 ξ k ξ j e -|ξ| 2 4s .
Thus in both cases we have

(1.30) p(∂ ξ )e -|ξ| 2 4s = 1 4s 2 p(ξ)e -|ξ| 2 4s ,
and (1.28) is proved.

Now we use the definition of the Gamma function

Γ(a) = ∞ 0 t a e -t dt t = |x| 2a ∞ 0
s a e -|x| 2 s ds s to write the kernel as an integral over Gaussians

(1.31) K(x) = p(x)|x| -d-2 = 1 Γ( d 2 + 1) ∞ 0 s d 2 +1 p(x)e -|x| 2 s ds s .
Taking Fourier transforms and using (1.28), we find with u

= |ξ| 2 4s (1.32) -K(ξ) = π d 2 4Γ( d 2 + 1) ∞ 0 s -1 p(ξ)e -|ξ| 2 4s ds s = p(ξ)|ξ| -2 π d 2 d 2 Γ( d 2 ) ∞ 0 e -u du = ν d p(ξ) |ξ| -2
as claimed.

1.6. Delta-delta discretization.

Let N ∈ N, fix some origin a N ∈ R d and define the cubic grid of meshwidth h = 1 N by

Σ N = {x N m = a N + m N | m ∈ Z d } .
We further define

ω N = {m ∈ Z d | x N m ∈ Ω} .
Then a very simple discretization of the strongly singular integral equation (λI

-A Ω )u = f (1.24) is the following (1.33) λu m -N -d n∈ω N ,m =n K(x N m -x N n )u n = f (x N m ) , (m ∈ ω N ) , or in shorthand (λI -T N )U = F .
The name "delta-delta discretization" points at the fact that this discretization formally looks like a Galerkin method for the integral equation (1.24) with Dirac deltas as both test and trial functions, except for the diagonal terms of T N , where we put zero, which is natural in view of the cancellation condition (1.21).

Our aim in this paper is to analyze the linear system (1.33), in particular its stability in 2 (ω N ), in the same way as we did above for the integral equation (1.24) in L 2 (Ω), and to compare the two.

We state a general result here, which we prove in the next section. More precise results will be given below in Section 3 for some examples, in particular those mentioned in Subsection 1.1.

Theorem 1.6. Let K be a strongly singular kernel satisfying (1. [START_REF] Yurkin | Convergence of the discrete dipole approximation. I. Theoretical analysis[END_REF]) and (1.21). Then there exists a compact convex set C ⊂ C such that for any λ ∈ C \ C, any N ∈ N for which ω N is non-empty, and for any F ∈ 2 (Ω N ), the system (1.33) has a unique solution, and there is a uniform estimate for the inverse in the 2 (Ω N ) operator norm

(1.34) (λI -T N ) -1 L( 2 (Ω N )) ≤ dist(λ, C) -1
. Furthermore, with the strongly singular integral operator A defined above in Section 1.5.1, there holds the inclusion

(1.35) W (A) ⊂ C .
Remark 1.7. Note that the inclusion W (A) ⊂ C implies that for λ ∈ C the singular integral equation is uniquely solvable, too, and provides the a priori estimate (1.25), for any domain Ω ⊂ R d . On the other hand, in order to guarantee stability for λ ∈ C \ C, the inclusion may need to be strict, as we shall see in the examples, and then there may be λ ∈ C \ W (A) for which the singular integral equation is well posed, but the delta-delta discretization is unstable.

THE DISCRETE SYSTEM

Let T N be the matrix representing the discretized integral operator in (1.33):

(2.1)

T N = (t N mn ) m,n∈ω N with t N mn = N -d K(x N m -x N n ) (m = n) 0 (m = n) .
Our aim is to bound the numerical range W (T N ) independently of N .

2.1. Toeplitz structure.

The matrix elements t N mn of T N do not depend on the choice of the origin a N , and since we assumed that K is homogeneous of degree -d, we have

N -d K(x N m -x N n ) = K(m -n) , hence T N is a finite section of a fixed infinite Toeplitz (discrete convolution) matrix (2.2) T = (t mn ) m,n∈Z d with t mn = K(m -n) (m = n) , 0 (m = n) .
Theorem 1.6 will be proved if we can show that T defines a bounded linear operator in 2 (Z d ) whose numerical range W (T ) contains W (A). We can then choose C as the closure of W (T ).

We use Fourier series and the convolution theorem to diagonalize the matrix T and to represent the sesquilinear form defined by the matrix T N , compare (1.15). For U = (u m ) m∈ω N , we find

(2.3) (U, T N U ) = (2π) -d Q F (τ ) |ũ(τ )| 2 dτ .
Here

ũ(τ ) = m∈ω N u m e im•τ and Q = [-π, π] d . F (τ )
is the symbol (characteristic function) of the Toeplitz matrix T :

(2.4)

F (τ ) = m∈Z d ,m =0 K(m) e im•τ .
The problem is now reduced to the study of the operator of multiplication by the function F in L 2 (Q).

Lemma 2.1. The operator T :

2 (Z d ) → 2 (Z d ) is bounded if and only if F ∈ L ∞ (Q).
The closure of W (T ) is the closed convex hull of the range im(F ) = {F (τ ) | τ ∈ Q} and is also equal to the closure of the union N ∈N W (T N ) .

Proof: This is immediate from (2.3).

Ewald method.

The problem that makes the statement F ∈ L ∞ (Q) non trivial is that the Fourier series (2.4) is not absolutely convergent. The sequence (K(m)) m∈Z d is of order O(|m| -d ) at infinity and therefore in p (Z d ) for all p > 1, but not for p = 1. Its membership in 2 (Z d ) implies, for example, that the series converges in the sense of L 2 (Q). The slow convergence of the Fourier series for F makes it also unsuitable for using it in numerical computations to find bounds for im(F ).

We will use a variant of a method introduced by P. P. Ewald [START_REF] Ewald | Die Berechnung optischer und elektrostatischer Gitterpotentiale[END_REF] in 1921 as a tool to compute slowly converging lattice sums. It has become a routine method for the computation of periodic and quasiperiodic Green functions, with application in numerical electrodynamics and other fields where periodic structures appear. Among the many presentations of the method: Appendix A of the article [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] or Section 2.13.3 in the book [START_REF] Ammari | Mathematical and computational methods in photonics and phononics[END_REF].

We use it here as a summation method for our slowly converging Fourier series. In our restricted setting it turns out to give surprisingly simple results.

The method introduces a decomposition K = K F + K P for the coefficients and correspondingly F = F F + F P for the Fourier series in such a way that both K F and the Fourier transform K P of K P are exponentially decreasing at infinity, so that both the Fourier series for F F (τ ) and the Poisson sum (compare (1.16)) for F P (τ ) are rapidly convergent, which not only proves the boundedness of F , but gives also a fast numerical algorithm for its computation.

In the literature one often labels the two terms in the decomposition "spatial" and "spectral" sums, but this is not pertinent to our situation, where the lattice sum runs over the Fourier variable, and the Fourier series runs overs spatial points. So we will use "Fourier" and "Poisson" sums as labels.

The idea of Ewald's method is to represent K(x) by an integral over Gaussians from 0 to ∞ as we did already in Section 1.5.2 above:

(2.5)

K(x) = p(x)|x| -d-2 = p(x) Γ( d 2 + 1) ∞ 0 s d 2 e -|x| 2 s ds
and then to split the integral at a point β 2 > 0:

K F (x) = p(x) Γ( d 2 + 1) ∞ β 2 s d 2 e -|x| 2 s ds , (2.6)
K P (x) = p(x) Γ( d 2 + 1) β 2 0 s d 2 e -|x| 2 s ds . (2.7)
We see that K F is simply the product of K by a function exponentially decreasing at infinity (2.8)

K F (x) = K(x) Γ( d 2 + 1, β 2 |x| 2 ) Γ( d 2 + 1) with the (upper) incomplete Gamma function (see [1, §6.5]) Γ(a, x) = ∞ x t a-1 e -t dt . Therefore K F (x) = O(|x| 2 e -β 2 |x| 2 )
as |x| → ∞ , and the Fourier series for F F (τ )

(2.9) F F (τ ) = m∈Z d ,m =0 K F (m) e im•τ
converges rapidly, implying that F F is an analytic function on R d /(2πZ) d .

Consequently, the Fourier series for F P (τ ) converges as slowly as the one for F (τ ), and we use instead the Poisson summation formula (1.16) and write (2.10)

F P (τ ) = n∈Z d K P (τ + 2πn) .
We can evaluate K P with the formulas used for K in Lemma 1.5. As in (1.32) we obtain (2.11)

K P (ξ) = -π d 2 4Γ( d 2 + 1) β 2 0 s -1 p(ξ)e -|ξ| 2 4s ds s = -p(ξ)|ξ| -2 π d 2 d 2 Γ( d 2 ) ∞ |ξ| 2 4β 2 e -u du = K(ξ)e - |ξ| 2
4β 2 .

Therefore we also obtain a very simple form for the Fourier transform, namely that K P is just the symbol of A cut off at infinity, and therefore the series (2.10) converges absolutely and uniformly. At most one term in the sum may be discontinuous, when τ + 2πn = 0, and for τ ∈ Q this is the term with n = 0. We can summarize the result.

Proposition 2.2. The symbol F (τ ) of the infinite Toeplitz matrix T is a bounded function given for any β > 0 by the exponentially convergent sums

(2.12)

F (τ ) = m∈Z d ,m =0 K(m) Γ( d 2 +1,β 2 |m| 2 ) Γ( d 2 +1) e im•τ + n∈Z d K(τ + 2πn) e - |τ +2πn| 2 4β 2
.

In the period cube

Q = [-π, π] d , it is C ∞ outside of 0,
and it has the form

(2.13) F (τ ) = K(τ ) + F 0 (τ ) where F 0 is analytic in Q and F 0 (0) = 0.
Proof. We have proved equation (2.12) above, except for one point: From Poisson's summation formula follows that the Poisson sum (2.10) equals the Fourier series with coefficients K P (m), m ∈ Z d , including m = 0. But in the Fourier series (2.4) defining F (t) as well as in (2.9) defining F F (t), we have excluded m = 0. So we should compensate for K P (0), which according to (2.7) equals

K P (0) = p(0)β d Γ( d 2 + 1)
. Now, since we assumed p(x) to be a homogeneous polynomial of degree 2, we have p(0) = 0 and hence no compensation is needed.

Representing F 0 as

F 0 (τ ) = m∈Z d ,m =0 K(m) Γ( d 2 +1,β 2 |m| 2 ) Γ( d 2 +1) e im•τ + n∈Z d ,n =0 K(τ + 2πn)e - |τ +2πn| 2 4β 2 + K(τ ) e - |τ | 2 4β 2 -1 ,
we see immediately that it is analytic. For finding F 0 (0), we can use the following observation.

Lemma 2.3. Let S ⊂ R d be a finite set that is cubically symmetric, i. e. invariant under reflections at coordinate planes and under permutations of the coordinates, and let p be a homogeneous polynomial of degree 2 satisfying the cancellation condition

S d-1 p = 0. Then x∈S p(x) = 0.
This is immediately clear when p is one of the a jk or b jk from (1.29), and it is therefore true for all p satisfying the (spherical) cancellation condition.

For any M ∈ R, the set {m ∈ Z d | |m| 2 = M } is either empty or cubically symmetric. Therefore for τ = 0, the two sums in the representation of F 0 (τ ) are 0. The last term

K(τ ) e - |τ | 2 4β 2 -1 = -ν d p(τ ) e - |τ | 2 4β 2 -1 |τ | 2
tends to ν D p(0)/(4β 2 ) = 0 as τ → 0, and hence F 0 (0) = 0.

Poposition 2.2 implies Theorem 1.6, where C is the closed convex hull of im(F ). The inclusion W (A) ⊂ C is easy to see from (2.13): Given > 0, let δ > 0 be such that for |τ | < δ we have |F 0 (τ )| < . Since K is homogeneous of degree zero, it takes all of its values already on the ball B δ (0) of radius δ. Thus

im( K) ⊂ F (B δ (0)) + B (0) ⊂ im(F ) + B (0) .
Taking convex hulls shows that

W (A) ⊂ C + B (0) for all > 0 .
Remark 2.4. The very simple form of the Ewald representation (2.12) comes from the very simple form of the Fourier transforms (1.27) and (1.28), which in turn rely on the cancellation condition (1.21). Now for the kernel K this condition is natural, because it is necessary in order to represent K as a homogeneous distribution and to have a bounded Fourier transform. But for the symbol K it is not as natural. We can add a constant and still have a function homogeneous of degree zero, which will then not satisfy the cancellation condition. An example is ξ j ξ k |ξ| -2 for all j, k, even for j = k.

On the other hand, the representation K(x) = p(x)|x| -d-2 may not be the most natural, one may come across cases (see Example 5 below) like

K jk (x) = δ jk |x| -d -d x j x k |x| -d-2 ,
where for j = k the two terms in the sum do not separately satisfy (1.21). This fits into our framework, however, because

K jk (x) = -d b jk (x)|x| -d-2 for j = k, and K kk (x) = d j=1 a jk (x)|x| -d-2 .
If one treats the two terms individually, one may get formulas for Fourier transforms and for the Ewald splitting that are less symmetric than what we presented above.

An integral representation.

We have another look at the Ewald splitting for the numerical symbol F (ξ) = F F (ξ) + F P (ξ) described in (2.6)-(2.12)

F F (ξ) = m∈Z d p(m) Γ( d 2 + 1) ∞ β 2 s d 2 e -|m| 2
s ds e im•ξ (2.14)

F P (ξ) = n∈Z d -π d 2 4Γ( d 2 + 1) β 2 0 s -2 p(ξ + 2πn)e -|ξ+2πn| 2 4s
ds . (2.15) These formulas are valid for any 0 < β < ∞. All the sums and integrals are converging absolutely here, and therefore we can interchange sums and integrals.

F F (ξ) = ∞ β 2 H F (ξ, s) ds with H F (ξ, s) = m∈Z d p(m) Γ( d 2 + 1) s d 2 e -|m| 2 s e im•ξ (2.16) F P (ξ) = β 2 0 H P (ξ, s) ds with H P (ξ, s) = n∈Z d -π d 2 p(ξ + 2πn) 4Γ( d 2 + 1) s -2 e -|ξ+2πn| 2 4s . (2.17)
From the definition (2.16) of H F and the fact that |m| ≥ 1 in the sum follows without difficulty that for any 0 < γ < 1 there exists a constant C such that

(2.18) |H F (ξ, s)| ≤ C e -γs for all s ≥ 1, ξ ∈ R d .
To see the behavior of H P (ξ, s) from (2.17), we decompose It remains to analyze the term with n = 0, i.e. H 1 . It is clear that it vanishes for ξ = 0, and for every ξ = 0 there exists a constant C ξ and 0 < γ < |ξ| 2 4 such that

H P (ξ, s) = H 0 (ξ, s) + H 1 (ξ, s) with H 0 (ξ, s) = -π d 2 4Γ( d 2 +1) n∈Z d ,n =0 p(ξ + 2πn)s -2 e -|ξ+2πn| 2 4s , (2.19) H 1 (ξ, s) = -π d 2 4Γ( d 2 +1) p(ξ)s -2 e -
(2.22) |H 1 (ξ, s)| ≤ C ξ min{s -2 , e -γ s } for all s ∈ (0, ∞) .
Thus H 1 (ξ, s) is integrable over s ∈ (0, ∞) for all ξ, but there is no uniform bound for C ξ :

Considering sup s>0 s -2 e -|ξ|
F 0 (ξ) = ∞ 0 H 0 (ξ, s) ds and K(ξ) = ∞ 0 H 1 (ξ, s) ds .
In these integrals, the functions s → H 0 (ξ, s), s → H 1 (ξ, s), and s → H(ξ, s) are integrable on (0, ∞) for any ξ ∈ Q, for any ξ ∈ R d \ {0}, and for any ξ ∈ Q \ {0}, respectively.

The integral representations (2.23) and ( 2.24) will be used below to get bounds for the function F (ξ) from estimates for H(ξ, s). The latter will be a consequence of the following observation that can be proved using Fourier representations (2.16) for H and (1.28) for H 1 .

Lemma 2.6. The functions

(ξ, s) → s -d 2 H 0 (ξ, s), (ξ, s) → s -d 2 H 1 (ξ, s), (ξ, s) → s -d 2 H(ξ, s
), are solutions of the heat equation

(∂ s -∆ ξ )u(ξ, s) = 0 in Q × (0, ∞).

Matrix-valued kernels.

Until now, we have considered kernel functions with values in C and integral operators acting on scalar functions. The generalization to vector-valued functions and matrix-valued kernels is simple and straightforward, and we do not find it necessary to introduce typographic distinctions for the vector-valued objects. The main difference is that in the general theory of Section 1.5, one has to use the numerical range W (K(x)) of the matrix K(x) instead of the value K(x) in statements such as Lemma 1.1 and Proposition 1.2. In particular (2.25) W (A) is the closed convex hull of ξ∈R d W ( K(ξ)) . Theorem 1.6 remains literally true, but for the construction of the set C one has once again to use the numerical range W (F (t)) of the matrix-valued function F . In Lemma 2.1, the characterization of the numerical range W (T ) is to be understood as follows.

Lemma 2.7. The closure of W (T ) is the closed convex hull of τ ∈Q W (F (τ )) and is also equal to the closure of the union N ∈N W (T N ).

The basic Parseval-convolution formula (1.15) now has to be written, instead of the scalar version (2.3), as

(2.26) (U, T N U ) = (U, T U ) = (2π) -d Q ũ(τ ) F (τ ) ũ(τ ) dτ .
Here ũ(τ ) = m∈ω N u m e im•τ , and F (τ ) is the matrix-valued symbol of the block Toeplitz matrix

T = K(m -n) m,n∈Z d : (2.27) F (τ ) = m∈Z d ,m =0 K(m) e im•τ .
From (2.26) one can immediately read the properties of the numerical range stated in Lemma 2.7.

In this paper, most considered examples of kernels are real-valued and the matrices symmetric, in which case the integral operators are selfadjoint, and the numerical ranges consist of intervals in the real line. We start with the simplest example of a strongly singular integral equation and show that the stability of its delta-delta approximation can be completely analyzed, resulting in a kind of ideal stability theorem.

-A Ω )u = f , (3.1) λu(x) - 1 iπ Ω u(y) x -y dy = f (x) , x ∈ Ω .
The integral is understood in the Cauchy principal value sense. The kernel function K(x) = 1 iπx has the Fourier transform K(ξ) = sign ξ.

The operator A of convolution with K on R is the Hilbert transformation. It satisfies A 2 = I, and its spectrum (in a large class of function spaces, for instance L p (R) with 1 < p < ∞) is {-1, 1}, consisting of two eigenvalues of infinite multiplicity.

The finite Hilbert transformation A Ω and its spectral theory are also well studied classical objects, see for example [START_REF] Koppelman | Spectral representations for finite Hilbert transformations[END_REF]. Here the spectrum depends on the function space; for L p (Ω) it is strongly dependent on p, but not on Ω, as long as Ω is a proper subinterval of R. For p = 2 one has the following description.

(3.2) (λI -A Ω ) -1 = dist(λ, C) -1 .
Explicit formulas for the resolvent are known. For the infinite Hilbert transformation this is trivially obtained by algebra:

(λI -A) -1 = 1 λ 2 -1 (λI + A) ,
and for the finite Hilbert transformation, formulas for the resolvent can be found for example in [START_REF] Öhngen | Zur Theorie der endlichen Hilbert-Transformation[END_REF] or [START_REF] Tricomi | Integral equations[END_REF].

3.1.2. The discrete system. We use the notation of Section 1.

2 with d = 1, in particular x N m = a N + m N and ω N = {m ∈ Z | x N m ∈ Ω}.
The simple delta-delta discretization of our singular integral equation (3.1) is

(3.3) λu m - 1 iπN n∈ω N ,m =n u n x N m -x N n = f (x N m ) , (m ∈ ω N ) .
The system matrix T N with matrix elements

1 iπN 1 x N m -x N n (m, n ∈ ω N ) is a finite section of the infinite Toeplitz matrix T = 1 iπ(m -n) m,
n∈Z with zero on the diagonal.

The symbol F (τ ) is now given by the Fourier series

(3.4) F (τ ) = m∈Z,m =0 e imτ iπm = ∞ m=1 2 sin mτ πm , τ ∈ Q = [-π, π] .
This series converges for all t ∈ Q to the well known saw-tooth function

(3.5) F (τ ) = sign τ - τ π (τ = 0) , F (0) = 0 .
The range of this function is the interval (-1, 1).

Properties of the matrix T follow immediately from this symbol F and can be summarized as follows. For λ ∈ C \ C, there is a uniform resolvent estimate in the 2 operator norm

(3.7) (λI -T N ) -1 ≤ dist(λ, C) -1 .
The converse is also true: If there is a uniform stability estimate For such λ, there is an estimate for the operator norms

(λI -T N ) -1 ≤ C for all N ,
(3.8) (λI -T N ) -1 L( 2 (ω N )) ≤ (λI -A Ω ) -1 L(L 2 (Ω)) . 3.2. Example 2. Dimension d = 2, kernel x 1 x 2 |x| -4 .
We consider now the simplest higher-dimensional example where in the notation of Section 1.5.2 d = 2 and p(x) = -1 π b 12 (x), see (1.29). We show that the stability of its delta-delta approximation follows a similar simple pattern as in the previous one-dimensional example, although the proof is non-trivial.

3.2.1.

The singular integral equation. The kernel and its Fourier transform are (3.9)

K(x) = - x 1 x 2 π|x| 4 , K(ξ) = ξ 1 ξ 2 |ξ| 2 .
For Ω ⊂ R 2 , we consider the singular integral equation (λI

-A Ω )u = f as in (1.24) (3.10) λu(x) -p.v. Ω K(x -y)u(y) dy = f (x).
Observing that the range of the function K is the interval [-1 2 , 1 2 ], we can formulate the result of Proposition 1.2 as follows Lemma 3.5.

Let C = [-1 2 , 1 2 ].
For Ω = R 2 , both the spectrum Sp(A Ω ) and the closure of the numerical range W (A Ω ) in L 2 (Ω) are equal to C. For any open subset Ω ⊂ R 2 , the closure of the numerical range in L 2 (Ω) satisfies W (A Ω ) ⊂ C, and there is a resolvent estimate in the L 2 (Ω) operator norm

(3.11) (λI -A Ω ) -1 ≤ dist(λ, C) -1 .
3.2.2. The discrete system. Let now Ω be a bounded domain in R 2 . In the notation of Section 1.6 with d = 2, the regular grid consists of the points

x N m = a N + m N , indexed by ω N = {m ∈ Z 2 | x N m ∈ Ω}. The simple delta-delta discretization of our singular integral equation (3.10) is (3.12) λu m + 1 π N -2 n∈ω N ,m =n (x N m,1 -x N n,1 )(x N m,2 -x N n,2 ) |x N m -x N n | 4 u n = f (x N m ) , (m ∈ ω N ) .
The system matrix T N is now a finite section of the infinite Toeplitz matrix

T = - 1 π (m 1 -n 1 )(m 2 -n 2 ) |m -n| 4 m,n∈Z 2
with zero on the diagonal.

Its symbol is therefore given by the Fourier series for τ

∈ Q = [-π, π] 2 (3.13) F (τ ) = - m∈Z 2 ,m =0 m 1 m 2 π|m| 4 e im•τ = 4 π ∞ m 1 ,m 2 =1 m 1 m 2 (m 2 1 + m 2 2 ) 2 sin(m 1 τ 1 ) sin(m 2 τ 2 ) .
Whereas we do not know an explicit closed form expression for this function, we know from the results of Section 2.2 using Ewald's method that it is bounded and that it can be written as in equation (2.13)

(3.14) F (τ ) = K(τ ) + F 0 (τ ) where F 0 is analytic in Q and F 0 (0) = 0.
In addition, we know from (3.13) that F vanishes on the boundary of Q, hence

(3.15) F 0 (τ ) = -K(τ ) (τ ∈ ∂Q) .
In the previous example, we used the explicit expression of F (τ ) for finding the range of F . In fact, the function F 0 in that case was just the linear interpolation between the two values of the symbol K on ∂Q, which implied that the closed convex hull of im(F ) was the same as the convex hull of im( K). In the present case, we do not have a simple formula, but we can still prove that the conclusion is true.

Lemma 3.6. Let F (τ ) be as defined in (3.13). Then for any τ ∈ Q

(3.16) F (τ ) ∈ C = [- 1 2 , 1 2 ].
The proof is not obvious, although the claim is numerically evident if we compute F using Ewald's method and plot its graph, see the contour plot in Figure 1. Before we give the proof, let us draw the conclusion for the stability of the numerical scheme (3.12).

Corollary 3.7.

Let C = [-1 2 , 1 2 ] and λ ∈ C \ C.
Then for any N the linear system (3.12) has a unique solution, and there is a uniform resolvent estimate

(3.17) (λI -T N ) -1 L( 2 (Ω N )) ≤ dist(λ, C) -1 .
Proof. For symmetry reasons, it is sufficient to prove (3.16) for τ ∈ Q ++ = (0, π) 2 . For the proof of Lemma 3.6, we will show the following:

(3.18) For any ξ ∈ Q ++ , F (ξ) ≥ 0 and F 0 (ξ) ≤ 0 .

This implies 0 ≤ F (ξ) ≤ K(ξ) ≤ 1 2 , hence (3.16) . We use the integral representations from Proposition 2.5

(3.19) F 0 (ξ) = ∞ 0 H 0 (ξ, s)ds , K(ξ) = ∞ 0 H 1 (ξ, s) ds , F (ξ) = ∞ 0 H(ξ, s) ds . According to (2.19), H 0 (ξ, s) = H(ξ, s) -H 1 (ξ, s) with (3.20) H 1 (ξ, s) = ξ 1 ξ 2 4s 2 e -|ξ| 2 4s , H(ξ, s) = n∈Z 2 H 1 (ξ + 2πn, s) .
Let 0 ≤ < T and Σ T = Q ++ × ( , T ). In Σ T , we want to use the maximum principle for the heat equation (see Lemma 2.6) for the functions H0 (ξ, s) = s -1 H 0 (ξ, s) and H(ξ, s) = s -1 H(ξ, s).

Since H1 (ξ, s) = s -1 H 1 (ξ, s) is continuous for (ξ, s) ∈ R 2 × [0, ∞) \ {0, 0} and the Poisson series H0 (ξ, s) = n∈Z 2 ,n =0 H1 (ξ + 2πn, s)
converges uniformly for (ξ, s) ∈ Σ T 0 for all T > 0, we see that H0 is continuous in Σ T 0 with initial value H0 (ξ, 0) = 0. On the lateral boundary we use the Fourier representation (see (2.16))

H(ξ, s) = 4 π ∞ m 1 ,m 2 =1
m 1 m 2 e -|m| 2 s sin(m 1 ξ 1 ) sin(m 2 ξ 2 ) .

If ξ 1 or ξ 2 is in {0, π}, this implies that H = 0 and therefore H0 (ξ, s) = -H1 (ξ, s) ≤ 0 for (ξ, s) ∈ ∂Q ++ × (0, T ] .

According to Lemma 2.6, H0 satisfies the heat equation (∂ s -∆ ξ ) H0 = 0 in Σ T 0 . Thus we can apply the maximum principle to H0 and obtain H0 (ξ, s) ≤ 0 in Σ T 0 , hence also H 0 (ξ, s) ≤ 0. Integrating over s ∈ (0, ∞) yields

F 0 (ξ) ≤ 0 for ξ ∈ Q ++ .
For H, we cannot apply the maximum principle directly in Σ T 0 , because H is not continuous at (0, 0) ∈ Σ T 0 , but we can apply it in Σ T for any 0 < < T . On the lateral boundary, H vanishes as seen above, and for the initial value at s = we have H(ξ, ) = H0 (ξ, ) + H1 (ξ, ) ≥ H0 (ξ, ) ≥ δ( ) with δ( ) = inf ξ∈Q ++ H0 (ξ, ). Hence by the maximum principle, in Σ T we have H(ξ, s) ≥ min{0, δ( )} .

Now, as we have seen above, H0 (•, s) tends to 0 uniformly as s → 0, hence δ( ) → 0 as → 0, which implies H(ξ, s) ≥ 0 for any s > 0 and ξ ∈ Q ++ . After integrating over s, we finally get F (ξ) ≥ 0 for ξ ∈ Q ++ , and the proof of the Lemma is complete.

Remark 3.8. In conclusion, for this example we find the same "ideal" stability estimate as in the previous one-dimensional example.

3.3. Example 3. Dimension d = 2, kernel (x 2 1 -x 2
2 )|x| -4 . We consider another two-dimensional example where in the notation of Section 1.5.2 d = 2 and p(x) = -1 2π a 12 (x), see (1.29). We show that the complement of the stability zone in this case is strictly larger than the image of the symbol of the integral operator.

3.3.1.

The singular integral equation. We use the same notation for analogous objects as in the preceding example. Therefore in this section, the letters K, K, T etc. are redefined to have new meanings. The kernel and its Fourier transform are now

(3.21) K(x) = x 2 2 -x 2 1 2π|x| 4 , K(ξ) = ξ 2 1 -ξ 2 2 2|ξ| 2 = ξ 2 1 |ξ| 2 - 1 2 .
The normalization is chosen so that the range of the function K is again the interval [-1 2 , 1 2 ]. In fact, this kernel is the same as in the previous example (3.9) after a 45 • rotation of the coordinate system. Therefore if we write the singular integral equation as in (3.10), we can copy verbatim the statement of the previous example concerning the numerical range of the integral operator A Ω (see Lemma 3.5) and the resolvent estimate (3.11). Lemma 3.9. Lemma 3.5 is true for the singular integral equation (3.10) defined with the kernel (3.21).

3.3.2. The discrete system. To the delta-delta discretization

(3.22) λu m -N -2 n∈ω N ,m =n K(x N m -x N n )u n = f (x N m ) (m ∈ ω N )
corresponds the finite section T N of the infinite Toeplitz matrix

T = 1 2π (m 2 -n 2 ) 2 -(m 1 -n 1 ) 2 |m -n| 4 m,n∈Z 2
with zero on the diagonal.

The numerical symbol (symbol of T ) is now defined as The conjecture (3.24) then corresponds to the claim that F attains its maximum on Q in the point τ = (π, 0).

(3.23) F (τ ) = m∈Z 2 ,m =0
To prove (3.25), we first transform the slowly converging double Fourier series

(3.26) F (π, 0) = m∈Z 2 ,m =0 (-1) m 1 m 2 2 -m 2 1 2π(m 2 1 + m 2 2 ) 2
into a rapidly convergent single series. One way to get this is to start with the Poisson summation formula applied to the function f (x) = (x -iy) -1 whose Fourier transform is f (ξ) = 2πi1 + (ξ)e -yξ for y > 0. The result is then valid for all y = 0. This can be used to reduce the double Fourier series for F (ξ) to a single rapidly convergent Fourier series. We are here only interested in the limit t → 0:

(3.29) n∈Z n 2 -y 2 (n 2 + y 2 ) 2 = -π 2 sinh 2 πy .
Hence, by decomposing the double sum m∈Z 2 \{0} as m 1 =0,m 2 ∈Z\{0} + m 1 ∈Z\{0},m 2 ∈Z and using ∞ n=1 1 n 2 = π 2 6 (which can also be obtained from (3.29) by looking at the pole in y = 0) we finally get

(3.30) Λ 0 = m∈Z 2 ,m =0 (-1) m 1 m 2 2 -m 2 1 2π(m 2 1 + m 2 2 ) 2 = π 6 - ∞ n=1 (-1) n π sinh 2 πn .
This series converges rapidly, with 5 terms giving 15 significant digits: Λ 0 = 0.547109903806619... The series is covered by the formulas for IX s with s = 2 in [START_REF] Zucker | The summation of series of hyperbolic functions[END_REF]. The explicit expression for Λ 0 given in the Lemma can be deduced from this.

Remark 3.12. The conformal radius (or logarithmic capacity) of the unit square is known to be [12, Tables]

R = Γ( 1 4 ) 2 4π 3 2
. This implies the remarkable relation (3.31)

πR 2 = 2Λ 0 .
The conjecture that Λ + = Λ 0 is clearly supported by numerical evidence. Here are the results of two different approaches for the approximation of Λ + :

In Table 1 we approximate the numerical symbol F from (3.23) using the Ewald method (2.12) from Proposition 2.2.

(3.32)

F (τ ) ≈ |m 1 |,|m 2 |≤M,m =0 K(m)Γ(2, π|m| 2 )e im•τ + |n 1 |,|n 2 |≤M K(τ + 2πn) e -|τ +2πn| 2 4π
.

We take the maximum of F (τ ) over a regular N × N grid discretizing the period square Q = [-π, π] For comparison, we show in Table 3 the analogous computations for the matrices from Example 2, where Λ + = 0.5.

In the previous Example 2, we were able to prove the equation Λ + = 0.5 using an argument involving the maximum principle for the heat equation, see the proof of Lemma 3.6, in particular (3.18). While we have no proof for the equation Λ + = Λ 0 here, it is possible to use an analogous argument to obtain an upper bound for Λ + . The square Q ++ (of area π 2 ) of the previous example now has to be turned by 45 • and to be replaced by the lozenge (a square of area 2π 2 ) Then one can see again that F (ξ) = 0 on ∂Q ♦ . But now ∂Q ♦ contains two points of discontinuity of F , the origin (0, 0) and the point (2π, 0). Therefore the decomposition F = F 0 + K has to be refined into

Q ♦ = ξ ∈ R 2 | 0 < ξ 1 < 2π; |ξ 2 | < min{ξ 1 , 2π -ξ 1 } . N λ max (T N ) extrap. diff
F (ξ) = F 00 (ξ) + K(ξ) + K(ξ -(2π, 0)) .
The function F 00 defined by this will then be continuous on the closure of Q ♦ . Now one can use the integral representation from Proposition 2.5 similarly to (3.19) and use the maximum principle for the heat equation as before to conclude that For any ξ ∈ Q ♦ , F (ξ) ≥ 0 and F 00 (ξ) ≤ 0 . Unfortunately, this estimate is much less sharp than the estimate by 1 2 in the previous example. To illustrate the behavior of the numerical symbol F (ξ), we present in Figure 2 a surface graph of F on the square Q and a contour plot on the lozenge Q ♦ . The parts exceeding the range of the symbol K are indicated in bright red hues. The maximum at the midpoint (π, 0) of Q ♦ is in clear evidence. Let us summarize the stability result obtained for this example.

This implies 0 ≤ F (ξ) ≤ K(ξ) + K(ξ -(2π, 0)) in Q ♦ ,
Corollary 3.13. Let C = [-Λ + , Λ + ] and λ ∈ C \ C. Then for any N the linear system (3.12) has a unique solution, and there is a uniform resolvent estimate

(3.34) (λI -T N ) -1 L( 2 (Ω N )) ≤ dist(λ, C) -1 . For λ ∈ [-Λ + , - 1
2 ) ∪ ( 1 2 , Λ + ] the integral equation (λI -A Ω )u = f with kernel (3.21) is wellposed in L 2 (Ω), but the corresponding delta-delta approximation scheme (3.22) is unstable.

Example 4. Dimension

d = 2, kernel (x 1 + ix 2 ) 2 |x| -4 .
Let d = 2 and p(x) = -1 1π (a 12 (x)+2ib 12 (x)). The corresponding kernel and its Fourier transform are

(3.35) K(x) = x 2 2 -x 2 1 -2ix 1 x 2 π|x| 4 , K(ξ) = (ξ 1 + iξ 2 ) 2 |ξ 1 + iξ 2 | 2 .
The normalization is chosen so that | K(ξ)| = 1 for ξ ∈ R 2 . We include this example, which has features combining those of the two preceding examples, mainly for purposes of illustration.

Because the singular integral operator and the system matrices of the corresponding delta-delta discretization in this case are non-selfadjoint, we expect to see less trivial relations between spectra and numerical ranges than in the selfadjoint case.

It is obvious from the definition (3.35) that the spectrum Sp(A) of the operator of convolution with K in L 2 (R 2 ) is the unit circle {ξ ∈ C | |ξ| = 1} and that its numerical range is the unit disk. Whereas we do not know the spectrum Sp(A Ω ) for a bounded domain Ω ⊂ R 2 , the numerical range is still the unit disk, compare (1.26), (3.36)

Sp(A Ω ) ⊂ W (A Ω ) = W (A) = {ξ ∈ C | |ξ| ≤ 1}.
For the system matrices T N of the delta-delta discretization scheme, Theorem 1.6 and Lemma 2.1 provide the following relations.

(3.37)

Sp(T N ) ⊂ W (T N ) ⊂ W (T ) = conv M ∈N W (T M ) and W (A) ⊂ W (T ) .
In Figure 3 we show for the case of a square domain Ω and two values of N the spectrum Sp(T N ) (red points), the boundary of the numerical range W (T N ) (red line), and the unit circle, which is the boundary of W (A) (green line). We can see the inclusions from (3.37) between Sp(T N ) and W (T N ), and we can perceive the asymptotic inclusion of W (A) in W (T N ) as N tends to infinity.

We can also see that the eigenvalues of the matrices, in contrast to the numerical range, will not fill the whole unit disk asymptotically. On the other hand, we clearly see the overshoot W (T ) \ W (A Ω ), that is the region of λ where the volume integral equation is uniquely solvable and the operator λI -A Ω is sectorial, so that every L 2 -conforming Galerkin method would converge, whereas the delta-delta scheme is unstable. It appears that the limits for the real part of this overshoot are the same (scaled by a factor 2) as in the previous example, that is ±2Λ 0 with Λ 0 defined in Lemma 3.10. In the quasi-static Maxwell volume integral equation (see Section 1.2), the right hand side and the solution are C d -valued functions, and the singular integral operator is defined as the matrix of second distributional derivatives of the convolution with the free-space Green function g for the Laplace operator, see equation (1.5). If we call this operator A 0 , then it is not the same as the operator A defined by the Cauchy principal value of the integral with the same kernel, but there is a simple relation: Let

A 0 u(x) = -∇ div R d g(x -y)u(y)dy = p.v. R d K(x -y)u(y)dy + 1 d u(x) = (A + 1 d I)u(x) .
This is most easily seen by first using the symmetries of the kernel with respect to reflections at coordinate axes and permutations of the variables in order to deduce that the distribution kernel of The function space is now L 2 (Ω; C d ).

Let us note the explicit form of the kernel, valid in any dimension d ≥ 2, where we consider points in R d as column vectors,

(3.41) K(x) = - 1 ν d (x x - 1 d I d |x| 2 ) |x| -d-2 , with ν d = 2π d 2 d Γ( d 2 )
.

The simplest way to see this is to first look at the symbol of the operator. For this we employ d-dimensional Fourier transformation and use the fact that g(ξ) = |ξ| -2 , hence

(3.42) F (-D 2 g)(ξ) = ξ ξ |ξ| 2 and K(ξ) = ξ ξ -1 d I d |ξ| 2 |ξ| 2 .
We check that tr K = 0 and that K satisfies the spherical cancellation condition. Indeed, in the notation of Lemma 1.5, the (matrix-valued) polynomial p(ξ) is given by

(3.43) p(ξ) = - 1 ν d (ξ ξ - 1 d I d |ξ| 2 ) .
Thus the off-diagonal elements of the matrix p(x) are given by

- 1 ν d x j x k = - b jk (x) ν d ,
and the diagonal elements by We recall the discussion of matrix-valued kernels in Section 2.4 above, in particular the properties of the numerical range stated in Lemma 2.7.

- 1 ν d (x 2 k - 1 d |x| 2 ) = 1 d ν d
The basic stability estimate follows. This justifies (3.47). All the other statements of the proposition are instances of the statements of Section 1.6, in particular Theorem 1.6, and their proofs in Section 2, based on Ewald's method.

What remains is to get information on the numbers Λ

+ ], the integral equation is well-posed in L 2 (Ω; C d ), but the delta-delta discretization scheme is unstable in 2 (Ω N ; C d ).

We will discuss this for the practically relevant cases d = 2, where we get rather precise information, and d = 3, which is the most important case because of its relevance for the DDA method in computational electromagnetics. Under this assumption, one can write simple series expansions for the numbers Λ ± . If all the coordinates of τ are 0 or π, then the off-diagonal elements of the matrix F (τ ) vanish and the 3 eigenvalues are the diagonal elements. Therefore the Fourier series for F (τ ) gives (3.52)

Λ (3) -= m∈Z 3 ,m =0 (-1) m 1 +m 2 4π m 2 1 + m 2 2 -2m 2 3 (m 2 1 + m 2 2 + m 2 3 ) 5 2
, Λ

+ = m∈Z 3 ,m =0 (-1) m 3 4π m 2 1 + m 2 2 -2m 2 3 (m 2 1 + m 2 2 + m 2 3 ) (3)
.

These sums, although not absolutely convergent, appear to converge quite well in the sense of partial sums over cubes, On the other hand, the corresponding DDA scheme is stable in 2 (Z 3 ) if and only if r is either nonreal or contained in the interval ( min , max ). For r ∈ (0, min ] ∪ [ max , ∞) the integral equation (and therefore the dielectric scattering problem) is well-posed, but the DDA scheme is unstable.

To conclude this discussion, we show in Table 4 the result of some computations for the spectrum of the system matrix T N for a cube in three dimensions. One can see convergence to the expected values (3.51), even for rather modest values of N . Compare also [START_REF] Yurkin | Application of the discrete dipole approximation to very large refractive indices: Filtered coupled dipoles revived[END_REF]FIG. 8

2 4

 2 the fact that for ξ ∈ Q and n = 0 we have |ξ + 2πn| ≥ π. Therefore for any δ < π there is a constant C such that (2.21) |H 0 (ξ, s)| ≤ C e -δ s for all 0 < s ≤ 1, ξ ∈ Q, and H 0 (ξ, s) is analytic in ξ for all s.

3 . EXAMPLES 3 . 1 .

 331 Example 1. Dimension d = 1. Finite Hilbert transformation.

3. 1 . 1 .

 11 The singular integral equation. Let a, b ∈ R with a < b. On the interval Ω = (a, b) we consider the singular integral equation, abbreviated as (λI

Lemma 3 . 1 .

 31 The finite Hilbert transformation A Ω is a bounded selfadjoint operator in L 2 (Ω), unitarily equivalent to the operator of multiplication by σ in L 2 (-1, 1) with σ(ξ) = ξ. Both the spectrum Sp(A Ω ) and the closure of the numerical range W (A Ω ) are equal to C = [-1, 1]. For all λ ∈ C \ C and any f ∈ L 2 (Ω), the integral equation (3.1) has a unique solution u ∈ L 2 (Ω), and for the resolvent one has in the L 2 (Ω) operator norm

Lemma 3 . 2 . 1 . 3 . 3 .

 32133 The infinite Toeplitz matrix T defines a bounded selfadjoint operator in 2 (Z), unitarily equivalent to the operator of multiplication by F in L 2 (-π, π) with F given in (3.5). Both the spectrum Sp(T ) and the closure of the numerical range W (T ) are equal to C = [-1, 1]. For all λ ∈ C \ C the operator λI -T is invertible in 2 (Z), and for the resolvent one has in the 2 (Z) operator norm (3.6) (λI -T ) -1 = dist(λ, C) -Corollary The matrix T N of the system (3.3) is selfadjoint with its eigenvalues in C = [-1, 1].

1 C

 1 then one also has (by a standard Galerkin argument) (λI -T ) -1 ≤ C, hence dist(λ, C) ≥ and λ ∈ C. Combining this with Lemma 3.1, we obtain the following description of the stability result for our delta-delta discretisation of the finite Hilbert transform. Theorem 3.4. For λ ∈ C the following are equivalent: (i) The singular integral equation (3.1) has a unique solution u ∈ L 2 (Ω) for any f ∈ L 2 (Ω). (ii) The discretization method (3.3) is stable in the 2 norm.(iii) λ ∈ C, where C = [-1, 1].

5 FIGURE 1 .

 51 FIGURE 1. Contour plot of F (ξ) on the quarter square Q ++ , Example 2.

2

 2 = 0.5471... . Let F (τ ) be as defined in(3.23). Then there existsΛ + ≥ Λ 0 such that F (Q) = C = [-Λ + , Λ + ].Conjecture 3.11. Numerical evidence suggests equality(3.24) Λ + = Λ 0 .Proof of Lemma 3.10. The function F is odd with respect to permutation of ξ 1 and ξ 2 . The decompositionF = F 0 + K with F 0 continuous on Q implies that F takes its maximum Λ + on Q. Therefore its image F (Q) is a closed symmetric interval C = [-Λ + , Λ + ].We are going to show that (3.25) F (π, 0) = Λ 0 .

  and hence by taking the maximum, (3.33) Λ + ≤ 1 .

FIGURE 2 .

 2 FIGURE 2. Numerical symbol for Example 3. Left: F on Q, right: F on Q ♦ .

FIGURE 3 .

 3 FIGURE 3. Spectrum and numerical range. Left: N = 8, right: N = 32.

( 3 .

 3 38) K(x) = -D 2 g(x) = -∂ i ∂ j g(x) i,j∈1,...,d for x = 0 . Then (3.39)

3 . 5 . 1 .

 351 A 0 -A must be a scalar multiple of the d × d identity matrix I d times the Dirac distribution δ 0 , and then determining this multiple by taking traces: tr(-D 2 g) = -∆g = δ 0 = tr( 1 d I d δ 0 ). The singular integral equation. We consider the strongly singular integral equation, still written as (λI -A Ω )u = f , (3.40) λu(x) -p.v. Ω K(x -y)u(y) dy = f (x) with K given in (3.38) .

a 3 . 5 . 2 .

 352 jk (x) , compare Remark 2.4.From our formulas of Section 1.5.1 we find the explicit form (3.41) for our kernel. For d = 2, we have ν d = π, and we recognize the kernels studied in the Examples 2 and 3.The matrix ξ ξ|ξ| 2 is an orthogonal projection matrix, hence its numerical range is the interval[0, 1]. Therefore W ( K(ξ)) = [-1 d , 1 -1 d ]for any ξ = 0. We immediately get the following instance of Proposition 1.2.Lemma 3.14. Let C = [-1 d , 1 -1 d ].Then for all λ ∈ C and any f ∈ L 2 (Ω; C d ), the integral equation (3.40) has a unique solution u ∈ L 2 (Ω; C d ), and there is a resolvent estimate in the L 2 (Ω; C d ) operator norm(3.44) (λI -A Ω ) -1 ≤ dist(λ, C) -1 .The discrete system. With the d × d matrix-valued kernel K and vector-valued functions u and f , we can write the delta-delta discretization (λI -T N )U = F of the integral equation(3.40) in the same form as in the scalar case(3.45) λu m -N -d n∈ω N ,m =n K(x N m -x N n )u n = f (x N m ) , (m ∈ ω N ) ,where now the system matrix T N is of size d|ω N | × d|ω N | and is considered as a linear operator in 2 (ω N ; C d ).

Proposition 3 . 15 .

 315 Let K be the kernel defined in (3.38),(3.41). Then there exist Λ following holds.(i) For τ ∈ Q = [-π, π] d , τ = 0, F (τ ) is a real symmetric matrix with eigenvalues contained in the interval C = [Λ (F (τ )) . (ii) For any N , the numerical range W (T N ) is contained in W (T ) = C. (T N )) .(iii) The delta-delta scheme (3.45) is stable if and only if λ ∈ C \ C, and one has the stability estimate in the 2 (ω N ; C d ) operator norm(3.49) (λI -T N ) -1 ≤ dist(λ, C) -1 .Proof. The matrix K(x) is symmetric for x = 0, implying that also F (τ ) is a symmetric matrix for τ = 0. The symmetry K(-x) = K(x) implies that the matrix elements of F (τ ) are real. Therefore the numerical range ofF (τ ) is the interval [λ -(τ ), λ + (τ )],where λ -(τ ) = min(Sp(F (τ )) , λ + (τ ) = max(Sp(F (τ )) .

±

  and to see whether the inequalities (3.46) are strict. In that case, for λ ∈ [Λ(d) -, -1 d ) ∪ (1 -1 d , Λ

3. 5 . 3 .

 53 Dimension d = 2. . Here the numerical symbol has the form F (τ ) = a(τ ) b(τ ) b(τ ) -a(τ ) with real-valued functions a and b. The eigenvalues are λ ± (τ ) = ± a(τ ) 2 + b(τ ) 2 , implying Λ (d) -= -Λ (d) + . The functions a and b have been studied in the previous examples, a in Example 3 and b in Example 2.

FIGURE 5 .

 5 FIGURE 5. d = 3. Eigenvalues of F (τ ) on middle line of face of Q.

. 3 ]

 3 We do not know whether explicit expressions for these sums exist.By means of the Clausius-Mossotti relation (1.8) one can express the stability results equivalently in terms of the relative permittivity r . of the quasi-static Maxwell volume integral operator corresponds to r ≤ 0. The volume integral equation is therefore well posed in L 2 (Ω) if the relative permittivity r is either non-real or positive.

  is given either by the Fourier series H F in (2.16) or, equivalently, by the lattice sum H P in (2.17). The decomposition F = F 0 + K in Proposition 2.2 corresponds to the decomposition H = H 0 + H 1 with H 0 and H 1 defined in (2.19) and (2.20), and there holds

	0 according to (2.21) and (2.22), but, because of Poisson's summation formula, they are in fact the same H (2.24)

2 4s , one sees that C ξ = O(|ξ| -2 ) as ξ → 0.

Thus we see that H F is integrable as s → ∞ according to

(2.18)

, and H P is integrable as s → F (ξ, s) = H P (ξ, s) , so we can use all of the above estimates for both of them. We can summarize Proposition 2.5. The symbol F (ξ) has the integral representation

(2.23) F (ξ) = ∞ 0 H(ξ, s) ds ,

where H(ξ, s)

  It can be written for t ∈ [-π, π] as

	(3.27)		n∈Z	e int n -iy	= iπ	e yσ(t) sinh(πy)	with σ(t) = -t + π sign t .
	Remark: Euler's formula (3.5) is a simple consequence of this.
	Taking the derivative with respect to y and subtracting the formulas for y and -y leads to
	(3.28)	n∈Z	n 2 -y 2 (n 2 + y 2 ) 2 e int = π	σ(t) sinh(σ(t)y) sinh(πy) -π cosh(σ(t)y) cosh(πy) sinh 2 πy	.

TABLE 1 .

 1 2 . Results are shown for N = 1001, so that the point (π, 0) is included. One sees the rapid convergence of the sums in the Ewald method. Computation of Λ +In Table2we show the maximum eigenvalue of the matrix T N where Ω is the unit square, together with an extrapolated value and its difference with Λ 0 .

		M	Maximum	diff with Λ 0
		1 0.5466820485568409	-0.00043
		2 0.5471099022284376	-1.578e-9
		3 0.5471099038066192	1.11e-16
		4 0.5471099038066192	1.11e-16
	N	λ max (T N )	extrap.	diff with Λ 0
	16 0.541802946417726
	24 0.544571778645890
	36 0.545922219922679 0.547207966733364	9.81e-5
	54 0.546562896841136 0.547141211191569	3.13e-5
	81 0.546860792009930 0.547119678405314	9.77e-6

TABLE 2 .

 2 Computation of max(Sp(T N )), Example 3

TABLE 3 .

 3 Computation of max(Sp(T N )), Example 2

	with 0.5

TABLE 4 .

 4 ]. N λ max (T N ) λ min (T N ) λ max (T N ) -λ min (T N ) Computation of max(Sp(T N ))and min(Sp(T N )), Example 5

	4 0.67730278666935 -0.3896455148525014	1.06694830152185
	8 0.73653727456221 -0.4130173055963489	1.14955458015856
	12 0.75323748914578 -0.4193953119966648	1.17263280114245
	16 0.76017444184544 -0.4220149407429199	1.18218938258836

Acknowledgment. This work was partially supported by a grant from the Niels Hendrik Abel Board. The authors acknowledge support of the Centre Henri Lebesgue ANR-11-LABX-0020-01.

In particular, b(τ ) = 0 for τ ∈ ∂Q, and therefore (3.50) for τ = (π, 0), λ + (τ ) = a(τ ) = Λ 0 with the number Λ 0 = 0.5471... encountered in Example 3, Lemma 3.10. This implies

+ ≥ Λ 0 , and we are in the same situation as in Example 3: Strong numerical evidence suggests that the function λ + attains its maximum in the point τ = (π, 0) and therefore Λ

(2) + = Λ 0 , but we do not have a formal proof for this. The positive eigenvalue λ + is plotted in Figure 4.

-= -0.4260241507272727 , Λ

This implies a length of W (T ) of Λ

+ -Λ

-= 1.1969263735019922 instead of 1, which is the length of W (A), thus an overshoot of almost 20%.