Solution of the inhomogeneous Maxwell ’ s equations using a Born series

An algorithm for the numerical solution of the inhomogeneous Maxwell’s equations is presented. The algorithm solves the inhomogeneous vector wave equation of the electric field by writing the solution as a convergent Born series. Compared to two dimensional finite difference time domain calculations, solutions showing the same accuracy can be calculated more than three orders of magnitude faster. c © 2017 Optical Society of America OCIS codes: (260.2110) Electromagnetic optics; (290.5825) Scattering theory; (290.4020) Mie theory. References and links 1. W. Drexler and J. G. Fujimoto, eds., Optical Coherence Tomography (Springer, 2008). 2. D. R. Lynch, K. D. Paulsen, and J. W. Strohbehn, “Finite element solution of Maxwell’s equations for hyperthermia treatment planning,” J. Comput. Phys. 58, 246–269 (1985). 3. R. L. Gordon and C. A. Mack, “Lithography simulation employing rigorous solutions to Maxwell’s equations,” Proc. SPIE 3334, 176–196 (1998). 4. D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science 314, 977–980 (2006). 5. A. Rottler, B. Krüger, D. Heitmann, D. Pfannkuche, and S. Mendach, “Route towards cylindrical cloaking at visible frequencies using an optimization algorithm,” Phys. Rev. B 86, 245120 (2012). 6. D. W. Mackowski and M. I. Mishchenko, “Calculation of the T matrix and the scattering matrix for ensembles of spheres,” J. Opt. Soc. Am. A 13, 2266–2278 (1996). 7. J. Schäfer and A. Kienle, “Scattering of light by multiple dielectric cylinders: comparison of radiative transfer and Maxwell theory,” Opt. Lett. 33, 2413–2415 (2008). 8. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley-VCH, 1998). 9. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). 10. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite Difference Time-Domain Method (Artech House, 1995). 11. Q. H. Liu, “The PSTD algorithm: A time-domain method requiring only two cells per wavelength,” Microw. Opt. Technol. Lett. 15, 158–165 (1997). 12. Q. H. Liu and G. Zhao, “Review of PSTD methods for transient electromagnetics,” Int. J. Numer. Modell. Electron. Networks Devices Fields 17, 299–323 (2004). 13. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705-714 (1973). 14. M. Yurkin and A. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” J. Quant. Spectrosc. Radiat. Transf. 106, 558–589 (2007). 15. F. Assous, P. Degond, E. Heintze, P. Raviart, and J. Segre, “On a finite-element method for solving the threedimensional Maxwell equations,” J. Comput. Phys. 109, 222–237 (1993). 16. G. Osnabrugge, S. Leedumrongwatthanakun, and I. M. Vellekoop, “A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media,” J. Comput. Phys. 322, 113–124 (2016). 17. Y. Saad and M. H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. and Stat. Comput. 7, 856–869 (1986). 18. J.-P. Schäfer, “http://www.mathworks.com/matlabcentral/fileexchange/36831-matscat,” Last accessed October 13, 2014, 3:47 pm. 19. J. Schäfer, S.-C. Lee, and A. Kienle, “Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence,” J. Quant. Spectrosc. Radiat. Transf. 113, 2113–2123 (2012). 20. S.-C. Lee, “Scattering by closely-spaced radially-stratified parallel cylinders,” J. Quant. Spectrosc. Radiat. Transf. 48, 119–130 (1992). 21. K. Moriya, “Light scattering from defects in crystals: Scattering by dislocations,” Philos. Mag. B 64, 425–445 (1991). Vol. 25, No. 21 | 16 Oct 2017 | OPTICS EXPRESS 25165 #303477 https://doi.org/10.1364/OE.25.025165 Journal © 2017 Received 2 Aug 2017; revised 11 Sep 2017; accepted 13 Sep 2017; published 3 Oct 2017 22. J. Stark, T. Rothe, S. Kieß, S. Simon, and A. Kienle, “Light scattering microscopy measurements of single nuclei compared with GPU-accelerated FDTD simulations,” Phys. Med. Biol. 61, 2749-2761 (2016).


Introduction
Today there are numerous applications where solutions of Maxwell's equations in inhomogeneous media are required. These applications include the propagation of electromagnetic waves in biological materials [1,2], in lithography [3], and in optical metamaterials [4,5].
While there are analytic solutions for simple geometries such as spheres and infinite cylinders, most real scattering problems need numerical calculations for a realistic description [6][7][8]. For this purpose several numerical methods have been developed. These methods include the finite difference time domain (FDTD) method [9,10] where the time integration of Maxwell's equations is performed on a discrete grid using finite differences to calculate the spatial derivatives. For this method, very small cell sizes have to be used. This restriction can partly be overcome by using the pseudo spectral time domain (PSTD) method [11,12] where the spatial derivatives are calculated using a Fourier transformation. A further method is the discrete dipole approximation (DDA), [13,14] where the system is discretized by small polarizable grains. Another method is the solution using finite element methods (FEM) [2,15].
Recently, Osnabrugge et al. [16] found a solution of the Helmholtz equation of the electric field that uses plane waves as trial functions. The resulting linear system can then be solved using a Born series. The advantage of this method is that only one electric-field distribution needs to be stored at a time. In contrast, Krylov subspace solvers like the generalized minimization of residuals [17] need to store several basis vectors of the Krylov subspace. The Helmholtz equation is an approximation of the wave equation that emerges from Maxwell's equations. The divergence of the electric field vanishes only in two dimensions for waves with electric field perpendicular to the scattering plane. Thus an accurate solution is possible in two dimensions with one distinct polarization only. In all other cases the Helmholtz equation is an approximation to Maxwell's equations. This approximation looses accuracy when the divergence of the electric field increases.
Here, we show how the method of Osnabrugge et al. can be extended to solve the full wave equation, i. e. Maxwell's equations. This extension allows for a better accuracy of the calculated field in systems with a strong spatial change of the refractive index for all polarizations.
This work is organized as follows. In the remaining part of this section the idea of Osnabrugge et al. is shortly sketched and the Helmholtz equation as an approximation of the wave equation that can be derived from Maxwell's equations is discussed. In the second part we describe the new algorithm with special attention to its comparison to the work of Osnabrugge et al. In the third and fourth part we present numerical results using the described algorithm. These results are compared with an analytical theory for scattering by multiple cylinders [7,[18][19][20]. Finally, a conclusion is given.

Maxwell's equations and the Helmholtz equation
We start with Maxwell's equations where the electric and magnetic fields are harmonic oscillations in time with a given frequency ω. The values of and μ are constant in time and no current or free charges are present. The fields then become E(t) = Ee −iωt and H (t) = He −iωt .
We get Taking the curl of Eq.
(3) and inserting Eq. (4) one gets assuming that μ is not spatially dependent, which is usually assumed, e. g. for biophotonics. This equation can be rewritten as Helmholtz equation (6) with the spatially dependent wave number k with k 2 = ω 2 μ . If the first term on the lefthand side can be neglected, we arrive at the Helmholtz equation. For the Helmholtz equation in Cartesian coordinates the components of the electric field do not mix. Therefore, each component has to fulfill the scalar Helmholtz equation where ψ is one of the three components of the electric field and s is a source term. Except the source term Eqs. (7) and (6) are the same if ∇( ∇ · E) = 0. This expression can be rewritten as − ∇ E · ∇ ln( ) = 0 (8) by making use of and Eq. (1). Plugging this expression into Eq. (6) one finds that the Helmholtz equation is identical to the wave equation if ∇ ln( ) = 0 or E ⊥ ∇ ln( ). The first condition is identical to the absence of any scatterer in the system. The second condition is for example fulfilled for two-dimensional scattering where the electric field is perpendicular to the scattering plane. Furthermore it is possible to use the Helmholtz equation as an approximation for some other systems where the term ∇( ∇ · E) is small compared to the other two terms in Eq. (6). This is the case if varies slowly on the length scale of a wavelength.

Numerical solution of the Helmholtz equation using a Born series
Equation (7) can be solved numerically using a Born series as shown by Osnabrugge et al [16]. They write Eq. (7) as where k 0 is a constant background wave number and δ is a constant real number. The left-hand side is independent of the position. It is therefore possible to derive a Green's function g( r) by solving where δ D is the Dirac function. After a Fourier transform to k space this equation becomes since a derivative in real space transforms to a multiplication in k space. The Fourier transform of the Green's function then readŝ For the calculation of the Green's function in real space one needs to distinguish between twodimensional and three-dimensional systems. The solutions can be found in Eqs. (35) and (39), respectively. Writing Eq. (10) on a discrete grid, we get with the matrices G and V representing the discretized Green's function g and potential v, respectively. The vectors Ψ and S contain the fields ψ and s at every grid point, respectively. This can be rewritten as or in the Born series Here, 1 1 is the identity matrix. Cutting the sum after the first term would give a result similar to the perturbative Born approximation [21]. The problem is that even with higher terms included this series converges only if the magnitudes of all eigenvalues of GV are smaller than one.
To overcome this problem, Osnabrugge et al. introduce the preconditioner [16] It can be seen that γ is a diagonal matrix with the entries Multiplying Eq. (14) with γ yields This can be solved for Ψ as and written as the series From this series the iterative solution with Ψ 0 = 0 can be derived.
for which this iteration will always converge if is a unitary matrix. However from Eq. (18) one can see that if the value of k in at least one cell fulfills the conditions Re(k 2 − k 2 0 ) = 0 and the matrix γ becomes singular. In this case the multiplication in Eq. (19) changes the solution of the linear system. Thus the restriction in Eq. (23) is replaced with the slightly stricter form that ensures that the preconditioner is always non-singular.

Algorithm
In this section our algorithm for the solution of the wave equation is described. For its derivation we follow the ansatz used by Osnabrugge et al. in their calculations for the Helmholtz equation [16]. For the wave equation we now have to solve Here s is a source term that depends on the incident light. This term will be discussed later and will be calculated in Eq. (61). From Eq. (28) the Green's function is calculated by solving where the Green's function g( r) is a 3 × 3 matrix. As shown above for the Helmholtz equation it is convenient to transform Eq. (29) into Fourier space. This yields employing the fact that a differentiation in real space becomes a multiplication in k space. In the following one needs to distinguish between the two-dimensional and three-dimensional cases.

Two-dimensional system
In the two-dimensional system there is no change in z direction. Therefore the derivatives in z direction vanish and Eq. (30) becomes It becomes readily visible that the system separates into one equation for polarization in the xy plane and another for the polarization perpendicular to it. The Green's function is calculated from solving Eq. (31) forĝ( k). For the matrix elements one findŝ where k 2 = k 2 x + k 2 y . For the transformation back to real space we use Here K n (ar) is the modified Bessel function of the second kind of order n with the parameter ar. The function in Fourier space depends on the square of a. Thus, a can assume two different values. Without loss of generality we assume that the real part of a is positive. Setting a 2 = −k 2 0 − iδ, the first part of Eq. (32) can be transformed. For the second part we use the expression from which we find the Green's function in real space. The first term is the Green's function for the Helmholtz equation.

Three-dimensional system
For the three-dimensional system the Green's function in k space is derived from the full Eq. (30) including k z . The matrix elements arê To calculate the Green's function in real space we use the identity where a is assumed to have a positive real part as for the two-dimensional case. Setting a 2 = −k 2 0 − iδ allows us to calculate the inverse Fourier transformation of the first part of Eq. (36). To calculate the inverse Fourier transform of the second part we employ This yields where a is the square root of −k 2 0 − iδ with positive real part. As for the two-dimensional case the first term is the Green's function for the Helmholtz equation.

Born series
For the calculation of our generalized Born series we follow the calculations of Osnabrugge et al. for the Helmholtz equation. With the aid of the Green's function for the wave equation one can write In the spatially discretized form this equation reads with a vector Ψ of size 3N containing the electric field at each of the N grid points. G is a dense 3N × 3N matrix. This matrix consists of N 2 3 × 3 blocks that are g( r i − r j ), where r i and r j are the positions of the according grid points. V is a diagonal 3N × 3N matrix. Similar to the calculations of Osnabrugge et al. we rewrite this expression as with the preconditioner and a value of It is worth noting that in contrast to the work of Osnabrugge et al. γ is now a 3N × 3N matrix. Equation (42) can be rewritten as or in the series expansion From this expression the iterative solution with Ψ 0 = 0 can be derived.
For the numerical evaluation one needs to keep in mind that the dense matrix G represents a convolution as can be seen from Eqs. (40) and (41). Since a convolution in real space becomes a local multiplication in Fourier space, one gets with a block diagonal matrixĜ This matrix has N 3 × 3 blocks containing the matrixĝ( k) for different values of k. The multiplication with G has a time complexity of O(N 2 ). Since the multiplication withĜ and the Fourier transforms have time complexities of O(N ) and O(N log(N )), respectively, it is faster to use Eq. (48) than Eq. (47) for computations.

Convergence
The iteration is not necessarily convergent. However, Osnabrugge et al. have proven that this series will always converge if Im(k 2 − k 2 0 ) ≥ 0 is fulfilled everywhere in the sample, meaning that there should be no gain in the medium, that the stricter condition Im(k 2 − k 2 0 ) > 0, i. e. absorption, is fulfilled somewhere in the sample, and that the matrix With this matrix we get It is worth noting that for the two-dimensional case k is always in the xy plane and thus sin(θ) = 1. We now calculateû This consists of two rotations and one unitary matrix. Therefore,û is unitary as well.

Truncating condition
In Eq. (46) the resulting vector Ψ is written as an infinite sum. Therefore, to obtain the exact solution one needs to apply the iteration in Eq. (48) an infinite number of times. Therefore, we need a condition that helps us to decide whether the iteration can be truncated at the current state or not. During the update of Eq. (48) one needs to calculate the value that is the residual of Eq. (41) for Ψ n . Thus the calculation of the new value Ψ n+1 gives us the current residual r n for free. This value can now be used to judge whether the iteration has been converged or further iterations are needed.
We have defined two conditions. The absolute error is given by the maximum of the residual over all cells. For the relative error we divide the residual by the field in the cell before finding the maximum. In the data shown in this paper we have used a maximum error of 0.001 and a relative error of 0.01. Stricter values do not affect the accuracies derived in Sec. 4.

Sources
In Eq. (28) we have defined a source term that includes the incident light. To calculate this term we write the electric field as the sum of an incident field E i that comes from external sources and a field E s that stems from the scattering in the system. Then Eq. (6) reads The incident field alone needs to fulfill the wave equation for the medium that surrounds the simulated system. Here k m is the constant wave number in the medium. In these two equations we have employed the fact that the divergence of E i vanishes as can be seen from Eq. (8), keeping in mind that is constant in the medium. By plugging Eq. (58) into Eq. (57) one gets In a shape that is more in line with Eq. (28) this reads By comparing this with Eq. (28) one easily obtains

Boundary conditions
As boundary conditions we use the ones introduced by Osnabrugge et al. The simulated system is surrounded by an absorbing region of width l, in which the wave number is given by with the polynomial and the distance d from the simulated system. The scatterer is a cylinder whose center is not necessarily aligned with the rectangular grid. Its circumference is denoted by the blue line. The grid cells (solid black lines) are divided into subcells (dashed black lines). For each subcell a black dot at its center denotes that the subcell is considered to be inside the cylinder. The effective refractive index of the grid cell (denoted by the gray scale) depends on the number of subcells that are considered to be inside the cylinder. (right) Our test system consists of five cylinders with a radius of 5 μm and a refractive index n 2 that are placed in a homogeneous medium with a refractive index of n 1 . The near field is calculated in a square that contains all cylinders. From the calculated near field the far field is calculated by a near to far field transformation on the boundary of this square.

Anti-aliasing
Before we start the calculation using our algorithm we need to assign a refractive index to each simulation cell of the grid. One possible method is to assign the refractive index at the center of each cell. To get a smoother surface and to reduce stair casing effects in calculations with round surfaces we decided to calculate the refractive index using an anti-aliasing routine. For this calculation we divide each cell in n s × n s subcells as depicted in Fig. 1 for n s = 2. We then assign to each subcell the refractive index at the center of the respective subcell. Subsequently, the average of the refractive indices in all subcells of one cell is calculated and assigned to the cell. For the results shown in this work we use n s = 10.

Runtime complexity
In this section we calculate the dependence of the simulation time on the size of the system. For simplicity we assume a system with the same size L and number of cells N in each spatial direction, i. e. the grid contains N 2 and N 3 cells for two and three dimensions, respectively.

FDTD
In the case of FDTD the time that is needed for the calculation of one simulation step scales as with the dimension D of the system. For the simulation of a time t with a time step of Δt this yields  To ensure the stability of the simulation the time step has to be given by where the Courant number S is a constant that depends on the dimension, c is the speed of light in the simulated medium, and Δ is the grid spacing. To describe the interaction with scatterers accurately we need to ensure that the light is able to travel through the grid n times. The maximum distance in the grid is given by its diagonal. The maximum time we need to Vol. 25,No. 21 | 16 Oct 2017 | OPTICS EXPRESS 25175 simulate is therefore This yields if n is assumed to be independent of the system size.
For the determination of a proper value for the grid spacing we note that for small systems the grid spacing needs to scale linearly with the wavelength to keep the error from the finite differences constant. In addition to this error we have to take into account that in FDTD simulations there is a numerical dispersion relation that is different to the analytical one. The error due to the occurring phase lag becomes maximal at the boundaries of the simulated area. This maximum dispersion determines the accuracy for large systems. The numerical wave number k n can be calculated as [10] with the analytical wave number k a = 2π/λ. We note that the value of k n depends on the propagation direction of the wave with respect to the grid axes. For simplicity we use the dispersion relation for a wave propagating parallel to one grid axis only. We now aim to keep the maximum phase lag φ max constant. For this we approximate φ max by the phase lag that occurs while the wave moves through the sample once. The maximum phase lag can thus be written as where λ is the wavelength. Solving this relation for Δ yields The value on the abscissa is the accuracy that is defined in Eq. (76) and is determined from a comparison with the respective analytical results. The results from the Born series and the FDTD are fitted with a power law for both sets of refractive indices. A comparison of these fits shows that for an accuracy of 0.1 the Born series is up to a factor of 7300 faster, depending on the simulated system.
leading to the conclusion that Δ has to scale with the inverse square root of L. Indeed, numerical experiments with a wave traveling along the x axis of an empty grid with S = 0.5 and a grid spacing of Δ = 0.2λ 1.5 /L 0.5 have shown that the phase lag over the whole grid is 18 degrees, irrespective of the grid size. We thus get a calculation time that scales with

This algorithm
In the calculation of one iteration of the present algorithm, the part with the highest runtime complexity per iteration is the calculation of the Fourier transformation. The runtime for a D dimensional Fourier transformation is of the order Since this algorithm contains no time integration there is no dispersion relation and N is proportional to L/λ. Osnabrugge et al [16]. show that for each additional iteration of the Born series the light travels a distance of about 2k 0 /δ. This expression scales linearly with the wavelength. The number of iterations then scales with L/λ. We find Vol. 25, No. 21 | 16 Oct 2017 | OPTICS EXPRESS 25179

Numerical tests
In this section we present a test of the accuracy and the speed of the present algorithm and compare the results to analytical calculations using an analytical solution [7,[18][19][20] and FDTD simulations.

Calculated system
For the test we use a square region with an edge length of 100 μm that is filled with five cylinders with a radius of 5 μm. The background has a refractive index of n 1 while the cylinders have a refractive index denoted by n 2 . The positions of the cylinders are given in Tab. 1. They are positioned so that there is a small space between them and the boundary. The full system is depicted in Fig. 1. Table 1. Positions of the five cylinders in the calculated system relative to the center of the simulation area. The system is depicted in Fig. 1. To get a deeper insight into the effect of the material parameters we performed calculations for two different sets of refractive indices. In the first set with n 1 = 1 and n 2 = 1.2 the cylinders are situated in vacuum. The second set with n 1 = 1.36 and n 2 = 1.39 has parameters that are typical for biological tissue [22]. In both cases the incident light is traveling in x direction and has a vacuum wavelength of 1 μm. The incident field is E i = E 0 (0, 1, 1) T .
For the boundary conditions in the Born series calculations we used the parameters α = 0.35/μm, l = 25 μm, and N = 4 as described in Sec. 2.7. A reduction of l, i. e. the size of the absorbing region, would make an increase of the damping α necessary to ensure that a wave is sufficiently damped and that no reflections at the edges of the grid occur. But a larger value for α leads to reflections in the transition region where the absorption sets in.
The area that is covered by the absorbing boundaries is (L + 2l) 2 − L 2 = 12500 μm 2 . This area is by a factor of 1.25 larger than the area covered by the region containing the scatterers. This means that for the following two-dimensional benchmarks the size of the whole simulation grid is about two times larger than that used for the finite difference time domain simulations. To determine how the simulation time is affected by this increase in size we start from the expression in Eq. (74) and approximate the calculation time T with This coarse approximation leads to an increase in simulation time by a factor of 3.7. Thus more efficient boundary conditions could significantly reduce the simulation time of the Born series method in the future.

Near field
The near field calculated using the present algorithm and the difference to the analytical calculation are shown in Figs. 2 and 3.
It can be seen that the cylinders act like converging lenses. Directly behind the cylinders there is a region where the fields are stronger than the fields of the incident wave. Staying behind the cylinder but increasing the distance beyond the focal length the field strongly decreases. This leads to a strong interaction between cylinder one (lower left corner) and cylinder two (lower right corner) governed by the fact that there is nearly no light directly impinging on the cylinder two due to shading from cylinder one.
In the data from the Born series as shown in Fig. 2, one can see that for the cylinders in vacuum the positions with maximum error are the surfaces of the cylinders, especially at positions where the electric field becomes discontinuous. These positions are for the x and y components at the left/right and at the top/bottom surface of the cylinders, respectively. The z component of the electric field is always continuous, leading to a much smaller maximum error. For the system with refractive indices of biological tissue the error at the surfaces of the cylinders is much smaller since the jump in the refractive index is much smaller. The error is about an order of magnitude smaller, even though the wavelength in the biological tissue is smaller and thus there are less cells per wavelength.
In Figs. 2(a) and 2(g) one finds that there is a sizable x component of the electric field even though the incident field has a vanishing x component. In the calculations using the Helmholtz equation the x component would vanish since in Cartesian coordinates the Helmholtz equation neglects the interaction between different components of the field as discussed in Sec. 1.1.
In the FDTD data in Fig. 3 there is also a region with large errors around the cylinders. In addition, there is an error at the left and right boundaries due to the numerical dispersion of the FDTD.

Far field
With the present algorithm and FDTD we calculate the far field at a distance of 1 m using a near to far field transformation. In the analytical calculations the far field is directly calculated using a series expansion for large distances. The results for both sets of refractive indices are shown in Figs. 4 and 5. In these figures the results from the present model and the results of FDTD are compared with the analytical results. Assuming that the analytical results are exact, it can be seen that the accuracy of the results from the FDTD calculations is much worse than the accuracy of the results from the present algorithm, even though the FDTD calculations have a much higher resolution. Figures 4 and 5 show results with different accuracies in the far field. For the comparison of the calculation time of the far fields we aim to compare results with the same quality. For this task we define the accuracy of the far field as where E and E ⊥ are the electric fields parallel and perpendicular to the scattering plane, respectively. The φ j are the angles at which the far field is calculated. Here, we use the median instead of the mean to be less sensitive to the large relative errors that occur at positions with very small field amplitude, i. e. the spikes that can be seen in Figs. 4 and 5.
The calculation time is measured on a system with two Intel(R) Xeon(R) E5-2640 v4 CPUs at 2.40 GHz with 10 cores each. The results of this measurement are depicted in Fig. 6. The measured values are fitted with a power law. It can be seen that for an accuracy of 0.1 the present algorithm is about a factor of 5000 faster than the FDTD simulations, depending on the simulated system.

Assets and drawbacks of the algorithm
As it can be seen from Fig. 6, compared to the FDTD the present algorithm performs better for the case of biomaterials than for the cylinders in vacuum. The reason for this is that the difference of the refractive indices in the cylinders and the medium is smaller in the biomaterials. This leads to the fact that smaller values of δ are possible (see Eq. (44)). From Eqs. (43) and (47) it can be seen that a lower value of δ directly leads to a larger value of γ and thus to an increasing change of the field in each iteration. In general the algorithm is best suited for systems with a small variation in the refractive index.

Conclusion
In this work we extend the solution of the Helmholtz equation of Osnabrugge et al. to the full vectorial wave equation derived from Maxwell's equations. The present algorithm calculates the solution of the wave equation in the frequency domain. It allows us to calculate the solution for arbitrary polarizations of the wave, scattering geometry, and refractive indices of the material while it shows its best performance for systems in which the difference between the refractive indices is small. Comparison of our algorithm to FDTD calculations shows better accuracy even for much larger cell sizes. For far fields of comparable quality we find that our algorithm is by almost four orders of magnitude faster.