Helmholtz equation and non-singular boundary elements applied to multidisciplinary physical problems

Evert Klaseboer and Qiang Sun

1 Institute of High Performance Computing，1 Fusionopolis Way，Singapore 138632，Singapore

2 Australian Research Council Centre of Excellence for Nanoscale BioPhotonics，School of Science，RMIT University，Melbourne，VIC 3001，Australia

Abstract 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 curlfree constraints.There is further a need for a numerical method that is accurate，reliable and takes into account radiation conditions at infinity.The classical boundary element method satisfies the last condition，yet one has to deal with singularities in the implementation.We review here how a recently developed singularity-free three-dimensional boundary element framework with superior accuracy can be used to tackle such problems only using one or a few Helmholtz equations with higher order(quadratic) elements which can tackle complex curved shapes.Examples are given for acoustics (a Helmholtz resonator among others) and electromagnetic scattering.

Keywords: acoustics，Helmholtz cavity，electromagnetics，scattering，boundary integral method

1.Introduction

Hermann von Helmholtz (1821–1894) [1，2] was born 200 years ago on 31 August 1821.The scalar Helmholtz equation，∇2φ+k2φ=0，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 sound waves)，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 (3D) 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(BEM) 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 section 2，with as examples，among others，the Helmholtz resonator in acoustics.In section 2.2，the recently developed non-singular boundary integral method(NS-BIM) to solve 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.After a brief discussion on curl-free vector fields in section 3，we discuss divergence-free vector fields in section 4，with particular attention to electromagnetic scattering problems.Then，the discussion and conclusions are given in sections 6 and 7，respectively.

2.Scalar Helmholtz equation

2.1.Acoustics

The wave equation describing the propagation of a quantity φ′is in the time domainwith c the wave speed and t time.Assuming a harmonic time dependencewith angular frequency ω and unit imaginary number‘i’，the wave equation transforms into the wellknown Helmholtz equation as:

with k=ω/c.This is the equation as it appeared in Helmholtz’s book[1]on‘acoustics in pipes with open ends’while studying organ pipes as equation (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，k is a real valued number，often referred to as a wave number in wave physics，while φ is a complex quantity.The boundary conditions are usually that φ is given (Dirichlet condition)，or its normal derivative dφ/dn=n·∇φ is given (Neumann condition)with n 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.

2.2.Non-singular BEM for the Helmholtz equation

Only for a very limited number of cases，an exact theoretical solution can be found for the Helmholtz equation.Therefore an efficient numerical method to solve it is required.A good candidate is the boundary integral method(or BEM when it is used in discretized form) as it can represent surface geometry accurately and reduces a 3D problem to a two-dimensional problem[3].

The Helmholtz equation is an elliptic equation.This means that the properties on the surface(s) 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 equation (1):

whereGk≡Gk(x，x0) =is the Green’s function，x0is a source point on the surface，x is the computation point on the surface S，and c(x0) is a constant related to the solid angle at x0.

Though the conventional boundary integral method(CBIM) in equation (2) is often used to solve the Helmholtz equation，equation (1)，there are two main problems worth mentioning.Firstly，the integrands of CBIM in equation (2)are singular when x →x0.Although these singular integrands can be integrated，their practical treatment in numerical implementations is not straightforward [4–7].Such singularities have no physical basis but are purely mathematical artifacts.Secondly，the computation of the solid angle c(x0)is also not straightforward for collocation nodes，which makes the use of surface elements other than planar constant elements rather tedious.

The treatment and even elimination of the singularities has been a subject of intense study from the beginning of the BEM.Well-known are the so-called constant potential subtraction methods that remove the singularity on the left-hand side of equation (2) (see for example [8])，yet the singular behavior on the right-hand side of equation (2)still has to be dealt with.Another way to avoid the singularity issue is to use the method of fundamental solutions or equivalent source method [9，10] where the sources are put outside of the domain.

We have developed an advanced non-singular BEM in which the solid angle and singularities in the integrals are removed analytically.This improvement has made the BEM much easier to implement，especially for higher order quadratic elements.

Recently，the NS-BIM or boundary regularised integral equation formulations have been developed for applications in fluid mechanics [11–14]，acoustics [15–17]，molecular electrostatics [18]，electromagnetics [19–21]，interactions between light and matter [22，23] and linear elastic waves[24].The objective of the NS-BIM is to analytically remove the singularities in Gk(x，x0) and ∂Gk(x，x0)/∂n as x →x0as well as the solid angle c(x0).

The singular behavior of Gk(x，x0) is the same as that of the free-space Green’s function of the Laplace equation:G0≡G0(x，x0)=1/|x-x0|with ∇2G0(x，x0)=-4πδ(x-x0)，since Gk(x，x0)≡G0(x，x0)+ΔG with

which is regular when x →x0[25].The same analysis and conclusion can be made for ∂Gk(x，x0)/∂n(see appendix A for more details).Using this fact，we start with a known function ψ(x)that satisfies the Laplace equation as，∇2ψ(x)=0，and the conventional boundary integral representation of it is

Let us assume that ψ(x) has the form

where φ(x0) and ∂φ/∂n(x0) are constants in this context，and g(x) and f(x) satisfy the following conditions

Introducing equation (5) into (4) and then subtracting the result from equation (2):

Equation (7) is the non-singular boundary integral equation for the Helmholtz equation，equation(1)，in which the integrands are all regular as x →x0and the term c(x0)φ(x0)with the solid angle is eliminated.It is worth mentioning that equation (7) is valid for either k=0，k being real，purely imaginary，or a random complex number.Note that each node x0has its own f and g function.Further proof that equation(7)no longer contains singular terms can be found in appendix A.

There are many possible choices for functions g(x) and f(x) that can satisfy the conditions in equation (6) which ensure equation (7) is free of singularities.For instance a constant and a linear function [18]:

Note that for external problems，the integrals with the above choice of φ(x) with g(x) and f(x) in equation (8) over the closed surface at infinity do not vanish.Nevertheless，the integral value can be found analytically as 4πφ(x0) and thus this value should be added to the left-hand side of equation (7) 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 for all nodes (including the previously singular ones).In all the simulation results that follow，we have used the desingularized BEM as described above with quadratic elements，except for figure 11(b)where we have used a desingularized version of the Burton–Miller BEM.The Sommerfeld radiation condition at infinity，i.e.only outgoing waves are allowed，is automatically satisfied (see also appendix B).

Assuming the surface is divided into surface elements with a total of N nodes，then all potentials and their normal derivatives given in equation (7) are related by the following matrix system:

The N×N matricesG andH are the numerical matrix equivalent of the boundary element integrals in equation (2)，andandare column vectors of length N in which each component φ(i) and ∂φ(i)/∂n belongs to φ and ∂φ/∂n，respectively，of the i’th node located at x0(i)with i=1，…，N.There are，by the way，other methods to solve the Helmholtz equation，such as with finite elements.

2.3.Acoustic transducer

To demonstrate some interesting wave phenomena in acoustics，in figure 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 figure 1(a)，the absolute pressure is shown，while in figure 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.

Figure 1.Sound transducer on top，receiving (stationary and rigid) rotated parabolic bowl below.(a) Absolute pressure amplitude and (b)instantaneous pressure profile.The graphs show the physical principles of generation and reflection of sound.Also note the focal areas of the two bowls.

2.4.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 inside a neck at the top.The Helmholtz resonator is assumed to be acoustically hard: ∂φtot/∂n=0 on its surface (thus ∂φsc/∂n=-∂φinc/∂n)，where the superscripts‘tot’，‘inc’ and ‘sc’ indicate the total，incoming and scattered wave，respectively.In figure 2 the Helmholtz cavity as used in our simulations is shown，both in full 3D view and in crosssectional view.The spherical part has an outer radius of 1.08a，an inner radius of 0.92a，a neck with length L=0.67a and the neck tapers down from a value of 0.32a to 0.168a at the top.This will give us a theoretical resonance that should occur at0.208 with A=(0.162a)2π the opening area of the neck，V=4/3(0.92a)3π the inner volume of the sphere [26].

Figure 2.Helmholtz cavity as used in the simulations consisting of a central spherical part and a neck with an opening (a) full 3D view.(b) Cross-sectional view revealing the interior of the cavity.The theoretical resonance occurs atka= with a the radius of the sphere，A the opening area of the neck，V the volume of the sphere and L the length of the neck.For the current cavity this would be around ka=0.225.

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

Figure 3.(a)ka=0.05，the pressure inside the cavity p=1.05 is slightly higher than the reference pressure outside p=1.00.(b)ka=0.23 near the Helmholtz cavity resonance frequency.(c) ka=0.5 The pressure p=0.25 is considerably lower than the reference pressure.(d) ka=π，near the first internal resonance frequency of a sphere.(e) ka=2.05π near the second resonance frequency of the sphere.(f) ka=4.3π，another internal resonance frequency.

The maximum pressure in the Helmholtz cavity as a function of ka is shown in figure 4.Starting at frequency ka=0 the pressure inside the cavity is 1.0 and thus equal to the reference pressure as it should be.It then reaches a very high peak near the Helmholtz cavity resonance around ka=0.225，this is close to the theoretically predicted value of ka=0.208.The difference can probably be attributed to the fact that the neck in our study is tapered and not straight.For large ka values the pressure inside the cavity drops below the reference pressure and becomes rather low to reach a minimum around ka=1.5.A second peak can be observed near ka=π，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 figure 4，but the pressure profiles are shown in figures 3(e) and (f) for some higher resonances near ka=2π and ka=4π).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.

Figure 4.The maximum pressure inside the Helmholtz cavity as a function of ka.A peak appears corresponding to the cavity resonance around ka=0.225 with a maximum of several 100's times the reference pressure.The peak at the internal resonance of the sphere is also visible around ka=π.

3.Curl free vector Helmholtz equations

For a curl-free vector field，∇×u=0 which also satisfies a Helmholtz equation ∇2u+k2u=0，such as the case for the velocity field of a sound wave，we can introduce a potential u=∇φ.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 [27]，for example the Poggio and Miller [28]，Chang and Harrington[29]，and Wu and Tsai [30] (PMCHWT) theory，although alternative methods exist，such as the method of fundamental solutions [31，32].In the frequency domain，if the medium is homogeneous and source free，the propagation of an electric field E obeys the following equations derived from Maxwell’s equations [33]

Equation (10b) is Gauss’ law for electricity (assuming no volume charges).Equation (10a) is clearly a 3D version of the Helmholtz equation and is essentially composed of three scalar Helmholtz equations.In equation(10a)(see pages 58–59 of[33])，the separate components of E only satisfy the scalar Helmholtz equation in a Cartesian coordinate system(i.e.Ex，Eyand Ez).The difficulty encountered in electromagnetic wave theory is that the divergence-free condition in equation (10b) must be satisfied simultaneously with the vector wave equation (10a).

Obviously，there is the conceptual advantage of dealing directly with the physical quantity of interest，the electric field vector E，instead of with surface currents or hypersingular dyadic Green’s functions [27].We will now introduce a recently developed field-only non-singular surface integral method to solve equation(10)straightforwardly.In a typical simulation，we solve equation (10) for the scattered electric field Escof an incoming field Eincin the external domain，since only the scattered field satisfies the Sommerfeld radiation condition at infinity.The total electric field would then be Etot=Einc+Esc，where Eincis a plane wave in this work.If the object from which the scatter occurs is dielectric in nature，we also need to solve for the transmitted electric field，Etrinside 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 the scatterer surface (since the boundary conditions are often given in those components) as:

Here the normal vector，n，is pointing out of the domain，and t1and t2are the two tangential vectors at the surface according to the convention n=t1×t2and (t1·t2)=0.The unit vector in the x-direction is ex(similar for y and z).Thus for example the x-component of the electric field can be expressed in terms of the normal and tangential components as:

with nx=n·ex，t1x=t1·exand t2x=t2·ex.

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 fields thus:

4.3.Divergence-free condition implementation

The key to solving equation (10) is how to deal with the divergence-free condition in equation (10b).An elegant way to ensure that the divergence-free 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:

with κ the curvature of the surface and ∂/∂n=n·∇the normal derivative，∂/∂t1=t1·∇，the tangential derivative in the tangential vector t1direction and similar for ∂/∂t2in the t2direction.En，Et1and Et2are the normal and tangential components of the electric field on the scatterer surface，respectively.A convenient way to prove equation (14) is to write the electric field into its normal and tangential components on the surface as E=Enn+Et1t1+Et2t2.Then the divergence can be written as

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，Exis a solution of the Helmholtz equation，then ∂Ex/∂x is as well.Then∂Ey/∂y and ∂Ez/∂z are solutions as well.The sum of several solutions of the Helmholtz equation will also obey the Helmholtz equation，then ∇·E will also obey the Helmholtz equation.Since ∇·E=0 on the surface，the normal derivative ∂(∇·E)/∂n must also be zero(since these two quantities are related by the boundary element framework).Because 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 2(∇·E)≡∇2(x·E)+k2(x·E)=0，with x the position vector.Thus one more Helmholtz equation for (x·E) will also guarantee(∇·E)=0，see [19，20] for more details.

The condition in equation(14)can replace equation(10b).On the surface of a perfect electric conductor，the total tangential fieldsandequation (14) becomes:

Thus for the scattered field，we find:

Similar to the decomposition of the electric field in equation(12)，a decomposition into normal and tangential components can be done for the normal derivatives as:

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 figure 5.A cube with the edges having length 2a is rotated in such a way to have one vertex pointing in the direction of an incoming wave with ka=10.A complex pattern of interaction between the incoming and scattered waves can be seen in the plane at y=0.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 timeindependent) as a scalar color plot.

Figure 5.Electromagnetic Scattering of a plane wave from a cube with rounded corners with lengths 2a and ka=10.The incoming plane wave travelling in the positive z-direction is polarized horizontally.The instantaneous electrical field vectors and amplitude are plotted on the plane y=0.The cube is rotated in such a way that one of its corners points towards the incoming wave.Note the strong reflection at the right side of the cube and the ‘shadow’ region behind the cube.

4.4.Dielectric formulation

Another type of scatterer that often appears is of dielectric nature.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 need to solve another set of Helmholtz equations for the dielectric with a different wave numberas opposed to the external wave numberTo simplify the equations，we assume that for the permeabilities: μin=μout(if these are different see Sun et al [35])，while the permittivities of both domains (∈inand ∈out) are different.

We seek again a solution of the external scattered electric field Escin terms of the incoming Eincand transmitted fieldEtrinto the dielectric object.The boundary conditions for such a system [36] are (with no free surface charges or currents):

Since equation (19a) was derived from the divergencefree electric field conditions on both sides of the boundary，we only have to ensure that ∇·Esc=0 while ∇·Etr=0will then be automatically satisfied (∇·Einc=0 already).The above boundary conditions are not all independent as equation (19b) guarantees that equation (19c) is satisfied.

The next task at hand is to convert the tangential continuity conditions on the magnetic field in equation (19d) in terms of the electric field.Since

this leads to t2·[n×H]=Ht1and t1·[n×H]=-Ht2(note the inversion of subscripts ‘1’ and ‘2’ here and the plus and minus signs).We can then get (by replacing H in [n×H]，using Faraday’s law of induction in the frequency domain:∇×E=iωμH):

Since

it will lead to

and results in equation (20).

The rest of the derivation for a dielectric scattering problem will result in a 6N×6N matrix system from which the quantitiesn·∂Esc/∂n，t1·∂Esc/∂n and t2·∂Esc/∂n 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 2a and the short axis is a.We first study the case when the incoming wave is at an angle of 45 degrees.For relatively low frequencies as shown in figure 6，we can see that the electric field inside the lens is still more or less parallel.No focusing effect is observed for kina=1.25 and kouta=1.00，neither for kina=2.5 and kouta=2.0，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 kin/koutis fixed).Another example is shown next with kina=5，kouta=4 in figure 7(a)，and we can see some disturbance in the electric field in the lower right corner.This disturbance becomes bigger in figure 7(b) where we have used kina=7.5 and kina=6.In figure 8，we show two more examples for the parameters sets kina=10，kouta=8 and kina=15，kouta=12.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 figure 9(a).We can now clearly see the focal area of the lens.In figure 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)，and the incoming and focused waves can clearly be observed.As a reference case，we have also plotted these parameters for waves travelling at a zero angle towards the lens，figures 9(c) and (d).

Figure 6.Instantaneous electric field vectors on the plane y=0 for a dielectric lens object with a/b=2 and the longest axis is 2a centered at(x，y，z)=(0，0，0) under irradiation at an angle of 45 degrees (from top left to bottom right) (a) kina=1.25，kouta=1.0;(b) kina=2.5，kouta=2.0.

Figure 7.As figure 6 but for higher ka numbers:(a)kina=5，kouta=4;(b)kina=7.5，kouta=6.Focusing effects start to appear behind the lens (in the lower right corner).

Figure 8.As figure 6 but for even higher ka numbers: (a) kina=10，kouta=8 and (b) kina=15，kouta=12 (see also figure 9).The ‘lens’gradually focuses the waves as the wavenumber gets higher.

In figure 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 current methods，our non-singular field-only surface integral framework has no issues when ka approaches zero (the long wavelength or electrostatic limit).We have chosen kina=0.002 and kouta=0.001 for this particular example.Another example is shown in figure 10(b)，but now with ∈in=∈out，thus the object essentially becomes invisible and the electric field lines are not being disturbed by the object.Here we have used kouta=kina=0.001.

Figure 9.Electric field for a dielectric lens with a/b=2 under irradiation at an angle of 45 degrees with kina=15 and kouta=12(from top left to bottom right): (a) absolute value of the electric field; (b) absolute instantaneous value of the electric field.Under zero degrees irradiation，(c) absolute value of the electric field; (d) absolute instantaneous value of the electric field.

Figure 10.Electric field lines for a near zero frequency case(ka →0)for an ellipsoid dielectric object with longest length 2a under irradiation at an angle of 45 degrees.(a) With kouta=0.001 and kina=0.002.Note that the internal electric field is parallel，but not parallel to the external field at infinity.(b)With kouta=kina=0.001，the ellipsoid essentially becomes a transparent object.The examples clearly show that there is no ‘zero frequency catastrophe’ for this method.

5.Other physical systems with Helmholtz equations

One physical phenomenon of special scientific interest concerns acoustic waves in solid materials since transverse and longitudinal waves，each with a different wave number，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.

with kL=ω/cLand kT=ω/cTthe longitudinal and transverse wave numbers，respectively.We can 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 u=uT+uLwhere ∇·uT=0 (transverse)and ∇×uL=0(longitudinal)[24]，then equation(23)becomes:

Here we have chosen a potential function uL=∇Φ to automatically satisfy the curl-free condition for uL.The above framework was used recently to find an analytical solution for a vibrating rigid sphere(i.e.moving periodically up and down in an infinite elastic material)[24]and for the same sphere but surrounded by an additional elastic shell [39].

The Helmholtz equation is classically used to describe acoustic waves，and the wave number k is often a real number.However，k can be complex，and equation(1)then represents an acoustic wave in a medium with absorption.As an example we can mention the description of acoustic boundary layers around a sphere in a viscous liquid[40]，which was based on the Navier equation (23)，but with complex k's.There are other types of variations of the Helmholtz equation as well.For example，

when k=iκDHis imaginary (thus k2is negative)，equation (1)becomes the Debye–Hückel model for the molecular electrostatic potential Ψ，satisfying ∇2Ψ -Ψ = 0，in colloidal systems in which κDHis the inverse of the Debye length [18].When k2is purely imaginary，equation(1)describes diffusion or heat transfer and no longer wave phenomena.Nevertheless，the same numerical boundary element framework (see section 2.2)can still be applied.

Figure 11.(a)Convergence and accuracy study for a standard scattering hard sound sphere at ka=1.01π(very close to the internal resonance frequency of ka=π).The incoming wave travels from bottom to top and the observation point is taken at r=1.2a above the sphere.The relative error between the numerical results by equation(7)and the analytical solution[49]reduces from 10-3 for 1000 nodes to about 10-5 for 8000 nodes for the scattered field.The computational time increases from a few seconds to 250 s.(b)Comparison between the analytical solution and the result with the non-singular Burton–Miller BEM of equation(54)for the case sketched in the inset of(a)when ka is swept from 1 to 30 with a step of Δka=0.01(note only 1 in every 25 points is plotted).The sphere surface is represented by 7842 nodes connected by 3920 quadratic triangular elements.Good agreement is found and no fictitious frequencies appear.

6.Discussion

All the above classic problems can，in principle，be solved just based on one or a set of Helmholtz equations.The following points are worth noting concerning the non-singular boundary element framework of section 2.2.

· 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[41]elements.Thus，anyone with a Helmholtz boundary element solver for sound waves，can in principle extend this to electromagnetic scattering simulations.Moreover，our method is stable for electromagnetic problems in the long wavelength limit.

· When solving vector wave problems，such as in electromagnetics，iterative solvers can also be used only requiring one N×N matrix to be solved，although careful attention needs to be paid to the convergence behavior.

· In figure 11(a)，the convergence behavior of the solution of the Helmholtz equation is given for scattering on a hard sphere for ka=1.01π.The solution at a sample point at r=1.2a (see the inset of the figure) is compared to the theoretical solution and the error is plotted as a function of node number.The relative error is about 10-3for 1000 nodes and goes down to 10-5for about 8000 nodes.The required CPU time is also indicated in the same graph.

· Concerning problems in the time domain: a timedependent 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 sound waves and for electromagnetic pulses in Klaseboer et al [16，21].

When using a BEM to numerically solve the Helmholtz equation for external problems，one challenge is that nonphysical solutions will show at certain discrete frequencies[42，43].This numerical issue is known as the occurrence of spurious solutions and occurs at fictitious frequencies.The two most popular methods to deal with fictitious frequencies are the CHIEF method proposed by Schenck [44] and the Burton–Miller method [45].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 which involves integration in parts and a principle value integration.Despite the fact that the Burton–Miller formulation results in hyper singular integrals，it is still possible to fully desingularize all the integrals involved.Since we believe this fully non-singular Burton–Miller formulation has not been shown before，the derivation is given in detail in appendix D.In figure 11(b)the scattering from a rigid sphere is shown where the desingulared Burton–Miller framework of equation(54)is used for ka=0–30.As expected，no fictitious frequencies appear here.

7.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 nonsingular 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 BEM 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 BEM 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 BEM

The expressions for ∂Gg/∂n and ∂G0/∂n are (with r=|x-x0|):

The terms in between brackets in equation (25) can thus be written with equation (26) as:

Here [ikr-1]eikr+1=2ikr+o(k2r2) and the term(x-x0)·n(x)behaves as 1/r2since in the limit of x going to x0，(x-x0)and n are perpendicular.Thus equation(27)does not contain any singular term.The term with ∇φ·(x-x0)in equation(25)goes as r and thus cancels out the 1/r behavior of ∂Gk/∂n.Similar for ∇g·(x-x0).

For the third term in equation (7) with，as long as ∂g/∂n approaches zero in a linear manner，it will cancel out the 1/r singularity of G0.

A similar analysis can be performed for the right-hand side of equation(7).For the first two terms on the right-hand side，we can write with again two similar Taylor expansions:∂φ(x)/∂n ≈∂φ(x0)/∂n+∇[∂φ/∂n]·(x-x0) and ∂f(x)/∂n ≈∂f(x0)/∂n+∇[∂f/∂n]·(x-x0)，and then

where we have used ∂f(x0)/∂n=1.The term with[Gk-G0]is regular，and the terms with the gradients contain both (x-x0)，which thus cancel out the 1/r singularity from Gkand G0，respectively.

Since f approaches zero as x →x0，it cancels out the third term with ∂G0/∂n on the right-hand side of equation (7).

In the Burton–Miller implementation of equation(51)(to be shown in appendix D) it is claimed that

when x →x0.The proof of this statement will be provided here.Due care has to be taken with the derivative ∂/∂n0，which now acts on x0and not on x，creating an additional minus sign.Applying this to equation (26) we get:

Using a Taylor expansion eikr≈1+ikr-k2r2/2，we get:

As explained before，when x →x0，term n0·(x-x0)and term n·(x-x0) both converge as r2and thus the last part of equation(30)is non-singular.The first term on the right hand side gives the required k2/(2r) in equation (51)，since n0·n →1.

Appendix B.Remarks on integrals at infinity and Sommerfeld condition

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 [46，47]and applies to the conventional BEM as given in equation(2)where the two integrals at infinity are:

Suppose we draw a very large sphere with radius R and consider this as ‘infinity’ then Gk=eikR/R and ∂Gk/∂n=eikRn·(x-x0)(ikR-1)/R3.Since n·(x-x0)→R and dS=4πR2，then:

The term with-1 vanishes as R →∞since φ decays at least as fast as 1/R.Then:

which can be shown to be zero with the help of the Sommerfeld radiation condition [46，47] expressing the fact that only outgoing waves are allowed (and realizing that at ∞，∂φ/∂n=∂φ/∂r) as

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

The second place where due care has to be taken with integrals at infinity is in section 2.2，where the desingularized BEM 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 equation (7):

With the particular choice of g=1 from equation(8)，it can be seen that the second integral is zero immediately(∂g/∂n=0).The first integral will become，with ∂G0/∂n=(x-x0)·n/r3→R/R3(as remarked earlier n·(x-x0)→R) and dS=4πR2，

Appendix C.Derivation for the dielectric case

We need to expand n·∂E/∂t1and n·∂E/∂t2in equation(20).First we write again E=Enn+Et1t1+Et2t2.Since ∂n/∂t1=-κ1t1and ∂t1/∂t1=κ1n and ∂t2/∂t1=0，we obtain n·∂E/∂t1=n·n∂En/∂t1-κ1n·t1En+n·t1∂Et1/∂t1+κ1n·nEt1+n·t2∂Et2/∂t1+0=∂En/∂t1+κ1Et1.We can then write:

where κ1is the curvature along the t1direction and κ2is the curvature along the t2direction.We assumed that μin=μout，in equation (20) and using the tangential boundary conditions for the tangential magnetic field，equation (19b)，results in

Thus equation (38) becomes

Due to the tangential continuity of the electric field in equation(19b)，the terms with κ1and κ2disappear in the above equation.Further using the normal conditions on the electric field in equation(19a)to replaceinto equation (39) gives:

Thus the two tangential components of the vectorare now known in terms of scattered and incoming electric fields.Applying the divergence-free condition of the scattered and total external field using equation (14) gives:

Subtracting these two equations and using the tangential boundary conditions in equations (19b) again，will eliminate the terms withandAlso using the normal conditions on the electric field to replace

Now we have all the components of the normal derivative of the transmitted electric field，both in Cartesian as in tangential and normal decomposition.Let us write

To reconstruct the Cartesian components，say the x-component of，we can use:

again with nx=(n·ex)，t1x=(t1·ex) and t2x=(t2·ex).If we replace the terms in brackets in equation (44) with equation (42) and equation (40) then:

We express the transmitted field in the x-directioninto normal and tangential components of the scattered and incoming field with the boundary conditions as:

with α=x，y，or z and

The first row of the matrix system just states(and similar for rows 2 and 3 for y and z)，but expressed in terms of normal and tangential components.

Appendix D.Desingularization of the Burton–Miller method

We demonstrate here a Burton–Miller boundary integral method for the Helmholtz equation without any singularities.To the best of our knowledge，a full desingularization of the Burton–Miller method has not appeared in literature yet，although hypersingular integrals have been studied extensively in the past [48].

The term on the right-hand side of equation(48)is hypersingular as x →x0.Performing the same operation on equation (4)，we have

If we choose ψ(x)=φ(x0)in the domain which indicates∇ψ(x)=0，introduce it into equation (49)，and subtract the resulting equation from equation (48) we have

The integrands in equation (50) are all weakly singular since (see appendix A)

To analytically remove the remaining singularities and the terms associated with the solid angle at x0in equation(50)，we follow the same procedure illustrated in section 2.2 to set up the following two functions that satisfy the Laplace equation in the domain where φ is valid:

Subtracting the conventional boundary integral equations corresponding to equation (52) from equation (50)，and considering the integrals over the surface at ∞for external problems，we get

The integrands in equation (53) are now all regular and this equation is now fully desingularized.

The Burton–Miller idea is to eliminate the spurious solutions at fictitious frequencies by combining equation(53)and equation (7) (multiplied by an imaginary parameter iβ)since the fictitious frequencies in equations (53) and (7) are always different from each other.To balance equations (53)and (7) in terms of physical dimensions (dimensional homogeneous in length)，the parameter β should be related to a characteristic length of the problem under consideration.In our simulations，shown in figure 11(b)，we used min [0 .5a，.The obvious choices are the size of the scattering object or the inverse of the wave number k.If we introduce equation (8) into (7)，and combine the resulting formulation with equation (53) according to the Burton–Miller idea，we obtain the non-singular Burton–Miller boundary integral method for the Helmholtz equation as

It is worth noting that in equation(54)，the integrals over the surface at ∞for external problems have been included，resulting in the terms with 4π.For internal problems，those terms disappear.Our goal to have a non-singular version of the Burton–Miller boundary integral method，as shown in equation (54) is now achieved and the fictitious frequencies are eliminated.