# Helmholtz equation and non-singular boundary elements applied to multi-disciplinary physical problems

The famous scientist Hermann von Helmholtz was born 200 years ago. Many complex physical wave phenomena in engineering can effectively be described using one or a set of equations named after him: the Helmholtz equation. Although this has been known for a long time from a theoretical point of view, the actual numerical implementation has often been hindered by divergence free and/or curl free constraints. There is further a need for a numerical method which is accurate, reliable and takes into account radiation conditions at infinity. The classical boundary element method (BEM) satisfies the last condition, yet one has to deal with singularities in the implementation. Since these singularities are mathematical in origin, they can actually be removed without losing accuracy by subtracting a carefully chosen theoretical solution with the same singular behavior. We review here how a recently developed singularity-free 3D boundary element framework with superior accuracy can be used to tackle such problems only using one or more Helmholtz equations with higher order (quadratic) elements which can tackle complex shapes. Examples are given for acoustics (a Helmholtz resonator among others) and electromagnetic scattering. We briefly touch on the Helmholtz decomposition for dynamic elastic waves as well.

## Authors

• 2 publications
• 32 publications
04/27/2021

### Including monopoles to a fully desingularized boundary element method for acoustics

The inclusion of domain (point) sources into a three dimensional boundar...
04/20/2022

### Isogeometric boundary element method for acoustic scattering by a submarine

Isogeometric analysis with the boundary element method (IGABEM) has rece...
02/02/2021

### The conforming virtual element method for polyharmonic and elastodynamics problems: a review

In this paper, we review recent results on the conforming virtual elemen...
07/10/2018

### A Numerical Comparison of an Isogeometric and a Classical Higher-Order Approach to the Electric Field Integral Equation

In this paper, we advocate a novel spline-based isogeometric approach fo...
03/16/2021

### Isogemetric Analysis and Symmetric Galerkin BEM: a 2D numerical study

Isogeometric approach applied to Boundary Element Methods is an emerging...
09/17/2020

### The Boundary Element Method of Peridynamics

The peridynamic theory reformulates the governing equation of continuum ...
12/20/2019

### Fast hybrid numerical-asymptotic boundary element methods for high frequency screen and aperture problems based on least-squares collocation

We present a hybrid numerical-asymptotic (HNA) boundary element method (...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Hermann von Helmholtz (1821-1894) [8, 7] was born 200 years ago on 31 August 1821. The scalar Helmholtz equation,

, which inherits its name from the famous German scientist, appears in many fields of science and engineering involving waves. The waves thus described can be longitudinal (as in soundwaves), transverse (as in electromagnetic waves) or both (as in dynamic elastic waves in solid materials). This appears to us as the perfect moment to revisit the Helmholtz equation and see how it can be solved most efficiently for real three dimensional engineering problems.

The main goal of this article is to demonstrate that many classical engineering problems such as in acoustics, electromagnetics and elasticity can, in principle, all be formulated using one or a set of scalar Helmholtz equations. Thus a robust and efficient method to solve the Helmholtz equation is most desirable. Having developed a boundary element method for the scalar Helmholtz equation which is free of singularities, it is now possible to tackle more problems in classical applied physics essentially using the same numerical framework. This method was developed after the realization that, ideally, if there are no singularities in the physical system, no such singularities should have to appear in the mathematical equivalent description of this physical system. We will show examples mainly for (decaying) scattering waves from an object situated in an infinite medium and focus on sound and electromagnetic waves. Especially for the field of electromagnetics, the framework presented here is substantially different from the industry standard implementation.

The structure of this paper proceeds as follows. After this short introduction, different classic physical problems that can be formulated based on the scalar Helmholtz equation are demonstrated in Sec. 2

, with as examples, among others, the Helmholtz resonator in acoustics. After a brief discussion on curl free vector fields in Sec.

3, we discuss divergence free vector fields in Sec. 4, with particular attention to electromagnetic scattering problems. In Sec. 6, the recently developed non-singular boundary integral method for solving the Helmholtz equation is derived in detail in which all the singularities in the integrands and the terms associated with the solid angle are eliminated analytically. Then, the discussion and conclusions are given in Secs. 7 and 8, respectively.

## 2 Scalar Helmholtz equation

### 2.1 Acoustics

The wave equation describing the propagation of a quantity is in the time domain , with the wave speed and time. Assuming a harmonic time dependence , with angular frequency and unit imaginary part ‘’, the wave equation transforms into the well-known Helmholtz equation as:

 ∇2ϕ+k2ϕ=0 (1)

with . This is the equation as it appeared in Helmholtz’ book [8] on “acoustics in pipes with open ends” while studying organ pipes as Eq.3b on page 18/19. The Helmholtz equation is often used to describe acoustic waves, in which represents the velocity potential or the pressure perturbation in the medium. Normally, is a real valued number, often referred to as wave number in wave physics, while is a complex quantity. The boundary conditions are usually that on a surface bounding the domain is given (Dirichlet condition), or its normal derivative on a surface is given (Neumann condition) with the normal vector pointing out of the domain, or a combination of those two (Robin conditions). Sound-waves are typical examples of so-called longitudinal waves.

The Helmholtz equation is an elliptic equation. This means that the properties on the surface(s) of the domain determine entirely the solution anywhere in the domain. This allows us to write the following surface integral equation which is equivalent to Eq. (1):

 c(x0)ϕ(x0)+∫Sϕ(x)∂Gk(x,x0)∂n d% S(x)=∫S∂ϕ(x)∂nGk(x,x0) dS(x) (2)

where is the Green’s function, is a source point on the surface, is the computation point on the surface , and is the solid angle at . At this moment, when assuming the surface is divided into elements with a total of nodes, it is enough to note that all potentials and their normal derivatives are related by the following matrix system:

 H⋅ϕ––=G⋅∂ϕ/∂n–––––––. (3)

The matrices and are the numerical matrix equivalents of the boundary element integrals in Eq. (2), and and are column vectors of length in which each component and belongs to and , respectively, of the ’th node located at with .

Note that Eq. (2) is the conventional version of the boundary element method (BEM) for the Helmholtz equation, in which the difficulties to deal with singular integrals (and solid angle calculations) have hindered the implementation of this method considerably. Nevertheless, as detailed in Sec. 6, we have developed an advanced non-singular boundary element method in which the solid angle and singularities in the integrals are removed analytically. This improvement has made the boundary element method much easier to use, especially for higher order quadratic elements. There are, by the way, other methods to solve the Helmholtz equation, such as with finite elements.

To demonstrate some interesting wave phenomena in acoustics, in Fig. 1, we simulated an acoustic transducer (essentially a parabolic disk) oscillating up and down to generate sound waves which are then reflected onto a similar but stationary rigid disk that is rotated at an angle. Clearly, the focal points of the transducers can be observed. In Fig. 1(a), the absolute pressure is shown, while in Fig. 1(b), the instantaneous pressure profile is shown When simulating this problem, only a surface mesh on the two disk boundaries is needed, and the pressure in the domain is obtained by postprocessing where the complex interference pattern is clearly visible.

### 2.2 The Helmholtz resonator

In this article, dedicated to Helmholtz, it seems appropriate to simulate an object that was named after him as well: a Helmholtz resonator. It consists of a thin spherical shell with an opening with a neck on top. The Helmholtz resonator is assumed to be acoustically hard: on its surface (thus ), with the superscripts , and indicating the total, incoming and scattered wave, respectively. In Fig. 2 the Helmholtz cavity as used in our simulations is shown, both in full 3D view and in cross sectional view. The spherical part has an outer radius of , an inner radius of , a neck with length and the neck tapers down from a value of to at the top. This will give us a theoretical resonance which should occur at with the opening area of the neck, the inner volume of the sphere [20].

The incident acoustic wave is coming from the left and travels to the right in Fig. 3. We start with a very low value in Fig. 3(a); the pressure inside the spherical part of the cavity is rather uniform, (slightly higher than the reference pressure of ). In Fig. 3(b) at near the Helmholtz cavity resonance frequency, the pressure amplitude inside the resonator reaches a value of (i.e. 20 times the reference pressure). In Fig. 3(c) we increase to a value of . The pressure is now and is considerably lower than the reference pressure. The first internal resonance frequency of the sphere is shown in Fig. 3(d) at . The maximum pressure is now situated in the neck and reaches a value of about . A further resonance is shown in Fig. 3(e) for , near the second resonance frequency of the sphere. The maximum pressure amplitude reaches with multiple maximum values inside the sphere. Finally, in Fig. 3(f) yet another resonance is shown, now at , with two horizontally placed pressure peaks of . Also note the interference pattern to the left of the Helmholtz cavity for higher values, which is caused by the interaction of the incoming and reflected wave on the external wall of the cavity.

The maximum pressure in the Helmholtz cavity as a function of is shown in Fig. 4. Starting at frequency the pressure inside the cavity is and thus equal to the reference pressure as it should be. It then reaches a very high peak near the Helmholtz cavity resonance around , this is close to the theoretically predicted value of

. The difference can probably be attributed to the fact that the neck in our study is tapered and not straight. For large

values the pressure inside the cavity drops below the reference pressure and becomes rather low to reach a minimum around . A second peak can be observed near , which corresponds to the internal resonance frequency of a sphere. The maximum pressure does not necessarily occur in the main spherical part of the cavity, but can also occur in the neck part. Further peaks in the spectrum can be found associated with higher order internal resonance frequencies (not shown in Fig. 4, but the pressure profiles are shown in Fig. 3(e) and (f) for some higher resonances near and ). The results clearly show that a Helmholtz cavity can be simulated, for both the low ”Helmholtz resonance” as well as for the higher frequencies associated with the inner resonances of the spherical part.

## 3 Curl free vector Helmholtz equations

For a curl-free vector field, which also satisfies a Helmholtz equation , such as the case for the velocity field of a sound wave, we can introduce a potential . Then the framework of the previous section can be applied again. Thus curl free (or longitudinal waves) are relatively easy to describe using a single Helmholtz equation for the potential . Such a simple approach is no longer available for divergence free vector fields as we will see in the next section.

## 4 Divergence free vector Helmholtz equations; Electromagnetism

### 4.1 Introduction to electromagnetic scattering

Electromagnetic waves are typical examples of transverse waves. These waves satisfy the vector Helmholtz equation, but must simultaneously obey the zero divergence condition (a solenoidal vector field). Traditionally this has hindered the development of a simple and consistent framework just using Helmholtz equations alone and recourse had to be taken to Dyadic Green’s functions or current based theories [9], for example the Poggio and Miller [18], Chang and Harrington [4], and Wu and Tsai [36]

(PMCHWT) theory. In the frequency domain, if the medium is homogeneous and source free, the propagation of an electric field

obeys the following equations derived from the Maxwell’s equations  [2] equationparentequation

 ∇2E+k2E = 0, (4a) ∇⋅E = 0. (4b)

Eq. (4b) is Gauss’ law for electricity (assuming no volume charges). Eq. (4a) is clearly a 3D version of the Helmholtz equation and is essentially composed of three scalar Helmholtz equations. In Eq. (4a) (see pages 58-59 of [2]), the separate components of (i.e. , and

) only satisfy the scalar Helmholtz equation in a Cartesian coordinate system. The difficulty encountered in electromagnetic wave theory is that the divergence free condition in Eq. (

4b) must be satisfied simultaneously with the vector wave equation (4a).

Obviously, there is the conceptual advantage of dealing directly with the physical quantity of interest, the electric field vector , instead of with surface currents or hypersingular dyadic Green’s functions [9]. We will now introduce a recently developed field-only non-singular surface integral method to solve Eq. (4) straightforwardly. In a typical simulation, we solve Eq. (4) for the scattered electric field of an incoming field in the external domain, since only the scattered field satisfies the Sommerfeld radiation condition at infinity. The total electric field would then be , where is a plane wave in this work. If the object from which the scatter occurs is dielectric in nature, we also need to to solve for the transmitted electric field, inside the object.

### 4.2 Boundary conditions for perfect electric conductors

For the development of the numerical framework, it is convenient to express the field not only in global Cartesian components, but also in the local normal and tangential components on scatterer surfaces (since the boundary conditions are often given in those components) as:

 E=exEx+eyEy+ezEz=nEn+t1Et1+t2Et2. (5)

Here the normal vector, , is pointing out of the domain, and and are the two tangential vectors at the surface according to the convention . The unit vector in the -direction is (similar for and ). Thus for example the -component of the electric field can be expressed in terms of the normal and tangential components as:

 Ex=(n⋅ex)En+(t1⋅ex)Et1+(t2⋅ex)Et2=nxEn+t1xEt1+t2xEt2 (6)

with , and .

Perfect electric conductors are often used as an idealization of a metallic object under the illumination of an electromagnetic wave. For a perfect electric conductor, the boundary condition is that the tangential component of the total electric field is zero. The total field is the sum of the incoming and scattered field thus: equationparentequation

 Esct1=−Einct1;Esct2=−Einct2(Perfect Electric Conductor) (7a)

### 4.3 Divergence free condition implementation

The key to solve Eq. (4) is how to deal with the divergence free condition in Eq. (4b). An elegant way to ensure that the divergence free condition condition of the electric field is satisfied in the domain, is to ensure it is satisfied on the surface of an object. There are two things we need to investigate; firstly how do we ensure that the divergence free condition is satisfied on the surface and secondly, does that guarantee that the divergence free condition is satisfied everywhere else in the domain? These questions were answered in Sun et al. [34] and will be summarized here.

Let us start with the divergence of the electric field on the boundary of the domain:

 ∇⋅E=n⋅∂E∂n−κEn+∂Et1∂t1+∂Et2∂t2=0 (8)

with the curvature of the surface and the normal derivative, the tangential derivative in the tangential vector direction and similar for in the direction. , and are the normal and tangential components of the electric field on scatterer surfaces, respectively. A convenient way to see this is to write the electric field into its normal and tangential components on the surface as . Then the divergence can be written as . With (the curvature), . Also and . Thus: which is the desired equation.

Setting the divergence of the electric field on the boundary to zero is important since this will ensure that the divergence is also zero in the domain outside the object. This can be shown by realizing that if, for example is a solution of the Helmholtz equation, then is as well. Then and are solutions as well. The sum of several solutions of the Helmholtz equation will also obey the Helmholtz equation, then will also obey the Helmholtz equation. Since on the surface, the normal derivative must also be zero (since these two quantities are related by the boundary element framework). Since the Helmholtz equation is elliptic in nature, the whole field must be divergence free.

There are other ways to ensure the divergence free condition, for example using the vector identity , with the position vector. Thus one more Helmholtz equation with will occur see [14, 28] for more details.

The condition in Eq. (8) can replace Eq. (4b). On the surface of a perfect electric conductor, the total tangential fields and , Eq. (8) becomes:

 n⋅∂Etot∂n=κEtotn. (9)

Thus for the scattered field we find:

 n⋅∂Esc∂n=κEscn+κEincn−n⋅∂Einc∂n. (10)

Similar to the decomposition of the electric field (Eq. 6), a decomposition into normal and tangential components can be done for the normal derivatives as:

 ∂Ex∂n=nx(n⋅∂Esc∂n)+t1x(t1⋅∂Esc∂n)+t2x(t2⋅∂Esc∂n) (11)

Thus the matrix equations , , can be combined and expressed in terms of , and for each node as a matrix system using Eq. 10:

 (12)

where , and similar for the and terms on the right hand side. It is also possible to work with the magnetic field instead of the electric field and get a similar matrix system.

The above formulation has been thoroughly tested against the Mie scattering solution for a sphere [34]. Another example of scattering of an electromagnetic wave off a perfect conducting cube with rounded corners is shown in Fig. 5. A cube with the edges having length is rotated in such a way to have one vertex pointing in the direction of an incoming wave with . A complex pattern of interaction between the incoming and scattered waves can be seen in the plane at . Since the right hand side exhibits a face, the perturbed electric field is rather chaotic there. The left hand side has an edge in the plane where we plot the electric field and the perturbation of the incoming field is much less there. Besides the instantaneous electric field (in red vectors), we have also plotted the absolute value of the electric field (which is time independent) as a scalar color plot.

### 4.4 Dielectric formulation

Another type of scatterers that often appear are dielectric objects. When a plane wave is scattered by a (closed) dielectric object, contrary to the perfect electric conductor case, we must also calculate the transmitted wave into the dielectric object. This means that we also need to solve another set of Helmholtz equations for the dielectric with a different wave number as opposed to the external wave number . To simplify the equations, we assume that for the permeabilities: (if these are different see Sun et al. [33]), while the permitivities of both domains ( and ) are different.

We seek again a solution of the external scattered electric field in terms of the incoming and transmitted field into the dielectric object. The boundary conditions for such a system [26] are (no surface charges or currents): equationparentequation

 ϵinEtrn=ϵout(Eincn+Escn);jump in normal % electric field (13a) Etrt1=Einct1+Esct1;Etrt2=Einct2+Esct2;Tang. continuity of electric field (13b) μinHtrn=μout(Hincn+Hscn);jump in normal magnetic field (13c) Htrt1=Hinct1+Hsct1;Htrt2=Hinct2+Hsct2Tang. continuity of magn. field. (13d)

Since Eq. (13a) was derived from the divergence free electric field conditions on both sides of the boundary, we only have to ensure that while will then be automatically guaranteed. The above boundary conditions are not all independent as Eq. (13b) guarantees that Eq. (13c) is satisfied.

The next task at hand is to convert the tangential continuity conditions of the magnetic field in Eq. (13d) in terms of the electric field. Since this leads to and (note the inversion of subscripts ‘1’ and ‘2’ here and the plus and minus signs). We can then get:

 (14)

Since , it will lead to and and gives Eq. 14 using Maxwell’s equations.

The rest of the derivation for a dielectric scattering problem will result in a matrix system from which the quantities , , , , and are solved. From this the Cartesian components of the electric field and its normal derivatives can easily be reconstructed. The derivation is relatively straightforward and is given in Appendix C.

As an illustrative example of a dielectric scattering problem, we consider an oblate spheroidal shaped object which functions as a lens, under illumination with a plane incoming wave at an angle with respect to the optical axis of the lens. The long axis of the lens has dimension and the short axis is . We first study the case when the incoming wave is at an angle of 45 degrees. For relatively low frequencies as shown in Fig 6, we can see that the electric field inside the lens is still more or less parallel. No focusing effect is observed for and , neither for and , which could be expected for the long wavelengths associated with these wavenumbers. We keep the material constants the same in all these examples (thus the ratio is fixed). Another example is shown next with , in Fig. 7(a), and we can see some disturbance in the electric field in the lower right corner. This disturbance becomes bigger in Fig. 7(b) where we have used and . In Fig. 8, we show two more examples for the parameters sets , and , . We can see that the electric field vectors in the region to the right and below the lens are getting larger. In order to see more clearly what is going on here, we have plotted the absolute value of the complex valued electric field in Fig. 9(a). We can now clearly see the focal area of the lens. In Fig. 9(b) we have also plotted the length of the instantaneous electric field (thus the absolute value of the real part of the electric field vector), the incoming and focused waves can clearly be observed. As a reference case, we have also plotted these parameter for waves travelling at a zero angle towards the lens, Figs. 9(c) and (d).

In Fig. 10(a), the electric field lines around a dielectric object near the zero frequency limit are shown. Contrary to other numerical methods, such as surface currents methods, our non-singular field only surface integral framework has no issues when approaches zero (the long wavelength or electrostatic limit). We have chosen and for this particular example. Another example is shown in Fig. 10(b), but now with , thus the object essentially becomes invisible and the electric field lines are not being disturbed by the object. Here we have used .

## 5 Other physical systems with Helmholtz equations

One physical phenomenon of special scientific interest is acoustic waves in solid materials since transverse and longitudinal waves, each with a different wave numbers, can occur simultaneously. The transverse waves satisfy a zero divergence, while for the longitudinal waves, the curl is zero. These waves travel at different speeds (the longitudinal waves always travel faster than the transverse waves). Assume a linear isotropic homogeneous elastic material, with zero body force. The material has density and the Lamé constants and ( is the shear modulus) as material elastic properties.

In the time domain waves occurring in the material can be described with the following elastodynamic equation where , with the displacement vector, and

the stress tensor. The superscript

indicates the transpose of the tensor and is the unit tensor. This equation essentially expresses the force balance on an infinitesimal volume element; with inertial forces on the left and elastic forces on the right. In the frequency domain this leads to the Navier equation (assuming again a harmonic time dependency [5, 1]: where and . and are the transverse and longitudinal wave speeds, also often referred to as shear wave velocity and dilatational wave velocity. Alternatively, using the vector identity , Navier equation can be written as:

 (k2Tk2L−1)∇∇⋅u+∇2u+k2Tu=0 (15)

with and the longitudinal and transverse wave numbers, respectively. We will now describe the waves occurring in a linear elastic medium with the help of Helmholtz equations alone, using the Helmholtz decomposition. The displacement field is decomposed into longitudinal and transverse components where (transverse) and (longitudinal) [15], then Eq. (15) becomes:

 ∇2uT+k2TuT=0 ;∇⋅uT=0 (16) ∇2Φ+k2LΦ=0 .

Here we have chosen a potential function to automatically satisfy the curl-free condition for .

We have demonstrated that the Helmholtz equation is classically used to describe acoustic waves, and the wave number is often a real number. However, can be complex, and Eq. (1) then represents an acoustic wave in a medium with absorption. There are other types of variations of the Helmholtz equation as well. For example, when is imaginary (thus is negative), Eq. (1) becomes the Debye-Hückel model for the molecular electrostatics potential , satisfying , in colloidal systems in which is the inverse of the Debye length [27]. When is purely imaginary, Eq. (1) describes diffusion or heat transfer and no longer wave phenomena. Nevertheless, the same numerical boundary element framework (see Sec. 6) can still be applied.

## 6 Non-singular boundary element method for Helmholtz equation

All the above classic problems can, in principle, be solved just based on one or a set of Helmholtz equations. Only for a very limited number of cases an exact theoretical solution can be found. Therefore an efficient numerical method to solve the Helmholtz equation is needed. A good candidate is the boundary integral method (or boundary element method when it is used in discretized form) as it can represent surface geometry accurately and reduces a three-dimensional (3D) problem to a 2D problem.

However, the conventional formulation of this method involves the presence of singularities originating from the Green’s function. Such singularities have no physical basis but are a purely mathematical artefact. They can cause mathematical complexities, such as divergent integrals in the classical sense. We describe an easy to implement way to remove this singular behavior, while retaining the elegance and accuracy of the boundary element method.

Though the conventional boundary integral method (CBIM) in Eq. (2) is often used to solve the Helmholtz equation, Eq. (1), there are two main problems worth mentioning. Firstly, the integrands of CBIM in Eq. (2) are singular when . Although these singular integrands can be integrated, their practical treatment in numerical implementations is not straightforward [6, 35, 21, 22]. Secondly, computation of the solid angle is also not straightforward for collocation nodes, which makes the use of surface elements other than planar constant elements rather tedious.

Recently, the non-singular boundary integral methods (NS-BIM) or boundary regularised integral equation formulations (BRIEF) have been developed for applications in fluid mechanics [12, 29, 30, 31], acoustics [32, 11], molecular electrostatics [27], electromagnetics [14, 28, 13] and linear elastic waves [15]. The objective of the non-singular boundary integral method is to analytically remove the singularities in and as as well as the solid angle .

The singular behavior of is the same as that of the free-space Green’s function of the Laplace equation: with , since with

 ΔG≡exp(ik|x−x0|)−1|x−x0| (17)

which is regular when [37]. The same analysis and conclusion can be made for (see A for more details). Using this fact, we start with a known function that satisfies the Laplace equation as, , and the conventional boundary integral representation of it is

 c(x0)ψ(x0)+∫Sψ(x)∂G0∂n dS(x)=∫S∂ψ(x)∂nG0 dS(x). (18)

Let us assume that has the form

 ψ(x)=g(x)ϕ(x0)+f(x)∂ϕ(x0)∂n (19)

where and are constants in this context, and and satisfy the following conditions equationparentequation

 ∇2g(x)=0,limx→x0g(x)=1,limx→x0∂g(x)∂n=0; (20a) ∇2f(x)=0,limx→x0f(x)=0,limx→x0∂f(x)∂n=1. (20b)

Introducing Eq. (19) into Eq. (18) and then subtracting the result from Eq. (2):

 ∫S[ϕ(x)∂Gk∂n−ϕ(x0)g(x)∂G0∂n+ϕ(x0)∂g(x)∂nG0] dS(x) =∫S[∂ϕ(x)∂nGk−∂ϕ(x0)∂n∂f(x)∂nG0+∂ϕ(x0)∂nf(x)∂G0∂n]d% S(x). (21)

Eq. (21) is the non-singular boundary integral equation for the Helmholtz equation (1) in which the integrands are all regular as and the term with the solid angle is eliminated. It is worth mentioning that (21) is valid for either , being real, purely imaginary, or a random complex number. Note that each node has its own and function. Further proof that Eq. (21) no longer contains singular terms can be found in Appendix A.

There are many possible choices for functions and that can satisfy the conditions in Eq. (20) which ensure Eq. (21) is free of singularities. For instance a constant and a linear function [27]: equationparentequation

 g(x)=1, ∂g(x)∂n=0; (22a) f(x)=n(x0)⋅(x−x0), ∂f(x)∂n=n(x0)⋅n(x). (22b)

Note that for external problems, the integrals with the above choice of with and in Eq. (22) over the closed surface at infinity do not vanish. Nevertheless, the integral value can be found analytically as and thus this value should be added to the left-hand side of Eq. (21) for external problems (see Appendix  B).

It is worth emphasizing that the above framework enables us to use higher order surface elements, such as quadratic elements, together with the standard Gauss integration methods in the numerical implementation.

## 7 Discussion

The following points are worth noting concerning the above described non-singular boundary element framework:

• Since we use the same numerical framework for acoustics and electromagnetics, we can use the same codes with quadratic elements for both. For the electromagnetic problem we do not use RWG [19] elements. Thus, anyone with a Helmholtz boundary element solver for soundwaves, can in principle extend this to electromagnetic scattering simulations. Moreover, our method is stable for electromagnetic problems in the long wavelength limit.

• When solving the vector wave problems, such as in electromagnetics, iterative solvers can also be used only requiring one matrix to be solved, although careful attention need to be paid to the convergence behavior.

• Concerning problems in the time domain: a time dependent problem can always be decomposed into its Fourier components. Then those components can be treated in the frequency domain and the solution can be reassembled in the time domain. This was actually implemented for soundwaves and for electromagnetic pulses in Klaseboer et al. [11, 13].

When using a boundary element method to numerically solve the Helmholtz equation for external problems, one challenge is that non-physical solutions will show at certain discrete frequencies [16, 17]. This numerical issue is known as the the occurrence of spurious solutions and occurs at fictitious frequencies. In fact with our desingularized boundary element method, the probability of hitting a spurious solution is actually very low. In fact, it is even difficult to find such a frequency when specifically looking for it (see the figures in [10]). The two most popular methods to deal with fictitious frequencies are the CHIEF method proposed by Schenck [23] and the Burton-Miller method [3]. The CHIEF method uses additional internal points creating an over-determined system. Therefore an additional computational technique, such as the least square method must be employed. Furthermore, the successful removal of spurious solutions at fictitious frequencies cannot be fully guaranteed by this method. The Burton-Miller method claims that the spurious solutions at the fictitious frequencies disappear. However, the drawback is that it contains an integral with a hyper-singularity.

## 8 Conclusions

In this article, in memory of the famous German scientist Hermann von Helmholtz for his 200-year birthday, we revisited the wide engineering applications of the elliptic partial differential equation named after him, the Helmholtz equation, and introduced a robust, efficient and simple non-singular boundary integral (element) method to solve the Helmholtz equation. We have shown that the Helmholtz equation, classically used to describe scalar quantities, such as the potential or pressure in sound wave simulations in the frequency domain, can also be used for other physical systems. We have given an example for the case of electromagnetic waves, which is essentially a set of three coupled Helmholtz equations, one for each Cartesian component of the electric field. The main difficulty that is faced historically, is how to enforce the divergence free condition of the electric field. We have shown that this condition can be satisfied relatively easily. Some examples are shown for perfect electric conductors and dielectric objects. For the scalar Helmholtz equation we have chosen an example with transducers and reflectors and a Helmholtz cavity.

A major advantage in solving actual problems is the development of an entirely desingularized boundary element method for the Helmholtz equation without compromising on accuracy or efficacy. The main idea behind this being that if the physical problem does not have any singularities, then the mathematical model also should not have any singular behavior. The singular behavior in the classical boundary element method originates from the Green’s function. It turns out that if a carefully chosen analytical solution, with the same singular behavior, is subtracted from the original equation, a totally non singular framework is created. All elements, including the previously singular ones, can be integrated with standard Gauss integration. It also enables us to implement higher order quadratic elements with quadratic shape functions with ease.

## Appendix A Notes on non-singular boundary element method

A closer inspection of the desingularized equation (21) will be performed. The first and second terms in Eq. (21), i.e. contain no singularities as can be shown by taking a Taylor expansion at as:

 ϕ(x)∂Gk∂n−ϕ(x0)g(x)∂G0∂n (23) ≈ ϕ(x0)∂Gk∂n+∇ϕ⋅(x−x0)∂Gk∂n−ϕ(x0)∂G0∂n−ϕ(x0)∇g⋅(x−x0)∂G0∂n = ϕ(x0)[∂Gk∂n−∂G0∂n]+∇ϕ⋅(x−x0)∂Gk∂n−−ϕ(x0)∇g⋅(x−x0)∂G0∂n.

The expressions for and are (with ):

 ∂Gk∂n =[ikr−1]eikrr3(x−x0)⋅n(x);∂G0∂n =−1r3(x−x0)⋅n(x) (24)

The terms in between brackets in Eq. (23) can thus be written with Eq. (24) as:

 [∂Gk∂n−∂G0∂n]=([ikr−1]eikr+1)1r3(x−x0)⋅n(x). (25)

Here and the term behaves as since in the limit of going to , and are perpendicular. Thus Eq. (25) does not contain any singular term. The term with in Eq. (23) goes as and thus cancels out the behavior of . Similar for .

For the third term in Eq. (21) with , as long as approaches zero in a linear manner, it will cancel out the singularity of .

A similar analysis can be performed for the right hand side of Eq. (21). For the first two terms on the right hand side, we can write with again two similar Taylor expansions: and , and then

 ∂ϕ(x)∂nGk−∂ϕ(x0)∂n∂f(x)∂nG0 (26) ≈ ∂ϕ(x0)∂n[Gk−G0]+∇∂ϕ∂n⋅(x−x0)Gk−∇∂f∂n⋅(x−x0)G0

where we have used . The term with is regular, and the terms with the gradients contain both , which thus cancel out the singularity from and , respectively.

Since approaches zero as , it cancels out the third term with on the right hand side of Eq. (21).

## Appendix B Remarks on integrals at infinity

There are two remarks that can be made about the integrals at infinity that occur for external Helmholtz problems in the boundary element framework. The first one refers to the application of the Sommerfeld radiation condition [25, 24] and applies to the conventional boundary element method as given in Eq. (2) where the two integrals at infinity are:

 ∫∞ϕ∂Gk∂ndS−∫∞∂ϕ∂nGkdS=∫∞[ϕ∂Gk∂n−∂ϕ∂nGk]dS. (27)

Suppose we draw a very large sphere with radius and consider this as ‘infinity’ then and . Since and , then:

 (28)

The term with vanishes as since decays at least as fast as . Then:

 (29)

which can be shown to be zero with the help of the Sommerfeld radiation condition [25, 24] expressing the fact that only outgoing waves are allowed (and realizing that at , ) as

 limr→∞r[∂ϕ∂r−ikϕ]=0. (30)

This is the formula as it appears in the original Sommerfeld article [25] (his Eq. (21)).

The second place where due care has to be taken with integrals at infinity is in Sec. 6, where the desingularized boundary element method is described. In this case, there is actually a term that appears from one of the integrals at infinity. There are four integrals that now need to be considered in Eq. (21):

 ∫∞ϕ(x0)g(x)∂G0∂n dS(x) ;∫∞ϕ(x0)∂g(x)∂nG0 dS(x) (31) ∫∞∂ϕ(x0)∂n∂f(x)∂nG0dS(x) ;∫∞∂ϕ(x0)∂nf(x)∂G0∂ndS(x).

With the particular choice of from Eq. (22), it can be seen that the second integral is zero immediately (). The first integral will become, with (as remarked earlier ) and

 ∫∞ϕ(x0)g(x)∂G0∂n dS(x)=ϕ(x0)1R3R4πR2=4πϕ(x0). (32)

This is the value mentioned in Sec. 6. The other two remaining integrals will give zeroes since and , with the angle between the normal vector on the surface at and the normal vector on the surface at infinity. This will cause the contributions of the integrals at infinity to cancel for these two last integrals. Thus:

 ∫∞∂ϕ(x0)∂n∂f(x)∂nG0dS(x)=0and∫∞∂ϕ(x0)∂nf(x)∂G0∂ndS(x)=0. (33)

## Appendix C Derivation for the dielectric case

We need to expand and