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

TL;DR

This paper analyzes the stability of a simple discretization method for strongly singular integral equations, including applications to electromagnetic scattering, revealing conditions under which the method is stable or unstable.

Contribution

It introduces a novel analysis of a discretization method for singular integral equations, connecting it to finite section methods and spectral properties using Ewald's technique.

Findings

The method can be viewed as a finite section of an infinite Toeplitz matrix.

The spectrum and numerical range of the discretization are characterized.

The method is stable for many parameters but can be unstable even when the continuous problem is well-posed.

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 $\lambda$, there are values of $\lambda$ for which the singular integral equation is well posed, but the discretization method is unstable.

Tables

Table 1. Table 1. Computation of Λ + subscript Λ \Lambda_{+}

$M$	Maximum	diff with ${\Lambda }_0$
1	0.5466820485568409	-0.00043
2	0.5471099022284376	-1.578e-9
3	0.5471099038066192	1.11e-16
4	0.5471099038066192	1.11e-16

Table 2. Table 2. Computation of max ⁡ ( Sp ( T N ) ) Sp superscript 𝑇 𝑁 \max(\mathop{\mathrm{Sp}}\nolimits(T^{N})) , Example 3.8

$N$	${\lambda }_{\max }\left(T^N\right)$	extrap.	diff with ${\Lambda }_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 3. Table 3. Computation of max ⁡ ( Sp ( T N ) ) Sp superscript 𝑇 𝑁 \max(\mathop{\mathrm{Sp}}\nolimits(T^{N})) , Example 3.4

$N$	${\lambda }_{\max }\left(T^N\right)$	extrap.	diff with $0.5$
16	0.4299869696672885
24	0.4526591158216325
36	0.4683227545642122	0.5033301483277116	3.33e-3
54	0.4789372344435390	0.5012512843991882	1.25e-3
81	0.4860451011088278	0.5004526691660266	4.53e-4

Table 4. Table 4. Computation of max ⁡ ( Sp ( T N ) ) Sp superscript 𝑇 𝑁 \max(\mathop{\mathrm{Sp}}\nolimits(T^{N})) and min ⁡ ( Sp ( T N ) ) Sp superscript 𝑇 𝑁 \min(\mathop{\mathrm{Sp}}\nolimits(T^{N})) , Example 3.4

$N$	${\lambda }_{\max }\left(T^N\right)$	${\lambda }_{\min }\left(T^N\right)$	${\lambda }_{\max }\left(T^N\right)-{\lambda }_{\min }\left(T^N\right)$
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

Equations

$$\lambda u(x)-{\int }_{\Omega }K(x,y)u(y)dy=f(x)(x\in \Omega )$$

$$\lambda u_m-\sum _{x_n\in \Omega ,n\ne m}{h}^{d}K(x_m,x_n)u_n=f(x_m)(x_m\in \Omega ).$$

$$\mathrm{curl}\mathrm{curl}u-{\kappa }^2ϵu=i\kappa J.$$

$$\mathrm{curl}\mathrm{curl}u-{\kappa }^{2}u=-{\kappa }^2(1-ϵ)u+i\kappa J.$$

$$u=-(\nabla \mathrm{div}+{\kappa }^2)g_{\kappa }\star (1-ϵ)u+u^{\mathrm{inc}}.$$

$$K(x,y)=-(D^2+{\kappa }^2)g_{\kappa }(x-y).$$

$$D^2{\int }_{{\mathbb{R}}^d}{g}_{\kappa }(x-y)u(y)dy=p.v.{\int }_{{\mathbb{R}}^d}{D}^2{g}_{\kappa }(x-y)u(y)dy-\frac{1}{d}u(x).$$

$$\lambda =\frac{1}{1-{ϵ}_r}-\frac{1}{d}=\frac{d-1+{ϵ}_r}{d(1-{ϵ}_r)}.$$

$$\hat{f}(\xi )=\mathcal{F}f(\xi )={\int }_{{\mathbb{R}}^d}f(x)e^{i\xi \cdot x}dx.$$

$$f(x)={\mathcal{F}}^{-1}\hat{f}(x)=(2\pi {)}^{-d}{\int }_{{\mathbb{R}}^d}\hat{f}(\xi )e^{-ix\cdot \xi }d\xi .$$

$$\tilde{a}(\tau )=\sum _{m\in {\mathbb{Z}}^d}a(m)e^{im\cdot \tau },\tau \in Q=[-\pi ,\pi {]}^d.$$

$$a(m)=(2\pi {)}^{-d}{\int }_Q\tilde{a}(\tau )e^{-im\cdot \tau }d\tau .$$

$$f\mapsto (2\pi {)}^{-\frac{d}{2}}\hat{f}:L^2({\mathbb{R}}^d)\to L^2({\mathbb{R}}^d)\text{mbox}anda\mapsto (2\pi {)}^{-\frac{d}{2}}\tilde{a}:{\ell }^2({\mathbb{Z}}^d)\to L^2(Q)$$

$${\int }_{{\mathbb{R}}^{2d}}\overline{u(x)}k(x-y)v(y)dydx=(2\pi {)}^{-d}{\int }_{{\mathbb{R}}^d}\overline{\hat{u}(\xi )}\hat{k}(\xi )\hat{v}(\xi )d\xi ,$$

$$\sum _{m.n\in {\mathbb{Z}}^d}\overline{a(m)}c(m-n)b(n)=(2\pi {)}^{-d}{\int }_Q\overline{\tilde{a}(\tau )}\tilde{c}(\tau )\tilde{b}(\tau )d\tau .$$

$$\sum _{m\in {\mathbb{Z}}^d}f(m)e^{im\cdot \tau }=\sum _{n\in {\mathbb{Z}}^d}\hat{f}(\tau +2\pi n).$$

$$f\raise-2.15277pt\hbox{\big{|}}_{{\mathbb{Z}}^{d}}\in\ell^{1}({\mathbb{Z}}^{d})\quad\mbox{ and }\;\widehat{f}\in L^{1}({\mathbb{R}}^{d}).$$

$$f(x)=e^{-|x{|}^{2}s}⟺\hat{f}(\xi )=(\frac{\pi }{s}{)}^{\frac{d}{2}}{e}^{-\frac{|\xi {|}^2}{4s}}.$$

$$\sum _{m\in {\mathbb{Z}}^d}{e}^{-|m{|}^{2}s}{e}^{im\cdot \tau }=\sum _{n\in {\mathbb{Z}}^d}(\frac{\pi }{s}{)}^{\frac{d}{2}}{e}^{-\frac{|\tau +2\pi n{|}^2}{4s}}.$$

$$K(x)=p(x)\,|x|^{-d-2}\quad\mbox{where $p$ is a homogeneous polynomial of degree $2$}\,.$$

$$⟨K_{ϵ},\varphi ⟩={\int }_{|x|<ϵ}K(x)(\varphi (x)-\varphi (0))dx+{\int }_{|x|>ϵ}K(x)\varphi (x)dx.$$

$${\int }_{{\text{SS}}^{d-1}}Kds=0.$$

$$⟨K,\varphi ⟩=p.v.\int K(x)\varphi (x)dx=\underset{ϵ\to 0}{\lim }{\int }_{|x|>ϵ}K(x)\varphi (x)dx.$$

$$\mathcal{F}Au=\hat{K}\hat{u}.$$

$$W(A)=\{(u,Au)\mid \Vert u\Vert =1\},$$

$$\Vert (\lambda \mathbb{I}-A{)}^{-1}\Vert \le \mathrm{dist}(\lambda ,W(A){)}^{-1}.$$

$$\lambda u(x)-p.v.{\int }_{\Omega }K(x-y)u(y)dy=f(x)(x\in \Omega ).$$

$$\Vert (\lambda \mathbb{I}-A_{\Omega }{)}^{-1}\Vert \le \mathrm{dist}(\lambda ,\mathcal{C}{)}^{-1}.$$

$$\overline{W(A_{\Omega})}=\overline{W(A)}\quad\mbox{ for any open subset $\Omega\subset{\mathbb{C}}$.}$$

$$\frac{(u,Au)}{(u,u)}=\frac{(u_{\rho ,a},A{u}_{\rho ,a})}{(u_{\rho ,a},u_{\rho ,a})}\in W(A_{\Omega }).$$

Full text

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 $\lambda$, there are values of $\lambda$ for which the singular integral equation is well posed, but the discretization method is unstable.

1. 1 Introduction

2. 1.1 Motivation

3. 1.2 The Discrete Dipole Approximation

4. 1.3 Outline of the paper

5. 1.4 Notation for Fourier transforms and Fourier series

6. 1.5 Kernels and their symbols

7. 1.5.1 Homogeneous kernels

8. 1.5.2 Special kernels

9. 1.6 Delta-delta discretization

10. 2 The discrete system

11. 2.1 Toeplitz structure

12. 2.2 Ewald method

13. 2.3 An integral representation

14. 2.4 Matrix-valued kernels

15. 3 Examples

16. 3.1 Example 1. Dimension $d=1$. Finite Hilbert transformation

17. 3.1.1 The singular integral equation

18. 3.1.2 The discrete system

19. 3.2 Example 2. Dimension $d=2$, kernel $x_{1}x_{2}|x|^{-4}$

20. 3.2.1 The singular integral equation

21. 3.2.2 The discrete system

22. 3.3 Example 3. Dimension $d=2$, kernel $(x_{1}^{2}-x_{2}^{2})|x|^{-4}$

23. 3.3.1 The singular integral equation

24. 3.3.2 The discrete system

25. 3.4 Example 4. Dimension $d=2$, kernel $(x_{1}+ix_{2})^{2}|x|^{-4}$

26. 3.5 Example 5. Dimension $d\geq 2$. Volume Integral Equation for the Quasi-static Maxwell system

27. 3.5.1 The singular integral equation

28. 3.5.2 The discrete system

29. 3.5.3 Dimension $d=2$

30. 3.5.4 Dimension $d=3$

1. Introduction

1.1. Motivation

Introduced almost 50 years ago by Purcell and Pennypacker [13], 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 [18, 3]), 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 [19], 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 non-zero 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 $\Omega\subset{\mathbb{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\geq 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 [14], 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 ${\mathbb{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 [4, 5]) that the same form of the essential spectrum is valid for more general bounded Lipschitz domains. In [14], 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 [14], 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 [20], 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 [20, 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 [15], 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 $\sim 0.093$ ) or very large (larger than $\sim 11.8$). We prove this here for the zero-frequency limit, but expect that it is also true for non-zero frequencies.

1.2. 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}\mid n\in{\mathbb{Z}}^{d}\}$ of meshwidth $h>0$.

Thus a linear integral equation on a bounded domain $\Omega\subset{\mathbb{R}}^{d}$

[TABLE]

will be approximated by the finite dimensional linear system

[TABLE]

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 [10] and, with special emphasis on the two-dimensional situation, in [5]. We write the time harmonic Maxwell equations with normalized frequency $\kappa\in{\mathbb{C}}$ as a second order system for the electric field $u$.

[TABLE]

Here it is assumed that the magnetic permeability is constant (normalized to $1$) in the whole space. If one further assumes that the permittivity $\epsilon$ 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

[TABLE]

Here the right hand side has compact support, and therefore convolution with the outgoing fundamental solution $g_{\kappa}$ of the Helmholtz equation and application of the operator $\nabla\operatorname{\rm div}+\kappa^{2}$ leads to the volume integral equation in distributional form

[TABLE]

Here the incoming field $u^{\rm 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\geq 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\times d$ matrix valued kernel

[TABLE]

The integral operator thus defined involves second order distributional derivatives of the weakly singular kernel $g_{\kappa}(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)

[TABLE]

If we further assume that the permittivity $\epsilon$ equals a constant $\epsilon_{r}\in{\mathbb{C}}\setminus\{1\}$ in $\Omega$, we can divide by $1-\epsilon_{r}$ 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 $\lambda$ defined by the relation

[TABLE]

For $d=3$, this relation $\lambda=\frac{2+\epsilon_{r}}{3(1-\epsilon_{r})}$ is known in the DDA literature as Clausius-Mossotti polarizability, referring to the fact that $\frac{1}{\lambda}$ 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 [9, Section 4.5].

The principal part of the volume integral operator is obtained by taking the limit $\kappa\to 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.

1.3. 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 $\Omega$.

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 $\lambda$ for which any Galerkin scheme of the integral equation would be stable.

In Section 3 we study in detail five representative examples.

Example 3 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 3.4 and 3.8 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 3.4 the ranges of the symbol of the integral operator and of the symbol of the infinite Toeplitz matrix are identical, whereas in Example 3.8 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 3.13, 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.8 and 3.4, respectively.

The kernels studied in Examples 3.4 and 3.8 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 3.4. We study this for dimensions $d\geq 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.8. 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 ${\mathbb{R}}^{d}$.

[TABLE]

Inverse:

[TABLE]

For Fourier series, we use the following notation. For a sequence $a:{\mathbb{Z}}^{d}\to{\mathbb{C}}$, its Fourier series is defined as

[TABLE]

Inverse:

[TABLE]

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

[TABLE]

are unitary (i.e. isometric Hilbert space isomorphisms).

Combining the Parseval formula and the convolution theorem gives

[TABLE]

[TABLE]

From these formulas follows immediately that the operators of convolution with $k$ in $L^{2}({\mathbb{R}}^{d})$ and of discrete convolution with $c$ in $\ell^{2}({\mathbb{Z}}^{d})$ are bounded if and only if the “symbols” $\widehat{k}$ and $\tilde{c}$ are bounded functions belonging to $L^{\infty}({\mathbb{R}}^{d})$ and $L^{\infty}(Q)$, respectively.

Sufficient conditions for this are that $k\in L^{1}({\mathbb{R}}^{d})$ and $c\in\ell^{1}({\mathbb{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

[TABLE]

A sufficient (but in no way necessary) condition for (1.16) to hold for all $\tau$ is that

[TABLE]

If we do not assume f\raise-2.15277pt\hbox{\big{|}}_{{\mathbb{Z}}^{d}}\in\ell^{1}, but only $\widehat{f}\in 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$.

[TABLE]

For this example the Poisson summation formula takes the form (for $\tau\in{\mathbb{R}}^{d}$)

[TABLE]

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 ${\mathbb{R}}^{d}$ that are convolutions with kernel functions of the form

[TABLE]

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 ${\mathbb{R}}^{d}$, positively homogeneous of degree $-d$ and smooth outside of the origin. For a given $\epsilon>0$, one can define a distribution $K_{\epsilon}\in{\mathscr{S}}^{\prime}({\mathbb{R}}^{d})$ that coincides with $K$ on ${\mathbb{R}}^{d}\setminus\{0\}$ by its action on a test function $\phi$ as

[TABLE]

This is independent of $\epsilon$ if and only if $K$ satisfies the cancellation condition on the unit sphere $\SS^{d-1}$

[TABLE]

In this case, we denote the distribution simply by $K$, and we can take the limit $\epsilon\to 0$, thus we get the Cauchy principal value.

[TABLE]

Another consequence of the cancellation condition (1.21) is that the Fourier transform $\widehat{K}$ of the homogeneous distribution $K$ is a bounded function homogeneous of degree [math], smooth outside of the origin and also satisfying the cancellation condition. The operator $A$ of convolution with $K$ is therefore bounded in $L^{2}({\mathbb{R}}^{d})$. Note that in the absence of condition (1.21), $\widehat{K}_{\epsilon}$ would have a logarithmic singularity at [math].

The operator $A$ is diagonalized by Fourier transformation:

[TABLE]

Therefore in $L^{2}({\mathbb{R}}^{d})$, we can obtain information about the spectrum $\mathop{\mathrm{Sp}}\nolimits(A)$ and about the numerical range $W(A)$ from the corresponding easily checked information about the operator of multiplication by the symbol $\widehat{K}$.

We recall that the numerical range of $A$ is defined by

[TABLE]

where $(\cdot,\cdot)$ denotes the Hilbert space inner product. It is convex by the Toeplitz-Hausdorff theorem and it contains the spectrum of $A$. Denote by $\operatorname{\rm im}(\widehat{K})=\widehat{K}({\mathbb{R}}^{d})$ the image (range) of $\widehat{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 $\mathop{\mathrm{Sp}}\nolimits(A)$ is the image $\operatorname{\rm im}(\widehat{K})$, and the closure $\overline{W(A)}$ of the numerical range of $A$ is the convex hull of $\,\operatorname{\rm im}(\widehat{K})$.

It is well known (and easy to prove) that the numerical range allows estimates for the operator norm of the resolvent: For any $\lambda\in{\mathbb{C}}\setminus\overline{W(A)}$,

[TABLE]

It is also monotone with respect to inclusions of subspaces, a property not shared by the spectrum. Given an open set $\Omega\subset{\mathbb{R}}^{d}$, we denote by $A_{\Omega}$ the restriction of the convolution operator $A$ to $L^{2}(\Omega)$ and consider the strongly singular integral equation $(\lambda\mathbb{I}-A_{\Omega})u=f$, or in detail

[TABLE]

From the definition of the numerical range follows immediately the inclusion $W(A_{\Omega})\subset W(A)$.

We can summarize this discussion:

Proposition 1.2**.**

Let ${\mathcal{C}}\subset{\mathbb{C}}$ be a closed convex set such that $\widehat{K}(\xi)\in{\mathcal{C}}$ for all $\xi\in\SS^{d-1}$. Then for all $\lambda\not\in{\mathcal{C}}$ and any $f\in L^{2}(\Omega)$, the integral equation (1.24) has a unique solution $u\in L^{2}(\Omega)$, and there is a resolvent estimate in the $L^{2}(\Omega)$ norm

[TABLE]

Remark 1.3*.*

The same argument implies stability for any Galerkin method: Let $X_{h}$ be any closed subspace of $L^{2}(\Omega)$, and let $A_{h}:X_{h}\to X_{h}$ be the operator defined by restricting the sesquilinear form $(u,Av)$ to $X_{h}\times X_{h}$. Then the statement of Proposition 1.2 remains true if we replace $A_{\Omega}$ by $A_{h}$.

Remark 1.4*.*

Whereas there is, in general, no simple relation between the spectra $\mathop{\mathrm{Sp}}\nolimits(A_{\Omega})$ and $\mathop{\mathrm{Sp}}\nolimits(A)$, for the numerical ranges of our convolution operators with homogeneous kernels we not only have the inclusion $W(A_{\Omega})\subset W(A)$, but also the converse. Namely there holds

[TABLE]

Proof.

The set of Rayleigh quotients $\frac{(u,Au)}{(u,u)}$, where $u\in L^{2}({\mathbb{R}}^{d})\setminus\{0\}$ has compact support, is a dense subset of $W(A)$. We show that it is a subset of $W(A_{\Omega})$: Indeed, let $u$ be such a function and let $\rho>0$ and $a\in{\mathbb{R}}^{d}$ be chosen such that the support of the function $u_{\rho,a}$ defined by $u_{\rho,a}(x)=u(\rho x+a)$ is contained in $\Omega$. Then

[TABLE]

∎

1.5.2. Special kernels

For $d=1$, there is essentially only one non-trivial kernel homogeneous of degree $-d$, namely $K(x)=\frac{1}{x}$.

In ${\mathbb{R}}^{d}$ for $d\geq 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 ${\mathbb{R}}^{d}$ according to (1.22), and we denote by $\widehat{K}$ its Fourier transform.

Lemma 1.5**.**

Let $K$ have the form (1.19) and satisfy (1.21). Then

[TABLE]

Proof.

We first compute the Fourier transform of $p(x)e^{-|x|^{2}s}$, using (1.17)

[TABLE]

The evaluation of these derivatives leads to the following simple result, as we will show:

[TABLE]

For $j,k\in\{1,\dots,d\}$ with $j\neq k$, let

[TABLE]

Any homogeneous polynomial of degree $2$ satisfying the cancellation condition $\int_{\SS^{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 $\partial_{\xi_{j}}e^{-\frac{|\xi|^{2}}{4s}}=-\frac{1}{2s}\xi_{j}e^{-\frac{|\xi|^{2}}{4s}}$ and \partial_{\xi_{j}}^{2}e^{-\frac{|\xi|^{2}}{4s}}=\big{(}-\frac{1}{2s}+\frac{1}{4s^{2}}\xi_{j}^{2}\big{)}e^{-\frac{|\xi|^{2}}{4s}}.

Then for $p=a_{jk}$, we see

[TABLE]

and for $p=b_{jk}$, we see

[TABLE]

Thus in both cases we have

[TABLE]

and (1.28) is proved.

Now we use the definition of the Gamma function

[TABLE]

to write the kernel as an integral over Gaussians

[TABLE]

Taking Fourier transforms and using (1.28), we find with $u=\frac{|\xi|^{2}}{4s}$

[TABLE]

as claimed. ∎

1.6. Delta-delta discretization

Let $N\in{\mathbb{N}}$, fix some origin $a^{N}\in{\mathbb{R}}^{d}$ and define the cubic grid of meshwidth $h=\frac{1}{N}$ by

[TABLE]

We further define

[TABLE]

Then a very simple discretization of the strongly singular integral equation $(\lambda\mathbb{I}-A_{\Omega})u=f$ (1.24) is the following

[TABLE]

or in shorthand $(\lambda\mathbb{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 $\ell^{2}(\omega^{N})$, in the same way as we did above for the integral equation (1.24) in $L^{2}(\Omega)$, 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.19) and (1.21). Then there exists a compact convex set ${\mathcal{C}}\subset{\mathbb{C}}$ such that for any $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$, any $N\in{\mathbb{N}}$ for which $\omega^{N}$ is non-empty, and for any $F\in\ell^{2}(\Omega^{N})$, the system (1.33) has a unique solution, and there is a uniform estimate for the inverse in the $\ell^{2}(\Omega^{N})$ operator norm

[TABLE]

Furthermore, with the strongly singular integral operator $A$ defined above in Section 1.5.1, there holds the inclusion

[TABLE]

Remark 1.7*.*

Note that the inclusion $W(A)\subset{\mathcal{C}}$ implies that for $\lambda\not\in{\mathcal{C}}$ the singular integral equation is uniquely solvable, too, and provides the a priori estimate (1.25), for any domain $\Omega\subset{\mathbb{R}}^{d}$. On the other hand, in order to guarantee stability for $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$, the inclusion may need to be strict, as we shall see in the examples, and then there may be $\lambda\in{\mathcal{C}}\setminus W(A)$ for which the singular integral equation is well posed, but the delta-delta discretization is unstable.

2. The discrete system

Let $T^{N}$ be the matrix representing the discretized integral operator in (1.33):

[TABLE]

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

[TABLE]

hence $T^{N}$ is a finite section of a fixed infinite Toeplitz (discrete convolution) matrix

[TABLE]

Theorem 1.6 will be proved if we can show that $T$ defines a bounded linear operator in $\ell^{2}({\mathbb{Z}}^{d})$ whose numerical range $W(T)$ contains $W(A)$. We can then choose ${\mathcal{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\in\omega^{N}}$, we find

[TABLE]

Here $\tilde{u}(\tau)=\sum_{m\in\omega^{N}}u_{m}e^{im\cdot\tau}$ and $Q=[-\pi,\pi]^{d}$. $F(\tau)$ is the symbol (characteristic function) of the Toeplitz matrix $T$:

[TABLE]

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:\ell^{2}({\mathbb{Z}}^{d})\to\ell^{2}({\mathbb{Z}}^{d})$ is bounded if and only if $F\in L^{\infty}(Q)$.

The closure of $W(T)$ is the closed convex hull of the range $\operatorname{\rm im}(F)=\{F(\tau)\mid\tau\in Q\}$ and is also equal to the closure of the union $\bigcup_{N\in{\mathbb{N}}}W(T^{N})$ .*

Proof: This is immediate from (2.3).

2.2. Ewald method

The problem that makes the statement $F\in L^{\infty}(Q)$ non trivial is that the Fourier series (2.4) is not absolutely convergent. The sequence $(K(m))_{m\in{\mathbb{Z}}^{d}}$ is of order $O(|m|^{-d})$ at infinity and therefore in $\ell^{p}({\mathbb{Z}}^{d})$ for all $p>1$, but not for $p=1$. Its membership in $\ell^{2}({\mathbb{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 $\operatorname{\rm im}(F)$.

We will use a variant of a method introduced by P. P. Ewald [7] in 1921 as a tool to compute slowly converging lattice sums. It has become a routine method for the computation of periodic and quasi-periodic 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 [6] or Section 2.13.3 in the book [2].

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 $\widehat{K}^{P}$ of $K^{P}$ are exponentially decreasing at infinity, so that both the Fourier series for $F^{F}(\tau)$ and the Poisson sum (compare (1.16)) for $F^{P}(\tau)$ 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 [math] to $\infty$ as we did already in Section 1.5.2 above:

[TABLE]

and then to split the integral at a point $\beta^{2}>0$:

[TABLE]

We see that $K^{F}$ is simply the product of $K$ by a function exponentially decreasing at infinity

[TABLE]

with the (upper) incomplete Gamma function (see [1, §6.5])

[TABLE]

Therefore $K^{F}(x)=O(|x|^{2}e^{-\beta^{2}|x|^{2}})$ as $|x|\to\infty$ , and the Fourier series for $F^{F}(\tau)$

[TABLE]

converges rapidly, implying that $F^{F}$ is an analytic function on ${\mathbb{R}}^{d}/(2\pi{\mathbb{Z}})^{d}$.

Consequently, the Fourier series for $F^{P}(\tau)$ converges as slowly as the one for $F(\tau)$, and we use instead the Poisson summation formula (1.16) and write

[TABLE]

We can evaluate $\widehat{K}^{P}$ with the formulas used for $\widehat{K}$ in Lemma 1.5. As in (1.32) we obtain

[TABLE]

Therefore we also obtain a very simple form for the Fourier transform, namely that $\widehat{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 $\tau+2\pi n=0$, and for $\tau\in Q$ this is the term with $n=0$. We can summarize the result.

Proposition 2.2**.**

The symbol $F(\tau)$ of the infinite Toeplitz matrix $T$ is a bounded function given for any $\beta>0$ by the exponentially convergent sums

[TABLE]

In the period cube $Q=[-\pi,\pi]^{d}$, it is $C^{\infty}$ outside of $\,0$, and it has the form

[TABLE]

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\in{\mathbb{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

[TABLE]

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

[TABLE]

we see immediately that it is analytic. For finding $F_{0}(0)$, we can use the following observation.

Lemma 2.3**.**

Let $S\subset{\mathbb{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 $\int_{\SS^{d-1}}p=0$. Then

[TABLE]

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\in{\mathbb{R}}$, the set $\{m\in{\mathbb{Z}}^{d}\mid|m|^{2}=M\}$ is either empty or cubically symmetric. Therefore for $\tau=0$, the two sums in the representation of $F_{0}(\tau)$ are [math]. The last term

[TABLE]

tends to $\nu_{D}p(0)/(4\beta^{2})=0$ as $\tau\to 0$, and hence $F_{0}(0)=0$. ∎

Poposition 2.2 implies Theorem 1.6, where ${\mathcal{C}}$ is the closed convex hull of $\operatorname{\rm im}(F)$. The inclusion $W(A)\subset{\mathcal{C}}$ is easy to see from (2.13):

Given $\epsilon>0$, let $\delta>0$ be such that for $|\tau|<\delta$ we have $|F_{0}(\tau)|<\epsilon$. Since $\widehat{K}$ is homogeneous of degree zero, it takes all of its values already on the ball $B_{\delta}(0)$ of radius $\delta$. Thus

[TABLE]

Taking convex hulls shows that

[TABLE]

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 $\widehat{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 $\xi_{j}\xi_{k}|\xi|^{-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 3.4 below) like

[TABLE]

where for $j=k$ the two terms in the sum do not separately satisfy (1.21). This fits into our framework, however, because

[TABLE]

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.

2.3. An integral representation

We have another look at the Ewald splitting for the numerical symbol $F(\xi)=F^{F}(\xi)+F^{P}(\xi)$ described in (2.6)–(2.12)

[TABLE]

These formulas are valid for any $0<\beta<\infty$. All the sums and integrals are converging absolutely here, and therefore we can interchange sums and integrals.

[TABLE]

From the definition (2.16) of $H^{F}$ and the fact that $|m|\geq 1$ in the sum follows without difficulty that for any $0<\gamma<1$ there exists a constant $C$ such that

[TABLE]

To see the behavior of $H^{P}(\xi,s)$ from (2.17), we decompose

[TABLE]

with

[TABLE]

Now we use the fact that for $\xi\in Q$ and $n\neq 0$ we have $|\xi+2\pi n|\geq\pi$. Therefore for any $\delta<\frac{\pi^{2}}{4}$ there is a constant $C$ such that

[TABLE]

and $H_{0}(\xi,s)$ is analytic in $\xi$ for all $s$.

It remains to analyze the term with $n=0$, i.e. $H_{1}$. It is clear that it vanishes for $\xi=0$, and for every $\xi\neq 0$ there exists a constant $C_{\xi}$ and $0<\gamma<\frac{|\xi|^{2}}{4}$ such that

[TABLE]

Thus $H_{1}(\xi,s)$ is integrable over $s\in(0,\infty)$ for all $\xi$, but there is no uniform bound for $C_{\xi}$: Considering $\sup_{s>0}s^{-2}e^{-\frac{|\xi|^{2}}{4s}}$, one sees that $C_{\xi}=O(|\xi|^{-2})$ as $\xi\to 0$.

Thus we see that $H^{F}$ is integrable as $s\to\infty$ according to (2.18), and $H^{P}$ is integrable as $s\to 0$ according to (2.21) and (2.22), but, because of Poisson’s summation formula, they are in fact the same

[TABLE]

so we can use all of the above estimates for both of them. We can summarize

Proposition 2.5**.**

The symbol $F(\xi)$ has the integral representation

[TABLE]

where $H(\xi,s)$ 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}+\widehat{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

[TABLE]

In these integrals, the functions $s\mapsto H_{0}(\xi,s)$, $s\mapsto H_{1}(\xi,s)$, and $s\mapsto H(\xi,s)$ are integrable on $(0,\infty)$ for any $\xi\in Q$, for any $\xi\in{\mathbb{R}}^{d}\setminus\{0\}$, and for any $\xi\in Q\setminus\{0\}$, respectively.

The integral representations (2.23) and (2.24) will be used below to get bounds for the function $F(\xi)$ from estimates for $H(\xi,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

[TABLE]

are solutions of the heat equation

[TABLE]

2.4. Matrix-valued kernels

Until now, we have considered kernel functions with values in ${\mathbb{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

[TABLE]

Theorem 1.6 remains literally true, but for the construction of the set ${\mathcal{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 $\bigcup_{\tau\in Q}W(F(\tau))$ and is also equal to the closure of the union $\bigcup_{N\in{\mathbb{N}}}W(T^{N})$.

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

[TABLE]

Here $\tilde{u}(\tau)=\sum_{m\in\omega^{N}}u_{m}e^{im\cdot\tau}$, and $F(\tau)$ is the matrix-valued symbol of the block Toeplitz matrix T=\big{(}K(m-n)\big{)}_{m,n\in{\mathbb{Z}}^{d}}:

[TABLE]

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.

3. Examples

3.1. Example 1. Dimension $d=1$. Finite Hilbert transformation

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.

3.1.1. The singular integral equation

Let $a,b\in{\mathbb{R}}$ with $a<b$. On the interval $\Omega=(a,b)$ we consider the singular integral equation, abbreviated as $(\lambda\mathbb{I}-A_{\Omega})u=f$,

[TABLE]

The integral is understood in the Cauchy principal value sense. The kernel function $K(x)=\frac{1}{i\pi x}$ has the Fourier transform

[TABLE]

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

The finite Hilbert transformation $A_{\Omega}$ and its spectral theory are also well studied classical objects, see for example [11]. Here the spectrum depends on the function space; for $L^{p}(\Omega)$ it is strongly dependent on $p$, but not on $\Omega$, as long as $\Omega$ is a proper subinterval of ${\mathbb{R}}$. For $p=2$ one has the following description.

Lemma 3.1**.**

The finite Hilbert transformation $A_{\Omega}$ is a bounded selfadjoint operator in $L^{2}(\Omega)$, unitarily equivalent to the operator of multiplication by $\sigma$ in $L^{2}(-1,1)$ with $\sigma(\xi)=\xi$. Both the spectrum $\mathop{\mathrm{Sp}}\nolimits(A_{\Omega})$ and the closure of the numerical range $\overline{W(A_{\Omega})}$ are equal to ${\mathcal{C}}=[-1,1]$. For all $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$ and any $f\in L^{2}(\Omega)$, the integral equation (3.1) has a unique solution $u\in L^{2}(\Omega)$, and for the resolvent one has in the $L^{2}(\Omega)$ operator norm

[TABLE]

Explicit formulas for the resolvent are known. For the infinite Hilbert transformation this is trivially obtained by algebra:

[TABLE]

and for the finite Hilbert transformation, formulas for the resolvent can be found for example in [16] or [17].

3.1.2. The discrete system

We use the notation of Section 1.2 with $d=1$, in particular $x^{N}_{m}=a^{N}+\frac{m}{N}$ and $\omega^{N}=\{m\in{\mathbb{Z}}\mid x^{N}_{m}\in\Omega\}$. The simple delta-delta discretization of our singular integral equation (3.1) is

[TABLE]

The system matrix $T^{N}$ with matrix elements $\frac{1}{i\pi N}\frac{1}{x^{N}_{m}-x^{N}_{n}}$ ($m,n\in\omega^{N}$) is a finite section of the infinite Toeplitz matrix

[TABLE]

The symbol $F(\tau)$ is now given by the Fourier series

[TABLE]

This series converges for all $t\in Q$ to the well known saw-tooth function

[TABLE]

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.

Lemma 3.2**.**

The infinite Toeplitz matrix $T$ defines a bounded selfadjoint operator in $\ell^{2}({\mathbb{Z}})$, unitarily equivalent to the operator of multiplication by $F$ in $L^{2}(-\pi,\pi)$ with $F$ given in (3.5). Both the spectrum $\mathop{\mathrm{Sp}}\nolimits(T)$ and the closure of the numerical range $\overline{W(T)}$ are equal to ${\mathcal{C}}=[-1,1]$. For all $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$ the operator $\lambda\mathbb{I}-T$ is invertible in $\ell^{2}({\mathbb{Z}})$, and for the resolvent one has in the $\ell^{2}({\mathbb{Z}})$ operator norm

[TABLE]

Corollary 3.3**.**

The matrix $T^{N}$ of the system (3.3) is selfadjoint with its eigenvalues in ${\mathcal{C}}=[-1,1]$. For $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$, there is a uniform resolvent estimate in the $\ell^{2}$ operator norm

[TABLE]

The converse is also true: If there is a uniform stability estimate

[TABLE]

then one also has (by a standard Galerkin argument) $\|(\lambda\mathbb{I}-T)^{-1}\|\leq C$, hence $\mathrm{dist}(\lambda,{\mathcal{C}})\geq\frac{1}{C}$ and $\lambda\not\in{\mathcal{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 $\lambda\in{\mathbb{C}}$ the following are equivalent:

(i) The singular integral equation (3.1) has a unique solution $u\in L^{2}(\Omega)$ for any $f\in L^{2}(\Omega)$.

(ii) The discretization method (3.3) is stable in the $\ell^{2}$ norm.

(iii) $\lambda\not\in{\mathcal{C}}$, where ${\mathcal{C}}=[-1,1]$.

For such $\lambda$, there is an estimate for the operator norms*

[TABLE]

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)=-\frac{1}{\pi}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

[TABLE]

For $\Omega\subset{\mathbb{R}}^{2}$, we consider the singular integral equation $(\lambda\mathbb{I}-A_{\Omega})u=f$ as in (1.24)

[TABLE]

Observing that the range of the function $\widehat{K}$ is the interval $[-\frac{1}{2},\frac{1}{2}]$, we can formulate the result of Proposition 1.2 as follows

Lemma 3.5**.**

Let ${\mathcal{C}}=[-\frac{1}{2},\frac{1}{2}]$. For $\Omega={\mathbb{R}}^{2}$, both the spectrum $\mathop{\mathrm{Sp}}\nolimits(A_{\Omega})$ and the closure of the numerical range $\overline{W(A_{\Omega})}$ in $L^{2}(\Omega)$ are equal to ${\mathcal{C}}$. For any open subset $\Omega\subset{\mathbb{R}}^{2}$, the closure of the numerical range in $L^{2}(\Omega)$ satisfies $\overline{W(A_{\Omega})}\subset{\mathcal{C}}$, and there is a resolvent estimate in the $L^{2}(\Omega)$ operator norm

[TABLE]

3.2.2. The discrete system

Let now $\Omega$ be a bounded domain in ${\mathbb{R}}^{2}$. In the notation of Section 1.6 with $d=2$, the regular grid consists of the points $x^{N}_{m}=a^{N}+\frac{m}{N}$, indexed by $\omega^{N}=\{m\in{\mathbb{Z}}^{2}\mid x^{N}_{m}\in\Omega\}$. The simple delta-delta discretization of our singular integral equation (3.10) is

[TABLE]

The system matrix $T^{N}$ is now a finite section of the infinite Toeplitz matrix

[TABLE]

Its symbol is therefore given by the Fourier series for $\tau\in Q=[-\pi,\pi]^{2}$

[TABLE]

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)

[TABLE]

In addition, we know from (3.13) that $F$ vanishes on the boundary of $Q$, hence

[TABLE]

In the previous example, we used the explicit expression of $F(\tau)$ 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 $\widehat{K}$ on $\partial Q$, which implied that the closed convex hull of $\operatorname{\rm im}(F)$ was the same as the convex hull of $\operatorname{\rm im}(\widehat{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(\tau)$ be as defined in (3.13). Then for any $\tau\in Q$

[TABLE]

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 ${\mathcal{C}}=[-\frac{1}{2},\frac{1}{2}]$ and $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$. Then for any $N$ the linear system (3.12) has a unique solution, and there is a uniform resolvent estimate

[TABLE]

Proof.

For symmetry reasons, it is sufficient to prove (3.16) for $\tau\in Q_{++}=(0,\pi)^{2}$. For the proof of Lemma 3.6, we will show the following:

[TABLE]

This implies $0\leq F(\xi)\leq\widehat{K}(\xi)\leq\frac{1}{2}$, hence (3.16) .

We use the integral representations from Proposition 2.5

[TABLE]

According to (2.19), $H_{0}(\xi,s)=H(\xi,s)-H_{1}(\xi,s)$ with

[TABLE]

Let $0\leq\epsilon<T$ and $\Sigma_{\epsilon}^{T}=Q_{++}\times(\epsilon,T)$. In $\Sigma_{\epsilon}^{T}$, we want to use the maximum principle for the heat equation (see Lemma 2.6) for the functions $\tilde{H}_{0}(\xi,s)=s^{-1}H_{0}(\xi,s)$ and $\tilde{H}(\xi,s)=s^{-1}H(\xi,s)$.

Since $\tilde{H}_{1}(\xi,s)=s^{-1}H_{1}(\xi,s)$ is continuous for $(\xi,s)\in{\mathbb{R}}^{2}\times[0,\infty)\setminus\{0,0\}$ and the Poisson series

[TABLE]

converges uniformly for $(\xi,s)\in\overline{\Sigma_{0}^{T}}$ for all $T>0$, we see that $\tilde{H}_{0}$ is continuous in $\overline{\Sigma_{0}^{T}}$ with initial value $\tilde{H}_{0}(\xi,0)=0$. On the lateral boundary we use the Fourier representation (see (2.16))

[TABLE]

If $\xi_{1}$ or $\xi_{2}$ is in $\{0,\pi\}$, this implies that $\tilde{H}=0$ and therefore

[TABLE]

According to Lemma 2.6, $\tilde{H}_{0}$ satisfies the heat equation $(\partial_{s}-\Delta_{\xi})\tilde{H}_{0}=0$ in $\Sigma_{0}^{T}$. Thus we can apply the maximum principle to $\tilde{H}_{0}$ and obtain $\tilde{H}_{0}(\xi,s)\leq 0$ in $\Sigma_{0}^{T}$, hence also $H_{0}(\xi,s)\leq 0$. Integrating over $s\in(0,\infty)$ yields

[TABLE]

For $\tilde{H}$, we cannot apply the maximum principle directly in $\Sigma_{0}^{T}$, because $\tilde{H}$ is not continuous at $(0,0)\in\overline{\Sigma_{0}^{T}}$, but we can apply it in $\Sigma_{\epsilon}^{T}$ for any $0<\epsilon<T$. On the lateral boundary, $\tilde{H}$ vanishes as seen above, and for the initial value at $s=\epsilon$ we have

[TABLE]

with $\delta(\epsilon)=\inf_{\xi\in Q_{++}}\tilde{H}_{0}(\xi,\epsilon)$. Hence by the maximum principle, in $\overline{\Sigma_{\epsilon}^{T}}$ we have

[TABLE]

Now, as we have seen above, $\tilde{H}_{0}(\cdot,s)$ tends to [math] uniformly as $s\to 0$, hence $\delta(\epsilon)\to 0$ as $\epsilon\to 0$, which implies $\tilde{H}(\xi,s)\geq 0$ for any $s>0$ and $\xi\in Q_{++}$. After integrating over $s$, we finally get $F(\xi)\geq 0$ for $\xi\in 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_{1}^{2}-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)=-\frac{1}{2\pi}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$, $\widehat{K}$, $T$ etc. are redefined to have new meanings. The kernel and its Fourier transform are now

[TABLE]

The normalization is chosen so that the range of the function $\widehat{K}$ is again the interval $[-\frac{1}{2},\frac{1}{2}]$.

In fact, this kernel is the same as in the previous example (3.9) after a $45^{\circ}$ 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_{\Omega}$ (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

[TABLE]

corresponds the finite section $T^{N}$ of the infinite Toeplitz matrix

[TABLE]

The numerical symbol (symbol of $T$) is now defined as

[TABLE]

Lemma 3.10**.**

Let

[TABLE]

Let $F(\tau)$ be as defined in (3.23). Then there exists $\Lambda_{+}\geq\Lambda_{0}$ such that $F(Q)={\mathcal{C}}=[-\Lambda_{+},\Lambda_{+}]$.

Conjecture 3.11**.**

Numerical evidence suggests equality

[TABLE]

Proof of Lemma 3.10. The function $F$ is odd with respect to permutation of $\xi_{1}$ and $\xi_{2}$. The decomposition $F=F_{0}+\widehat{K}$ with $F_{0}$ continuous on $Q$ implies that $F$ takes its maximum $\Lambda_{+}$ on $Q$. Therefore its image $F(Q)$ is a closed symmetric interval ${\mathcal{C}}=[-\Lambda_{+},\Lambda_{+}]$. We are going to show that

[TABLE]

The conjecture (3.24) then corresponds to the claim that $F$ attains its maximum on $Q$ in the point $\tau=(\pi,0)$.

To prove (3.25), we first transform the slowly converging double Fourier series

[TABLE]

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 $\widehat{f}(\xi)=2\pi i\mathds{1}_{+}(\xi)e^{-y\xi}$ for $y>0$. The result is then valid for all $y\neq 0$. It can be written for $t\in[-\pi,\pi]$ as

[TABLE]

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

[TABLE]

This can be used to reduce the double Fourier series for $F(\xi)$ to a single rapidly convergent Fourier series. We are here only interested in the limit $t\to 0$:

[TABLE]

Hence, by decomposing the double sum $\sum_{m\in{\mathbb{Z}}^{2}\setminus\{0\}}$ as $\sum_{m_{1}=0,m_{2}\in{\mathbb{Z}}\setminus\{0\}}+\sum_{m_{1}\in{\mathbb{Z}}\setminus\{0\},m_{2}\in{\mathbb{Z}}}$ and using $\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\frac{\pi^{2}}{6}$ (which can also be obtained from (3.29) by looking at the pole in $y=0$) we finally get

[TABLE]

This series converges rapidly, with $5$ terms giving $15$ significant digits: $\Lambda_{0}=0.547109903806619...$ The series is covered by the formulas for $\mathrm{IX}_{s}$ with $s=2$ in [21]. The explicit expression for $\Lambda_{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]

[TABLE]

This implies the remarkable relation

[TABLE]

The conjecture that $\Lambda_{+}=\Lambda_{0}$ is clearly supported by numerical evidence. Here are the results of two different approaches for the approximation of $\Lambda_{+}$:

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

[TABLE]

We take the maximum of $F(\tau)$ over a regular $N\times N$ grid discretizing the period square $Q=[-\pi,\pi]^{2}$. Results are shown for $N=1001$, so that the point $(\pi,0)$ is included. One sees the rapid convergence of the sums in the Ewald method.

In Table 2 we show the maximum eigenvalue of the matrix $T^{N}$ where $\Omega$ is the unit square, together with an extrapolated value and its difference with $\Lambda_{0}$.

For comparison, we show in Table 3 the analogous computations for the matrices from Example 3.4, where $\Lambda_{+}=0.5$.

In the previous Example 3.4, we were able to prove the equation $\Lambda_{+}=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 $\Lambda_{+}=\Lambda_{0}$ here, it is possible to use an analogous argument to obtain an upper bound for $\Lambda_{+}$. The square $Q_{++}$ (of area $\pi^{2}$) of the previous example now has to be turned by $45^{\circ}$ and to be replaced by the lozenge (a square of area $2\pi^{2})$

[TABLE]

Then one can see again that $F(\xi)=0$ on $\partial Q_{\lozenge}$. But now $\partial Q_{\lozenge}$ contains two points of discontinuity of $F$, the origin $(0,0)$ and the point $(2\pi,0)$. Therefore the decomposition $F=F_{0}+\widehat{K}$ has to be refined into

[TABLE]

The function $F_{00}$ defined by this will then be continuous on the closure of $Q_{\lozenge}$. 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

[TABLE]

This implies $0\leq F(\xi)\leq\widehat{K}(\xi)+\widehat{K}(\xi-(2\pi,0))$ in $Q_{\lozenge}$, and hence by taking the maximum,

[TABLE]

Unfortunately, this estimate is much less sharp than the estimate by $\frac{1}{2}$ in the previous example.

To illustrate the behavior of the numerical symbol $F(\xi)$, we present in Figure 2 a surface graph of $F$ on the square $Q$ and a contour plot on the lozenge $Q_{\lozenge}$. The parts exceeding the range of the symbol $\widehat{K}$ are indicated in bright red hues. The maximum at the midpoint $(\pi,0)$ of $Q_{\lozenge}$ is in clear evidence.

Let us summarize the stability result obtained for this example.

Corollary 3.13**.**

Let ${\mathcal{C}}=[-\Lambda_{+},\Lambda_{+}]$ and $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$. Then for any $N$ the linear system (3.12) has a unique solution, and there is a uniform resolvent estimate

[TABLE]

For $\lambda\in[-\Lambda_{+},-\frac{1}{2})\cup(\frac{1}{2},\Lambda_{+}]$ the integral equation $(\lambda\mathbb{I}-A_{\Omega})u=f$ with kernel (3.21) is well-posed in $L^{2}(\Omega)$, but the corresponding delta-delta approximation scheme (3.22) is unstable.

3.4. Example 4. Dimension $d=2$, kernel $(x_{1}+ix_{2})^{2}|x|^{-4}$

Let $d=2$ and $p(x)=-\frac{1}{1\pi}(a_{12}(x)+2ib_{12}(x))$. The corresponding kernel and its Fourier transform are

[TABLE]

The normalization is chosen so that $|\widehat{K}(\xi)|=1$ for $\xi\in{\mathbb{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 $\mathop{\mathrm{Sp}}\nolimits(A)$ of the operator of convolution with $K$ in $L^{2}({\mathbb{R}}^{2})$ is the unit circle $\{\xi\in{\mathbb{C}}\mid|\xi|=1\}$ and that its numerical range is the unit disk. Whereas we do not know the spectrum $\mathop{\mathrm{Sp}}\nolimits(A_{\Omega})$ for a bounded domain $\Omega\subset{\mathbb{R}}^{2}$, the numerical range is still the unit disk, compare (1.26),

[TABLE]

For the system matrices $T^{N}$ of the delta-delta discretization scheme, Theorem 1.6 and Lemma 2.1 provide the following relations.

[TABLE]

In Figure 3 we show for the case of a square domain $\Omega$ and two values of $N$ the spectrum $\mathop{\mathrm{Sp}}\nolimits(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 $\mathop{\mathrm{Sp}}\nolimits(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)\setminus W(A_{\Omega})$, that is the region of $\lambda$ where the volume integral equation is uniquely solvable and the operator $\lambda\mathbb{I}-A_{\Omega}$ 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 $\pm 2\Lambda_{0}$ with $\Lambda_{0}$ defined in Lemma 3.10.

3.5. Example 5. Dimension $d\geq 2$. Volume Integral Equation for the Quasi-static Maxwell system

In the quasi-static Maxwell volume integral equation (see Section 1.2), the right hand side and the solution are ${\mathbb{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

[TABLE]

Then

[TABLE]

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 $A^{0}-A$ must be a scalar multiple of the $d\times d$ identity matrix $\mathbb{I}_{d}$ times the Dirac distribution $\delta_{0}$, and then determining this multiple by taking traces: $\mathop{\mathrm{tr}}(-D^{2}g)=-\Delta g=\delta_{0}=\mathop{\mathrm{tr}}(\frac{1}{d}\mathbb{I}_{d}\delta_{0})$.

3.5.1. The singular integral equation

We consider the strongly singular integral equation, still written as $(\lambda\mathbb{I}-A_{\Omega})u=f$,

[TABLE]

The function space is now $L^{2}(\Omega;{\mathbb{C}}^{d})$.

Let us note the explicit form of the kernel, valid in any dimension $d\geq 2$, where we consider points in ${\mathbb{R}}^{d}$ as column vectors,

[TABLE]

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 $\widehat{g}(\xi)=|\xi|^{-2}$, hence

[TABLE]

We check that $\mathop{\mathrm{tr}}\widehat{K}=0$ and that $\widehat{K}$ satisfies the spherical cancellation condition. Indeed, in the notation of Lemma 1.5, the (matrix-valued) polynomial $p(\xi)$ is given by

[TABLE]

Thus the off-diagonal elements of the matrix $p(x)$ are given by

[TABLE]

and the diagonal elements by

[TABLE]

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 $\nu_{d}=\pi$, and we recognize the kernels studied in the Examples 3.4 and 3.8.

The matrix $\frac{\xi\,\xi^{\top}}{|\xi|^{2}}$ is an orthogonal projection matrix, hence its numerical range is the interval $[0,1]$. Therefore $W(\widehat{K}(\xi))=[-\frac{1}{d},1-\frac{1}{d}]$ for any $\xi\neq 0$. We immediately get the following instance of Proposition 1.2.

Lemma 3.14**.**

Let ${\mathcal{C}}=[-\frac{1}{d},1-\frac{1}{d}]$. Then for all $\lambda\not\in{\mathcal{C}}$ and any $f\in L^{2}(\Omega;{\mathbb{C}}^{d})$, the integral equation (3.40) has a unique solution $u\in L^{2}(\Omega;{\mathbb{C}}^{d})$, and there is a resolvent estimate in the $L^{2}(\Omega;{\mathbb{C}}^{d})$ operator norm

[TABLE]

3.5.2. The discrete system

With the $d\times d$ matrix-valued kernel $K$ and vector-valued functions $u$ and $f$, we can write the delta-delta discretization $(\lambda\mathbb{I}-T^{N})U=F$ of the integral equation (3.40) in the same form as in the scalar case

[TABLE]

where now the system matrix $T^{N}$ is of size $d|\omega^{N}|\times d|\omega^{N}|$ and is considered as a linear operator in $\ell^{2}(\omega^{N};{\mathbb{C}}^{d})$.

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.

The basic stability estimate follows.

Proposition 3.15**.**

Let $K$ be the kernel defined in (3.38), (3.41). Then there exist $\Lambda_{-}^{(d)},\Lambda_{+}^{(d)}\in{\mathbb{R}}$ with

[TABLE]

*such that the following holds.

(i) For $\tau\in Q=[-\pi,\pi]^{d}$, $\tau\neq 0$, $F(\tau)$ is a real symmetric matrix with eigenvalues contained in the interval ${\mathcal{C}}=[\Lambda_{-}^{(d)},\Lambda_{+}^{(d)}]$,*

[TABLE]

*(ii) For any $N$, the numerical range $W(T^{N})$ is contained in $W(T)={\mathcal{C}}$.

[TABLE]

(iii) The delta-delta scheme (3.45) is stable if and only if $\lambda\in{\mathbb{C}}\setminus{\mathcal{C}}$, and one has the stability estimate in the $\ell^{2}(\omega^{N};{\mathbb{C}}^{d})$ operator norm

[TABLE]

Proof.

The matrix $K(x)$ is symmetric for $x\neq 0$, implying that also $F(\tau)$ is a symmetric matrix for $\tau\neq 0$. The symmetry $K(-x)=K(x)$ implies that the matrix elements of $F(\tau)$ are real. Therefore the numerical range of $F(\tau)$ is the interval $[\lambda_{-}(\tau),\lambda_{+}(\tau)]$, where

[TABLE]

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 $\Lambda_{\pm}^{(d)}$ and to see whether the inequalities (3.46) are strict. In that case, for $\lambda\in[\Lambda_{-}^{(d)},-\frac{1}{d})\cup(1-\frac{1}{d},\Lambda_{+}^{(d)}]$, the integral equation is well-posed in $L^{2}(\Omega;{\mathbb{C}}^{d})$, but the delta-delta discretization scheme is unstable in $\ell^{2}(\Omega^{N};{\mathbb{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.

3.5.3. Dimension $d=2$

. Here the numerical symbol has the form

[TABLE]

with real-valued functions $a$ and $b$. The eigenvalues are $\lambda_{\pm}(\tau)=\pm\sqrt{a(\tau)^{2}+b(\tau)^{2}}$, implying $\Lambda_{-}^{(d)}=-\Lambda_{+}^{(d)}$. The functions $a$ and $b$ have been studied in the previous examples, $a$ in Example 3.8 and $b$ in Example 3.4.

In particular, $b(\tau)=0$ for $\tau\in\partial Q$, and therefore

[TABLE]

with the number $\Lambda_{0}=0.5471...$ encountered in Example 3.8, Lemma 3.10. This implies

[TABLE]

and we are in the same situation as in Example 3.8: Strong numerical evidence suggests that the function $\lambda_{+}$ attains its maximum in the point $\tau=(\pi,0)$ and therefore $\Lambda_{+}^{(2)}=\Lambda_{0}$, but we do not have a formal proof for this. The positive eigenvalue $\lambda_{+}$ is plotted in Figure 4.

In any case, we have proved that in dimension $d=2$ for any bounded open set $\Omega\subset{\mathbb{R}}^{2}$ and any $\lambda\in(-0.5471,-0.5)\cup(0.5,0.5471)$ the delta-delta scheme $(\lambda\mathbb{I}-T^{N})U=F$ is unstable in $\ell^{2}(\Omega^{N};{\mathbb{C}}^{2})$ as $N\to\infty$, whereas the integral equation $(\lambda\mathbb{I}-A_{\Omega})u=f$ is well posed in $L^{2}(\Omega;{\mathbb{C}}^{2})$.

3.5.4. Dimension $d=3$

. The three eigenvalues $\lambda_{j}$ of $F(\tau)$ satisfy $\lambda_{1}+\lambda_{2}+\lambda_{3}=0$. Numerically, one sees that the minimal and maximal values are attained on the intersection of the boundary of $Q=[-\pi,\pi]^{3}$ with the coordinate planes. In Figure 5 we show a graph of the three eigenvalues on the line $\{(\pi,y,0)\mid y\in[-\pi,\pi]\}$. The values $-\frac{1}{3}$ and $\frac{2}{3}$ are shown as dashed lines.

This suggests $\Lambda_{-}^{(3)}=\min\mathop{\mathrm{Sp}}\nolimits(F((\pi,\pi,0)))$ and $\Lambda_{+}^{(3)}=\max\mathop{\mathrm{Sp}}\nolimits(F((\pi,0,0)))$.

The computed values are

[TABLE]

This implies a length of $W(T)$ of $\Lambda_{+}^{(3)}-\Lambda_{-}^{(3)}=1.1969263735019922$ instead of $1$, which is the length of $W(A)$, thus an overshoot of almost $20\%$.

Under this assumption, one can write simple series expansions for the numbers $\Lambda_{\pm}^{(3)}$. If all the coordinates of $\tau$ are [math] or $\pi$, then the off-diagonal elements of the matrix $F(\tau)$ vanish and the $3$ eigenvalues are the diagonal elements. Therefore the Fourier series for $F(\tau)$ gives

[TABLE]

These sums, although not absolutely convergent, appear to converge quite well in the sense of partial sums over cubes,

[TABLE]

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 $\epsilon_{r}$. Let

[TABLE]

The numerical range $\lambda\in[-\frac{1}{3},\frac{2}{3}]$ of the quasi-static Maxwell volume integral operator corresponds to $\epsilon_{r}\leq 0$. The volume integral equation is therefore well posed in $L^{2}(\Omega)$ if the relative permittivity $\epsilon_{r}$ is either non-real or positive.

On the other hand, the corresponding DDA scheme is stable in $\ell^{2}({\mathbb{Z}}^{3})$ if and only if $\epsilon_{r}$ is either non-real or contained in the interval $(\epsilon_{\min},\epsilon_{max})$. For $\epsilon_{r}\in(0,\epsilon_{\min}]\cup[\epsilon_{\max},\infty)$ 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 [20, FIG. 8].

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.