The Wave-Matching Boundary Integral Equation — An energy approach to Galerkin BEM for acoustic wave propagation problems

In this paper, a new Boundary Integral Equation (BIE) is proposed for solution of the scalar Helmholtz equation. Applications include acoustic scattering problems, as occur in room acoustics and outdoor and underwater sound propagation. It draws together ideas from the study of time-harmonic and transient BIEs and spatial audio sensing and rendering, to produce an energy-inspired Galerkin BEM that is intended for use with oscillatory basis functions. Pivotal is the idea that waves at a boundary may be decomposed into incoming and outgoing components. When written in its admittance form, it can be thought of setting the Burton–Miller coupling parameter differently for each basis function based on its oscillation; this is a discrete form of the Dirichlet-to-Neumann map. It is also naturally expressed in a reflectance form, which can be solved by matrix inversion or by marching on in reflection order. Consideration of this leads to an orthogonality relation between the incoming and outgoing waves, which makes the scheme immune to interior cavity eigenmodes. Moreover, the scheme is seen to have two remarkable properties when solution is performed over an entire obstacle: (i) it has a condition number of 1 for all positive-real wavenumber k on any closed geometry when a specific choice of cylindrical basis functions are used; (ii) when modelling two domains separated by a barrier domain, the two problems are numerical uncoupled when plane wave basis functions are used — this is the case in reality but is not achieved by any other BIE representation that the authors are aware of. Normalisation and envelope functions, as would be required to build a Partition-of-Unity or Hybrid-Numerical-Asymptotic scheme, are introduced and the above properties are seen to become approximate. The modified scheme is applied successfully to a cylinder test case: accuracy of the solution is maintained and the BIE is still immune to interior cavity eigenmodes, gives similar conditioning to the Burton–Miller method and iterative solution is stable. It is seen that for this test case themajority of values in the interaction matrices are extremely small and may be set to zero without affecting ∗ Corresponding author. E-mail address: j.a.hargreaves@salford.ac.uk (J.A. Hargreaves). https://doi.org/10.1016/j.wavemoti.2018.07.003 0165-2125/© 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/ 4.0/). J.A. Hargreaves and Y.W. Lam / Wave Motion 87 (2019) 4–36 5 conditioning or accuracy, thus the linear system become sparse a property uncommon in BEM formulations. ©2018 TheAuthors. Published by Elsevier B.V. This is an open access article under the CCBY license (http://creativecommons.org/licenses/by/4.0/).

Th e w a v e-m a t c hi n g b o u n d a r y in t e g r al e q u a tio n -a n e n e r g y a Pl e a s e c h e c k t h e m a n u s c ri p t fo r a n y fu r t h e r c o py ri g h t r e s t ri c tio n s.

a b s t r a c t
In this paper, a new Boundary Integral Equation (BIE) is proposed for solution of the scalar Helmholtz equation. Applications include acoustic scattering problems, as occur in room acoustics and outdoor and underwater sound propagation. It draws together ideas from the study of time-harmonic and transient BIEs and spatial audio sensing and rendering, to produce an energy-inspired Galerkin BEM that is intended for use with oscillatory basis functions. Pivotal is the idea that waves at a boundary may be decomposed into incoming and outgoing components. When written in its admittance form, it can be thought of setting the Burton-Miller coupling parameter differently for each basis function based on its oscillation; this is a discrete form of the Dirichlet-to-Neumann map. It is also naturally expressed in a reflectance form, which can be solved by matrix inversion or by marching on in reflection order. Consideration of this leads to an orthogonality relation between the incoming and outgoing waves, which makes the scheme immune to interior cavity eigenmodes. Moreover, the scheme is seen to have two remarkable properties when solution is performed over an entire obstacle: (i) it has a condition number of 1 for all positive-real wavenumber k on any closed geometry when a specific choice of cylindrical basis functions are used; (ii) when modelling two domains separated by a barrier domain, the two problems are numerical uncoupled when plane wave basis functions are used -this is the case in reality but is not achieved by any other BIE representation that the authors are aware of. Normalisation and envelope functions, as would be required to build a Partition-of-Unity or Hybrid-Numerical-Asymptotic scheme, are introduced and the above properties are seen to become approximate. The modified scheme is applied successfully to a cylinder test case: accuracy of the solution is maintained and the BIE is still immune to interior cavity eigenmodes, gives similar conditioning to the Burton-Miller method and iterative solution is stable. It is seen that for this test case the majority of values in the interaction matrices are extremely small and may be set to zero without affecting

Introduction
Boundary Integral Equations (BIEs) can be an effective technique for solving wave scattering problems in homogeneous media. The partial differential equations that hold in the medium may be reformulated as an integral equation that holds on the boundary of the obstacle, thereby reducing the dimensionality of the problem and allowing unbounded homogeneous media to be modelled in a straightforward way. In this paper we will predominantly consider boundary integrals that represent acoustic wave phenomena, in which case the unknown field quantity is scalar.
The Boundary Element Method (BEM), being the most common manifestation of BIEs as a computational method, represents the boundary quantities as weighted sums of piecewise polynomial functions defined on 'elements'. These boundary quantities would be pressure and the surface-normal component of particle velocity, or equivalently pressure gradient, for an acoustics problem. Despite the advantages of BEM -more straightforward meshing, a reduction in the number of degrees of freedom relative to methods that discretise the quantities in the volume and ease of modelling unbounded problems -its use in acoustic engineering has mostly been in a small set of specific cases. These include, for example: prediction and optimisation of the sound scattering properties of acoustic objects and treatments [1,2] and noise barriers [3,4]; estimation of head-related transfer functions from head geometry [5,6].
Despite these advantages and the success of BEM in these areas, it is undeniable that almost all acoustic practitioners, particularly in room acoustics, select geometrical acoustic modelling techniques [7], such as ray-tracing, beam-tracing and image source, over approaches such as BEM that solve the Helmholtz or scalar wave equations directly. Partly this is likely due to accumulated expertise in applying these methods and the availability of material data, though the latter is almost always estimated and is widely acknowledged to be a source of inaccuracy [8]. The more significant reason is the computational cost of algorithms such as BEM, in particular with respect to its scaling with wavenumber k. This occurs because the meshes must be refined to keep the elements small with respect to wavelength; for a D dimensional problem, with D − 1 dimensional boundary, the number of elements, and therefore the number of degrees of freedom, must grow with O ( k D−1 ) to maintain accuracy as k increases. BEM produces full interaction matrices, linking every basis function to every other basis function, leading to computational cost and storage requirements that scale O ( k 2(D−1) ) . Matrix compression techniques such as fast-multipole [9] and adaptive-cross-approximation [10] are capable of significantly reducing these requirements, but the scaling of computation cost and storage with frequency is still sufficiently unfavourable so as to preclude full audiblebandwidth simulation for most realistic sized room acoustics problems of interest.

BEM with oscillatory basis functions
An alternative strategy is to attempt to build the oscillation expected in the solution into the approximation space; the piecewise-polynomial basis functions are replaced with a smaller, ideally k-independent, number of appropriately chosen oscillatory basis functions, whose support may span many wavelengths. Typically, these are then multiplied with a standard piecewise-polynomial approximation space -these will be referred to as 'envelope' functions in what follows -but these will be defined on a coarse mesh that is not refined with frequency.
Such algorithms may be broadly grouped into two categories: Hybrid Numerical Asymptotic BEM (HNA-BEM) [11] and Partition-of-Unity BEM (POU-BEM) [12]. HNA-BEM schemes attempt to pick an optimal set of oscillatory basis functions a priori, leading to a reduction of the scaling of the number of degrees of freedom required to maintain accuracy with k; O (log k) [13] or even k-independent [14] have been reported. The disadvantage with this approach is that schemes are only applicable to certain classes of obstacle, and they typically require detailed analysis of the problem in question. POU-BEM on the other hand, aims to construct a general-purpose set of oscillatory basis functions that produce a scheme requiring fewer degrees of freedom per wavelength than conventional piecewise-polynomial scheme. This may be significant, e.g. reduction to 2.5-3 degrees of freedom per wavelength for 'engineering accuracy' from the more usual 8-10, but the scaling with k is unchanged. The contributions included in this paper are intended that they might be useful for either of these classes of algorithms, but for simplicity we will construct an approximation space that is more consistent with POU-BEM principles.
When using oscillatory basis functions, a factor that offsets the gains described above is the cost scaling of the numerical integration. This is typically invariant of k in standard polynomial approximation spaces, since the elements get smaller as the wavenumber of the kernel increases. In both POU-BEM and HNA-BEM, the elements are large and invariant of wavelength, and the basis functions are oscillatory too, meaning integration cost for a single interaction matrix coefficient will scale at least O ( k (D−1) ) for a collocation scheme or O ( k 2(D−1) ) for a Galerkin scheme. For such schemes to achieve their potential, there is therefore a need for algorithms that can efficiently evaluate the types of oscillatory integrals that arise. There has been significant progress towards achieving this in two-dimensions using Filon type methods (see e.g. the review in section 4 of [11]). In three dimensions, methods for example using variations of the method of steepest descent [15], stationary phase [14] and Stokes' theorem [16] have been proposed.
An interesting aspect of BEM with oscillatory basis function is the way that many of the oscillatory functions, and the waves they scatter [16], begin to resemble terms from a geometric phased beam-tracing algorithm as k is increased. Plane waves are a common choice of oscillatory function leading to the idea that they have 'propagation directions'; this is a concept that many authors use informally in their prose even if they do not formally acknowledge it in the mathematical analysis, and geometric phased beam-tracing has even been used to find optimal wave directions and support [17]. These plane waves are windowed spatially by the standard piecewise-polynomial approximation space, so naturally are 'phasedbeams' as they depart the boundary. At finite k diffraction will occur due to finite aperture effects and the scattered wave will spread out, but as k → ∞ they become beams with geometrically-defined edges. These observations suggest that BEM with oscillatory basis function might, in addition to the reasons above, be worthy of further study because it provides a potential pathway to unification between standard BEM techniques, which are accurate, efficient and well established for small ka (where a is a length typical of the problem size), and geometric acoustic algorithms, which are accurate, efficient and well established for large ka. Such a unification would be extremely impactful in acoustic modelling, where the full audible bandwidth covers many octave and wavelength changes from being large with respect to typical obstacle size to being very small, hence neither class of algorithm can presently meet all requirements. A parallel development is the application of boundary integral methods to model acoustic energy flow at high frequencies [18][19][20], again emphasising that BIEs are a framework that can be applied over a wide range of ka.
One recurring feature of geometric algorithms is that they are solved by marching on in reflection order. It will be seen that the BIE proposed herein also lends itself to solution in this manner, a property that has the potential to be interesting since it could in principal allow optimal wave directions (in the sense of HNA-BEM) to be found adaptively for each reflection order.
Another key characteristic of geometric algorithms is that they only consider the local boundary orientation, curvature and conditions when computing reflections. This is rather like considering the obstacle to be a collection of 'canonical' scattering problems, each of which encapsulate a certain scattering behaviour, that have been spatially windowed and joined together as illustrated in Fig. 1. The scattering predicted by a geometric algorithm would be the solution to each of those canonical problems, spatially windowed by the extent of the boundary that the problem matched. Typically, most geometric algorithms are only concerned with reflections from planar boundary sections (centre), though it is also possible to handle low-order reflections from simply curved boundary sections, e.g. a circular cylinder (right), and algorithms that compute edge diffraction by truncating the infinite wedge canonical problem (left) also exist [21,22]. Some HNA-BEM algorithms can be considered to follow a similar principal in that they use different oscillatory functions on different boundary sections [13,14,23].
This point is emphasised because it justifies the consideration in this paper of two separate, simple canonical scattering problems using oscillatory basis functions; the infinite plane and the circular cylinder. Clearly more interesting geometries would require the joining of segments 'cut' from these canonical problems, but that is identified as future work; this preliminary study is limited to considering these in isolation.

Non-uniqueness
It is well known that a BEM model of a scattering problem simultaneously solves two coupled problems: the intended exterior problem and a complementary interior problem. When the frequency hits an eigenfrequency of the complementary interior problem, the solution of the boundary integral equation becomes non-unique and conditioning of the resulting numerical system breaks down. This issue has been studied extensively, and two dominant mitigation techniques exist; the Combined Helmholtz Integral Equation Formulation (CHIEF) method [24] and that of Burton and Miller [25]. In the former, additional field points from the interior domain are added to the linear system of equations, making it overdetermined and excluding the internal field from the solution, so long as the points chosen do not happen to lie on nodes of the internal eigenmode. Burton and Miller showed that using a linear combination of the BIE and its surface-normal derivative gives a unique solution for all k, so long as the 'coupling parameter' that defines the relative weight of the two constituent BIEs is appropriately chosen. Choosing an optimal value of this parameter has since become an area of study in its own right. Generally the objective is to minimise the condition number of the resulting linear system, with earlier work typically concentrating on idealised obstacles such as spheres and cylinders [26][27][28] and latter work extending results to more complicated classes of obstacle [29,30]. Terai [31] commented that ''it is clearly from physical considerations'' that the sign of coupling parameter should be inverted depending on whether e +iωt or e −iωt time dependence is assumed, a detail that follows naturally from the time domain implementation discussed below, and Marburg showed that this can affect the solution in some cases [32]. It has also been shown that different values of the coupling parameter can be chosen on different sections of the obstacle so long as its imaginary part is one-signed [33]; the ideas presented in this paper can equivalently be thought of as applying different values for the coupling parameter to different basis functions.
Another research area where this effectively is done is in the group of algorithms that apply a Dirichlet-to-Neumann map (or vice-versa) to regularise the BIE and ensure a unique solution [10,[34][35][36][37]. These are typically stated as being motivated by the work on On-Surface Radiation Conditions (OSRCs) [38], though appear to most commonly be applied to generate matrix pre-conditioners [10,39]. The connection between these approaches and the Burton-Miller formulation is well-established, and it will be seen that the BIE proposed herein has strong similarities to some of them too. Despite this, the ideas presented herein remain novel in that: (i) the BIE, while similar to some that have previously been published, is not identical to any and is derived differently; (ii) the scheme considers admittance boundary conditions instead of just Neumann or Dirichlet; (iii) the scheme uses oscillatory basis functions and states the Dirichlet-to-Neumann map separately for each of them.
The CHIEF and Burton-Miller approaches were compared for POU-BEM in [40]. Here the authors concluded that CHIEF outperformed Burton-Miller, though implementation details such as regularising coordinate transforms appear to have had a significant effect. It is worth pointing out here that regularisation of the hypersingular operator, which occurs in the Burton-Miller formulation, is rather challenging for collocation schemes, as was used in [40]. In contrast, Galerkin schemes allow the hypersingular operator to be rewritten in a form that only contains a weak singularity [41,42].
Issues due to the complementary interior problem also affect time domain BEM, though here they manifest as solver instability. An approach equivalent to the Burton-Miller method has been applied for both scalar acoustic problems [43] and vector electromagnetic problems [44] and was shown to be effective, though instability may still arise due to discretisation error and is more likely when lightly-damped external features are present [45]. Chappell et al. [46] analysed these schemes using a similar technique to that which was used to analyse frequency domain formulations, but the analysis in [43] focuses instead on the trapping of energy in the fictitious internal cavity. Such ideas were later taken further by deriving the Galerkin BEM formulation directly from an energy identity [47,48]; a similar approach will be adopted in this paper, but for frequency domain problems with oscillatory basis functions.

Application of BIEs in spatial audio rendering
The same BIEs that appear in acoustic scattering problems also appear in the study of sound-field rendering [49,50]. Here the Kirchhoff-Helmholtz boundary integral equation (KHBIE), which is often interpreted as boundary source densities, would represent a physical loudspeaker array [51]. Often a second boundary integral is set up to evaluate the accuracy of the rendering by the source array; this has the physical interpretation of being an array of microphones, but mathematically it resembles the inner-product 'testing' integral of a Galerkin BEM scheme. Equivalently a collocation scheme would represent a collection of equally-weighted, spatially discrete microphones.
It is known that the same non-uniqueness issue that affects BEM also affects microphone arrays [49,50,52]. The configuration in [50] is particularly interesting since the 'testing' integral is physically separated from the boundary on which the source densities are defined, making it immediately clear that non-uniqueness is associated with the 'testing' integral rather than the KHBIE. This fact is also known within the scattering community, manifesting in the standard BIE and its surface-normal derivative version being non-unique at pressure-release and rigid cavity eigenfrequencies respectively [53], but it is not so immediately obvious due to the two integrals being defined on the same boundary. It is also known that the non-uniqueness problem in microphone arrays may be overcome by using directional microphones e.g. Cardioid microphones as considered in [52]. These measure a linear combination of pressure and surface-normal pressure gradient, hence are in fact implementing a physical manifestation of the Burton-Miller method.
In the context of loudspeaker and microphone arrays, the presence of oscillatory basis functions in the integrals introduces the possibility of beam-steering of the array's radiation or sensitivity pattern. This also introduces the possibility that the relative weighting of the sensing of pressure and its surface-normal derivative might be optimised for a particular wave direction. This is done in [54] to eradicate 'front-back confusion' in planar arrays, and in [54][55][56] with spherical arrays to distinguish between waves arriving from outside or emanating from inside the array. Green's theorem is used to extend this to arbitrary shaped arrays in [57]; a similar formulation to this will be derived in two dimensions herein.

Overview of this paper
This paper proposes a new boundary integral equation for the scalar Helmholtz equation, motivated by the state of the art in time-domain BEM and microphone array design. This will be developed for two canonical scattering problems -an infinite plane and a cylinder -with a view that these may somehow later be stitched together to allow modelling of more realistic problem geometries. A particular aim is to ensure that fictitious interior cavity resonances are avoided; this is identified as a prerequisite for avoiding non-uniqueness in the frequency domain and achieving stability if solving the system iteratively by marching on in time or reflection order. Parallels to sound-field rendering and sensing will be made where possible to give physical context to the equations presented. Section 2 will present the standard boundary integral operators, discretisation with oscillatory basis functions and solution by Galerkin BEM. This material is well-known, but is included because the new scheme will build on it, and because the schemes will be used as a reference in Section 5.
Section 3 presents the new BIE. First the integral equation itself is derived from an energy identity, then the key concept of incoming and outgoing waves is articulated, moving later to specific forms for planar and cylindrical geometries. Section 4 uses this numerical machinery to solve the two canonical problems, without the inclusion of envelope functions. Some remarkable results are derived regarding wave orthogonality and linear system conditioning, and the extent to which these properties apply for broader classes of boundary geometry is also considered.
Section 5 introduces envelope functions, thereby testing the scheme in a way that is more representative of how it is expected to ultimately be used.
Finally, Section 6 draws conclusions and makes suggestions for future research priorities.

Standard Boundary Integral Formulations
In this section some standard boundary integral formulations will be presented, along with their numerical implementation as Galerkin BEM. In all that follows, we will consider scalar problems that can be modelled by the Helmholtz equation in two dimensions: Here wavenumber k is assumed to be positive and real and x is a point in 2D Cartesian space. is the scalar potential, which in acoustic applications would typically be the pressure of an acoustic wave; it will hereafter be referred to as pressure.

Standard BIE formulations for scattering problems
This study will mostly be restricted to exterior scattering problems; the configuration is depicted in Fig. 2. The obstacle boundary is assumed to be a Lipschitz contour, and it separates the space into interior and exterior domains − and + respectively, the latter of which contains the acoustic medium of interest. Total pressure (x) = ϒ (x) + (x), where ϒ is a wave incident on the scattering obstacle (defined in the absence of the obstacle) and is the wave that is scattered in response to this. The scattered wave satisfies Sommerfeld's radiation condition in the far-field of the obstacle. Under these conditions, (x) may be found from on the obstacle boundary by the KHBIE: where: is the fundamental solution (free-space Green's function) of the Helmholtz equation in 2D, with H + 0 being a zeroth order Hankel function of the first (outgoing) kind. The notation ∂/∂n y is shorthand forn y · ∇ y , wheren y is a unit vector pointing normal to and into + at point y ∈ , and subscript y on ∇ means ''with respect to y''. dl y in (2) denotes a line integral where y is the integration variable. It is common to write boundary integral equations in terms of two standard operators: the single-layer potential S and the double-layer potential D. These are defined as: Here φ is some scalar field defined on the boundary. It follows from (2) that: Often it will also be necessary to evaluate ∂ /∂n x (x), the surface-normal derivative of for x ∈ , in which case two further standard operators are introduced; the adjoint double-layer potential A and the hypersingular operator H: Taking the surface-normal derivative of (6) gives: It will be assumed that the obstacle has a locally-reacting admittance boundary condition, so (y) and ∂ /∂n y (y) are related by: Here Y is the value of locally-reacting specific admittance, which for simplicity is considered to be uniform over the obstacle.
To solve for and in a scattering problem, x must lie either limiting close to in + or somewhere in − ; in the latter case, that will be used herein, this means that (x) = 0. In practice x will still be chosen to be limitingly close to , lying on a surface − ⊂ − that is conformal to . The only common exception to this is CHIEF points [24], where an additional set of x ∈ − are chosen to make the resulting system of linear equations overdetermined. When x ∈ − , the traces of the operators D and A mean that (6) and (9) become: Here the admittance boundary condition of (10) has also been substituted and the identity operator I {φ} (x) = φ (x) has been introduced.
Since (x) = 0 for all x ∈ − implies that ∂ /∂n x (x) = 0 for x ∈ − , it is also possible to solve the surface-normal derivative form of (13): This BIE admits interior fields that satisfy ∂ /∂n x (x) = 0 for x ∈ − , sometimes called hard cavity modes. The Burton-Miller formulation [25] uses a linear combination of (13) and (14), thereby solving ∂ /∂n x (x) − α (x) = 0 for x ∈ − . Here α is the coupling parameter discussed in Section 1.2. This produces the BIE: Burton and Miller showed that (15) has a unique solution for all real values of k, so long as the imaginary part of α is non-zero. We follow common practice [28,32] and choose α = ik.

Interpretation of the Burton-Miller formulation in terms of microphone directivity
Microphone directivity is usually assessed in a far-field sense by evaluating the response of the sensor to plane waves. The interpretation of the fundamental solution in (3) as an outgoing wave implies e −iωt time dependence; with this a plane wave defined (x) = e ikd·x propagates in the direction of the unit vectord as time progresses. As stated above, the Burton-Miller formulation attempts to satisfyn x · ∇ x (x) − α (x) = 0. If (x) is a plane wave as defined above, this becomes The termn x ·d − 1 is recognised as a Cardioid directivity with its maximum sensitivity orientated in the directionn x and its 'null' direction of zero sensitivity being −n x ; this is the orientation suggested to overcome the shortcomings of spherical microphone arrays in [52]. It can be concluded from this that the Burton-Miller formulation will 'react' more strongly to waves arriving from outside and less strongly to waves arriving from inside, with plane waves propagating in exactly the directionn x automatically satisfying This observation is inextricably intertwined with the energy trapping argument in section 1D of [43] and the fact that (15) does not support time-invariant trapped interior eigenmodes.

Discretisation and solution using the Galerkin method
on is approximated by a function˜ ∈ R N×Q defined as a weighted sum of basis functions b n,q : Here a double summation over n and q has been used because it will simplify the introduction of oscillatory basis functions. The ⌊· · ·⌋ brackets denote indicate rounding down to the nearest integer and Q must be odd. ϕ ∈ C N×Q is the set of discretisation coefficients and R N×Q ⊂ L 2 ( ) is the N × Q -dimensional approximation space spanned by˜ . The basis functions b n,q (y) will be taken to be a product of an envelope function e n (y) and an oscillatory function o n,q (y).
q will be used to index wave direction, hence the envelope functions are independent of it. The oscillatory function includes the index n to allow for the possibility that oscillatory functions might differ for different envelope functions, as would be required in the 'windowed canonical problem' concept illustrated in Fig. 1. The envelope functions e n (y) are defined only for y ∈ , and would typically be standard piecewise polynomial shape functions defined on a coarse mesh of elements in most existing HNA-BEM or POU-BEM schemes. We assume that the oscillatory functions o n,q (y) are solutions of the Helmholtz equation, so should also be defined for y / ∈ . (16) can represent a non-oscillatory scheme if Q = 1 and o n,0 (y) = 1 are chosen for all n.
Solution by the Galerkin method of (13) involves finding˜ ∈ R N×Q such that for all T ∈ R N×Q : Equivalently for the Burton-Miller formulation in (15): Here a straight bar over a quantity indicates conjugation, and T is some testing function that is also a member of the approximation space R N×Q . A common choice is to take the set of testing functions to equal the set of basis functions, hence the conjugation of the b m,p terms that appear below. Note that the integration domain − has been replaced by since the explicit inclusion of the traces of D and A (being the ± 1 2 I terms) mean that x ∈ will be evaluated as x ∈ − .
Eqs. (17) and (18) are respectively solved numerically by the following matrix equations: [ The entries of the matrices introduced here are defined as: Here it is assumed that the pairs of subscripts m, p and n, q are each collapsed into a single index so that S, D, A, H and M are matrices and υ and υ ′ are vectors. The matrix M corresponds to the operator I but has been named M to avoid confusion with the standard identity matrix I and because its definition matches the well-known 'mass matrix' that is found in Finite Element Method.
Evaluation of these matrices is complicated by the fact that G (x, y) → ∞ as |x − y| → 0. For a smooth boundary, as is considered herein, (22) and (23) can be shown to be non-singular, so can be computed using standard quadrature techniques. Eq. (24) includes a stronger singularity, but in practice the matrix H will be evaluated using the following statement [41,42], which only contains a weak singularity: This leaves only weak singularities in (21) and (28); these were regularised using the Sato transform [58] of order 5, as has been shown to be optimal for weakly-singular oscillatory kernels [59].
This concludes the presentation of existing formulations. The following section will build-upon this and present the new BIE that is the key contribution of this paper.

The proposed 'wave-matching' BIE
Here we propose a new BIE that is, like the time domain BIEs proposed in [47,48], derived from energy principles. The approach here will however differ because we are considering a frequency domain formulation. In particular, the formulations in [47] and [48] begin by considering the time-derivative of the total system energy, but this is necessarily zero in a time-harmonic scheme, including for the cavity eigenmodes the aim is to prevent. Ergin et al. [43] showed that their time-domain version of the Burton-Miller method allowed energy to escape from the cavity. We will instead begin by requiring that energy cannot enter the cavity. This is a natural extension of the conditions that (x) = 0 and ∂ /∂n x (x) = 0 for x ∈ − , as discussed in Section 2.1. It will be seen in Section 4 that this produces a BIE that optimally (i.e. entirely) both excludes energy from entering the cavity and entirely allows it to escape, for both finite obstacle and half-space problems. We focus, as in the rest of the paper, on the case where k is purely real, but remark that this same approach should also be directly applicable to lossy media where I (k) ̸ = 0.
Regarding the name proposed for the BIE, Aimi and her co-workers [48,60] have called their BIE various names usually involving the terms ''energy'' or ''energetic''. The BIE proposed herein could be described in similar language, but it is felt that this would cause confusion with the energy BEM described in [18][19][20]. Hence the name 'Wave-Matching', which reflects some properties of the BIE that are demonstrated in Section 4 and was used for the same ideas in [61], is proposed instead.
Time averaged acoustic intensity I is given by [62]: Here the first statement is the standard definition, with U being acoustic particle velocity. The second statement has the real operator expanded as half the sum of the argument and its conjugate. In the third step Euler's equation −ρ 0 du/dt = ∇ϕ has been substituted, which for e −iωt time dependence becomes iωρ 0 U = ∇ . Here ρ 0 is the equilibrium density of air. Acoustic intensity is a measure of power flux density in Watts/metre 2 . It and acoustic energy density E (y) are related by the well-known energy-flux relation dE/dt = −∇ · I . Applying the divergence theorem over the interior volume − and then noting that dE/dt = 0 in a time-harmonic problem gives: This equation has the physical interpretation that, because the KHBIE models − as being a lossless medium (like + ) that is source-free, energy is not created or destroyed in − and any energy that enters through the boundary will also leave again. It does not however preclude trapped cavity eigenmodes, for which there is no energy flow across .
The key ingredient that allows us to distinguish between energy entering and leaving − is the idea that (x) may be decomposed into a summation of two waves, − (x) and + (x), that are incoming (−) and outgoing (+) at the boundary . Here + and -have been chosen to reflect the sign of the dot product between their nominal direction vectors and the (outward pointing) boundary normal vectorn y . This idea has previously been used with piece-wise constant interpolation on planar elements for time-domain BEM [63], the primary reason being to ensure causal boundary conditions, and is proposed as a more generally valid decomposition in [64]. Here we do not attempt to consider the mathematical details of how this might apply to different classes of boundaries, but instead in Section 3.1 demonstrate that the decomposition is straightforward for the two canonical problems we are considering; a cylinder and an infinite plane.
We assume that − (x) and + (x) can both be approximated by (16) but with differing sets of coefficients ϕ − and ϕ + and differing sets of basis functions b − n,q and b + n,q . For the latter, it will be assumed that the envelope functions are common to both, but that the oscillatory functions differ; this enables construction of a reflectance-based scheme. It will be assumed that − (x) and + (x) are represented by a weighted sum of o − n,q (x) and o + n,q (x), but that these are then windowed by the envelope functions e n (x). We also assume that ∂ − /∂n x (x) and ∂ + /∂n x (x) are represented by ∂o + n,q /∂n x (x) and ∂o − n,q /∂n x (x) windowed in the same way; that is the basis functions possess a surface-normal derivatives ∂b ± n,q /∂n x (x) = e n (x) ∂o ± n,q /∂n x (x). This is a straightforward and valid choice but may not be optimal is terms of achieving only inward or outward propagation; this is discussed more in Section 5. It is observed that (30) already takes the form of a pair of inner-product integrals. We therefore use it to set up a Galerkin BEM scheme, where the testing functions will be substituted in place of the conjugated terms involving ; these are drawn solely from the set of incoming waves that can be represented by (16) , where˜ is the result of substituting the approximation˜ to into (11), giving: This is the 'Wave-Matching' BIE for scattering problems. It is enforced for all T ∈ R − N×Q .

Incoming and outgoing waves
In this section the decomposition into incoming and outgoing waves for two boundary geometries is presented: an infinite planar boundary and a cylindrical obstacle submerged in a medium. In both cases the waves will be parameterised by their oscillation along the boundary, expressed as a tangential wavenumber k t . k t will be real for problems with real k, and may take both positive and negative values. The values it takes will be defined from q differently for different boundary geometries; this is detailed in the following sections. Its usage draws an interesting parallel with the energy BEM scheme in [19], which parameterises wave direction in the same way. The main difference is that we will allow evanescent waves for which |k t | > k, whereas the scheme in [19] uses a geometric propagation operator, so cannot include evanescent waves and must restrict |k t | ≤ k.

Infinite planar boundaries
Consider a planar boundary, so a straight line in two-dimensions. We continue to assume that the surface-normal vector n, which is position invariant for a plane, points into the domain + that contains the acoustic medium.
A pressure field (x) = e iktt·[x−x 0 ] , wheret is the tangential unit vector that is perpendicular ton and x 0 ∈ , could arise on the boundary from either (or some linear combination) of the following two waves: Assuming e −iωt time dependence, + can be seen to be propagating into + , so is outgoing with respect to the boundary, whereas − is propagating in a direction incoming from + .
These are the outgoing and incoming wave pair for this value of k t . It can be seen that ∂ + /∂n x (x) = ik n + (x) and ∂ − /∂n x (x) = −ik n − (x) for all x ∈ . We will denote these ratios α ± (k t ) = ±ik n . We remark that multiplication by α ± (k t ) implements a Dirichlet-to-Neumann map, since ∂ ± /∂n x (x) = α ± (k t ) ± (x), and note that expressions like that for k n appear in the OSRC literature e.g. following (10) in [36]. In the case where k t > k, k n becomes imaginary and the wave becomes evanescent, obeying a non-oscillatory exponential taper in the direction ±n. The choice of whether k n follows the positive or negative imaginary branch of the square root is in some senses arbitrary, but we choose the positive branch because this makes physical sense. This means an outgoing evanescent wave, arising from some property of the obstacle boundary that varies with k t > k, will decay away from the obstacle as real evanescent waves do. Similarly, an incoming evanescent wave will be decaying away from whatever body radiated it.

Cylindrical boundaries
Consider a cylindrical boundary, so a closed circular arc in two-dimensions, of radius a. We define a polar coordinate system (r, θ) that has its origin at the centre of the cylinder. The domain outside the cylinder is deemed to contain the acoustic medium and is termed + and the interior is − . The surface-normal vectorn points into + and equalsr, the radial unit vector of the coordinate system.
A pressure field (a, θ) = e iaktθ could arise on the boundary from either (or some linear combination) of the following two waves: Here H + akt and H − akt are Hankel functions of order ak t that are outgoing and incoming respectively; assuming e −iωt time dependence these are respectively of the first and second kind.
as a Dirichlet-to-Neumann map.

Comparison of planar and cylindrical waves
As the radius of a cylinder increases the curvature of its boundary reduces until it becomes planar as a → ∞. We therefore expect that for the large a the incoming and outgoing cylindrical waves described in Section 3.1.2 should begin to match the planar versions in Section 3.1.1. We also note that the definitions for Dirichlet-to-Neumann maps in the OSRC literature factor in boundary curvature, so dependence of α ± (k t ) on a is expected.
It is clear from their formulae that all the waves match in tangential oscillation, as this is governed by the tangential wavenumber k t in both cases. It is however also clear that they will not match in magnitude, since the plane waves in Section 3.1.1 are unitary and the cylindrical waves in Section 3.1.2 have a magnitude that depends on radius; this can be overcome by simply scaling the latter, as will be done in Section 4.2.6. What is more important is the ratio of the surface-normal derivatives, as captured in the quantities α ± (k t ).
These are plotted for the outgoing case in Fig. 3. The heavy-black line termed a = ∞ is the planar case. It can clearly be seen that the values for cylinder tend to this curve as ka is increased. At the ends of the k t scale when ka is large, α + (k t ) → ik as k t → 0 and α + (k t ) → −k t as k t becomes large. That these trends apply in ka rather than a simply reflects that it is the radius of curvature compared to the wavelength that is important.
Regarding α − (k t ), a peculiarity is that α − (k t ) = −α + (k t ) for the planar whereas α − (k t ) = α + (k t ) for the cylindrical case. It will be seen in Section 4.2.6 that it will be necessary to take α − (k t ) = −α + (k t ) in the cylindrical case too when using a reflectance boundary condition and envelope functions.
An important point about both these incoming and outgoing wave decompositions is that they are not limited to the boundary geometries that suggested them. They are however much easier to use with this geometry that with arbitraryshaped boundaries. This is because the waves are 'separable' on such boundaries -that is, they can be written as the product of a boundary-normal term and a boundary-tangential term. This makes the α coefficients invariant over the boundary, simplifying the definition of reflectance boundary conditions (below) and meaning the Wave-Matching scheme matrices can be readily calculated from the standard Galerkin BEM scheme matrices; these latter relations are presented in Section 4.2.6 and subsequently used for the schemes in Section 5.

Formulation using a reflectance boundary condition
The decomposition into incoming and outgoing waves suggests the use of a reflectance boundary condition. In this the discretisation weights ϕ + and ϕ − , for − and + respectively, are related by: Here R is a reflectance matrix, the entries in which will be discussed further below. This form makes it easy to check that the boundary condition is causal and passive (i.e. does not add energy).
Inserting the approximate statements for = − + + into (31) gives the following matrix equation, the terms of which will be defined below:W + ϕ + +W − ϕ − = −υ. This cannot be solved until the boundary condition is substituted because it contains too many unknowns. Doing so gives: This may be solved to find ϕ − , then ϕ + may be found from ϕ − by (32). Note that a 'breve' over a matrix or vector will be used to indicate that it is the Wave-Matching scheme version.
In (31), the testing functions T (x) were drawn from the set of incoming basis functions b − m,p (x). These will be one of the incoming waves discussed in Sections 3.1.1 and 3.1.2 multiplied by an envelope function; exact definitions will be given in the sections that follow. The right-hand side vectorυ of (33) is defined as: The matricesW ± are most easily expressed as the following sum of four terms: Here the definitions for ∂ /∂n x (x) and (x) in (6) and (9) have been applied to the approximate form of ± following (16). Values of these for specific scenarios will be given in Sections 4 and 5. The reflection coefficient R at a point y may be found by applying the admittance boundary condition (10) to a pair of basis function: When the obstacle has a geometry on which the incoming and outgoing waves are 'separable' (see note in previous section), this becomes position invariant.

Iterative solution by marching on in reflection order
The reflectance boundary condition leads naturally to a process of iterative solution by marching on in reflection order. In this, the incoming and outgoing waves are expressed a sum of reflections ± (x) = ∑ ∞ r=0 ± r (x), where r indexes reflection order and all the terms ± r (x) are approximated by (16) with discretisation weight vectors ϕ ± r . (33) is solved iteratively by first solving for ϕ − r at each reflection order using: ϕ + r is then found by applying the reflectance boundary condition (32) applied to ϕ − r . This scheme converges if the magnitude of all eigenvalues of the update matrix inv (W − )W + R are smaller than one. This is examined numerically for a cylinder in Section 5.2.3.

Formulation using an admittance boundary condition
It is not essential to perform discretisation using the incoming and outgoing waves and reflectance in a Wave-Matching scheme. In this formulation the incoming wave is still used for 'testing' purposes, but the discretisation of surface quantities is performed using the basis functions b m,p from Section 2.2, with ∂ /∂n y (y) related to (y) by the admittance boundary condition (10). The testing functions T (x) in (31) will again be the set of incoming basis functions b − m,p (x). We now solve: The right-hand side vectorυ is unchanged from (34). The four matrices, now lacking ± exponents since they are not associated with incoming or outgoing waves, are defined: This concludes the definition of the components of the scheme. It should be noted that these have all been stated without placing restrictions on the geometry of , though limitations will in practice apply because methods for construction of appropriate basis functions are only so far proposed for some specific classes of boundary geometry. In Sections 4 and 5 they will be applied to model scattering from the cylinder and plane canonical problems, without and with windowing by envelope functions respectively. This will allow various simplifications to be made and features observed.

Scheme for entire canonical problems
In this section we will investigate the properties of the new Wave-Matching BIE on two entire canonical scattering problems; an infinite plane in Section 4.1 and a cylinder in 4.2. By entire we mean that we will not use any envelope functions.
Or more precisely, and following the discretisation scheme in (16), we will take N = 1 and define e 0 (x) = 1 for all x ∈ .
Since n always takes the value 0, it will be dropped throughout this section.

Planar boundaries
Here we consider the infinite planar boundary defined in Section 3.1.1. For a given surface-tangential wavenumber k t , the incoming and outgoing basis functions used to discretise + and − are: Here k n = √ k 2 − k 2 t is the surface-normal component of wavenumber and x 0 ∈ is some phase reference point. Wave direction has been indexed by the real-valued tangential wavenumber k t , rather than integer-valued q as used elsewhere in the paper; this is because an inner product on an infinite boundary will give rise to a Fourier transform whereas on finite boundaries it will produce Fourier series. The ratios α ± (k t ) are defined: These can be used to find the Wave-Matching scheme matrices in Sections 3.2 and 3.3 from the standard Galerkin BEM matrices. We would now regard them as entries from the discrete Dirichlet-to-Neumann map, since they apply to individual basis functions within the approximation space rather than the total solution. They can also be used to simplify (40), again using k t in place of q, to find the reflection coefficient that would occur if the plane had admittance Y : For non-evanescent plane waves k n = k cos θ, where θ is the angle the wave direction makes with the boundary normal vector. (49) then simplifies to the plane wave pressure reflection coefficient widely used in acoustics: For rigid boundaries R = 1 for all k t .

Radiation of the incoming and outgoing waves by the boundary operators
Here we consider the KHBIE (6) in its interpretation as a sound-field rendering operator, as is physically implemented in Wave Field Synthesis [51]. In particular it is well known that for radiation problems, for which ϒ (x) = 0, (x) = (x) if x ∈ + and (x) = 0 if x ∈ − . Care must however be taken in scattering problems, since + is bounded by both and another distant boundary over which the Sommerfeld radiation condition was applied; because of this the behaviour is not always as expected. Radiation of + is straightforward; it can be understood directly in a radiation sense since energy is leaving . We find Radiation of − is more easily understood in a scattering sense by assuming ϒ (x) = − (x). In this case ϒ (x) can physically be thought of as the 'shadow' wave, that cancels out whatever wave is incident within − . Applying these relations to the basis functions in (47) gives: Note that if is planar then b + kt (y) = b − kt (y) for y ∈ and ∂b ± kt /∂n x is just a linear scaling of b ± kt according to (48). In this case, statements for D can if desired be found from the sum and difference of (50) and (51). Images of these waves with envelope functions applied are included later in the paper in Fig. 10. The waves + kt and − kt are the pressure scattered by the individual basis functions b + kt and b − kt . Since the boundary quantities and ∂ /∂n y are approximated by a weighted sum of b + kt and b − kt , it follows that the total scattered wave will be approximated by a weighted sum of + kt and − kt , with the same weights as occurred in the boundary representation of . Examination of these is therefore interesting since they are the building blocks of the scattered wave. Moreover, they are the waves that appear inside the inner products of the various Galerkin integrals in Sections 2 and 3, hence are interesting since these are what are spatially integrated with a testing function to obtain entries in the interaction matrices. If they are zero over the support of the testing function, or if the testing function does not successfully demodulate their oscillatory behaviour, then that matrix coefficient will be very small or zero. This property will be investigated in subsequent sections.

Properties of the incident wave term
In this section we consider the evaluation of the right-hand side terms given by (26), (27) and (34). Note that since is infinite, all the boundary integrals over it will have to be evaluated in a limiting and normalised sense in order to yield finite results, as Fourier transform integrals typically are i.e.
dl. We continue to index wave direction by k t instead of q.
Let us assume that the incident wave ϒ (x) is a plane wave of complex amplitude A propagating in directiond, so ϒ (x) = Ae ikd·[x−x 0 ] and ∂ϒ/∂n x (x) = ikd ·n ϒ (x). Here the x 0 term in the exponent is not essential, but simplifies the statements below.
For comparison purposes, we first evaluate the right-hand side vector υ of the standard scheme (26).
Substituting these into a limiting and normalised version of (26) gives: Here δ (· · ·) is a Kronecker delta, which equals 1 when its argument equals zero and zero otherwise. It is clear from this that υ kt is only non-zero whend ·t = k t /k. This is satisfied when ϒ ( . This means that the inner product statement with an oscillatory test function is, as expected, performing beam-steering. It cannot however distinguish between o + kt and o − kt ; it is equally sensitive to both. In microphone-array terminology, (52) would be said to suffer from 'front-back confusion'.
Consider now the right-hand side of the Burton-Miller scheme (20). The second term υ ′ kt (27) is identical to υ kt except that it involves ∂ϒ/∂n . Hence with α = ik, the right-hand side of the Burton-Miller formulation equals: This is also only non-zero whend ·t = k t /k, so beam-steering is again being performed. Consider first the case wherê d ·t = 0, sod = ±n. Under these conditions, (53) is only sensitive to waves arriving surface normal from + ; waves arriving surface-normal from − cause it to equal zero. However, as the incident direction swings away from being entirely surface normal, the term 1 −d ·n becomes non-zero and the statement becomes sensitive to waves arriving from − , though it remains more sensitive to waves from + . This behaviour is consistent with the Cardioid microphone array analogy given in Section 2.1.1.
Finally, attention is turned to the Wave-Matching right-hand side vectorυ, as defined in (34). This statement also involves ∂b − kt /∂n x (x) = −ik n b − kt (x). Using the same process as above gives: As with (52) and (53), this statement performs beamforming and is only non-zero whend ·t = k t /k. Moreover, d ·n = √ 1 −d ·t becausen andt are perpendicular, which withd ·t = k t /k gives ⏐ ⏐ ⏐d ·n ⏐ ⏐ ⏐ = k n /k. This means that if d ·n > 0, which is true for any wave arriving from − , the bracketed term equals zero; it is therefore only sensitive to waves arriving from + . Hulsebos et al. give a physical interpretation of this in Fig. 7 of [54]. It is equivalent to choosing a microphone sensitivity pattern that has a null at the 'front-back confusion' angle, being the angle you are aiming to sense reflected in the plane of the array.
These observations allow (54) to be simplified:

Properties of the interaction matrices: self-interaction
We will now consider the integral that produces the entries in the matricesW ± . These are most easily analysed by considering the left-hand side of (31) rather than the four constituent terms from Section 3.2. For this we will require the surface-normal derivatives of (50) and (51): First, we consider same-plane (self) interaction. Modifying (31) to find the coefficient inW ± that would correspond to radiation from the basis function with tangential wavenumber k ty , being tested by an incoming basis function with tangential wavenumber k tx , gives: We again must evaluate this in a normalised limiting sense, like a Fourier transform, to get finite results. Since x ∈ − ⊂ − , we state it as x = x 0 + lt − ϵn, where x 0 ∈ as before and ϵ is some limitingly small positive real constant; we will take the limit as ϵ → 0. Noting from (50) and (56) (47) with (51) and (57) gives: This statement clearly includes a beam-steering effect, like the statements found in Section 4.1.2 did;W − kt x ,kt y = 0 unless k ty = k tx . Under this condition k nx = k ny , so for self-interactionW − kt x ,kt y = 2ik n when k ty = k tx , and zero otherwise.
These two results show that plane waves are orthogonal under the Wave-Matching BIE metric. This orthogonality applies both in terms of tangential wavenumber and whether the waves are incoming or outgoing. It means that for an approximation space including any set of incoming and outgoing waves with differing k t ,W + would equal zero andW − would be diagonal. Note that this does not imply that ϕ + is necessarily zero; it is found from ϕ − by (32). What it does however mean is that the iterative solver in Section 3.2.1 converges fully in a single reflection order r = 0; the outgoing wave + causes no higher-order reflection. Physically this manifests as Snell's law; at a plane boundary a plane wave reflects a single plane wave with equal tangential wavenumber and amplitude define by the reflection coefficient. No higher-order reflections occur.
Actually, the standard and Burton-Miller formulations would produce diagonal matrices for this self-interaction case too -it occurs due to the inner-product of the oscillatory basis functions taking the form of a Fourier transform. The advantages of the Wave-Matching operator will however become clear when non-self-interaction is considered; this is done in the following section.

Properties of the interaction matrices: non-self-interaction
Consider now two parallel planes with oppositely oriented normal vectors separated by distance d as depicted in Fig. 4. Fig. 4a depicts the case where they separate two acoustic media (white) by an internal domain (grey); this has the geometry of a classic transmission problem, but since the current scheme specifies no model describing propagation through 1− ∪ 2− it should form an impenetrable barrier. From a physical perspective, we therefore expect the two scattering problems in 1+ ∩ 2− and 1− ∩ 2+ to be decoupled. Fig. 4b in contrast shows the case where the two planes enclose an interior cavity; in this case they are expected to strongly interact.
Consider first the cavity in Fig. 4b. We assume that the boundary operators D, S, A and H are evaluated over 1 and the Galerkin testing integral is evaluated over 2 , though equal results would be achieved the other way around due to the symmetry of the problem. From (51) and (57) it is immediately apparent that the block ofW − representing interaction across the cavity equals 0 because 2 ⊂ 1+ . Substituting (47), (50) and (56) into (58) gives the following statement for the block ofW + representing interaction across the cavity. It is similar to (59), but some terms are negated becausen 2 = −n 1 : This features the beam-steering behaviour that has been widely seen, meaning this block ofW + will be diagonal. An iterative scheme would see the entries of the iteration matrix inv (W − )W + R equal to −e ikn y d R(k ny ). This, as expected, contains the phase-lag across the cavity and the reflection coefficient. The convergence of the scheme would follow the decay of reflection amplitude with reflection order so is, as expected, dictated by the absorption of the cavity. The standard and Burton-Miller formulations would work in an equivalent way for this scenario too. The inner product integral with the oscillatory basis functions would again produce diagonal matrices, and both schemes are sensitive to all waves arriving from above 2+ , which includes reflections from 1 , so reflections across the cavity would be solved for correctly. Differences will however occur for the scenario in Fig. 4a, where the objective is to avoid interactions across the fictitious interior cavity and to decouple the two problems on 1 and 2 .
Consider the blocks of the matrices of the Wave-Matching reflectance scheme that represent propagation across the interior cavity in the scenario in Fig. 4a. We now have thatW + = 0, as was the case for self-interaction, because 2 ⊂ 1− .
The entries ofW − are given by: Because of the relationship between k t and k n the bracket term will always equal zero whenever the Kronecker delta term is non-zero. This block ofW − therefore also equals zero, meaning the problems on 1 and 2 are completely decoupled. This property applies for the admittance version of the Wave-Matching formulation (Section 3.3) too since any radiation from a uniform admittance boundary can be converted to a weighted sum of + kt y and − kt y terms, for which these properties apply.
The situation for the standard and Burton-Miller formulation is, on the other hand, rather different for this scenario. Assuming the scattering from 1 is some weighted sum of + kt y and − kt y , then only − kt y will be present at 2 ⊂ 1− . For the same reasons discussed in Section 4.1.2, the interaction coefficient of the standard formulation will equal e ikn y d δ ( k ty − k tx ) . This will produce a diagonal matrix, but its diagonal will be unitary magnitude for all k t . Importantly there is no loss in the interaction, and for frequencies where e ikn y d = ±1 the standard scheme will not have a unique solution.
The coefficient for the Burton-Miller scheme will equal . This is again diagonal and is zero for perfectly normal waves with k t = 0, but for all other angles there will be some interaction. Importantly however, Here the denominator term arises in the corresponding Burton-Miller self-interaction coefficient, and this fraction represents the proportion of energy that is lost on each reflection across the cavity. This is smaller than one for all k and k t , so represents a mode that is decaying with reflection order. Such a mode is stable for a solver that marches on in time or reflection-order. For a time-harmonic model, it invalidates the time-invariant condition and should therefore be somehow rejected by the BIE. This is, from a marching on in reflection-order perspective, why the Burton-Miller formulation always gives a unique solution.
As a side note, it was observed in [32] that the Burton-Miller formulation is solvable and gives a unique solutions when α = −ik is chosen, despite going against all the microphone array analogy and energy arguments made in this paper. An explanation for this is that choosing α = −ik makes the interaction coefficient ratio described always greater than one i.e. growing with respect to time or reflection order. Such a mode would cause instability in a solver that marches on in time or reflection-order, but for a time-harmonic model, it invalidates the time-invariant condition just the same as a decaying mode would, so should therefore also be somehow rejected by BIE. Clearly this argument is less robust than the more mathematical ones given elsewhere, but it is felt by the authors that its intuitive physical nature gives it some merit. This section has demonstrated the superiority of the Wave-Matching BIE, from the perspective of avoiding internal resonances and mathematically decoupling problems that, from a physical standpoint, should be decoupled. The next section will demonstrate that these results apply to arbitrary-shaped open boundaries also.

Properties for arbitrary non-closed boundaries
Here it will be shown that the various results derived for the Wave-Matching scheme in Sections 4.1.1 to 4.1.4 also hold for arbitrary-shaped non-closed boundaries, subject to a few conditions. This is based on deformation of the canonical planar boundary. The associated plane wave basis functions, which were already defined for points off the original boundary in (47) and (48), are retained. The statement for reflection coefficient in (49) will however no longer apply due to the change of boundary geometry.
The first step is to acknowledge that the KHBIE, on which (50), (51), (56) and (57) are based, was derived from Green's theorem, so the result given by integration over would also be given by any other boundary so long as the two fields involved satisfy the Helmholtz equation for every point contained in the domain between and the new boundary. Here we are using the oscillatory functions o ± (y) to represent (y) on the boundary, so these must satisfy this condition; for plane wave o ± this will always be true. The Green's function is singular and does not satisfy the Helmholtz equation when x = y, so the domain between and the new boundary cannot contain x. This in practice means that the values of ± kt (x) will change for any x that is in-between the old and new boundaries.
The second step is to notice that the integrals in (31), on which all the coefficients in the linear system are based, also take the form required to apply Green's theorem. These will therefore also be invariant of a change to the boundary so long as both fields satisfy the Helmholtz equation in the domain in-between. o − (y) is used for testing so again this condition must hold for that. ± kt (x) are discontinuous across the radiating surface, so the testing surface cannot cross it, though it can get limitingly close.
The upshot of this is that all the results above still hold for non-planar surfaces: • For self-interactionW + is zero andW − is diagonal with entries equal to 2ik n .
• For interaction across an external cavityW − is zero andW + is diagonal.
• For interaction across an internal cavityW + andW − are both zero.
Note that in the first point k n still refers to the component of wavenumber perpendicular to some plane surface that the plane wave is defined with respect to; the oscillatory function must be kept the same for the above to hold i.e. the plane wave definition does not move with the boundary. The relations hold for evanescent waves too, but only when their oscillatory components are co-planar. A significant practical limitation that applies to this however, is definition of the reflectance matrix. This is straightforward for geometries on which the incoming and outgoing oscillatory functions are separable, but is not straightforward on arbitrary boundaries. This issue is likely to define how useful this extrapolation of operators to other boundaries process is in practice. To find it one would likely have to solve the scattering problem with the incident wave ϒ(x) set to match each incoming basis function, and then somehow find the component of each outgoing basis function present in the scattered field (x). Such a process appears to render the proposed method entirely useless, since the problem has to be solved by another method in order for it to be solved by this one. The data required is however rather similar to that which is measured by holographic techniques [65] and it is also likely, given the setup of these measurement techniques, that a smaller, geometrically local simulation of just the boundary section of interest could compute the local reflectance matrix approximately. Hence a set of smaller, geometrically-detailed models could feed scattering data into a larger, geometricallyidealised model based on the proposed method. This approach of replacing a geometrically complicated boundary with a simpler one and some scattering data is well established within geometrical acoustic methods [7]. Envelope functions have been notably absent from this argument, as they are from all of this section. We will see in Section 5.1 that they compromise many of these orthogonality type properties. One can however argue that if the oscillatory functions are chosen globally, then the properties above hold for the BIE globally; how the envelope functions partition this is just a question of how the BIE is 'chopped up'. It follows that the global solution should stay more or less the same if different envelope functions are defined, but that individual matrix coefficients might differ quite significantly.
This concludes the analysis of schemes with plane wave oscillatory functions but no windowing by envelope functions. Some numerical results will be presented for this case with envelope functions in Section 5.1. The next section concentrates on solution of the entire cylindrical canonical problem without envelope function.

Cylindrical boundaries
Here we consider a cylindrical boundary as defined in Section 3.1.2, and the polar coordinates (r, θ) will be used interchangeably with (x). The incoming and incoming and outgoing basis functions used to discretise + and − are: Here the surface-tangential wavenumber k t can be indexed by the integer parameter q as proposed in (16), because the boundary is of finite length 2π a; it will take values k t = q/a. The n subscript used in (16) is still dropped because envelope functions are not used in this section either.
The ratios α ± q , being entries from the discrete Dirichlet-to-Neumann map, are defined: These can be used with (49) to find the reflection coefficient that would occur if the boundary had admittance Y . They can also be used to find the Wave-Matching scheme matrices in Sections 3.2 and 3.3 from the standard Galerkin BEM matrices. This will be done in Section 4.2.6, when a normalised form of (62) will be considered. For now, we must recall that a standard BEM scheme with oscillatory functions would not include the radial term in (62) and would instead use b q = e iqθ .

Radiation of the incoming and outgoing waves by the boundary operators
Here we again consider the KHBIE in its interpretation as a sound-field rendering operator. We will see again that radiation of + is straightforward, whereas − requires more detailed consideration. Eq. (50) can be applied directly, though we restate it here with subscript q: This wave is depicted for ka = 30 and q = 2 in Fig. 5a.
Radiation of − is complicated by two factors: (i) the wave is radiated with negative amplitude, as for the plane wave in Section 4.1.1; (ii) because is closed, − will pass through − and then emerge into + . To understand what is happening it is helpful to first consider a regular oscillatory function o r Bessel function of order q. For the same reasons given when discussing the incoming plane wave, this radiates as: This wave is depicted for ka = 30 and q = 2 in Fig. 5b.
This wave is depicted for ka = 30 and q = 2 in Fig. 5c.

Orthogonality of inward and outward waves
Before we consider the values this choice of basis gives for the actual left and right-hand side statements of the Wave-Matching BIE, we will derive an orthogonality relation that will simplify their evaluation. Consider the following quantity, where (p, q) ∈ I and ( †, ‡) ∈ (+, −): HereȊ {· · · , . . .} can be thought of as the identity operator of the Wave-Matching scheme. It evaluates an integral of the same form as those in (31), but the quantity being 'tested' is another basis function. First the definition for o † p and o ‡ q has been substituted, and then the integral has been evaluated analytically, producing a Kronecker delta that is only non-zero when p = q. In the final statement it has been recognised that the bracketed terms cancel if † ̸ = ‡, or that they are a Wronskian (Eq. 10.5.5 of [66]) if † = ‡, which equals †4i/π ka.
(67) is the cylindrical wave equivalent of the relation for plane waves that was derived in (59). It shows that cylindrical waves are also orthogonal under the Wave-Matching BIE metric. This orthogonality applies both in order q, which also represents surface-tangential wavenumber, and whether the waves are incoming or outgoing. An equivalent property was demonstrated for spherical waves in 3D in [57]. Substitution of o r This latter term may also be proven by Green's theorem, since o r p (x) satisfies the Helmholtz equation for all x ∈ − and all p ∈ I.
These expressions will now be used to simplify analytical evaluation of the incident and self-interaction statements for the Wave-Matching scheme.

Properties of the incident wave term
In this section we consider evaluation of the right-hand side term defined in (34) as was done for the plane wave basis in Section 4.1.2, though here the integrals can be evaluated without modification since the boundary is finite. The incident wave ϒ (x) is assumed to satisfy the Helmholtz equation everywhere in the neighbourhood of the obstacles. This means it may be expressed as a weighted sum of the regular oscillatory functions o r q : Substituting (70) into in (34) produces a collection of termsȊ , from which it is apparent that υ p = −4iυ r p ; the Wave-Matching right-hand side integral has 'matched' each term of the expansion in (70) and extracted its coefficient υ r q , albeit multiplied by a constant factor −4i. This is the 'Wave-Matching' property that the BIE is named after applied to cylindrical waves. It also holds for the circular microphone array designs given in [54].
A potentially useful aspect is that this means the excitation vectorυ may be found by other means if that is easier than directly evaluating (34). For example if ϒ (x) were a cylindrical wave centred on a point in + , then rotation and translation formulae may be used to obtainυ directly. Similarly, if ϒ (x) is a plane wave thenυ can be found by the Jacobi-Anger expansion; for a plane-wave ϒ (r, θ) = e ikr cos θ travelling in the θ = 0 direction, this states that υ r q = i q , henceυ p = −4i p+1 .

Properties of the interaction matrices
We now consider the integral on the left-hand side of (31) that produces the entries in the matricesW ± . Here we are only considering one cylinder so these are effectively all self-interaction matrices. Noting that (31) is evaluated over − ⊂ − , (65) and its surface-normal derivative immediately tell us thatW + = 0. Substituting (66) into (31), it is apparent that This analysis is equivalent to the Fourier analysis performed in [26,28]. There a Fourier series was used as a convenient orthonormal basis; here it arises from the choice of approximation space. Indeed, these analytical results forW ± can also be derived using the operator identities given in those papers. In [28] Amini discusses how the choice α = ik is ''almost optimal'' for the Burton-Miller in that it ''almost minimises'' the condition number of the linear system produced. In terms of that criteria, we can state that the Wave-Matching BIE applied to the approximation scheme in (16) with N = 1, e 0 (x) = 1 for all x ∈ and oscillatory functions defined by (62), is exactly optimal in that it has a condition number equal to 1 for all k and any choice of boundary condition.
It is also interesting to compare how the scheme computes scattering from a cylinder compared to how an analytical model would formulate it. Such a scheme would typically use o r q to represent the incident field in the form of (70) and then o + q to represent the scattered wave. SinceW + = 0,W − = 8iI andυ p = −4i q+1 , we find that solving (33) or applying the iterative scheme in Section 3.2.1 immediately gives ϕ − =υ p (8i) −1 = i q /2. This is half the value expected in an analytical scheme using coefficients υ r q because o r The differing reflection coefficient definitions compensate for this in the values of ϕ + they produce. Specifically, in an analytical scheme the reflection coefficient takes the values − , whereas substituting (62) into (40) gives . We therefore see that, aside from some fixed scaling terms, the Wave-Matching scheme produces the same results that a simple analytical model would do, when no envelope functions are used. This may seem trivial, but it should be reassuring. As we proceed to introduce other complexities such as more complicated geometry and envelope functions, we know that we are at least starting from an algorithm that solves the problem in the way that one would choose to do it by hand.

Properties for arbitrary closed boundaries
Here it will be shown that the various results derived for the Wave-Matching scheme in Sections 4.2.1 to 4.2.4 also hold for arbitrary-shaped closed boundaries, subject to a few conditions. This is based on deformation of the canonical cylindrical boundary. The associated basis functions, which were already defined for points off the original boundary in (62) and (63), are retained. The statement for reflection coefficient in (49) will however no longer apply due to the change of boundary geometry.
This follows exactly the same logic as Section 4.1.5. The KHBIE, which computes the scattered wave terms, will give the same result for a difference choice of so long as and G satisfy the Helmholtz equation at every point in the domain in-between; this is true because it was derived from Green's theorem. The integrals of the Wave-Matching scheme (31) also take the form required for the application of Green's theorem, so will also produce results that are invariant of a change of boundary, subject to similar conditions. Since we only have one boundary, the restrictions due to G only mean that the boundary must remain Lipschitz and cannot cross itself. The cylindrical basis functions are singular at r = 0, so the origin of their polar coordinate system must remain inside the boundary, otherwise the values of the integrals will change. Actually, it is easily shown via Green's theorem that (67) and (68) equal zero if the origin of the polar coordinate system is not in − . (67) and (68) therefore still hold so long as this condition is added to the non-zero branches. In summary, cylindrical waves are orthogonal under the Wave-Matching BIE metric for any closed boundary that contains the origin of the polar coordinate system in which they are based.
As mentioned in Section 4.1.5 however, a significant practical limitation that applies to the use of this is the definition of the reflectance matrix. This is straightforward for geometries on which the incoming and outgoing oscillatory functions are separable, but is not straightforward on arbitrary boundaries; finding this would likely amount to solving another BIE instead. This issue is likely to define how useful this extrapolation of operators to other boundaries process is in practice. Having said that, the ability to replace an integral over a complicated boundary with an integral over a simpler one may be attractive in some circumstances.

Normalisation of basis functions
The choice of basis functions given in (62) has produced some remarkable results when used with the Wave-Matching BIE. Unfortunately however, there are some numerical issues hiding in the detail. These are associated with the magnitudes of the basis functions, which are scaled by H ± q (ka). The magnitude of these terms is plotted in Fig. 6, and can clearly be see to explode for q > ka. Returning to the derivation of the orthogonality relation in (67), it is clear that these should cancel out in the Wronskian. However for q even quite modestly above ka subtraction error begins to occur, meaning the matrix coefficients become inaccurate and finally the scheme collapses entirely. This limits the useable values of Q and therefore the accuracy that can be achieved.
Another issue with the definition in (62) is connected to the idea that we aim to 'join up' sections of canonical problems, as illustrated in Fig. 1. Having the plane-wave basis functions unitary and the cylindrical basis functions non-unitary adds unnecessary complication; they could be 'joined up' more easily if they were both unitary and equal on the boundary. In addition, a more serious problem to do with the orthogonality relation above is also aggravated when the basis functions are not normalised. Eq. (67) also relies on the orthogonality of the terms e iqθ when integrated over a range 0 to 2π . The introduction of envelope functions compromises this orthogonality, meaning there will be off-diagonal terms when p ̸ = q. For these the terms of the Wronskian do not cancel, so very large basis function magnitudes will produce very large coefficients.
Normalising the basis functions also makes sense from the point of view that the Wave-Matching BIE is attempting to satisfy an energy metric. Waves with q > ka are largely evanescent and propagate relatively little energy; it can be seen from Fig. 3 that their pressure and particle velocity will be out-of-phase with one another. However in (67), they produce the same interaction matrix coefficients as the non-evanescent waves; this is because they are larger amplitude for the same value of discretisation weight. Hence a normalised scheme makes more sense from an energy perspective too; evanescent waves have smaller coefficients because they cause less energy flow. The normalisation is straightforward. Instead of (62), the basis functions are defined: This makes the basis functions unitary and equal on the boundary; that is b + The definition for α + q in (63) still holds with this revised definition. Numerical experiments have however shown that the definition of α − q must be changed to align with that which occurs for plane waves, otherwise extremely poor conditioning occurs. We no longer use α − q = α + q ; instead we take α − q = −α + q . The reason for this poor-conditioning is visible in Fig. 3.
As ka and q increase, the imaginary part of α + q decays away so α + q ≈ α + q . This produces pair of incoming and outgoing wave terms that are almost identical, hence the scheme is ill-conditioned. More mysterious perhaps is why this problem does not manifest for the non-normalised scheme; this is a question that has not yet been answered and should ideally be by future research. Taking α − q = −α + q also makes sense from an energy flow perspective, since negating the particle velocity for the same pressure will negate the power flux through the boundary.
With regard to how changing to choose α − q = −α + q affects − q (x), this shown by Fig. 7. Here ka = 30 and q = 32; a value of q > ka is necessary to make the differences visible. It is apparent that the old choice in Fig. 7a is radiating most energy outwards, which is unsurprising since it tends towards equalling + q (x) as ka and q increase. In contrast, the new choice in Fig. 7b is radiating more energy inwards than outwards. + q (x) is shown in Fig. 7c; there is no question that this is radiating purely outwards and that the value of α + q is correct.
A peculiarity of this modification is that α − p = α + p must however still be used for the testing scheme. Again, the exact reason for this is an open question. One possible explanation comes from the time-reversed back-extrapolation presented in [54]. This suggests that the wave radiated by the KHBIE with the Green's function conjugated gives an impression of the incident waves that it is sensitive to. This idea requires more detailed consideration, but it is certainly true that a wave coming from + , as a time-reversed version of + would be, is what the testing integral is designed to 'look' for.
With these modifications, and defining a diagonal matrix α with α qq = α + q , the expressions for the reflectance scheme in Section 3.2 and the admittance scheme in 3.3 can be found from the standard Galerkin BEM matrices in Section 2.2 by: Here the matrices with the superscript ± are for the reflectance scheme and those without are for the admittance scheme (which uses the standard Galerkin BEM basis functions for representation of the boundary quantities and only uses b − q for testing). These statements hold for both cylindrical and plane wave basis function, subject to the differing definitions of α + q , and will be applied in Section 5 when envelope functions are introduced too. The statement for reflectance of plane wave basis in (49) holds with these normalised cylindrical basis functions too. It follows that the Wave-Matching admittance scheme is equal to the reflectance version when Y = 0. We remark that α is the discrete form of the Dirichlet-to-Neumann map and that similar uses occur in the OSRC literature [34], though there it is more common to invert the map and apply it to H.
With regard to what effect these changes have on the matricesW ± , the new choice of α − p means the identities in [28] are required. The details will not be presented here, but it can be shown thatW + = 0 still holds andW − is still diagonal, but has changed fromW − = 8iI toW − = 4π α. This means that all the most significant properties from the previous section still hold, and that the matrixW − is still trivial to invert, but its condition number will no longer equal one. Note that since these values are predictable, it would be possible to use preconditioning to compensate for them. This would however invalidate the energy argument made above for why normalisation of the basis functions makes sense, and it may also cause numerical problems since α + (k t ) = 0 for k t = k on a planar boundary (we note that this issue also arises in the OSRC literature -see e.g. (12) of [39]). The extrapolation to non-cylindrical boundaries also still holds -we are still working with all the same incoming and outgoing wave building blocks and have simply rescaled them -but the basis functions will not be unitary on a boundary where r ̸ = a.
Conditioning results are shown for these schemes in Fig. 8. Results are shown for Y = 0 (a-d) and Y = 1 (e-h), to show how the boundary condition can affect conditioning and accuracy, and are presented for 1.5k ≤ Q ≤ 2.5k, to cover the region where k t approaches, equals and passes k in the medium; this is the point at which convergence begins. The conditioning of the Burton-Miller method is shown in a & e. It gives small condition numbers that vary with the boundary condition. The conditioning of the Wave-Matching reflectance scheme is shown in b & f. It has the highest condition numbers of those  (72) and (73), but gives different results for Y = 1; these are the lowest of any found for this test case. Finally, the bottom row shows the normalised L 2 error for the two boundary conditions. Only one figure is presented for each case because all the schemes gave the same error; it is obviously dictated by the approximation space. Conditioning results were also calculated for the standard operator but they are not included here because they were heavily corrupted by interior resonances, so it was impossible to see any trends. Condition numbers were however much larger -often three orders of magnitude so more than the other schemes.
This demonstrates that the two versions of the Wave-Matching scheme both give correct results and are both immune to interior resonances, as the Burton-Miller method is. In more cases than not, the Burton-Miller method gave slightly lower condition numbers; clearly the increased complexity of the Wave-Matching formulation and the range of values on the diagonal ofW − are not quite delivering the benefits expected based on the analysis in Section 4.2.4. We will now move on to considering the performance of these schemes when envelope functions are introduced.

Scheme for windowed canonical problems
In this section we examine the properties of the new BIE when envelope functions are introduced. The spatial windowing that they enable is an essential ingredient of both POU-BEM [12], where it means a set of oscillatory functions that do not exactly match the incident wave or its scattering may still approximate it, and HNA-BEM, where differing oscillatory functions will be defined on sections of the boundary that are 'illuminated' differently [13,14,23].
The same canonical scattering problems will be considered -a planar barrier between two domains and a cylinder -but they will be spatially windowed. The same set of pre-defined oscillatory functions will also be retained; hence a POU-BEM type scheme arises. Regarding this choice, we remark that, unlike many of the HNA-BEM papers, choice of approximation space is not the key focus of this paper. Clearly a properly considered HNA-BEM approximation space has many advantages, but for simplicity we choose to use something predefined.
The first choice that must be made is what envelope function to use. POU-BEM requires that the envelope functions be C 0 continuous [12], but Peake et al. [67] showed that versions with C 1 continuity produce more accurate results. The original intention here was to use staggered copies of Peake's N 2 shape function, which is a raised cosine function; this is capable of forming a POU scheme on its own without the corresponding N 1 and N 3 . It was however discovered that this caused the numerical integration of the hypersingular operator H to converge very slowly, since (28) includes db/dl terms that, for a C 1 envelope, are only C 0 . To mitigate this, a piecewise-polynomial C 2 POU envelope function was designed; it is depicted in Fig. 9 and given by: where f (µ) = 8µ 3 − 8µ 4 . µ is a normalised edge coordinate, defined relative to the centre of the basis function's associated boundary element. This element is considered to occupy the region − 1 2 ≤ µ ≤ 1 2 , so it is obvious from Fig. 9 that the basis functions will overlap and span multiple elements; this is an inherent feature of a POU scheme.
The parameter γ scales the width of the transition region. This was included so that its optimal length could be investigated. Certainly the idea of 'chopping' up canonical problems in Fig. 1 implies an distinct transition without overlap, but some degree of smoothness is necessary from both a numerical and theoretical perspective. This issue is of course ignored by geometric methods operating under asymptotically high ka assumptions.
With γ = 1 the envelope function has no flat-top -the transitions are as gradual as the support of the basis function allows. With γ = 0 the transition is instantaneous, so the envelope is not even C 0 . Such discontinuities are acceptable for the standard scheme -use of piecewise-constant discretisation is common -but any scheme that uses the hypersingular operator H requires basis functions with C 1,α continuity [68]. Moreover, the weak-singularity form in (28) assumes a continuous field in its derivation [41], and γ = 0 means db/dl would be infinite for µ = ± 1 2 . We therefore make the restriction that γ > 0.
It is also important to consider that multiplication in space is equivalent to convolution in k t . This means we expect some degree of spectral 'smudging' due to multiplication by envelope functions. The sharper their transitions, the 'whiter' their wavenumber spectrum, and the greater the amount of 'smudging' that will occur. This will be visible in the results in Section 5.1. The same envelope functions are currently used for pressure and pressure gradient, which, as briefly mentioned in Section 3, may not be the optimal choice. In particular the relationship between the envelope for o ± n,q (x) and ∂o ± n,q /∂n x (x) will affect the approximation of the Dirichlet-to-Neumann map, and in the OSRC literature this usually includes surfacetangential derivatives of the basis functions (often via the Laplace-Beltrami operator [10]). This suggests that the surfacetangential derivative of the envelope function should be taken into account somehow too, but this is not attempted in this initial set of case studies and it remains an open research question.
It is also visible in the pressure field the basis functions radiate; this is depicted in Fig. 10. Here, γ = 0 has been used for subplots a-d, with γ = 0.2 and γ = 1 being used in e and f respectively. Note that in these figures the envelope has been stretched so that its support always fills the entirety of the black line. The first four plots (a-d) show how the radiation from the double-layer and single-layer potentials is respectively a sum and difference of the radiation from the inward and outward basis functions. This point has until now only been articulated through equations; it is discussed in greater detail in [16]. The last three (d-f) aim to show the effect of the envelope function on the radiated wave. In d there is strong diffraction visibly emanating from the limits of the element, where the hard transitions occur. In e, introducing a smooth transition has already made quite a significant difference to the amount of corner diffraction, even though the transition region is still quite narrow. In f, distinct corner diffraction is no longer visible; there is still spreading in the beam, but it occurs further from the element. The discussion in Section 4.3 of [16] points out that this makes it more similar to a geometric beam than d is. In wavenumber spectrum terms, f has been subjected to less spectral 'smudging' than d or e.
Having discussed the principles and implications of windowing, we now move onto the two case study problems. Here the process is that the standard Galerkin BEM matrices in Section 2.2 will first be calculated, which also makes those schemes available for comparison, and then the matrices of the Wave-Matching scheme will be found from them using the statement in Section 4.2.6. Finally, L 2 error is computed versus an analytical model of scattering from a cylinder. The integration was done numerically by brute force. First, the oscillation rate of the integration kernels was calculated allowing the integration domain to be split into regions smaller than half a kernel cycle, then an order 11 Gauss-Legendre scheme was applied. This is high enough to ensure integration accuracy does not significantly affect results, but gives this demonstration algorithm a total quadrature cost that scale O , which, since Q must in practice scale with k (see Fig. 8), corresponds also to O Clearly this is not scalable to large problems; as Section 1.1 already identified, efficient integration of these kernels is a key challenge that needs to be addressed.

Planar boundary case study
Here we revisit the infinite barrier problem that was depicted in Fig. 4a and was analysed in Section 4.1.4. Again, we are looking for the interaction across the interior cavity − to be minimised, both to avoid internal eigenmodes and to ensure the two exterior problems are decoupled.
Meshing two infinite boundaries is obviously not possible, so rather than attempt to numerically solve the problem, we will just examine the interactions between some nearby elements on opposite sides of the boundary. Results have been computed for one receiving element and nine radiating elements; all are of length L and the width of the interior domain (d in Fig. 4a) is also L. The elements were aligned so that the middle radiating element was directly below the receiving element. The envelope function in (79) was used and the oscillatory basis functions were defined as: Here is has been assumed that there is a unique reversible mapping between x and µ that depends on the index n of the element that it is centred on. The oscillatory part of the basis function is implicitly extrapolated for x / ∈ following the form of (47). Noting thatt · [x − x 0 ] = µL if µ = 0 corresponds to x = x 0 , we therefore have that k t = 2π q/L.
The magnitudes of the interaction matrix coefficients are shown in Fig. 11. The entries in the Burton-Miller and Wave-Matching matrices have been normalised by 2ik in line with (53) and (59) when k t = 0. This compensates for the scaling caused by the p = q = 0 self-interaction coefficient when the matrix is inverted.
Comparison between a-c, for which γ = 0.2, and d-f, for which γ = 1, makes clear the spectral 'smudging' effect discussed above. It is particularly apparent for the standard formulation in a and d, since the interactions are much stronger.
It might seem surprising that the strongest interactions lie on the p = −q diagonal instead of p = q; this is simply because the elements are rotated opposite ways since they are on opposite sides of the obstacle, so have opposing normal and tangent vectors. The standard operator gives values close to 1 2 down the diagonal. It is no coincidence that this is the value of its boundary trace -it is also the magnitude of its 'beam' in the high kL limit, as discussed in [16]. The Burton-Miller scheme successfully suppresses all energy arriving from the −n direction as expected from the discussion in Sections 2.1.1 and 4.1.2; this manifests as the coefficients being zero around p = q = 0. However away from that zone it is sensitive to energy arriving accross the interior cavity. This is most notable for the non-adjacent case h, where a 'beam' from the offset radiating element For (a-f), the receiving element is directly above the radiating element. For (g-i), they are offset by one element length. Plots (a, d, g) show data for the     intersects the receiving element at an angle the Burton-Miller formulation is sensitive to. Finally it is clear that the Wave-Matching scheme (c,f,i) gives the smallest magnitude interaction coefficients for all cases. This further reinforced by Fig. 12a, which takes the maximum over all radiating elements and values of p and q for a range of kL. The maximum coefficient in the standard formulation tends to 1 2 as k increases as expected. The Burton-Miller scheme offers a significant improvement but it is clearly the Wave-Matching formulation that is most immune to energy crossing the internal cavity. Fig. 12b shows that it is not just the maximum values that is small in magnitude; the vast majority of coefficients are far smaller. This is the sparse diagonal form ofW ± in section 4 manifesting in an approximate sense for windowed basis functions too.
These results have shown the superiority of the Wave-Matching BIE when used with plane wave basis functions on planar geometry. A cylindrical test-case will now be investigated.

Cylindrical boundary case study
The circular boundary will be divided into eight equal sized elements, so N = 8 in (16) and L = 2π a/N, where a is still the radius of the cylinder. The basis function definition from (79) and (80) will be retained, though the oscillatory part of the basis function is implicitly extrapolated for x / ∈ following the form of (71). The characteristics of the waves that these basis functions radiate is shown in Fig. 13 for ka = 30 and q = 2; for N = 8 this corresponds to k t a = Nq = 16. Note that the black line indicates the entire support of the basis function. It is seen that the incoming (a) and outgoing (b) waves propagate in the directions intended; note that definition α − q = −α + q from Section 4.2.6 is used here because the oscillatory functions are unitary on the boundary. For the incoming wave it is clearly seen how the radiated wave enters − and then leaves to enter + , with some focusing occurring in − due to the curvature of . This is rather more intuitive to understand than Fig. 5c, where the waves in − and + appear unrelated. It is also noted that for the incoming wave there is some 'leakage' outwards through the boundary; it is thought that this might be due to using the same envelope functions for pressure gradient as for pressure.
It should be noted that the approximation space created by this definition is incomplete; it cannot represent the incident field or solution for some incident wave directions. Using 1i in the exponent of (80) resolves this issue and was demonstrated to be accurate, but produces an approximation space that is so badly conditioned that it is impossible to see any difference in the condition numbers due to the different boundary operators. Very high condition numbers are the norm in POU-BEM and specialist matrix solvers are required to overcome them [40], but that does not help with this analysis. Using 2i, as is instead done here, gives better conditioning since all the basis functions with a particular value of m are orthogonal when multiplied by a single envelope function. Non-orthogonality arises from the squaring of the envelope function in self-interaction, and the overlap of basis functions on adjacent elements; this is why the condition numbers increase as γ increases. Clearly using such an approximation space is extremely limiting, but for the purpose of characterising the boundary operators in this case study it is satisfactory.

Convergence and conditioning
Convergence versus Q is shown for ka = 25 and Y = 0 in Fig. 14. Part (a) shows the condition number for the identity (mass) matrix M. It is seen to be significantly higher with γ = 1 compared with γ = 0.2, and increases with Q ; both of these aspects are expected based on the preceding discussion regarding the approximation space. This causes the conditioning of the boundary operators in (b) to also be much higher than was seen in Fig. 8 for the unwindowed but normalised scheme in Section 4.2.6. The condition number is slightly higher for the Wave-Matching scheme than for the Burton-Miller scheme with both values of γ . In (c), the conditioning trends are seen to be offset against the solution accuracy that can be achieved; smoother approximation spaces appear to give better convergence. Intermediate values of γ were also computed for all results in this section but are not shown since they lie predictably between the γ = 0.2 and γ = 1 trends.
It was demonstrated in Section 4 that the Wave-Matching scheme without envelope functions is immune to the nonuniqueness problem that affects standard (13) and surface-normal derivative (14) forms of BEM. Arguments were made in Section 5.1 that this should also be the case for the windowed version with plane wave basis.   16 presents the conditioning and accuracy of the scheme versus ka, with Q = 2 ⌈ka/2⌉ + 1. Results are presented for γ = 0.2 and γ = 1 and, unlike most results in this section, for Y = 1 as well as Y = 0. For the Y = 1 case conditioning is shown for the reflectance and admittance versions of the Wave-Matching scheme, since they differ slightly. Normalised L 2 error is shown only for the Wave-Matching scheme since all schemes achieved almost exactly the same accuracy, except for the occasional glitch from the standard scheme; clearly this is limited by the approximation space. Condition numbers are seen to increase with ka, so Q also, and γ for all schemes, whereas accuracy improves. It appears from c and f that Q = 2 ⌈ka/2⌉ + 1 might be increasing with ka faster than is required to maintain a predetermined accuracy. Presently as ka increases, L 2 error decreases.

Interaction matrix sparsity
It was seen in Section 4 that all the interaction matrices for the non-windowed scheme are diagonal and therefore extremely sparse; they have Q 2 entries but only Q of these are non-zero. In Section 5.1 it was seen that a windowed scheme   Population nnz(W ) that must remain to achieve normalised L 2 error <0.01 versus ka, plus trendlines. Both plots show data for Y = 0 and γ = 1. with plane wave oscillatory functions also produces matrices that are 'approximately sparse', in that the vast majority of their entries are extremely small. Fig. 17 shows that this is also the case for a cylinder with cylindrical oscillatory functions, and that it holds more strongly for larger ka and for the envelope function with the more gradual transition; this latter property it related to the diffraction and spectral 'smearing' effects discussed in Section 5.1. Results are only shown for the Wave-Matching scheme, but the Burton-Miller scheme actually gave almost identical results, with those from the standard scheme only slightly higher; clearly on a cylinder the beamforming effect of the basis functions is a more significant factor in dictating this, compared to the ratio of sensing of pressure and its surface-normal derivative.
We are most interested in the size of the dark area at the bottom of the graph. Bigger is better, as this means there are relatively few larger magnitude coefficients representing significant interactions. This leads naturally to the question of whether the smaller coefficients can be discarded to make the matrices truly sparse without adversely affecting accuracy of conditioning. Fig. 18 shows the results of applying this process to the interaction matrices used to calculate the Wave-Matching results in Fig. 16. In Fig. 18a, all matrix entries below some threshold are discarded, then the % of entries remaining is plotted again the condition number of this 'culled' matrix and the solution accuracy it achieves. This is shown for ka = 25, and it is clearly seen that around 75% of entries may be discarded before any adverse effects occur. Fig. 18b is the result of computing results like Fig. 18a for a range of ka, and then tracking the location of the black dot -the point past which further culling would cause normalised L 2 error to rise above some threshold, here chosen to be 1%. It is seen that while the total number of entries (numel) of the matrix scales by definition with Q 2 and therefore (ka) 2 , the number of non-zero entries (nnz) that must be retained instead follows a O(k) trend -an order better. This shows that the propagation operator is approximately diagonal, which is consistent with the beamforming interpretation that has been given and with what geometric algorithms assume at large ka. How to find these significant interactions a priori, so only they need be calculated, is an open question. One approach might be to use geometric algorithms and/or an approximate model of diffraction to estimate them.

Iterative solution by marching on in reflection order
Finally, we come to results for the iterative solver that, as suggested in the introduction, is interesting because it may allow adaptive finding of optimal wave directions at each reflection order. As stated in Section 3.2.1, it is the largest eigenvalues of the update matrix inv (W − )W + R that defines the stability and convergence of the series. This is plotted for Y = 0 in Fig. 19, and it is seen to be comfortably less than one for all k studied. This is also not significantly affected by the choice of γ , though γ = 0.2 actually converges slightly faster. Figures for error versus iteration number are not shown here, but the solver typically reaches the accuracy achieved by matrix inversion within 15 iterations.

Conclusions and further work
This paper has proposed a new boundary integral equation for solution of the scalar Helmholtz equation, termed the 'Wave-Matching' BIE. It was derived from energy flow principles adapted from time-domain BEM, but uses a decomposition of the field into incoming and outgoing waves, something that was shown to be straightforward with the oscillatory basis functions chosen. The properties of the new scheme were first analysed on two canonical scattering problems without introducing envelope functions. Here the combination of inner product terms used in the testing integral was shown to have a 'Wave-Matching' property; given a trial wave (the testing function), it will return a coefficient proportional to the amplitude of that wave component within whatever field it receives. This form was also identified for use in microphone arrays in [54] and has previously been presented by the authors of this paper in three dimensions in [57]. This led to a set of orthogonality relations between the incoming and outgoing waves, which makes the interaction matrices diagonal. It was also shown to entirely eliminate interactions across the interior cavity of a planar barrier, decoupling the problems in the two exterior domains and eradicating the possibility of internal cavity resonances leading to non-unique frequencies. Moreover, it was seen that: (i) the interaction matrices for the cylinder with non-normalised basis functions is a scaled version of the identity matrix, so has condition number equal to one; (ii) these properties can be extrapolated to a much wider class of boundary geometry using Green's theorem.
Following this, envelope functions were introduced and the same canonical problems reconsidered. A full BEM model of an infinite planar boundary was of course not possible, but a test case examining the interaction matrices across a planar barrier was performed, and the Wave-Matching BIE was seen to outperform both the standard and Burton-Miller formulations in terms of its rejection of cross-cavity interactions. A full model of a cylinder was then examined in more detail. The new scheme was shown to achieve the same accuracy as the standard and Burton-Miller formulations, with conditioning that was almost as good as Burton-Miller and which was unaffected by internal cavity eigenfrequencies. The matrices produced were seen to be sparse, and discarding small coefficients was seen to be reduce the number of entries from O ( k 2 ) to O (k) while maintaining a required accuracy threshold. Finally, a marching-on-in-reflection-order iterative solver was implemented and seen to be accurate and stable for all cases considered. In terms of future work, there is much to do. This has been an expansive exposition of the new BIE and its properties, but it is still just an initial study and a versatile and practical algorithm is some way off. Mathematicians studying such problems would typically bring more advanced analysis techniques to bear on a new formulation e.g. checking for existence and uniqueness of solutions [35] and for coercivity [48,69]; this has yet to be done. There are also a few open questions within this paper require answers e.g. an explanation for the choice of α − q found to be necessary in Section 4.2.6. It is also possible that using different envelopes for pressure and pressure gradient, as the OSRC literature suggests, may give better results. Coupling to other technologies from this field, e.g. OSRC preconditioners, would be interesting to investigate, as would coupling to fast-multipole.
Engineers will want to see the algorithm applied to more realistic geometries, for which 'joining up' of the canonical problems will be necessary as discussed in Section 1.1. Whether this is possible and how best to do it is currently unknown; a first step might be to join up the two problems studied herein to make the torpedo problem in [40]. Following that, consideration of how to model wedge diffraction in this framework, as depicted in Fig. 1, would be a high priority, since it would enable the modelling of polygons, as has been successfully done with HNA-BEM [13,23]. The algorithm would ideally also be taken into three-dimensions. Some parts of the framework for this are already in place [16,57], but many challenges remain, hence this paper has restricted itself to two-dimensional problems. Efficient integration is a challenge that restricts the effectiveness of all these applications, but much work is being done to address that. How scattering by geometrically complex boundaries can best be represented using this method is also an open question.
There are however many interesting qualities that the new scheme has that have great potential and justify the research outlined above. The idea of adaptively finding optimal wave directions while marching on in reflection order has the potential to achieve the gains of HNA-BEM while avoiding the limitation that classes of obstacle must be studied in great detail by hand. The similarity of this approach to geometric algorithms in the large ka limit suggests that fully-integrated wide-bandwidth hybrid algorithms may be possible. This is of great interest in acoustics for example, where the audible bandwidth covers many octaves and applications such as auralisation require compatible modelling of all of this. Marching on a full wave solution for a few reflection orders until this scattering becomes chaotic, and then transferring that energy into a phase-less energy BEM such as [18][19][20], has the potential to be an efficient solution.