On hybrid convolution quadrature approaches for modeling time-domain wave problems with broadband frequency content

We propose two hybrid convolution quadrature based discretizations of the wave equation on interior domains with broadband Neumann boundary data or source terms. The convolution quadrature method transforms the time-domain wave problem into a series of Helmholtz problems with complex-valued wavenumbers, in which the boundary data and solutions are connected to those of the original problem through the  -transform. The hybrid method terminology refers specifically to the use of different approximations of these Helmholtz problems, depending on the frequency. For lower frequencies, we employ the boundary element method, while for more oscillatory problems, we develop two alternative high frequency approximations based on plane wave decompositions of the acoustic field on the boundary. In the first approach, we apply dynamical energy analysis to numerically approximate the plane wave amplitudes. The phases will then be reconstructed using a novel approach based on matching the boundary element solution to the plane wave ansatz in the frequency region where we switch between the low and high frequency methods. The second high frequency method is based on applying the Neumann-to-Dirichlet map for plane waves to the given boundary data. Finally, we investigate the effectiveness of both hybrid approaches across a range of numerical experiments.

The convolution quadrature (CQ) approach was developed by Lubich [2][3][4] and provides a simple way to obtain a stable time stepping scheme using the Laplace transform of the boundary integral kernel.Since the initial work of Lubich, significant advances have been made in the CQ method for time-dependent wave modeling including high-order Runge-Kutta implementations for a variety of wave equations. 5,6CQ possesses favorable stability properties due to an implicit regularization in time.For the discretization of boundary integral equations (BIEs), one of its main advantages is that it avoids having to evaluate the convolution kernel in the time domain and instead involves solving a simplified system of frequency domain BIEs in the spatial region.[10] We employ a direct boundary element method (BEM) for the spatial discretization of the lower frequency Helmholtz problems since the BEM requires us to discretize along the boundary only, whereas other methods such as finite element method and finite difference method, require us to discretize over the whole domain.This reduces the dimensionality of the problem and hence the complexity of the calculation.Other examples of work on the CQBEM for solving a range of problems can be found in References 11-18.Of particular relevance to the present study is the work of Reference 18, which proposes a similar hybrid approach to the one introduced here.The methods proposed in Reference 18 are for exterior wave scattering problems where typically multiple wave reflections do not have to be considered, particularly for convex wave scatterers.In the work here we apply a high frequency approximation for interior wave problems that also includes multiple reflections, and compare with a simpler approach that ignores them.Rather than wave scattering, the applications of the study are therefore mainly in architectural acoustics with natural extensions to modeling wave fields in enclosed electromagnetic environments or in linearly elastic structures.
For transient broadband wave sources or boundary data, the CQBEM requires large numbers of spatial boundary elements and very small time-steps to model the high frequency part of the signal, and thus can quickly become infeasible due to a large computational cost to obtain an appropriate accuracy.A number of statistical or ray based methods have been developed specifically for modeling wave problems at high frequencies, including statistical energy analysis (SEA), ray tracing and dynamical energy analysis (DEA).SEA is the most well-known approach for modeling high frequency noise and vibrations, examples can be found in References 1, 19, and 20.In SEA, a structure is split into subsystems and ergodicity of the underlying ray dynamics as well as quasi-equilibrium conditions in each subsystem are assumed.That is, the energy density in each subsystem is assumed to be approximately constant with Lambertian reflections at boundaries.These assumptions result in a simplified set of SEA equations based only on flow rates between subsystems. 21However, the particular sub-division of the structure is critical to the validity of the underlying assumptions, which are often hard to verify.In recent years, SEA has also been developed for transient problems 22,23 where the success of the approach relies on the success of the corresponding steady-state analysis.
5][26] At the heart of the method is a linear integral operator based model for transporting densities along ray trajectories in phase-space, between intersections with the boundary of a domain or sub-domain.In recent years the capability of the method has been extended to large-scale problems from industry through an efficient implementation on finite element type meshes (in the position variable). 24,25In this study, we make use of the recently proposed Petrov-Galerkin approach to DEA introduced in Reference 27, where the directivity is modeled using a finite global set of Dirac delta distributions that can model highly directive problems to high accuracy. 27The hybrid method developed in this work will be more computationally efficient than CQBEM, since at high-frequencies we can solve in terms of a frequency independent phase-space density.However, this phase-space density does not include phase information and so here we propose a novel approach to reconstruct the phase by matching up the solutions from the BEM and DEA in the frequency region where we switch between these two methods.We additionally propose the use of an alternative simple high frequency approximation (SHFA) formula based on a physical optics approximation that omits reflected wave contributions.This is similar to the high frequency approximation employed in Reference 18, but here is formulated for interior boundary value problems (BVPs), as opposed to exterior wave scattering.
The remainder of this article is structured as follows.In Section 2, we detail the initial-BVP for the wave equation under consideration and its reformulation as a second-kind BIE.Section 3 then presents the CQ time discretization procedure, and in particular, the reformulation of the time-domain BIE introduced in Section 2 as a set of Helmholtz BIEs at complex-valued wavenumbers k.In Section 4, we then outline the methodologies employed to solve these Helmholtz problems, depending on the magnitude of Re(k), or equivalently, the frequency.A BEM scheme for low frequencies is detailed in Section 4.1, while two alternative high frequency approximations are described in Section 4.2 and then the reconstruction of the solution in the time-domain from the set of Helmholtz solutions is reported in Section 4.3.We then present numerical results for a variety of polygonal domains in Section 5, where the problem is driven by either nonhomogeneous boundary data or and interior point source, before drawing our conclusions in Section 6.

INTEGRAL FORMULATION OF THE WAVE EQUATION
Let Ω ⊂ R 2 be a finite domain with regular boundary Γ = Ω.We note that the models and discretization methods proposed in this work all have natural extensions to three-dimensional applications, albeit with an associated increase of computational cost that could be reduced through, for example, the use of fast BEMs in conjunction with iterative matrix solvers-see Reference 18.We will consider the following Neumann BVP for the inhomogeneous wave equation with initial conditions and boundary condition for some T > 0. Here, we assume f and g are scalar and real-valued functions of space and time, c > 0 is the wave speed and n is the unit outward normal to the boundary.We will consider problems where internal wave sources described by g undergo reflections at rigid boundaries (f ≡ 0), or where g ≡ 0 and the boundary Γ corresponds to an interface with a vibrating structure that generates an inhomogeneous boundary condition f .In order to treat the inhomogeneous wave equation using a boundary integral formulation, we introduce an incident solution v and a reverberant solution u in the whole of R 2 , and write  = (u + v)| Ω .We assume that v solves the free-space problem with initial conditions and ĝ is an extension of g to the whole of R 2 .In fact, we will consider only spatially localized sources of the form for some prescribed temporal source profile g 0 ∶ (0, T) → R and a spatial source point x 0 ∈ Ω.Hence, the extension of g from Ω to the whole of R 2 is simply an extension by zero.The solution to (2a) and (2b) may be conveniently expressed in terms of the fundamental solution of the wave equation G via the convolution theorem as where and H is the Heaviside step-function.Note that we have made use of the form of g specified in (3) to derive the solution (4) and the temporal integral may be numerically evaluated using an appropriate quadrature rule.With the incident solution v in hand, then it is easy to see that the reverberant solution u must satisfy the homogeneous wave equation, and since f ≡ 0 in (1c), then the boundary condition for u is given by The BVP for u therefore takes the same form as (1a)-(1c) when g ≡ 0 and f = −v∕ n.
In the remainder of this section, we outline a boundary integral formulation for the wave equation (1a) with g ≡ 0, together with initial conditions (1b), and a general inhomogeneous Neumann boundary condition (1c).The solution of this initial-BVP may be expressed by the following direct boundary integral representation formula for x ∈ Ω, which results from Green's second identity. 28Here,  and  are, respectively, the single and double layer potential operators Considering now the limit as the solution point x approaches a point on the boundary where Γ is locally differentiable, then the representation formula (7) becomes a second-kind BIE for  on Γ.Here we have used V and K to denote the traces of  and  onto Γ, respectively, and I is the identity operator.In the next section, we initiate the discretization process for the BIE (9) by considering its semi-discretization in time.

TIME DISCRETIZATION VIA CQ
In this section, we outline the semi-discretization of the BIE (9) in time using Lubich's CQ method. 2,4In particular, through the application of CQ we are able to transform the space-time BIE (9) into a series of frequency domain problems with complex wavenumbers.We do not recall the theoretical framework here but summarize the application of the method.First, we note that boundary integral operators V and K in (9) are time convolution operators and make use of the notation to emphasize this.Note that Ṽ( t )f is standard notation for V * f in the CQ literature. 4,12Here Ṽ(s) and K(s) are the Laplace transforms of V and K, respectively, and s is the Laplace transform parameter.We note that Ṽ(s) and K(s) correspond to the traces onto Γ of the single and double layer potentials, respectively, for the Helmholtz equation with wavenumber k = is∕c.Explicitly, for x ∈ Γ the operators Ṽ(s) and K(s) are given by where F and Φ denote the Laplace transforms of f and , respectively.In addition, G k is the fundamental solution to the Helmholtz equation in two dimensions given by with H (1) 0 being the zeroth order Hankel function of the first kind.To discretize the time convolution ( K( t ) ) , we split the time interval [0, T] into N steps of equal length Δt = T∕N and compute an approximate solution at the discrete time steps t n = nΔt.The continuous convolution operator K( t ) at the discrete times t n is replaced by the discrete convolution operator for n = 0, … , N − 1, where  j =  Δt (⋅, t j ) and the superscript Δt is used to denote quantities that have been semi-discretized in time.The convolution weights are defined by their -transform The function (z) is the quotient of the generating polynomials of the multistep method used to discretize in time.In our numerical examples, we consider only the second order backward difference formula (BDF2) where (z) = 1 2 (z 2 − 4z + 3).The convolution weights can be calculated using the inverse -transform via an approximation of Cauchy's integral formula using the trapezoidal rule.The approximate convolution weights are then given by a scaled inverse discrete Fourier transform as where C is taken as a circular contour centered at the origin of radius  < 1 and In this work, we allow the choices of N and Ñ to be decoupled in order to potentially over-resolve in the Laplace domain for better accuracy as proposed in Reference 16, although typically we choose Ñ = 2N as recommended in Reference 3. We now show that applying the semi-discretization ( 14), together with an analogous time discretization procedure for ( Ṽ( t )f ) , in the integral equation ( 9) leads to a system of integral equations for the Helmholtz equation with complex wavenumbers k l = is l ∕c for l = 0, 1, … , Ñ − 1.By extending the sum in (14) to j = N − 1, substituting in the approximate weights (16) and then applying in the BIE (9), we obtain a new system of equations for  Δt, (x, t n ) =   n (x): Substituting the definition of w Δt, j (16) into (17), then multiplying by  n and applying a discrete Fourier transform with respect to n gives where are the -transforms of  Δt, , f , respectively.We have thereby reduced the problem of numerically solving the wave equation to numerically solving a system of Helmholtz equations with complex wavenumbers k l = is l ∕c, for l = 0, 1, … , Ñ − 1.An example of the range of wavenumbers k l is shown in Figure 1.In the next section, we describe the spatial discretization of these Helmholtz problems via a conventional BEM as well as using a high frequency approximation for cases when Re(k l ) is large relative to a typical length scale of Ω.

SPATIAL DISCRETIZATION IN THE FREQUENCY DOMAIN
For the spatial discretization, we either employ a piecewise constant collocation BEM, or for a high frequency region to be specified below (in terms of Re(k l )) we employ a high frequency approximation.Note that due to the reflective symmetry of the wavenumbers k l , l = 0, 1, … , Ñ − 1 in the imaginary axis (see Figure 1), then only Ñ∕2 + 1 Helmholtz problems are approximated (when Re(k l ) ⩽ 0) and the remaining Helmholtz (approximate) solutions are obtained through simple complex conjugation (when Re(k l ) > 0).We specify a threshold k * for which we employ the BEM when |Re(k l )| ⩽ k * and let  ⩽ Ñ∕2 be the minimal integer valued index of the minimal |Re(k l )| > k * , which is the region for which we apply the high frequency approximation.Considering the structure of the wavenumbers k l , l = 0, 1, … , Ñ − 1 shown in Figure 1, and noting that the indexing l = 0, 1, … , Ñ − 1 starts where Re(k 0 ) = 0 at the bottom of the loop and runs clockwise, then we specify the location of the wavenumber k  to be within the lower left quarter of the loop of possible k l values-see Figure 1.This choice is a consequence of the fact that the criteria for (heuristically) choosing the threshold k * is based only on the modulus of Re(k l ), which controls how oscillatory the solution of the Helmholtz problem is.The high-frequency approximations will only be viable when |Re(k l )| is sufficiently large and if one continues around the loop of wavenumbers k l beyond the lower left quarter then |Re(k l )| begins to decrease again.The reason to choose the lower half of the loop instead of the upper half is heuristic and is discussed further at the end of this section.We also note that the choice between positive or negative Re(k l ) for the location of k  is entirely arbitrary since, as we mention above, the wavenumbers k l , l = 0, 1, … , Ñ − 1 are symmetric about the imaginary axis and the solutions are simply complex conjugates.Here we choose negative real parts simply for convenience, since those wavenumbers are listed first as l is increased from zero.
The choice of the value for the threshold k * is essentially a trade-off between computational cost and accuracy.Choosing the threshold too low will result in faster calculations with poor accuracy due to applying high frequency approximations at low frequencies where they introduce significant errors.On the other hand, a very large choice of k * will mean that most of the calculation is performed using the CQBEM without any noticeable reduction in accuracy and only a minimal reduction in the computational time.A specific choice for the k * value could therefore be made in much the same way as you might choose a discretization parameter in a numerical method, for example, N or M in the CQBEM scheme here.That is, start with a relatively low value and increase until reaching an accuracy level that you deem acceptable.The effect of changing the value of k * on both the accuracy and computational cost will be investigated numerically in Section 5, where a range of k * values between 25 and 350 will be considered.It is also worth noting that the values of k l , l = 0, 1, … , Ñ − 1 are inversely proportional to the time-step Δt and so for a fixed threshold k * , then decreasing Δt has the effect of increasing the number of Helmholtz problems approximated by a high-frequency approximation, while the number approximated using the BEM remains approximately fixed. 18he high frequency approximation will be provided either by a numerical method based on a plane-wave approximation, or by an incident illumination approximation where only the direct contribution of the source term on the boundary is considered and reflected contributions are assumed to play an insignificant role.The latter approximation will be reasonable when the imaginary part of the wavenumber is large and the ray trajectories in the corresponding geometrical optics model are sufficiently long to allow for significant decay before any reflections occur.The numerical approach based on a plane-wave approximation will be able to go beyond the incident illumination model in terms of the reflection order, but will introduce additional sources of error owing the numerical discretization and phase reconstruction procedures.The plane wave amplitudes in a finite set of directions are calculated using DEA, with a recently proposed Petrov-Galerkin discretization. 27The phases are constructed by matching the high frequency approximation with the BEM results in the frequency region where we switch from BEM to a high frequency approximation.In relation to the location of the wavenumber k  mentioned above, it is advisable to perform this matching using Helmholtz problems with wavenumbers having smaller imaginary parts and hence use values from the lower half of the loop of wavenumbers shown in Figure 1.The reason for this is heuristic and based on the fact that Helmholtz solutions for wavenumbers with larger imaginary parts will consequently include greater levels of decay that could potentially give rise to conditioning issues.

BEM in the low frequency region
The low frequency BEM discretization is performed by dividing Γ into M subintervals (or elements) E m of approximately equal size.We assume that the -transformed numerical solution Φ l and boundary data F l can be approximated by We substitute (19) into the integral equation (18) and solve for the transformed solution coefficients Φ l,m .The resulting system of equations may be written explicitly in the form for each

High frequency approximation
We will consider two alternative high frequency approximations.First, let us consider a plane wave superposition of the form where  l = Re(ck l ) and In the first high frequency approach, the amplitude terms A  are approximated using DEA as described in the Section 4.2.1.In fact, the DEA method computes a stationary phase-space density  on Γ, which is then related to the amplitude A  and the phase S  via: where p ∈ R 2 is the slowness (or momentum) vector satisfying |p| = c −1 .In particular, for a plane wave whose direction relative to the x-axis is defined by the angle Θ, then c p = (cos Θ, sin Θ).The approximation of the phase terms S  will then be performed by matching the BEM solution (the first of ( 19)) to the expression ( 22) at l =  and l =  + 1 using the amplitude terms calculated via DEA and leaving only the phase terms to be determined.This procedure will be described in Section 4.2.2.A second high-frequency approach based on a simple approximation formula will then be described in Section 4.2.3.

DEA approximation of the amplitudes A 𝛼
In DEA, phase-space densities are transported through Ω using a boundary integral form of the Frobenius-Perron (FP) operator, 26 which is expressed in terms of boundary phase-space coordinates X = (q, p), where the arclength parameter q ∈ [0, |Γ|) corresponds to the position x ∈ Γ and |Γ| is used to denote the total length of the boundary.The momentum coordinate p is used to denote the tangential component of the outgoing slowness vector p at the point q.A density  is then transported from the phase-space boundary coordinate Y = (q ′ , p ′ ) to the next intersection with the boundary X via the phase-space boundary integral operator where  is the Dirac delta distribution and the map  is the boundary map for Γ.In particular,  transports the boundary coordinate Y = (q ′ , p ′ ) along a ray path starting from the position specified by q ′ and traveling in the direction specified by p ′ until reaching Γ again-see Figure 2. The result (Y ) = ( q (Y ),  p (Y )) is a new boundary phase-space coordinate whose position  q (Y ) is determined from the intersection of the ray with the boundary and whose tangential slowness  p (Y ) corresponds to a specular reflection of the incoming ray once it reaches the boundary.The factor w in (24) contains a damping term and, if necessary, any reflection/transmission coefficients.For the damping only case, where d is the Euclidean distance metric between a pair of positions on Γ specified in terms of arclength parameters.The boundary map (Y ) is relatively simple to determine for convex domains since for each ray direction p ′ , oriented from q ′ into Ω, there is a unique intersection point where the ray arrives at the boundary again.To treat non-convex domains directly using DEA would require one to omit any rays that travel outside Ω, possibly using a visibility function.Often a simple alternative, particularly for complex geometries, is to decompose Ω into a set of convex sub-domains whose union is Ω as described in References 27 and 29.
For frequency domain problems, we need to calculate the stationary boundary density  induced by an initial boundary density  0 via The boundary map  acting on the boundary phase-space coordinates Y = (q ′ , p ′ ) ∈ Γ × (−c −1 , c −1 ).The coordinate Y denotes a ray leaving the boundary at position q ′ and in direction  ′ = arcsin(cp ′ ) with respect to the normal vector to the boundary Γ at the position q ′ .The result of the mapping is denoted , which corresponds to the next intersection with a boundary edge where the ray undergoes a specular reflection.The arrows represent the vector p at the boundary positions q ′ and  q (Y ), which is represented by only its tangential component, p ′ and  p (Y ), respectively, in the boundary coordinate system as annotated on the figure where  () refers to the th iterate of the operator .The initial density term  0 is related to the amplitude of the prescribed boundary condition F l for the Helmholtz problem under consideration-specific details will be provided in Section 5 for each type of boundary data considered.We calculate an approximation to the density  by solving a finite dimensional approximation of (26).A number of techniques for this based on a Galerkin projections have been proposed, which have the advantage of circumventing any difficulties associated with the Dirac  in ( 24) by introducing additional integrals with respect to X, that are then evaluated simply using the properties of .Here, we adopt the Petrov-Galerkin scheme proposed recently in Reference 27, which uses a standard piecewise constant Galerkin method in the position variable q and a Petrov-Galerkin approximation in the tangential slowness variable p.The piecewise constant spatial approximation of the amplitudes is then consistent with the one given by the BEM in the low frequency range and the same spatial meshes will be used in each case.The Petrov-Galerkin approximation in the p variable uses Dirac delta distributions for the basis, which correspond to a finite set of globally consistent directions, with a test function space of piecewise constants.We now give brief details of this discretization process, but for full details the interested reader is referred to Reference 27.
Consider a set of global ray directions Θ l ∈ [0, 2), where l = 1, 2, … , L, defined anti-clockwise relative to the positive x-axis.Here we choose equispaced angles Θ l = 2(l − 1)∕L, but they could also be chosen in an example dependent way to include any known dominant transport directions.Now define   (q) ∈ (−∕2, ∕2),  = 1, 2, … , P m to be the local ray directions at q ∈ Γ.The local directions correspond to the subset of the global directions that are directed into Ω at q and have been relabeled according to the angle they make with the interior normal vector at q-see Figure 3 for a simple example case.We approximate  on Γ × (−c −1 , c −1 ) using a finite-dimensional approximation of the form where p (q) = sin (  (q)) ∕c and bm (q) = |E m | −1∕2 for q ∈ E m , and zero elsewhere, with |E m | = diam(E m ).For m = 1, 2, … , M, bm defines an orthonormal basis of piecewise constant functions with respect to the standard L 2 inner product, onto which we apply a standard (Bubnov) Galerkin projection in the position variable q.However, in the tangential slowness variable p we choose a set of test functions   that are orthogonal (in fact orthonormal) in the L 2 inner product to (p − p (q)) for  = 1, 2, … , P m , that is , where   =  − Θ +1 for  = 1, 2, 3, and the global coordinate index has been shifted to run over the correct subset of directions.The constant  is the global direction of the inward normal vector at q, and so  = ∕2 here unless  =  ′ , when the inner product is one.We do this by taking   (p) = χ (arcsin(cp)) = χ () = 1 if  ∈ I  , or zero otherwise.Here I  ,  = 1, 2, … , P m are a set of disjoint subintervals satisfying In particular, ) and A Petrov-Galerkin projection of the operator  on to the basis and test function combination described above leads to a matrix representation B with entries given by where the multi-indices J = (m, ) and J ′ = (m ′ ,  ′ ).The four-dimensional integral in the Galerkin projection has been reduced to a single integral over the boundary element E m ′ ⊂ Γ due to the local support of the spatial basis and the properties of the Dirac  distributions arising in Equations ( 24) and (27).The calculation of B J,J ′ is relatively simple because bm and   are locally constant and will be zero unless the direction  = arcsin(cp) ∈ I  and s ∈ E m , meaning that most entries to the matrix will be zero.The calculation when the matrix elements are nonzero only involves the integral of the exponential term (25) over the element E m ′ and multiplication by the pre-factor (|E m ||E m ′ |) −1∕2 .Furthermore, for polygonal boundaries the Euclidean distance function d(q ′ ,  q (q ′ , p ′ (q ′ ))) is linear in q ′ ∈ E m ′ and hence the integral in the third line of ( 29) can be performed analytically (after subdivision when either of the inputs correspond to corners).The coefficients of the expansion ( 27) can be found by solving the linear system  = (I − B) −1  0 , which corresponds to the discretized form of Equation (26).Here  0 and  represent the coefficients of the expansions of  0 and , respectively, when projected onto the finite dimensional basis.Note that the invertibility of I − B is a consequence of the fact that for A-stable CQ time discretizations such as the BDF2 scheme employed here, we have Im(k l ) > 0 and hence the damping term (25) will ensure that the leading eigenvalue of the operator  is smaller than 1 in absolute value.The entries of the source vector  0 corresponding to an initial density  0 are given using the property (28) and orthonormality of the spatial basis via We note that while the total number of matrix entries in B scales as (L 2 ), a particularly appealing feature of the Petrov-Galerkin discretization introduced in Reference 27 and summarized above is that its natural sparsity means the number of nonzero matrix entries, and hence the assembly cost, only scales linearly in the number of ray directions L. This sparsity was not present in earlier DEA schemes, such as those detailed in References 24, 26, and 29.For the hybrid methods described in this work, we consider the number of spatial boundary elements M to be fixed sufficiently large so that the BEM solution of the Helmholtz problems for all |Re(k l )| < k * can be determined accurately in some sense.In the calculations presented in Section 5, we adopt the commonly applied rule of thumb 30 and specify good accuracy as given by choosing around six elements per wavelength.This choice of M should also be more than sufficient for the DEA calculations at higher frequencies since the unknown is a relatively slowly varying amplitude in comparison to the oscillating Helmholtz solutions determined using the BEM.We will see in the next section that fixing M also provides limitations on the number of directions that can be included in the DEA-based hybrid scheme.Nevertheless, the combination of linear complexity in L with fixed M provides the realistic prospect of extending the proposed hybrid methods to three-dimensional cases where the global direction set would be naturally described using spherical polar coordinates.Assuming (L) directions subdividing both the polar and azimuthal angle would give a matrix with (L 4 ) total entries, but only (L 2 ) assembly cost.The solution of the resulting large sparse nonsymmetric linear systems would rely on the availability (or development) of an efficient iterative solver as mentioned at the start of Section 2.

Wave matching approximation of the phase terms
When Γ is subdivided using the same set of M piecewise constant elements in both the BEM and DEA, we can find the amplitude A  for each  directly from (26) by combining ( 23) and (27), and choosing the directions in the sum over  in (23) to correspond to those of the discretization represented in (27).However, estimating the corresponding phase term S  requires the development of new methodology.We propose a method to reconstruct the phase terms from a full wave solution at both the maximal frequency BEM calculation performed before we switch to the high frequency formulation and the lowest frequency at which we apply the high frequency formulation.Denoting these frequencies as  −1 = Re(ck −1 ) and   = Re(ck  ), respectively, we apply both BEM and DEA to obtain to a set of equations of the form for l =  − 1 and l = .We note that the left side of ( 31) is provided by the BEM solution and the linear form of the phase terms S  (x) = sin(  )q∕c +  l  expressed in ( 31) is a consequence of the fact that the wave speed c is taken to be constant.The constants  l  ,  = 1, 2, … , P, l =  − 1,  are the unknowns to be determined by imposing (31) at a set of points q on Γ.The variables   ,  = 1, 2, … , P are the local angles with respect to the unit normal introduced in the previous section.To calculate the undetermined constants  l  , we calculate the solutions to the Helmholtz equations using BEM at the frequencies  −1 and   corresponding to the pair of frequencies where we switch from using the BEM for the low frequency content to using DEA for the high frequency content.
The procedure must be performed at more than one frequency due to the nonuniqueness of the phase solution at a single frequency owing to the periodicity of the plane waves.The phase terms at  −1 and   may then be related by Solving (32) for  ∈ Z allows us to recover a unique set of phase constants via for either l =  − 1 or l = .Once the   values are known, we approximate the Helmholtz solutions Φ l for all frequencies with absolute value larger than | −1 | via (22) with S  (x) = sin (  ) q∕c +   .The specific details of the calculation of the  −1  (and    ) values will be example dependent.For the polygonal domains considered here the values will need to be calculated for each edge of the polygon separately, since each edge will have a different subset of the global directions Θ l , l = 1, 2, … , L associated to it (i.e., those pointing into Ω) and so the interpretation of the    values will be different on each edge.We generate a system M equations of the form (31) by choosing x = x i for i = 1, 2, … , M as the collocation points from the BEM approximation in (31).The next task is the determine how many of the amplitudes A  are nonzero at every collocation point x i on a given edge, since this provides a reduction in the number of phase constants    that we need to recover.The system of equations (31) can then be solved as a linear system for i = 1, 2, … , M and where q i is the arc-length parameter for the coordinate x i .The unknowns v l  = exp ) may be determined using the Moore-Penrose pseudo-inverse to obtain the least squares solution and from the result one can directly calculate the phase constants  l  = −i log Since this matching process is only performed for two frequencies, the computational cost of this part of the scheme is relatively small in comparison to the total cost.It involves solving a linear system with P unknowns for each boundary edge.Note that the number of local directions P pointing into Ω from a given edge is linearly related to the total number of global directions and typically L ≈ 2P, depending on the alignment of the boundary edge relative to the ray directions.Since we fix the total number of spatial boundary elements M as discussed in the previous section, this imposes a practical limitation on the number of directions that can be employed; taking P significantly larger than the number of boundary elements along the edge, that is the number of equations, will adversely affect the accuracy.

Incident illumination high frequency approximation
In this section, we describe a simple high frequency approximation (SHFA) based on the observation that the wavenumbers k l typically have large imaginary part-see Figure 1.The dissipative factor (25) appearing in the operator  will then damp out all contributions except those from very short ray trajectories meaning the approximation  ≈  0 will, in many cases, be reasonably good.There is no need to perform any DEA calculations for this approach and the solution for a wavenumber k l with a large enough imaginary part can be approximated by simply rescaling of the -transformed boundary data F l , corresponding to the Neumann-to-Dirichlet map for plane waves.In particular, we set , where  0 (x) defines the direction of the source term at x ∈ Γ relative to the normal direction.For example, for a point source from For a BVP with boundary data related to a plane wave entering the domain from one or more edges, then the angle  0 can be found directly from the plane wave direction.

Construction of the time-domain solution
Once the -transformed boundary solution Φ l , l = 0, 1, … , Ñ − 1 has been computed using the methods described in Sections 4.1 and 4.2, the space-time discretized solution   j can then be approximated via the trapezoidal rule for computing the inverse -transform as The interior solution is also calculated by applying the same time and spatial discretization to (7).We derive the Laplace domain interior solution Φ Ω l , at any point x inside the domain as An analogous inverse transform to (35) is used to approximate the interior solution (x, t) from the Laplace domain interior solution (36).We will present numerical results for both the interior and the boundary solutions for a selection of examples in the next section.

NUMERICAL RESULTS
In this section, we consider the numerical solution of inhomogeneous Neumann BVPs where the boundary data corresponds to a plane wave traveling into the domain, as well as homogeneous Neumann problems where the excitation of the system is driven by an internal point source.

Plane wave boundary data
We consider a plane wave traveling from left to right in a direction specified (in global coordinates) by Θ 0 ∈ [0, ∕2) relative to the positive x 1 -axis.For accuracy reasons, we limit our choice of Θ 0 to one of the global directions (from the DEA discretization) within the first quadrant, that is Θ 0 = Θ l for one of l = 1, 2, … with l < L∕4.We note that the directions Θ l can be chosen in a problem specific fashion and so this choice is not a limitation of the method.We apply this boundary condition on a unit square and L-shaped domain as shown in Figure 4. Explicitly, the corresponding boundary conditions in the original wave problem (1c) are taken as Here, x 1 and x 2 are the entries of the vector x, so the first condition specifies the left edge of the domain and the second specifies the lower edge.The function W is given by the normal derivative of a Gaussian pulse of the form ) , where n 1 and n 2 are the entries of the vector n and x ∈ R is given by either x 2 sin(Θ 0 ) − ct or x 1 cos (Θ 0 ) − ct as specified in (37).The parameters t 0 > 0 and  > 0 control the position of the peak of the Gaussian pulse and its bandwidth, respectively.Note that these parameters should be chosen carefully to ensure that the pulse has decayed sufficiently at t = 0 for the initial conditions to be approximately satisfied.
The combination of regular geometric features and the uni-directional boundary condition mean that throughout this section we will only need to use L = 8 global directions in the DEA implementation as illustrated in Figure 3.The boundary condition (37) leads to a DEA initial density  0 of the form where p 0 = sin( 0 )∕c and  0 is the local direction of the plane wave inducing the boundary data, relative to the normal vector to Γ at position q.That is, on the left edge (x 1 = 0) we have  0 = Θ 0 and on the lower edge (x 2 = 0) then  0 = −Θ 0 + ∕2.Recall that F l is the -transform of f and k l = is l ∕c.The numerator of the initial density (38) therefore represents the square amplitude of the -transformed boundary data transported along the direction of the incident plane wave.The denominator in (38) results from applying the Neumann-to-Dirichlet map to F l (and then taking the square amplitude) so that the phase-space densities in DEA relate to the square amplitudes of the wave functions rather than their normal derivatives-see (23).

Unit square domain
We first consider the case when Ω is a unit square and choose Θ 0 = 0 for simplicity.In this case, for early times (before any reflections) we have a simple exact solution that we can use to test the accuracy and convergence of the proposed methods.For simplicity, we fix c = 1 and choose  = 1∕64 in order to obtain a broadband signal.Note that max | (x, 0) | and max | t  (x, 0) | for x ∈ Ω ∪ Γ are obtained at x 1 = 0. Choosing t 0 = 0.1 then gives that to four digits, max | t (x, 0)| = 1.627e − 18 and max | t  (x, 0))| = 1.333e − 15, meaning that the initial conditions are approximately satisfied throughout Ω.A similar analysis can be performed to find the maximal time for which the exact solution is valid by considering the behavior when x 1 = 1, since the solution will only be valid during the early time period for which  ≈ 0 at x 1 = 1 and hence no reflections have occurred.Using the symmetry of the solution about t = t 0 one can see that the equivalent limitation for the maximal time would restrict the validity to t ∈ [0, 1].
Figure 5 shows a comparison between the exact and numerical interior solutions at x = (0.2, 0.25).We apply a high frequency approximation whenever |Re(k l )| > 350 and apply M = 1024 boundary elements to provide a good level of accuracy up to the BEM cut-off frequency  −1 .The plots compare the results of using the SHFA with the DEA based plane wave approximation up to t = 2, as well as comparing both these approximations with the exact solution up to t = 1.In this case, both high frequency approximations produce identical looking results, and one observes an improvement in the resolution of the reflected peak at t = 1.9 when using N = 8192 or N = 4096 time-steps compared to N = 2048.For earlier times t ⩽ 1, it is difficult to distinguish the seven curves in the plot.Calculating the relative errors for t ⩽ 1 via we obtain to four digits Error(x) = 0.06229 when N = 2048, Error(x) = 0.01609 when N = 4096 and Error(x) = 0.004101 when N = 8192.The errors for both DEA and the SHFA are identical up to the number of digits stated and we notice that the method is approximately converging at second order as would be expected for the BDF2 scheme, meaning that M = 1024 has been fixed large enough for the temporal discretization error to dominate.
In Figure 6, we show the solution along the boundary using both the DEA and SHFA high frequency approximations with the same parameter choices as before in the case N = 4096.We notice that the two solutions are visually identical over the entire boundary (as well as at the interior point x = (0.2, 0.25) investigated previously).The plots depict the Gaussian pulse moving over the left edge of the square (3 < q < 4) at around t = t 0 = 0.1 and then traveling along the upper and lower edges until reflecting on the right edge at approximately t = 1.1 and then returning back along the upper and lower edges.
Figure 7 shows the -transformed boundary solution Φ l .The left subplot (A) shows the absolute value at l = , where we switch from using BEM to calculate Φ l and instead use a high frequency approximation.The values of |Φ  | calculated using the BEM are compared to both of the proposed high frequency approximations.We observe that the main differences between the three methods are along the upper (2 < q < 3) and lower (0 < q < 1) edges where the high frequency approximations give zero and the BEM does not, since both high frequency approximations of Φ l omit waves tangential to the boundary.In the zoomed panel, we highlight that the main difference between DEA and the SHFA is apparent on the right hand edge (1 < q < 2) where DEA includes reflected wave contributions but the SHFA does not.The right subplot (B) shows Re(Φ l ) for wavenumbers close to Re(k) = 350.In this plot, Re(Φ l ) is computed using the DEA based approximation for Re(k) > 350, and using the BEM otherwise.Along the top edge (2 < q < 3), the value of Re(Φ l ) jumps to zero when we switch from the BEM to the DEA approximation as would be expected from the left subplot (A).Along the left edge (3 < q < 4), the solutions from the two methods match up very well demonstrating the success of the phase reconstruction process outlined in Section 4.2.2.
In the results so far we have demonstrated that the developed hybrid approaches can both provide accurate results, but owing to the choice of wavenumber threshold k * = 350, the computational cost saving relative to using BEM for all frequencies is relatively modest.The corresponding full CQBEM calculation was found to take approximately five times longer than the DEA based hybrid calculation and around six times longer than the SHFA based calculation when N = 4096, but the cost saving improves as N increases, and for the SHFA in particular since the cost is approximately (1) in N as for the high frequency approximation introduced in Reference 18.In Table 1, we investigate the effect of lowering the wavenumber threshold k * with final time T = 1 in the range where the exact solution is valid.We test the idea that as the maximal wavenumber in the BEM calculation decreases then so can the number of boundary elements employed for its solution, since the minimal wavelength that we need to resolve in the BEM increases and we use this minimal wavelength to choose the spatial element size as described at the end of Section 4.2.1.We find that in order to retain a reasonable accuracy level it is necessary to over-resolve in the Laplace domain by increasing the total number of Helmholtz problems to solve Ñ while fixing the number of time-steps N as proposed in Reference 16 and discussed earlier in Section 3. We notice from Table 1 that this has the effect of keeping the number of Helmholtz problems solved using the BEM reasonably stable, but cost savings are made through the reduction in the number of boundary elements M. We find that reducing the wavenumber threshold to k * = 175 has only a minimal effect on the error but gives a significant decrease in computational time and decreasing further gives more moderate improvements in computational efficiency while causing more significant increases in the error.We note also that the phase reconstruction process involved in the DEA based solution fails when the wavenumber threshold is reduced to k * = 45, leading to significantly larger errors than for the SHFA.For the larger choices of wavenumber threshold considered, the performance of the DEA and SHFA based schemes was very similar, but with the DEA approach requiring more computational resources as would be expected.Note that the computational times quoted are for the entire calculation using non-optimized MATLAB codes and significant speed-up would be expected for an optimized implementation in a C++ or Python, for example.
Figure 8 shows the results of performing the same numerical experiments again but with the wave direction changed to Θ 0 = ∕4, although in this case we only show the DEA based result for the boundary solution versus time in subplot (B), since both the DEA and SHFA results have an identical appearance as before.In subplot (A), we notice that the interior solutions at x = (0.2, 0.25) are visually identical, even for the N = 2048 case that was visibly less accurate when Note: As the maximal BEM wavenumber is decreased, there is also a proportionate decrease in the number of boundary elements together with an increase in the number of Helmholtz problems solved overall.Note that the number of Helmholtz problems solved using the BEM remains approximately unchanged. 0 = 0.The boundary solution in subplot (B) depicts the Gaussian pulse entering along the lower and left edges of the square (0 < q < 1 and 3 < q < 4) at the origin (q = 0 ≡ 4 (mod 4)) around t = t 0 = 0.1 and then traveling along those edges until reaching the vertices of the right edge (q = 1 and q = 3) at around t = 0.1 + √ 2∕2 ≈ 0.81 where a small part of the wave reflects, but mostly the propagation then continues along the right and top edges.The two boundary contributions then meet at the top-right vertex (q = 2) at around t = 0.1 + √ 2 ≈ 1.5 where there is a highly localized peak in the amplitude before the wave reflects to return along the right and top edges.The bottom row of subplots shows the -transformed boundary solution Φ l and corresponds to Figure 7 for Θ 0 = 0.The left subplot (C) shows the absolute value at l =  as before.In this case, there are no edges where the wave direction is tangential to the boundary and hence the agreement between the BEM solution and the high frequency approximations is much better than for Θ 0 = 0.The most significant divergence between three models occurs along the upper and right edges (1 < q < 3) where the SHFA prediction is zero due to not accounting for the reflected wave contributions.The DEA result shows a smooth continuation of the solution across all edges, whereas the BEM solution exhibits a small jump and oscillation close to the bottom-right (q = 1) and top-left (q = 3) vertices.The right subplot (D) shows a plot of Re(Φ l ) for wavenumbers close to Re(k) = 350.The solution is only shown for one of the dominant regions of Φ l near the lower-left vertex.The solution is computed using the DEA based approximation for Re(k) > 350, and using the BEM otherwise.We observe that the solutions from the two methods match up very well, again demonstrating the success of the phase reconstruction process outlined in Section 4.2.2.

L-shaped domain
We now repeat the above analysis for the L-shaped domain shown in Figure 4B and again initially focus on the case Θ 0 = 0 where we can make use of the same exact solution (39) for early times (before the wave reaches the re-entrant corner).As discussed previously in Sect.4.2.1, the DEA approximation process needs to be modified for non-convex domains such as an L-shape and the calculation is actually performed on a sub-divided version of the domain where each of the (two) sub-domains is convex.In this case, the sub-division was implemented by introducing an (artificial) internal interface connecting the vertices at q = 1 and q = 3 (see Figure 4B) to form two convex quadrilateral sub-domains.The extension of DEA to multi-domains is a straightforward process and we omit the details here for brevity-the interested reader can find details in References 27 and 29.We also mention for completeness that we must omit any amplitudes associated with the internal interface from the DEA result, and reorder the degrees of freedom to be consistent with the low frequency BEM calculations before implementing into the CQ algorithm.
Figure 9A shows a comparison between the exact and numerical interior solutions at x = (0.2, 0.25) and for T = 2.We again apply a high frequency approximation whenever |Re(k l )| > 350 and use M = 1024 boundary elements.As for the unit square, we observe that both high frequency approximations produce identical looking results, as well as an improvement in the resolution of the reflected peak at t = 1.9 when using N = 8192 or N = 4096 time-steps compared to N = 2048.Calculating the relative errors for t ⩽ 0.5 using (40) (but with the sums only until N∕4) we obtain to four digits Error(x) = 0.06229 when N = 2048, Error(x) = 0.01611 when N = 4096, and Error(x) = 0.004221 when N = 8192.The errors for both DEA and the SHFA are identical up to the number of digits stated and the values themselves are either the same as or very slightly larger than for the unit square.The method is again approximately converging at second order as expected for BDF2, meaning that M = 1024 has been fixed large enough for the temporal discretization error to dominate.In Figure 9B, we show the solution along the boundary computed using the DEA based high frequency approximation with the same parameter choices as before in the case N = 4096.The DEA and SHFA based approaches give visually identical solutions over the entire boundary and so just one of the plots is shown for brevity.The plot in Figure 9B shows the Gaussian pulse moving over the left edge of the L-shape (3.5 < q < 4) at around t = t 0 = 0.1 and then traveling along the upper and lower edges.At around t = 0.6, the pulse along the upper edge partially reflects at the re-entrant corner, whereas along the lower edge the pulse continues until approximately t = 1.1 where it reflects from the lower part of the right hand edge.After this point, the wave field becomes more complex, but the most dominant features are the pulse reflected back along the lower edge (0 < q < 1) for t > 1.1 and the pulse traveling back along the upper edge of the left half of the domain (3 < q < 3.5-see Figure 4) for t > 1.6.
Figure 10 shows the -transformed boundary solution Φ l .The left subplot (A) shows the absolute value at l = .The values of |Φ  | calculated using the BEM are compared to both of the proposed high frequency approximations.The main difference between the three methods can be observed along the left-upper (3 < q < 3.5) and lower (0 < q < 1) edges where the high frequency approximations give zero and the BEM does not, since both high frequency approximations of Φ l omit waves tangential to the boundary.In the zoomed panel, we highlight the main difference between the DEA and the SHFA results is apparent on the lower part of the right hand edge (1 < q < 1.5) where DEA includes reflected wave contributions but the SHFA does not.These results are consistent with those for the unit square, however there is also an optically shaded region for this problem for x 2 ⩾ 0.5 where the high frequency approximations are zero (between 1.5 < q < 3).This leads to an additional discrepancy between the BEM and the high frequency approximations that is particularly evident close to the edges of the shadow region (see the zoomed panel close to q = 1.5 and q = 3) where diffracted contributions are most dominant.The right subplot (B) shows a plot of Re(Φ l ) for wavenumbers close to Re(k) = 350.In this plot, Re(Φ l ) is computed using the DEA based approximation for Re(k) > 350, and using the BEM otherwise.Along the left-upper edge (3 < q < 3.5), the value of Re(Φ l ) jumps to zero when we switch from the BEM to the DEA approximation as would be expected from the left subplot (A).Along the left edge (3.5 < q < 4), the solutions from the two methods match up very well demonstrating the success of the phase reconstruction process.
Table 2 demonstrates the effect of lowering the wavenumber threshold k * with final time T = 0.5 up to which the exact solution is valid.As for the unit square, to retain a reasonable accuracy level it is necessary to over-resolve in the Laplace domain.We find again that reducing the wavenumber threshold to k * = 175 has only a minimal effect on the error but gives a significant decrease in computational time and decreasing further gives more moderate improvements Note: As the maximal BEM wavenumber is decreased, there is also a proportionate decrease in the number of boundary elements together with an increase in the number of Helmholtz problems solved overall.Note that the number of Helmholtz problems solved using the BEM remains approximately unchanged.in computational efficiency while causing more significant increases in the error.In this case, the phase reconstruction process does not fail until we reduce the wavenumber threshold to k * = 25 (compared to k * = 45 for the unit square), leading to significantly larger errors than for the SHFA.For the larger choices of wavenumber threshold considered, the performance of the DEA and SHFA based schemes was again very similar, but with the DEA approach requiring more computational resources as before.
Figure 11 shows the results of performing the same numerical experiments again but with the wave direction changed to Θ 0 = ∕4.In subplot (A), we notice that the interior solutions at x = (0.2, 0.25) are visually identical, even for the N = 2048 case that was visibly less accurate when  0 = 0, and consistent with the observations for the unit square.The boundary solution in subplot (B) shows the Gaussian pulse entering along the lower and left edges of the square (0 < q < 1 and 3.5 < q < 4) at the origin (q = 0 ≡ 4 (mod 4)) around t = t 0 = 0.1 and then traveling along those edges until reaching a vertex (q = 1 on the bottom-right and q = 3.5 at the top of the left-most edge) where a small part of the wave reflects, but mostly the propagation then continues along the adjoining edges.The wave field then becomes more complex and the main other observation is the localized peak at the vertex where q = 2.5 and at around t = 0.1 + 5 √ 2∕4 ≈ 1.87.This peak is the result of the reflected pulse from the right-hand edge that starts from the lower right vertex (q = 1) at around t = 0.1 + √ 2∕2 ≈ 0.81 and meets the vertices at q = 2 and q = 3 at the same time around t = 0.1 + √ 2 ≈ 1.5.One then sees the pulse move along the edges adjoining the vertex at q = 2.5 between t = 0.1 + √ 2 and t = 0.1 + 5 √ 2∕4 ≈ 1.87, when the pulses traveling along both of these edges meet and constructively interfere at q = 2.5.
The bottom row of subplots show the -transformed boundary solution Φ l and corresponds to Figure 10 for Θ 0 = 0. Owing to the change of incident wave direction, there is no optically shaded region in this case, but one would expect to see a diffracted contribution originating from the re-entrant corner (q = 3) that is not captured by the high frequency approximations.The left subplot (C) shows the absolute value at l =  as before.In this case, there are no edges where the wave direction is tangential to the boundary and hence the agreement between the BEM solution and the high frequency approximations is better than for Θ 0 = 0.The most significant divergence between three models occurs along the left-upper edge (3 < q < 3.5) and the lower part of the right edge (1 < q < 1.5) where the SHFA prediction is zero due to not accounting for and reflected wave contributions.The DEA result shows a smooth continuation of the solution across all edges whereas the BEM solution exhibits a small jump and oscillation both close to the bottom-right vertex (q = 1) and along the entire left-upper edge (3 < q < 3.5).The right subplot (D) shows a plot of Re(Φ l ) for wavenumbers close to Re(k) = 350.The solution is only shown for one of the dominant regions of Φ l near the lower-left vertex as for the unit square case.The solution is computed using the DEA based approximation for Re(k) > 350, and using the BEM otherwise.We observe that the solutions from the two methods match up very well, once again demonstrating the success of the phase reconstruction process.

Excitation by an interior point source
In this section, we consider the case when the wave problem is driven by a point source excitation of the form (3) and choose the temporal source profile to be a Gaussian pulse As before the parameters t 0 > 0 and  > 0 control the peak position and bandwidth of the Gaussian pulse, respectively.We perform numerical tests on a unit square domain and an irregular polygon as shown in Figure 12, which also depicts the location of the source point x 0 in each case.The fact that we are now approximating a circular wave as a superposition of plane waves in the DEA based approach means that using a small number of wave directions will no longer suffice and throughout this section we use L = 360 global directions in the DEA implementation.The boundary condition (6) related to the point source excitation gives rise to a DEA initial boundary density  0 of the form where ṽl is the -transform of the convolution integral (4), k l = is l ∕c, as before, and p 0 = sin( 0 )∕c.The angle  0 defines the direction of the outgoing plane wave (reflected specularly after arriving from x 0 ) at position q on Γ relative to the normal direction at q.The initial density (42) therefore represents the square amplitude of the -transformed Dirichlet boundary data |ṽ l | 2 transported along the direction corresponding to a specular reflection at q on Γ after arriving from x 0 .The exponential term represents the energy decay between leaving x 0 and arriving at q.The c cos ( 0 ) pre-factor is a consequence of projecting the density in Ω induced by the point source onto the boundary Γ-see Reference 31 for details.
Applying the CQ discretization to the convolution (4) provides an expression for |ṽ l | 2 as follows where gl 0 denotes the -transform of g 0 .Note that the approximation in the second line represents a high frequency asymptotic approximation of the Hankel function in the definition of G k l .

5.2.1
Unit square domain We first consider the case when Ω is a unit square with the source point located in the center at x 0 = (0.5, 0.5) and so we expect the resulting boundary solution to be the same on each edge by symmetry.As before we fix c = 1 for simplicity and choose  = 1∕64 in order to obtain a broadband signal.A choice of t 0 = 0.1 then gives that g 0 (0) is zero to double precision.Figure 13A shows the numerical interior solution at x = (0.2, 0.25), comparing the results of using the SHFA with the DEA based plane wave approximation up to t = 2.Most of the plotted lines represent the reverberant solution u, which is equivalent to the full solution  with the contribution of the direct source excitation corresponding to the free-space Green's function removed.The solid line shows a plot of  for completeness.As before, we apply a high frequency approximation whenever |Re(k l )| > 350 and employ M = 1024 boundary elements to provide a good level of accuracy up to the BEM cut-off frequency  −1 .We observe that both high frequency approximations produce identical looking results, but one observes small differences in the vicinity of the extrema for t > 1.3 when using N = 4096 time-steps compared to N = 2048.Figure 13B shows the solution  along the boundary computed using the DEA based high frequency approximation with the same parameter choices as before in the case N = 4096.The DEA and SHFA based approaches give visually identical solutions over the entire boundary and so just one of the plots is shown for brevity.The result in Figure 13B shows the expected symmetry, with identical solutions along all four edges of the square.One can also observe the expected physical behavior since the Gaussian pulse first arrives at the center of each edge (q = 0.5, 1.5, 2.5, 3.5) and then spreads across to the corners, after which time the peaks represent reflected waves propagating in a regular and symmetric fashion.
Figure 14 shows the -transformed boundary solution Φ l corresponding to the same parameter choices as for Figure 13.The left subplot (A) shows the absolute value at l =  calculated using the BEM as well as both of the proposed high frequency approximations.The result is shown only along the lower edge (0 < q < 1) to display the result more clearly instead of showing all four symmetrically repeated copies.In this case, the BEM solution has a noisy appearance owing to the presence of highly oscillatory interference patterns, particularly close to the vertices (also evident in the right subplots).In contrast to this, the SHFA result has a smooth slowly varying appearance corresponding to the fact that the SHFA does not include reflected wave contributions.The DEA result, on the other hand, does include the reflected contributions but loses accuracy for two interrelated reasons.First, the error resulting from the approximation of circular waves by a plane wave superposition and second, the challenge of phase reconstruction when the amplitudes approximated by the full wave BEM simulation are rapidly oscillating.The DEA result provides an underestimate of the oscillation amplitude near the corners and an overestimate either side the peak at the center of the edge.The issues present in both the DEA and SHFA approximations are also shown in the right-hand subplots (B, C), respectively, which show plots of Re(Φ l ) for wavenumbers close to Re(k) = 350.That is, one can observe a change in the plots when Re(k) > 350, which in subplot (B) for DEA appears as spurious noise most evident around the center of the edge at q = 0.5, and in subplot (C) for the SHFA appears as the omission of wave interference patterns most evident close to the corners at q = 0 and q = 1.Despite the limitations of the high frequency approximations observed in Figure 14, there did not appear to be a related noticeable issue in the time-dependent solution shown in Figure 13.To investigate this further, we repeated the above analysis with a lower threshold k * = 175 for switching to the high frequency approximations.We also amended the discretization parameters consistent with the results in Tables 1 and 2 in the previous section and use M = 512 boundary elements with N = Ñ∕4 = 4096.The results are shown in Figure 15 and we again observe the same limitations with the high frequency approximations in subplots (B-D), but now the difference in the result either side of the (lowered) threshold k * = 175 is even more noticeable in the lower subplots (C, D).Furthermore, the consequence of these accuracy limitations in the high frequency approximations is now clearly evident in the time-dependent interior solution shown in subplot (A).Note that we have only plotted the latter part of the time window 1.25 < t < 2, where the effect of the high frequency approximations can be observed as spurious small amplitude oscillations.Similar results in terms of the manifestation of inaccuracies from a physical optics high frequency approximation in the corresponding time-domain wave solution have also been observed for the exterior scattering models considered in Reference 18.

Irregular polygon domain
We now repeat the above analysis for the irregular polygon domain shown in Figure 12B, taking the source point to be x 0 = (−0.4,0.5).We choose c = 1 as before and the same parameters in (41) with  = 1∕64 and t 0 = 0.1.Figure 16A shows the numerical interior solution at x = (−0.75,0.5), comparing the results of using the SHFA with the DEA based plane wave approximation up to t = 2.As before, most of the plotted lines represent the reverberant solution u except the solid line, which shows a plot of the full solution  for completeness.We apply a high frequency approximation whenever |Re(k l )| > 350 with M = 1024 boundary elements.We again observe that both high frequency approximations produce identical looking results, as well as small differences in the vicinity of the extrema for the latter part of the time range when using N = 4096 time-steps compared to N = 2048.Figure 16B shows the solution  along the boundary computed using the DEA based high frequency approximation (here, identical in appearance to the SHFA result) with the same parameter choices as before in the case N = 4096.In comparison to the result for the unit square shown in Figure 13B, the plot here shows the expected loss of symmetry and lack of regularity corresponding to the irregular geometry.One can relate the position of the initial wavefront reaching Γ shown in Figure 16B to the geometry shown in Figure 12B, where the vertices appear as cusp points in the initial wavefront at q = 0, 0.96, 1.8, 2.6, and 3.1.
Figure 17 shows the -transformed boundary solution Φ l corresponding to the same parameter choices as for Figure 16.The upper subplot (A) shows the absolute value at l =  calculated using the BEM as well as both of the proposed high frequency approximations.As for the case of the unit square, the BEM solution has a noisy appearance owing to the presence of highly oscillatory interference patterns, particularly close to the vertices and the two shorter edges (q > 2.6).However, the SHFA result again has a smooth and slowly varying appearance corresponding to the fact that the SHFA does not include reflected wave contributions.The DEA result also corresponds to the previous observations for the unit square and includes some of the oscillations seen in the BEM result, but to a poor level of accuracy with the same root causes as before.The issues present in both the DEA and SHFA approximations are also shown in the lower subplots (B, C), respectively, which show plots of Re(Φ l ) for wavenumbers close to Re(k) = 350.In subplot (B) for DEA, we see the influence of the spurious noise in the DEA result, particularly for the shorter edges when q > 2.6.Interestingly, the result close to the longer edges is significantly better in terms of retaining its structure than the result for the unit square.This is unexpected because the DEA approximation on the square is direction preserving and all reflected ray directions are within the global direction basis set Θ l , l = 1, 2, … , 360.This is not true for the irregular polygon and there will be some inevitable directional approximation by mapping to the nearest direction in the basis set.In subplot (C) for SHFA, the omission of wave interference patterns can be observed, particularly in the vicinity of the vertex at q = 0.96, but this deficiency also appears to be less noticeable than for the unit square.
We have also repeated the study presented in Figures 16 and 17 with a lower threshold k * = 175 for switching to the high frequency approximations and make similar observations to the case of the unit square (the result is omitted for brevity).That is, the mismatch to the high frequency approximations becomes more severe and this, in turn, means that the inaccuracy also becomes evident in the time-domain solution plots.We remark that there is little to be gained from trying to improve the accuracy of the DEA result by increasing the number of directions in the global direction basis set since, in the present set up, the DEA calculations with L = 360 global directions used in this section only provide a modest time saving on using the BEM for all frequencies.This saving would be lost by adding more degrees of freedom to the model.

CONCLUSION
We have introduced two hybrid CQ based discretizations of the wave equation for interior acoustic Neumann problems with broadband boundary data or source terms.The CQ method was applied to reformulate the time-domain wave problem into a series of frequency domain Helmholtz problems with complex-valued wavenumbers, in which the boundary data and solutions are connected to those of the original time-dependent problem through the -transform.The idea behind the hybrid method is to use different methodologies for the approximation of these Helmholtz problems, dependent on the magnitude of the real part of the wavenumber, or equivalently, the frequency.For lower frequencies, we applied a conventional BEM for the discretization of the Helmholtz problems, while for more oscillatory problems, we developed two high frequency approximation methods, both of which are based on a plane wave decomposition of the acoustic field on the boundary.In the first approach, we apply DEA to numerically approximate the plane-wave amplitudes.The phases were reconstructed by matching the BEM solution to the plane wave ansatz at the maximal frequency where we apply the BEM and the minimal frequency where we apply the plane wave ansatz.We also proposed a SHFA method based on applying the Neumann-to-Dirichlet map for plane waves to the given boundary data.This approach is only valid if reflected wave contributions are negligible, which can often be the case for CQ discretizations if the imaginary part of the wavenumber is relatively large leading to high levels of dissipation.We then performed numerical experiments to demonstrate the effectiveness of both hybrid approaches for the cases of plane wave boundary data and excitation by an interior point source.In both situations, the hybrid methods were able to provide faster computations than using the BEM for all frequencies.However, the DEA based scheme for point sources suffered from poor accuracy owing to its inability to accurately represent circular waves without using a large number of ray directions that would remove any efficiency gains.A promising direction for further work in this case could be to use an image-source type method for the high frequency approximation, in which the solution would be expressed more efficiently as a series of circular waves.
Finally, we compare the DEA and SHFA high frequency approximations.In terms of computational efficiency and simplicity of implementation then the SHFA is highly preferable since it can be implemented at low cost within a single line of code.It also has the advantageous property that when refining Δt alone and fixing constant the maximal frequency at which a BEM calculation is performed, then the computational cost of the method overall then scales as (1), the same as in the related scheme proposed for exterior wave scattering in Reference 18.The DEA approach is able to provide better approximations that include reflected wave contributions.However, these improvements did not typically make any significant improvement to the time-domain result and when using lower threshold frequencies for switching from BEM to the high frequency approximation, DEA often performed worse owing to inaccuracies introduced in the phase reconstruction process.Given that the DEA approach also incurs a greater computational overhead, particularly when larger models including more discretized ray directions are needed, then we conclude that the accuracy gains from including higher order reflections in the hybrid CQ schemes here do not provide sufficient improvements to choose DEA over the SHFA based scheme.One significant reason for this appears to be the large imaginary parts of the wavenumbers within the BDF2 based CQ method that provide a strong damping effect within the DEA calculations, making the reflected wave/ray contributions relatively small.

respectively, where x m , m = 1 , 2 ,
… , M, are a set of collocation points located at the centers of the corresponding boundary elements E m .The basis functions b m , m = 1, 2, … , M, are piecewise constant functions defined by b

) 4
Problem set-up for the numerical experiments conducted in Section 5.1 showing the domains considered, the value of the boundary arclength q at each vertex, the propagation direction Θ 0 for the plane wave boundary data, and the interior evaluation point "•."Bold boundary sections indicate positions where the boundary data may be inhomogeneous

1 F I G U R E 5 6
Interior solution to the wave equation at x = (0.2, 0.25) inside a unit square with boundary data (37) and parameters Θ 0 = 0, t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 time-steps for T = 2 as specified in the legend, and  = 10 −8∕Ñ .High frequency approximations are applied whenever |Re(k l )| > 350.The exact solution is only displayed up to t = 1 since it is only valid prior to the first reflection.The zero solution until t = 0.2 is omitted from the plot Solution to the wave equation along the boundary of a unit square domain with boundary data (37) and parameters Θ 0 = 0, t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 = 4096 time-steps, and  = 10 −8∕Ñ .The high frequency approximations (A) DEA and (B) SHFA are applied when |Re(k l )| > 350 .

F I G U R E 7
The -transformed solution Φ l on the boundary of a unit square with boundary data (37) and parameters Θ 0 = 0, t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 = 4096, and  = 10 −8∕Ñ .The left subplot (A) shows |Φ l | for l = , which corresponds to the wavenumber k  where we stop using the BEM and instead use a high frequency approximation.In this plot, we compare the BEM result with the results from using both of the high frequency approximations.The right subplot (B) shows Re(Φ l ) computed using the DEA based hybrid CQ scheme on the upper and left edges of the square for a range of wavenumbers in the vicinity of l = , which approximately corresponds to Re(k) = 350 U R E 8 A repeat of the numerical results from Figures 5 to 7 to show the effect of changing the wave direction in the boundary condition (37) from Θ 0 = 0 to Θ 0 = ∕4, with all other parameters unchanged.The four subplots show: (A) The interior solution at x = (0.2, 0.25) (cf.Figure 5), (B) the solution along the boundary computed using the DEA based hybrid CQ scheme (cf.

Figure 6 )
, (C) the absolute value -transformed boundary solution Φ l at l =  (cf.Figure7A), and (D) the real part of the -transformed boundary solution Φ l for a range of wavenumbers in the vicinity of Re(k) = 350 computed using the DEA based hybrid CQ scheme (cf.Figure7B).
Solution to the wave equation for an L-shaped domain with boundary data (37) and parameters Θ 0 = 0, t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, T = 2, and  = 10 −8∕Ñ .The left-hand plot (A) shows the interior solution at x = (0.2, 0.25) for N = Ñ∕2 time-steps for as specified in the legend.The exact solution is only displayed up to t = 0.5 since it is only valid prior to the pulse arriving at the re-entrant corner.The zero solution until t = 0.2 is omitted from the plot.The right-hand plot (B) shows the solution along the boundary of the L-shaped domain for N = Ñ∕2 = 4096 using the DEA based hybrid CQ scheme.In both cases, high frequency approximations are applied whenever |Re(k l )| > 350 U R E 10The -transformed solution Φ l on the boundary of an L-shaped domain with boundary data (37) and parameters Θ 0 = 0, t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 = 4096, and  = 10 −8∕Ñ .The left subplot (A) shows |Φ l | for l = , which corresponds to the wavenumber k  where we stop using the BEM and instead use a high frequency approximation.In this plot, we compare the BEM result with the results from using both of the high frequency approximations.The right subplot (B) shows Re(Φ l ) computed using the DEA based hybrid CQ scheme on the upper and left edges of the left part of the L-shape (see Figure4B) for a range of wavenumbers in the vicinity of l = , which approximately corresponds to Re(k) = 350 U R E 11 A repeat of the numerical results from Figures 9 and 10 to show the effect of changing the wave direction in the boundary condition (37) from Θ 0 = 0 to Θ 0 = ∕4, with all other parameters unchanged.The four subplots show: (A) The interior solution at x = (0.2, 0.25) (cf. Figure 9A), (B) the solution along the boundary computed using the DEA based hybrid CQ scheme (cf. Figure 9B), (C) the absolute value -transformed boundary solution Φ l at l =  (cf. Figure 10A), and (D) the real part of the -transformed boundary solution Φ l for a range of wavenumbers in the vicinity of Re(k) = 350 computed using the DEA based hybrid CQ scheme (cf.Figure 10B)

12 13
Problem set-up for the numerical experiments conducted in Section 5.2 showing the domains considered, the value of the boundary arclength q at each vertex (to 2 digits if not integer valued), the source point x 0 and the interior evaluation point "•" Solution to the wave equation for a unit square domain driven by a point source located at the center x 0 = (0.5, 0.5) with temporal profile (41) and parameters t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, T = 2, and  = 10 −8∕Ñ .The left-hand plot (A) shows either the interior solution  or the reverberant solution u (see Section 2) at x = (0.2, 0.5) for N = Ñ∕2 time-steps for as specified in the legend.The zero solution until t = 0.2 is omitted from the plot.The right-hand plot (B) shows the solution along the boundary of the unit square domain for N = Ñ∕2 = 4096 using the DEA based hybrid CQ scheme.In both cases, high frequency approximations are applied whenever |Re(k l )| > 350 The -transformed solution Φ l on the boundary of a unit square driven by a point source located at the center x 0 = (0.5, 0.5) with temporal profile (41) and parameters t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 = 4096, and  = 10 −8∕Ñ .All subplots show results along the lower edge of the square only, since the results on all four edges are identical by symmetry.The left subplot (A) shows |Φ l | for l = , which corresponds to the wavenumber k  where we stop using the BEM and instead use a high frequency approximation.The right subplots (B, C) show Re(Φ l ) computed using the DEA and SHFA based hybrid CQ schemes, respectively.The results are shown for a range of wavenumbers in the vicinity of l = , which approximately corresponds to Re(k) = 350 U R E 15 (A) Solution to the wave equation for a unit square domain driven by a point source located at the center x 0 = (0.5, 0.5) with temporal profile (41) and parameters t 0 = 0.1, c = 1,  = 1∕64.Plot (A) shows either the interior solution  or the reverberant solution u (see Section 2) at x = (0.2, 0.5).The solution is shown for 1.25 < t < T = 2 to focus on the late time behavior.(B-D) The -transformed solution Φ l on the boundary of the unit square driven by the point source excitation used for plot (A).Subplots (B-D) show results along the lower edge of the square only, since the results on all four edges are identical by symmetry.Subplot (B) shows |Φ l | for l = , which corresponds to the wavenumber k  where we stop using the BEM and instead use a high frequency approximation.The lower subplots (C, D) show Re(Φ l ) computed using the DEA and SHFA based hybrid CQ schemes, respectively.The results are shown for a range of wavenumbers in the vicinity of l = , which approximately corresponds to Re(k) = 175.Discretization parameters: M = 512 boundary elements, N = Ñ∕4 = 4096, and  = 10 −8∕Ñ

16
Solution to the wave equation for an irregular polygon domain driven by a point source located at x 0 = (−0.4,0.5) with temporal profile (41) and parameters t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, T = 2, and  = 10 −8∕Ñ .The left-hand plot (A) shows either the interior solution  or the reverberant solution u (see Section 2) at x = (−0.75,0.5) for N = Ñ∕2 time-steps for as specified in the legend.The zero solution until t = 0.2 is omitted from the plot.The right-hand plot (B) shows the solution along the boundary of the unit square domain for N = Ñ∕2 = 4096 using the DEA based hybrid CQ scheme.In both cases, high frequency approximations are applied whenever |Re(k l )| > 350 U R E 17 The -transformed solution Φ l on the boundary of an irregular polygon driven by a point source located at x 0 = (−0.4,0.5) with temporal profile (41) and parameters t 0 = 0.1, c = 1,  = 1∕64.The discretization parameters are chosen as M = 1024 boundary elements, N = Ñ∕2 = 4096, and  = 10 −8∕Ñ .The upper subplot (A) shows |Φ l | for l = , which corresponds to the wavenumber k  where we stop using the BEM and instead use a high frequency approximation.The lower subplots (B, C) show Re(Φ l ) computed using the DEA and SHFA based hybrid CQ schemes, respectively.The results are shown for a range of wavenumbers in the vicinity of l = , which approximately corresponds to Re(k) = 350 Errors and computation times for the DEA and SHFA based hybrid CQ schemes on the unit square domain for T = 1 when lowering the choice of wavenumber threshold k * TA B L E 1 Errors and computation times for the DEA and SHFA based hybrid CQ schemes on the L-shaped domain for T = 0.5 when lowering the choice of wavenumber threshold k * TA B L E 2