On the Ewald Oseen scattering formulation for light scattering: stability, singularity and parallelization

In this paper we discuss some of the mathematical and numerical issues that have to be addressed when calculating wave scattering using the Ewald Oseen scattering (EOS) formulation which is newly developed for solving linear and nonlinear scattering problems. The discussion is framed in context of light scattering by objects whose optical response can be of a nonlinear and/or inhomogeneous nature. The discussions address two issues that, more likely than not, will be part of any investigation of wave scattering using the EOS approach.


Introduction
A new hybrid numerical approach for solving linear and nonlinear scattering problems, the Ewald Oseen scattering (EOS) formulation, has recently been introduced and applied to the cases of 1D transient wave scattering [1] and 3D light scattering [2]. The approach combines a domain-based method and a boundary integral representation in such a way that the wave fields inside the scattering objects are updated in time using the domain-based method, while the integral representation is used to update the boundary values of the fields, which are required by the inside domain-based method. In such a way, for the numerical implementations, no numerical grids outside the scattering objects are needed. This greatly reduces the computational complexity and cost compared to fully domain based methods like the finite difference time domain (FDTD) method or the finite element methods. The method can handle inhomogeneous and/or nonlinear optical response, and include the time dependent boundary element method (TBEM), as a special case.
For the case of 1D transient wave scattering [1], we illustrate our approach by implementing the EOS formulation for two different toy models for scattering electromagnetic waves which are our particular interest. The method solves the model equations accurately and efficiently. This does not mean that only models that in some way are related to electromagnetic scattering can be subject to our EOS approach. The only one requirement for the EOS approach to be applicable is that it must be possible to derive an explicit integral formulation for the PDEs of interest. However, we do not expect the 1D case to be fully representative for the problems and issues that need to be resolved, while using the EOS formulation to calculate wave scattering. We do, however, expect the case of 3D light scattering [2] to be fairly representative with respect to which problems arise, and also the computational and mathematical severity of these problems. We have seen three types of mathematical and computational issues arise for the case of light scattering which we believe are to be found in any nontrivial application of the EOS formulation to wave scattering.
Firstly, we have the issue of numerical stability. Instabilities in numerical implementations of the EOS formulation can arise from discretization of the domain part of the algorithm but also from discretization of the boundary update part of the algorithm. The numerical instability arising from the boundary part of the algorithm has been noted earlier in the context of transient light scattering from objects that has a linear homogeneous optical response. For this situation, realized for example in antenna theory, the boundary part of the EOS algorithm can be disconnected from the domain part of the algorithm, which in this case can be discarded. The EOS formulation becomes a pure boundary update algorithm which is solving a set integro-differential equations located on the boundary of the scattering objects. These integro-differential equations, which are the defining equations for TBEM, are subject to an instability that, in many common situations, strikes at late times. This late time instability is a major nuisance, and has prevented TBEM from being more widely applied than it is today. The sources of these instabilities are not yet fully understood, but we believe that our investigation of light scattering using the EOS approach, gives some new insight into the origin of these instabilities.
Even without a true understanding of the underlying causes of the late time instability, efforts have been made and several techniques have, over the last several decades, been developed with the goal of improving the stabilities of the numerical schemes designed to solve the integro-differential equations underlying TBEM.
Broadly speaking, there are two different directions that have been pursued. One direction is to delay or remove the late time instability by applying increasing accurate spatial integration schemes [3][4][5][6][7][8]. For instance Weile and his co-authors have published a series of articles focused on illustrating the dependence of the stability on the different numerical integration schemes [3][4][5]. The other direction is aimed at designing more stable time discretization schemes. Bluck and his co-authors developed a stable, but implicit numerical method, [7,8] for the integro-differential equations underlying TBEM, for the case when the magnetic response is the dominating one. These are the so called magnetic field integral equations. Some authors have reported some success in mitigating the instability by both making better approximations to the integrals and also applying improved algorithms for the time derivatives [9,10].
Our work has not been directly aimed at contributing to this discussion, but, as already noted above, the integro-differential equations discussed by these authors can be seen as a special case of our general EOS approach, and we therefore believe that the insights we have gained on how this long time instability depend on the different pieces of the EOS algorithm, in particular how it depends on the material parameters describing the optical response of the scattering object, do have some relevance to the discussion described above.
Secondly, there is the issue of the singular integrals that appear when the integral part of the EOS algorithm is discretized. This issue is very much present in BEM and in TBEM [11][12][13][14], but they are more prevalent and severe for the EOS formulation, where we have to tackle both surface integrals and volume integrals. We believe that the type of singular integrals, and how to treat them for the case of light scattering, are fairly representative for the level of complexity one will encounter, while applying the EOS approach to wave scattering problems. For this reason we find it appropriate to include a section in this paper, where we discuss relevant types of integrals, and how to treat them.
Thirdly, the fundamental equations underlying both the TBEM and our more general EOS approach to transient wave scattering, are retarded in time. This retardation is unavoidable since their underlying equations can only be derived using spacetime Green's functions. Thus the solutions at a certain time depend on a values of the solutions from a potentially very long previous interval of time. Computationally this means that the method can be very demanding with respect to memory, and it also means that the updating of the boundary values of the fields, which is done by the boundary part of the EOS algorithm, can be very costly. Parallel processing, either using a computational cluster or a shared memory machine can take on these computational tasks. However, whenever large scale parallel processing is needed, the issue of appropriate partitioning of the problem and load balancing inevitably comes into play. In our work the EOS algorithm was implemented on a large cluster, but we will not in this paper report on any of the parallel issues that our EOS approach for light scattering gave rise to. These kind of considerations, which are important in practical terms, but typically have fairly low generality, are somewhat distinct from the mathematical and numerical issues that are the focus of the current paper, and will therefore be reported elsewhere at a later time.
However, the high memory requirement of the EOS approach to light scattering, is something that should be addressed at this point. On the one hand, the EOS approach represents a large, potentially very large, reduction in memory use, as compared to fully domain based methods, since only the surface and inside of the scattering objects has to be discretized. On the other hand, because of the retardation, there is a large, potentially very large increase in memory use compared to the memory usage needed by the domain part of the algorithm. It is appropriate to ask if anything has been gained with respect to memory usage compared to a fully domain based method like the FDTD method? We do not, as of yet, know the answer to this question, and the answer is almost certainly not going to be a simple one. It will probably depend on the detailed structure of the problems like the nature of the source, the number, shape and distribution of scattering objects etc. However, even if the memory usage for purely domain based methods and our EOS approach are roughly the same for many problems of interest, our approach avoid many of the sources of problems that need to be taken into account while using purely domain based methods. These are problems like stair-casing at sharp interfaces defining the scattering objects, issues of accuracy, stability and complexity associated with the use of multiple grids in order to accommodate the possibly different geometric shapes of the scattering objects and the need to minimize the reflection from the boundary of the finite computational box. The EOS approach is not subject to any of these problems.
In this paper our effort are aimed towards testing the EOS formulations of light scattering with respect to implementation complexity and numerical stability. Thus we illustrate the method by the simplest situation where we have single scattering object in the form of a rectangular box.
In section 2 we analyze the numerical stability of our EOS scheme for light scattering by using eigenvalues of the matrix defining the linearized version of the scheme exactly like for the case of 1D wave scattering [1]. We find, just like for the 1D case, that the internal numerical scheme, Lax-Wendroff for our case determines a stability interval for the time step. In the 1D case, the stability interval of the EOS formulation is purely determined by the internal numerical scheme. However for the 3D case, there is another lower limit of the stability interval determined by the integral part of the scheme which leads to the situation where the lower limit of the stability interval is determined by the integral equations, and the upper limit is determined by the internal numerical scheme. We find that the late time instability is highly depended on the features of the scattering materials and specifically, it is directly related to the values of the relative magnetic permeability μ 1 and the relative electric permittivity ε 1 . Using this we prove that, for the relative permeability and permittivity in a certain range, the numerical scheme for our EOS formulation of light scattering, works well and is without any late time instabilities. The late time instability is only observed for high relative electric permittivity or high relative magnetic permeability. We also observe that the lower limit of the stability interval for the time step is more sensitive to relative differences in magnetic permeability μ 1 than electric permittivity ε 1 between the inside and outside of the scattering objects.
In section 3 we present the singular integrals that appear in our EOS formulation for light scattering and the techniques we use to reduce their calculation to a singular core, which we calculate exactly, and a regular part which we calculate numerically.

Stability
In this section we discuss instabilities showing up at late times when we discretize the EOS formulation for light scattering. Whether or not the late time instabilities show up, depends on the values of the material parameters defining the problem. The overall method is far to complex for an analytical investigation of the stability to be feasible, but using numerical calculation of the eigenvalues of a linearization of the system of difference equations defining the numerical implementation of the EOS formulation, supplemented by running of the full algorithm, we find that the domain part and the boundary part of the algorithm contribute to the instability separately and in different ways. The focus of this section is to disentangle these two contributions to the instability. For the domain part of the algorithm we use Lax-Wendroff, which is an explicit method. The discrete grid inside the scattering object must, for the EOS formulation of light scattering, support both discrete versions of the partial derivatives, and also discretizations of the integrals defining the boundary update part of the algorithm. For this reason the grid is nonuniform close to the boundary. The discretization of the domain part of the algorithm takes the form of a vector iteration where Q is a vector containing the components of the electric field and the magnetic field at all points of the grid with a size´´Ń N N 6 x y z , where N x , N y and N z are the numbers of grid points in the x, y and z directions. The entries of the matrix M are presented in appendix. In order to get a stable numerical solution, as discussed in [1], the largest eigenvalue of the matrix M must have a norm smaller than 1. For the non-uniform grids and the discretizations in [2], we find that the vector where τ=c 1 Δt/Δx. Figure 1 illustrates the intensity of the electric field at a specific point inside the object, as a function of time, for different values of τ. The instability, which in the TBEM literature is called the late time instability, is illustrated in the second panel of figure 1. As we mentioned in the introduction in the paper, the term late time instability has been much used in the community that is focused on time dependent boundary element method. We believe that in their domain of application, like antenna theory, the physical parameters are such that the largest eigenvalue for the iteration is always only slightly bigger than 1, like it is in panel two of figure 1. That is why the instability always shows up at late times. In panel three of the figures we are deeper into the unstable domain for τ, and the larges eigenvalue is now so large that it destroys the whole calculation. The late time instability has thus been transformed into an early time instability. Note that the outside source in figure 1 is the same as in [2].
In our numerical experiments, we found that the stable range of the EOS formulations is not only restricted by the eigenvalues of the matrix M, but is also restricted by the boundary integral identities through the relative electric permittivity ε 1 and the relative magnetic permeability μ 1 . Figure 2 shows how the stability depends on the values of ε 1 , and figure 3 shows how it depends on the values of μ 1 . Together, they tell us that increasing the electric permittivity or the magnetic permeability narrows the stable range. Figure 3 also tells us that μ 1 and ε 1 do not affect the stability of the full scheme in the same way. It seems that the method is more sensitive to μ 1 than ε 1 . After a series of numerical experiments, our conclusion is that, for an explicit numerical method like the one we are using, the lower limit of the stable range of the EOS formulation is restricted by the electric permittivity ε 1 and the magnetic permeability μ 1 while the upper limit of the stable range is determined by the inside domain-based method. This conjecture is verified by the following two tests.

Instabilities coming from the domain-based method
For the first test we consider a homogeneous model without current and charge inside the object which implies μ 1 =μ 0 , ε 1 =ε 0 , J 1 =0 and ρ 1 =0. Under these assumptions, the electric field and the magnetic field are continuous across the surfaces where + + E B , and --E B , are the integral representations of the solutions on the surface by taking the limit from the inside and the outside of the object respectively. The electric field inside the object can be calculated by the outside sources 1 Also from [2] we have the boundary integral identity On the other hand, [2] gives the integral representations for the inside domain by  Figure 4 compare the solutions calculated in these three ways, where μ 1 , ν 1 and τ have been fixed in the stable range. Both Method 2 and Method 3 are stable and give solutions that agree with the exact solution to high accuracy. In figure 5, τ has been set to be 0.49, and is thus larger than the upper limit of the stable range. The figure shows that Method 2 is now unstable but Method 3 is still stable and equal to the exact solution to high accuracy. The outside sources in

Instabilities coming from the boundary integral identities
In order to investigate the dependence of the stability on μ 1 and ε 1 , we set up a test based on the use of artificial sources as in [2]. The idea is to chose functional forms for an electromagnetic field, and then calculate the sources, charge density and current density, needed for making the chosen field solutions to Maxwell's equations driven by the calculated sources We now calculate the electromagnetic field inside the scattering object in two different ways. Method 1 uses the EOS formulation which combines the Lax-Wendroff method inside and the integral representations of the boundary fields. Method 2 is to calculate the inside field values by only using the Lax-Wendroff method supplemented by the exact boundary values of the electromagnetic field which are the ones we chose while setting up the artificial sources. Figure 6 is the numerical result where the upper limit of the stable range is kept while the values of μ 1 and ε 1 have been chosen to break the lower limit of the stable range of the EOS formulations. It shows that even though the lower limit of the stable range has been broken, Method 2, which only involves the Lax-Wendroff method works perfectly. 5 and 6 tell us that the changing of the lower limit does not effect the stability of the Lax-Wendroff method and the changing of upper limit does not effect the stability of the surface integrals. For a general application where the source is located outside the object and there are current density and electric density inside the scattering object, the EOS formulations does have a range for a stable numerical implementation. The upper limit of the range is determined by the Lax-Wendroff method due to the non-uniform grids and the lower limit is determined by the changing μ 1 and ε 1 . The setting up of the artificial sources and the values of the parameters in figure 6 are the same as the artificial sources in [2]. From figures 5 and 6, we can also see that before the instabilities show up, both the EOS formulations and the Lax-Wendroff method solve the equations accurately.

Calculations of the singular integrals
In this section we introduce a technique to accurately calculate integrals with singularities which can be applied for both the singular volume integrals and the singular surface integrals occurring in the EOS formulations of the 3D Maxwell's equations. Here we illustrate the technique by calculating one

type of singular volume integral
where the integral domain V i,j,k is adjacent to the surfaces of the scattering object and given by with surfaces S , m m=1,2, L,6. Here, Δx, Δy and Δz are the grid parameters in x, y and z directions respectively.
The point x p is centered on one of the surfaces of the scattering object. The geometry is illustrated in figure 7, where n m is the unit normal vector on surface S m pointing out of V i,j,k . The components of the integration variable in (3.1) are given by and let us introduce the quantity r . We want to apply the divergence theorem on (3.1), and therefore need to find a function j(r) that satisfies or equivalently ( ) ( ) j j + ¢ = r r r r 3 1 .
Solving the above equation, we get Because of the singularity on S 1 , we can not apply the divergence theorem directly, however we can write f 1 as (3.2) is not singular any more and can be calculated by 2D Gaussian quadrature. However we will compute f 1 by reducing the surface integral into a line integral, which is also the approach we use to calculate the singular surface integrals appearing in the implementation discussed in this paper.
We first consider the integral over S 2 . The geometry is shown in figure 8. As shown in figure 8, the surface S 2 is bounded by the union of four straight lines = L n ,  On this surface we have and the unit normal is The goal is to use the divergence theorem on this surface integral and thereby reduce it to line integrals over the four lines that forms the boundary of S 2 . We therefore seek a function (¯) j r that satisfies · (¯(¯))j  = + D r r y r 1 , Using the divergence theorem and taking into account of the singularity at¯( where L ò is a semicircle with radius ò centered at pointx and L Ω is the rest of L 21 . Heren n is the unit normal of L 2n , pointing out of S , 2 and¯ n is the unit normal of L ò , pointing out of S 2 .
For the integral over L Ω , we havē For the integral over L ò , using the polar coordinates, we havē  Due to the symmetry of r in V i,j,k along y direction, we have = s s .

2
The calculation of s 3 is similar to the one of s 2 with the final Also due to the symmetry of r in V i,j,k along z direction, we have = s s . 6 3 The only surface integral remaining to be calculated is the one over S . 4 On this surface we have we seek a function (¯) j r that satisfies · (¯(¯))j  = + D r r x r 1 .

2
This equation can be written in the form (¯)¯(¯)j j + ¢ = + D r r r r x 2 1 .

2
Solving the above equation gives  All the line integrals l 21 etc are non-singular and can be calculated accurately using numerical integration.

Summary
In this paper we have, by considering 3D light scattering, discussed some important issues that we believe will be generic for numerical implementations of the EOS formulation for wave scattering. We have shown that the numerical instabilities can be thought as arising separately from the domain part and the boundary update part of the algorithm. We have argued that the instability arising from the boundary part of the algorithm is strongly related to the late time instability noted earlier while solving antenna problems using TBEM. We find that our version of the late time instability can be completely removed by suitably chosen material values, in particular the jump in material values at the boundary of the scattering object should not be too severe.
In the limit where the material parameters simulate the properties of highly conductive metallic surfaces, we observe that our version of the late time instability is always present. Thus the instability interval vanishes in this limit. We take this as an indicator that for situations like in antenna theory, the late time instability should always be present, which it is. We are now aware of work where it has been noted that the instability can be removed by manipulating the material parameters defining the scattering objects. The EOS formulation gives thus different window into the late time instability that might be useful.
We have in our discretization used explicit methods. It would not be easy, but we believe that it is possible to do a fully implicit method for the EOS formulation, such an approach might remove all instabilities, which is the ultimate goal both for TBEM and for our EOS formulation.
In this paper we have also discussed how to calculate singular volume and surface integrals for light scattering. The reason for including this discussion is that we think the type of singular integrals we discuss are generic for the singular integrals that will arise while calculating wave scattering using the EOS approach.

Appendix. Matrix elements
In this section we detail the entries of the updating matrix M in (2.1) where Q is a vector containing the components of the electric field and the magnetic field at all points of the grid with a size´´Ń N N 6 x y z , where N x , N y and N z are the number of grid points in the x, y and z directions. To simplify the writing, we denote L =´Ń N N ,