The Burton and Miller method: Unlocking another mystery of its coupling parameter

The phenomenon of irregular frequencies or spurious modes when solving the Kirchhoff–Helmholtz integral equation has been extensively studied over the last six or seven decades. A class of common methods to overcome this phenomenon uses the linear combination of the Kirchhoff–Helmholtz integral equation and its normal derivative. When solving the Neumann problem, this method is usually referred to as the Burton and Miller method. This method uses a coupling parameter which, theoretically, should be complex with nonvanishing imaginary part. In practice, it is usually chosen proportional or even equal to i/k. A literature review of papers about the Burton and Miller method and its implementations revealed that, in some cases, it is better to use -i/k as coupling parameter. The better choice depends on the specific formulation, in particular, on the harmonic time dependence and on the fundamental solution or Green's function, respectively. Surprisingly, an unexpectedly large number of studies is based on the wrong choice of the sign in the coupling parameter. Herein, it is described which sign of the coupling parameter should be used for different configurations. Furthermore, it will be shown that the wrong sign does not just make the solution process inefficient but can lead to completely wrong results in some cases.


Introduction
The Burton and Miller formulation 1 for exterior acoustic problems is well-known since it is free of fictitious resonances, see also Ref. 2. The only parameter of the Burton and Miller formulation is known as the coupling parameter. When reading the literature applications of the method, it seems clear to choose this parameter to be i/k,  or at least positively proportional [30][31][32] or asymptotically proportional 33,34 to this value. (While i is the imaginary unit, the wave-number k = ω/c is the quotient of the angular frequency ω = 2πf with f denoting frequency and the speed of sound c.) Some authors did not explicitly mention their choice of coupling parameter. 35,36 Interestingly, some authors apply i/k as a negative value. 37, 38 Terai 39 has even clarified that the positive value is valid for a harmonic time dependence of e +iωt whereas in case of using e −iωt , the coupling parameter should be negative. Especially Terai's work is essential in this context but seems to have been ignored by most authors who have used e −iωt . However, Kress 14 and Amini 3 were both using this kind of harmonic time dependence and clearly found that i/k is a very good and even close to optimal choice to minimize the condition number for a sphere at high frequencies. Apparently, it is not completely clear, what the optimal coupling parameter for the Burton and Miller formulation is. Furthermore, it will even be shown that many authors are using a coupling parameter which is not optimal and which can lead to wrong results.

BEM Formulation of the Burton and Miller Method
Derivation of the wave equation, discussion of boundary conditions, weak formulation and discretization process are presented in a reduced way. A more detailed presentation is found in the concept chapter. 40

Helmholtz equation and boundary conditions
Herein, linear acoustic problems are defined in the infinite domain Ω with the complement finite domain of Ω c and Γ is representing the closed boundary between Ω and Ω c . The normal vector is pointing into the complementary domain Ω c . The wave equation is valid for the sound pressurep. Alternatively, a velocity potential may be used. To complete a solution, the differential equation requires boundary conditions and initial conditions, which will be specified when used. For time-harmonic problems, a time dependence is introduced asp with α = ±1. Applying the time-harmonic dependence of p to Eq. (1) leads to the Helmholtz equation for the sound pressure as This result is independent of α. Neumann boundary conditions are assumed. For them, the normal particle velocities of the fluid v f equal the (prescribed) particle velocity of the underlying radiator v s as Hence, the normal fluid particle velocity v f is related to the normal derivative of the sound pressure p by means of the Euler equation in frequency domain. The constant s is given by where 0 is the ambient density of the fluid. The vector n( x) represents the outward normal at the surface point x and ∂/∂n( x) is the normal derivative.

Weak formulation
A weak formulation is based on introducing the weight function χ( x) and "testing" it with the Helmholtz operator such that Integrating by parts twice gives The second part of Eq. (7) consists of two boundary integrals and one domain integral. This domain integral can be transformed into an integral-free term by the using fundamental solution G( x, y) in the sense of distributions. G is known as free-space Green's function. In terms of physics, G( x, y) can be understood as the sound pressure distribution according to a point source (monopole) in y. Together with the harmonic time-dependence of e −iωt , it represents an outgoing wave. In the literature, the author has found this fundamental solution for the Helmholtz equation as with r as the Euclidean distance between field point x and source point y as r( x, y) = | x− y|. Function G should represent the solution of the equation but it does not. (The function δ( x, y) is the Dirac or delta function with the origin at y.) where, similar to α in Eq. (2), β = α. However, there are many examples in literature where β = −α is chosen, e.g. see Refs. 5,10,11,15,16,[18][19][20][21][23][24][25][26][27][28][29] Usually, this choice is not obvious at this point but becomes clear with the formulation of the integral equation and the integral-free term. Therefore, this derivation is continued with the assumption that α = ±1 and β = ±1 are independent of each other. Then, the correct Green's function G should be written as However, as mentioned above, the author has found only the choice of G as given in Eq. (8) in the literature in the context of the Burton and Miller method. For solution of Eq. (10), the sign of the exponent of e −iαkr( x, y) is arbitrary and, thus, independent of α. It becomes determined when requiring the Green's function to fulfill the Sommerfeld radiation condition which depends on α as well.
Applying the property of the fundamental solution and the delta function of Eq. (8) results in With this result and application of the boundary condition (4), Eq. (7) is rewritten as Equation (13) is known as representation formula. For y ∈ Γ, it is known as the Kirchhoff-Helmholtz (boundary) integral equation. It is also referred to as the first integral equation. The term βζ( y)p( y) will be called the integral-free term for this first integral equation. The value of ζ (with 0 < ζ < 1) is determined by the surface smoothness. At a smooth surface point, ζ = 0.5. Note that plus and minus signs of either the first term or the second and the third term may be different if the direction of the normal vector is chosen in opposite direction.

Approximation and discretization by collocation
Independent of the discretization method, approximations of the physical quantities are formulated. First of all, the sound pressure p( x) is approximated as where p l represents the discrete sound pressure at point x l and φ l is the lth basis function for approximation. Further, it is assumed that a similar approximation is formulated for the particle velocity of the underlying structure If v s is explicitly known, these approximations are not necessary for evaluation of the boundary integrals in Eq. (13). However, there are many practical cases where the structural particle velocity is the result of a finite element simulation and available only as piecewise defined function. The number of basis functions φ l andφ j is given by N andN , respectively. If the particle velocity of the structure is a known function, N accounts for the degree of freedom. It is common thatN = N .
The collocation method requires testing Eq. (13) with the Dirac function δ( y, z). This integration is known analytically, cf. Eq. (12). It yields which is basically the same expression as shown in Eq. (13). The major difference between Eqs. (13) and (16) is that the former is actually a continuous integral equation whereas the latter is valid just for the discrete point z. This means that the integral equation is fulfilled at a number of discrete points, i.e. the collocation points z l . It is common practice that the collocation points coincide with the nodes of the piecewise formulated approximation of the sound pressure as shown in Eq. (14). For further considerations it is assumed that φ l ( z k ) = δ lk where δ lk is the Kronecker symbol with δ lk = 0 for l = k and δ lk = 1 for l = k. Then, applying the approximation of Eqs. (14) and (15) yields the matrix equation as Matrix G is the system matrix of the single layer potential as and matrix H contains the integral-free term and the contribution of the double layer potential as For the Neumann problem, Eq. (17) is solved for the sound pressure at the nodes p.
It is a well-studied fact that the system matrix H becomes ill-conditioned for frequencies close to the eigenvalues of the Dirichlet problem of the complementary domain, i.e. the eigenvalues of G. To overcome this problem on behalf of the Burton and Miller formulation, the normal derivative of Eq. (13) is used.

Normal derivative integral equation
The normal derivative of Eq. (13) is given as A similar discretization process as in the previous subsection yields the matrix equation as with Matrix F as the system matrix of the integral-free term and the adjoint double layer potential as and matrix E with the contribution of the hypersingular operator as For the Neumann problem, Eq. (21) is solved for the sound pressure at the nodes p. However, this is quite uncommon since the eigenvalues of matrix E are the same as for the interior Neumann problem. This means that the lowest eigenvalue and, with it, the lowest fictitious resonance occurs at the frequency of 0 Hz which makes the single use of this equation rather impractical.

The Burton and Miller method
The Burton and Miller method is based on a linear combination of Eqs. (13) and (20) for which the coupling parameter η is introduced. The coupled equation reads as follows The matrix form is easily derived from this equation and equations (17) and (21) as which is, again for the Neumann problem, solved for p. Many authors have discussed the choice of the coupling parameter η. [1][2][3]14,39 It is clear that the requirement for η is such that it should be complex with nonzero imaginary part to yield a unique solution for Eq. (24). This regards the continuous integral equation. The requirements are stronger for the discretized integral Eq. (25). Usually, the coupling parameter is chosen purely imaginary. Amini 3 found that the parameter i/k almost minimizes the condition number for a sphere at high frequencies. Terai 39 had provided a physically motivated explanation for this choice 10 years before.

Rigorous mathematical studies of the Dirichlet problem
Mathematical investigations in recent years focused at a similar problem, i.e. the Dirichlet problem of the Kirchhoff-Helmholtz integral equation for which a similar linear combination of the integral equation and its normal derivative is used. This method has been well-known since the first papers were published approximately 50 years ago 41,42 and these papers were known and mentioned by Burton and Miller. 1 Recent investigations confirmed for a large group of domains, so-called nontrapping domains, that choosing |η| proportional to 1/k minimizes the condition number of the system of equations for which the sign of η is not relevant. 43-47 According to Spence, 43 the only property of the system matrix which strongly depends on the sign of η is coercivity. The system matrix for the Dirichlet problem (discretized by a Galerkin method) is coercive for η = i/k but not for η = −i/k, 48 of course assuming αβ = 1. It is worth mentioning that the mathematical studies referenced here are all consistent insofar, that they use αβ = 1.

Test Problems
It has been explained in the previous sections that the choice of a positive or negative coupling parameter η depends on the formulation of the integral equations, in particular on the choice of the parameters α and β. In what follows, the correct choice of η will be referred as η = i/k, i.e. the correct choice for αβ = 1.

Spherical scatterer
The first example comprises a spherical scatterer in a field of plane waves. For the Neumann problem, the sphere is assumed to be rigid. Material data of air are used such that the speed of sound c = 340 m/s and the density ρ 0 = 1.3 kg/m 3 . The analytical solution for the total sound pressure p is well-known as sum of incident and scattered sound pressures, p i and p s , respectively, see for example Ref. 49 In Eq. (26), p 0 represents the sound pressure amplitude of the incident wave, j n and h n are the spherical Bessel and Hankel functions of the first kind, respectively. P n denotes the Legendre polynomials and R = 1 m is the radius of the spherical scatterer. Since the problem is axisymmetric, the two parameters r and ϑ allow a complete evaluation in space, r is the distance from the center of the sphere and ϑ is the angle such that the shadow zone is located at ϑ = 0 and the illuminated zone at ϑ = π. A sound pressure amplitude for the incident wave of p 0 = 1 Pa is chosen. The numerical model is a model which uses super-parametric boundary elements for which the geometry is approximated by quadratic quadrilateral elements (nine nodes per element) and the sound pressure is approximated by linear discontinuous boundary elements. 17 is approximately 0.1 m, i.e. 64 elements along the diameter of the sphere. Selecting a maximum frequency of 1700 Hz for the analysis, i.e. kR = 10π, results in a mesh for which 3.4 elements per wavelength are counted. According to a former study of the author, 50 this should be sufficient. Figures 1 and 2 show the sound pressure magnitude at r = 2R and ϑ = 0 for different configurations. The comparison in Fig. 1 confirms that a method to suppress the irregular modes is really necessary. In Fig. 2, the analytical solution at the same position is compared with the numeric solution of the Burton and Miller formulation using the coupling parameters i/k and −i/k, respectively. It can be recognized very clearly that the choice of η = −i/k is significantly changing the result. Above 1 kHz, the error of the scattered pressure p s reaches more than 25%, whereas it remains much smaller if i/k is chosen. Interestingly, a clear difference between the analytical and the i/k-solution can be recognized in the lower frequency range around 500 Hz. This is somewhat unexpected and remains a (marginal) open question of this study. Figure 3 presents the sound pressure magnitude at r = 2R and ϑ = π for the analytical solution and for the Burton and Miller formulation using the coupling parameters i/k and −i/k, respectively. For this position, the numeric solutions follow the analytical solution rather well. It becomes obvious that selection of the wrong coupling parameter is again spoiling the result. A difference between the numeric solution for η = i/k and the analytical solution is hardly recognizable except for the highest frequencies.
A practical aspect of numeric solutions consists in the computation time. The setup of the system of equations and the solution of the system of equations account for the main contributors for boundary element solutions. While the setup time is independent of the coupling parameter and in the current implementation even quite similar for solution of the Kirchhoff-Helmholtz integral equation alone and the Burton and Miller formulation, the convergence of the iterative solution of the system of equations is investigated. A GMRes algorithm 51,52 without preconditioning is used for this. The solution is assumed to have converged after the residual has reached a certain prescribed value. Herein the residual is evaluated as where A and b are the system matrix and the right hand side, respectively. The current iteration is denoted by x n . For all computations in this article, the iteration is terminated after R n ≤ 10 −10 .
The results are depicted in Fig. 4. Obviously, the simple Kirchhoff-Helmholtz integral equation which requires inversion of H only is very efficiently solved in the low frequency range. However, matrix H becomes the worse conditioned the higher the frequency is. Among other reasons, the increasing density of irregular frequencies leads to a rising condition number which leads to more iterations and, thus, more matrix vector products for

Cat's eye radiator
The second example is already known from a detailed previous study. 17 There, the cat's eye radiator was chosen to investigate the occurrence of irregular frequencies for a more complex structure than the sphere. The cat's eye model is a sphere with the positive octant cut out. In the radiation problem, it is assumed that the vibrating surface coincides with the spherical part of the surface. The plane surface areas of the missing octant remain calm. This implies that the noise transfer function, i.e. sound pressure divided by surface velocity in terms of frequency, returns a very smooth function similar to the solution of the radiating sphere where it is easy to identify frequencies where the solution fails. Consequently, in Ref. 17 it could be shown that some methods which should suppress the irregular frequencies are still failing if higher frequencies are investigated. The Burton and Miller method, even in a reduced form where the hypersingular integral equation was only applied to collocation at the center point of a continuous quadratic element, supplied smooth results and was believed to be a reliable method to provide a solution free of fictitious modes. Similar to the spherical scatterer, material data of air as density ρ 0 = 1.3 kg/m 3 and speed of sound c = 340 m/s are assumed. The spherical radius is again taken to be R = 1.0 m and the particle velocity for the vibrating surface is v = 1 m/s. Again, linear discontinuous elements are used. Applying a very similar mesh as for the sphere in the previous example results in a mesh of 1920 quadrilateral super-parametric elements. Again, the sound pressure is approximated by linear discontinuous boundary elements. Hence, the same frequency range as for the spherical scatterer is analyzed. An analytical solution is not known for this particular problem of the cat's eye. However, Mechel 53 discussed analytical series solutions of the cat's eye model in particular for scattering.
For the current investigation, it is assumed that the positive octant has been cut out of a sphere, i.e. the octant which encompasses positive coordinates in all three directions. Hence, the point is located above the back side of the cat's eye. Furthermore, the points P 2 = 0.01R[10, 10, 1] and P 3 = 0.05R [10,10,1] are defined within the cat's eye, i.e. within the missing octant. Figure 5 depicts the sound pressure at P 1 . The result is very similar to the figure provided in the former paper 17 and the curves for both coupling parameters are very similar to each other. Both curves are smooth and do not give any indication of an irregular frequency. The sound pressure curves for P 2 and P 3 are presented in Fig. 6. While the general behavior of two curves for the same position is the same, there are at least some obvious deviations between the curves. This time, it is not possible to decide which of the curves is more accurate. However, the author assumes that the curves for η = i/k provide the more accurate solution here.
Analogously to Fig. 4, Fig. 7 shows the number of iterations required for the GMRes to converge. The picture looks similar to what was experienced for the sphere. Just worse for the wrong coupling parameter and the pure Kirchhoff-Helmholtz integral equation. The number of iterations is, again, decreasing for the coupling parameter i/k and, for high frequencies, it is surprisingly (< 50) low whereas for the wrong coupling parameter, the number of iterations reaches 900 for some frequencies.

The Radiatterer
The Radiatterer is a recently created benchmark problem of the European Acoustics Association. 54 This model is the result of searching for a new shape with many edges and corners while still being able to easily create hierarchical meshes. It contains a number of resonators  Again, material data of air as density ρ 0 = 1.3 kg/m 3 and speed of sound c = 340 m/s are assumed. For the current simple test case, a unit particle velocity is applied to all surfaces. This unit particle velocity is set to be 1 mm/s. While the sound pressure has been evaluated at the surface and at numerous other positions in the near field and in the far field of the Radiatterer, in this study, only the radiated sound power P 55,56 where it is unusual to use the magnitude |P |. However, this is required here since in some places when using the wrong coupling parameter η = −i/k, a negative radiated sound power is evaluated. Obviously, this is physically insufficient and must be wrong. These values of negative sound power occur mainly at some resonances and they do not appear in case of η = i/k.
While the results in the previous examples were always similar for both coupling parameters, a substantial error is generated by the wrong coupling parameter when a resonator is analyzed. Interestingly, this effect is not observed for all resonances but only a few. Moreover, the negative sound power at resonances is mainly associated to standing waves within the Helmholtz resonator shown in the lower right corner of the lower right subfigure of Fig. 8. At these resonances, the main contribution to the radiated sound power stems from inside the resonator. Therefore, it is worth to check the sound pressure distribution inside the resonator. Only the real part of the sound pressure is relevant since a real valued particle velocity is used. Figure 10 displays the real part of the sound pressure at the surface of the Radiatterer for the two different values of αβ for the frequency of 262 Hz. While, from the outside, the Radiatterer's surface appears green, i.e. the real component of the sound pressure is small compared to the same quantity at other locations, the solution becomes colorful inside the resonator. Although looking quite similar at first glance, the maxima and minima are inverted when changing the value of αβ from −1 to 1.   low frequency range but diverge from each other rather quickly. Different from the results in the previous examples, the number of required iterations is generally increasing with frequency even for the correct coupling parameter. Above approximately 650 Hz, solution of the Kirchhoff-Helmholtz integral equation takes more iterations than the Burton and Miller formulation with the coupling parameter i/k. Finally, the Radiatterer is investigated for a finer tuning of the coupling parameter. The motivation behind this test consists in the fact that some authors are using a coupling parameter of η = Ci/k with a real value C fulfilling the condition 0 < C < 1. 30-32 Figure 12 reveals that, at least for the Radiatterer example, values of C < 1 can be advantageous in the low frequency range. In the higher frequency range, the curve for C = 1/ √ 10 approaches the C = 1-curve whereas the C = 1/10-curve is clearly above the other two. Hence, an optimal choice of the coupling parameter depends on the frequency and, most likely, on the problem. Therefore, the search for an optimal coupling parameter appears to be very complex and would be beyond the scope of this paper.

Conclusion
The main motivation for this investigation has been the apparent unawareness of many researchers, mainly engineers, that they are using the wrong coupling parameter for the Burton and Miller formulation. Scanning the literature showed that a coupling parameter with the correct sign was chosen in Refs. 3, 4, 6, 7, 9, 12-14, 17, 22, 30-34, 37, 39 while a coupling parameter with the wrong sign was chosen in Refs. 5,10,11,15,16,[18][19][20][21][23][24][25][26][27][28][29]. For some papers, it was impossible to decide whether the authors chose the better solution or not. 8,35,36,38,57 In these cases, some information for the decision has been missing in the paper. Although being the linear combination of two zeros, cf. Eq. (24), the sign of the coupling parameter is not just influencing the efficiency of the iterative solution of the system of equations. It also manipulates the numeric results and can end up in solutions which are completely wrong and unphysical as shown. To avoid this and to decide whether the parameter is correctly chosen, it is recommended that α in the harmonic time dependence and β as the factor in the integral-free term are always supplied.