## Abstract

In this paper, we develop an asymptotic scheme to approximate the trapped mode solutions to the time harmonic wave equation in a three-dimensional waveguide with a smooth but otherwise arbitrarily shaped cross section and a single, slowly varying ‘bulge’, symmetric in the longitudinal direction. Extending previous research carried out in the two-dimensional case, we first use a WKBJ-type ansatz to identify the possible *quasi-mode* solutions that propagate only in the thicker region, and hence find a finite *cut-on* region of oscillatory behaviour and asymptotic decay elsewhere. The WKBJ expansions are used to identify a *turning point* between the cut-on and *cut-off* regions. We note that the expansions are non-uniform in an interior layer centred on this point, and we use the method of matched asymptotic expansions to connect the cut-on and cut-off regions within this layer. The behaviour of the expansions within the interior layer then motivates the construction of a uniformly valid asymptotic expansion. Finally, we use this expansion and the symmetry of the waveguide around the longitudinal centre, *x*=0, to extract trapped mode wavenumbers, which are compared with those found using a numerical scheme and seen to be extremely accurate, even to relatively large values of the small parameter.

## 1. Introduction

Acoustic propagation in waveguides and similar domains with infinite boundaries is a class of problems that is of broad interest, and there is a long history of research into various methods of calculating their acoustic modes [1,2]. By *trapped modes*, we are referring to acoustic modes wherein we find a localized region of finite oscillatory energy, and which decay elsewhere. The existence of trapped modes in various types of topographically slowly varying waveguides has been demonstrated previously, in the cases of elastic plates and rods [3,4], in two-dimensional acoustic waveguides [5,6] and in weakly curved quantum waveguides [7]. Additionally, there has been previous research into the calculation of *quasi-modes*—a term borrowed from Gridin & Craster [8,9] in order to denote perturbed modes of a slowly varying waveguide—in both two- and three-dimensional acoustic waveguides with a slow curvature [8,10], and in two-dimensional elastic plates [9], among others.

In particular, in Biggs [5], an asymptotic scheme is developed to approximate the trapped mode solutions of the time-harmonic wave equation in a two-dimensional waveguide with a slowly varying symmetric bulge around the longitudinal centre *x*=0. This slowness is incorporated into the scheme as a large longitudinal length scale, *L*_{x}∼1/*ϵ*, where *ϵ*≪1 is taken as the small parameter used to construct an asymptotic approximation to the solution. By expanding both the amplitude and phase in powers of *ϵ*, the cut-on and cut-off regions are determined by the relative magnitudes of the wavenumber, *k*, and the eigenvalues of the one-dimensional duct cross section at each *x*. A turning point between the two regions is identified, and the solution behaviour around this motivates the construction of a uniformly valid approximation, which is used to calculate the trapped mode eigensolutions.

In this paper, we extend this work to create an asymptotic approximation to the trapped mode solutions within a three-dimensional waveguide with a smooth and simply connected, but otherwise arbitrarily shaped cross section with a slow, symmetric variation centred on the plane *x*=0. In §2, the slowly varying nature of the waveguide boundary is introduced in terms of the small parameter, 0<*ϵ*≪1, and the geometry is non-dimensionalized and mapped onto a stretched coordinate system, such that this parameter is incorporated into the differential operator.

In §3, we first determine the possible oscillatory or evanescent quasi-mode solutions [5,8,10] that may be present within the waveguide by expanding the solution using a WKBJ-type ansatz of the form , wherein the amplitude and phase terms are expanded as *A* = *A*_{0} + *ϵA*_{1} + ⋯ and *P*=*ϵ*^{−1}*P*_{−1}+*ϵP*_{1}+⋯ , respectively. A hierarchy of equations in the *A*_{i} and *P*_{i} terms is obtained, which we solve iteratively to obtain a leading-order quasi-mode approximation to *ϕ*. It will be seen that the behaviour of the solution depends upon the relative sizes of the wavenumber, *k*, and the transverse eigenvalues of the waveguide cross section, which we denote *λ*(*x*).

In §4, the existence of a turning point is discussed, around which the solution is either oscillatory or evanescent. Motivated by the quasi-mode solutions found in §3, we find this to be a point *x*_{*}, for which *λ*(*x*_{*})=*k*. Taking established results in domain perturbation theory [11,12], we are able to assert that a positive turning point will be found in a region for which the cross-sectional area is decreasing, such that we expect a solution to be oscillatory only in the bulge, up until the turning point, and decaying towards infinity.

An interior layer is then introduced around the turning point, wherein an interior solution is found and asymptotically matched with the oscillatory and evanescent solutions on either side of the region. We find that the interior solution takes the form of an Airy function in the *x*-direction, which motivates the uniformly valid expansion introduced in §6. Using this expansion, we are able to extract trapped mode wavenumbers by considering the symmetry of the waveguide about the *x*=0 plane. As such, we find that, to leading order, the wavenumbers of the trapped modes that are either symmetric or antisymmetric in the *x*-direction are those for which *ϕ*_{x}(0,*y*,*z*)=0 or *ϕ*(0,*y*,*z*)=0, respectively. We find there to be a sequence of trapped mode wavenumbers in each case of the form for the Laplacian eigenvalues *λ*_{n}(*x*) of the waveguide cross section at *x*.

Lastly, in §8, both the asymptotic scheme and a spectral collocation method are applied to calculate the trapped mode wavenumbers of two example problems, to demonstrate the accuracy of the former. The first test case considered is that of an azimuthally constant ‘cylindrical’ waveguide with a slow variation. Owing to this azimuthal invariance, the problem in this particular geometry may be reduced to a sequence of two-dimensional problems, so that wavenumbers may be quickly extracted via the spectral method, and we see an excellent agreement between these and those obtained using the asymptotic scheme. Further benefits of the asymptotic method are shown in the second example problem, wherein no such reduction is possible and a three-dimensional spectral collocation or similar numerical method proves to be prohibitively expensive. In such a case, we see that the asymptotic scheme can still be used to efficiently approximate the trapped mode wavenumbers to a high degree of accuracy.

## 2. Formulation

We are interested in time harmonic solutions to the wave equation, and hence omit the common factor throughout. The problem is then equivalent to that of finding trapped mode solutions to the Helmholtz equation within a three-dimensional waveguide of arbitrary cross section and smooth boundary. Using standard Cartesian coordinates with *x* oriented along the waveguide, we take the duct to be infinite in the *x*-direction and symmetric about the plane *x*=0. We will take the cross section at *x* to be a simply connected domain , with smooth boundary . We seek eigensolutions of
2.1We will consider individually both the sound soft and sound hard waveguide, i.e. that with either a Dirichlet or Neumann condition on the waveguide boundary, respectively, given by
2.2aand
2.2bwith denoting the outward normal vector to at (*x*,*y*,*z*). If we define a non-dimensionalization constant *L*_{yz}, which characterizes the length scale of the waveguide in the *y*- and *z*-directions, by
and similarly take *L*_{x} to characterize the length scale of the perturbation in the *x*-direction, we therefore take the waveguide to be slowly varying in the sense that *L*_{yz}/*L*_{x}=*ϵ*, where 0<*ϵ*≪1. The exact nature of *L*_{x} will depend upon the geometry under consideration; in those considered in §8, for example, for which the perturbation does not have compact support and instead decays exponentially with |*x*|, *L*_{x} is equal to the e-folding length scale of the decay. We thus non-dimensionalize by introducing , where *ξ*=*x*/*L*_{x}, *η*=*y*/*L*_{yz} and *ζ*=*z*/*L*_{yz} on the stretched domain . Transforming (2.1) into these stretched coordinates, we have the non-dimensional system
2.3where is the non-dimensionalized wavenumber. The boundary conditions (2.2a) and (2.2b) are then equivalent to
2.4aand
2.4brespectively, where ∇_{ϵ}=(*ϵ*^{2}∂_{ξ},∂_{η},∂_{ζ}) is a skewed del operator and
is the transformed normal vector.

## 3. WKBJ ansatz

We first seek quasi-mode solutions valid within the waveguide by expanding *ϕ* using a WKBJ-type ansatz [13], of the form
3.1where
We note that there is no *P*_{0} term, because we observe that the *O*(1) phase contribution can be subsumed into the leading-order amplitude term. In fact, as will be seen, it is possible to obtain a good leading-order approximation with just the leading-order phase component. We include the phase expansion in order to preserve generality and to permit the possibility of calculating higher order refinements. We also require the phase terms *P*_{j} to be either strictly real or strictly imaginary, because we are seeking modes that are either oscillatory or evanescent in different regions within the waveguide. By substituting the above ansatz into (2.3) and equating powers of *ϵ*, we extract a hierarchy of equations in *A*_{j} and *P*_{j}, from which we can determine the leading-order terms of an asymptotic expansion for *ϕ*. At leading order, *O*(*ϵ*^{−2}), we obtain
3.2which, because we have specified that either or and we require a non-zero leading-order amplitude for a non-trivial solution, implies that *P*_{−1}(*ξ*,*η*,*ζ*)≡*f*(*ξ*), where *f* is a function to be determined. We find that the *O*(*ϵ*^{−1}) terms are trivially satisfied and hence we turn to the *O*(1) terms, at which point we consider separately the case of the sound soft boundary, as in (2.4a), and the sound hard boundary, as in (2.4b).

### (a) Sound soft boundary

Note that, to satisfy (2.4a), it is sufficient to choose
such that *A*_{j}=0 on the boundary ∂*D*(*ξ*) for all . We denote ∇_{ηζ}=(∂_{η},∂_{ζ}) to represent the two-dimensional del operator in *η* and *ζ* only. Then, considering the *O*(1) terms in the asymptotic expansion, we have that
3.3on the domain *D*(*ξ*), with boundary condition *A*_{0}=0 on ∂*D*(*ξ*), where we write *λ*^{2}(*ξ*)= *k*^{2}+*f*^{′2}(*ξ*) for notational convenience. We consider a Dirichlet eigenvalue problem over an arbitrarily shaped domain, *Ω*(*ξ*),
3.4with eigensolutions *λ*_{n}(*Ω*), *E*_{n}(*Ω*) for all and , where *E*_{n} are normalized such that
We will henceforth omit the subscript notation on the del operator for brevity, such that we refer to ∇_{ηζ} as ∇, unless otherwise specified. We may then write the solutions of (3.3) as *A*_{0;n}=*α*_{n}(*ξ*)*E*_{n} with associated eigenvalues *λ*_{n}(*ξ*). Now, briefly disregarding the subscript *n* notation for convenience, at *O*(*ϵ*) we have
3.5which we multiply by *A*_{0} and integrate over *D*(*ξ*) to obtain the solvability condition
We apply equation (3.3) and use Green's identities to yield
where denotes the outward unit normal to ∂*D*(*ξ*). Then because the *A*_{j} terms vanish on ∂*D*, we retain only the condition
We turn now to the general form of the Reynolds transport theorem, which states that for a region *Ω*(*t*), with boundary ∂*Ω*(*t*),
3.6where ** v** denotes the rate of change of the boundary with

*t*. We apply this to the two-dimensional domain

*D*(

*ξ*) and see that, because the Dirichlet condition requires

*A*

_{0}to vanish on the boundary, the derivative may simply be taken outside of the integral and the solvability condition reduces to so that we have 3.7for some constant . Therefore, reapplying the subscript notation, we find the leading-order quasi-mode solutions to be 3.8where 3.9for constants .

### (b) Sound hard boundary

In the case of the sound hard waveguide, the boundary condition (2.4b) applied to the WKBJ ansatz tells us that, at leading order, we require
3.10on ∂*D*(*ξ*). The *O*(1) terms then provide the two-dimensional equation (3.3) as before, albeit now with the boundary condition upon *A*_{0} given by
3.11where ∇≡∇_{ηζ}, as specified above, and denotes the (*η*,*ζ*)-direction components of the outward normal, ** n**, to the boundary ∂

*D*(

*ξ*).

Similar to the sound soft waveguide, we consider a general eigenvalue problem over a domain *Ω*(*ξ*),
3.12where ** n** is the outward normal to the domain, and write its eigenvalues and normalized eigenfunctions as , , respectively, for . Now given that the leading-order amplitude in the sound hard duct satisfies (3.3) and (3.11), it admits the eigensolutions , with eigenvalues , for . We again briefly omit the subscript notation, and hence the solvability condition from the

*O*(

*ϵ*) terms can be written as 3.13which, through an application of Green's identities and (3.3), is equivalent to 3.14where denotes the outward unit normal to ∂

*D*(

*ξ*). The

*O*(

*ϵ*) terms of the Neumann boundary condition require that which can be written as 3.15on ∂

*D*(

*ξ*), where

*n*_{jk}is as before, and similarly

*n*

_{i}denotes the component of the normal vector acting in the

*ξ*-direction. Thus, applying both this and condition (3.11) to the solvability condition yields 3.16We denote by

**the rate of change of the waveguide boundary with respect to**

*v**ξ*and see from appendix A that (3.16) is equivalent to We may thus apply the Reynolds transport theorem in the form presented in (3.6) to obtain, as in the Dirichlet case, for constant . Hence, reintroducing the subscript notation, we have the sound hard quasi-mode solutions 3.17where 3.18for constants .

These solutions thus imply that the presence of oscillations in the *ξ*-direction depend upon the relative magnitudes of the wavenumber, *k*, and the eigenvalues, *λ*_{n}(*ξ*), of the Laplacian over the domain *D*(*ξ*). In the case presented here, in which we have a waveguide symmetric about the plane *ξ*=0, we require, for a trapped mode solution, a region of oscillatory energy up to a particular point, 0≤|*ξ*|≤|*ξ*_{*}|, and decay for |*ξ*|≥|*ξ*_{*}|. In fact, we may restrict our attention to just the *ξ*≥0 half of the waveguide, since it is symmetric, and we see from (3.9) and (3.18) that we hence require a turning point *ξ*_{*}>0 satisfying *λ*_{n}(*ξ*_{*})=*k* and *λ*_{n}′(*ξ*_{*})>0 for such a solution to exist.

## 4. Turning point

We briefly discuss a result in domain perturbation theory for a necessary condition for a suitable turning point to exist in the Dirichlet waveguide. Given an arbitrary, piecewise smooth domain and a domain , with a small (i.e. *O*(*δ*) for 0<*δ*≪1), but not necessarily localized, perturbation from *Ω*(*x*), it follows from the work of Kato [11] and Rellich [12] that the perturbed Dirichlet eigenvalues can be expressed as a uniform asymptotic expansion in *δ*. In particular, if we denote by *λ*_{j} a particular Dirichlet eigenvalue on *Ω*(*x*), and by the associated perturbed eigenvalue on *Ω*^{δ}(*x*), then we can expand
We then find that, for the leading-order variation, , wherein *f* is the normal distance between the boundaries ∂*Ω* and ∂*Ω*^{δ}. We may say, then, that a decrease in the size of the domain coincides with an increase in the magnitude of the eigenvalues of the Dirichlet Laplacian on it. In terms of our sound soft waveguide, a necessary condition for the existence of a turning point *ξ*_{*}>0, with the properties described above, is for the size of the waveguide cross section to be decreasing at *ξ*_{*}.

In the example waveguide geometry that we will return to in §8*a*, for instance, we consider a rotationally symmetric waveguide with a circular cross section of radius *h*(*ξ*). It is straightforward to solve the eigenproblem over the cross section analytically and show that the eigenvalues are given by *λ*_{mn}(*ξ*)=*β*_{mn}/*h*(*ξ*), where *β*_{mn} denotes the *n*th root of the Bessel *J*_{m} function. Hence, we see directly in this case that the magnitude of the eigenvalues is inversely proportional to the radius of the cross section, as expected.

Unfortunately, no such result holds in general for the Neumann boundary and thus we have no simple sufficient condition on the waveguide geometry for a trapped mode solution to exist within the sound hard duct. However, we may be able to make an assertion about the behaviour of the eigenvalues of the Neumann Laplacian over its cross section, particularly for certain ‘simple’ geometries. Returning to the example of the waveguide of circular cross section, as mentioned earlier, we find that the Neumann eigenvalues over the circle of radius *h*(*ξ*) are given similarly by *λ*_{mn}(*ξ*)=*β*′_{mn}/*h*(*ξ*), where *β*′_{mn} denotes the *n*th root of the derivative of the Bessel *J*_{m} function.

## 5. Asymptotic matching

Given the quasi-mode solutions (3.8), with (3.9), and the existence of a turning point, *ξ*_{*}>0, we now discuss the behaviour of the solution near this turning point, and show that the asymptotic expansions exhibit non-uniformity within a small layer around it. To do so, we consider the related reflection–transmission problem formed by the half-waveguide *D*(*ξ*) in and the constant duct *D*(*ξ*)≡*D*(0) for .

### (a) Region of non-uniformity

We consider for clarity just the sound soft waveguide, but note that the analysis shown also holds for the sound hard duct. Given that the bulge maximum is at *ξ*=0 and the waveguide tapers with increasing *ξ*, we can write the general quasi-mode solutions in the cut-on and cut-off regions as
5.1aand
5.1bfor *k*≥*λ*(*ξ*) and *k*≤*λ*(*ξ*), respectively, where the oscillatory solution is composed of an incident and reflected component. We expect there to be a layer of non-uniformity between these two regions within which these expansions break down around a turning point *ξ*_{*}, where *λ*(*ξ*_{*})=*k*, and consider the behaviour of the quasi-mode solutions as . Taking the Taylor series of *λ*^{2}(*ξ*) near *ξ*=*ξ*_{*}, we have that
Writing *Γ*=2*λ*′(*ξ*_{*})*λ*(*ξ*_{*}), so that |*k*^{2}−*λ*^{2}(*ξ*)|∼|*Γ*(*ξ*−*ξ*_{*})|, we can consider (5.1b), for example, to see that
from which it is evident that the quasi-mode expansions break down in the region where (*ξ*−*ξ*_{*})=*O*(*ϵ*^{2/3}).

### (b) Interior layer

Having found the region within which the WKBJ expansion breaks down, we place an interior layer of width *O*(*ϵ*^{2/3}) in this region by considering a solution in terms of the stretched variable *χ*=*ϵ*^{−2/3}(*ξ*−*ξ*_{*}), which is *O*(1) in the region of interest. We suppose that *ϕ*=*X*(*χ*)*E*, where *E*_{ηη}+*E*_{ζζ}+*λ*^{2}*E*=0, and hence reformulate the original equation (2.3) such that it becomes
Then by taking the Taylor expansions of *E* and its derivatives near *ξ*_{*}, and noting that the *O*(1) terms vanish, we are left with
Hence, at leading order, we require
5.2whose solution is given by
5.3where Ai and Bi are the Airy functions of the first and second kind, and *F* and *G* are constants. We note that we retain only the Ai component of (5.3), because the Airy Bi function is not appropriate to the trapping problem, owing to its exponential growth. Given (5.2), we see that *Γ*=2*λ*′(*ξ*_{*})*λ*(*ξ*_{*})>0, where *λ*(*ξ*) is an eigenvalue for the domain *D*(*ξ*), which corresponds to a cut-on region *ξ*≤*ξ*_{*}, in the thicker part of the waveguide.

### (c) Matching coefficients

Given the tapering geometry of the waveguide, we seek solutions that are oscillatory for *ξ*≤*ξ*_{*} and decay in *ξ*≥*ξ*_{*}. Following this assumption, we rewrite the expansion (5.1a) and (5.1b) as
5.4aand
5.4bfor *ξ*≤*ξ*_{*} and *ξ*≥*ξ*_{*}, respectively. Then in the oscillatory region, we have that, as ,
whereas the Airy function in (5.3) has the asymptotic expansion [14] as
and we hence require the constants *F*, *R* and *I* to satisfy
5.5in order to match these outer and inner expansions. Similarly, if we look at expansion (5.4b) as , we see that
which matches with the interior expansion for
provided that
5.6It is then straightforward to substitute (5.5)–(5.6) into (5.4a)–(5.4b) and subsequently construct two one-term composite expansions, valid across either side of the interior layer.

## 6. Uniform expansion

Having observed the appearance of an Airy function of the first kind in the interior layer connecting the oscillatory and evanescent quasi-mode solutions, and indeed motivated by its use in constructing a uniformly valid expansion for the two-dimensional waveguide in Biggs [5], we now seek a uniformly valid expression for *ϕ* by expanding in terms of Airy functions. We introduce the new ansatz
6.1where
and we restrict *g*(*ξ*) to be real. The Ai′ function is included to form an orthogonal basis, but we note that only the Ai function appears to leading order in the interior layer expansion (5.3), and so omit an *O*(1) Ai′ coefficient. We also note here that the presence of an interior layer of width *O*(*ϵ*^{2/3}) suggests that our ansatz should be expressed in powers of *ϵ*^{2/3}. Inserting this ansatz into the governing equation (2.3) yields a new hierarchy of equations in powers of *ϵ*, each in coefficients of the linearly independent functions Ai and Ai′. As with the quasi-mode expansions, we will consider first the sound soft waveguide, with a Dirichlet condition on its surface, followed by the sound hard duct with a Neumann condition.

### (a) Sound soft boundary

In the sound soft case, the boundary condition (2.4a) under this ansatz resolves simply to the requirement that
6.2So looking first at the *O*(1) coefficients of Ai, we have
6.3where we now define . We denote the eigensolutions of (6.3) over *D*(*ξ*), with boundary condition *B*_{0}(*ξ*,*η*,*ζ*)=0 for (*η*,*ζ*)∈∂*D*(*ξ*), by *B*_{0;n}=*β*_{n}(*ξ*)*E*_{n}, with associated eigenvalues *λ*_{n}(*ξ*), where *E*_{n} and *λ*_{n} are the eigensolutions of (3.4) for .

For the sake of notational convenience, we take it as understood for the time being that we are dealing with an arbitrary *n*th eigensolution and omit the subscript *n* notation. Following the convention in §3, we also remove the subscript *ηζ* from ∇_{ηζ}. Because the eigenvalue *λ*(*ξ*) is now known from (6.3), we have an expression for *g*_{0}(*ξ*), taking care to note the separation between cases where . As in the notation above, if we assume there to be a turning point *ξ*_{*}>0 such that *λ*(*ξ*_{*})=*k* and *λ*′(*ξ*_{*})>0, then we may write
6.4We then find an expression identical to (6.3) in *C*_{1} from the *O*(*ϵ*^{2/3}) coefficients of Ai′, such that *C*_{1}=*γ*(*ξ*)*E* with eigenvalue *λ*. Now the *O*(*ϵ*^{2/3}) coefficients of Ai are
which we multiply through by *B*_{0} and integrate over *D*(*ξ*) to yield the solvability condition
6.5The Dirichlet boundary condition causes the first integral to disappear upon an application of Green's theorem, which along with the orthonormality of *E* leaves just
6.6for some constant . Turning to the *O*(*ϵ*^{4/3}) coefficients of Ai′, we have that
6.7which we once again integrate with *B*_{0} and, given that *C*_{2} vanishes on the boundary, find
Because *B*_{0} vanishes on the boundary ∂*D*(*ξ*), we may take the derivative outside of the integral, such that
6.8for constant . Hence, to leading order, we have the uniformly valid asymptotic expansions
6.9for constants , where in a slight amendment to the notation in (6.4), to highlight the dependence of *ξ*_{*} upon *k*,
6.10

### (b) Sound hard boundary

In the case of the sound hard boundary, we return to (6.1) under boundary condition (2.4b). Expanding this boundary condition in powers of *ϵ* and collecting coefficients of Ai and Ai′, we once again extract a hierarchy of equations that must be satisfied. Firstly, taking the *O*(1) coefficients of the Ai terms, we have the condition
on ∂*D*(*ξ*). Thus, at leading order, we have the PDE (6.3) with a Neumann boundary condition on ∂*D*(*ξ*). Analogous to the sound soft case, we will denote the eigensolutions of such as , with eigenvalues , as in (3.12). Once again, we also extract a similar PDE for *C*_{1} from the *O*(*ϵ*^{2/3}) coefficients of Ai′ with a Neumann boundary condition on ∂*D*, and thus we find the eigensolutions of *C*_{1} to be , with eigenvalues . We shall briefly remove the subscript notation for clarity. Then, the Neumann solvability condition is found to be identical to (6.6), and so we turn to the *O*(*ϵ*^{4/3}) coefficients of Ai′, wherein we find the next solvability condition
Returning to the hierarchy of equations extracted from the boundary condition, we find that the *O*(*ϵ*^{4/3}) coefficients of Ai′ yield the condition
6.11on ∂*D*(*ξ*), such that the solvability condition becomes
Hence using the property given in Appendix A, the Reynolds transport theorem and the orthonormality of , we have that
as in (6.8) in the sound soft case. Thus, reintroducing the subscript notation, we have the uniformly valid solution for the sound hard duct
6.12with
6.13

## 7. Trapped modes

Given the leading-order uniform expansions (6.9) and (6.12), valid for , we notice that, by considering a reflection of this semi-infinite duct through the (*η*,*ζ*)-plane, we may extract leading-order trapped mode solutions by considering those solutions that are either antisymmetric or symmetric about *ξ*=0. We hence choose appropriate wavenumbers, *k*, for which either *ϕ*(0,*η*,*ζ*)=0 or *ϕ*_{ξ}(0,*η*,*ζ*)=0, respectively. Taking just the Dirichlet waveguide, for example, this is equivalent to requiring that the trapped mode wavenumbers *k* are solutions of either
7.1aor
7.1bwhere we take *z*_{p} and *z*′_{p} to be the zeros of Ai and Ai′, respectively, ordered such that 0>*z*_{1}>*z*_{2}…, and similarly for *z*′_{p}.

## 8. Numerical verification

We note that the full numerical solution of (2.1) is, in an arbitrarily shaped waveguide, far from straightforward and, in general, computationally expensive. In the case of a waveguide with circular cross section, however, we may use the azimuthal periodicity of the solution to reduce the problem to a two-dimensional one, to which we may apply a spectral collocation method [15,16] to demonstrate the accuracy of the asymptotic scheme. We find spectral methods to be particularly appropriate in the numerical solution of waveguide problems [3,10]. This is due in part to the simplicity with which one may incorporate an infinite boundary and exponential decay into the solution by expanding in Laguerre polynomials in the longitudinal direction.

### (a) Waveguide with circular cross section

We consider first the case of the rotationally symmetric, cylindrical waveguide with a slow, uniform bulge around *x*=0. It is natural to describe this geometry in cylindrical coordinates and we are thus seeking trapped mode solutions of
8.1where we will consider the example waveguide profile *h*(*ξ*,*θ*)≡*h*(*ξ*)=1+(*h*_{1}−1) sech *ξ*, such that we are considering a duct with unit radius as *ξ* approaches and which bulges to a maximum radius of *h*_{1} at *ξ*=0. We will present the numerics for the sound soft duct, so that we also require *ϕ*(*ξ*,*h*(*ξ*),*θ*)=0, although note that it is a straightforward adjustment to instead consider the sound hard case, wherein *ϕ*_{r}(*ξ*,*h*(*ξ*),*θ*)=0. In this case, the waveguide cross section at any *ξ* is merely a circle of radius *h*(*ξ*) and, as such, we are able to derive the transverse eigenvalues *λ*(*ξ*) and uniform solution analytically. Following the methods of §6, we obtain the uniformly valid, leading-order eigensolutions
8.2with associated eigenvalues *λ*_{mn}(*ξ*)=*β*_{mn}/*h*(*ξ*), where *β*_{mn} denotes the *n*th root of the Bessel *J*_{m} function. We note here that, in fact, the Dirichlet Laplacian on the circular cross section has, for all except the zeroth azimuthal mode, double eigenvalues, with either a cosine or sine eigenfunction. In this unforced problem we choose, without loss of generality, the cosine solution to preserve the zeroth azimuthal mode. Given the waveguide topography under consideration, it is clear that, for all positive *ξ*, *h*′(*ξ*)≤0, and hence *λ*′_{mn}(*ξ*)>0. Thus, as discussed in §4, we may reasonably anticipate the existence of one or more turning points per wavenumber, and hence associated trapped mode solutions.

We notice that the separability and behaviour of (8.2) in the *θ*-direction permits the rewriting of (8.1) as
8.3such that the *θ*-derivative has been removed and hence the three-dimensional problem has been reduced to *m* two-dimensional problems. In general, this modification will not be possible, and to apply a full three-dimensional spectral collocation method with *M*, *N* and *P* interpolation points in the spatial directions will require the construction and calculation of the eigenvalues of an (*M*×*N*×*P*)×(*M*×*N*×*P*) differentiation matrix—not a computationally cheap exercise.

In order to apply the spectral method, we map the domain onto a semi-infinite rectangular domain via the transformation *ρ*=2*r*/*h*(*ξ*)−1, such that it follows from (8.3) that the function *ψ*(*ξ*,*ρ*)=*ϕ*(*ξ*,*r*) satisfies
for *ρ*∈[−1,1), , where
and with boundary conditions as and *ψ*(*ξ*,1)=0. We may simplify the problem further by noting that we are seeking either symmetric or antisymmetric modes, and may consider only the half plane , with the respective boundary conditions *ψ*_{ξ}(0,*ρ*)=0 or *ψ*(0,*ρ*)=0 and as .

We then apply a Chebyshev interpolant in the *ρ*-direction and a Laguerre interpolant in the *ξ*-direction to discretize (8.3) and approximate the solution.

Taking *ϵ*=0.1 as our small parameter, figure 1 depicts the comparisons between the asymptotic scheme and spectral approximation for the first azimuthal and radial modes, i.e. for *m*=*n*=1, using a spectral grid of 60×15 points in the *ξ*-, *ρ*-directions, respectively. This is sufficient to give the ‘tenth’ longitudinal wavenumber, , to eight decimal places of accuracy. It can be seen that the leading-order asymptotic approximation is in excellent agreement with the numerical scheme. Indeed, turning to table 1, we see that, if we fix *h*_{1}=1.5 as a sample bulge radius, the first wavenumber carries the greatest error at ≈ 0.136 per cent, and in fact the average error over the first eight calculated wavenumbers is just ≈ 0.032 per cent.

It is also interesting to note that, while we have exploited the ‘small’ parameter *ϵ* in order to construct an asymptotic scheme, we now find that the accuracy of the approximation remains fairly high with much larger *ϵ*. Table 2 shows the relative error for the first symmetric and antisymmetric trapped mode wavenumbers for the *m*=*n*=1 modes with varying *ϵ*, wherein we see that, even with a value of *ϵ*=1, the relative error for the first symmetric trapped mode is only ≈ 1.156 per cent.

### (b) Elliptical waveguide

The method presented here can be applied to any waveguide geometry featuring a smooth cross section and slow, axially symmetric perturbation, and is of the most interest in any case where it is extremely difficult or computationally expensive to numerically solve the full three-dimensional problem. We demonstrate this by considering a waveguide that is elliptical at large |*ξ*| and has a localized, one-sided bulge around *ξ*=0. Once more in cylindrical coordinates, we consider the problem described by (8.1), with the waveguide profile given instead by
8.4and a bulge height *h*_{1}=1.5, as shown in figure 2. We map this geometry onto a circular cross section cylinder of unit radius in order to apply a spectral collocation method, using the transformation *ρ*=*r*/*h*(*ξ*,*θ*). Thus (8.1) becomes
where, similar to the previous case,
and such that *ρ*∈[0,1), *θ*∈[0,2*π*] and . It is of note here that, following the methods of Trefethen [15], we map *ρ* onto the interval [0,1), as opposed to the ‘full’ Chebyshev interval [−1,1), and later discard the negative Chebyshev interpolation points. This ensures a much more even spacing of Chebyshev points in the *ρ*-direction and thus far greater accuracy for a marginal computational cost. We note, however, that, even with such a modification, a full spectral collocation method with *M*, *N* and *P* interpolation points in the *ξ*-, *r*- and *θ*-directions, respectively, involves calculating the eigenvalues of an (*M*×*N*×*P*)×(*M*×*N*×*P*) matrix. For the purposes of this verification, we apply this method with an interpolation grid of (18×18×20) points, such that the problem is discretized into a 6480×6480 matrix.

Comparing the results of this numerical scheme with our leading-order asymptotic results, we see from table 3 that, for the first eight wavenumbers, they are once again in excellent agreement. Furthermore, once the wavenumbers have been found, the full leading-order solution is immediately obtained from substitution into (6.9) and (6.10). Figure 3 shows the cross section of the first four leading-order antisymmetric trapped mode solutions (*k*_{np} for *n*=1,…,4) inside this particular waveguide geometry, taken slightly away from the bulge centre, at *ξ*=0.2.

## 9. Conclusion

In this paper, we have presented an extension of the asymptotic scheme derived in Biggs [5] for calculating the trapped mode solutions inside a two-dimensional sound soft or sound hard waveguide, and shown that it is applicable and numerically accurate in the case of a three-dimensional duct. In doing so, we have replaced a computationally expensive and cumbersome numerical solution with a method involving only the solution of a nonlinear algebraic expression.

We also note that, despite the construction of the scheme hinging on the slowly varying nature of the waveguide geometry (i.e. a variation over a length scale 1/*ϵ* for small *ϵ*), the leading-order approximations obtained are still accurate in some cases to a relative error of ≈1 per cent when the ‘slow’ parameter is as large as *O*(1).

Furthermore, we have introduced no restrictions on the overall size or shape of the perturbation, other than the smoothness of the cross section and symmetry around *x*=0, such that the method presented here may be applied to a variety of waveguide geometries in which we have an axially symmetric perturbation. We propose to extend this work to encompass other waveguide geometries with relaxed conditions, wherein the cross section need not be smooth, and to consider the case of trapped modes in an elastic plate.

## Acknowledgements

The authors thank Prof. Michael Levitin and Prof. Michael Baines for their helpful discussions and interest in this work. We are also grateful to the EPSRC, which financially supports this research.

## Appendix A

Given a smooth surface ** S**=(

*x*,

*y*,

*z*), parametrized by it is clear that the rate of change of the cross section with

*x*is given by

**=(**

*v**η*

_{x},

*ζ*

_{x}). We also have that the surface normal is given by the cross product and hence

- Received June 27, 2012.
- Accepted September 7, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.