Calculating coherent light-wave propagation in large heterogeneous media

TL;DR

This paper introduces an efficient computational method for modeling light-wave propagation in large, complex, and anisotropic media, enabling detailed analysis of scattering and birefringence in biological tissues.

Contribution

The paper presents a novel algorithm that accurately solves Maxwell's equations in arbitrary anisotropic materials, facilitating large-scale simulations of light propagation in complex media.

Findings

Enables detailed modeling of light in birefringent and chiral materials

Supports calculations in media with negative refractive index

Improves understanding of light-tissue interactions for microscopy

Abstract

Understanding the interaction of light with a highly scattering material is essential for optical microscopy of optically thick and heterogeneous biological tissues. Ensemble-averaged analytic solutions cannot provide more than general predictions for relatively simple cases. Yet, biological tissues contain chiral organic molecules and many of the cells' structures are birefringent, a property exploited by polarization microscopy for label-free imaging. Solving Maxwell's equations in such materials is a notoriously hard problem. Here we present an efficient method to determine the propagation of electro-magnetic waves in arbitrary anisotropic materials. We demonstrate how the algorithm enables large scale calculations of the scattered light field in complex birefringent materials, chiral media, and even materials with a negative refractive index.

Tables

Table 1. Table S1: Storage requirements for the algorithm.

spatial variability	isotropic	anisotropic
$ϵ$	$11P$	$19P$
$ϵ$, and $𝝁$	$12P$	$28P$
$ϵ$, $𝝁$, $𝝃$, and $𝜻$	$14P$	$46P$

Equations

$$\mathbf{E}$$

$$\mathbf{M}$$

$$\mathbf{\Gamma }$$

$$\nabla \times \nabla \times \mathbf{E}(\mathbf{x})-k_0^2ϵ(\mathbf{x})\mathbf{E}(\mathbf{x})$$

$$\left(\nabla \times \nabla \times -\alpha k_0^2\right)\mathbf{E}-k_0^2(ϵ-\alpha )\mathbf{E}$$

$$\left({\mathcal{O}}_h+{\mathcal{O}}_i\right)\mathbf{E}$$

$${\mathcal{O}}_h\overleftrightarrow{G}=\nabla \times \nabla \times \overleftrightarrow{G}(\mathbf{x},{\mathbf{x}}^{\prime })-k_0^2\alpha \overleftrightarrow{G}(\mathbf{x},{\mathbf{x}}^{\prime })=1_3{\delta }^{(3)}(\mathbf{x}-{\mathbf{x}}^{\prime })$$

$$\overleftrightarrow{G}(\mathbf{x},{\mathbf{x}}^{\prime })=\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}\tilde{\mathbf{G}}(\mathbf{k}){\mathrm{e}}^{\mathrm{i}\mathbf{k}(\mathbf{x}-{\mathbf{x}}^{\prime })}$$

$$\overleftrightarrow{G}={\mathcal{F}}^{-1}\tilde{\mathbf{G}}\mathcal{F},$$

$$\mathbf{k}\times \mathbf{k}\times \tilde{\mathbf{G}}(\mathbf{k})+k_0^2\alpha \tilde{\mathbf{G}}(\mathbf{k})=-1_3.$$

$${\mathbf{\Pi }}_L=\frac{\mathbf{k}\otimes \mathbf{k}}{{\mathbf{k}}^2}.$$

$$k_0^2\alpha g_L(\mathbf{k})$$

$$\left({\mathbf{k}}^2-\alpha k_0^2\right)g_T(\mathbf{k})$$

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

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

$$\frac{1}{z}=\frac{1}{2\mathrm{i}\mathrm{Im}[z]}\left(1-\frac{z^{\star }}{z}\right)=\frac{1}{2\mathrm{i}\mathrm{Im}[z]}\left(1-{\mathrm{e}}^{-2\mathrm{i}\phase{z}}\right),$$

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

$$\mathbf{U}={\mathcal{F}}^{-1}\left(\frac{k^2-{\alpha }^{\star }{k}_0^2}{k^2-\alpha k_0^2}{\mathbf{\Pi }}_T+\frac{{\alpha }^{\star }}{\alpha }{\mathbf{\Pi }}_L\right)\mathcal{F},$$

$$\begin{gathered} \mathbf{U}{\mathbf{U}}^{†}= \\ {\mathcal{F}}^{-1}\left(\frac{k^2-{\alpha }^{\star }{k}_0^2}{k^2-\alpha k_0^2}{\mathbf{\Pi }}_T+\frac{{\alpha }^{\star }}{\alpha }{\mathbf{\Pi }}_L\right)\left(\frac{k^2-\alpha k_0^2}{k^2-{\alpha }^{\star }{k}_0^2}{\mathbf{\Pi }}_T+\frac{\alpha }{{\alpha }^{\star }}{\mathbf{\Pi }}_L\right)\mathcal{F} \\ =\text{mathbbold}{1}_3, \end{gathered}$$

$$\mathbf{\Delta }\tilde{\mathbf{E}}$$

$${\mathbf{E}}_N$$

$$\Rightarrow \Delta {\mathbf{E}}_N$$

$${\mathbf{E}}_N$$

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

$$\overleftrightarrow{G}$$

$$\mathbf{M}$$

$$\mathrm{Re}\left\{{\mathbf{E}}_{\max }^{†}\mathbf{ϵ}{\mathbf{E}}_{\max }\right\}$$

$$\mathrm{Im}\left\{{\mathbf{E}}_{\max }^{†}\mathbf{ϵ}{\mathbf{E}}_{\max }\right\}$$

$$\begin{gathered} {\mathbf{E}}_{\max }^{†}\left[\mathbf{\Gamma }\left(\frac{\text{mathbbold}{1}_3-\mathbf{U}}{2}\frac{\mathrm{i}}{{\alpha }_{\mathrm{i}}}\mathbf{\chi }-\text{mathbbold}{1}_3\right)+\text{mathbbold}{1}_3\right]{\mathbf{E}}_{\max }= \\ {\mathbf{E}}_{\max }^{†}\left[\mathbf{\Gamma }\left(\frac{1}{2}\frac{\mathrm{i}}{{\alpha }_{\mathrm{i}}}\mathbf{\chi }-\text{mathbbold}{1}_3\right)+\text{mathbbold}{1}_3\right]{\mathbf{E}}_{\max }-{\mathbf{E}}_{\max }^{†}\mathbf{\Gamma }\frac{\mathbf{U}}{2}\frac{\mathrm{i}}{{\alpha }_{\mathrm{i}}}\mathbf{\chi }{\mathbf{E}}_{\max }. \end{gathered}$$

$$\begin{gathered} {\mathbf{E}}_{\max }^{†}\mathbf{M}{\mathbf{E}}_{\max }= \\ \sum _i|c_i{|}^2\left[\frac{{\lambda }_i^2}{2}-{\lambda }_i+1\right]-\frac{1}{2}{\left({\mathbf{\Gamma }}^{†}{\mathbf{E}}_{\max }\right)}^{†}\mathbf{U}\frac{\mathrm{i}}{{\alpha }_{\mathrm{i}}}\mathbf{\chi }{\mathbf{E}}_{\max }. \end{gathered}$$

$$\begin{gathered} \left|\sum _i|c_i{|}^2\left[\frac{{\lambda }_i^2}{2}-{\lambda }_i+1\right]\right|\le \sum _i|c_i{|}^2\left[1-\frac{|{\lambda }_i{|}^2}{2}\right]= \\ \sum _i|c_i{|}^2-\sum _i|c_i{|}^2\frac{|{\lambda }_i{|}^2}{2}=1-\frac{1}{2}\sum _i|c_i{|}^2|{\lambda }_i{|}^2. \end{gathered}$$

$${\alpha }_{\mathrm{i}}>\max \left\{A_1,A_2\right\},$$

$$A_1$$

$$A_2$$

Full text

**Supplementary Information:

Calculating coherent light-wave propagation

in large heterogeneous media.**

S1   The algorithm

S1A   Detailed description of the general algorithm

The iteration loop of the algorithm is relatively short and identical for magnetic and non-magnetic materials (Algorithm S1, lines 16-24). Only the preconditioning steps differ. The algorithm starts by checking whether the material has magnetic properties and it determines the background permittivity, $\alpha=\alpha_{\mathrm{r}}+\mathrm{i}\alpha_{\mathrm{i}}$, and if required, the permeability scale, $\beta$. This enables the definition of the susceptibility, $\boldsymbol{\chi}$, the dyadic Green function, \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}, and the source distribution, $\boldsymbol{S}$.

Before starting the iteration loop, the electric field is either initialized to all zero, or to an approximate solution if available. On line 17, the loop starts by calculating the next term in the series, $\Delta\boldsymbol{E}$, using the operators $\boldsymbol{\chi}$ and \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}\ast on the source distribution, $\boldsymbol{S}$, and the current estimate of the field, $\boldsymbol{E}$. The next term is added to the current estimate, $\boldsymbol{E}$, under the condition that the l2-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, $\alpha$, 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 $\boldsymbol{\chi}\equiv\left(\boldsymbol{\epsilon}-\boldsymbol{\xi}\boldsymbol{\mu}^{-1}\boldsymbol{\zeta}\right)/\beta-\mathbb{1}_{3}\alpha-\mathrm{i}\boldsymbol{\xi}\boldsymbol{\mu}^{-1}\boldsymbol{\mathcal{D}}/\beta+\mathrm{i}\boldsymbol{\mathcal{D}}\boldsymbol{\mu}^{-1}\boldsymbol{\zeta}/\beta+\boldsymbol{\mathcal{D}}\left(\mathbb{1}_{3}-\boldsymbol{\mu}^{-1}/\beta\right)\boldsymbol{\mathcal{D}}$, where $\boldsymbol{\mathcal{D}}\equiv k_{0}^{-1}\boldsymbol{\nabla}\times=k_{0}^{-1}\mathcal{F}^{-1}\boldsymbol{k}\otimes\mathcal{F}$, with $\otimes$ the outer product, while the forward and inverse Fourier transforms are represented by $\mathcal{F}$ and $\mathcal{F}^{-1}$, 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\times 3$-dot products for each application of the operator $\boldsymbol{\chi}$ on the electric field, $\boldsymbol{E}$. The source is defined by $\boldsymbol{S}\equiv\mathrm{i}\omega\mu_{0}\;\boldsymbol{j}/\left(\beta k_{0}^{2}\right)$. The dyadic Green’s function, \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}\ast, is discussed in section S1B 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 \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}} to that of $\boldsymbol{S}$.

S1B   The dyadic Green’s function

The dyadic Green function, \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}, is integral to the calculation of the modified Born series:

[TABLE]

Although well established in the literature of classical electromagnetism 1, the form of the dyadic Green function that we have presented in equation (6) of the main text may not be familiar to the reader. In this appendix we justify this equation, and show a useful representation of the Green function in terms of a unitary operator, $\boldsymbol{U}$.

The vector Helmholtz equation can be separated into a homogeneous, $\boldsymbol{\mathcal{O}}_{h}$, and an inhomogeneous part, $\boldsymbol{\mathcal{O}}_{i}$, as follows

[TABLE]

The dyadic Green function, \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}, is defined as the impulse response solution to the homogeneous part:

[TABLE]

First, we take the Fourier transform of equation (S5), recognizing that—due to the homogeneity of the medium—\overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}} must be a function of $\boldsymbol{x}-\boldsymbol{x}^{\prime}$, and thus of a single variable $\boldsymbol{k}$ in Fourier space

[TABLE]

or equivalently in the operator notation used in the main text (where integrals are subsumed into our ‘’ product)

[TABLE]

where $\tilde{\boldsymbol{G}}$ is diagonal in Fourier space. Substituting expression (S6) into (S5) we obtain the following equation for the Fourier components of the Green function

[TABLE]

At this point we decompose the identity matrix on the right into two parts $\mathbb{1}_{3}=(\mathbb{1}_{3}-\boldsymbol{\Pi}_{L})+\boldsymbol{\Pi}_{L}=\boldsymbol{\Pi}_{T}+\boldsymbol{\Pi}_{L}$, where

[TABLE]

These two operators can be physically understood as projecting out the longitudinal ($\boldsymbol{\Pi}_{L}$) and transverse ($\boldsymbol{\Pi}_{T}$) parts of the electromagnetic field, associated with electrostatic and radiative contributions respectively. We similarly decompose the Green function as $\tilde{\boldsymbol{G}}(\boldsymbol{k})=g_{L}(\boldsymbol{k})\boldsymbol{\Pi}_{L}+g_{T}(\boldsymbol{k})\boldsymbol{\Pi}_{T}$, finding that equation (S8) separates into two parts

[TABLE]

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

[TABLE]

which in our operator notation we write as

[TABLE]

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 $\boldsymbol{k}$ indicates that this matrix is diagonal in the Fourier indices, with each diagonal entry corresponding to a different value of $\boldsymbol{k}$. This completes the demonstration of equation (S12).

The Green function operator (S12) 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 (S1) and its convergence. This representation of the Green function relies on the background permittivity $\alpha$ being chosen as a complex number. To find this representation we first note the following identity for any complex number $z$

[TABLE]

where $\phase{\,z}$ indicates the complex argument of $z$. Applying identity (S13) to (S12), the Green function operator (S12) is reduced to the following form

[TABLE]

where $\alpha_{\mathrm{i}}$ is the imaginary part of $\alpha$ and the unitary matrix $\boldsymbol{U}$ is given by

[TABLE]

where $k=\|\boldsymbol{k}\|$ is the $l_{2}$-norm of $\boldsymbol{k}$, while $\alpha^{\ast}$ is the complex conjugate of $\alpha$. The operator (S15) is unitary, due to the fact that our Fourier transform operators can be chosen such that $\mathcal{F}^{\dagger}=\mathcal{F}^{-1}$, and thus

[TABLE]

where we used the fact that $\boldsymbol{\Pi}_{T}\boldsymbol{\Pi}_{L}=0$, $\boldsymbol{\Pi}_{T}\boldsymbol{\Pi}_{T}=\boldsymbol{\Pi}_{T}$, and $\boldsymbol{\Pi}_{L}\boldsymbol{\Pi}_{L}=\boldsymbol{\Pi}_{L}$.

Note that the unitary operator only changes the complex argument of the transverse and longitudinal projected components. It can thus be seen from equation (S14) that the eigenvalue of \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}} with the largest absolute value equals $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}k_{0}^{2}}$, and that the real part of the eigenvalue must be no larger than $\frac{1}{2\alpha_{\mathrm{i}}k_{0}^{2}}$ in absolute value, while its imaginary part cannot be negative. More specifically, all eigenvalues are contained within a disk of radius $\frac{1}{2\alpha_{\mathrm{i}}k_{0}^{2}}$ around the point $\frac{\mathrm{i}}{2\alpha_{\mathrm{i}}k_{0}^{2}}$ in the complex plane.

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, $|\boldsymbol{k}|^{2}$, for normalizing the k-vectors must be stored.

S1C   Computation and memory efficiency

The computational efficiency of the modified Born series iteration has been discussed previously by Osnabrugge et al. 2. 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\times 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 2, and $3$ times slower than the isotropic vector algorithm 3. It should also be noted that simulating inhomogeneous magnetic properties introduces the discretized differential operator, $\boldsymbol{\mathcal{D}}$. 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 2.

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, $\boldsymbol{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, $\Delta\boldsymbol{E}$, thus bringing the total to $6P$.

Anisotropic permittivity can be represented by a $3\times 3$ matrix for each point in space, to make up a block-diagonal matrix of dimension $3P\times 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 (S1B), though it is necessary to determine the normalization factor $\|\boldsymbol{k}\|$ for all k-vectors. This requires storage for $P$ values.

Anisotropy in the permeability, $\boldsymbol{\mu}$, or the coupling factors, $\boldsymbol{\xi}$ and $\boldsymbol{\zeta}$, will increase the memory requirements accordingly. If the source current $\boldsymbol{S}$ is not sparse, a further $3P$ complex values need to be stored. The memory requirements are summarized in Table S1 for the case of a spatially variant $\boldsymbol{\epsilon}$, and optionally, a spatially variant $\boldsymbol{\mu}$, $\boldsymbol{\xi}$, and $\boldsymbol{\zeta}$. The memory requirements are listed in different columns for the isotropic and the anisotropic cases. Mixed cases are omitted for brevity.

S1D   Sampling and prevention of aliasing

The sampling grid size and density are important considerations for the calculation accuracy. In Algorithm S1, every iteration requires a multiplication by the preconditioned general susceptibility, $\boldsymbol{\chi}$, an addition of the source, $\boldsymbol{S}$, a convolution with the dyadic Green’s function, and another multiplication by $\boldsymbol{\chi}$. Multiplications and additions cannot increase the spatial extent of the field $\boldsymbol{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:

[TABLE]

In Fourier space, each step requires a convolution, $\ast$, with the Fourier transform of the susceptibility, $\widetilde{\boldsymbol{\chi}}$, followed by an addition with the Fourier transform of the source, $\widetilde{\boldsymbol{S}}$, a multiplication with the Fourier transform of the dyadic Green’s function, $\widetilde{\boldsymbol{G}}$, and another convolution with $\widetilde{\boldsymbol{\chi}}$. Although the multiplication by $\widetilde{\boldsymbol{G}}$, tends to suppress high spatial frequencies, it does not strictly limit the bandwidth support of the product. On the other hand, convolutions with $\widetilde{\boldsymbol{\chi}}$ do extend the bandwidth support twice per iteration. If we define $W_{\boldsymbol{\chi}}$, $W_{\boldsymbol{E}}$, and $W_{\boldsymbol{S}}\leq W_{\boldsymbol{\chi}}+W_{\boldsymbol{E}}$ as the spatial-frequency band-widths of $\boldsymbol{\chi}$, $\boldsymbol{E}$, and $\boldsymbol{S}$, respectively; it can seen that the bandwidth of the iteration update, $\boldsymbol{\Delta\boldsymbol{E}}$, must be $W_{\boldsymbol{\Delta\boldsymbol{E}}}\leq 2W_{\boldsymbol{\boldsymbol{\chi}}}+W_{\boldsymbol{\boldsymbol{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_{\boldsymbol{\chi}}$ with every iteration step.

In practice, the calculation must be performed on a sample grid with finite bandwidth, $W$, and sample spacing, $W^{-1}$. 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\geq W_{\boldsymbol{\chi}}+W_{\boldsymbol{E}}$, in other words, with a sampling density that is no smaller than the sum of the Nyquist rates for $\boldsymbol{E}$ and $\boldsymbol{\chi}$. The two convolutions with $\widetilde{\boldsymbol{\chi}}$ expand the support in frequency space to $2W_{\boldsymbol{\chi}}+W_{\boldsymbol{E}}$, thereby causing aliasing in the upper $W_{\boldsymbol{\chi}}$-part of the band. However, the lowest spatial frequencies in the $W_{\boldsymbol{E}}$-band are not affected. All aliasing artefacts can thus be avoided by eliminating all but the spatial frequencies within the lower $W_{\boldsymbol{E}}$-band of the iteration update, $\boldsymbol{\Delta\boldsymbol{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 4, 5.

S1E   Robustness to errors

Unlike time-stepped methods such as FDTD, the iteration presented here is robust to numerical errors. The $N^{\mathrm{th}}$-correction term, $\Delta\boldsymbol{E}_{N}$, is obtained from recursive equation (S1):

[TABLE]

For the initial field estimate, If we begin with a null field estimate $\boldsymbol{E}_{0}=\boldsymbol{0}$, this iteration is equivalent to the modified series (S1). However, it should be noted that when $\boldsymbol{E}_{0}\neq\boldsymbol{0}$, the iteration corresponds to a different series, yet one that converges under the same conditions and to the same limit:

[TABLE]

It can be seen that the $N^{\mathrm{th}}$-term of the series differs by $\boldsymbol{M}^{N}\boldsymbol{E}_{0}$ from that of series (S1), and that its corresponding residue, $\left\|\boldsymbol{E}-\boldsymbol{E}_{N}\right\|$, is

[TABLE]

Provided that the absolute largest eigenvalue of $\boldsymbol{M}$ is smaller than one, the upper bound on the residue is tightened with every iteration, independently of the choice of the initial field. However, it can be seen from equation (S19) that starting from an approximate solution has a similar effect as skipping the first iterations, and allowing the algorithm to reach convergence in less time. The corrections in consecutive terms also prevent the accumulation of numerical errors, a common issue with techniques such as the finite-difference time-domain method.

Since the iteration does not have to be started from the all-zero field, convergence may be reached faster if an approximate solution is provided as the initial field of the iteration. It is straightforward to start the algorithm with an approximate solution to the corresponding isotropic or scalar problem. Techniques such as the beam-propagation method or two-dimensional fast marching 6, could be leveraged as an initialization of the anisotropic algorithm so to reduce the total run-time.

S2   Convergence of the modified Born series

S2A   Geometrical Interpretation

The iteration as given by equation (S1) calculates a series of correction terms of the form $\boldsymbol{M}^{p}\boldsymbol{E}_{0}$. Independently of the initial value $\boldsymbol{E}_{0}$, this series is guaranteed to converge when all eigenvalues of $\boldsymbol{M}$ are less than one in absolute value. Recall that the preconditioned iteration operator, $\boldsymbol{M}$, is given by definition (S2).

In the previous section it is shown with equation (S14) that the dyadic Green function, \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}, can be decomposed in the scaled sum of an identity and a unitary transformation, $\boldsymbol{U}$:

[TABLE]

Substitution in equation (S2) allows it to written in terms of just the preconditioner, $\boldsymbol{\Gamma}$, and the susceptibility, $\boldsymbol{\chi}$:

[TABLE]

The operator $\boldsymbol{\Gamma}$ must be chosen so that the largest eigenvalue of the preconditioned iteration operator, $\boldsymbol{M}$, is less than $1$ in absolute value. In other words, $\left|\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{\max}\right|<1$, where $\boldsymbol{E}_{\max}$ is the eigenvector of the largest eigenfunction of $\boldsymbol{M}$. We will now investigate the effects of these operations using a geometrical construction of the complex values that $\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{\max}$ can take and show under what conditions these fall within the unit disk around the origin in the complex plane.

The expression $\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{\max}$ is constructed step-by-step, starting from the permittivity, $\boldsymbol{\epsilon}\equiv\boldsymbol{\chi}+\alpha\mathbbold{1}_{3}$. The real and imaginary parts of $\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{\epsilon}\boldsymbol{E}_{\max}$ can be seen to be:

[TABLE]

where gain-free media have a positive definite imaginary part. Assuming that the medium is free from gain and finite-valued, then there must exist a semi-circle as shown in blue in Figure S1(a) that contains all values $\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{\epsilon}\boldsymbol{E}_{\max}$ in the upper half of the complex plane. The values $\boldsymbol{E}_{\max}^{\dagger}\left(\boldsymbol{\epsilon}-\alpha_{\mathrm{r}}\mathbbold{1}_{3}\right)\boldsymbol{E}_{\max}$ are contained within the same semi-circle, now centered at the origin as shown in Figure S1(b). Although $\boldsymbol{E}_{\max}$ is not necessarily also an eigenfunction of $\boldsymbol{\chi}$, it can be seen that when $\boldsymbol{\epsilon}$ is positive definite, such semi-circle must exist with radius, $\alpha_{\mathrm{i}}$, equal to the largest singular value of $\boldsymbol{\epsilon}-\alpha_{\mathrm{r}}\mathbbold{1}_{3}$. The values $\boldsymbol{E}_{\max}^{\dagger}\left(\boldsymbol{\epsilon}-\alpha\mathbbold{1}_{3}\right)\boldsymbol{E}_{\max}=\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{\chi}\boldsymbol{E}_{\max}$ can be seen to be contained within the semi-circle centered at $-\mathrm{i}\alpha_{\mathrm{i}}$ and with its circumference touching the origin, as shown in Figure S1(c). The multiplication by $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}$ normalizes and rotates the semi-circle $90^{\circ}$ around the origin for the values $\boldsymbol{E}_{\max}^{\dagger}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{\max}$ as shown in Figure S1(d). The radius of the semi-circle is now $1$, and its base diagonal parallel to the imaginary axis through the point $1+0\mathrm{i}$ on the real axis.

When the dyadic Green function tensor is applied to a point in Hilbert space, the result is contained within the sphere with as diagonal the line-segment between the origin and the original point, scaled by $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}k_{0}^{2}}$. The scaling factor has already been treated together with the susceptibility. To understand the effect of the dyadic Green function, we decompose $\boldsymbol{E}_{\max}$ in the eigenfunction basis of $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$ as $\boldsymbol{E}_{\max}=\sum_{i}c_{i}\boldsymbol{E}_{i}$, and first consider its effect on the eigenfunctions $\boldsymbol{E}_{i}$ of $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$ with eigenvalues $\lambda_{i}$. Figure S1(e) shows how the circle’s center point is constructed by the identity term in the dyadic Green function $\boldsymbol{E}_{i}^{\dagger}\frac{\mathbbold{1}_{3}}{2}\left(\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}\right)$, while the inclusion of the unitary term, $\boldsymbol{E}_{i}^{\dagger}\frac{\mathbbold{1}_{3}-\boldsymbol{U}}{2}\left(\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}\right)$, outlines the circle’s circumference as shown in Figure S1(f). Note that the unitary transformation may reduce the magnitude of the final projection, hence the circle only indicates the bounds of a disk of possible values.

The combined effects of the susceptibility and the dyadic Green function allow us to study the bounds on the eigenvalues of the preconditioned iteration operation $\boldsymbol{M}$. The value of $\boldsymbol{E}_{i}^{\dagger}\frac{\mathbbold{1}_{3}-\boldsymbol{U}}{2}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$ can be inside any disk between the origin and any point within the previously described semi-circle. Figure S2(a) shows the semi-circle in blue and the radii and circumference of a disk for one specific example in black dotted line. It can be noted that some values of $\boldsymbol{E}_{i}^{\dagger}\frac{\mathbbold{1}_{3}-\boldsymbol{U}}{2}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$ still fall outwith the unit circle in the complex plane, marked as a green background. This lies at the basis of the divergence and emphasizes the need for the preconditioner. The final three operations to construct $\boldsymbol{E}_{i}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{i}$ as defined in equation (S21) are a subtraction of the identity, a multiplication by the preconditioner $\boldsymbol{\Gamma}$, and the re-addition of the identity. Since $\boldsymbol{\Gamma}$ is a linear operation, it can be understood as a scaling and rotation around the origin in the complex plane, and the three combined operations are the same scaling and rotation, but around the point $1+0\mathrm{i}$. We will now show that this ensures that all values fall within the unit disk. The black dotted circle shown in Figure S2(a) goes through the point $1+0\mathrm{i}$ when $\boldsymbol{E}_{i}^{\dagger}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$ coincides with the base of the blue semi-circle because the base and the real axis form a right-angle triangle with the diagonal of the black dotted circle. To ensure that all points of the black dotted circle shown in Figure S2(a) are on the unit disk, one has to rotate it around the point $1+0\mathrm{i}$ and make it tangent to the unit disk. It can be seen that this can only be achieved when the operation $\boldsymbol{\Gamma}$ introduces a phase $\theta$ that is equal to that introduced by $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$ in the previous step. In other words, $\theta$ is the complex argument of the previously studied expression $\boldsymbol{E}_{i}^{\dagger}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$. This operation rotates the origin of the black dotted circle onto the real axis. The multiplication by $\boldsymbol{\Gamma}$ can also scale the result. The scaling factor should be positive and sufficiently small so that the resulting radius remains below one. When the same scaling factor as $\boldsymbol{E}_{i}^{\dagger}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$ is used, the equilateral triangle is scaled to touch the imaginary axis as shown in Figure S2(c). The base of equilateral triangle (shown in red), can be seen to trace out the values $\boldsymbol{E}_{i}^{\dagger}\left(\mathbbold{1}_{3}-\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\right)\boldsymbol{E}_{i}$ inside the unit disk. Although a larger scaling factor would be permissible in the particular case shown in Figure S1(c), it is convenient to define the preconditioner as $\boldsymbol{\Gamma}\equiv\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$, to ensure that all the eigenvalues are contained within the unit disk.

Similarly one can see that the values $\boldsymbol{E}_{i}^{\dagger}\frac{\mathbbold{1}_{3}-\boldsymbol{U}}{2}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$ on the circumference of the blue semi-circle are mapped to fall into the unit disk when using the same rotation angle, $\theta$, and scaling factor as $\boldsymbol{E}_{i}^{\dagger}\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{i}$. Figure S2(d) shows how a different equilateral triangle can be constructed between the diagonal of the black dotted circle and the real axis. When the rotation caused by $\boldsymbol{\Gamma}$ is equal to that in the previous step, it rotates the apex of the equilateral triangle onto the real axis as shown in Figure S2(e). The black dotted circle now falls within the bounds of the unit disk. When the preconditioner $\boldsymbol{\Gamma}$ scales the result as in the previous step, the apex translates to coincide with the origin and the diagonal of the black dotted circle goes through the origin. The points connecting the base of the equilateral triangle must therefore lie on the same circle. Since one of the endpoints lies on the circumference of the unit disk, the diagonal line of the black dotted circle must fall entirely within the unit disk.

Values that are in between the two extreme cases discussed above can be written as a weighted sum. The linearity of the operation and the convexity of the unit sphere guarantees that also non-extreme values have eigenvalues that are no greater than those encountered at the boundaries. While this shows that the $\left|\boldsymbol{E}_{i}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{i}\right|\leq 1$ for all eigenfunctions, $\boldsymbol{E}_{i}$, of $\boldsymbol{\chi}$; the same is not necessarily true for, $\boldsymbol{E}_{\max}$, an eigenfunction of $\boldsymbol{M}$ that may be a linear combinations of the eigenfunctions $\boldsymbol{E}_{i}$ of $\boldsymbol{\chi}$.

The general expression $\boldsymbol{E}_{\max}^{\dagger}\boldsymbol{M}\boldsymbol{E}_{\max}$ can be written as the difference

[TABLE]

In the special case that the eigenfunctions of $\boldsymbol{\chi}$ are orthonormal and using \overset{\text{\tiny\leftrightarrow}}{\boldsymbol{G}}\equiv\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}, the first term can be rewritten as a weighted sum over the complex values corresponding to each eigenfunction of the system using $\boldsymbol{E}_{\max}=\sum_{i}c_{i}\boldsymbol{E}_{i}$ with $c_{i}=\boldsymbol{E}_{i}^{\dagger}\boldsymbol{E}_{\max}$ and $\sum_{i}|c_{i}|^{2}=1$:

[TABLE]

Here, the summation term defines the center position of the circle for $\boldsymbol{E}_{\max}$ as a weighted average of the center positions of the circles corresponding to the eigenfunctions $\boldsymbol{E}_{i}$. Since for each eigenfunction, a circle with radius $\frac{1}{2}|\lambda_{i}|^{2}$ fits within the unit disk of the complex plane so the absolute value of the summation term is limited by the inequality

[TABLE]

This leaves sufficient space to fit a circle of radius $\frac{1}{2}\sum_{i}|c_{i}|^{2}|\lambda_{i}|^{2}$ within the unit disk of the complex plane. The second term in equation (S24) is one-half of the dot product of a unitary operation and two terms with identical $l_{2}$-norm $\left\|\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}\boldsymbol{E}_{\max}\right\|=\sqrt{\sum_{i}|c_{i}|^{2}|\lambda_{i}|^{2}}=\left\|\boldsymbol{\Gamma}^{\dagger}\boldsymbol{E}_{\max}\right\|$, when $\boldsymbol{\Gamma}$ is defined as $\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$. Hence, by defining $\boldsymbol{\Gamma}\equiv\frac{\mathrm{i}}{\alpha_{\mathrm{i}}}\boldsymbol{\chi}$, it can be seen that equation (S24) cannot be larger than $1$, where $\alpha_{\mathrm{i}}$ is larger than the largest singular value of $\|\boldsymbol{\Delta}\|$. The series must thus converge when the eigenfunctions of the susceptibility distribution, $\boldsymbol{\chi}$, are orthogonal, a very common situation. Yet, its eigenfunctions will not be orthogonal when the reactive, $\boldsymbol{\Delta}_{\mathrm{r}}$, and dissipative parts, $\boldsymbol{\Delta}_{\mathrm{i}}$, of $\boldsymbol{\Delta}=\boldsymbol{\chi}+{\mathrm{i}}\alpha_{\mathrm{i}}\mathbbold{1}_{3}=\boldsymbol{\Delta}_{\mathrm{r}}+\mathrm{i}\boldsymbol{\Delta}_{\mathrm{i}}$ 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 analyzed.

S2B   Numerical demonstration of the convergence

In the main text we showed that the following choice of $\alpha_{\mathrm{i}}$ guarantees convergence of the modified Born series

[TABLE]

where $A_{1}$ and $A_{2}$ are defined as:

[TABLE]

It the main text it was noted that without constraining $\boldsymbol{\Delta}_{\mathrm{i}}$ to be positive definite, $A_{2}$ can be arbitrarily large, when $\langle\boldsymbol{n}|\boldsymbol{\Delta}_{\mathrm{i}}|\boldsymbol{n}\rangle$ is close to zero. We begin this Supplementary Section by showing that this divergence can be avoided so long as $\boldsymbol{\Delta}_{\mathrm{i}}$ has an empty kernel (no eigenvectors with zero eigenvalue). Suppose that $\boldsymbol{\Delta}_{\mathrm{i}}$ has an eigenvector $|\boldsymbol{n}_{0}\rangle$ with eigenvalue $\lambda$, and consider $|\boldsymbol{n}\rangle=|\boldsymbol{n}_{0}\rangle+\eta|\boldsymbol{n}_{\perp}\rangle$ where $\eta\ll 1$. The expression for $A_{2}$ inside the maximization is then given by

[TABLE]

If $\lambda=0$ (and the kernel of $\boldsymbol{\Delta}_{\mathrm{i}}$ thus contains $|\boldsymbol{n}_{0}\rangle$ then the above quantity diverges as $1/\eta$ as we take $\eta$ to zero. The quantity $A_{2}$ thus becomes infinite. However, if $\lambda$ is non-zero (but arbitrarily small) then as $\eta\to 0$, (S28) tends to $|\langle\boldsymbol{n}_{0}|\boldsymbol{\Delta}|\boldsymbol{n}_{0}\rangle|$ which is finite, and we note, smaller than or equal to the largest singular value of $\boldsymbol{\Delta}$.

Assuming that $\langle\boldsymbol{n}|\boldsymbol{\Delta}_{\mathrm{i}}|\boldsymbol{n}\rangle$ is never zero, the quantity $A_{2}$ can be rewritten in a form that is easier to compute. First we write the positive definite Hermitian matrix $\boldsymbol{\Delta}_{\mathrm{i}}$ as the square of another Hermitian matrix $\boldsymbol{a}$

[TABLE]

Defining $|\boldsymbol{n}^{\prime}\rangle=\boldsymbol{a}|\boldsymbol{n}\rangle$, the quantity $A_{2}$ can be written as

[TABLE]

which is simply the numerical radius $r(\boldsymbol{m})$ of the matrix

[TABLE]

The numerical radius of a matrix is always less than or equal to its norm 8 $r(\boldsymbol{m})\leq\|\boldsymbol{m}\|$. In addition, the norm of a sum is never larger than the sum of the norms, and we can therefore estimate $A_{2}$ as

[TABLE]

where $\|\cdot\|$ indicates the $l_{2}$-norm. While the spectra of $\boldsymbol{\Delta}$ and e.g. $\boldsymbol{a}\boldsymbol{\Delta}\boldsymbol{a}^{-1}$ 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_{2}$ is the magnitude of the largest eigenvalue $|\lambda_{\rm max}|$ of $\boldsymbol{\Delta}$. Meanwhile, $A_{1}$, being the square root of the numerical radius of a Hermitian operator, is equal to the square root of the operator norm $(1/2)\|\boldsymbol{\Delta}\boldsymbol{\Delta}^{\dagger}+\boldsymbol{\Delta}^{\dagger}\boldsymbol{\Delta}\|$, which is bounded by

[TABLE]

where the estimate for $A_{2}$ is here given by equation (S32). We can therefore use this estimate of $A_{2}$ to find an upper bound for the value of $\alpha_{\mathrm{i}}$. Figure S3b shows a numerical test, where we chose $A_{2}$ according to (S32) and repeatedly evaluated the inner product $\langle\boldsymbol{n}|\boldsymbol{M}|\boldsymbol{n}\rangle$ for different choices of random complex $40\times 40$ matrices $\boldsymbol{\Delta}$ and $\boldsymbol{U}$. These tests represent $10^{6}\times 40=4\times 10^{7}$ evaluations of the inner product and appear to indicate that $\alpha_{\mathrm{i}}>\left\|\boldsymbol{\Delta}\right\|$ is a tighter bound (Fig. S3a).

S3   Impedance matching and Bianisotropy

By eliminating back reflections, impedance matching ensures efficient energy transfer and communication through waveguides. Supplementary Figure S4(a-d) introduces the concept of impedance matching in one dimension. Supplementary Figures S4(a,b) shows two samples with isotropic permittivity, $\epsilon$, and permeability, $\mu$. Both samples contain objects (at $10\,\mathrm{\upmu m}<z<20\,\mathrm{\upmu 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$. Supplementary Figure S4(b) thus represents an object that is impedance matched with the surrounding medium. As can been seen from Supplementary Figure S4(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. Supplementary Figure S4(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.

Impedance matching also has important practical applications for simulations of infinite volumes in a finite space. Supplementary Figure S4(e) shows the electric dipole field in a box with regular absorbing boundaries and homogeneous permeability. As in the one-dimensional case, plotted in Supplementary Figure S4(a,c), significant reflections can be noted from the boundaries. Supplementary Figure S4(f) shows a dipole in a box with impedance matched absorbing boundaries. It can be seen that the interfering reflections are suppressed. In higher dimensions, impedance matching is only an approximation to perfectly matching layers. This may explain the weak, low-frequency, beating that can be noticed in the top-right corner.

The ability to account for the coupling factors enables the calculation of electromagnetic waves in materials with chiral properties. Supplementary Figures S4(g,h) show how a chiral substance slowly rotates the linear polarization of a wave that propagates through it. Supplementary Figure S4(g) shows the intensity in the horizontal and vertical polarization component, while Supplementary Figure S4(h) shows how the polarization angle changes linearly over the simulation volume’s width of $10\;\mathrm{mm}$.

Supplementary References

• 1

B. Hecht L. Novotny.

Principles of Nano-Optics.

Cambridge University Press, June 2012.

• 2

Gerwin Osnabrugge, Saroch Leedumrongwatthanakun, and Ivo M Vellekoop.

A convergent born series for solving the inhomogeneous helmholtz equation in arbitrarily large media.

Journal of computational physics, 322:113–124, 2016.

• 3

Benjamin Krüger, Thomas Brenner, and Alwin Kienle.

Solution of the inhomogeneous maxwell’s equations using a born series.

Opt. Express, 25(21):25165–25182, October 2017.

• 4

Alok Dutt and Vladimir Rokhlin.

Fast fourier transforms for nonequispaced data.

SIAM Journal on Scientific computing, 14(6):1368–1393, 1993.

• 5

Jeffrey A Fessler and Bradley P Sutton.

Nonuniform fast fourier transforms using min-max interpolation.

IEEE Transactions on Signal Processing, 51(2):560–574, 2003.

• 6

Amedeo Capozzoli, Claudio Curcio, Angelo Liseno, and Salvatore Savarese.

Two-dimensional fast marching for geometrical optics.

Optics express, 22(22):26680–26695, 2014.

• 7

T. E. Oliphant.

A guide to NumPy.

Trelgol Publishing, 2006.

• 8

S. S. Dragomir.

Inequalities for the Numerical Radius of Linear Operators in Hilbert Spaces.

Springer, 2013.