## Abstract

An improved version of the multimodal admittance method in acoustic waveguides with varying cross sections is presented. This method aims at a better convergence with respect to the number of transverse modes that are taken into account. It is based on an enriched modal expansion of the pressure: the *N* first modes are the local transverse modes and a supplementary (*N*+1)th mode, called boundary mode, is a well-chosen transverse function orthogonal to the *N* first modes. This expansion leads to the classical form of the coupled mode equations where the component of the boundary mode is of evanescent character. Under this form, the multimodal admittance method based on the Riccati equation on the admittance matrix (the Dirichlet-to-Neumann operator) is straightforwardly implemented. With this supplementary mode, in addition to the improvement of the convergence of the pressure field, results show a superconvergence of the scattered field outside of the varying cross sections region.

## 1. Introduction

Coupled mode theory is a classical method to deal with the propagation of inhomogeneous waveguides [1–3]. Starting from the second-order coupled mode equations, the multimodal admittance method is based on the introduction of the admittance matrix, which is the modal representation of the Dirichlet-to-Neumann operator on the transverse section of the waveguide [4–6]. The admittance matrix is governed by a matrix Riccati differential equation that can be integrated from the radiation condition at the outlet of the waveguide towards the source position. Then, it remains only to account for the source to get an initial condition on the pressure at the inlet that can be used to calculate the field in the whole waveguide by integrating a first-order differential equation. This method allows us to avoid the problem of numerical divergence owing to the presence of evanescent modes. Besides, when only the scattering matrix is needed, one advantage of the whole method is that the admittance matrix does not have to be stored during the integration of the initial value problem starting from the radiation condition. As any method used to solve the coupled mode equations, the multimodal admittance method has a rate of convergence with respect to the number of local transverse modes in the series expansion of the pressure at each location on the axis of the waveguide. In the literature, several propositions have been made to improve this rate of convergence [7–11]. Similar to the classical attachment mode used in structural mechanics [12–14], all these techniques use a boundary mode which is not a local transverse mode but that encapsulates the less convergent part of the series. Nevertheless, these techniques are not well suited to implement the multimodal admittance method for varying cross section waveguides. This is because it is not possible to define a radiation condition on the boundary mode.

Following the results presented in a recent paper [15], an improved version of the multimodal admittance method is presented, using the idea of a boundary mode. It is based on an enriched modal expansion of the pressure, as in [7,10,11,15]. The pressure is projected on *N*+1 transverse functions associated with *N*+1 components. The first *N* transverse functions are chosen as the local transverse modes, corresponding to the solutions in a waveguide with a constant cross section equal to the local cross section. The choice of the supplementary (*N*+1)th transverse function is the key point of our improved method. As in [7,10,11,15], it is chosen to absorb the less convergent part of the series representing the pressure, but the novelty is that we take it to be orthogonal to the first *N* local transverse modes. When projecting the wave equation on these (*N*+1) transverse functions, this latter property ensures that the coupled mode equations take the following classical form: the identity matrix appears in front of the higher order differential terms, and more importantly the boundary mode appears as an evanescent mode. Regarding the multimodal admittance method, this allows us to follow the usual formalism by defining a radiation condition at the outlet for this boundary mode; as this mode is evanescent, it is sufficient to impose exponentially decreasing behaviour (as with any evanescent mode). Using this property, it is then straightforward to implement the multimodal admittance method.

The paper is organized as follows. In §2, we present the formulation of the improved multimodal admittance method for waveguides with varying cross sections. It is shown that, when the boundary mode is introduced, the coupled mode equations with a boundary mode, equation (2.7), keep the same form as the coupled mode equations without boundary mode. Thus, the whole formalism of the multimodal admittance method is conserved. Both the two- and the three-dimensional axisymmetric cases are considered. Section 3 focuses on the convergence of the improved modal method compared with the classical modal method. We show a much better improvement in the convergence of the scattering coefficients than on the pressure field, leading to what may be called a superconvergence of the scattering. Preceding the concluding remarks, a discussion on previous works implementing a modal method with a boundary mode is presented in §4.

## 2. Improved admittance method

The propagation of linear acoustic waves in a waveguide with rigid walls is considered in the harmonic regime. The time dependance e^{−iωt} is omitted in the following. For the sake of comprehensiveness, the two- and three-dimensional axisymmetric cases are detailed in the following.

### (a) Two-dimensional case

We start from the Helmholtz equation
2.1with Neumann boundary condition on the lower straight wall ∂_{y}*p*(*x*,0)=0 and on the varying upper wall ∂_{n}*p*(*x*,*h*(*x*))=0 (figure 1). Defining *u*≡∂_{x}*p*, it can be written as a first-order equation along the *x*-axis of the waveguide
2.2with boundary conditions
2.3

In the spirit of coupled mode equations, expansions of *p* and *u* onto transverse functions *φ*_{n}(*y*;*x*) are written as
2.4The *N* functions *φ*_{n}, 0≤*n*≤(*N*−1), are the usual local transverse modes. They correspond to the eigenfunctions of the transverse eigenproblem , with homogeneous Neumann boundary condition on the walls ∂_{y}*φ*_{n}(0;*x*)=∂_{y}*φ*_{n}(*h*(*x*);*x*)=0 (the associated eigenvalue is *γ*_{n}=*nπ*/*h*)
2.5Here, for the improved method, the (*N*+1)th function *φ*_{−1} is not a local transverse mode but it will be chosen as a function orthogonal to the modes *φ*_{n}, 0≤*n*≤(*N*−1), with the particular boundary conditions ∂_{y}*φ*_{−1}(*h*(*x*);*x*)≠0 and ∂_{y}*φ*_{−1}(0;*x*)=0. This choice is meant to better account for boundary condition (2.3) at *y*=*h* that cannot be satisfied by the local transverse modes. It is very easy to write such a boundary mode for any order of truncation *N*
2.6where denotes the scalar product (where overbar means conjugate). Any *χ* with ∂_{y}*χ*(*h*(*x*);*x*)≠0 and ∂_{y}*χ*(0;*x*)=0 is possible, and in the following we have chosen
where *a*_{N} is a normalization factor such that ∥*φ*_{−1}∥=1. Then, to obtain the coupled mode equations, we just have to project wave equation (2.2) on the (*N*+1) modes *φ*_{n}. Details of the calculations are presented in appendices A and B. Eventually, the coupled mode equations take the form
2.7with **p**≡(*p*_{n}), **u**≡(*u*_{n}), and where I is the identity matrix and K is a diagonal matrix with diagonal element given by
2.8and the coupling matrix A has components
2.9The square root is chosen with positive real part or with positive imaginary part. For 0≤*n*≤(*N*−1), correspond to the usual horizontal wavenumbers of the modes in the straight part of the waveguide. is found negative (assuming of course that the truncation includes all the propagative modes), *k*_{−1} is imaginary, representing an evanescent part of the wavefield. As *N* is increased, it is more and more evanescent, as can be seen in the explicit expression obtained from (2.8)
2.10which implies that for large *N* (for details, see [15]).

Coupled mode equations (2.7) with the boundary mode (as well as equations (2.8) and (2.9)) have exactly the same form as the coupled mode equations without the boundary mode [6]. Consequently, from this point, it is straightforward to implement the admittance method following the same steps as in [6]. We first define the admittance matrix Y and the propagator matrix G
2.11where *x*_{f} is the abscissa of the outlet and *x*≤*x*_{f}. The matrix Y satisfies a Ricatti equation and G a differential equation coupled to Y,
2.12Both equations are integrated from the radiation condition at the right part of the waveguide *x*=*x*_{f} with the initial conditions Y(*x*_{f})=Y_{c} (Y_{c} is the characteristic admittance defined by Y_{c}≡*i*K(*x*_{f})) and G(*x*_{f})=I. The above Riccati equation is used because of its stability [4–6], and the numerical integration by Magnus method is detailed in [6]. Once equations (2.12) are integrated up to the inlet *x*=*x*_{i}, the reflection and transmission matrices can be calculated. Indeed, the wave can be split into right- and left-going waves **p**=**p**^{+}+**p**^{−} and **u**=Y_{c}(**p**^{+}−**p**^{−}), so that the reflection matrix is defined by **p**^{−}(*x*_{i})=R**p**^{+}(*x*_{i}) and the transmission matrix is defined by **p**(*x*_{f})=T**p**^{+}(*x*_{i}). Using from equation (2.11) **u**(*x*_{i})=Y(*x*_{i})**p**(*x*_{i}) and **p**^{+}(*x*_{f})=**p**(*x*_{f})=G(*x*_{i})**p**(*x*_{i}), we get
2.13

Note that the great advantage of the method is that the matrices Y and G do not have to be stored during the integration of differential equations (2.12). If needed, the solution inside the scattering region, **p**(*x*_{i}≤*x*≤*x*_{f}), can be obtained from the first equation of coupled mode equations (2.7)
2.14with the incident wave condition at the left part of the waveguide *x*=*x*_{i}.

It is important to remark that the introduction of boundary mode does not change the formalism of the admittance method: equations (2.7), (2.12) and (2.14) have the same form as for the usual case without the boundary mode and, from the structure of coupled mode equations (2.7), the energy conservation is preserved (for details on energy conservation in the multimodal method, see [16]).

### (b) The three-dimensional axisymmetric case

This case corresponds to the acoustic propagation in rigid ducts, with varying circular cross sections (without mean flow). It is similar to the two-dimensional case, with a difference lying only in the form of the transverse Laplacian. We give below only the expressions that are different from the two-dimensional case and all the above remarks for the two-dimensional case remain valid for the three-dimensional axisymmetric case. The Helmholtz equation is
2.15with Neumann boundary condition on the varying wall ∂_{n}*p*(*z*,*h*(*z*))=0. It can be written as a first-order equation along the *z*-axis of the waveguide
2.16with boundary condition
2.17The same form of decomposition as in (2.4) are chosen
2.18where *φ*_{n}(*r*;*z*) (*n*≥0), defined by
2.19are the eigenfunctions of the transverse eigenproblem with Neumann boundary condition at *r*=*h*(*z*), ∂_{r}*f*_{n}(*h*)=0. Here, *γ*_{n}=*β*_{n}/*h* and *β*_{n} is the (*n*+1)th zero of *J*_{1} (and *J*_{n} denotes the *n*th Bessel function). The boundary mode, as in two dimensions, is defined by equation (2.6), with a difference in the definition of the scalar product and *χ* is chosen such that ∂_{r}*χ*(*h*(*z*);*z*)≠0. The function *χ* has been chosen as follows:
where *β* is the first zero of *J*_{0} (with five digits: *β*=2.4048), and, as in the two-dimensional case, *a*_{N} is the normalization factor such that ∥*φ*_{−1}∥=1.

Then, following the same steps as in two dimensions in the previous section, one obtains the same form of coupled mode equations as in (2.7). Only the expressions of the matrices occurring in this equation change: K is a diagonal matrix with diagonal element given by and A_{nm}≡(*φ*_{n},∂_{z}*φ*_{m}). Equations (2.11)–(2.14) remain valid. Details of the calculations are presented in appendix C.

## 3. One-dimensional toy model

In this section, for the sake of clarity, we consider very simple one-dimensional models to examine two aspects that will be useful in the following. The first is concerned with the improved convergence using boundary modes (a technique which is also called polynomial subtraction [17–19]). The second allows us to illustrate the different types of errors that occur when solving an ordinary differential equation with truncated series.

### (a) Improved convergence in one dimension

Let us take a function *f* with sufficient regularity defined on [0,1]. We want to project *f* on the Neumann modes *φ*_{n} defined by with *φ*_{n}′(0)=*φ*_{n}′(1)=0. These modes *φ*_{n} are similar to the ones considered in the waveguide problem. The explicit expression of the normalized *φ*_{n} is with *γ*_{n}≡*nπ*. Here, the scalar product is defined by . The expansion of *f* is
3.1Integrations by parts give
3.2with , leading to
3.3Hence, successive integrations by part give the behaviour of the coefficients *p*_{n}
3.4With *γ*_{n}∝*n*, the convergence of the coefficients *p*_{n} is 1/*n*^{2} if *f*′ does not vanish at the boundaries and it is at least 1/*n*^{4} if it vanishes. The convergence of the truncated series
3.5is given by the rest
3.6whose norm is
3.7If *f*′ does not vanish at the boundaries, we have
3.8

More generally, from equation (3.4), it is seen that the convergence of the series is governed by the values of the successive derivatives of the function *f* at the boundaries. Therefore, it is possible to improve the convergence by introducing a new function *g*
3.9with *χ*′_{0}(1)=*χ*_{1}′(0)=0 and *χ*′_{0}(0)=*χ*_{1}′(1)=1. The function *g* satisfies by construction *g*′(0)=*g*′(1)=0, so that, from (3.4), its coefficients *q*_{n}≡(*φ*_{n},*g*) behave as 1/*n*^{4}. This polynomial subtraction technique (because *χ*_{0} and *χ*_{1} can be chosen as polynomials) improves the convergence of the series. Indeed, considering the truncated series
3.10the following behaviours are obtained:
3.11It follows that a truncated series approximation for *f* can be built, as
3.12that satisfies
3.13It is the basic idea of the boundary mode method: in equation (3.12) is a much better approximation of *f* than *f*_{N} in equation (3.5) owing to the introduction of the functions *χ*_{0} and *χ*_{1} (known explicitly).

Let us take the concrete example of on [0,1], satisfying *f*′(0)=0 and *f*′(1)≠0. From (3.4), one knows that *p*_{n}∝1/*n*^{2}. Here, the explicit expression of *p*_{n} is
3.14and one recovers that *p*_{n}∝1/*n*^{2}. To improve the series expansion, the following choice for *χ* is done: *χ*(*y*)≡*χ*_{1}(*y*)=*y*^{2}/2 (here, *χ*_{0} is useless for *f*′(0)=0) and we define *g*≡*f*−*f*′(1)*χ*, as in (3.9). The coefficients *q*_{n}=(*φ*_{n},*g*) are explicit
3.15for *n*>0, where
3.16Eventually, the truncated series
3.17converges as
3.18We see that it is very simple to improve the convergence by the technique of polynomial subtraction when *f*′(1) is known. In this paper, the improved multimodal method uses the same kind of ideas in a more involved situation, where one has to solve a partial differential equation with *f*′(1) unknown. The procedure consists in introducing an extra degree of freedom *p*_{−1}, replacing *f*′(1) in (3.17), and supposed to approximate *f*′(1). As another remark, note that decomposition (3.17) can be rewritten as
3.19This choice has the advantage of dealing with orthogonal functions, the *φ*_{n} (*n*≤*N*) and the function (what we call loosely a boundary mode). This is the same form as the expansion chosen in equation (2.4). This makes the projection easier and this is the choice made for the improved multimodal method in the waveguide.

### (b) Boundary value problem: remainder error and coefficient error

Let us consider the following boundary value problem:
3.20with
3.21where and , such that is the known exact solution. The numerical solution is sought as an *N*th partial sum approximation
3.22with *φ*_{n}(*y*) the same Neumann modes as in the previous section. To obtain the equations on , one uses the Galerkin method by projecting equation (3.20) on the *φ*_{n} (for *n*≤*N*). One gets
3.23where the following properties have been used: *f*′(0)=0, (*φ*_{n},*f*)=*p*_{n}, and the integration by parts
Then, it is sufficient to use the boundary condition at *y*=1 to obtain the linear system governing the
3.24for 0≤*n*≤*N*. It is a system of (*N*+1) equations with (*N*+1) unknowns that can be solved by a standard inversion technique. The error between the exact solution *f*(*y*)
3.25and the numerical Galerkin solution (3.22) can be written as the sum of two parts
3.26The second part in the left-hand side, denoted *R*_{N}(*y*), is the remainder error corresponding to the distance between *f* and the set spanned by the functions *φ*_{0},…,*φ*_{N}. The first part in the left-hand side is the coefficient error *r*_{N}: let us see where does it comes from. Equation (3.23) is exact but the elimination of *f*′(1) using is not exact. It is possible to write an exact system by using . More precisely, one writes
3.27Thus, from (3.23) and (3.27), the exact coefficients *p*_{n} (0≤*n*≤*N*) satisfy the exact linear system of (*N*+1) equations
3.28It appears that the coefficients and *p*_{n} in the coefficient error (3.26) obey two different systems of equations (3.24) and (3.28). These two systems can be written formally as
3.29with , *S*_{n}≡−*βφ*_{n}(1), *T*_{n}≡−*αφ*_{n}(1)*R*_{N}(1) and with **p**=(*p*_{0},…,*p*_{N})^{T}, . Thus, the coefficient error is
3.30In this specific example, the exact solution is , such that *T* is known and the coefficient error in (3.26) can be calculated by numerical inversion of *M*. The convergence is calculated and it can be summarized as follows:
3.31where the coefficient error and the remainder error have been calculated with the norm *L*^{2}. In this simple example, the coefficient error *r*_{N} is smaller than the remainder error *R*_{N}. For the waveguide problem, it will be shown that this is not always the case.

## 4. Results

In this section, the improved admittance multimodal method is applied for the three-dimensional axisymmetric case. The computation performed following the steps presented in §2. We have checked that the same conclusions are drawn for the two-dimensional case.

To begin with, it is shown that the use of the impedance matrix method implemented with the boundary mode allows us to reach a high accuracy with very few degrees of freedom. A typical example is shown in figure 2 at a frequency (*kh*_{0}=7, *h*_{0}=1 being the height of the right lead) just below the cut-on frequency of the third mode in the right lead, so that there are two propagating modes in the right lead and one propagating mode in the left lead. The method is performed with *N*_{m}=*N*+1=4 degrees of freedom, achieving an accuracy better than 1% on the scattering coefficient. The incident wave is mode 0 from the left and the outgoing boundary condition is imposed at the right.

In the following, the convergence of the improved method is examined with respect to the wave field and with respect to the scattering coefficients, where a superconvergence will be shown. To do that, the geometry of a cosine waveguide is considered, namely
4.1Except for the result shown in figure 7, the frequency is *kh*_{0}=0.9*π* with *h*_{0}=1, which corresponds to one propagating mode in the leads. Moreover, the mode 0 is incident from the left and the outgoing boundary condition is imposed at the right.

Let us first examine the asymptotic behaviour of the modal coefficients. Equations (2.4) and (2.6) show that
4.2with
4.3Expansions (4.2) and (2.4) correspond, respectively, to expansions (3.17) and (3.19). As seen in the toy model section, the coefficients *q*_{n} are expected to decrease more quickly than the coefficients *p*_{n} (equations (3.11) and (3.8)). In our impedance multimodal method, the system is solved on the coefficients *p*_{n} that are found to follow *p*_{n}∝1/*n*^{2} (figure 3). Nevertheless, when *q*_{n} is computed from equation (4.3), the behaviour *q*_{n}∝1/*n*^{4} is actually recovered [7,11,15].

Next, the local and global errors on the field are defined by
4.4and
4.5where *p*(*z*,*r*) is the solution calculated at order *N* (equation (2.4), i.e. with *N*_{m}=*N*+1 coefficients) and *p*^{ex}(*z*,*r*) the converged solution considered to be exact. In order to assess the advantage of the boundary mode, the usual multimodal method without the boundary mode is considered, using the expansion
4.6whose corresponding local and global errors are defined as in equations (4.4) and (4.5). Note that above expression equation (4.6) has (*N*+1) components on the classical transverse modes *φ*_{n}; it has been chosen as such in order to have the same number of degrees of freedom *N*_{m}=*N*+1.

The rates of convergence of the global errors are shown in figure 4 as a function of the total number of transverse functions taken into account (number of modes plus the boundary mode in the case of the improved method and number of modes in the case of the usual multimodal method). The rate of convergence is 1/*N* for the usual method and 1/*N*^{3.5} for the improved method.

In order to get further insights on the non-uniformity of the convergence, it is interesting to inspect the pattern of the local error in space. The real part of the total fields as well as these local errors in space *ϵ*_{N}(*z*,*r*) are shown in figure 5 both for the usual and the improved methods. The total number of modes is *N*_{m}=10 in both cases. Both results are close to the reference solution but one can note spurious wiggles near the boundaries for the computation without the boundary mode (of course, the amplitudes of these wiggles decrease with the order of truncation [4]). When plotting the solution (figure 5*a*), the wiggles appear to be located close to the inclined wall. It is interesting also to inspect the pattern of the local error for the computation with boundary mode. Indeed, in this case, the pattern of the error presents a characteristic length that corresponds to the first omitted mode. This corresponds to the intuitive idea of the error owing to a truncation for series of functions of one variable, namely if , the error is . A contrario, the pattern of the error without the boundary mode does not follow this simple behaviour. This means that the coefficient error (the modal coefficient depends on *N* and *p*_{n} converges towards its asymptotic value ) has a non-negligible weight, as already discussed in the toy model section.

The next results concern the convergence of the elements of the scattering matrix. These are the most important quantities because usually, one is more interested in obtaining the transmission and reflection coefficients owing to the inhomogeneity of the waveguide rather than the detailed pattern of the field inside the scattering region. Note that, to compute the scattering matrix, only the integrations of the equations on Y and G are needed. They are evolution equations posed as an initial value problem which means that they do not imply large memory storage. In our computation example, with *kh*_{0}=0.9*π*, only one mode is propagative in the straight terminations of the waveguide so that the scattering is characterized by one transmission coefficient *T*=*T*_{00} and one reflection coefficient *R*=*R*_{00}. Figure 6 shows the rate of convergence of the transmission coefficient as a function of the total number of transverse modes. Without the boundary mode, the convergence is |*T*−*T*^{ex}|∼*N*^{−1}, and with the boundary mode the convergence is |*T*−*T*^{ex}|∼*N*^{−4.5}. Moreover, it can be seen that two modes and the boundary mode achieve an excellent accuracy that is reached for a truncation around 40 modes using the classical modal method. This is what we call the superconvergence. As a rule of thumb, this excellent accuracy appears to be obtained as soon as the truncation includes the first evanescent mode. In our example, with one propagative mode, this means that the computation includes the propagative mode 0, the first evanescent mode 1 and the boundary mode. Figure 7 illustrates this accuracy in the whole range of frequency below the first cut-on frequency of the termination *k*_{c}=*γ*_{1}/*h*_{0} (*γ*_{1}≃3.83). The error |*T*−*T*^{ex}| is plotted as a function of the dimensionless frequency *kh*_{0} in three cases: with and without the boundary mode using *N*=2 in (4.2) (modes 0, 1 and the boundary mode) and in (4.6) (modes 0, 1 and 2). When comparing the usual and improved multimodal methods with the same number of transverse modes (a total of *N*_{m}=3 modes in each case), the accuracy is improved by one order of magnitude. Thus, it is observed that the superconvergence of the improved method works uniformly with respect to the frequency.

### (a) On the convergence and the superconvergence

Following the idea of the remainder error and of the coefficient error presented in the toy model section, we examine the partition of the error when using the usual or the improved multimodal method. Let us consider the expressions of the pressure field without the boundary mode (usual method)
4.7with the boundary mode (improved method)
4.8and the converged field, considered as an exact solution
4.9To analyse the convergence of the pressure field, the error |*p*−*p*^{ex}| is split into three parts
4.10and with the boundary mode
4.11with *N*_{p}+1 the number of propagating modes. *ϵ*_{1} refers to the error on the propagative field and *ϵ*_{2}+*ϵ*_{3} refers to the error on the evanescent field. The rates of convergence indicated for each term have been obtained numerically. For *ϵ*_{3}, these rates of convergence can be obtained analytically by using the behaviour of for large *n*, namely and . Also, it appears that without the boundary mode, *ϵ*_{1} and *ϵ*_{2} are dominant, whereas with the boundary mode *ϵ*_{2} and *ϵ*_{3} are dominant.

These rates of convergence are observed in the varying part of the waveguide but in the straight parts of the waveguide, only *ϵ*_{1} subsists there because *ϵ*_{2} and *ϵ*_{3} become exponentially small (they are associated with evanescent modes). Consequently, the error in the straight parts is given by *ϵ*_{1} and the rate of convergence can be different inside and outside the scattering region. Importantly, it turns out that the addition of the boundary mode improves more significantly the rate of convergence of the propagative part of the field (from *N*^{−1} to *N*^{−4.5}) than that of the evanescent part of the field (from *N*^{−1} to *N*^{−3.5}). This explains the superconvergence observed for the scattering coefficients as they are defined from the propagative modes outside the scattering region.

## 5. Coupled formulation of the improved method

Our improved multimodal decomposition is
5.1where *p*_{−1} is found to be an evanescent mode, with an imaginary wavenumber *k*_{−1} that tends to infinity as . Also, it can be noted that an uncoupled system of equations on the (*p*_{n}) is obtained in the straight parts of the waveguide (namely , for *n*=−1,…,*N* that can be written **p**′′+*K*^{2}**p**=**0**).

Previous works on improved multimodal methods used also a boundary mode [7,11]. Nevertheless, the decomposition was written differently, with a boundary mode non-orthogonal to the first *N* ones
5.2For the same function *χ*(*r*;*z*), the two modal decompositions coincide owing to equation (4.3) and *p*_{−1}=*q*_{−1}. The decomposition, equation (5.2), leads to coupled mode equations on **q**≡(*q*_{n})
5.3In the straight parts of the waveguide, C=0 but B and D remain non-diagonal matrices, inducing a coupling that makes the implementation of the admittance method not straightforward, if not problematic. The approach chosen in [7,11] is to pose the coupled mode equations as a boundary value problem with Dirichlet boundary condition on *q*_{−1}. In view of our work, we propose an interpretation of this boundary condition based on an analogy with the perfectly matched layers (PML). Indeed, as *q*_{−1} is evanescent in the straight parts of the waveguide, the straight terminations of the waveguide play the role of PML with *q*_{−1}∼e^{−|k−1|L} decreasing with the size *L* of the PML.

If *L* is large enough, the choice *q*_{−1}=0 at the end of the PML is similar to the outgoing radiation condition for *q*_{−1} that is used in this paper. Also, because |*k*_{−1}|∝*N*, increasing the order of the truncation *N* decreases the necessary size *L* of the PML.

Results in [7,11]—boundary value problem with Dirichlet boundary conditions on *q*_{−1}—have shown an improvement of the convergence inside the scattering region of the same kind as the one we obtain in this paper. It could be interesting to know if, for small *N*, the accuracy is the same as in our method and if the superconvergence for scattering coefficient that we obtain might be also observed with Dirichlet boundary condition on *q*_{−1}.

## 6. Concluding remarks

The introduction of a boundary mode in the multimodal admittance method has been shown to improve significantly the efficiency of the method for acoustic waveguides with varying cross sections. The convergence of the pressure field inside the varying part of the waveguide increased but, more importantly, a superconvergence is obtained on the scattering coefficients, from 1/*N* to 1/*N*^{4.5}. Moreover, this asymptotic rate of convergence is accompanied by an error that is already small at small *N*. It means that a great accuracy can be reached by just integrating the Riccati equation on the admittance and the evolution equation on the propagator—without the need of memory storage along the axis—with a small number of modes. Typically, at low frequency, taking three degrees of freedom achieves an accuracy better than 1% on the reflection and transmission coefficients (e.g. figure 7). The presented method can be applied to other classes of problems where the usual modes do not satisfy the right boundary condition; this is notably the case for the waveguide with Robin boundary condition.

## Funding statement

This work was supported by the Agence Nationale de la Recherche through grant ANR ProCoMedia, project ANR-10-INTB-0914.

## Acknowledgements

We acknowledge Eric Luneville and Christophe Hazard for fruitful discussions.

## Appendix A. Coupled mode equations

To obtain the coupled mode equations, equation (2.2) is projected on the (*N*+1) modes *φ*_{n}. For the first equation ∂_{x}*p*=*u*, the projection is That yields the first coupled mode equation
For the second equation ∂_{x}*u*=−(∂_{y2}+*k*^{2})*p*, the calculation is as follows. Integrating by part leads to
where the boundary conditions ∂_{y}*p*(*x*,0)=0 and ∂_{y}*p*(*x*,*h*)=*h*′*u*(*x*,*h*) have been used. Then for 0≤*n*≤*N*−1 (and −1≤*m*≤*N*−1),
By using the symmetry of (∂_{y}*φ*_{n},∂_{y}*φ*_{m}), it appears that for any −1≤*n*≤*N* and −1≤*m*≤*N*. Inserting this result into the projection of the second equation gives
This equation can be simplified owing to the normalization of the modes (*φ*_{n},*φ*_{m})=*δ*_{nm} that when differentiated with respect to *x* gives (*φ*_{n},∂_{x}*φ*_{m})+(∂_{x}*φ*_{n},*φ*_{m})+*h*′*φ*_{n}(*h*;*x*)*φ*_{m}(*h*;*x*)=0. Eventually, the second coupled mode equation is

## Appendix B. Expressions for the two-dimensional case

**(a) The boundary mode**

The function *χ* is chosen such that , where *a*_{N} is a normalization factor to be found and where
The projection of on the eigenmodes *φ*_{n} ( with ∂_{y}*φ*_{n}=0 at *y*=0 and *y*=*h* and where *ϵ*_{n}=2−*δ*_{n0}) is obtained by an integration by parts
With and (*γ*=*π*/2*h*), one gets . Thus, the explicit expression for is
A function orthogonal to *φ*_{n} (0≤*n*≤*N*) is with . The boundary mode *φ*_{−1} is normalized

**(b) Transverse wavenumber ∥∂ _{y}φ_{n}∥^{2}**

For the usual modes, 0≤*n*≤(*N*−1), we have .

To compute ∥∂_{y}*φ*_{−1}∥^{2} an integration by parts is used
so that . The explicit expression is

**(c) The coupling matrix A**

The elements of the coupling matrix are defined by
where *Q*_{nm} are constant, and with, for 0≤*n*≤(*N*−1), 0≤*m*≤(*N*−1)
B1In these expressions, *ϵ*_{n}≡2−*δ*_{n0} and
To derive these expressions, the following relations have been used. For diagonal elements, differentiating the normalization condition (∥*φ*_{n}∥^{2}=1) with respect to *x* gives . For non-diagonal elements, *n*≠*m*, successive integrations by parts have been used. For non-diagonal elements with *n*=−1 and 0≤*m*≤(*N*−1) (the row number −1), successive integrations by parts yield can be used. For instance, one finds an intermediate result
To obtain *A*_{n,−1}, it is sufficient to use the differentiation w.r.t. *x* of the mode orthogonality (*φ*_{n},*φ*_{−1})=0 to get (*φ*_{n},∂_{x}*φ*_{−1})+(*φ*_{−1},∂_{x}*φ*_{n})+*h*′*φ*_{n}(*h*;*x*)*φ*_{−1}(*h*;*x*)=0, which implies

## Appendix C. Expressions for the axisymmetric case

**(a) The boundary mode**

The transverse modes *φ*_{n}(*r*;*z*) (0≤*n*≤(*N*−1)) are defined by
where *β*_{n} is the (*n*+1)th zero of *J*_{1} (with five digits: *β*_{0}=0, *β*_{1}=3.8317, *β*_{2}=7.0156, *β*_{3}=10.1735…), *J*_{n} denotes the *n*th Bessel function and *f*_{n}=*J*_{0}(*γ*_{n}*r*) with *γ*_{n}=*β*_{n}/*h* are eigenfunctions of the transverse eigenproblem . The definition of the scalar product is . To build the boundary mode *φ*_{−1}, *χ* is chosen such as ∂_{r}*χ*(*r*=*h*;*z*)≠0 and in the following:
where *β* is the first zero of *J*_{0} (with five digits: *β*=2.4048) and where *a*_{N} is the normalization factor such that ∥*φ*_{−1}∥=1. The projection of on the eigenmodes *φ*_{n} is obtained as in the two-dimensional case by an integration by parts: , where *γ*_{n}=*β*_{n}/*h* and *γ*=*β*/*h*. The explicit expression for is thus
A function orthogonal to the *φ*_{n} (0≤*n*≤*N*) is with . Therefore, we choose to normalize the boundary mode

**(b) Transverse wavenumber ∥∂ _{r}φ_{n}∥^{2}**

For 0≤*n*≤(*N*−1),
To compute ∥∂_{r}*φ*_{−1}∥^{2}, integrations by parts give . With and , the explicit expression is

**(c) The coupling matrix A**

As in the two-dimensional case, the coupling matrix can be calculated using integration by parts and the elements of *A* are defined by
with, for 0≤*m*,*n*≤(*N*−1)
C1In these expressions,

- Received July 11, 2013.
- Accepted January 3, 2014.

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