## Abstract

The coupled-mode equations of Shevchenko are extended to three-dimensional irregular acoustic waveguides. These equations are solved for a model of a cylindrically symmetric seamount and a time-harmonic point source in a horizontally stratified ocean. An algebraic formula for the scattered field beyond the seamount is obtained from this solution. The formula is characterized by a set of infinite-dimensional matrices that are independent of the source. A stable numerical procedure is developed to compute finite-dimensional approximations to these matrices for use in truncated versions of the formula.

## 1. Introduction

This paper uses coupled-mode theory to solve a scattering problem for a simple model of a seamount. The seamount is cylindrically symmetric, and its surface has a finite slope at every point. Sound from a time-harmonic point source that is located beyond the seamount is incident on it, and the problem is to compute the scattered field beyond the seamount.

A recent paper [1] that deals with this problem is based on a ‘stepwise’ formulation of coupled-mode theory. This formulation replaces the ‘ideal’ seamount of the actual problem with a ‘step’ structure made up of nested cylindrical shells, and the scattering problem is solved for this step structure instead of the ideal seamount. The step structure casts an acoustic shadow that is almost identical to the shadow cast by the ideal seamount. However, the fields outside the shadows differ because the step structure does not emit lateral waves from the same locations or in the same directions as the ideal seamount.

This paper introduces a formulation of coupled-mode theory that deals directly with the ideal seamount. It is based on partial differential equations that govern propagation of coupled modes in three-dimensional waveguides. A somewhat lengthy analysis is required to convert these equations into systems of ordinary differential equations that are suitable for numerical work. However, these systems are no harder to program than the algebraic systems on which the stepwise formulation is based. They can be integrated by standard numerical methods that are stable and accurate if the integration steps are small enough to resolve the coupling processes that are associated with emission of lateral waves and radiation into the bottom. These aspects of numerical implementation, as well as a computational validation of the theory, are discussed in the electronic supplementary material.

If the source is located beyond the seamount, then the scattered field beyond the seamount can be computed from a simple algebraic formula. This formula is characterized by a set of matrices that are independent of the source. These scattering matrices are infinite-dimensional and there are infinitely many of them. However, numerical implementations of the coupled mode solution can be adapted to compute finite-dimensional approximations to these matrices for truncated versions of the formula that are suitable for practical work.

This paper is organized as follows. Partial differential equations that govern propagation of coupled modes in three-dimensional waveguides are derived in §2. The derivation is done in rectilinear coordinates and does not depend on any symmetries of the waveguides. In §3, these equations are reformulated in cylindrical coordinates for waveguides with cylindrical symmetry, and systems of ordinary differential equations for the scattered field of the seamount are obtained from them by Fourier analysis. In §4, these systems are converted to equivalent systems that can be used as a basis for numerical work, and a Riccati transformation is proposed as a practical solution method. In §5, the solutions of these systems are used to construct the algebraic formula for the scattered field and the scattering matrices. It is demonstrated in the appendix that the scattering matrices are symmetric.^{1}

## 2. First-order coupled-mode equations for three-dimensional waveguides

Consider a three-dimensional plane-parallel waveguide. Let the coordinates *x*, *y* and *z* indicate range, cross-range and height, respectively. The surface of the waveguide coincides with the horizontal plane in which *z*=0, and the bottom of the waveguide coincides with the horizontal plane in which *z*=−*L*. An internal interface in the waveguide is described by the equation *z*=−*H*(*x*,*y*), where the function *H*(*x*,*y*) is differentiable and 0<*H*(*x*,*y*)<*L* for all *x* and *y*. That is, the interface is smooth, having a finite slope at every point, and it does not make contact with the surface or the bottom. The waveguide is filled with fluid above the interface and below it. The local sound speed *c*(*x*,*y*,*z*) and the local mass density *ρ*(*x*,*y*,*z*) are differentiable in each layer, but they may be discontinuous at the interface. For simplicity, loss is omitted from this model, but it can be included with little change to the theory [2].

Sound in the waveguide depends harmonically on time *t* with angular frequency *ω*. The acoustic pressure is the real part of the product of the complex pressure *p*(*x*,*y*,*z*) and the phase factor *p*(*x*,*y*,*z*) satisfies the reduced wave equation
*x*,∂/∂*y*,∂/∂*z*). The pressure vanishes at the surface and at the bottom, and it is continuous at the interface. The component of fluid velocity normal to the interface also is continuous at the interface. Now the vector **n**=(∂*H*/∂*x*,∂*H*/∂*y*,1) is normal at the interface, and fluid velocity associated with sound in the waveguide is proportional to *ρ*^{−1}∇*p*. Therefore, the quantity *ρ*^{−1}**n**⋅∇*p*, which equals

The coupled-mode equations derived in this section are extensions of the equations that were originally derived in [3] for a two-dimensional waveguide. The method in [4] is used to cope with the divergence of certain expansions at the interface.^{2} It involves taking moments of equation (2.1) with respect to the eigenfunctions that correspond to the local modes in the waveguide.

The eigenfunctions *f*_{n} and the corresponding local wavenumbers *ξ*_{n} are solutions of the eigenvalue problem
*f*_{n} vanish at the surface and at the bottom and that *f*_{n} and *ρ*^{−1}∂*f*_{n}/∂*z* be continuous at the interface. The squared wavenumbers *n* increases through the integers 1,2,3,…. The eigenfunctions are orthonormal in the sense that *f*_{n} and *ξ*_{n} are implicitly functions of *x* and *y*. It may be assumed that *f*_{n} is real-valued and that the sign of ∂*f*_{n}/∂*z* at the surface is independent of *x* and *y*. Under the assumptions about the interface and the material properties of the fluids in the waveguide, this is enough to make *f*_{n} a differentiable function of *x* and *y*.

The functions *p*, ∂*p*/∂*x* and ∂*p*/∂*y* have the following eigenfunction expansions:
*x* and *y*. Since *p* is continuous everywhere, the expansion for *p* converges uniformly in the interval −*L*≤*z*≤0. Since ∂*p*/∂*x* and ∂*p*/∂*y* may be discontinuous at the interface, the expansions for them may diverge there. However, these expansions always can be integrated termwise against piecewise continuous functions.

Now consider the moments of equation (2.1). Multiply this equation by *f*_{n} and integrate the result from *z*=−*L* to *z*=0:
*F*(*z*) let
*F*(*z*) at *z*=−*H*. The partial derivative ∂/∂*x* can be moved out of the first integral in equation (2.10) as follows:
*y* can be moved out of the second integral the same way. Integrating the third integral by parts twice and using the conditions on *p* and *f*_{n} at the boundaries and equation (2.3) to simplify the result, one finds that
_{n} combines the contributions from all discontinuities at the interface:
_{n}=0 because the quantity defined in equation (2.2) and *f*_{n} are continuous at the interface. Now one can use the definitions of the expansion coefficients in equations (2.7)–(2.9) to write equation (2.14) as follows:
*n*th moment of equation (2.1) becomes the following first-order equation:

This equation must be supplemented with equations for the partial derivatives ∂*a*_{n}/∂*x* and ∂*a*_{n}/∂*y*. They are obtained by differentiating equation (2.7) and the normalization integral *x* yields
*p* from equation (2.4) can be substituted in the jump term as well as in the integral in this equation because it is uniformly convergent. Using equation (2.8) and collecting terms, one finds that
*x* yields

Equations (2.17), (2.21) and (2.22) are the first-order coupled-mode equations for the three-dimensional waveguide. It is convenient to state them in terms of the following coupling coefficients:
*n*=1,2,…,

## 3. A scattering problem with cylindrical symmetry

A cylindrically symmetric waveguide is used to model the seamount and the ocean environment beyond it (figure 1). The axis of symmetry coincides with the vertical coordinate axis. The interface in the waveguide represents the surface of the seamount and the sea floor. The bottom of the waveguide has no physical counterpart. Its only purpose is to make the spectrum of wavenumbers be entirely discrete; for this reason, it is called a ‘false’ bottom.

The time-harmonic point source is suspended in the ocean but not above the seamount, and the ocean environment beyond the seamount is assumed to be horizontally stratified. The seamount is assumed to have a plateau. This feature prevents a singularity that would appear in the mathematical formulation of the scattering problem if the seamount had a pointed peak. Under these assumptions, mode coupling can occur only over the base of the seamount, and the source-free equations (2.25)–(2.27) govern it.

Coupled-mode equations in cylindrical coordinates are obtained from equations (2.25)–(2.27) with the aid of some two-dimensional vectors that are constructed from expansion coefficients and coupling coefficients. Let **w**_{n}=(*b*_{n},*c*_{n}) and **e**_{mn}=(*B*_{mn},*C*_{mn}) for *m*,*n*=1,2,…. Note that
_{∥}=(∂/∂*x*,∂/∂*y*) is the horizontal gradient operator.

Let *r* and *θ* be the radial and angular coordinates, and let *u*_{n} and *v*_{n} are the modal expansion coefficients of the radial and angular components of the horizontal pressure gradient ∇_{∥}*p*. Since the eigenfunctions in a cylindrically symmetric waveguide have no angular dependence, it follows from equation (3.2) that

Equations (2.25)–(2.27) can be expressed in terms of these vectors as follows:
_{∥}⋅**w**_{n}=∂*u*_{n}/∂*r*+*r*^{−1}*u*_{n}+*r*^{−1}∂*v*_{n}/∂*θ* and **w**_{m}⋅**e**_{mn}=*u*_{m}*E*_{mn}, equation (3.4) is equivalent to
*v*_{n} in equation (3.6):

The first step towards a solution is to separate variables in these equations. The source is brought into the analysis at this stage (figure 2). Let the angle *θ* be measured with respect to the horizontal line from the axis of symmetry to the source. Since the symmetry of the problem makes the pressure and its radial derivative be even functions of *θ*, their modal expansion coefficients can be represented by Fourier cosine series in this angle:
*a*_{ln}(*r*) and *u*_{ln}(*r*) are called the ‘azimuthal harmonics’. A system of coupled ordinary differential equations governs the azimuthal harmonics with a given harmonic index *l*. They are derived from equations (3.7) and (3.9) by multiplying these equations by *θ*=0 to *θ*=2*π*. The identity *l*=0,1,…,

The second step is to construct a set of azimuthal harmonics for the incident field from the source. Beyond the seamount, the eigenfunctions depend only on depth and the wavenumbers are constant:
*R*_{1} is the radius of the base of the seamount. Let *R*_{S} be the horizontal distance of the source from the axis of symmetry and let *z*_{S} be the depth of the source. If the seamount were absent, then the acoustic pressure at any point in the waveguide would have the following eigenfunction representation [5]:
*ρ*^{0}(*z*_{S}) is the density at the source. The Fourier cosine series for the factor *ϵ*_{0}=1 and *ϵ*_{l}=2 for *l*=1,2,…. This formula is absolutely convergent if *r*<*R*_{S} and can be differentiated termwise with respect to *r*. Therefore,

The third step is to define another set of azimuthal harmonics for the scattered field:
*r*>*R*_{1}, then

The scattered field beyond the seamount consists of waves that propagate or decay in the radially outward direction. Therefore, the residuals must have the following analytic forms if *r*>*R*_{1}:
*σ*_{ln} is a constant. A boundary value at the perimeter of the base of the seamount determines *σ*_{ln}:

A system of coupled ordinary differential equations governs the residuals with a given harmonic index. These equations are derived from equations (3.12) and (3.13) by substitution of equations (3.25) and (3.26). Equations (3.23) and (3.24) are used to simplify the results. Thus, for *l*=0,1,…,
*ϕ*_{ln}(*r*) and *ψ*_{ln}(*r*) are driven by the references:

Now the scattering problem for the seamount is formulated as a set of two-point boundary value problems for the residuals. For each harmonic index, equations (3.35) and (3.36) are to be integrated over the interval 0<*r*<*R*_{1} subject to radiation conditions at *r*=*R*_{1} from equation (3.34) and to ‘regularity’ conditions that require the residuals to be finite at *r*=0. Using the calculated values of the residuals at the perimeter of the base, one can evaluate the scattered field beyond the seamount with equations (3.27)–(3.33).

Direct numerical integration of equations (3.35) and (3.36) usually is unstable due to extraneous solutions that are singular at the axis of symmetry. Although the regularity conditions are supposed to suppress them, they are excited in discrete calculations by unavoidable numerical errors, which in turn are amplified by them. Therefore, it is necessary to convert these two-point boundary value problems into equivalent problems that can be solved numerically in a stable way.

## 4. Conversion of the two-point boundary value problems

The conversion process discussed in this section is based on ‘variation of parameters’ [7]. Special functions are defined that describe radiating and standing wave solutions of equations (3.35) and (3.36) with the coupling and forcing terms omitted. The residuals are represented as combinations of these functions multiplied by some scaling factors and ‘parameters’ that replace the residuals as unknown functions. New systems of differential equations are derived that govern the parameters. A Riccati transformation is proposed as a stable method for decoupling finite-dimensional truncations of the two-point boundary value problems for the new systems.

### (a) The functions *H*_{ln} and *J*_{ln}

The special functions that are used in the conversion are denoted *H*_{ln}(*r*) and *J*_{ln}(*r*), where *l* is a harmonic index and *n* is a modal index. The function *H*_{ln}(*r*) is the solution of the ordinary differential equation
*r*<*R*_{1} that satisfies the terminal conditions
*J*_{ln}(*r*) is a solution of the same equation
*r*<*R*_{1} that is regular at *r*=0. It should be kept in mind that the local wavenumber *ξ*_{n} appearing in equations (4.1) and (4.4) is implicitly a function of the radius *r* that generally is not constant. For this reason, these equations generally are not equivalent to a Bessel equation, and their solutions generally are not equivalent to cylinder functions.

The plateau enters the analysis at this stage. It is assumed that the wavenumbers are constant over the plateau. Let *R*_{0} be the radius of the plateau, and let
*r*<*R*_{0}, and the initial condition on *J*_{ln}(*r*) implies that
*γ*_{ln} is a constant. With no loss of generality, this constant is assumed to be a real number different from zero, which makes *J*_{ln}(*r*) real-valued for all *r*. Let *H*_{ln}(*r*). This function also is a solution of equation (4.1) because *H*_{ln}(*r*) and *H*_{ln}(*r*) and *J*_{ln}(*r*) is a linear combination of these functions, which can be expressed in the following form because *J*_{ln}(*r*) is real-valued:
*A*_{ln} is a complex constant. This constant is determined up to a multiplicative factor of ±1 by the constraint that |*A*_{ln}|=1.

The reciprocal of the Wronskian of *J*_{ln}(*r*) and *H*_{ln}(*r*) is required in the next stage of the conversion process. The Wronskian is found from equations (4.7), (4.8). Since |*A*_{ln}|=1, the reciprocal is

Computation of the functions *H*_{ln}(*r*) and *J*_{ln}(*r*) is discussed briefly below. It is assumed in this discussion that the factor *r*_{ln} if *l*≠0. This assumption is satisfied if the squared wavenumbers are non-decreasing functions of radius, which is true for many physically reasonable models. It also is assumed that *r*_{ln}<*R*_{1}, which is equivalent to the restriction that *l*<*k*_{n}*R*_{1} if

The function *H*_{ln}(*r*) is computed by integrating equation (4.1). The integration is based on the following representation of *H*_{ln}(*r*). Let
*P*_{ln}(*r*) and *Q*_{ln}(*r*) are real-valued functions. This representation is possible because *H*_{ln}(*r*) does not vanish, which is a consequence of equation (4.7). These functions are solutions of a pair of coupled first-order differential equations,
*r*<*R*_{1}. They satisfy the terminal condition
*P*_{ln}+*iQ*_{ln}=(*dH*_{ln}/*dr*)/*H*_{ln}.

In practice, equations (4.11) and (4.12) are integrated simultaneously, and *H*_{ln}(*r*) is computed from equation (4.10) as the integrations progress. If *l*≠0, then the behaviour of *H*_{ln}(*r*) changes from oscillatory for *r*>*r*_{ln} to exponential for *r*<*r*_{ln}. Numerical integration of equation (4.12) becomes inaccurate as *r* decreases below *r*_{ln}, but this equation can be integrated analytically. The integral is obtained by substituting equation (4.10) in equation (4.7):
*r* decreases, *Q*_{ln}(*r*) rapidly becomes negligible, but *P*_{ln}(*r*) is singular at the axis of symmetry:

If *l*=0, or if *l*≠0 and *r*_{ln}<*R*_{0}, then there is no need to integrate equation (4.4) to compute *J*_{ln}(*r*) for *r*>*R*_{0}. It can be computed from equation (4.8) after the function *H*_{ln}(*r*) has been computed and the constant *A*_{ln} has been determined. The basis for determining *A*_{ln} is the identity
*R*_{0}. In the usual case that *J*_{l}(*χ*_{n}*R*_{0})≠0, it follows from equation (4.16) that |*A*_{ln}|=1 if
*S*_{ln}(*R*_{0}) is found from equation (4.17) with equation (4.6):
*J*_{ln}(*R*_{0}) also specifies the sign of *A*_{ln}:
*A*_{ln}=±*i*|*H*_{ln}(*R*_{0})|/*H*_{ln}(*R*_{0}) in the exceptional case that *J*_{l}(*χ*_{n}*R*_{0})=0.

If *l*≠0 and *R*_{0}<*r*_{ln}, then *J*_{ln}(*r*) cannot be computed reliably for *R*_{0}<*r*<*r*_{ln} from equation (4.8). In this range, *J*_{ln}(*r*) rapidly tends to zero as *r* decreases, but |*H*_{ln}(*r*)| rapidly increases. All significant digits in computed values of the expression Re{*A*_{ln}*H*_{ln}(*r*)} rapidly are lost in finite-precision arithmetic for this reason. However, the procedure described above can be modified to compute *J*_{ln}(*r*) for *R*_{0}<*r*<*r*_{ln}.

This modification is based on the fact that *S*_{ln}(*r*) is a solution of the differential equation
*S*_{ln}(*R*_{0}) from equation (4.19) as the initial value, this equation is integrated from *R*_{0} to *r*_{ln} to compute *S*_{ln}(*r*_{ln}). *J*_{ln}(*r*_{ln}) and *A*_{ln} are computed with *r*_{ln} in place of *R*_{0} in equations (4.18) and (4.20):
*r*>*r*_{ln}, *J*_{ln}(*r*) is computed from equation (4.8), but for *R*_{0}<*r*<*r*_{ln} *J*_{ln}(*r*) is computed from the identity
*S*_{ln}(*r*) that are stored while equation (4.21) is integrated.

### (b) Variation of parameters with scaling

The next stage of the conversion process introduces pairs of ‘parameters’ *p*_{ln} and *q*_{ln} that are used to represent the residuals *H*_{ln}| and its reciprocal do not appear in the usual formulation of this method, but they are essential in the present application of it.^{3} The reason for this is given below.

Equations for *dp*_{ln}/*dr* and *dq*_{ln}/*dr* are obtained by substituting equations (4.25) and (4.26) in equations (3.35) and (3.36). These equations are simplified with equations (4.1) and (4.4) and are easily solved for the derivatives. Using equation (4.9) and the fact that *d*|*H*_{ln}|/*dr*=*P*_{ln}|*H*_{ln}|, one can write the solutions as follows:

The quantities defined in equations (4.29)–(4.33) are regular at *r*=0 because the factors
*r*=0. (The scaling factors were included in equations (4.25) and (4.26) to accomplish this.) However, the factor *P*_{ln} in equations (4.27) and (4.28) is singular at *r*=0, as indicated by equation (4.15). This factor prevents these equations from being integrated in a neighbourhood of *r*=0 by standard numerical methods.

The plateau enters the analysis again at this stage. If the coupling coefficients *E*_{ij} vanish over the plateau, then equations (4.27) and (4.28) can be integrated analytically for 0<*r*<*R*_{0}. This assumption, and the previous assumption that the wavenumbers are constant over the plateau, are valid if the fluids in the waveguide are horizontally stratified above the plateau and below it. They also are approximately satisfied in most other cases where *R*_{0} is small, for example, if *R*_{0} is less than a wavelength.

Under these assumptions, equations (4.27) and (4.28) have the following general solutions for 0<*r*<*R*_{0}:
*α*_{ln} and *β*_{ln} are constants, and
*r*<*R*_{0}:

The constant *α*_{ln} is determined by the condition that the residuals be regular at *r*=0. The only terms in equations (4.38) and (4.39) that are singular at *r*=0 are the terms that are proportional to *α*_{ln}. Therefore, the residuals are regular at *r*=0 if and only if *α*_{ln}=0.

The integral in equation (4.35) can be evaluated in closed form because the factors *J*_{ln}(*τ*) and *ϕ*_{ln}(*τ*) in the integrand are proportional to *J*_{l}(*χ*_{n}*τ*) and *J*_{l}(*k*_{n}*τ*), respectively. In particular, using equations (3.19), (4.6) and (4.37), one obtains the following formula from equation (4.35) for the value of *p*_{ln} at the perimeter of the plateau:

The value of *q*_{ln} at the perimeter of the base of the seamount is determined by the radiation condition on the residuals. It follows from equations (4.2), (4.3) and equations (4.25), (4.26) that equation (3.34) is satisfied if and only if *q*_{ln}(*R*_{1})=0.

Therefore, the two-point boundary value problems for the residuals are equivalent to the following two-point boundary value problems for the parameters. For each harmonic index *l*, the system of equations (4.27) and (4.28) (*n*=1,2,…) is to be integrated over the interval *R*_{0}<*r*<*R*_{1} subject to initial values *p*_{ln}(*R*_{0}) from equation (4.40) and to terminal values *q*_{ln}(*R*_{1})=0.

### (c) Decoupling by the Riccati transformation

The two-point boundary value problems for the parameters involve an infinite number of modes. In practice, only a finite number of modes can be retained in these problems. Let *N* be the number of modes that are retained. Let **p**_{l}, **q**_{l}, **f**_{l} and **g**_{l} be the *N*-dimensional column vectors with components *p*_{ln}, *q*_{ln}, *f*_{ln} and *g*_{ln} in the *n*th rows, respectively. Let **M**_{l}, **B**_{l} and **C**_{l} be the *N*×*N* matrices with elements *M*_{lnm}, *B*_{lnm} and *C*_{lnm} in the *n*th rows and in the *m*th columns, respectively. Let **A**_{l} be the *N*×*N* diagonal matrix that has *A*_{ln} as the *n*th diagonal element, and let **P**_{l} be the *N*×*N* diagonal matrix that has *P*_{ln} as the *n*th diagonal element. Note that these vectors and matrices are functions of the radius, except **A**_{l}, which is a constant matrix.

For each harmonic index *l*, the system of equations that consists of the truncated versions of equations (4.27) and (4.28) for *n*=1,…,*N* can be written in vector form as follows:
*R*_{0}<*r*<*R*_{1}. The components of **p**_{l}(*R*_{0}) are specified by equation (4.40) and **q**_{l}(*R*_{1})=**0**.

The coupling between equations (4.42) and (4.43) complicates this task. However, these equations can be decoupled by a Riccati transformation [8].^{4} This transformation uses new *N*-dimensional column vectors **x**_{l} and **y**_{l} to represent **p**_{l} and **q**_{l} as follows:

The *N*×*N* matrices **T**_{l} and **S**_{l} in these equations satisfy differential equations that are designed to make the vectors **x**_{l} and **y**_{l} be independent of one another. **T**_{l} is the solution of the terminal value problem
**T**_{l}(*R*_{1})=**0**. **S**_{l} is the solution of the initial value problem
**S**_{l}(*R*_{0})=**0**.

Substituting equations (4.44) and (4.45) in equations (4.42) and (4.43) yields equations for the vectors **x**_{l} and **y**_{l} that are uncoupled. **x**_{l} is the solution of the initial value problem
**x**_{l}(*R*_{0})=**p**_{l}(*R*_{0}). This initial condition is implied by equation (4.44) because **S**_{l}(*R*_{0})=**0**. **y**_{l} is the solution of the terminal value problem
**y**_{l}(*R*_{1})=**0**. This terminal condition is implied by equation (4.45) because **q**_{l}(*R*_{1})=**0** and **T**_{l}(*R*_{1})=**0**.

Equations (4.46)–(4.49) are the final results of the conversion process. They have no extraneous solutions that can destabilize numerical integration of the boundary value problems. However, it is useful to make a simple change of variables in equations (4.46) and (4.47) that reduces the computational cost of integrating them. Let
**B**_{l} and **C**_{l} are symmetric, it follows from equation (4.46) that *Θ*_{l} is the solution of the terminal value problem
*Θ*_{l}(*R*_{1})=**0**. Because of the boundary condition and the symmetry of **K**_{l} and **C**_{l}, *Θ*_{l} is symmetric. Similarly, it follows from equation (4.47) that *Π*_{l} is the solution of the initial value problem
*Π*_{l}(*R*_{0})=**0**. Because of the boundary condition and the symmetry of **K**_{l}, *Π*_{l} also is symmetric. Since *Θ*_{l} and *Π*_{l} are symmetric, equations (4.52) and (4.53) can be integrated with about half as many arithmetic operations as equations (4.46) and (4.47).

It is convenient to write equation (4.48) in terms of these matrices as follows:

The integrations are done in the following order. First, equation (4.52) is integrated from *r*=*R*_{1} down to *r*=*R*_{0}. The values of *Θ*_{l} that are computed at intermediate integration steps are stored for later use. Next, equation (4.53) is integrated from *r*=*R*_{0} up to *r*=*R*_{1}, and equation (4.54) is integrated in parallel with it. These integrations use the stored values of *Θ*_{l}.

Integration of equation (4.54) yields the value of **x**_{l}(*R*_{1}). Since *q*_{ln}(*R*_{1})=0, equations (4.2) and (4.25) imply that
**p**_{l}(*R*_{1})=**x**_{l}(*R*_{1}) because **y**_{l}(*R*_{1})=**0**. Therefore, if *x*_{ln} is the component in the *n*th row of **x**_{l}, then

The solution of equation (4.54) depends on the source. If the source is moved, then this equation must be integrated again in order to evaluate the scattered field. The following section describes an alternative formulation of the coupled mode solution that requires exactly *N* integrations of equation (4.54) but allows the scattered field to be evaluated for an arbitrary number of source positions. This is the most efficient formulation if more than *N* source positions are considered.

## 5. Algebraic formula for the scattered field and the scattering matrices

The references *j* is the index of the mode that is associated with the standing wave. The references are expanded as follows:

Let

The operations in equations (3.27), (3.28), (3.31), (3.33) and (5.5) can be rearranged to produce the following formula for the scattered field beyond the seamount. If *r*>*R*_{1}, then
*w*_{lnj}=*w*_{ljn}. Therefore, the coefficient *w*_{lnj} can be regarded as the element in the *n*th row and in the *j*th column of a symmetric ‘scattering matrix’ **W**_{l}. The symmetry of **W**_{l} is consistent with the requirement that equations (5.6) and (5.7) satisfy the reciprocity principle.

In practice only a finite number of terms are retained in equations (5.6) and (5.7). Let *N* be the number of modes that are retained. Equation (4.54) is used to compute the boundary values *w*_{lnj} from equation (5.8). This is done as follows. First, for *n*,*j*=1,…,*N*, let

Next, let **F**_{l} be the *N*×*N* matrix with the element *f*_{lnj} in the *n*th row and in the *j*th column, and let **G**_{l} be the *N*×*N* matrix with the element *g*_{lnj} in the *n*th row and in the *j*th column. These matrices appear in the following differential equation for another *N*×*N* matrix **X**_{l}, which is equivalent to *N* copies of equation (4.54) with the columns of **X**_{l}, **F**_{l} and **G**_{l} in place of the vectors **x**_{l}, **f**_{l} and **g**_{l}:
*r*=*R*_{0} to *r*=*R*_{1} subject to initial conditions on the matrix **X**_{l} that are derived from equation (4.40):
*x*_{lnj} is the element in the *n*th row and in the *j*th column of **X**_{l}. Integration of equation (5.13) yields **X**_{l}(*R*_{1}), and by the same reasoning that leads to equation (4.56), one finds that

For each harmonic index, equations (4.52), (4.53) and (5.13) need to be integrated only once because they are independent of the source. The computed coefficients *w*_{lnj} are stored for later use with truncated versions of equations (5.6) and (5.7). These equations are evaluated efficiently with Clenshaw's algorithm for the cosine series and recurrence relations for the Hankel functions.

Equations (5.13) and (5.16) have no intrinsic symmetries to ensure that the computed finite-dimensional approximations to the scattering matrices are symmetric. However, the size of their asymmetric parts provides a good error estimate for the approximations.

## 6. Conclusion

Mode coupling in three-dimensional irregular acoustic waveguides is governed by a system of first-order partial differential equations. These equations can be solved analytically to find the scattered field from a cylindrically symmetric seamount with a smooth surface. There is a simple algebraic formula for the scattered field beyond the seamount. It is characterized by a set of infinite-dimensional matrices that are independent of the source. Finite-dimensional approximations to these matrices can be computed with a stable numerical procedure that is based on the analytical coupled mode solution. This procedure is well suited for parallel implementation because the scattering matrices are computed independently of one another.

## Competing interests

I declare I have no competing interests.

## Funding

I received no funding for this study.

## Appendix A. Symmetry of the scattering matrices

Consider the response of the seamount to a standing wave in the waveguide. This field has the form *l* is the index of an angular harmonic, *k*_{j} is a modal wavenumber and *r*>*R*_{1}:
*w*_{lnj}=*w*_{ljn}. Since *w*_{lnj} is the element in the *n*th row and in the *j*th column of the scattering matrix **W**_{l}, this relation implies that **W**_{l} is symmetric.

Let *V*_{R} be the vertical cylinder of radius *R* that is centred on the axis of symmetry and extends from the false bottom up to the surface of the waveguide. Since the fields *S*_{R} be the vertical surface of the cylinder *V*_{R}. Since the fields vanish at the false bottom and at the surface, the divergence theorem implies that
*j*≠*n* and *R*>*R*_{1}. The surface integral in this equation can be evaluated by substituting the far-field representations of *φ*_{j}=*k*_{j}*R* and *φ*_{n}=*k*_{n}*R*, one finds that
*w*_{lnj}=*w*_{ljn}.

## Footnotes

↵1 The original version of this paper also had a section with examples of scattering matrices and a field calculation. This section now appears as in electronic supplementary material, §8.

↵2 In this paper, coupled-mode equations for electromagnetic waveguides are derived from Maxwell's equations. The abstract explains the basis of the method as follows. ‘The essential point is that a function, which for practical purposes is sufficiently arbitrary, may be represented in numerous ways by a series of orthogonal functions; and that when some such series are non-differentiable, the required relations between the coefficients of the series may be obtained from Maxwell's equations by integration rather than by conventional substitution followed by differentiation’.

↵3 These scaling factors were suggested by the form of the ‘normalized range functions’ that are used in [1]. The notation for the functions

*H*_{ln}and*J*_{ln}was chosen in recognition of this fact.↵4 Additional references to the application of the Riccati transformation to decoupling of two-point boundary value problems can be found in [2].

- Received July 7, 2015.
- Accepted December 22, 2015.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.