1. Introduction

Determining how light propagates in heterogeneous media is a notoriously hard problem [1]. Unless the system of interest has symmetries such as periodicity [2], one needs to solve the Maxwell equations ab initio with appropriate boundary conditions. While this may be feasible for relatively simple systems such as Mie scattering [3], multiple scattering of light can lead to many subtle effects [4–6]. The fractal propagation method can simulate light in biological tissues [7]; however, calculating the exact light field distribution in arbitrary large heterogeneous materials remains out of reach for the current generation of computational methods [1]. A further complication is that the media are generally not isotropic, meaning that the refractive index is different depending on the orientation of the field polarization. Such birefringence is common in many samples of interest such as TiO₂, the lipid bilayer cell-membrane [8–10], or muscle fibers [11]. Indeed, birefringence is a valuable, label-free, contrast method [11–16]. Furthermore, biological tissues contain chiral organic molecules such as glucose, a property that is linked to the electro-magnetic coupling and cannot be modeled by anisotropic permittivity. In this paper we investigate the application of the modified Born series [17,18] to solve Maxwell’s equations in large heterogeneous electromagnetic media, characterized by arbitrary linear constitutive relations.

Although finite-element and finite-difference time-domain (FDTD) methods can in principle handle general electromagnetic problems of an arbitrary size, such methods do not scale well to the dimensions relevant in microscopy. Only recently it has become possible to represent the complete field distribution in computer memory for larger samples. Conventional methods require amounts of high-speed access storage that are considerably larger. For P sample points the finite-element method typically requires the representation of a sparse 3P × 3P matrix in memory [19], which is then iteratively applied in order to reach the desired solution. Meanwhile, the FDTD method requires far higher sampling densities to limit the accumulation of errors. A quite different, and more physically meaningful iterative method is the Born series, which calculates the scattered field, with each iteration including the effect of the next highest order of scattering events. This approach is commonly used for the analytical theory of multiple scattering [20,21]. When implemented numerically, all operations can be performed in a constant space (∝ P) and using only fast-Fourier transforms and operations on diagonal matrices. The high computational efficiency and modest memory requirements make the Born series an ideal candidate for solving large-scale electromagnetic problems. However, a severe limitation is that the basic form of the series only converges in the limit of weak scattering. A modified Born series for scalar waves, which converges for any value of the scattering strength, was recently proposed by Osnabrugge et al. [17], and generalized to vector waves by Krüger et al. [18]. Yet, as it stands, this approach is limited to media with an isotropic permittivity and without magnetic properties. The modified Born series method thus excludes a large class of materials such as biological tissues that exhibit birefringence or contain chiral organic compounds, as well as any potential applications to metamaterial structures. Here we generalize the numerical Born series method to arbitrary linear materials, including those with heterogeneous magnetic properties and bi-anisotropy. The modification to the Born series is found to be notably more subtle in these materials. We demonstrate that the algorithm enables large scale calculations of the scattered light field in complex birefringent materials, chiral media, and even materials with negative refractive index. Furthermore, we show that the iteration is robust to numerical errors.

To start, we review the simplest application of the Born series and its modified form, as given in [17,18]. At a fixed frequency ω, the electric and magnetic fields satisfy Maxwell’s equations, ∇ × E(x) = iωB(x), and ∇ × H(x) = j(x) − iωD(x), where j(x) is the electric current density as a function of the spatial coordinate, x. In first instance, we consider the case that the medium is isotropic and non-magnetic, and can thus be characterized in terms of a scalar relative permittivity ∊(x). This quantity connects the displacement field D and magnetic flux density B fields to the electric field E and magnetic field H in the constitutive relations:

The Born series [22] is a physically motivated version of this iterative procedure that was used as a mathematical technique long before the current numerical approaches [20]. It involves splitting the Helmholtz operator, 𝒪, into two parts 𝒪 = 𝒪_i + 𝒪_h, where the inverse of the homogeneous part 𝒪_h is known, expanding the rest as what is known as the Born series:

For our particular case we take the homogeneous medium of 𝒪_h to have a constant permittivity α = α_r + iα_i. The inverse of this operator is given by the Green function for this medium, ${\mathit{𝒪}}_h^{-1}\equiv \overleftrightarrow{\mathit{G}}$, which is the two-index object (dyadic) that is the solution to

Osnabrugge et al. proved the convergence of the modified Born series for the scalar Helmholtz equation [17], the only proviso being that the system exhibits solely dissipation (no gain), and that the permittivity is nowhere singular. This was soon extended to the Born series (5) for Maxwell’s equations in isotropic materials [18]. To ensure convergence, two ingenious steps were originally taken by Osnabrugge et al. [17]. The first is to notice that the complex background permittivity α can be chosen freely, changing the eigenvalues of the operator $k_0^2\overleftrightarrow{\mathit{G}}\chi$, without affecting the solution. The second is that one can apply a pre-conditioner, Γ, to both sides of recursion relation (5) as follows: $\mathrm{\Gamma }\mathit{E}=k_0^2\mathrm{\Gamma }\overleftrightarrow{\mathit{G}}\chi \mathit{E}+\mathrm{\Gamma }\overleftrightarrow{\mathit{G}}\mathit{S}\to \mathit{E}=\left(k_0^2\mathrm{\Gamma }\overleftrightarrow{\mathit{G}}\chi +{\mathrm{𝟙}}_3-\mathrm{\Gamma }\right)\mathit{E}+\mathrm{\Gamma }\overleftrightarrow{\mathit{G}}\mathit{S}$.

The corresponding so-called ‘modified Born series’, which is the counterpart of series (5), is then given by

The paper is organized as follows. In Section 2 we describe how the modified Born series method can be generalized to anisotropic dielectrics. The iterative algorithm is introduced and we discuss how the preconditioner values can be calculated to ensure convergence. The algorithm is demonstrated with the polarization-dependent propagation through homogeneous and highly heterogeneous birefringent materials. Next, we show that the method is not limited to Hermitian permittivity matrix-functions and that it can be extended to non-Hermitian permittivity. In Section 3 we show how the preconditioning step of the algorithm can be adapted so that the same iteration can be used for magnetic materials. Accounting for electric-magnetic coupling also enables chiral, Tellegen [26], and general bi-anisotropic media. This is demonstrated with the propagation of a linearly polarized wave of visible light through 10 mm of a highly chiral substance. Finally we show how the method naturally handles more esoteric materials such as a solid with a negative refractive index, thus demonstrating the capability of simulating light propagation in such metamaterials.

2. A convergent Born series for anisotropic dielectrics

As a first step towards a Born series for arbitrary electromagnetic materials, we consider non-magnetic birefringent materials. Optical elements such as waveplates and Wollaston prisms are birefringent and characterized by anisotropic permittivity that cannot be represented by a scalar function, ∊(x). In this section we show how the modified Born series can be extended to permittivity tensors. Later it will be shown that arbitrary electromagnetic problems can be brought into the same form and solved using the iterative method introduced in this section.

To account for anisotropy, the scalar permittivity ∊(x) is replaced by the 3 × 3 matrix ∊(x) function in the spatial coordinates, x. In this more general case the constitutive relation for the electric displacement field (1) becomes

The question is whether there still exists a choice of α_i such that all the eigenvalues of M have a magnitude less than unity. To show that there is such a choice of α_i we consider the numerical radius, max_n |〈n|M|n〉|, where 〈n|n〉 = 1. Evidently the numerical radius will always be at least as large as the largest magnitude eigenvalue of M and thus we can enforce the convergence of the series (7) through the requirement

The dyadic Green function, G⃡, in convergence condition (11) can be written in terms of the identity and a unitary operator U as $\overleftrightarrow{\mathit{G}}\equiv \left(\mathit{U}-{\mathrm{𝟙}}_3\right)/\left(2\text{i}{k}_0^2{\alpha }_i\right)$ as discussed in Appendix 5.1.2. Condition (11) can then be rewritten using the triangle inequality |a + b| ≤ |a| + |b| as

For this particular bound it is very important that Δ_i is positive definite, rather than non-negative. If Δ_i has a kernel containing even one eigenvector |n₀〉 then for vectors |n〉 = |n₀〉 + η|n_⊥〉, A₂ diverges as ∼ 1/η, as η → 0. This is because A₂ is the analogue of a weak value [27] of Δ_r,i with respect to the two vectors |n〉 and Δ_i|n〉, a value which is well known to potentially lie outside the spectrum of the operator when the two vectors are close to orthogonal (so-called superweak values [28]). In Appendix 5.2.2 it is shown that when Δ_i has an empty kernel, no such divergence arises, and the Born series can be made to converge for any anisotropic medium. As it stands, it is difficult to make use of condition α_i > max{A₁, A₂}, because the value of A₂ cannot be readily estimated. Extensive numerical tests (App. 5.2.2) have also indicate the larges singular value ‖Δ‖ is a tighter bound for α_i that is valid for general gain-free materials, although we have not been able to show this analytically. As a diverging series is straightforward to detect, we have set α_i = ‖Δ‖ in all our simulations. Note that in the limiting case where the medium is isotropic and non-magnetic, Δ = Δ𝟙₃, and the basis |n〉 is the position basis, convergence condition α_i > max{A₁, A₂} reduces to α_i > max_x |Δ(x)|. This reproduces the convergence condition found by both Osnabrugge et al. [17] and Krüger et al. [18].

To demonstrate that this series can be used to calculate the propagation of light through an anisotropic material, we consider a birefringent crystal. Perhaps the most common example is calcite (CaCO₃, n_o = 2.776, n_e = 2.219 at λ = 500 nm), which splits an incident beam into two orthogonally-polarized beams that travel along different paths [29]. Fig. 1(a) shows that the modified Born series reproduces this effect for a circularly-polarized Gaussian beam (wavelength of 500 nm) traversing the crystal from left to right. The crystal displaces the extra-ordinary polarization laterally (e, magenta), while the ordinary-polarized component (o, cyan) travels along the original optical axis unaffected by the anisotropy.

 

Fig. 1 (a) Demonstration of anisotropic permittivity. Diagonally polarized light propagates from left to right through a calcite crystal (light gray box) cut at 45° with respect to its optic axis (indicated by the arrow). It can be seen that, as expected, the in-plane polarized extraordinary ray (e, magenta) is displaced from the ray that is polarized perpendicular to the plane (o, cyan). Some interference can be noticed between the incoming wave and its back-reflection at both the entrance and exit surface of the crystal. (b–e) A circularly polarized Gaussian beam incident from the left on a birefringent vaterite (CaCO₃) microrod with a diameter of 20 μm forms a complex scattering pattern instead of a single focus. Although the volume is homogeneous CaCO₃, complex, seemingly random, scattering occurs due to its subdivision in crystals of approximately 1 μm in cross-section for which the fast axis is oriented randomly with angles θ shown as hue in panel (e). Panels (b–d) show the field components E_x, E_y, and E_z, respectively. The darkness and hue indicate the field amplitude and phase, respectively, as indicated by the legend in panel (e). An overlaid gray grid outlines the crystal areas for reference, and the inset shows a 4× magnified detail of the field at the exit surface.

Download Full Size | PDF

This method of computing the solution to Maxwell’s equation is quite unusual compared to, for instance, a finite difference calculation. Firstly, each iteration can be performed in constant space and involves only a product of diagonal matrices and Fourier transforms, both of which have relatively low computational demands. Secondly, and perhaps even more interesting, the iterative corrections to the field predominantly depend on the local field. This is because the modified Green function, G⃡(x, x′), is the solution to Eq. (4), which is that for a dissipative medium characterized by a positive α_i. Both of these features can be advantageous for the efficient simulation of wave propagation in large heterogeneous materials. To illustrate this we simulated (λ₀ = 500 nm, 2D grid, pixel spacing 31.25 nm) the electric field within and behind a heterogeneous calcite rod of diameter 20 μm made up of crystal grains with a variable diameter of approximately 1 μm and with a random orientation (Figs. 1(b)–1(e)). An isotropic material would have focused both polarizations. Instead, a highly irregular speckle pattern can be seen to emerge, elongated along its propagation direction. The insets show close-ups of the high irregularity in the three field-components as the beam exits at the CaCO₃-air interface.

In the case of ordinary birefringence, the 3 × 3-matrix that represents the permittivity at a specific point in space is Hermitian. More general materials may have a permittivity that is non-Hermitian. To illustrate the application of the modified Born series to lossy anisotropic media we calculated the transmission through the simplest such system: an infinitely wide (one-dimensional, grid spacing Δz = 15.625 nm) absorbing polarizer. Fig. 2 shows a series of polarizers with varying alignment and an anisotropic extinction coefficient of κ = 0.1. For clarity, back reflections are avoided by setting the refractive index of the polarizers to the same value as that of the embedding medium (n = 1). The two geometries shown are two cross-polarizers (a,c,e), the second of which blocks all transmission. When a third polarizer is inserted in-between the same two polarizers, and with its transmission axis at 45° to the first two polarization axes (b,d,f), this results in a non-zero transmission in agreement with Malus’ law.

 

Fig. 2 Demonstration of non-Hermitian anisotropy. An E_x-polarized wave traverses two (a,c,e), or three (b,d,f), polarizers from the left to the right. The first and the last polarizer are in cross-diagonal-orientation, preventing transmission through the first system shown in (a). The intensity (black line) and the real part of the electric field (green line) are shown for the vertical and horizontal components in (c,e) and (d,f), respectively. Arbitrary units are used so that the maximum value of the intensity and field match for clarity of display.

Download Full Size | PDF

Algorithm 1 shows the pseudo code to calculate the electric field in anisotropic non-magnetic materials. Before starting the iteration loop, the algorithm must choose a background permittivity, α that will lead to a fast convergence. This is done by numerically finding the value α_r that minimizes the largest singular value of Δ = ∊(x) − α_r𝟙₃, using the Nelder-Mead algorithm [30]. The imaginary part of the background permittivity, α_i, is then set to equal the minimized upper bound for the singular value ‖Δ‖. Although the strict equality does not necessarily guarantee convergence, in practice we have never encountered divergence. It should be noted that calculating the largest singular value is impractical for large problems, so an upper bound must be calculated instead and α_i will likely be larger than the largest singular value, ‖Δ‖. For good measure, the complete algorithm (2) presented in Appendix 5.1, does include a check for divergence, though in practice it has never been encountered. The algorithms’ description is followed by a discussion of its efficiency, robustness to errors, and the issues of sampling and aliasing.

Algorithm 1. Calculation of the electric field in materials with anisotropic permittivity, ∊.

View Table | View all tables in this article

The iteration loop starts on line 6 and repeats until the condition on line 9 is met. In each loop a correction term, ΔE, is calculated for the electric field estimate, E. Unless divergence is detected, the correction term is added to the current electric field estimate. In the unlikely event that divergence is detected, the current iteration is repeated for a more conservative, larger, imaginary value for the permittivity bias, α_i. This procedure is repeated until the correction term becomes smaller than a preset tolerance, r_max. Alternative stopping criteria may be considered to improve the algorithm’s performance [18].

3. Magnetic and chiral materials

As it stands, Algorithm 1 can only be used on non-magnetic materials. However, the preconditioning step can be generalized to account for general magnetic media, without requiring changes to the iteration loop. A general algorithm is given in Appendix 5.1.1, which makes use of a generalized susceptibility, χ, to account for permeability and any electro-magnetic coupling constants.

The constitutive relations for general linear electromagnetic materials are given by (e.g. [31]):

Substituting Eqs. (16) and (17) into the Maxwell equations and eliminating the magnetic field, H = (iωµ₀)⁻¹μ⁻¹[∇ × E − ik₀ζE], leads to the same form of vector Helmholtz Eq. (9). Moreover, it can be noted that both sides of Eq. (9) can be divided by any positive constant, β ∈ ℝ_≥0. Eq. (9) thus remains valid for general materials when the source and susceptibility are generalized as

Although the modified Born series is now formally identical for these general constitutive relations (16) and (17), to that for the anisotropic dielectric discussed in Section 2, it is important to note that the susceptibility (19) is no longer a diagonal matrix in the position variable. This is because it now contains the differential operators 𝒟. Although this affects neither the form of the modified Born series (i.e. as a product of diagonal matrices and Fourier transforms), nor the proof of its convergence (which did not make any assumptions about the form of χ), it can make the convergence of the series less rapid. This is because the largest singular value of 𝒟 has a magnitude, π/(k₀Δr), which is inversely dependent on the grid spacing, Δr. Therefore, the required number of iterations increases quadratically with the density of the sampling grid. The impact on the convergence rate is highly dependent on the permeability distribution, and for magnetic materials with a highly variable permeability distribution such as metamaterials, the propagation of field corrections may be on the order of a wavelength or less per iteration. However, it follows from Eq. (19) that ‖𝟙₃ − μ⁻¹/β‖ can be minimized with an appropriate choice of the permeability scale, β, so to keep the convergence rate at a maximum.

In principle, the optimal value of α for general materials may be determined using a similar method as for dielectric case. Instead, albeit not critical to ensure convergence, the additional parameter β introduced in Eqs. (18) and (19) is leveraged to maximize the convergence rate for the magnetic case. Algorithm 2 in Appendix 5.1.1 generalizes Algorithm 1 to magnetic materials by redefining σ̂ as a function of two variables, α_r and β, and minimizing it in both variables simultaneously. Using the usual inequalities for matrix norms [35], the function σ̂(β, α_r) determines an upper bound for the singular values of the separate terms in the sum (19), and using the largest singular value, σ_𝒟 = π/(k₀Δr), of the discretized differential operator, 𝒟. The Nelder-Mead method is used to numerically minimize the product σ̂(β, α)β, which we found to correlate well with the inverse of the convergence rate [17]. The minimization provides values for the permeability scale, β, and the real part of the background permittivity α_r. The imaginary part, α_i, is set to σ̂(β, α), an upper bound estimate for the largest singular value of Δ = χ + iα_i𝟙₃.

Algorithm 2. A function that implements the general algorithm for arbitrary materials.

View Table | View all tables in this article

At optical frequencies, it is a common approximation to assume that the permeability, μ, equals μ₀𝟙₃ in the entire volume [34,36]. However, this precludes calculations with a broad class of materials of theoretical and practical interest. Maxwell’s equations are scale invariant so that the same solutions are obtained when the spatial and temporal frequencies, k and ω, are scaled by the same factor. This observation lies at the basis of impedance matching, transformation media [37], and many other interesting wave effects (e.g. [38]). In what follows we demonstrate the ability to accurately model variations in permeability with various examples.

By eliminating back reflections, impedance matching ensures efficient energy transfer and communication through waveguides. Figs. 3(a)–3(d) introduces the concept of impedance matching in one dimension. Figs. 3(a) and 3(b) show two samples with isotropic permittivity, ∊, and permeability, μ. Both samples contain objects (at 10 μm < z < 20 μm) with identical permittivity, larger than the surrounding medium. In the left-hand panel the permeability equals the background, while in the right-hand panel it equals the permittivity for all z. Fig. 3(b) thus represents an object that is impedance matched with the surrounding medium. As can been seen from Fig. 3(c), reflections from the front and back surface interfere with the incoming wave. This is most clearly visible in the ‘beating’ in the intensity where the incoming and the reflected wave interfere. Fig. 3(d) shows how impedance matching successfully suppresses back reflections at both interfaces. The absence of oscillations in the field amplitude indicates an absence of back reflections. Figs. 3(a) to 3(f) use a grid spacing of 32.25 nm.

 

Fig. 3 Demonstration of impedance matching (a–f), and propagation in a chiral medium (g,h). A plane wave in free space (∊ = μ = 1) with a wavelength λ₀ = 500 nm enters at x = 10 μm from the left into dielectric slabs of thickness 10 μm, with μ = 1 (a,c) and μ = 1.5 (b,d). In both cases the permittivity is 1.5 (green line, panels a and b). The interference between the incoming and reflected wave is clearly visible as oscillations in intensity (|E|², black line, c). It can be seen that a fraction of the wave is reflected from the slab without impedance matching (μ(z) ∝̸∊(z) in panel a). In contrast, a constant intensity is seen in panel (d), indicative of the absence of back-reflection for the impedance matched slab (μ ∝ ∊ in panel b). To facilitate comparison, both the intensity and field are normalized to their respective maximum value in panels (c) and (d). Panels (e,f) show the (truncated) electric field amplitude for a dipole with absorbing (e) and impedance matched (f) boundary layers. The interference with the back reflected wave, visible as beating in panel (e), is suppressed by the impedance matched layers as seen in panel (f). (g,h) Linear polarization rotates upon propagation in a chiral medium with a high chirality that is 100 times of that of saturated glucose (n = 1.45, specific rotation ${\left[\alpha \right]}_{500\hspace{0.17em}\text{nm}}^{T}52.7°\text{mL}\hspace{0.17em}{\text{g}}^{-1}{\text{dm}}^{-1}$, at 909 g/L). The constitutive relations can thus be seen to be $∊=\sqrt{1.45}$, μ = 1, and $\xi =\zeta =52.7\frac{909{\lambda }_0}{360}\text{i}=66.53\times {10}^{-6}\text{i}$. The intensity transfer between the E_x-polarization (solid green) and E_y-polarization (dashed red) can be seen to occur several times over a propagation distance of 10 mm (g). The local angle, θ, of the linear polarization is shown in panel (h). Note the significantly larger length scale for panels (g) and (h).

Download Full Size | PDF

Impedance matching also has important practical applications for simulations of infinite volumes in a finite space. Fig. 3(e) shows the electric dipole field in a box with regular absorbing boundaries and homogeneous permeability. As in the one-dimensional case, plotted in Figs. 3(a) and 3(c), significant reflections can be noted from the absorbing boundaries. Fig. 3(f) shows a dipole in a box with impedance matched absorbing boundaries. Although the boundaries both have a linearly varying permittivity from 1 to 1 + 1i and an identical thickness of 4 μm, it is clear from Fig. 3(f) that impedance matching suppresses the interfering reflections.

Since the constitutive relations (16) and (17) include the electric-magnetic coupling tensors, ζ and ξ, one can also use the modified Born series to treat bi-isotropic and more generally bi-anisotropic media. Many organic molecules such as glucose have a chiral asymmetry that leads to bi-isotropy. This chirality causes a rotation of linear polarization around the axis of propagation. As a second example of our generalized approach, Figs. 3(g) and 3(h) show simulations where linearly polarized light is slowly rotated in a chiral solution (Δz = 86.2 nm). As the rotation takes several millimeters to complete, the calculation was performed for a wavelength of 344.8 nm in glucose and over a propagation length of 10 mm.

As a final example we consider negative index materials. Simultaneously negative values for ∊ and μ can be obtained by engineering a material at the sub-wavelength scale. Such materials are characterized by a negative refractive index [39,40]. The general algorithm presented here naturally handles negative values for both the permittivity, ∊, and the permeability, μ. Fig. 4 demonstrates this for a block of transparent material, yet with a negative refractive index of −1.5 (∊ = −1.5, μ = −1; Δx = Δy = 125 nm). A Gaussian beam is incident at 30° to the surface normal and refracts backwards into the material with an angle of −19.5°, opposite to what is expected from Snell’s law for regular glass. The wave appears to travel backwards from left to right inside the metamaterial. At the exit surface the beam couples out at 30°, again opposing Snell’s law. Because there is no impedance matching in this case, some interference can be seen between the incident and reflected beam in Fig. 4.

 

Fig. 4 Light-wave propagation through a material with negative refractive index (∊ = −1.5, μ = −1, gray). The field amplitude (brightness) and phase (hue) are shown for a Gaussian beam that enters the surface at 30° to the normal. The beam can be seen to refract at −19.5°, backwards, into the metamaterial at an angle opposite to that for normal glass.

Download Full Size | PDF

4. Conclusion

We have extended the computationally efficient modified Born series [17,18], to arbitrary linear electromagnetic media. This enables the accurate calculation of light propagation through birefringent, chiral, and even magnetic materials. The general algorithm is listed in Algorithm 2, and implemented as a numerical library with working examples [25]. While there is no doubt that finite element simulations continue to be a very versatile and effective to invert the Helmholtz operator, the modified Born series has features that may make it more suitable for simulating complex samples as those found in microscopy. As the finite element method, it calculates the field due to a fixed frequency source; however, it does this through computing a series of terms that iteratively propagate the field out from the source, not entirely unlike a time-dependent calculation such as FDTD (see for example [17]). Yet, unlike an FDTD calculation, the solution at a given iteration is not considered a boundary condition for the next iteration, rather an estimate of the solution to Maxwell’s equations that is improved with each iteration. As such, numerical errors do not accumulate over time. The modified Born series—as it was originally intended [17]—is best suited for large scale problems with a limited range of susceptibilities. The modest memory requirements of the anisotropic modified Born iteration make it an ideal candidate to study light propagation in large heterogeneous tissues or materials (Appendix 5.1.3). Our extended method allows for a wide range of material properties, including birefringence, polarization, variable permeability, chirality, and negative refractive index. This enables it to account for the birefringence and chiral effects of light traversing biological tissue. The algorithm is also expected to find use in scattering studies, where laboratory experiments are often performed on highly scattering powders that typically consist of birefringent particles such as rutile powder (TiO₂). The presented numerical method offers a bridge between approximate analytic predictions and experiments. In addition, such large-scale problems are not exclusive to optics and occur across wave physics, including problems involving metamaterials, to which the modified Born series can now be applied.

5. Appendix

5.1. The algorithm

5.1.1. Detailed description of the general algorithm

The iteration loop of the algorithm is relatively short and identical for magnetic and non-magnetic materials (Algorithm 2, lines 15–23). Only the preconditioning steps differ. The algorithm starts by checking whether the material has magnetic properties and it determines the background permittivity, α = α_r + iα_i, and if required, the permeability scale, β. This enables the definition of the susceptibility, χ, the dyadic Green function, G⃡, and the source distribution, S.

Before starting the iteration loop, the electric field is either initialized to all zero, or to an approximate solution if available. On line 16, the loop starts by calculating the next term in the series, ΔE, using the operators χ and G⃡* on the source distribution, S, and the current estimate of the field, E. The next term is added to the current estimate, E, under the condition that the l₂-norm of the new term is less than that of the previous term. Otherwise, the series must be divergent for the current choice of the background permittivity, α, so its imaginary part is increased by 50%. The iteration continues until the updates to the field are deemed sufficiently small.

The susceptibility is defined by χ ≡ (∊ − ξµ⁻¹ζ)/β − 𝟙₃α − iξµ⁻¹𝒟/β + i𝒟μ⁻¹ζ/β + 𝒟 (𝟙₃ − μ⁻¹/β)𝒟, where $𝒟\equiv k_0^{-1}\mathbf{\nabla }\times =k_0^{-1}{ℱ}^{-1}\mathit{k}\otimes ℱ$, with ⊗ the outer product, while the forward and inverse Fourier transforms are represented by ℱ and ℱ⁻¹, respectively. Although all operations can be represented as large matrix operations, it is more space and time efficient to use fast-Fourier transforms and point-wise 3 × 3-dot products for each application of the operator χ on the electric field, E. The source is defined by $\mathit{S}\equiv \text{i}\omega {\mu }_0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{j}/\left(\beta k_0^2\right)$. The dyadic Green’s function, G⃡*, is discussed in Appendix 5.1.2 in what follows. Note that to prevent issues with numerical precision, in the implementation the factor $k_0^{-2}$ is moved from the definition of G⃡ to that of S.

5.1.2. The dyadic Green’s function

The dyadic Green function, G⃡, is integral to the calculation of the modified Born series:

(20)$$\mathit{E}=\left[\sum _{p=0}^{\infty }{\mathit{M}}^p\right]\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}\mathit{S},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{where}$$

(21)$$\mathit{M}\equiv k_0^2\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}\mathit{\chi }-\mathbf{\Gamma }+{\mathrm{𝟙}}_3,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}$$

(22)$$\mathbf{\Gamma }\equiv \frac{\text{i}}{{\alpha }_{\text{i}}}\mathit{\chi }.$$

Although well established in the literature of classical electromagnetism [41] and derived explicitly by Krüger et al. [18], the form of the dyadic Green function presented in Eq. (6) may not be familiar to the reader. Here we justify this equation, and show a useful representation of the Green function in terms of a unitary operator, U.

The vector Helmholtz equation can be separated into a homogeneous, 𝒪_h, and an inhomogeneous part, 𝒪_i, as follows

(23)$$\mathbf{\nabla }\times \mathbf{\nabla }\times \mathit{E}(\mathit{x})-k_0^2∊(\mathit{x})\mathit{E}(\mathit{x})=\text{i}\omega {\mu }_0\mathit{j}(\mathit{x})$$

(24)$$\left(\mathbf{\nabla }\times \mathbf{\nabla }\times -\alpha k_0^2\right)\mathit{E}-k_0^2(∊-\alpha )\mathit{E}=\mathit{S}$$

(25)$$\left({\mathit{𝒪}}_h+{\mathit{𝒪}}_i\right)\mathit{E}=\mathit{S}$$

The dyadic Green function, G⃡, is defined as the impulse response solution to the homogeneous part:

(26)$${\mathit{𝒪}}_h\overleftrightarrow{\mathit{G}}=\mathbf{\nabla }\times \mathbf{\nabla }\times \overleftrightarrow{\mathit{G}}(\mathit{x},{\mathit{x}}^{\prime })-k_0^2\alpha \overleftrightarrow{\mathit{G}}(\mathit{x},{\mathit{x}}^{\prime })={\mathrm{𝟙}}_3{\delta }^{(3)}(\mathit{x}-{\mathit{x}}^{\prime })$$

First, we take the Fourier transform of Eq. (26), recognizing that—due to the homogeneity of the medium—G⃡ must be a function of x − x′, and thus of a single variable k in Fourier space

(27)$$\overleftrightarrow{\mathit{G}}(\mathit{x},{\mathit{x}}^{\prime })=\int \frac{d^3\mathit{k}}{{\left(2\pi \right)}^3}\tilde{\mathit{G}}(\mathit{k}){\text{e}}^{\text{i}\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}$$

or equivalently in the operator notation used in the previous sections (where integrals are subsumed into the product)

(28)$$\overleftrightarrow{\mathit{G}}={ℱ}^{-1}\tilde{\mathit{G}}ℱ,$$

where G̃ is diagonal in Fourier space. Substituting expression (27) into (26) we obtain the following equation for the Fourier components of the Green function

(29)$$\mathit{k}\times \mathit{k}\times \tilde{\mathit{G}}(\mathit{k})+k_0^2\alpha \tilde{\mathit{G}}(\mathit{k})=-{\mathrm{𝟙}}_3.$$

At this point we decompose the identity matrix on the right into two parts 𝟙₃ = Π_T + Π_L, where

(30)$${\mathbf{\Pi }}_L=\frac{\mathit{k}\otimes \mathit{k}}{k^2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathbf{\Pi }}_T\equiv {\mathrm{𝟙}}_3-{\mathbf{\Pi }}_L,$$

and k ≡ ‖k‖ is the l₂-norm of k. These two operators can be physically understood as projecting out the longitudinal (Π_L) and transverse (Π_T) parts of the electromagnetic field, associated with electrostatic and radiative contributions respectively. We similarly decompose the Green function as G̃(k) = g_L(k)Π_L + g_T (k)Π_T, finding that Eq. (29) separates into two parts

(31)$$k_0^2\alpha g_L(\mathit{k})=-1$$

(32)$$\left(k^2-\alpha k_0^2\right)g_T(\mathit{k})=1,$$

from which we can deduce that the Green function is given by

(33)$$\overleftrightarrow{\mathit{G}}(\mathit{x},{\mathit{x}}^{\prime })=\int \frac{d^3\mathit{k}}{{\left(2\pi \right)}^3}\left[\frac{{\mathbf{\Pi }}_T}{k^2-\alpha k_0^2}-\frac{{\mathbf{\Pi }}_L}{\alpha k_0^2}\right]{\text{e}}^{\text{i}\mathit{k}(\mathit{x}-{\mathit{x}}^{\prime })}$$

which in our operator notation we write as

(34)$$\overleftrightarrow{\mathit{G}}={ℱ}^{-1}\left(\frac{{\mathbf{\Pi }}_T}{k^2-\alpha k_0^2}-\frac{{\mathbf{\Pi }}_L}{\alpha k_0^2}\right)ℱ$$

where the quantity in the rounded brackets must be understood as a matrix in both vector and Fourier indices. The dependence of the bracketed quantity on a single Fourier space variable k indicates that this matrix is diagonal in the Fourier indices, with each diagonal entry corresponding to a different value of k. This completes the demonstration of Eq. (34).

The Green function operator (34) can be written as a linear combination of the identity operator and a unitary operator, an observation that is instrumental in the study of the modified Born series (20) and its convergence. This representation of the Green function relies on the background permittivity α being chosen as a complex number. To find this representation we first note the following identity for any complex number z

(35)$$\frac{1}{z}=\frac{1}{2\text{i}\hspace{0.17em}\hspace{0.17em}\text{Im}\left[z\right]}\left(1-\frac{z^{*}}{z}\right)=\frac{1}{2\text{i}\hspace{0.17em}\hspace{0.17em}\text{Im}\left[z\right]}\left(1-{\text{e}}^{-2\text{i}\angle z}\right),$$

where ∠z indicates the complex argument of z. Applying identity (35) to (34), the Green function operator (34) is reduced to the following form

(36)$$\overleftrightarrow{\mathit{G}}=-\frac{1}{2\text{i}{\alpha }_{\text{i}}{k}_0^2}\left({\mathrm{𝟙}}_3-\mathit{U}\right)$$

where α_i is the imaginary part of α and the unitary matrix U is given by

(37)$$\mathit{U}={ℱ}^{-1}\left(\frac{k^2-{\alpha }^{*}{k}_0^2}{k^2-\alpha k_0^2}{\mathbf{\Pi }}_T+\frac{{\alpha }^{*}}{\alpha }{\mathbf{\Pi }}_L\right)ℱ,$$

where α^* is the complex conjugate of α. The operator (37) is unitary, due to the fact that our Fourier transform operators can be chosen such that ℱ^† = ℱ⁻¹, and thus

(38)$$\mathit{U}{\mathit{U}}^{†}={ℱ}^{-1}\left(\frac{k^2-{\alpha }^{*}{k}_0^2}{k^2-\alpha k_0^2}{\mathbf{\Pi }}_T+\frac{{\alpha }^{*}}{\alpha }{\mathbf{\Pi }}_L\right)ℱ{ℱ}^{-1}\left(\frac{k^2-\alpha k_0^2}{k^2-{\alpha }^{*}{k}_0^2}{\mathbf{\Pi }}_T+\frac{\alpha }{{\alpha }^{*}}{\mathbf{\Pi }}_L\right)ℱ={\mathrm{𝟙}}_3,$$

where we used the fact that Π_TΠ_L = 0, Π_TΠ_T = Π_T, and Π_LΠ_L = Π_L.

We note that, for numerical purposes, it is space and time-efficient to implement the dyadic Green function operation in k-space. The resulting multiplication only requires space to store the Fourier-transformed input, output, and Green function. The vector operations of the latter can be implemented without constructing the full, multi-dimensional, matrix for the dyadic Green function, though at least a scalar array, |k|², for normalizing the k-vectors must be stored.

5.1.3. Computation and memory efficiency

The computational efficiency of the modified Born series iteration has been discussed previously by Osnabrugge et al. [17]. In general the convergence of the iteration is approximately inversely proportional to the range of susceptibilities in the calculation volume. The method is therefore not efficient for calculations that involve metals. The anisotropic algorithm is no different in this respect. It is most efficient for heterogeneous dielectric materials such as biological tissue. However, the anisotropic version does require the addition of 3-vectors and multiplications of 3 × 3 matrices instead of scalar operations. The algorithm for anisotropic permittivity can thus be expected to be 9 times slower than the scalar wave algorithm [17], and 3 times slower than the isotropic vector algorithm [18]. It should also be noted that simulating inhomogeneous magnetic properties introduces the discretized differential operator, 𝒟. This largest singular value of this operator depends on the sampling density of the computation volume. This is equivalent to having a large variation in optical properties in the sample, which is known to lead to significantly slower convergence of the modified Born series [17].

The main limitation of any large scale electro-magnetic calculation is computer memory. The advantage of the presented algorithm is that the required memory scales with the number of sample points, P. This is important, considering that the calculation volumes of interest can be several orders of magnitude larger than the wavelength in all three dimensions. At a very minimum the electro-magnetic field in the calculation volume must be stored. The magnetic field can be calculated from the electric field, hence it is sufficient to store only the electric vector field, E. This can be done using 3P complex floating point numbers to represent a single frequency field in the calculation volume. Each iteration calculates a correction term for the field, ΔE, thus bringing the total to 6P.

Anisotropic permittivity can be represented by a 3 × 3 matrix for each point in space, to make up a block-diagonal matrix of dimension 3P × 3P. However, this can be efficiently stored using 9P complex numbers, while all matrix operations, as well as the fast-Fourier transforms can be performed in-place. It is not necessary to explicitly calculate the matrix for the dyadic Green’s function (App. 5.1.2), though it is necessary to determine the normalization factor ‖k‖ for all k-vectors. This requires storage for P values.

Anisotropy in the permeability, μ, or the coupling factors, ξ and ζ, will increase the memory requirements accordingly. If the source current S is not sparse, a further 3P complex values need to be stored. The memory requirements are summarized in Table 1 for the case of a spatially variant ∊, and optionally, a spatially variant μ, ξ, and ζ. The memory requirements are listed in different columns for the isotropic and the anisotropic cases. Mixed cases are omitted for brevity.

Table 1. Storage requirements for the algorithm.

View Table | View all tables in this article

5.1.4. Sampling and prevention of aliasing

The sampling grid size and density are important considerations for the calculation accuracy. In Algorithm 2, every iteration requires a multiplication by the preconditioned general susceptibility, χ, an addition of the source, S, a convolution with the dyadic Green’s function, and another multiplication by χ. Multiplications and additions cannot increase the spatial extent of the field E beyond the union of that of its left and right hand sides. However, the Green function has an infinite extent and therefore also the result of its convolution with a finite field. When the convolution operation is implemented as a multiplication in Fourier space, periodic boundary conditions are implicitly assumed. Alternative boundary conditions can be readily simulated by defining layers of non-transmitting material at the boundary. The sampling grid size must therefore be sufficiently large to fit both the volume of interest and multi-cell boundaries that adequately absorb the field.

To discuss the sampling density, we consider the Fourier transform of the iteration update calculation:

(39)$$\mathbf{\Delta }\tilde{\mathit{E}}=\frac{\text{i}}{{\alpha }_{\text{i}}}\tilde{\mathit{\chi }}*\left[\tilde{\mathit{G}}\left(\tilde{\mathit{\chi }}*\tilde{\mathit{E}}+\tilde{\mathit{S}}\right)-\tilde{\mathit{E}}\right].$$

In Fourier space, each step requires a convolution, *, with the Fourier transform of the susceptibility, χ̃, followed by an addition with the Fourier transform of the source, S̃, a multiplication with the Fourier transform of the dyadic Green’s function, G̃, and another convolution with χ̃. Although the multiplication by G̃, tends to suppress high spatial frequencies, it does not strictly limit the bandwidth support of the product. On the other hand, convolutions with χ̃ do extend the bandwidth support twice per iteration. If we define W_χ, W_E, and W_S ≤ W_χ + W_E as the spatial-frequency band-widths of χ, E, and S, respectively; it can seen that the bandwidth of the iteration update, ΔE, must be W_ΔE ≤ 2W_χ + W_E. Therefore, even when both the material properties and the source are smooth functions with a finite bandwidth, the sampling density of the calculation must, in principle, be increased by 2W_χ with every iteration step.

In practice, the calculation must be performed on a sample grid with finite bandwidth, W, and sample spacing, W⁻¹. With the notable exception of superoscillations, as long as the material properties and source field are spatially band-limited, the solution for the electric field can be expected to be concentrated around the Ewald sphere. We found that the smoothing effect of the Green function is generally sufficient to suppress the highest spatial frequencies that are affected by aliasing. When we accept the approximation that the solution must be band-limited, aliasing can be completely eliminated by low-pass filtering the field after each iteration step. To achieve this, the calculation must be performed with a sampling band-width W ≥ W_χ + W_E, in other words, with a sampling density that is no smaller than the sum of the Nyquist rates for E and χ. The two convolutions with χ̃ expand the support in frequency space to 2W_χ + W_E, thereby causing aliasing in the upper W_χ-part of the band. However, the lowest spatial frequencies in the W_E-band are not affected. All aliasing artefacts can thus be avoided by eliminating all but the spatial frequencies within the lower W_E-band of the iteration update, ΔE. Since the suppression of the highest spatial frequencies can be applied at the end of each iteration as a projection onto a subspace, it can be seen that convergence must also hold in this sub-space.

It should also be noted that the algorithm as it is presented here requires the material properties to be sampled on a regular grid. This enables efficient convolutions using the fast-Fourier transform and it simplifies the implementation in general. However, one could imagine structures that require a higher sampling density in specific regions. To address such need, we envision that the method can be extended to irregular grids by using non-uniform Fourier transforms to perform the convolutions [42,43].

5.1.5. Robustness to errors

Unlike time-stepped methods such as FDTD, the iteration presented here is robust to numerical errors. The N^th-correction term, ΔE_N, is obtained from recursive Eq. (20):

(40)$${\mathit{E}}_N=\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*k_0^2\mathit{\chi }{\mathit{E}}_{N-1}-\mathbf{\Gamma }{\mathit{E}}_{N-1}+{\mathit{E}}_{N-1}+\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*\mathit{S}$$

(41)$$\Rightarrow \mathrm{\Delta }{\mathit{E}}_N={\mathit{E}}_N-{\mathit{E}}_{N-1}=\mathbf{\Gamma }\left[\overleftrightarrow{\mathit{G}}*\left(k_0^2\mathit{\chi }{\mathit{E}}_{N-1}+\mathit{S}\right)-{\mathit{E}}_{N-1}\right].$$

For the initial field estimate, If we begin with a null field estimate E₀ = 0, this iteration is equivalent to the modified series (20). However, it should be noted that when E₀ ≠ 0, the iteration corresponds to a different series, yet one that converges under the same conditions and to the same limit:

(42)$${\mathit{E}}_N=\mathbf{\Gamma }\left[\overleftrightarrow{\mathit{G}}*\left(k_0^2\mathit{\chi }{\mathit{E}}_{N-1}+\mathit{S}\right)-{\mathit{E}}_{N-1}\right]+{\mathit{E}}_{N-1}$$

(43)$$=\mathit{M}{\mathit{E}}_{N-1}+\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*\mathit{S}$$

(44)$$={\mathit{M}}^N{\mathit{E}}_0+\left[\sum _{p=0}^{N-1}{\mathit{M}}^p\right]\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*\mathit{S}$$

It can be seen that the N^th-term of the series differs by M^N E₀ from that of series (20), and that its corresponding residue, ‖E − E_N‖, is

(45)$$\Vert \left[\sum _{p=0}^{\infty }{\mathit{M}}^p\right]\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*\mathit{S}-{\mathit{M}}^N{\mathit{E}}_0-\left[\sum _{p=0}^{N-1}{\mathit{M}}^p\right]\mathbf{\Gamma }\overleftrightarrow{\mathit{G}}*\mathit{S}\Vert =\Vert {\mathit{M}}^N(\mathit{E}-{\mathit{E}}_0)\Vert .$$

Provided that the absolute largest eigenvalue of M is smaller than one, the upper bound on the residue is tightened with every iteration, independently of the choice of the initial field. The corrections in consecutive terms prevent the accumulation of numerical errors, a common issue with techniques such as the finite-difference time-domain method.

5.2. Convergence of the modified Born series

5.2.1. Normal susceptibility matrices

The iteration as given by Eq. (20) calculates a series of correction terms of the form M^pE₀. Independently of the initial value E₀, this series is guaranteed to converge when all eigenvalues of M are less than 1 in absolute value. By substituting Eq. (36) in Eq. (21), the iteration operator can be written as a polynomial of the preconditioner, Γ, and a unitary transformation, U:

(46)$$\mathit{M}\equiv \mathbf{\Gamma }\left[\frac{1}{2}\left({\mathrm{𝟙}}_3-\mathit{U}\right)\frac{\text{i}}{{\alpha }_{\text{i}}}\mathit{\chi }-{\mathrm{𝟙}}_3\right]+{\mathrm{𝟙}}_3.$$

When the preconditioner is chosen to be $\mathbf{\Gamma }\equiv \frac{\text{i}}{{\alpha }_{\text{i}}}\mathit{\chi }={\mathrm{𝟙}}_3-\mathit{P}$, the expression simplifies to

(47)$$\mathit{M}\equiv \mathbf{\Gamma }\left[\frac{1}{2}\left({\mathrm{𝟙}}_3-\mathit{U}\right)\mathbf{\Gamma }-{\mathrm{𝟙}}_3\right]+{\mathrm{𝟙}}_3=\frac{1}{2}\left({\mathrm{𝟙}}_3+{\mathit{P}}^2\right)-\frac{1}{2}\mathbf{\Gamma }\mathit{U}\mathbf{\Gamma }.$$

The matrix $\mathit{P}\equiv {\mathrm{𝟙}}_3-\frac{\text{i}}{{\alpha }_{\text{i}}}\mathit{\chi }=-\frac{\text{i}}{{\alpha }_{\text{i}}}\mathbf{\Delta }$ can be seen to have a positive semidefinite real part for gain-free materials, and its largest singular value ‖P‖ < 1 when α_i is chosen to be larger than the largest singular value of Δ ≡ χ + iα_i. To prove that the eigenvalues of M are less than 1, it is sufficient to show that its numerical radius is less than 1. Using Cauchy-Schwartz to eliminate the unitary matrix, U, the numerical radius can be seen to be bound as follows

(48)$$\left|{\mathit{n}}^{†}\mathit{Mn}\right|\le \frac{1}{2}\left|1+{\mathit{n}}^{†}{\mathit{P}}^2\mathit{n}\right|+\frac{1}{2}\Vert {\mathbf{\Gamma }}^{†}\mathit{n}\Vert \Vert \mathbf{\Gamma }\mathit{n}\Vert ,$$

where Γ ≡ 𝟙₃ − P has the same eigenvectors as P and χ.

Often the susceptibility, χ, can be represented by a normal matrix, which is a matrix with orthonormal eigenvectors. This implies that the eigenvectors of the complex transpose P^† and Γ^† are also the same as those of χ. If we consider an eigenvector e_i of P with eigenvalue λ, then the corresponding eigenvalue of P^† is its complex conjugate λ^*, that of Γ is 1 − λ, and that of Γ^† is 1 − λ^*. When $\mathit{n}\equiv {\sum }_{i=1}^N{c}_i{\mathit{e}}_i$ is a linear combination of orthonormal eigenvectors e_i and complex coefficients c_i, the numerical radius can be seen to be bound as

(49)$$\left|{\mathit{n}}^{†}\mathit{Mn}\right|\le \sum _{i=1}^N{\left|c_i\right|}^2\left|{\mathit{e}}_i^{†}\mathit{M}{\mathit{e}}_i\right|\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\forall \mathit{n}:\mathit{n}\equiv \sum _{i=1}^{\mathit{N}}{c}_i{\mathit{e}}_i\wedge \Vert \mathit{n}\Vert =1.$$

Since ${\sum }_{i=1}^N{\left|c_i\right|}^2=1$ for orthonormal eigenvectors, to show that the convergence condition holds it is sufficient that we prove that $\left|{\mathit{e}}_i^{†}{\mathit{Me}}_i\right|<1$ for all eigenvectors, e_i.

Substituting Pe_i = λ, P^†e_i = λ^*, Γe_i = 1 − λ, and Γ^†e_i = 1 − λ^* in Eq. (48) enables us to rewrite the convergence condition as the scalar inequality:

(50)$$\left|{\mathit{e}}_i^{†}{\mathit{Me}}_i\right|\le \frac{\left|1+{\lambda }^2\right|}{2}+\frac{\left|1-{\lambda }^{*}\right|\left|1-\lambda \right|}{2}<1.$$

In terms of the real, λ_r, and imaginary, λ_i, parts of the eigenvalue λ ≡ λ_r + iλ_i, this yields:

(51)$$\frac{\sqrt{{\left|1+{\lambda }_{\text{r}}^2-{\lambda }_{\text{i}}^2\right|}^2+{\left|2{\lambda }_{\text{r}}{\lambda }_{\text{i}}\right|}^2}}{2}+\frac{1+{\lambda }_{\text{r}}^2+{\lambda }_{\text{i}}^2-2{\lambda }_{\text{r}}}{2}<1.$$

The square root can be eliminated by rearranging the terms and squaring both sides:

(52)$${\left|1+{\lambda }_{\text{r}}^2-{\lambda }_{\text{i}}^2\right|}^2+{\left|2{\lambda }_{\text{r}}{\lambda }_{\text{i}}\right|}^2<{\left(1-{\lambda }_{\text{r}}^2-{\lambda }_{\text{i}}^2+2{\lambda }_{\text{r}}\right)}^2.$$

For gain-free materials, the eigenvalues of P have a non-negative real part, λ_r ≥ 0. This allows us to rearrange the terms and simplifies the condition to:

(53)$$\left|1-{\lambda }_{\text{i}}^2\right|{\lambda }_{\text{r}}^2<\left(1-{\lambda }_{\text{i}}^2\right)\left(2{\lambda }_{\text{r}}-{\lambda }_{\text{r}}^2\right)+2{\lambda }_{\text{r}}^2\left(1-{\lambda }_{\text{i}}^2\right)-2{\lambda }_{\text{r}}^3.$$

Since the singular values of P are less than 1 (‖P‖ < 1), so must be its eigenvalues. Therefore, $1-{\lambda }_{\text{i}}^2\ge 0$, so we can subtract the left-hand from both sides and simplify the expression to get:

(54)$$0<\left(1-{\lambda }_{\text{i}}^2\right){\lambda }_{\text{r}}-{\lambda }_{\text{r}}^3.$$

It can be seen that the strict inequality does not hold for λ_r = 0. However, the condition does hold when there is a non-zero loss, λ_r > 0, which allows us to divide the expression by λ_r to yield ${\lambda }_{\text{i}}^2+{\lambda }_{\text{r}}^2<1$. Since the eigenvalues are less than 1, we have ${\lambda }_{\text{i}}^2+{\lambda }_{\text{r}}^2={\left|\lambda \right|}^2<1$, proving that $\left|{\mathit{e}}_i^{†}{\mathit{Me}}_i\right|\le 1$ for e_i any eigenvector of χ. Although convergence cannot be guaranteed in lossless materials, the strict inequality does not tend to be a limitation in practice since any non-zero loss at the boundary is sufficient to ensure convergence.

We thus conclude that the algorithm converges when the eigenvectors of the susceptibility distribution, χ, are orthogonal, a very common case. Yet, its eigenvectors will not be orthogonal when the reactive and dissipative parts of the susceptibility do not commute. This would occur when a birefringent crystal also has a polarization dependent absorption, yet with a different axis. In what follows, such more general susceptibilities are considered.

5.2.2. Numerical demonstration of the convergence

In Section 2 we showed that the following choice of α_i > max{A₁, A₂} guarantees convergence of the modified Born series, where A₁ and A₂ are defined by Eq. (15). It was noted that without constraining Δ_i to be positive definite, A₂ can be arbitrarily large, when 〈n|Δ_i|n〉 is close to zero. We begin this Section by showing that this divergence can be avoided so long as Δ_i has an empty kernel (no eigenvectors with zero eigenvalue). Suppose that Δ_i has an eigenvector |n₀〉 with eigenvalue λ, and consider that n can be written as the sum of two orthogonal vectors |n〉 = |n₀〉 + η|n_⊥〉 where η ≪ 1 and 〈n₀|n_⊥〉 = 0. The expression for A₂ inside the maximization is then given by

(55)$$\frac{\left|〈\mathit{n}\left|\mathbf{\Delta }{\mathbf{\Delta }}_{\text{i}}+{\mathbf{\Delta }}_{\text{i}}\mathbf{\Delta }\right|\mathit{n}〉\right|}{2〈\mathit{n}\left|{\mathbf{\Delta }}_{\text{i}}\right|\mathit{n}〉}=\frac{\left|\begin{array}{l}2\lambda 〈{\mathit{n}}_0\left|\mathbf{\Delta }\right|{\mathit{n}}_0〉+\eta \left(〈{\mathit{n}}_{\perp }\left|{\mathbf{\Delta }}_{\text{i}}\mathbf{\Delta }\right|{\mathit{n}}_0〉+〈{\mathit{n}}_0\left|\mathbf{\Delta }{\mathbf{\Delta }}_{\text{i}}\right|{\mathit{n}}_{\perp }〉\right)\\ +\eta \lambda \left(〈{\mathit{n}}_{\perp }\left|\mathbf{\Delta }\right|{\mathit{n}}_0〉+〈{\mathit{n}}_0\left|\mathbf{\Delta }\right|{\mathit{n}}_{\perp }〉\right)+{\eta }^2〈{\mathit{n}}_{\perp }\left|\mathbf{\Delta }{\mathbf{\Delta }}_{\text{i}}+{\mathbf{\Delta }}_{\text{i}}\mathbf{\Delta }\right|{\mathit{n}}_{\perp }〉|\end{array}\right|}{2\left(\lambda +{\eta }^2〈{\mathit{n}}_{\perp }\left|{\mathbf{\Delta }}_{\text{i}}\right|{\mathit{n}}_{\perp }〉\right)}$$

If λ = 0 (and the kernel of Δ_i thus contains |n₀〉 then the above quantity diverges as 1/η as we take η to zero. The quantity A₂ thus becomes infinite. However, if λ is non-zero (but arbitrarily small) then as η → 0, (55) tends to |〈n₀|Δ|n₀〉| which is finite, and we note, smaller than or equal to the largest singular value of Δ.

Assuming that 〈n|Δ_i|n〉 is never zero, the quantity A₂ can be rewritten in a form that is easier to compute. First, we write the positive definite Hermitian matrix, Δ_i, as the square of another Hermitian matrix, a, so that it equals: Δ_i = a². Defining |n′〉 = a|n〉, the quantity A₂ can be written as

(56)$$A_2=\underset{{\mathit{n}}^{\prime }}{\text{max}}\frac{\left|〈{\mathit{n}}^{\prime }\left|{\mathit{a}}^{-1}\mathbf{\Delta }\mathit{a}+\mathit{a}\mathbf{\Delta }{\mathit{a}}^{-1}\right|{\mathit{n}}^{\prime }〉\right|}{2〈{\mathit{n}}^{\prime }|{\mathit{n}}^{\prime }〉}$$

which is simply the numerical radius r(m) of the matrix

(57)$$\mathit{m}=\frac{1}{2}\left({\mathit{a}}^{-1}\mathbf{\Delta }\mathit{a}+\mathit{a}\mathbf{\Delta }{\mathit{a}}^{-1}\right).$$

The numerical radius of a matrix is always less than or equal to its norm [45] r(m) ≤ ‖m‖. In addition, the norm of a sum is never larger than the sum of the norms, and we can therefore estimate A₂ as

(58)$$A_2=\frac{1}{2}\left[\Vert {\mathit{a}}^{-1}\mathbf{\Delta }\mathit{a}\Vert +\Vert \mathit{a}\mathbf{\Delta }{\mathit{a}}^{-1}\Vert \right]$$

where ‖ · ‖ indicates the l₂-norm. While the spectra of Δ and e.g. aΔa⁻¹ are the same, their norms are not. Given that the spectral radius of an operator is always less than or equal to its norm, the lowest possible value of A₂ is the magnitude of the largest eigenvalue |λ_max| of Δ. Meanwhile, A₁, being the square root of the numerical radius of a Hermitian operator, is equal to the square root of the operator norm (1/2)‖ΔΔ^† + Δ^†Δ‖, which is bounded by

(59)$$A_1\le \frac{1}{2}\Vert \mathbf{\Delta }{\mathbf{\Delta }}^{†}\Vert +\frac{1}{2}\Vert {\mathbf{\Delta }}^{†}\mathbf{\Delta }\Vert =\left|{\lambda }_{\text{max}}\right|\le A_2$$

where the estimate for A₂ is here given by Eq. (58). We can therefore use this estimate of A₂ to find an upper bound for the value of α_i. Fig. 5(b) shows a numerical test, where we chose A₂ according to (58) and repeatedly evaluated the inner product 〈n|M|n〉 for different choices of random complex 40 × 40 matrices Δ and U. These tests represent 10⁶ × 40 = 4 × 10⁷ evaluations of the inner product and appear to indicate that α_i > ‖Δ‖ is a tighter bound (Fig. 5(a)).

 

Fig. 5 Numerical check that condition α_i > max {A₁, A₂} ensures that the numerical radius of M, given by Eq. (21), is less than unity. To produce this figure we used the condition on the numerical radius, |〈n|M|n〉| < 〈n|n〉, with 40 × 40 matrices. We generated 1 × 10⁶ random matrices Δ (positive definite Δ_i) and unitary matrices U, along with the corresponding values of α_i calculated as indicated in each panel (using Python’s numpy library for random matrix generation [44]). For each matrix we calculated values for |〈n|M|n〉| for a set of 40 random but orthogonal complex vectors, n. The largest magnitude of these values is one of the 10⁶ points plotted in each panel. In (b) (largest magnitude 0.999997) we used the condition (58), which we know analytically to be sufficient to move the eigenvalues of M within the unit circle. In (a) (largest magnitude 0.999999) we show that α_i > ‖Δ‖ appears to guarantee also that the eigenvalues are within the unit circle

Download Full Size | PDF

Funding

Philip Leverhulme Prize (PLP-2015-216); Royal Society and TATA (RPG-2016-186).

References

1. O. G. Ernst and M. J. Gander, “Why it is difficult to solve Helmholtz problems with classical iterative methods,” in Numerical Analysis of Multiscale Problems, I. G. Graham, T. Y. Hou, O. Lakkis, and R. Scheichl, eds. (Springer, 2012).

2. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light - Second Edition (Princeton University, 2008).

3. H. C. van de Hulst, Light Scattering by Small Particles (Dover Publications Inc., 2003).

4. S. Rotter and S. Gigan, “Light fields in complex media: Mesoscopic scattering meets wave control,” Rev. Mod. Phys. 89, 015005 (2017). [CrossRef]  

5. I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. 32(16), 2309–2311 (2007). [CrossRef]   [PubMed]  

6. I. M. Vellekoop, A. Lagendijk, and A. P. Mosk, “Exploiting disorder for perfect focusing,” Nat. Photonics 4, 320–322 (2010). [CrossRef]  

7. A. K. Glaser, Y. Chen, and J. T. Liu, “Fractal propagation method enables realistic optical microscopy simulations in biological tissues,” Optica 3(8), 861–869 (2016). [CrossRef]   [PubMed]  

8. Z. Salamon and G. Tollin, “Optical anisotropy in lipid bilayer membranes: coupled plasmon-waveguide resonance measurements of molecular orientation, polarizability, and shape,” Biophys. J. 80(3), 1557–1567 (2001). [CrossRef]   [PubMed]  

9. M. Everett, K. Schoenenberger, B. Colston, and L. Da Silva, “Birefringence characterization of biological tissue by use of optical coherence tomography,” Opt. Lett. 23(3), 228–230 (1998). [CrossRef]  

10. W.-C. Kuo, S.-C. Lin, and C.-Y. Chuang, “Birefringence measurement in polarization-sensitive optical coherence tomography using differential-envelope detection method,” Rev. Sci. Instrum. 81, 053705 (2010). [CrossRef]   [PubMed]  

11. G. Jarry, F. Henry, and R. Kaiser, “Anisotropy and multiple scattering in thick mammalian tissues,” J. Opt. Soc. Am. A 17(1), 149–153 (2000). [CrossRef]  

12. J. F. De Boer, T. E. Milner, M. J. van Gemert, and J. S. Nelson, “Two-dimensional birefringence imaging in biological tissue by polarization-sensitive optical coherence tomography,” Opt. Lett. 22(12), 934–936 (1997). [CrossRef]   [PubMed]  

13. S. Alali, K. J. Aitken, A. Schröder, D. J. Bagli, and I. A. Vitkin, “Optical assessment of tissue anisotropy in ex vivo distended rat bladders,” J. Biomed. Opt. 17(8), 0860101 (2012). [CrossRef]  

14. A. Pierangelo, A. Nazac, A. Benali, P. Validire, H. Cohen, T. Novikova, B. H. Ibrahim, S. Manhas, C. Fallet, M.-R. Antonelli, and A.-D. Martino, “Polarimetric imaging of uterine cervix: a case study,” Opt. Express 21(12), 14120–14130 (2013). [CrossRef]   [PubMed]  

15. M. Koike-Tani, T. Tani, S. B. Mehta, A. Verma, and R. Oldenbourgh, “Polarized light microscopy in reproductive and developmental biology,” Mol. Reprod. Dev. 82(7–8), 548–562 (2013). [CrossRef]   [PubMed]  

16. D. C. Adams, L. P. Hariri, A. J. Miller, Y. Wang, J. L. Cho, M. Villiger, J. A. Holz, M. V. Szabari, D. L. Hamilos, R. Scott Harris, J. W. Griffith, B. E. Bouma, A. D. Luster, B. D. Medoff, and M. J. Suter, “Birefringence microscopy platform for assessing airway smooth muscle structure and function in vivo,” Sci. Transl. Med. 8(359), 359ra131 (2016). [CrossRef]   [PubMed]  

17. 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). [CrossRef]  

18. B. Krüger, T. Brenner, and A. Kienle, “Solution of the inhomogeneous Maxwell’s equations using a Born series,” Opt. Express 25(21), 25165–25182 (2017). [CrossRef]   [PubMed]  

19. I. Koutromanos, Fundamentals of Finite Element Analysis: Linear Finite Element Analysis (John Wiley & Sons, 2018).

20. M. Born, “Quantenmechanik der stoßvorgänge,” Z. Phys. 38, 803 (1926). [CrossRef]  

21. E. Akkermans and G. Montambaux, Mesoscopic Physics of Electrons and Photons (Cambridge University, 2007), chap. Coherent backscattering of light. [CrossRef]  

22. R. G. Newton, Scattering Theory of Waves and Particles (Springer, 2013).

23. S. A. R. Horsley, M. Artoni, and G. C. L. Rocca, “Spatial kramers-kronig relations and the reflection of waves,” Nat. Photonics 9, 436 (2015). [CrossRef]  

24. J. S. Walker, Fast Fourier Transforms (CRC Press, 2017).

25. T. Vettenburg, “Macromax: A python 3 library for solving macroscopic maxwell’s equations for electromagnetic waves in gain-free heterogeneous (bi-)(an)isotropic (non)magnetic materials,” figshare (2019). [Retrieved 17 Jan 2019], https://doi.org/10.6084/m9.figshare.7600100.

26. B. D. H. Tellegen, “The gyrator, a new electric network element,” Philips Res. Rep. 3(2), 81–101 (1948).

27. Y. Aharonov, D. Z. Albert, and L. Vaidman, “How the result of a measurement of a component of the spin of a spin-1/2 particle can turn out to be 100,” Phys. Rev. Lett. 60(14), 1351–1354 (1988). [CrossRef]   [PubMed]  

28. M. V. Berry and P. Shukla, “Typical weak and superweak values,” J. Phys. A 43, 354024 (2010). [CrossRef]  

29. G. Waldman, Introduction to Light: The Physics of Light, Vision, and Color (Dover, 2002).

30. J. A. Nelder and R. Mead, “A simplex method for function minimization,” Comput. J. 7(4), 308–313 (1965). [CrossRef]  

31. A. Priou, A. Sihvola, S. Tretyakov, and A. Vinogradov, Advances in Complex Electromagnetic Materials (Springer, 2012).

32. D. R. Smith, “Calculation and measurement of bianisotropy in a split ring resonator metamaterial,” J. Appl. Phys. 100, 024507 (2006). [CrossRef]  

33. A. Lakhtakia, Beltrami Fields in Chiral Media (World Scientific, 1994). [CrossRef]  

34. L. D. Landau, L. P. Pitaevskii, and E. M. Lifshitz, Electrodynamics of Continuous Media (Butterworth-Heinemann, 2004).

35. I. S. Gradsheyn and I. M. Ryzhik, Table of Integrals Series and Products (Academic Press, 2000).

36. R. Merlin, “Metamaterials and the Landau–Lifshitz permeability argument: Large permittivity begets high-frequency magnetism,” Proc. Natl. Acad. Sci. U.S.A. 106(6), 1693–1698 (2009). [CrossRef]   [PubMed]  

37. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science. 312(5781), 1780–1782 (2006). [CrossRef]   [PubMed]  

38. B. Vial, Y. Liu, S. A. R. Horsley, T. G. Philbin, and Y. Hao, “A class of invisible inhomogeneous media and the control of electromagnetic waves,” Phys. Rev. B 94, 245119 (2016). [CrossRef]  

39. V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of ∊ and μ,” Sov. Phys. Uspekhi 10(4), 509–514 (1967). [CrossRef]  

40. J. B. Pendry, “Negative refraction makes a perfect lens,” Phys. Rev. Lett. 85(18), 3966–3969 (2000). [CrossRef]   [PubMed]  

41. B. H. L. Novotny, Principles of Nano-Optics (Cambridge University, 2012). [CrossRef]  

42. A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14(6), 1368–1393 (1993). [CrossRef]  

43. J. A. Fessler and B. P. Sutton, “Nonuniform fast fourier transforms using min-max interpolation,” IEEE Trans. Signal Process. 51(2), 560–574 (2003). [CrossRef]  

44. T. E. Oliphant, A guide to NumPy (Trelgol Publishing, 2006).

45. S. S. Dragomir, Inequalities for the Numerical Radius of Linear Operators in Hilbert Spaces (Springer, 2013). [CrossRef]