{'fraction_non_alphanumeric': 0.029935693695025502, 'fraction_numerical': 0.009793776332323157, 'mean_word_length': 4.4576442113755546, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A 3D Nonlinear Maxwell's Equations Solver Based On A Hybrid Numerical Method

In this paper we explore the possibility for solving the 3D Maxwell's equations in the presence of nonlinear and/or inhomogeneous material response. We propose using a hybrid approach which combines a bound- ary integral representation with a domain-based method. This hybrid approach has previously been successfully applied to 1D linear and non- linear transient wave scattering problems. The basic idea of the approach is to propagate the Maxwell's equations inside the scattering objects for- ward in time by using a domain-based method, while a boundary integral representation of the electromagnetic field is used to supply the domain- based method with the required surface values. Thus no grids outside the scattering objects are needed and this greatly reduces the computational cost and complexity.


Introduction
Boundary Element method (BEM), as a tool for solving scattering problems, has several attractive features. First and foremost, BEM is well suited to treating scattering problems in unbounded domains because the boundary integral equations are located on the surfaces of the scattering objects and thus one whole dimension is taken out of the problem. Also, the scattering objects are usually defined by sharp material boundaries and thus a domain-based method must seek to resolve the fast variation in the corresponding solutions generated by the boundaries. This problem has been more or less solved by using perfectly matched layers [1,2], but the solution comes with additional cost and complexity.
Secondly, the surface localization of boundary integral equations also means that the boundary discretization, which leads to the BEM equations, can be optimized with respect to the geometry of each surface separately when there are more than one scattering object, which there usually is. For domain-based methods like Finite Element method (FEM) [3][4][5] and Finite Difference Time Domain method (FDTD), such an optimization can only be achieved by using non-uniform or multiple grids tailored to the geometrical shape of the objects. These kinds of efforts have met with some degree of success, but it does add new layers of complexity.
Thirdly, BEM, is exceptionally well suited for modeling scattering problems where the sources are slow on relevant timescales. In this setting the boundary integral equations are derived in the spectral domain and the discretized equations defining BEM only needs to be solved for the small set of discrete frequencies that are required for an accurate representation of the time dependent source. For a domain-based method like FDTD, near-stationary sources are the worst possible case since FDTD is, as the name indicate, a time domain method and slow sources mean that the Maxwell's equations have to be solved for a long interval of time, which is costly in terms of computational resources.
Given all these attractive features of BEM for solving scattering problems in electromagnetics, it is somewhat surprising that in a popularity contest, FEM and FDTD beat BEM hands down [6]. The reason for this is that in addition to all its attractive features, BEM has some real drawbacks too.
Firstly, reformulating the scattering problem for a system of PDEs in terms of boundary integral equations, typically will involve classical but fairly intricate mathematical tools, as will the discretization of the resulting integral equations. In particular, one will need to content with singularities which appear in the limits that are always taken while deriving the boundary equations from the PDEs. Accurately representing these singularities in the ensuing discretization leading to BEM, is a major issue. FEM has some of the same issues, whereas FDTD is much simpler to implement using mathematical tools that are straight forward and known to all.
Secondly, BEM relies on Green's functions, and such functions can only be defined for linear systems of PDEs. In our opinion this is the major reason why BEM is significantly less popular than the major contenders FEM and FDTD. Nonlinearities are common in most areas of pure and applied science, and having to change your whole computational approach when nonlinearities are added to your model is a major nuisance. In some areas of application nonlinearities dominate and hence a computational approach based on Green's functions is out of the question, fluid mechanics is such an area. However in other areas one can frequently disregard nonlinearities and a computational approach based on Green's functions is feasible. Scattering of electromagnetic waves, where nonlinearities only come into play at very high field intensities, is such an area.
Thirdly, BEM, which is simplest to formulate in the spectral domain, is at it's best when we are looking at a near stationary situation with narrow band sources. For transient scattering with broad band sources, one needs to formulate BEM in the time domain, and this involves space-time Green's functions that in general tends to be more singular than the frequency domain ones. However, for the specific case of electromagnetic scattering an integral formulation has been derived [7] whose singularities are fairly weak. This integral formulation of electromagnetics is the foundation of all time dependent BEM formulations to date, and is also the basis for our approach. However, the progress in developing these time dependent BEM schemes has been slow due to several drawbacks.
The time-domain integral equations are retarded and this means that in order to compute the solutions at a certain time, one needs to retain the solutions for an interval of previous times that in some cases can be very large. This leads to a large memory requirement that needs to be met using parallel processing. In today's computational environment parallel implementation of time dependent BEM is fairly standard, but the possibly limited efficiency due to the problem of load balancing is something that always must be contended with.
However the major obstacle that has prevented this method from being widely applied for electromagnetic scattering is the occurrence of numerical instabilities. These instabilities, whose source is not fully understood, occur not at early times, but at later times, and have become known as the late time instability. Many efforts have been made and several techniques have been developed in order to improve the stabilities of BEM schemes for time dependent electromagnetics in the last several decades [8][9][10][11][12][13][14]. Our work has not been aimed at joining in or improving on any of the efforts pursued by these research groups. Our major aim has been to generalize boundary approaches to electromagnetic scattering for cases where the scattering objects has an inhomogeneous and/or nonlinear response. Most work in the area of time dependent BEM has been focused on metallic objects with linear response whose dimensions are large with respect to the wave length, antennas is a major example of the kind of structure one has been interested in. The kind of problem we have in mind is scattering from very small objects, from micron to nanometer scale, objects who might have engineered inhomogeneities in their structure and strong, also engineered, nonlinear response. At this scale, the standard simplifying assumption of disregarding the inside of the scattering objects and modeling them using surface charges and currents, is not applicable. The skin depth at this microscopic scale can easily be as large as the scattering objects themselves, which is a marked difference to what is true for macroscopic antenna theory. For this reason, and also because of the complex inner structure of these microscopic scattering objects, a different approach is needed. The traditional BEM can not be used here.
Our approach is based on the same integral representation of the electromagnetic field [7] as the traditional BEM, however, we use the integral representation in a way that is different from what one does in BEM. We solve the Maxwell's equations on the inside of each scattering object, as an initial boundary value problem, and use the integral identities to supply the boundary values needed in order to make the initial boundary value problem for Maxwell, well posed. This kind of approach for solving electromagnetic scattering problems was first proposed in 1972 by E. Wolf and D. N. Pattanyak [15] in the context of stationary linear scattering and was based on the Ewald-Oseen optical extinction theorem. For this reason we call this particular way of reformulating the electromagnetic scattering problem for the Ewald Oseen Scattering(EOS) formulation. This reformulation can be applied to any kind of wave scattering situation. It has previously been applied to two toy models of 1D linear and nonlinear transient wave scattering by the authors [16] where the EOS formulations work perfectly well with high accuracy and low computational load and without any instabilities, even at very late times.
In section 2 we show some of the details of the derivation of our EOS formulation for Maxwell's equations and in section 3 we discuss some tests we have run on our numerical implementation of the EOS formulation of Maxwell's equations. In this paper we do not describe the numerical implementational details, like how we handle the singular integrals and issues of numerical stability.
These, mainly very technical considerations, would cloud the main message of the current paper, which is that our EOS formulation of scattering problems works. The technical details pertaining to our choise of implementation, some of which are probably relevant for most numerical implementations of the EOS formulation, will be reported elsewhere in a paper soon to appear. Here we will just note that, just like for the 1D case, the internal numerical scheme, Lax-Wendroff for our case, determines a stability interval for the time step. The difference is that, in the 1D case, the stability interval is purely determined by the internal numerical scheme while in 3D case, there is another lower limit of the stability interval determined by the integral part of the scheme. We also find that the late time instability is highly depended on the features of the scattering materials. Section 4 summarizes what we have achieved and discuss extensions of our work that could be of interest to pursue.

EOS formulations of the 3D Maxwell's equations
In this paper, we investigate an electromagnetic scattering problem described by the 3D Maxwell's equations where J and ρ are the current density and the charge density of free charges. Bound charges and currents determine D and H as functionals of E and B, In the simplest situation, where the response from the bound charges and currents is linear, isotropic, homogeneous and instantaneous, we have where 1 εµ = c is the speed of light in the material. For this particular situation, we have We now rewrite the Maxwell's equations into a form that is a suitable starting point for our EOS formulation of the elctromagnetic wave scattering problem.
First, observe that All fields we consider will be driven by the source that will operate for some finite time interval. This means that at some time in the past t = t 0 , we have and this together with (2.2) and (2.3) imply that for any t ∇ · B(x, t) = 0 holds true. If we now use the equation of charge conservation also holds true at all time. Taking these considerations into account, Maxwell's equations can be written in the following equivalent form In order to complete the model, we must supply an equation of motion for the current The specific form for the functional F is determined by what kind of material response we are considering. In order for the system to lead to an efficient numerical method it is important that the sources ρ, J are confined to some small region. In this paper, in order to be specific, we look at the case of a small metallic object interacting with light. We are not seeking to make a detailed computational investigation of this system, but is rather focused on testing our computational approach with respect to implementational complexity and numerical stability. For this reason we choose the following simple nonlinear model for the metal response of such a system where α, β and γ are constants.
Following the usual approach, it is easy to show that the electric field satisfy the following equation (2.6) Each vector component of equation (2.6) is an inhomogeneous wave equation. Let's suppose the scattering object is confined in a compact homogeneous domain denoted by V 1 while the light source is located in an unbounded domain outside the object which is denoted by V 0 . µ, ε are the magnetic permeability and the electric permittivity with their values µ 1 , ε 1 inside and µ 0 , ε 0 outside respectively. c represents the speed of light, with value c 1 inside and c 0 outside the scattering object V 1 . The sources J 0 and ρ 0 are given and J 1 , ρ 1 are the response sources generated by the metallic object interacting with the light field.
We are now ready to start the construction of the EOS formulations of this scattering problem.
Applying the integral relation for the wave equation (A.6) derived in Appendix A on equation (2.6) in domain V 0 and V 1 respectively, we get with j = 0 representing the outside domain V 0 and j = 1 representing the inside domain V 1 . Here x ∈ V j and n is the unit normal to the boundary, S of V 1 , at the point x ∈ S, pointing out of the domain V 1 . The upper sign applies to the case j = 0 and the lower sign for the case j = 1. The same convention applies to all the following expressions in this section. After a series of algebraic manipulations, starting with (2.7), we obtain (2.8) These manipulations are detailed in Appendix B.
Like the electric field, the magnetic field also satisfies a wave equation After a set of algebraic manipulations, similar to the ones we did for the electric field, we obtain (2.10) The identities (2.8) and (2.10), for the electric and magnetic field, are our version of the general integral identities for the electromagnetic field derived by D.S.Jones [7]. In addition to these two identities we get, in a very similar way, two additional integral identities [17], and for x ∈ V j , j = 0, 1. In the above expressions, where x ∈ V j , j = 0, 1. In the end, we have a full set of the integral identities of the inside and the outside fields expressed by (2.8), (2.10), (2.11) and (2.12) which can be written compactly as We will now derive the boundary integral identities of (2.4) by letting x approach the surface from the inside and the outside of the scattering object V 1 , separately. We observe that, in this limit, strong singularities only appear in the last term of the integrals in (2.13) and (2.14). Hence we are faced with a singular term which takes the form of where A(x , T ) is a vector function with x located on the surface S which is an small disk of radius . If we let x approach a surface point ξ, from the inside of V 1 , along a direction where n is the unit normal vector pointing out of V 1 , at the point ξ, and β is a unit vector along the direction x − ξ, tangential to S, at the same point ξ, we have where η = x − ξ + β. Using spherical coordinates, (2.16) turns into lim →0ˆ2 π 0ˆ 0 ρ (ρ cos θ, ρ sin θ, α) So in the end, where I + and I − denote the limit of equation (2.15) by letting x approach S from the inside and the outside of V 1 respectively. After taking these inside and outside limits, we get the following set of equations where E + is the limit of E 1 with x approaching the surface from the inside of the object while E − is the limit of E 0 with x approaching the surface from the outside of the object. These equations, because of the limits taken, contains singular integrals that must be interpreted as Cauchy principal integrals. Adding the first two equations of (2.17) gives us and adding the last two equations of (2.17) gives us Repeating the derivations we just did for the electric field, give us, in a similar way, the following set of equations for the magnetic field where B + is the limit of B 1 with x approaching the surface from the inside of the object while B − is the limit of B 0 with x approaching the surface from the outside of the object. Also in these equations the singular integrals that occur must be interpreted as Cauchy principal value integrals. So far, we have two outer equations for the outer limit fields E − , B − and two inner equations for the inner limit fields E + , B + . We also have the usual electromagnetic boundary conditions at the surface S which separate regions with different susceptibilities and permittivities It might appear that we have more equations than we need here. The very same problem was encountered earlier while deriving the EOS formulation for two 1D toy models [16]. It appears as if we can use the two outer equations to solve for E − and B − and then use the boundary conditions to find E + and B + . But these field values inside the scattering object cannot in general be consistent with the field values derived directly from the two inner equations for E + and B + . For example, if there is a source inside of V 1 and no source outside of V 1 , the first approach would give vanishing electric and magnetic field whereas the second approach certainly would not. On the other hand, for a given source, the Maxwell equations has a unique solution, which by construction also satisfy all the integral identities.
In order to understand what the problem is, and how to fix it, we must just realize that, from an abstract point of view, we have the following formal situation where A and B are singular but where we know that (2.22) has a unique solution.
In this situation, let us assume that αA + B is nonsingular for some choice of α. Any solution of (2.22) is a solution of (αA + B)X = αb + c, (2.23) and since αA + B is nonsingular the unique solution of the singular system (2.22) must in fact be the unique solution of the nonsingular system (2.23). We know that the solution of Maxwell is unique for a given source, so since the integral equations are equivalent to Maxwell, our four integral equations for the two unknown fields on S must have a unique solution. This happens only if they are singular. Thus in our situation, we can simply add (2.18) and (2.19), which gives and add (2.20) and (2.21), which gives Observe that for any vector A, the following identities hold true n × (n × A) = (n · A)n − (n · n)A, where x ∈ S, n is the unit normal vector pointing out of V 1 at the point x, I is the identity matrix and (2.30) our scheme, which in general will involve multiple, arbitrarily shaped, scattering objects, that include linear and nonlinear optical response. We expect however, that the general nature of these issues will reveal themselves already in the simplest possible setting, where we have one scattering object of rectangular shape. The numerical implementation consists of a domain method for the model (2.4) and (2.5), determining the evolution of the fields inside the scattering object, and a scheme for updating the boundary values of the fields using the integral identities (2.27a) and (2.27b). For the internal domain method we choose to use a combination of Lax-Wendroff on (2.4) and modified Euler's method on (2.5), this is similar to what we did for the simple 1D case [16], previously.
For the boundary part of the scheme we use the mid-point rule to the nonsingular integrals appearing in (2.27). The treatment of the singular integrals is technical and rather lengthy and will therefore be reported elsewhere [18].
Here it is enough to note that we calculate the singular integrals by reducing them to a singular core, which we calculate exactly, and nonsingular surface and line integrals, that we calculate numerically. The reductions proceed through a nontrivial use of well known integral identities. For the inside of the scattering object we will use a rectangular grid. This grid is however not uniform close to the boundary. This is because the grid has to support both the discrete approximations to the partial derivatives and discrete approximations to the integrals, used to update the boundary values based on the current and previous internal values of the fields. The fact that in our scheme the boundary values are dynamical variables, enforce some special difference rules that applies close to the boundary. This is an extra complication for our scheme, but they are manageable, and will be part of any scheme that implements the EOS formulation introduced in this paper. Details are given in [18].
What we do in this section is to report on some tests that we have run on our scheme. A usual approach to testing of numerical implementations involve finding exact special solutions corresponding to special source functions. In this section we do not use this approach, but rather use an artificial source test to verify the correctness of our EOS formulations. The basic idea behind the artificial source test, of some numerical scheme designed for a system of PDEs, L ψ = 0 , is to slightly modify the system by adding an arbitrary source to all the equations in the system, creating a new modified system L ψ = g . This modification typically lead to minimal modifications to the numerical scheme, where most of the effort and complexity are usually spent on the derivatives and nonlinear terms. For the equations, however, the presence of the sources change the situation completely. This is because the presence of the added sources implies that any function is a solution to the equations for some choice of sources. Thus we can pick a function ψ 0 and insert it into the model and calculate the source function g 0 = L ψ 0 so that our chosen function is a solution to the extended equation. Finally we run the numerical scheme with the calculated source function and compare the numerical solution with the exact solution ψ 0 .
A modified model of (2.4) and (2.5) with artificial sources is generally con-structed by where ϕ 1 , ϕ 2 , ϕ 4 and ϕ 5 are a set of vector functions and ϕ 3 is a scale function.
Based on this we suppose ϕ 1 = 0 and ϕ 4 = 0 which can simplify the choice of the exact solutions. We also observe that if ϕ 2 and ϕ 3 are set to be both 0, then the continuity equation ∂ t ρ + ∇ · J = 0 is automatically satisfied. So in the end, the source extended model is given by For model (3.1), any choice ofẼ,B can be a solution if the artificial source is given by whereJ andρ are given respectively bỹ Due to (3.1a), we can simply choose a vector function φ, such that For a general case where the electromagnetic fields inside the object are produced by the outside source, we set up the outside source J 0 and ρ 0 to be a combination of a bump function in time and a delta function in space which is easily integrated in space. In order to satisfy the continuity equation we can choose a vector function ϕ such that   with x 0 = (x 0 , y 0 , z 0 ). Our numerical scheme, being explicit, is not unconditionally stable. There is however a stable range, τ 1 < τ < τ 2 , for the time step ∆t ∆t = τ c 1 M in{∆x, ∆y, ∆z}, where τ 1 and τ 2 determine the lower and upper boundaries of the stability range for the sceheme. We have carefully investigated the source of the upper and lower bound of the range and how the width of the stability range depends on material parameters. It is not appropriate to include these fairly techincal numerical investigations here, a full discussion will be presented elsewhere in [18]. Here it is enough to note that the source of the lower stability bound is the numerical implementation of the boundary part of the algoritm, and the source of the upper bound is the numerical implementation of the domain part of the algorithm, for our case this is a combination of Lax-Wendroff for the electromagnetic fields, and modified Euler for the current .

Conclusions
In this paper we have showed that our EOS formulation of electromagnetic scattering can be accurately and stably implemented using one particular choice of numerical scheme for the inside of the objects and for the integral representation of the boundary values required by the inside scheme. For a stable numerical solution, the time step needs to be confined in some range, where we have found that this range is not only determined by the internal domain-base method due to the non-uniform grids but also determined by the boundary integral representations. Discussions on how the internal non-uniform grids and the boundary integral representation effect the time stable range is reported in [16] which will be published elsewhere. It is worth stressing that the existence of the stability range and its width depends not only on the material parameters but certainly on choices made for the numerical implementation of the boundary part and domain part of the algorithm. In principle, any numerical scheme can be used for the domain part of the algorithm, also the extremely well established FDTD method. It would be interesting to see how this method would perform with respect to stability. There is also the question of going fully implicit, both for the boundary part and the domain part of the algorithm. One would think that this would have a chance of producing an unconditionally stable algorithm for our EOS formulation for Maxwell's equations.

Appendices
A The integral identity for a 3D wave equation We will start by considering a wave equation in 3D where x = (x, y, z), c is the propagation speed and Let D × T be a given space-time domain. We will assume that the source ρ(x, t) is entirely contained in D × T . The operator is formally self adjoint. Observe that for any pair of functions defined in D × T we have This is the fundamental integral identity for the wave equation in 3D. Next we will need the advanced Green's function for L which is given by where Q(y, s) satisfies 1 c 2 ∂ ss Q(y, s) − ∇ 2 Q(y, s) = δ(y)δ(s), s < 0, (A .3) and Q(y, s) = 0, s > 0.
Because of translational invariance, we can take the Fourier transform of (A.3) and get Applying the inverse Fourier transform on (A.4) gives and C ε is a contour slightly below the real axis. See figure ??. Observe that the integrand has simple poles at z = ±cζ, so if s > 0, we must close the contour in the lower half-plane. By Cauchy, q(ξ, s) = 0.
If s < 0, we must now close the contour in the upper half plane and Cauchy's theorem gives where r = |y|. Similarly, In the end, Since if s < 0, δ(r − sc) = 0, finally Q(y, s) = 1 4πr δ(s + r c ).
Thus our advanced Green's function is We now apply (A.5) to the fundamental integral identity (A.2). Let ϕ(x, t) be a solution to the wave equation (A.1) Lϕ(x, t) = ρ(x, t), Since t 2 > t, for the advance Green's function And since the field is entirely driven by the source ϕ(x, t 1 ) = 0, we get