## Abstract

The scattering of acoustic waves by irregular structures plays an important role in a wide range of problems of scientific and technological interest. This contribution focuses on the rapid and highly accurate numerical approximation of solutions of Helmholtz equations coupled across irregular periodic interfaces meant to model acoustic waves incident upon a multi-layered medium. We describe not only a novel surface formulation for the problem in terms of boundary integral operators (Dirichlet–Neumann operators), but also a Boundary Perturbation methodology (the Method of Operator Expansions) for its numerical simulation. The method requires only the discretization of the layer interfaces (so that the number of unknowns is an order of magnitude smaller than volumetric approaches), while it avoids not only the need for specialized quadrature rules but also the dense linear systems characteristic of Boundary Integral/Element Methods. The approach is a generalization to multiple layers of Malcolm & Nicholls' Operator Expansions algorithm for dielectric structures with two layers. As with this precursor, this approach is efficient and spectrally accurate.

## 1. Introduction

The interaction of acoustic waves with irregular structures plays an important role in a wide range of problems of scientific and technological interest. From remote sensing (Tsang *et al.* 1985) to underwater acoustics (Brekhovskikh & Lysanov 1982), the ability to simulate, in a robust and accurate way, the fields generated by such structures is of crucial importance to researchers from many disciplines. This contribution focuses upon the rapid and highly accurate numerical approximation of solutions of Helmholtz equations coupled across irregular periodic interfaces meant to model acoustic waves incident upon a multi-layered medium. We describe not only a novel surface formulation for the problem in terms of boundary integral operators (Dirichlet–Neumann operators, DNOs), but also a Boundary Perturbation methodology (the Method of Operator Expansions) for their numerical simulation.

A wide array of numerical algorithms have been devised for the simulation of problems akin to the one we consider here. The classical Finite Difference (Pratt 1990), Finite Element (Zienkiewicz 1977) and Spectral Element (Komatitsch & Tromp 2002) methods are available but suffer from the requirement that they discretize the full volume of the problem domain which not only introduces a huge number of degrees of freedom, but also raises the difficult question of appropriately specifying a far-field boundary condition explicitly.

The surface integral methods are popular and appealing alternatives (Sanchez-Sesma *et al.* 1989; e.g. Boundary Integral Methods—BIM—or Boundary Element Methods—BEM), which only require a discretization of the layer *interfaces* (rather than the whole structure) and which, owing to the choice of the Green function, satisfy the far-field boundary condition exactly. While these methods can deliver high-accuracy simulations with greatly reduced operation counts, there are several difficulties which need to be addressed (Reitich & Tamma 2004). First, high-order simulations can only be realized with specially designed quadrature rules which respect the singularities in the Green function (and its derivative, in certain formulations). Additionally, BIM/BEM typically give rise to dense linear systems to be solved which require carefully designed preconditioned iterative methods (with accelerated matrix–vector products; e.g. by the Fast–Multipole Method; Greengard & Rokhlin 1987) for configurations of engineering interest.

Boundary Perturbation Methods (BPMs) have recently received attention as an alternative strategy that maintain the reduced numbers of degrees of freedom of BIM/BEM while avoiding the need for special quadrature formulas or preconditioned iterative solution procedures for dense systems. Bruno & Reitich introduced the Method of Field Expansions (FE) (Bruno & Reitich 1993*a*,*b*,*c*) for doubly layered media, and (Malcolm & Nicholls 2011*b*) generalized this to an arbitrary number of layers, greatly improving upon the prohibitive operation counts of Dinesen & Hesthaven's (2000) extension based upon iterated two-layer solvers.

A closely related BPM was devised by Milder (1991*a*,*b*) (see also the improvements of Coifman *et al.* 1999) for the simulation of scattering by impenetrable gratings. This Operator Expansions (OE) approach was recently generalized by Malcolm & Nicholls (2011*a*) in the case of doubly layered media for the purpose of devising an algorithm for the inverse problem of identifying the interface shape based upon far-field measurements. This approach enjoys all of the speed and accuracy of Bruno & Reitich's FE method, while only requiring *one* surface unknown rather than the *two*, which the implementation of Malcolm & Nicholls (2011*b*) currently requires.

In this contribution, we generalize this OE method of Milder and Malcolm & Nicholls to the case of (*M*+1)-many layers. In this formulation, we realize the *minimum* number of problem unknowns (*M* surface functions, also used by BIM/BEM) which halves the number currently mandated by the approach of Malcolm & Nicholls (2011*b*). While we fully anticipate that the FE algorithm of Malcolm & Nicholls (2011*b*) can be reformulated to involve only one unknown per interface, our new approach has the additional benefit of having rather *explicit* dependence upon the interface shapes, which make their application towards an inverse problem algorithm very appealing (cf. Malcolm & Nicholls 2011*a*). However, we leave this latter observation for future work and focus instead upon this novel forward solver.

The organization of the paper is as follows. In §2, we recall the governing equations of acoustic scattering in a triply layered medium, and in §2*a*, we present a novel boundary formulation of the problem in terms of DNOs. In §§2*b*–*d*, we outline an OE methodology for approximating solutions to our new boundary formulation. In §3, we repeat these developments in the general case of (*M*+1) layers. We outline a boundary formulation in §3*a* with an OE implementation in §3*b*. Finally, in §4, we provide details of our numerical experiments. In §4*a*, we discuss a class of exact solutions and in §4*b*, we give details of our numerical approach, including error measurement in §4*c* and Padé approximation in §4*d*. In §§4*e*,*f*, we display results in two and three dimensions, respectively, and brief concluding remarks are given in §5.

## 2. Governing equations: three layers

For ease of notation, we begin our developments with a detailed description of scattering by a triply layered media. We will return to the general setting in the following section. It is well-known (Petit 1980; Ihlenburg 1998) that the (reduced) scattered pressure inside a periodic structure satisfies the Helmholtz equation with isonification conditions at the interface, and outgoing wave conditions at positive and negative infinity. More precisely, for a triply layered medium, we define the domains
with *d*_{x}×*d*_{y} periodic interfaces *g* and *h*,
and (upward pointing) normals
All three domains are constant-density acoustic media with velocities *c*_{j} (*j*=*u*,*v*,*w*); we assume that plane-wave radiation is incident upon the structure from above:
2.1
With these specifications, we can define in each layer the parameter *k*_{j}=*ω*/*c*_{j} which characterizes both the properties of the material and the frequency of radiation in the structure. If the reduced scattered fields (i.e. the full scattered fields with the periodic time dependence factored out) in *S*_{u}, *S*_{v} and *S*_{w} are, respectively, denoted {*u*,*v*,*w*}={*u*(*x*,*y*,*z*),*v*(*x*,*y*,*z*),*w*(*x*,*y*,*z*)}, then these functions will be quasiperiodic (Petit 1980)
and the system of partial differential equations (PDEs) to be solved are
2.2a
2.2b
2.2c
2.2d
2.2e
2.2f
and
2.2g
where
2.2h
and
2.2i
If continuity is enforced inside the structure, then *θ*=*μ*≡0; however, as we shall see, it is no impediment to the method if we set it to any non-zero function. In these equations, the operator enforces the condition that scattered solutions must either be ‘outgoing’ (upward in *S*_{u} and downward in *S*_{w}) if they are propagating, or ‘decaying’ if they are evanescent. We make this ‘Outgoing Wave Condition’ (OWC; Petit 1980) more precise in the Fourier series expression for the exact solution (see (2.3) below).

The quasiperiodic solutions of the Helmholtz equations—(2.2*a*, *c* and *e*)—and the OWCs—(2.2*b*,*f*)—are given by Petit (1980)
2.3a
2.3b
and
2.3c
where
*j*=*u*,*v*,*w*. The OWC mandates that we choose the positive sign in front of *γ*_{u,p,q} in (2.3a) and the negative sign in front of *γ*_{w,p,q} in (2.3c). These formulas are valid provided that (*x*,*y*,*z*) are *outside* the grooves, i.e.

### (a) Boundary formulation

The key to our solver is the realization that if we recover the Dirichlet and Neumann traces of *u*, *v* and *w*,
then integral formulas will tell us *u*, *v* and *w* everywhere. We point out that the choice of signs on prescribe normal derivatives *exterior* to the domain of definition of the corresponding fields {*u*,*v*,*w*}. Furthermore, if we define the DNOs
2.4a
2.4b
and
2.4c
then it suffices to simply find the Dirichlet traces {*U*,*V* ^{g},*V* ^{h},*W*}. As the DNOs encapsulate the solution of the Helmholtz equations and the OWCs, it is not difficult to see that (2.2) are equivalent to the *surface* equations
We can further simplify by using the first and third equations to express
and then insert these into the second and fourth equations:
or
2.5

### Remark 2.1

At this point, we indicate that the formulation (2.5) is independent of algorithm choice and that the method we advocate below (the Method of Operator Expansions) is not the only approach which can be used to solve these equations. For instance, a surface integral approach could also be considered. The crucial consideration is how are the operators *H*^{gg}, *H*^{hg} and *J* to be computed and, more importantly, how is the inverse of the operator
to be approximated?

We choose the OE approach because of its favourable operation counts (which match integral equation methods), rapid execution times and remarkable accuracy (please see Bruno & Reitich (1993*a*,*b*,*c*) and (Milder 1991*a*,*b*) for a full verification of all of these claims). We note, however, that for very large and/or rough interfaces, the OE algorithm does face conditioning challenges. We refer the interested reader to the work of the author and F. Reitich for a comprehensive discussion of the relevant height-to-period ratios, smoothness requirements and frequency dependence of this phenomena (Nicholls & Reitich 2004*a*), together with potential remedies (Nicholls & Reitich 2004*b*).

### (b) An operator expansions method

We propose a perturbative approach to the solution of (2.5) based upon the assumption *g*(*x*,*y*)=*εf*(*x*,*y*) and *h*(*x*,*y*)=*εs*(*x*,*y*) where, *a priori*, *ε* is assumed small. If this is the case, then it can be shown that the data {*ζ*,*ψ*,*θ*,*μ*} and operators {*G*,*H*^{gg},*H*^{gh},*H*^{hg},*H*^{hh},*J*} depend analytically upon *ε* so that
Furthermore, the scattered fields can also be shown to be analytical so that
Inserting these into (2.5), we see that
At order we find
2.6
where
and
while at order
where
and
Solving for ,
2.7
Note that at *every* perturbation order in this approach, we repeatedly invert the same operator, *A*_{0}, which, as we shall see, is *block diagonal* in Fourier space and can, therefore, be accomplished very rapidly.

### (c) Expansions: surface data and the Dirichlet–Neumann operators *G* and *J*

The key to our algorithm is convenient, high-order formulas for the functions {*ζ*_{n},*ψ*_{n},*θ*_{n},*μ*_{n}}, and the operators . While many of these formulas have appeared in the literature (particularly for two-dimensional configurations), we collect all of them here for convenience.

We begin with *ζ*
where *F*_{n}(*x*,*y*):=*f*(*x*,*y*)^{n}/*n*! Thus
2.8
Similarly, for *ψ* we have
So
2.9a
and
2.9b
for *n*>0. In the current model, we have included *θ* and *μ* as forcing terms at the lower interface as they introduce no significant complications to the algorithm. However, for our present purposes, we simply enforce continuity at this interface and therefore set *θ*=*μ*≡0 so that *θ*_{n}=*μ*_{n}=0 for all *n*.

The operators
are a bit more involved and we use the method of OE (Milder 1991*a*; Craig & Sulem 1993; Nicholls & Reitich 2004*a*) to find the action of
In previous work (Nicholls & Reitich 2004*a*,*b*), it was shown that
where we have used Fourier multiplier notation
Furthermore,
2.10
In an exactly analogous fashion, we can derive
and
2.11

### (d) Expansions: the operator *H*

The matrix-valued operator *H* is more complicated than the DNOs *G* and *J* presented above; however, the same OE philosophy can be brought to bear upon this operator as well. To begin, we define the function
2.12
where
which satisfy Helmholtz's equation in the middle layer and
since
We insert this into the definition of *H*
or
2.13
With an eye towards simplifying these expressions, we define **S**=**S**^{g}+**S**^{h}, where
2.14
and
2.15
and **C**=**C**^{g}+**C**^{h}, where
2.16
and
2.17
We note that
and
so that, by factoring out (i*γ*_{v,p,q}), we realize an operator **C**(0) which is bounded as the wavenumbers increase. With these definitions, equation (2.13) reads
2.18

At this point, we see a further difference with the operator *H*: it depends not only on two Dirichlet data, but also two boundaries. While we view these two boundaries as of the same order of magnitude, *g*(*x*)=*εf*(*x*) and *h*(*x*)=*εs*(*x*), where both , we note that this need not be the case and, in fact, the variations can be quite independent. Expanding the operators *H*, **S** and **C** in Taylor series delivers
2.19

#### (i) Zeroth order

At order in (2.19), we find Noting that and we have 2.20

#### (ii) General domain

At order , *n*>0 in (2.19), we have
Using and, with a simplification in mind, introducing two cancelling factors of (i*γ*_{v,p,q}), we find
Using , we have
or, since,
for a general function (*ξ*,*ν*),
2.21

#### (iii) Taylor coefficients for the operators S and C

All that remains is to specify formulas for the **S**_{n} and **C**_{n}. We shall show that, in fact, these coefficients are quite closely related thus enabling a significant simplification of formula (2.21). For these, it is most convenient to invent the following notation
and we note that
With this notation in hand, we have from (2.14)
and (2.15)
so that
2.22a
and
2.22b
Similarly, (2.16) and (2.17) give
2.23a
and
2.23b

To make the simplifications we advertised at the beginning of this section, we follow the lead of our previous work on these types of problems (e.g. Nicholls & Reitich 2001*a*,*b*) and consider the sum of the third and fifth terms on the right-hand side (r.h.s.) of (2.21)
We can show by a direct invocation of the product rule that
which is reminiscent of , but not evidently identical. However, from (2.23b)
where the last step comes from (2.22a), so that
Using similar manipulations for the fourth and sixth terms on the r.h.s., we can simplify (2.21) to
2.24

## 3. Governing equations: (*M*+1)-many layers

We now consider a multiple-layered material with *M* interfaces at *z*=*a*^{(m)}+*g*^{(m)}(*x*,*y*) (1≤*m*≤*M*) separating (*M*+1)-many layers which define the domains
with (upward pointing) normals *N*_{m}:=(−∂_{x}*g*^{(m)},−∂_{y}*g*^{(m)},1)^{T}. The (*M*+1) domains are all constant-density acoustic media with velocities *c*_{m} (*m*=0,…,*M*) and we again assume that plane-wave radiation (2.1) is incident upon the structure from above. In each layer, the parameter *k*_{m}=*ω*/*c*_{m} characterizes both the properties of the material and the frequency of radiation in the structure. We denote the quasiperiodic reduced scattered fields in *S*^{(m)} by *v*^{(m)}(*x*,*y*,*z*), which satisfy the Helmholtz equations
3.1
These are coupled through the Dirichlet and Neumann boundary conditions
3.2a
and
3.2b
Again, outgoing wave conditions are enforced on *v*^{(0)} and *v*^{(M)} at positive and negative infinity, respectively.

### (a) Boundary formulation

We once again make the reduction to surface quantities, and define the (lower and upper) Dirichlet traces
and their (exterior, lower and upper) Neumann analogues
In terms of these, the boundary conditions become
Again, we introduce the DNOs
so that the boundary conditions become
The first of these can be used to eliminate *V* ^{(m),u}
so that the latter equations become
or
Stated more compactly, these read
3.3
where
where

### (b) An Operator Expansions method

Following the developments of §2*b*, we now describe an OE approach to solving (3.3), once again based upon the property *g*^{(m)}(*x*,*y*)=*εf*^{(m)}(*x*,*y*) where, *a priori*, *ε* is small. Again, the data {*ζ*,*ψ*}, and operators {*G*,*H*(*m*),*J*} can be shown to depend analytically on *ε* so that
Furthermore, the fields are also analytical implying that
Upon insertion of these into (3.3), we find
where
and
where
At order , we solve
3.4
while at order
or, solving for the one unknown **V**^{l}_{n},
3.5
Once again, we point out that at order *n*, each of the are known, and that the only inversion performed is that of **A**_{0}, which is repeated at *every* perturbation order.

## 4. Numerical results

In this section, we provide detailed results of numerical simulations of scattering quantities compared with *exact* solutions. We are able to show that our numerical method is not only fast and efficient, but also accurate and applicable to configurations featuring *large* interface deformations.

### (a) Exact solutions

Of course, in the case of non-trivial interface shapes there are no known exact solutions for plane-wave incidence. To carry out a convergence study for our algorithm, we use the following principle: in building a numerical solver for the homogeneous PDE and boundary conditions
it is often, as it is here, no more difficult to construct an algorithm for the corresponding inhomogeneous problem:
Selecting an arbitrary function *w*, we can compute
and now have an *exact* solution to the problem
namely *u*=*w*. In this way, we can rigorously test our inhomogeneous solver for which the homogeneous solver is a special case. However, one does need to be careful to consider *w*, which has the same ‘behaviour’ as solutions *u* of the inhomogeneous problem and here we find *w* such that . We point out though that our exact solution does *not* correspond to plane-wave incidence (but rather to plane-wave reflection).

To be more specific, we consider the functions
4.1
with *A*^{(M)}=*B*^{(0)}=0, which satisfy the Helmholtz equations (3.1) and the outgoing wave conditions so that in the notation above. However, the boundary conditions satisfied by these functions are *not* those satisfied by an incident plane wave. With the construction of the in mind, we compute the surface data
(cf. (3.2)). We now have a family of exact solutions against which to test our numerical algorithm for *any* choice of deformations {*g*^{(1)},…,*g*^{(m)}}.

### (b) Numerical implementation

The description of our numerical scheme is rather straightforward and is one of the powerful aspects of the approach. In the (*M*+1)-layer case, we are charged with finding, e.g. the Dirichlet traces at the lower domain boundaries
from which all other quantities (Dirichlet or Neumann traces, near or far-field data or even the full field everywhere) can be computed by suitable integral formulas. Our Boundary Perturbation approach posits an expansion of the form
and we seek as an approximate solution, the truncation of this Taylor series after *N* terms
Now, *without approximation*, we can recover the **V**^{l}_{n} from the formulas (3.4) and (3.5) at orders zero and *n*>0, respectively. However, the functions which appear in these formulas generally involve Fourier series with an infinite number of non-zero coefficients. Thus, we make a spectral approximation to each of the **V**^{l}_{n}(*x*,*y*) by
Products appearing in (3.5) are computed by fast convolutions via the fast Fourier transform algorithm (Gottlieb & Orszag 1977) and our final Fourier/Taylor approximation is
4.2

### (c) Error measurement

With these approximations in hand, we can make any number of error measurements versus the exact solutions (4.1). For definiteness, we choose to measure the defect in the Dirichlet traces, one of the most difficult we can pick owing to the fact that these data are posed on the very surfaces around which we perturb. For the results described in §4*e* and 4*f* we measure
4.3a
and
4.3b

### (d) Padé approximation

Before leaving our description of the numerical procedure, we mention that there are a number of choices for summing the Taylor series which appear in (4.2). To avoid an avalanche of impenetrable notation, we focus on the generic problem of approximating the analytical function
by its truncated Taylor series
It is a classical result that if *ε*_{0} is in the disc of convergence of *A*(*ε*), say {|*ε*|<*ρ*}, *A*^{N}(*ε*_{0}) will converge to *A*(*ε*_{0}) exponentially fast as . However, it is possible for *ε*_{0} to be a point of analyticity *outside* the disc of convergence of the Taylor series and for *A*^{N} to produce meaningless results. The classical numerical analytical continuation technique of Padé approximation (Baker & Graves-Morris 1996) has been successfully brought to bear upon BPMs in the past (e.g. Bruno & Reitich 1993*b*; Nicholls & Reitich 2003) and we use this here as well. In short, Padé approximation seeks to simulate the truncated Taylor series *A*^{N} by the rational function
where *L*+*M*=*N* and
well-known formulas for the coefficients {*a*_{l},*b*_{m}} can be found in Baker & Graves-Morris (1996). This approximant has the remarkable properties that, for a wide class of functions, not only is the convergence of [*L*/*M*] to *A* at *ε*=*ε*_{0} *faster* than that of *A*^{N} for |*ε*_{0}|<*ρ*, but also that [*L*/*M*] may converge to *A* for points of analyticity *ε*_{0} for which |*ε*_{0}|>*ρ*. We refer the interested reader to §2.2 of Baker & Graves-Morris (1996) and the insightful calculations of §8.3 of Bender & Orszag (1978) for a thorough discussion of the capabilities and limitations of Padé approximants.

### (e) Numerical tests: two dimensions

To begin, we consider the two-dimensional and 2*π*-periodic case where the profiles are *independent* of the *y*-variable; we will return to the fully three-dimensional case shortly. To fix upon a reasonable class of deformations to test, we recall the three profiles introduced by Nicholls & Reitich (2001*b*) for precisely this purpose: the sinusoid
4.4a
the ‘rough’ (*C*^{4} but not *C*^{5}) profile
4.4b
and the Lipschitz boundary
4.4c
We point out that all three profiles have zero mean, approximate amplitude 2 and maximum slope of roughly 1. The Fourier series representations of *f*_{r} and *f*_{L} are listed in Nicholls & Reitich (2001*b*) and in order to minimize aliasing errors, we approximate these by their truncated *P*-term Fourier series.

We begin with two three-layer configurations.

Two smooth interfaces (figure 1

*a*): physical and numerical parameters: 4.5Rough and Lipschitz interfaces (figure 1

*b*): physical and numerical parameters: 4.6

In these, the wavelengths of propagation, given by *λ*_{j}=2*π*/*k*_{j}, are
In the first configuration (which took just over 1 s to realize with a simple Matlab script on the author's dual quad-core PowerMac), (4.5), we show that only a small number of Fourier coefficients (*N*_{x}=32) and Taylor orders (*N*=8) are required to realize machine precision (up to the conditioning of our algorithm) for small, smooth profiles, (4.4a), which displays the spectral accuracy of the scheme. In simulation (4.6), we demonstrate that the algorithm performs well if the lower and upper interfaces are replaced by the rough, (4.4b), and Lipschitz, (4.4c), profiles, respectively (both truncated after *P*=40 Fourier series terms) provided that *N*_{x} and *N* are chosen sufficiently large.

Of course, there are an infinite quantity of number/profile combinations we could consider. To give a flavour for the performance of our algorithm, we select two more in the two-dimensional setting.

Six layer (figure 2

*a*): physical and numerical parameters: 4.7201 layer (figure 2

*b*): physical and numerical parameters: 4.8

Once again, we can see that in all cases, our algorithm provides highly accurate solutions in a stable and rapid manner provided that the problem is properly resolved in space and perturbation order.

### (f) Numerical tests: three dimensions

We now consider the general case of (2*π*)×(2*π*) periodic three-dimensional interfaces. Again, we follow the lead of Nicholls & Reitich (2001*b*) and select the following interface shapes: the sinusoid
4.9a
the ‘rough’ (*C*^{2} but not *C*^{3}) profile
4.9b
and the Lipschitz boundary
4.9c
Again, all three profiles have zero mean, approximate amplitude 2 and maximum slope of roughly 1. The Fourier series representations of and are given in Nicholls & Reitich (2001*b*) and in order to minimize aliasing errors, we approximate these by their truncation after *P*=20 coefficients.

In three dimensions, despite the immediate applicability of our recursions (e.g. single convolutions are simply replaced by double convolutions), the numerical simulations become much more involved. Therefore, we focus upon the four three-layer configurations outlined below.

Two smooth interfaces (figure 3

*a*): physical and numerical parameters: 4.10Smooth and rough interfaces (figure 3

*b*): physical and numerical parameters: 4.11Smooth and Lipschitz interfaces (figure 4

*a*): physical and numerical parameters: 4.12Rough and Lipschitz interfaces (figure 4

*b*): physical and numerical parameters: 4.13

We once again notice the rapid and stable convergence of our numerical scheme in all four instances. The behaviour is independent of interface shape provided that a sufficient number of discretization points are used to properly resolve all features of importance in the simulation.

## 5. Conclusion

In this paper, we have presented a novel surface formulation of acoustic scattering by layered media. We have also generalized the OE method of Milder and Malcolm & Nicholls to the case of (*M*+1)-many layers. In this formulation, we realize the *minimum* number of problem unknowns (*M* surface functions, also used by BIM/BEM) which halves the number currently mandated by the approach of Malcolm & Nicholls (2011*b*). Our new approach has rather *explicit* dependence upon the interface shapes which makes their application towards an inverse problem algorithm very appealing, which we leave to future work.

- Received September 13, 2011.
- Accepted October 17, 2011.

- This journal is © 2011 The Royal Society