Field-only integral equation method for time domain scattering of electromagnetic pulses
Evert Klaseboer Qiang Sun Derek Y. C. Chan

Sales Forecast in E-commerce using 
 Convolutional Neural Network
    
Kui Zhao

       College of Computer Science

       Zhejiang University

       Hangzhou, China

       zhaokui@zju.edu.cn
Can Wang

       College of Computer Science

       Zhejiang University

       Hangzhou, China

       wcan@zju.edu.cn

    Received: December 30, 2023/ Accepted: December 30, 2023
==========================================================================================================================================================================================================================================================

fancy

shortarticlesinglecolumn



§ INTRODUCTION


Current theoretical formulations and computational algorithms for time domain electromagnetic problems are well-developed for applications that range from transient spectroscopy <cit.> to telemetry involving high speed spacecrafts <cit.>.
Under conditions relevant to micro-photonics in which the interest in regimes of wavelengths can be small or large compared to the characteristic dimensions of the scatterers or the need to have accurate values of the phases and amplitudes of the field near boundaries, current methodologies encounter challenges.

This paper addresses the time domain electromagnetic problem of the scattering of electromagnetic pulses by working directly in terms of the components of the electric field using a non-singular surface integral formulation in the frequency domain. 
Fourier transform is then used to give the complete space-time behavior. 

This approach retains the advantages of reduction in spatial dimension of surface integral methods, and because of the non-singular nature of the surface integrals, it has the added ability to handle surface geometric intricacies that often arise in micro-photonics problems in which certain length scales in the problem many be small compared to the wavelength or where there is the need to obtain accurate results for field values on or near surfaces. 
The use of the Fourier transform to give the time evolution also avoids possible instabilities that can arise with march on time algorithms.
Examples of the space-time dependence of scattering of an incident wave pulse by conducting and dielectric scatterers are provided to illustrate the transition between wave and geometric optics. 
The structure of the field near the surfaces is used to demonstrate field focusing and multiple scattering effects.



§ BACKGROUND


One established approach to time domain electromagnetics is based on extending the surface integral formulations, the electric field integral equation (EFIE) or magnetic field integral equation (MFIE) method, to the time domain using a march on time method, see for example, <cit.>,  <cit.> and <cit.> for reviews. 
This widely used approach involves solving surface integral equations for the surface currents.
The electric and magnetic field are then obtained by post-processing the surface current values.
In present formulations, the surface integral equations that need to be solved for the induced surface current densities contain singular kernels that originate from the Green's function. This mathematical feature that does not have a physical basis means a loss of precision in the calculation of field values close to boundaries between different media. 
Furthermore, march on time algorithms can also lose precision as time progresses.

Another approach to time domain electromagnetics is the finite difference time domain (FDTD) method of Yee <cit.> where the space- and time-dependent Maxwell's equations are discretized in the 3D space variables and time stepping is used to track the space-time evolution from given initial conditions.
Although the FDTD algorithm is simple conceptually, there are a number of technical issues that require care in implementation. 
Convergence constraints impose limitations on the step sizes in time and space. 
Numerical dispersion effects associated with the relative orientation of the spatial grid and the direction of propagation can arise and unphysical reflections can occur at the boundaries (that conform to the stepwise nature of the grid) between regions of the different grid densities.
If the problem domain is infinite, an outer boundary needs to be constructed with suitable boundary conditions to satisfy the radiation condition at infinity so as to avoid unphysical reflections back into the solution domain <cit.>.

Earlier works on pulses have been based on extending the analytic Mie theory for scattering by a sphere <cit.> and as such are not readily applicable to consider problems involving general scatterers.





§ THEORY

The present theoretical development is motivated by the Mie theory <cit.> in which the solutions of Maxwell's equations in the frequency domain are constructed from 2 scalar Debye potentials that obey the Helmholtz equation. 
Although the Debye potentials can be found by solving boundary integral equations, they are not suitable for numerical implementation because the boundary conditions on the electric, E and magnetic, H fields involve the first and second derivatives of these scalar potentials.

The following is a summary of the field only, non-singular boundary integral formulation of the solution of Maxwell's equation in the frequency domain. Detailed development and numerical implementation issues have been given in the literature <cit.>.

For a vector field that satisfies wave equation

    ∇^2 E + k^2 E = 0

the divergence free condition on E can be written as

    2 ∇·E≡∇^2 (x·E) + k^2 (x·E) = 0

whereby the scalar function (x·E), where x is the position vector, also satisfies the scalar wave equation or the Helmholtz equation.

Therefore the solution of Maxwell's equation for the scattered electric field, E can be expressed as the simultaneous solution of 4 coupled scalar Helmholtz equations of the general form

    ∇^2 p_i(x) + k^2 p_i(x) = 0,       i = 1 .. 4
 
where p_i(x) denotes (x·E) or one of the 3 Cartesian components of E.
These 4 equations are coupled by the usual boundary conditions on the normal and tangential components of E at material boundaries.
The scalar function (x·E) is the result of the application of the angular momentum operator on one of the Debye potentials <cit.>. 

Writing the solution of Eq. <ref> as surface integral equations furnishes 3 relations between the 6 unknowns: E_α and ∂ E_α / ∂ n, (α =x,y,z), where ∂/∂ n ≡n·∇ and n is the outward unit normal of the surface, S of the solution domain. The boundary integral solution of Eq. <ref> for the quantity (x·E) provides one more relation between E_α and ∂ E_α / ∂ n since: ∂ (x·E) / ∂ n = n·E + x·∂E / ∂ n. The electromagnetic boundary conditions on the continuity of the tangential components of E provide the remaining 2 equations to determine E and ∂E / ∂ n completely. At the surface of a perfect electrical conductor, the tangential component of the incident field must be canceled by that of the scattered field.

The boundary integral solution of Eq. <ref> for the scattered field is based on Green's Second Identity that gives a relation between p_i(x) and its normal derivative ∂p_i/∂n at points x and x_0 on the boundary, S.  All singularities associated with the Green's function G ≡ G(x,x_0) = exp(ikr) / r, as r ≡ |x - x_0|→ 0 can be removed analytically by using the following non-singular formulation of the boundary integral equations <cit.>

    ∫_S^[p_i (x) - p_i (x_0) g(x) - ∂p_i(x_0)/∂n  f(x)] ∂G/∂n dS(x) =     
    ∫_S^G [ ∂p_i (x)/∂n - p_i (x_0) ∂g (x)/∂n - ∂p_i(x_0)/∂n∂f (x)/∂n ]dS(x).
 
The requirement on the functions f(x) and g(x) is that they satisfy the Helmholtz equation and the following conditions at the point x = x_0 on surface, S: <cit.>. 


    f(x_0)     = 0 
    n·∇ f(x_0)     = 1 
    
    g(x_0)     = 1 
    n·∇ g(x_0)     = 0



Therefore, if p_i (or ∂p_i/∂n) is given, then Eq. <ref> can be solved for ∂p_i/∂n (or p_i) in a straightforward manner. The reason is that for f(x) and g(x) that obey the above conditions, Eq. <ref>, the terms that multiply G and ∂ G/∂ n vanish at the same rate as the rate of divergence of G or ∂ G/∂ n as x→x_0 and consequently both integrals have non-singular integrands and can thus be evaluated accurately by quadrature, see <cit.> for details. 
Consequently, this approach affords higher numerical precision with fewer degrees of freedom and confers numerically robustness that enables, in particular, the accurate calculation of field values on or near boundaries without adverse numerical issues.
Note that the solid angle at x_0 that occurs in the traditional boundary integral equation has also been eliminated in Eq. <ref>. Suitable choices for f(x) and g(x) can be found in <cit.>.
Similar coupled equations also hold for the magnetic field H, see <cit.> for details. 
 
Furthermore, the present field only formulation is not affected by the so called zero frequency catastrophe as k → 0 that limits the accuracy of the familiar electric and magnetic field integral equations for the surface current density <cit.>.

As to be expected, the system of equations to be solved are simpler for perfect electrical conductor (PEC) scatterers <cit.> than for dielectric scatterers <cit.>, but the underlying physical concepts are the same. 





§ RESULTS


A demonstration of our field only approach for small and large dimensionless wave number, ka where a is a characteristic length scale of the problem, is now given. 
To this end we consider variations of the magnitude of the scattered field around:

a) a perfect electrical conductor (PEC) sphere of radius, a,

b) an axisymmetric bowl shaped PEC object obtained by rotating the parametric closed curve in the body coordinates, see Fig. <ref>

    (X, Z) = a ( 2 sinθ, βsin^2 θ + γ [cosθ - 1]),    0 ≤θ < 2 π,

about the Z-axis. 
The bowl rim has radius, 2a and the inner concave bottom of the bowl can be fitted to a parabola Z = X^2/(4f) where the focal length, f = R/2 is related to the radius of curvature, R = 4a/(2β-γ) at the inner center of the bowl. 
Here we choose β = 0.6 and γ = 0.5 that gives f = 2.86 a and R = 5.71 a as shown in Fig.<ref>. 

c) an axisymmetric thin nano-rod with length, L and maximum cross-section diameter, b whose surface is defined by rotating an analytic curve about the long axis <cit.>. The aspect ratio of the rod is L/b = 10. The cases of the nano-rod being a PEC and a dielectric are studied.

In Fig. <ref> we show the scattered field amplitudes by perfect electrical conductors (PEC): a sphere and a bowl, excited by an incident plane wave that propagates in the positive z-direction and with E polarized in the y-direction with wave number, k.
The results shown in the yz-plane demonstrate the transition from the electrostatic limit at small wave numbers, ka, to approaching geometric or physical optics as ka increases. 
Note also the high induced local field strengths associated with the corona effect that is evident around the high curvature rim of the bowl in the long wavelength (small ka) limit. 

 

To demonstrate time dependent effects, we consider the scattering of an incident plane wave pulse: E^inc = (0, E^inc_y, 0) that is polarized in the y-direction and propagates in the positive z-direction with the following form in a window of width, w 

    E^inc_y = {[      0,                   -4 π N_c ≤τ  < -2 π N_c; sin(τ) exp (-α |τ|),         -2 π N_c ≤τ≤ 2 π N_c;         0,                   2 π N_c < τ≤ 4 π N_c;                                                   ].

where τ≡ k_0 (z-vt). 
Examples of such pulse functions that comprise of 2N_c oscillatory cycles and their Fourier components calculated using a 128 point discrete Fourier transform are given in Fig. <ref> for different values of N_c and α <cit.>. 
Since the Fourier amplitudes are small at high wave numbers, using say the first 20 Fourier components is generally sufficient to give an accurate representation of the pulse.

For each of these Fourier components of the incident pulse we solve for E and (x·E) at the corresponding value of k and then the time domain behavior can be found by inverse discrete fast Fourier transform <cit.>.  
As in earlier work on the scattering of acoustic pulses <cit.>, we set the total window width of the incident pulse to be w = 20.1a, of which the central oscillatory portion has width w/2. 
So for a pulse with 2N_c oscillatory cycles over the width w/2, the wave number, k_0 in Eq. <ref> is determined by the relation (2N_c)λ_0 ≡ (2N_c) (2π/k_0) = w/2, giving k_0 a = 8π N_c/20.1. The results are presented in such a way, that the pulse just reaches the bottom of the scatterer at the first frame. To account for this, each frequency component of the incoming wave must be multiplied with a phase factor e^iβ with β=-k_ma[w/(4a)-1] and k_m the m^th wave number of the Fourier spectrum. 

In Fig. <ref> we show both the total and scattered E_y electric field amplitudes as well as the electric field vectors excited by the plane wave pulse given by Eq. <ref>. Results in the yz-plane around a PEC sphere are shown at selected time steps out of a total of 128 steps (1002 nodes and 500 quadratic elements were used). The field strengths are indicated by the color scale and their strengths and directions by arrows. These results have been validated against available analytic solutions <cit.>.







In Fig. <ref>, we show the magnitudes of the total field excited by a broader incident plane wave pulse given by Eq. <ref> with N_c = 1, α = 0.5 and hence k_0 a = 8π/20.1 in the neighborhood of a PEC bowl with 1002 nodes and 500 quadratic elements.
Here the high field strengths associated with the corona effect are more evident around the high curvature rim of the bowl although the maximum field amplitude is only about half that in Fig. <ref>.

In Fig. <ref> we show the total amplitudes and field vectors of E excited by the same incident plane wave pulse as in Fig. <ref> in the neighborhood of a PEC bowl with rim radius 2a as given in Eq. <ref>. 
The electric field strengths are indicated by the color scale and the field strengths and directions by arrows.
The focusing effects due to the curvature of the bowl are clearly evident with regions of high field amplitudes located near the position of the focal point discussed above. 
The focusing effect increases the field amplitude by up to about 50%. 
Complex structures in the vector E are evident in the enlarged sub-figures near the concave surface of the bowl.




To demonstrate the ability of the present field only formulation to handle multiple scattering effects in which spatial configurations of the scatterer can cause significant field enhancements, we consider the scattering by the combined effect of a PEC bowl that has a rim radius, 2a with a sphere of diameter a, that is located near the focal point of the bowl. The bowl and the sphere each have 362 nodes and 180 quadratic elements.
In Fig. <ref>, we show the space-time variation of the total field amplitudes and directions, for an incident pulse given by Eq. <ref> with N_c = 2, α = 0.1 and k_0 a = 16π/20.1 – the same as that in Fig. <ref>. Here we see that this configuration of a bowl with a sphere placed near its focus can enhance the local field amplitude by about a factor of 3, twice that achieved by a bowl on its own. The complex field structures in the region between the bowl and the sphere are shown in the enlarged sub-figures.



The scattering from an object with high aspect ratio is tackled next, both for a PEC and a dielectric object, a task which is rather challenging for conventional methods. 
In Fig. <ref> and Fig. <ref>, we compare the values of the E_y and E_z components of the total electrical field excited by the incident pulse given in Eq. <ref> by a PEC and by a dielectic nano-rod with refractive index n_p=3 <cit.>, the surrounding medium has unit refractive index. 
The y-component of the electric field, E_y is shown on the left and the z-component, E_z on the right at three time instances in both figures. 
In both cases, the long axis of the nano-rod is oriented at an angle of 45^∘ to the propagation direction of the incident pulse. 
The surface of the axial symmetric nano-rod is defined by rotating an analytic closed curve about the long axis <cit.>. 
The rod has length, L and maximum cross-section diameter, b at the middle of the rod, with L/b = 10. 
In both figures, the incident pulse is defined by Eq. <ref> with N_c=2 and α = 0.1. 
This corresponds to, k_0L = 32π/20.1 and k_0b = 3.2π/20.1 so that the length of the nano-rod is approximately the same wavelength, λ_0 ≡ 2π/k_0 of the incident pulse and the width of the nano-rod is about λ_0/10. 
As expected, the PEC nano-rod scatters much more strongly than its dielectric counterpart. 
This is evident by comparing the snapshots at the same time step in Fig. <ref> and Fig. <ref>. 
At frame 31, when nano-rod is near the center of the pulse, the scattering due to the PEC nano-rod is larger and consequently the perturbation of the total field is more severe. 
Towards the tail of the pulse at frame 53, the scattered field around the PEC is again much more pronounced than that of the dielectric nano-rod. After the pulse has passed, the scattered field from the dielectric nano-rod is essentially zero, whereas that due to the PEC is still creating a scattering effect which slowly dies out. 
This effect is even more obvious in the supplementary Visualizations 5 and 6). 
As a check of our numerical implementation, the scattering goes to zero smoothly, as expected, when the nano-rod is made transparent by taking the limit of its refractive index to be the same as that of the surrounding medium, n_p→ 1. 






§ CONCLUSIONS

We have leveraged on the direct nature and numerical stability of the field only formulation of the solution of Maxwell's equations in the frequency domain <cit.> to extend it to find time domain solutions using the Fourier transform. 
The field only formulation solves directly for the field components instead of first finding the surface current and then computing the field by post-processing. It is therefore particularly suitable for situations in which one is interested in the field amplitudes, phases and directions, especially near material boundaries - often desired in the design of microphotonic devices or in surface enhanced Raman spectroscopy.

In contrast, because of the singular nature of the surface integrals for the surface currents, the post processing step to find field components from the surface currents becomes unstable numerically as the field point becomes close to a boundary.
Although only results for the field in the yz-plane of symmetry are given for simpler visualization purposes, the complex 3D structure of the field around the scatterers are readily available. 

Although the finite difference time domain (FDTD) method also solves directly for the field components, it is necessary to discretise the 3D spatial domain as opposed to only to focus on 3D boundaries in surface integral formulations. This becomes especially beneficial for spatial regions with different characteristic length scales such as near the rim of the PEC bowl and around the sharp ends of the nano-rods in our examples. The 3D discretization problem can be challenging for conventional methods. Furthermore, when the natural orientation of the scatterer is at an angle to the direction of field propagation, care is needed in the construction of the 3D spatial mesh to avoid spurious dispersion effects <cit.>. 

The present non-singular boundary integral formulation <cit.> for the solutions of the Helmholtz equation for the field components, is unaffected by numerical instabilities at the low frequency regime or in the regime when the characteristic length scale of the problem is small compare to the wavelength. In this low frequency or long wavelength regime, the conventional boundary integral equation approach will exhibit numerically troublesome near-singular behavior sometimes referred to as the zero frequency catastrophe that limits the precision that can be attained in calculations that involve first having to find induced surface currents <cit.>. In contrast, the present approach is numerically robust for wavelengths that span the electrostatic to approaching the geometric optics regimes without the need to modify the algorithm.

The absence of singularities in the integral equations also means numerical instabilities do not arise in multiple scattering effects due to scatterers in close proximity or in the evaluation of field values close to material boundaries. And in the absence of singular kernels, it is easy to use quadratic surface elements with consistent quadratic interpolation of the continuous integrands <cit.>, as we have done so in our numerical examples, to evaluate the surface integrals. This means higher order precision can be obtained with fewer number of nodes compared to using planar surface elements 

For a fixed set of scatterers, the boundary integral equations only need to be solved once from which the space-time variation of the total field for different incident pulses can be found by inverse Fourier transform.

In summary, the present approach involves solving non-singular scalar surface integral equations directly for the electric field components on the boundary of scatterers. In contrast, the conventional surface integral equation approach involves solving integral equations with singular kernels for surface currents. Although there are well developed methods to deal with such issues in the conventional surface integral approach, the integral equations become numerically ill-conditioned when the characteristic length scales of the scatterers are small compared to the wavelength. Since the electric field is obtained from the surface current by post processing, the presence of singularities in the integral equations means that the precision of field values at or near surfaces is compromised. The present approach avoids such issues by having to solve 4 scalar Helmholtz equations by a non-singular integral equation formulation. Thus, in handling the time evolution by Fourier transform, the robustness of our method at all frequencies is a distinct advantage.

Compared to the finite difference time domain (FDTD) method, the ability to achieve a reduction in dimension by only having to solve the problem on the boundary of scatterers also removes the challenge of having to discretize a complex  3D domain when there are vastly different characteristic length scales and spatial directions in a problem. In the FDTD approach, changing incident field for the same set of scatterers, requires repeating the solution process with each new incident field. In contrast, with the present approach, once the frequency domain problem has been found for a fixed set of scatterers, the solution can be used to find the response to different incident fields by simply taking the inverse Fourier transform.  



 This work is supported by the Australian Research Council (ARC) through a Discovery Project Grant (DP170100376) to DYCC and a Discovery Early Career Researcher Award (DE150100169) to QS.






Main
 

