A spectrally accurate direct solution technique for frequency-domain scattering problems with variable media

A. Gillman, A. H. Barnett, and P. G. Martinsson

Department of Mathematics, Dartmouth College Department of Applied Mathematics, University of Colorado at Boulder

Abstract: This paper presents a direct solution technique for the scattering of time-harmonic waves from a bounded region of the plane in which the wavenumber varies smoothly in space. The method constructs the interior Dirichlet-to-Neumann (DtN) map for the bounded region via bottom-up recursive merges of (discretization of) certain boundary operators on a quadtree of boxes. These operators take the form of impedance-to-impedance (ItI) maps. Since ItI maps are unitary, this formulation is inherently numerically stable, and is immune to problems of artificial internal resonances. The ItI maps on the smallest (leaf) boxes are built by spectral collocation on tensor-product grids of Chebyshev nodes. At the top level the DtN map is recovered from the ItI map and coupled to a boundary integral formulation of the free space exterior problem, to give a provably second kind equation. Numerical results indicate that the scheme can solve challenging problems wavelengths on a side to 9-digit accuracy with 4 million unknowns, in under 5 minutes on a desktop workstation. Each additional solve corresponding to a different incident wave (right-hand side) then requires only seconds.

## 1. Introduction

### 1.1. Problem formulation

Consider time-harmonic waves propagating in a medium where the wave speed varies smoothly, but is constant outside of a bounded domain . This manuscript presents a technique for numerically solving the scattering problem in such a medium. Specifically, we seek to compute the scattered wave that results when a given incident wave (which satisfies the free space Helmholtz equation) impinges upon the region with variable wave speed, as in Figure 1. Mathematically, the scattered field satisfies the variable coefficient Helmholtz equation

(1) |

and the outgoing Sommerfeld radiation condition

(2) |

uniformly in angle. The real number in (1) and (2) is the free space wavenumber (or frequency), and the so called “scattering potential” is a given smooth function that specifies how the wave speed (phase velocity) at the point deviates from the free space wave speed ,

(3) |

One may interpret as a spatially-varying refractive index. Observe that is identically zero outside . Together, equations (1) and (2) completely specify the problem. When is real and positive for all , the problem is known to have a unique solution for each positive [10, Thm. 8.7].

The transmission problem (1)-(2), and its generalizations, have applications in acoustics, electromagnetics, optics, and quantum mechanics. Some specific applications include underwater acoustics [3], ultrasound and microwave tomography [14, 35], wave propagation in metamaterials and photonic crystals, and seismology [36]. The solution technique in this paper is high-order accurate, robust and computationally highly efficient. It is based on a direct (as opposed to iterative) solver, and thus is particularly effective when the response of a given potential to multiple incident waves is desired, as arises in optical device characterization, or computing radar scattering cross-sections. The complexity of the method is where is the number of discretization points in . Additional solves with the same scattering potential and wavenumber require only operations. (Further reductions in asymptotic complexity can sometimes be attained; see section 4.2.) For simplicity of presentation, the solution technique is presented in ; however, the method can be directly extended to .

### 1.2. Outline of proposed method

We solve (1)-(2) by splitting the problem into two parts, namely a variable-coefficient problem on the bounded domain , and a constant coefficient problem on the exterior domain . For each of the two domains, a solution operator in the form of a Dirichlet-to-Neumann (DtN) map on the boundary is constructed. These operators are then “glued” together on to form a solution operator for the full problem. The end result is a discretized boundary integral operator that takes as input the restriction of the incoming field (and its normal derivative) to , and constructs the restriction of the scattered field on (and its normal derivative). Once these quantities are known, the total field can rapidly be computed at any point .

#### Solution technique for the variable-coefficient problem on

For the interior domain , we construct a solution operator for the following homogeneous variable-coefficient boundary value problem where the unknown total field satisfies

(5) | ||||||

(6) |

Note that, for now, we specify the Dirichlet data ; when we consider the full problem will become an unknown that will also be solved for. It is known that, for all but a countable set of wavenumbers, the interior Dirichlet BVP (5)-(6) has a unique solution [26, Thm. 4.10]. The values at which the solution is not unique are (the square roots of) the interior Dirichlet eigenvalues of the penetrable domain ; we will call them resonant wavenumbers.

###### Definition 1.1 (interior Dirichlet-to-Neumann map).

###### Remark 1.2.

In this paper, a discrete approximation to is constructed via a variation of recent composite spectral methods in [16, 25] (which are also similar to [9]). These methods partition into a collection of small “leaf” boxes and construct approximate DtN operators for each box via a brute force calculation on a local spectral grid. The DtN operator for is then constructed via a hierarchical merge process. Unfortunately, at any given , each of the many leaves and merging subdomains may hit a resonance as described above, causing its DtN to fail to exist. As approaches any such resonance the norm of the DtN grows without bound. Thus a technique based on the DtN alone is not robust.

###### Remark 1.3.

We remind the reader that any such “box” resonance is purely artificial and is caused by the introduction of the solution regions. It is important to distinguish these from resonances that the physical scattering problem (1)-(2) itself might possess (e.g. due to nearly trapped rays), whose effect of course cannot be avoided in any accurate numerical method.

One contribution of the present work is to present
a robust improvement to the methods of [16, 25],
built upon hierarchical merges of impedance-to-impedance (ItI) rather than
DtN operators; see section 2.
The idea of using impedance coupling builds upon the work of
Kirsch–Monk [22].
ItI operators are inherently stable, with condition number ,
and thus exclude the
possibility of inverting arbitrarily ill-conditioned matrices as in
the original DtN formulation.
For instance, in the *lens* experiment of section 5,
the DtN method has condition numbers as large as , while for the new ItI
method the condition number never grows larger than .
The DtN of the whole domain is still needed; however, if has a resonance,
the size of can be changed slightly to avoid the resonance (see remark 2.6).

#### Solution technique for the constant coefficient problem on

On the exterior domain, we build a solution operator for the constant coefficient problem

(8) | |||||

(9) | |||||

(10) |

obtained by restricting (1)-(2) to . (Again, the Dirichlet data later will become an unknown that is solved for.) It is known that (8)-(10) has a unique solution for every wavenumber [10, Ch. 3]. This means that the following DtN for the exterior domain is always well-defined.

###### Definition 1.2 (exterior Dirichlet-to-Neumann map).

#### Combining the two solution operators

Once the DtN operators and have been determined (as described in sections 1.2.1 and 1.2.2), and the restriction of the incident field to is given, it is possible to determine the scattered field on as follows. First observe that the total field must satisfy

(12) |

We also know that the scattered field satisfies

(13) |

Combining (12) and (13) to eliminate results in the equation (analogous to [22, Eq. (2.12)]),

(14) |

As discussed in Remark 1.2, both and have order . Lamentably, this behavior adds rather than cancels in (14), so that also has order , and is therefore unbounded. This makes any numerical discretization of (14) ill-conditioned, with condition number growing linearly in the number of boundary nodes. To remedy this, we present in section 3 a new method for combining and to give a provably second kind integral equation, which thus gives a well-conditioned linear system.

Once the scattered field is known on the boundary, the field at any exterior point may be found via Green’s representation formula; see section 3.3. The interior transmitted wave may be reconstructed anywhere in by applying solution operators which were built as part of the composite spectral method.

### 1.3. Prior work

Perhaps the most common technique for solving the scattering problem stated in Section 1.1 is to discretize the variable coefficient PDE (1) via a finite element or finite difference method, while approximating the radiation condition in one of many ways, including perfectly matched layers (PML) [20], absorbing boundary conditions (ABC) [12], separation of variables or their perturbations [27], local impedance conditions [7], or a Nyström method [22] (as in the present work). However, the accuracy of finite element and finite difference schemes for the Helmholtz equation is limited by so-called “pollution” (dispersion) error [2, 4], demanding an increasing number of degrees of freedom per wavelength in order to maintain fixed accuracy as wavenumber grows. In addition, while the resulting linear system is sparse, it is also large and is often ill-conditioned in such a way that standard pre-conditioning techniques fail, although hybrid direct-iterative solvers such as [13] have proven effective in certain environments. While there do exist fast direct solvers for such linear systems (for low wavenumbers ) [32, 31, 37, 24], the accuracy of the solution is limited by the discretization. The performance of the solver worsens when increasing the order of the discretization—thus it is not feasible to use a high order discretization that would overcome the above-mentioned pollution error.

Scattering problems on infinite domains are also commonly handled by rewriting them as volume integral equations (e.g. the Lippmann–Schwinger equation) defined on a domain (such as ) that contains the support of the scattering potential [1, 8]. This approach is appealing in that the Sommerfeld condition (2) is enforced analytically, and in that high-order discretizations can be implemented without loss of stability [11]. Principal drawbacks are that the resulting linear systems have dense coefficient matrices, and tend to be challenging to solve using iterative solvers [11].

### 1.4. Outline

Section 2 describes in detail the stable hierarchical procedure for constructing an approximation to the DtN map for the interior problem (5)-(6). Section 3 describes how boundary integral equation techniques are used to approximate the DtN map for the exterior problem (8)-(10), how to couple the DtN maps and to solve the full problem (1)-(2), and the proof (Theorem 3.1) that the formulation is second kind. Section 4 details the computational cost of the method and explains the reduced cost for multiple incident waves. Finally, section 5 illustrates the performance of the method in several challenging scattering potential configurations.

## 2. Constructing and merging impedance-to-impedance maps

This section describes a technique for building a discrete approximation to the Dirichlet-to-Neumann (DtN) operator for the interior variable coefficient BVP (5)-(6) on a square domain . It relies on the hierarchical construction of impedance-to-impedance (ItI) maps; these are defined in section 2.1. Section 2.2 defines a hierarchical tree on the domain . Section 2.3 explains how the ItI maps are built on the (small) leaf boxes in the tree. Section 2.4 describes the merge procedure whereby the global ItI map is built, and then how the global DtN map is recovered from the global ItI map.

### 2.1. The impedance-to-impedance map

We start by defining the ItI map on a general Lipschitz domain, and giving some of its properties. (In this section only, will refer to such a general domain.)

###### Proposition 2.1.

Let be a bounded Lipschitz domain, and be real. Let , and . Then the interior Robin BVP,

(15) | |||||

(16) |

has a unique solution for all real .

###### Proof.

We first prove uniqueness. Consider a solution to the homogeneous problem . Then using Green’s 1st identity and (15), (16),

Taking the imaginary part shows that , and hence , vanishes on , hence in by unique continuation of the Cauchy data. Existence of now follows for data from the Fredholm alternative, as explained in this context by McLean [26, Thm. 4.11]. ∎

###### Definition 2.1 (interior impedance-to-impedance map).

We choose the impedance parameter (on dimensional grounds) to be . Numerically, in what follows, we observe very little sensitivity to the exact value or sign of .

For the following, we need the result that the DtN map is self-adjoint for real and . This holds since, for any functions and satisfying in and in ,

by Green’s second identity. This allows us to prove the following property of the ItI map that will be the key to the numerical stability of the method.

###### Proposition 2.2.

Let be a bounded Lipschitz domain, let be real, and let and . Then the ItI map for exists for all real frequencies , and is unitary whenever is also real.

###### Proof.

Existence of for all real follows from Proposition 2.1. To prove is unitary, we insert the definitions of and into (19) and use the definition of the DtN to rewrite , giving

(20) |

which holds for any data . Thus the ItI map is given in operator form by

(21) |

Since is self-adjoint and is real, this formula shows that is unitary. ∎

As a unitary operator, has unit operator -norm, pseudo-differential order , and eigenvalues lying on the unit circle. From (21) and the pseudo-differential order of one may see that the eigenvalues of accumulate only at .

### 2.2. Partitioning of into hierarchical tree of boxes

Recall that is the square domain containing the support of . We partition into a collection of equally-sized square boxes called leaf boxes, where sets the number of levels; see Figure 2. We place Gauss–Legendre interpolation nodes on each edge of each leaf, which will serve to discretize all interactions of this leaf with its neighbors; see Figure 3(a). The size of the leaf boxes, and the parameter , should be chosen so that any potential transmitted wave , as well as its first derivatives, can be accurately interpolated on each box edge from their values at these nodes.

Next we construct a binary tree over the collection of leaf boxes. This is achieved by forming the union of adjacent pairs boxes (forming rectangular boxes), then forming the pairwise union of the rectangular boxes. The result is a collection of squares with twice the side length of a leaf box. The process is continued until the only box is itself, as in Figure 2. The boxes should be ordered so that if is a parent of a box , then . We also assume that the root of the tree (i.e. the full box ) has index . Let denote the domain associated with box .

###### Remark 2.1.

The method easily generalizes to rectangular boxes, and to more complicated domains in the same manner as [25].

### 2.3. Spectral approximation of the ItI map on a leaf box

Let denote a single leaf box, and let and be a pair of vectors of associated incoming and outgoing impedance data, sampled at the Gauss–Legendre boundary nodes, with entries ordered in a counter-clockwise fashion starting from the leftmost node of the bottom edge of the box, as in Figure 3(a). In this section, we describe a technique for constructing a matrix approximation to the ItI operator on this leaf box. Namely, we build a matrix such that

holds to high-order accuracy, for all incoming data vectors corresponding to smooth transmitted wave solutions .

First, we discretize the PDE (15) on the square leaf box using a spectral method on a tensor product Chebyshev grid filling the box, as in Figure 3(b), comprised of the nodes at locations

where , are the Chebyshev points on . We label the Chebyshev node locations , for . For notational purposes, we order these nodes in the following fashion: the indices correspond to the Chebyshev nodes lying on the boundary of , ordered counter-clockwise starting from the node located at the south-west corner . The remaining interior nodes have indices and may be ordered arbitrarily (a Cartesian ordering is convenient).

Let be the standard spectral differentiation matrices constructed on the full set of Chebyshev nodes, which approximate the (horizontal) and (vertical) derivative operators, respectively. As explained in [33, Ch. 7], these are constructed from Kronecker products of the identity matrix and , where is the usual one-dimensional differentiation matrix on the nodes . Namely has entries where is the vector of barycentric weights for the Chebyshev nodes (see [33, Ch. 6] and [30, Eqn.(8)].) Let the matrix be the spectral discretization of the operator on the product Chebyshev grid, namely

where “diag S” indicates the diagonal matrix whose entries are the elements of the ordered set .

###### Remark 2.2.

The matrices , , and must have rows and columns ordered as explained above (i.e. boundary then interior) for the Chebyshev nodes; this requires permuting rows and columns of the matrices constructed by Kronecker products. For example, the structure of is

where the notation indicates the submatrix block , etc.

We now break the boundary Chebyshev nodes into four sets , denoting the south, east, north, and west edges, as in Figure 3(b). Note that includes the south-western corner but not the south-eastern corner (which in turn is the first element of ), etc.

We are now ready to derive the linear system required for constructing the approximate ItI operator. We first build a matrix which maps values of at all Chebyshev nodes to the outgoing normal derivatives at the boundary Chebyshev nodes, as follows,

(22) |

where (as is standard in MATLAB) the notation denotes the matrix formed from the subset of rows of a matrix given by the index set . Then, recalling (16), the matrix which maps the values of at all Chebyshev nodes to incoming impedance data on the boundary Chebyshev nodes is

(23) |

where denotes the identity matrix of size . Using for the vector of values at all Chebyshev nodes, the linear system for the unknown that imposes the spectral discretization of the PDE at all interior nodes, and incoming impedance data at the boundary Chebyshev nodes, is

(24) |

where is an appropriate column vector of zeros, and the square size- system matrix is

###### Remark 2.3.

At each of the four corner nodes, only one boundary condition is imposed, namely the one associated with the edge lying in the counter-clockwise direction. This results in a square linear system, which we observe is around twice as fast to solve as a similar-sized rectangular one.

To construct the “solution matrix” for the linear system, we solve (24) for each unit vector in , namely

In practice, is found using the backwards-stable solver available via MATLAB’s mldivide command. If desired, the tabulated solution can now be found at all the Chebyshev nodes by applying to the right hand side of (24).

Recall that the goal is to construct matrices that act on boundary data on Gauss (as opposed to Chebyshev) nodes. With this in mind, let be the matrix which performs Lagrange polynomial interpolation [34, Ch. 12] from Gauss to Chebyshev points on a single edge, and let be the matrix from Chebyshev to Gauss points. Let be with the last row omitted. For example, maps from Gauss points on the south edge to the Chebyshev points .

Then the solution matrix which takes incoming impedance data on Gauss nodes to the values at all Chebyshev nodes is

(25) |

Finally, we must extract outgoing impedance data on Gauss nodes
from the vector ,
to construct an approximation to the full ItI map on the Gauss nodes.
This is done by extracting (as in (22)) the relevant rows of the spectral differentiation matrices,
then interpolating back to Gauss points.
Let be the indices of all Chebyshev nodes on the
south edge, and correspondingly for the other three edges.
Then the index set
counts each corner twice.^{2}

Then, in terms of (25),

is the desired spectral approximation to the ItI map on the leaf box.

The computation time is dominated by the solution step for , which takes effort . We observe empirically that one must choose in order that not acquire a spurious numerical null space. We typically choose and .

### 2.4. Merging ItI maps

Once the approximate ItI maps are constructed on the boundary Gauss nodes on the leaf boxes, the ItI map defined on is constructed by merging two boxes at a time, moving up the binary tree as described in section 2.2. This section demonstrates the purely local construction of an ItI operator for a box from the ItI operators of its children.

We begin by introducing some notation. Let denote a box with children and where . For concreteness, consider the case where and share a vertical edge. As shown in Figure 4, the Gauss points on and are partitioned into three sets: : Boundary nodes of that are not boundary nodes of . : Boundary nodes of that are not boundary nodes of . : Boundary nodes of both and that are not boundary nodes of the union box . Define interior and exterior outgoing data via

The incoming vectors and are defined similarly. The goal is to obtain an equation mapping to . Since the ItI operators for and have previously been constructed, we have the following two equations

(26) |

Since the normals of the two leaf boxes are opposed on the interior “3” edge, and (using the definitions (17)-(18)). This allows the bottom row equations to be rewritten using only on the interior edge, namely

(27) |

and

(28) |

Let . Plugging (27) into (28) results in the following equation

By collecting like terms and solving for , we find

(29) |

Note that the matrix maps the incoming impedance data on to the incoming (with respect to ) impedance data on the interior edge. By combining the relationship between the impedance boundary data on neighbor boxes and (27), the matrix computes the impedance data .

The outgoing impedance data is found by plugging into equation (27). Now the top row equations of (26) can be rewritten without reference to the interior edge. The top row equation of (26) from is now

and the top row equation of (26) from is

Writing these equations as a system, we find and satisfy

(30) |

Thus the block matrix on the left hand side of (30) is , the ItI operator for .

###### Remark 2.4.

In practice, the matrix products and should be computed once per merge.

###### Remark 2.5.

Algorithm 1 outlines the implementation of the hierarchical construction of the impedance operators, by repeated application of the above merge operation. Algorithm 2 illustrates the downwards sweep to construct from incoming impedance data the solution at all Chebyshev discretization points in . Note that the latter requires the solution matrices at each level, and for each of the leaf boxes, that were precomputed by Algorithm 1.

The resulting approximation to the top-level ItI operator is a square matrix which acts on incoming impedance data living on the total of composite Gauss nodes on . An approximation to the interior DtN map on these same nodes now comes from inverting equation (20) for , to give

(31) |

where is the identity matrix of size . The need for conversion from the ItI to the DtN for the domain will become clear in the next section.

###### Remark 2.6.

Due to the pseudo-differential order of , we expect the norm of to grow linearly in the number of boundary nodes. However, it is also possible that falls close enough to a resonant wavenumber of that the norm of is actually much larger, resulting in a loss of accuracy due to the inversion of the nearly-singular matrix in (31). In our extensive numerical experiments, this latter problem has never happened. However, it is important to include a condition number check when formula (31) is evaluated. Should there be a problem, it would be a simple matter to modify slightly the domain to avoid a resonance. For instance, one can add a column of leaf boxes to the side of , and then inexpensively update the computed ItI operator for to the enlarged domain.

## 3. Well-conditioned boundary integral formulation for scattering

In this section, we present an improved boundary integral equation alternative to the scattering formulation (14) from the introduction.

### 3.1. Formula for the exterior DtN operator

We first construct the exterior DtN operator via potential theory. The starting point is Green’s exterior representation formula [10, Thm. 2.5], which states that any radiative solution to the Helmholtz equation in may be written,

(32) |

where and are respectively the frequency- Helmholtz double- and single-layer potentials with density , with the outgoing Hankel function of order zero. , The vector denotes the outward unit normal at ,

Letting approach in (32), one finds via the jump relations,

(33) |

where and are the double- and single-layer boundary integral operators on . See [10, Ch. 3.1] for an introduction to these representations and operators. Rearranging (33) gives , thus the exterior DtN operator is given in terms of the operators of potential theory by

(34) |

### 3.2. The new integral formulation

We apply from the left the single layer integral operator to both sides of (14), and use (34), to obtain

(35) |

a linear equation for , the restriction of the scattered wave to the domain boundary . Let

(36) |

be the boundary integral operator appearing in the above formulation. In the trivial case (no scattering potential) it is easy to check that , by using which can be derived in this case similarly to (34). Now we prove that introducing a general scattering potential perturbs only compactly, that is, our left-regularization of the original ill-conditioned (14) has produced a well-conditioned equation.

###### Theorem 3.1.

###### Proof.

Let satisfy (5) in , then by Green’s interior representation formula (third identity) [10, Eq. (2.4)],

(37) |

where denotes the Helmholtz volume potential [10, Sec. 8.2] acting on a function with support in . Define to be the solution operator for the interior Dirichlet problem (5)-(6), i.e. for , and let denote the operator that multiplies a function pointwise by , we can write the last term in (37) as . Taking to from inside in (37), the jump relations give

(38) |

where is restricted to evaluation on . Recall . Plugging this definition into (38), we find

When substituted into (36) results in cancellation of the terms, giving

Now is bounded [