On the stability of time-domain integral equations for acoustic wave propagation

We give a principled approach for the selection of a boundary integral, retarded potential representation for the solution of scattering problems for the wave equation in an exterior domain.

Many numerical approaches to solving this problem express the solution in terms of the retarded potentials derived from the fundamental solution to the wave equation. By analogy with the time harmonic case, we define the single and double layer potentials by the formulae for x ∈ Γ c . Here,μ(y, s) = ∂ s µ(y, s). Suppose now that we represent the solution in the form To enforce the boundary condition on Γ × [0, ∞), we take the limit of this representation from the exterior of D, yielding the integral equation µ(x, t) 2 + N µ(x, t) + K(a(y)μ + b(y)µ)(x, t) = g(x, t), for (x, t) ∈ Γ × [0, ∞). Here N is the weakly singular principal part of the double layer restricted to Γ × Γ, and K is the single layer on Γ × Γ.
There is a substantial literature on time-domain integral equations for scattering problems, which we do not seek to review here. (See, for example, [3,4,11,5,6,10] and the references therein.) We simply note that, while prior work has considered the use of the single layer alone, the double layer alone, as well as linear combinations of the two, there is relatively little discussion of a principled analytic approach for selecting one representation over another. Using classical results from scattering theory of Lax, Morawetz, and Phillips [8], we provide criteria for choosing the coefficients (a(y) and b(y)) above, at least when Γ is a nontrapping obstacle. This is a hypothesis about the behavior of reflected light rays in the exterior of D, which is always satisfied if D is convex or star-shaped. These results imply that if the disturbance g(y, t) has compact support in time, then the solution to the mixed Cauchy problem in equation (1) will eventually exhibit local exponential decay.
The principle that guides our selection of representation is that, ideally, the source µ(y, t) should display the same rate of exponential decay in time as the solution itself. If this is not the case, then the mapping from µ → u(x, t) in (3) must, after a certain time, involve some catastrophic cancellation, and hence can be expected to lose relative accuracy as time progresses. This goal can be achieved exactly for round spheres. For more general non-trapping obstacles, we show (modulo a technical hypothesis) that certain choices of (a, b) lead to sources that decay exponentially, though perhaps not at the same rate as u itself. As is shown in Section 5, in the smooth, strictly convex case, microlocal analysis indicates that the choice a ≡ 1 should be be optimal.
Our analysis relies upon taking the Fourier transform in the time variable and using recent estimates for the kernel of the CFIE done by Chandler-Wilde and Monk, and Baskin, Spence, and Wunsch, see [1,2]. The results of [1] and our own also rely on several modern improvements to the classical results in scattering theory alluded to above.

Time-harmonic Analysis
Recall that is the outgoing fundamental solution for (∆ + k 2 ), with the corresponding single and double layer potentials defined by: With respect to the complex linear pairing, u, v = ∂D uvdS, the single layer S k is self-dual and If we take the Fourier transform in time: then equation (4) becomes As before the limit is taken from the exterior of D. Here and in the sequel we assume that g(x, t) = 0 for t ≤ 0; by causality µ(x, t) can also be taken to vanish for t ≤ 0. The combined field operators, acting on functions on Γ, given by are well known to be Fredholm operators of second kind, provided that Γ is at least C 1 . For functions a(y) > 0 and b(y) > 0 for all y ∈ Γ, it can be shown, by a simple integration by parts argument, that A(a, b; k) is invertible for k in the closed upper half plane. To prove this we use that the dual of A(a, b; k) with respect to the complex linear pairing is Suppose that for some Im k ≥ 0, there is a non-trivial solution to A(a, b; k)f = 0. By the Fredholm theory this means that there is also a non-trivial solution h to We will have occasion to consider the function u(x) in both the interior of D and its exterior. To make this distinction clear, we will refer to the corresponding restriction of u by u − and u + , respectively. Observe now that u − = 0. If it were, then u + ↾ ∂D = 0 as well. In this case u + would be an solution to (∆ + k 2 )u + = 0 with vanishing Dirichlet data, satisfying the Sommerfeld radiation condition and it would therefore also be zero, which is impossible. Thus, when Im k ≥ 0, this gives a non-trivial solution in D to The boundary condition is equivalent to the fact that A ′ (a, b; k)h = 0.
Using the boundary condition, and integration by parts we obtain: Let k = k 1 + ik 2 , and take the real and imaginary parts to see that this implies The first relation shows that if k 1 = 0 and k 2 > 0, then u − = 0 as well. If k 1 = 0 and k 2 ≥ 0, then the second equation shows that u − = 0. If k 2 = 0, but k 1 = 0, then the first relation implies that u − ↾ bD = 0; the boundary condition in (13) then implies that ∂ n u − = 0 as well. It then follows from Green's formula that u − = 0. This proves the following basic theorem: If a(x) and b(x) are positive functions defined on ∂D, a compact C 1 -surface, then the operators A(a, b; k) and A ′ (a, b; k) are invertible for k with Im k ≥ 0.
Remark 1. If b = 0, then A(a, b; 0) has a non-trivial nullspace. Such a nullspace would generically destroy any possibility for exponential decay in the source function µ(x, t). In Kress' early paper [7] on the combined field operator, he shows that the optimal result (for the disk) when k is close to zero results from taking b = − i 2 + O(k 2 ) in our notation. This is somewhat at odds with our choice to take b real and positive.
Careful examination shows that Kress' choice only works if k 1 ≥ 0 : In the set k 2 < |k 1 | the real part of the quadratic form in (14) is indefinite (regardless of the values that a and b take) and so to obtain the desired result we need to employ the imaginary part of the quadratic form. If b = b 1 + ib 2 , then this form would become In order for the term in the brackets to be definite where k 2 > 0, we need to take a > 0. In order for this expression to be definite in both components of the set 0 < k 2 < |k 1 | it is clearly necessary to take b 2 = 0. Of course one could use functions of k more complicated than ak + ib as the coefficient of S k , but this would considerably complicate the relationship between the solutions in the frequency and the time domains.

Scattering Theory
Let (a(x), b(x)) be positive functions defined on ∂D. In the recent paper of Baskin, Spence, and Wunsch, see Theorem 1.10 of [1], it is shown that there is a positive number β 1 so that the operator (∆ + k 2 )u acting on data in D that satisfies the boundary condition This result does not assume that ∂D is non-trapping. Using this result, along with modern refinements of theorems of Lax and Phillips and Lax, Morawetz and Phillips [8,9], we can prove the following result: Proof. The facts we use, in addition to the result of [1] are 1. If D is a non-trapping region then the generator, B, of the compressed wavesemigroup Z(t) of Lax and Phillips has its spectrum, σ(B), in a half plane of the form Im k ≤ −β 0 < 0.

A number k belongs to σ(B) if and only if there is an "eventually outgoing
By definition (see [9] pp. 126-7) a solution is eventually outgoing (even for Im k < 0) if and only if it can be represented in the form with ν y = −n y , the outer normal to D c . See Theorem IV-4.3 in [9]. Here g k (x) is the "outgoing" fundamental solution for ∆ + k 2 in R 2n+1 .
We argue as before: if for some k, with Im k < 0, the operator A(a, b; k) has a non-trivial null-space, then so does A ′ (a, b; k). Let ϕ = 0 be in this latter nullspace and set u = Unlike the case where Im k ≥ 0, we do not know, a priori, that u − = 0. Indeed, there are just the two cases to consider: u − = 0 or not. If u − = 0, then u − is a non-trivial solution to the boundary value problem considered in [1], and therefore The other possibility is that u − = 0. This implies that u + ↾ ∂D = 0 as well, since u is continuous across ∂D. However u + need not be zero, since it is exponentially growing at infinity. Indeed, if ϕ = 0, and u − = 0, then the jump conditions for ∂ n u ± show that u + cannot be zero. Because u + ↾ ∂D = 0, and ∂ n u + = −ϕ, it follows that u + satisfies (17) and therefore that u + is an eventually outgoing solution in the sense of Lax and Phillips. Hence k must belong to the spectrum of B, and therefore Im k ≤ −β 0 . This shows that if A(a, b; k) is non-invertible, then The converse is also true.

Theorem 3.
If k is a scattering pole for the Dirichlet Laplacian, then A(a, b; k) has a non-trivial null-space.
Proof. The theorem of Lax and Phillips states that k is a scattering pole if and only there is an outgoing solution to the BVP the outgoing condition in this case is equivalent to If we let ϕ = ∂ n u, then S k ϕ vanishes on bD and from which it is immediate that which proves the claim. Thus every scattering resonance occurs among the set {k : ker A(a, b; k) = 0}, for any choice of positive a and b.
Remark 2. These arguments apply generally to identify the frequencies for which the null-space of any combination of single and double layer potentials is nontrivial. The only cases not explicitly covered are those of the single and double layers alone. One easily establishes that S k fails to be invertible for k 2 in the Dirichlet spectrum of D union with the scattering poles of the Dirichlet operator in Ω. The (exterior) double layer Id /2 + D k fails to be invertible for k 2 in the Neumann spectrum of D union with the scattering poles of the Dirichlet operator in Ω.
Spence et al. show that that the resolvent kernel for the interior impedance problem satisfies the norm estimate R a,b (k 2 ) ≤ C 1+|k| for k in a strip around the real axis. It is not clear what the analogous bound is for the analytic continuation of the exterior Dirichlet problem R Dir (k 2 ) into the lower half plane. This operator maps compactly supported data into functions that have exponential growth at infinity. For later applications we will need to determine bounds on A(a, b; k) −1 for k in the lower half plane. In particular it will be important to establish a bound of the form A(a, b; k) −1 ≤ M (1 + |k|) m .
Suppose that g(x, t) is data for the wave equation on the boundary of D with support for t ∈ [0, T ], and let The source term on the boundary as a function of t will be given by the inverse Fourier transform If this is the case, then the corresponding time-domain representation of a solution to the wave-equation is with S t and D t the retarded potential single and double layers for the wave equation.
Assume that for any ǫ > 0 and Im k > −(α − ǫ), we have an estimate like For g with sufficient smoothness in t we can deform the contour to conclude that, for σ < α, we have the representation: and therefore µ(·, t) ≤ Ce −σt g m,T .
For conveniently defined norms, this is essentially what we wanted to prove. Of course this may not give the "optimal" rate of decay. It should be noted that the rate of decay α is directly related to two spectral invariants, one being the genuine scattering resonances, and the other the spectrum of the BVP for ∆ on L 2 (D) with impedance BC At least near the negative imaginary axis, we can push this down into the lower half plane so that the first eigenvalue encountered is a genuine scattering resonance. This follows from a careful examination of the integration by parts argument given above.
Let a and b be positive constants. We assume that (∆u + k 2 u) = 0, and If we integrate by parts we see that, with k = k 1 + ik 2 : If k 1 = 0, the we can substitute from the imaginary part the fact that into the real part to see that Evidently the quantity in the brackets must be positive, which implies that Thus the impedance boundary value problem has a pole free region lying in the disk minus the imaginary axis Along the imaginary axis its clear that k 2 < −b/a is needed for (30) to hold. A somewhat better result can be obtained by using the fact that there is an estimate of the form which implies that, if k 1 = 0, then Taking b/a large, we can arrange to have the pole-free region encompass as much of the lower half plane as we want.

The Round Sphere
In this section we consider the case that Γ is equal to the unit sphere in R 3 . We take a to be a constant α and b a constant β. Since everything commutes with the action of the rotation group, we can analyze this problem one spherical harmonic subspace at a time.

Integral Equations on Spherical Harmonic Subspaces
We need to analyze the action of the operator where τ = t − |x − y|, on data of the form µ(y, t) = Y m n (y)f (t). Notationally it's easier to study the general class of operators of the form: Up to taking time derivatives of f, all our operators take this form.
The key observation is that if we let x ⊥ be a unit vector orthogonal to x, then we can define a coordinate system (θ, φ) on the unit sphere by setting Here, R φ is a rotation through angle φ about x. Substituting, we see that A simple calculation shows that Here P n (z) is the degree n Legendre polynomials normalized so that P n (±1) = (±1) n . Therefore we have that (42) Let s = 2(1 − cos θ), andk(s) = 2(1 − s)k(s), to obtain: To represent G α,β acting on data of this type we need to determine whatk is for each of the terms in (37). The kernels are Acting on Y m n f we get Note that we have replaced ∂ t with −∂ s . Integrating by parts with respect to s gives Denote the 1-d operator in the brackets by G α,β n f. If we let α = β = 1, we get a very simple integral equation of the second kind to solve for f : If n = 0, then this is simply f 0 0 (t) = g 0 0 (t).

The Roles of α and β
The equation is quite informative as regards the rate of decay and regularity of the solution. Assuming that 0 < α < 1, we simplify this equation to obtain If g(t) = 0 for t < 0, then we can let f (t) = 0, for t < 0 as well. The solution to this equation is formally given by the infinite series For data supported in [0, ∞) this sum is finite for any t.
In particular: If g is supported in [0, 2M ] then the sum can be made more explicit: Eventually f satisfies an estimate of the form This shows how the choice of α affects the decay of the source term on the boundary. Furthermore, if g is smooth but its support does not lie in an interval of the form [2l, 2(l + 1)], then the solution f typically has a jump discontinuity at every positive even integer. The n = 0, β = 0 case can be solved by using Laplace transform: where γ α 0 (τ ) = Note, however, that no matter what value α takes, this multiplier has a root at s = 0. This explains why we need to take β > 0 in order to get exponential decay.

Numerical Illustrations
To provide concrete examples of the implications of our analysis, we have carried out representative numerical simulations in the case n = 0. We consider two choices for the Dirichlet data: We also consider three choices for the parameters in the integral equation: (Note that the optimal choice for the sphere, a = b = 1, leads to a trivial equation when n = 0 so we avoid it.) Our numerical method is based on the standard Adams predictor-corrector idea. Recall that for this special case the equation to be solved is Introducing a time step ∆t we rewrite in correction form: In the examples which follow we interpolate past solution data on an interval (j∆t, (j + 1)∆t) by the Lagrange polynomial of degree 5 using the approximate solution at t = k∆t, k = j − 3, . . . , j + 2. Precisely the interpolant is used to calculate µ(t − 2), µ(t + ∆t − 2), t+∆t−2 t−2 µ(s)ds appearing on the right hand side of (57). The future integral uses the 6th order Adams-Bashforth-Moulton predictor-corrector method. That is, we compute a predicted value µ (p) (t + ∆t) using the degree 5 Lagrange interpolant of µ(t−k∆t), k = 0, . . . , 5 to approximate t+∆t t µ(s)ds. We then compute µ(t+∆t) using the interpolant of µ (p) (t+∆t) and µ(t − k∆t), k = 0, . . . , 4 inside the integral. We have verified that convergence at 6th order is generally obtained (in some experiments the observed rate was reduced to 5 when α = β = 0) and we have also cross-checked the results with those computed using 2nd and 4th order solvers constructed the same way. The solutions displayed in the graphs were calculated with a time step ∆t = 97/6400 for the non-oscillatory data and ∆t = 97/12800 for the oscillatory data. Comparisons with a coarsened computation suggest that the solutions are accurate to at least 7 digits in the non-oscillatory case and 4-6 digits when the data was oscillatory. Here the maximum errors in the oscillatory case are two orders of magnitude larger when α = β = 0; looking at the solutions below we attribute the difference to dispersion effects.
Solutions of the integral equation for the choice α = β = 0 are displayed in Figure 1. We see that in the case where the data is oscillatory µ does not decay but oscillates at high frequency. When the data is non-oscillatory we observe linear growth. In each case we note that the solution of the wave equation itself decays exponentially in time, so, as our analysis shows, the long time behavior of the density µ(t) is quite different than that of the function it represents.
Solutions of the integral equation for the choice α = 1, β = 0 are displayed in Figure 2. We see that in the case where the data is oscillatory µ apparently decays in time. However a closer look indicates that it in fact approaches a very small steady state. When the data is non-oscillatory µ clearly approaches an O(1) steady state, again in contrast with the solution of the wave equation itself. This behavior corresponds to the pole at the origin discussed above.  Lastly we consider the solutions which arise when α = 1, β = 1 2 ; these are displayed in Figure 3. Here in each case we have exponential decay of µ in time.

High Frequency Asymptotics for Convex Bodies
In the previous section we saw that, at least for round spheres, a principal determinant of the decay rate for the solution of the integral equations on the boundary, ]. Choosing α = 1 removes the sharp contribution from the antipodal point, and gives the optimal rate of decay available for this case. In this section we study the high frequency behavior of D k − i(ka + ib)S k for a smooth surface Γ. In general one cannot expect to have the sort of exact cancellation attained for a round sphere, however we demonstrate the rather remarkable fact that for convex bodies, the leading order part of the delay term is optimally canceled by taking a = 1. It is less clear how to choose b, but this is related to insuring invertibility of A(a, b; k) at zero frequency.
We begin by recalling the formulae for the relevant kernels: As is well known, these are weakly singular integral operators on a C 1 surface in R 3 .
Since we are interested in the high frequency asymptotics, the contribution of the diagonal singularity is of less interest to us than the that of the other critical points of the phase function φ x (y) = |x − y|.
For completeness we state the contribution from the diagonal: where H(x) = 1 2 (κ 1 (x) + κ 2 (x)) is the mean curvature of Γ at x with respect to the outer normal vector n y . Observe that the diagonal asymptotics of the combined field operator are given by: Now we turn to the other critical points of φ x : Even for a strictly convex body, these critical points can be degenerate. For example, if Γ contains an open subset, U, of a sphere, whose center x 0 lies on Γ, then U ⊂ C x 0 . To use the standard techniques of stationary phase we restrict our attention to the case of where Γ is C ∞ and C x consists entirely of non-degenerate critical points. There are, for example, many boundaries for which C x consists of a global maximum for every x ∈ Γ. A simple analysis of this property, using the characterization of C x in (62), shows that it is stable under C 2 small perturbations. If x is a non-degenerate critical point of φ x , then x lies along the line { x+tn x }. We choose orthogonal coordinates so that x = 0, n x = (0, 0, 1), so that x = (0, 0, ±d), with d = |x − x|. For a convex body we always have −d. In these coordinates we represent the surface Γ near to 0 as a graph over its tangent plane by taking y = (z, h(z)) with h(0) = ∇h(0) = 0.
Let ψ be a smooth function equal to 1 in a neighborhood of 0, with a single critical point of φ x in its support. With x = (0, 0, ±d), the asymptotic contribution of x to S k f (x), which we denote by S x k f (x), is given by To apply the stationary phase formula (at least to leading order) we need only compute the Hessian of the phase function at 0; it is H x x (0) = 1 d 1 ∓ dh 11 (0) ∓dh 12 (0) ∓dh 21 (0) 1 ∓ dh 22 (0) .