A new approach to vector scattering: the 3s boundary source method

This paper describes a novel Boundary Source Method (BSM) applied to the vector calculation of electromagnetic fields from a surface defined by the interface between homogenous, isotropic media. In this way, the reflected and transmitted fields are represented as an expansion of the electric fields generated by a basis of orthogonal electric and magnetic dipole sources that are tangential to, and evenly distributed over the surface of interest. The dipole moments required to generate these fields are then calculated according to the extinction theorem of Ewald and Oseen applied at control points situated at either side of the boundary. It is shown that the sources are essentially vector-equivalent Huygens’ wavelets applied at discrete points at the boundary and special attention is given to their placement and the corresponding placement of control points according to the Nyquist sampling criteria. The central result of this paper is that the extinction theorem should be applied at control points situated at a distance d = 3s (where s is the separation of the sources) and consequently we refer to the method as 3sBSM. The method is applied to reflection at a plane dielectric surface and a spherical dielectric sphere and good agreement is demonstrated in comparison with the Fresnel equations and Mie series expansion respectively (even at resonance). We conclude that 3sBSM provides an accurate solution to electromagnetic scattering from a bandlimited surface and efficiently avoids the singular surface integrals and special basis functions proposed by others. Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal


Introduction
Traditionally, the theory of electromagnetic waves has been applied most frequently in the field of electrical engineering, where it is central to the understanding of antennas and the design of radar imaging systems. With the growth of low-cost civil radar applications, wide-band, wireless communications and the continued miniaturization of information technology, rigorous modelling of electromagnetic scattering in the built environment is now of increasing interest. At optical frequencies, the application of rigorous scattering models is also gaining importance. To realize the considerable potential of phonic-bandgap and integrated fiber-optic devices our knowledge of electromagnetic scattering gained at radio frequencies must be transferred (through suitable computational tools) to the significantly higher frequencies and shorter wavelengths of the optical spectrum.
In many optical applications, the characteristic dimensions of structures of interest (e.g. lens apertures) are often several orders of magnitude greater than the wavelength. In these cases, the design and function of optical systems are best understood using the simplified models of physical optics [1]. From the early days of optics, ray tracing has been used to describe the propagation of light reflected or refracted by slowly varying surfaces and lies at the heart of the modern computational tools design packages that are routinely used to design complex lens systems [2]. More recently Monte Carlo methods have been used to model random scattering from rough surfaces and turbid media and are central to photorealistic rendering of computer-generated images [3]. These are powerful computational tools, however, they typically neglect diffraction effects and consequently are not (without significant modification) suited to the study of coherent scattering phenomena.
Huygens' principle provides us with a useful understanding of diffraction phenomena by explaining the propagation of light from a boundary surface as the sum of appropriately phased wavelets [4]. With the introduction of the obliquity factor, Fresnel showed that Huygens' principle follows directly from the scalar wave equation when expressed in the integral form [4]. Under the assumption that multiple scattering is negligible, the Huygens-Fresnel principle allows us to represent coherent scattering as a simple convolution operation and the optical output of measuring instruments to be written as linear filtering operations applied to a "foil" representation of the surface form [5,6]. In this way, we have a useful and conceptionally simple means to characterize the performance of optical surface profiling instruments including coherence scanning interferometers [7] and focus variation microscopes [8] when in normal use. It is noted that in this context normal use suggests an application to surfaces that are slowly varying on the scale of a wavelength (i.e. within the regime of physical optics) and with slope angles such that polarization effects and multiple scattering can be neglected. In practice, however, the effects of multiple scattering (e.g. from edges and v-grooves) are known to confound optical profilometers [9]. Making use of a-priori knowledge and inverse scattering models, multiple scattering can sometimes reveal features (e.g. undercuts) that otherwise would be invisible [10]. Furthermore, there is growing evidence that condition monitoring of machine tools and additive manufacturing processes is possible using straightforward measurements of scattering distribution and artificial intelligence (AI) algorithms [11].
The motivation of the work described in this paper was to create a computational tool that provides a rigorous solution to 3D scattering from an arbitrary interface between homogenous media for use in inverse scattering problems and the creation of synthetic training sets for AI condition monitoring applications. The method that we describe, properly accounts for polarization and the resonant behavior that is known to be problematic in coherent scattering. Using modest desktop computing it can be applied directly to irregular objects with a characteristic dimension of up to a few tens of wavelengths, open surfaces of equivalent area, or to arbitrarily small objects (with appropriate scaling). Moreover, the method is conceptionally simple and provides further insight into electromagnetic scattering from sub-wavelength surface features.

Background
In general terms the propagation of electromagnetic waves though in inhomogeneous media is governed by Maxwell's equations and these can be solved numerically according to well-defined illumination conditions. For the case of 3D media, Finite Element Methods (FEM) provide an efficient solution for monochromatic scattering behavior [12] while related Finite Difference Time-Domain (FDTD) analysis is typically applied to model wide-band or partially coherent illumination [13]. For the case of inhomogeneous and near-planar surface structures such as those forming semiconductor devices, Rigorous Coupled Wave Analysis (RCWA) provides an elegant solution that is partially expressed in the (spatial) frequency domain. RCWA is a particularly efficient way to analyse the performance of periodic structures such as blazed and coated diffraction gratings [14].
Electromagnetic scattering at the surface defined by the interface between homogenous media allows further simplification. In this case, the field at any point in a given space can be calculated from the field component on any surface that bounds that space according to the Kirchhoff's diffraction integral [15]. This approach can be formulated in several different ways. The Stratton-Chu formulation of the scattered field is defined in terms of all 6 components of the electric and magnetic fields the boundary surface [16]. It is noted, however, that Maxwell's equations relate these components and formulations due to Kottler [17] and later Franz [18] require only the 4 tangential components of these fields. In this way, if the illuminating field is known and appropriate boundary conditions applied, the surface field components can be calculated.
In computational terms the process can be treated as a Boundary Element Method (BEM) in which elemental areas act as finite sources and radiate into the volume of interest [19]. This approach presents some significant practical problems, however. First, the boundary elements that define typical surfaces have different shapes and areas and have varying radiation patterns. Second, the calculation of these radiation patterns is non-trivial as it requires integration of high-order singular functions. The Method of Moments (MoM) provides one solution to this problem [20] but requires significant computational overhead. The Method of Auxiliary Sources (MAS) [21] obviates the need for integration by expanding the fields in terms of a finite number of sources located on either side of the boundary. MAS avoids the singularities, but its accuracy is found to be sensitive to source position (this is discussed further in Section 3.3). Finally, we note that the Kirchhoff diffraction integral applies only to closed surfaces and many of the methods stated cannot be applied to open boundaries. With appropriate illumination conditions, however, these restrictions can be overcome (this is discussed further in Section 4.1).
The following section describes a new formulation of a Boundary Source Method (BSM) based on the Franz formulae of surface scattering. The method is similar in concept to the MAS [21] method, however, sources are now applied at the boundary and the extinction theorem (effectively an alternative boundary condition) is applied at a finite distance above and below the boundary. For the case of spherical particles and planar surfaces, we show that this method provides solutions that are almost identical to the exact solutions provided by Mie scattering and Fresnel equations respectively. Furthermore, we show that the method properly account for resonances.

Theory
In the following section we first discuss the electric field generated by the Franz surface integrals and relate these equations to an expansion of the electric field in terms of basis functions representing electric and magnetic dipoles. Calculation of the scattered electric field that is observed when a given object surface is illuminated by a known, monochromatic electric field is then discussed by applying appropriate field constraints. The relationship between this approach and the Huygens-Fresnel principle is then explored and the proper positioning of sources and constraints is discussed. Finally, the electromagnetic scattering problem is formulated as a straightforward matrix inversion.

Scattered field as a dipole expansion
Let us consider scattering from an object medium A within an ambient medium B as shown in Fig. 1. Let A and B be linear isotropic and homogenous media characterized by electric permittivities ε A and ε B , and magnetic permeabilities µ A and µ B , respectively. Definingn as the unit outward surface normal, according to the Franz formulation [18] a monochromatic electric field of angular frequency, ω, propagating, in medium A is given by, where ∇× denotes the curl of the vector field with respect to the space variable r; and G A (r , r) = exp(−i2πk A |r−r |) 4π |r−r | and k A = 1/λ A are the Green's function and wave number in medium A. It is also noted that the first term in Eq. (1) can be identified as the combined field due to magnetic dipoles that are tangential to the surface and have a strength (per unit area) that is proportional to the amplitude of the tangential component of the electric field. Similarly, the second term can be identified as electric dipoles that are also tangential to the surface and have strength (per unit area) that is proportional to the amplitude of the tangential component of the magnetic field. Noting that the ambient field, E B (r), is the sum of the incident field, E r (r), and the field scattered from the surface, we can write the field in medium B, where the Green's function and the wavenumber are appropriate to medium B and the terms can be interpreted as magnetic dipoles and electric dipole sources as discussed above. It is also noted that it is possible to write the magnetic field in each medium, H A (r) and H B (r) as similar superposition integrals [18], however, these fields are completely defined by the electric fields E A (r) and E B (r) and for simplicity are omitted here. Finally, we can simplify these formulae since the tangential components of the electric field and magnetic field are continuous across the boundary [22] such that,n (4) For computational purposes, it is convenient to combine the dipoles in each elemental area into discrete sources that are tangential to the surface and have magnetic and electric dipole moments where N is the number of source points. In this way Eqs. (5) and (6) define a multipole expansion of the electric fields in terms of 4N complex variables that define the phase and amplitude of each of the 4 electric and magnetic dipoles applied at each source point.
For a given incident field we can calculate the dipole moments required to satisfy appropriate boundary conditions. In the following, we make use of the extinction theorem of Ewald and Oseen [23] which can be written as the following constraints, in the ambient medium B; and in the object medium A. Here it is important to note the reversal of the Green's functions such that the exterior Green's function is used for the interior constraint and vice-versa. In order to solve Eqs. (7) and (8) we need to impose at least 4N constraints to the components of the field at the appropriate interior and exterior points as will be discussed in detail in Section 3.3. Before doing so, however, it is instructive to consider the electric field due to co-located orthogonal magnetic and electric dipoles as follows.

Vector Huygens' wavelets
Let us consider the field in medium A due to a magnetic dipole with moment m = m o k and an electric dipole of moment p = p o j co-located at the origin where i, j, k are the unit vectors of the coordinate system. If we apply the constraint of Eq. (7) at the position r = −Di where It is noted that by applying this constraint we have constructed a dipole source pair that does not radiate in the −i direction. The radiation pattern, E H (r), of this source pair is given by and is illustrated in Fig. 2. In this figure the red, blue and black arrows show the direction of the electric field, magnetic field and Poynting vector respectively while the shading corresponds to the magnitude of the Poynting Vector. The radiation pattern is strongest in the i direction and weakest (tends to zero as distance increases) in -i direction. It is interesting to note that if A 0 (r), is the amplitude at a distance r in the positive i direction, further analysis shows that the amplitude distribution of the field is given by, A(r, θ) = A 0 (r)(1 + cosθ)/2, where cosθ is the direction cosine measured from the x-axis. The term (1 + cosθ)/2 can be recognized as the well-known "obliquity factor" of the Huygens-Fresnel principle that results from a similar scalar analysis [24]. For these reasons we refer to such a forward propagating dipole source pair as a vector Huygens' wavelet and we note a similar "hypothetical vector Huygens' secondary point source" is described by Marathay [25].
From this analysis, the extinction condition of Eq. (7) requires that the tangential dipole source pairs must correspond to vector Huygens' sources that only radiate into medium A according to the interior Green's function G A (r). The extinction condition of Eq. (8) requires that the same tangential dipole source pairs exactly cancel the illumination field within medium A according to exterior Green's function G B (r). We are now in a position to consider the placement of control points where these conditions are imposed.

Optimised control point location
It can be concluded from the previous discussion that it is possible to generate electromagnetic fields internal and external to a closed body using tangential electric and magnetic dipole sources. Furthermore, for a given illumination field it is possible to calculate the dipole moments such that they satisfy the extinction theorem and thereby the transmitted and scattered field components can be calculated. It is well known, however, that dipole sources are mathematically intractable as they display third order singularity that can only be removed through volume integration [26]. The discrete dipole sources proposed in the previous section result in a rapidly changing field that is characterized by non-physical (infinitely) high spatial frequencies at the boundary. Propagation is known to be equivalent to linear filtering or smoothing, however, and after a short distance the high spatial frequency content is attenuated. By carefully considering the filtering, and by applying the constraints of Eqs. (7) and (8) at control points that are appropriate to the source spacing, the method of solution can be greatly simplified relative to conventional BEM.
As described in Section 3.1 we propose to generate internal and external electromagnetic fields using a basis of tangential electric and magnetic dipole sources that are placed at N source points on the surface of interest. In total, the basis is defined by 4N complex variables and in order to find these variables we require at least 4N constraints. In this work we constrain the tangential components of the electric field at 2N control points situated a small distance, d, above and below the boundary as shown in Fig. 3(a). The control points effectively sample the field generated by the sources on surfaces above and below the boundary and it is important that the field is properly sampled according to Nyquist-Shannon sampling theory [24].
To illustrate the process, let us consider a forward propagating plane wave generated by equal amplitude vector Huygens' sources placed on a regular grid of spacing, s, as shown in Fig. 3(b). In this case, the vector Huygens' source (as presented in Section 3.2) can be regarded as the convolution kernel that effectively smooths the granularity in the field imposed by the grid. We place a similar grid of control points in a parallel plane separated from the source plane. We require the separation distance, d, that is required to assure proper sampling. Figure 4(a) shows the spectrum of the spatial frequencies that are observed in the E x component of the electric field due to a vector Huygens' source at a wavelength λ = 1 [A.U.], at a distance, d = s = λ/100. The plot shown in Fig. 4(b) are sections through the distribution in the k x and k y directions and the green lines show the maximum spatial frequency according to the Nyquist limit defined by the control point spacing s (k x , k y < 1/2s). The Figs. 4(c) and 4(d) are similar to Figs. 4(a) and 4(b) but for a distance d = 3s = 3λ/100. At a distance d = 3s, the spatial bandwidth of the E x component of vector Huygens' wavelet predominantly falls within the bounds of the Nyquist limit and the effects of aliasing are therefore removed. Although this conclusion has been drawn for a specific component and wavelength, we note here that the spectra for other components are similar and a similar "rule of thumb" is drawn for all components for λ /1000 < s < λ /10. Finally, we note that if we increase the distance such that d >> 3s then several local sources will have a similar influence on a given control point and the system behaves in a similar manner to the underdetermined case where there are more sources than constraints, no unique solution exists and more complex and less efficient matrix procedures (such as pseudo-inverse) must be exploited. In the following section, the matrix solution of Eqs. (7) and (8) is considered in more detail.

Matrix formulation
To apply the theory of Section 3.1 it is necessary to write Eqs. just below and just above the boundary surface such that, Accordingly, Eqs. (7) and (8) can be written as the matrix equation, where S and C are vectors that represent the sources and constraints and W is a matrix that defines the influence of each source to each constraint such that, where where In this way the top row sub-matrices in W relate the sources to theû i component of the electric field at the control points in medium A using the Green's function for medium B such that, In a similar manner, the sub-matrices in the second row relate the sources to thev j component of the electric field at the control points in medium A using the Green's function for medium B. The third and fourth rows relate theû j andv j components of the electric field at the control points in medium B using the Green's function for medium A respectively, such that the fourth row is, We note that if the medium A is a perfect conductor then the field expansion requires only tangential electric dipoles and their strength can be found by cancelling the illumination field at every internal point such that, and Finally, we note that W is a square matrix with either 16N 2 elements (for the case of dielectric/dielectric interface) or 4N 2 elements (for the case of a dielectric/perfect conductor interface). If d = 3s and s<λ min /2 (where λ min is smallest of λ A or λ B ) we have found that this matrix is well conditioned and can always be inverted efficiently to find the source coefficients, such that, S = W −1 C, as we show in the following examples.

Application of 3sBSM
In this section, we apply the 3sBSM to calculate the scattering from planar and spherical boundaries comparing the results are compared with the well-known analytic formulae of Fresnel and Mie.

Scattering from a plane dielectric/dielectric interface
The aim of this section is to demonstrate that the 3sBSM can replicate the well-known Fresnel equations that define amplitude reflection and transmission coefficients at a plane surface between two dielectrics. The Fresnel equations are strictly only valid, however, for the case of infinite plane waves at an infinite plane (i.e. an open boundary) and the comparison presents a problem, as the theory developed in Section 3. is derived from the closed boundary integrals of Franz. However, we note here that the 3sBSM is applicable to open boundaries in the sense that each of the dipole basis functions obeys Maxwell's equations and secondly, regularly spaced tangential electric and magnetic dipoles form a complete basis (note this is not the case when the Stratton-Chu formulae are applied to an open boundary [27]). We can therefore apply the method with confidence if i) the expression used to describe the illumination, E r (r), obeys Maxwell's equations and ii) a sufficient surface area is covered by the basis such that the contribution due to omitted dipoles can be assumed negligible.
In this example, we consider reflection from an air/glass interface over a surface area of 20λ B × 20λ B with illumination by a Gaussian beam with a waist radius 2λ B (@1/e) at the surface according to the rigorous 3D vector Gaussian beam formulation [28]]. The refractive indexes of the upper medium B and lower media A are n B = (ε B µ B )/(ε 0 µ 0 ) = 1 and n A = (ε A µ A )/(ε 0 µ 0 ) = 1.5, where ε 0 and µ 0 are the free space electric permittivity and magnetic permeability. The angle of incidence is nominally the Brewster (polarising) angle (approx. 57 degrees). The surface is sampled such that s = λ B /5. Figure 5(a) shows extinction in the plane of incidence. The E x , E y and E z components of the sum of the incident field, E r (r), and that radiated from the surface dipoles at k B = 1/λ B are shown in the top row showing complete extinction in the lower medium A and reflection of the E y component in accordance with incidence at the Brewster angle. The lower row shows E x , E y and E z components of the electric field that is radiated from the surface dipoles at k A = 1/λ A illustrating extinction in the upper medium B and the transmitted field in the lower medium A. Figures 5(b) and 5(c) show a comparison of the reflection and transmission coefficients for different plane-wave components using the 3sBSM and Fresnel formulae. It is noted here that due to diffraction, the incident Gaussian beam can be decomposed into plane polarized, plane-wave components cover over a range of incident angles from about 49-66 degrees. The reflected and transmitted fields can be similarly decomposed. The phase and amplitude of the reflection coefficients in Fig. 5(b) and transmission coefficients Fig. 5(c) are in good agreement with those predicted by the Fresnel equations.
In Fig. 5(a) the surface is highlighted by the yellow pixels, while the ± 3s zone between the control points is highlighted by the red pixels. We note here that the field within the ± 3s zone is not actually an accurate representation of a physical electric field but is merely one solution to an infinite number of electric fields that satisfy the extinction theorem in the remaining domain as will be discussed further in Section 5.

Scattering from a dielectric sphere
The work in this section provides a comparison between 3sBSM and Mie formulas for the case of scattering of plane waves from dielectric spheres. In principle this closed boundary is an easier problem to solve than open boundary in the previous section, however, closed surfaces resonate at certain frequencies which can lead to inaccuracy [29].
The matrix equations of Section 3.4 have been solved for the case of a dielectric sphere of refractive index n A = (ε A µ A )/(ε 0 µ 0 ) = 1.5 in a vacuum n B = (ε B µ B )/(ε 0 µ 0 ) = 1 with radius r = λ B and s = λ B /10. The incident light is assumed to be a plane wave which propagates in the k direction (i.e. along z-axis). Figure 6(a) demonstrates the extinction of the sum of the scattered field and incident field inside of the sphere and the extinction of the transmitted field outside the sphere while the surface position and ± 3s zones are highlighted by yellow and red pixels as before. In Fig. 6(b) the comparison with Mie series is shown at a radius, r = 1.4λ B , is presented.
Although a very good agreement with Mie scattering is observed, it is well known that resonant cavities present additional complexity and difficult to model properly [30], [29]. When resonances occur the Scattering cross section C scat of spherical particles increases sharply [31] such that, where a n , b n are the coefficients in the expansion of the scattered fields in vector spherical harmonics (for details, please see [31] Chapter 4.4). The scattering cross section C scat for a dielectric sphere with refractive index n A = 2, n B = 1 and radius r = λ B is shown as a function of the wavelength [A.U.] of the incident plane wave with field components in Fig. 7(a). Resonant peaks are identified A-D. The 3sBSM has been used to model the same system and a comparison of the scattered fields with Mie predictions are shown in the lower half of this figure. The incident light is assumed to be a plane wave which propagates in the k direction (i.e. along z-axis). The plots shown in Fig. 7. show an excellent agreement between the predictions of 3sBSM and those of Mie series. We also note that scattered fields calculated for a very small sphere with radius r = λ B /100, s = r/10 = λ B /1000, n B = 1, n A = 1.5 also coincide nicely with Mie series. The reasons for this are discussed further in the following section.

Discussion and conclusions
In this paper we have presented a straightforward method to compute the electric fields that are scattered and transmitted by the surface at the interface between homogenous media. The electric fields are expressed as an expansion of sources in the form of electric and magnetic dipole pairs that are i) tangential to the surface and ii) spaced at approximately even intervals. The relative phase and amplitude of the dipole sources are found by applying the extinction theorem to the tangential components of the electric field at suitable control points in each medium. We have found that control points are optimally positioned at a perpendicular distance of 3s on either side of the interface where s is the nominal source spacing the source, resulting in a robust and efficient solution based on regular matrix inversion and for this reason we refer to the method as 3sBSM.
The choice of sample spacing, s, effectively defines resolution of the technique. It might be expected from the results presented that the 3sBSM provides an exact calculation of any smooth surface generated by a point cloud of this resolution (e.g. generated using a cubic-spline interpolation). It should be remembered, however, that it is the electric field components at the control points that are constrained, and these have propagated a minimum distance, 3s, from the sources. As mentioned in Section 4.1 this electric field could in principle be generated by an infinite number of sources positioned anywhere in the region between surfaces defined by the inner and outer control points. It is more reasonable therefore to expect the method to correctly represent a smooth surface that lies somewhere between these bounds, or equivalently a surface form that is defined by a bandlimited function and is adequately sampled (by the Nyquist criteria) at a resolution of 3s. We also note here that more advanced sampling strategies such as those described in the context of MAS might also be possible [32,33]. Finally, it is worth stating that in order to sample the propagating electric fields correctly, the maximum sample spacing is half of the wavelength in the medium of greatest refractive index. There is no apparent restriction to the minimum sample spacing.
With these restrictions the 3sBSM provides a straightforward and easy to implement method to solve electromagnetic scattering problems. The results in this paper show a good agreement between 3sBSM and analytical results obtained from Mie series expansion and Fresnel formulae.  6. Scattering from a sphere with radius r=λ B , n B = 1, n A = 1.5, at discretization s = λ B /10; (a) extinction of the sum of the scattered field and incident field inside of the sphere and the extinction of the transmitted field outside the sphere, (b) comparison of the scattered fields with Mie series. Fields are plotted in xz plane, outside the sphere at radius a = 1.4λ B as a function of the angle θ between radius vector of the point and the x-axis (such that in forward scatter θ = π/2). Furthermore, the method appears to provide accurate results for resonant cases where problems have been reported elsewhere. We believe this is because i) the basis functions used to generate the electric fields at the control points are complete and ii) the fields at the control points are correctly sampled according to Nyquist criteria. With these two conditions we believe that the 3sBSM provides a single, unique solution to the electric fields defined by the extinction theorem.

Disclosures
The authors declare no conflicts of interest.