Calculating coherent light-wave propagation in large heterogeneous media

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.


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][5][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 2 , the lipid bilayer cell-membrane [8][9][10] , or muscle fibers 11 . Indeed, birefringence is a valuable, label-free, contrast method [11][12][13][14][15][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: D(x) = 0 (x)E(x), and B(x) = µ 0 H(x) (1) where 0 and µ 0 are the vacuum permittivity and permeability, respectively. By substituting both constitutive relations into Maxwell's equations we obtain the vector Helmholtz equation for the electric field where k 0 = ω/c is the free space wavenumber and S(x) represents the radiation source. In order to solve equation (2) for the electric field E(x), one must invert the operator O = −k 2 0 (x) + ∇ × ∇×. We could attempt to do this directly, representing O as a matrix in some basis of eigenfunctions and then applying a numerical inversion algorithm. However, for moderately large systems the direct representation and inversion of such matrix is infeasible, due to demands of both memory and time. Iterative solutions are thus required where, starting from some initial guess for the electric field, one repeatedly applies a matrix until the result is arbitrarily close to the one that would be obtained from a direct inversion, as for instance in finite element simulations (see e.g. 19 , Chapter 19).
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 where 1 3 is the identity operator in three-dimensional space. The operator O h is typically associated with propagation through a homogeneous medium, and O i is due to the inhomogeneity of the material. We can thus understand the Born series (3) as representing a sequence of scattering events, where as our initial guess the source S generates the same field as if the material were homogeneous E 0 (x) = O −1 h S, the first iteration adds to this events where a single scattering event , and so on. The main problem with this intuitive expansion is that in many cases it will not converge. Physically, this is connected to the existence of bound states 22 , which are due to the constructive interference of an infinite number of scattering events. On the other hand, the mathematical origin of this divergence is simple, it is the same as the divergence of the scalar series (1 − q) −1 = ∑ ∞ p=0 q p , ∀q : q ∈ C: for convergence we must require that |q| < 1. Correspondingly, the Born series, as given by (3), will only converge if all of the eigenvalues of O −1 h O i have a magnitude less than unity (for discussion see 22 ). Although there are cases where exact results can be deduced from this series 23 , the application of the Born series is generally restricted to weakly scattering media.
For our particular case we take the homogeneous medium of O h to have a constant permittivity α = α r + iα i . The inverse of this operator is given by the Green function for this medium, The remaining part of our operator O (2), is the spatially- is the isotropic susceptibility with respect to our choice of background permittivity, α ∈ C. The Born series, as defined in (3), for the electric field E then becomes Although we have not yet addressed the problem of convergence, the numerical advantage of this iterative series is that it can be done in a constant space and with a combination of fast-Fourier Transform algorithm 24 ) and operations on matrices that are diagonal in their spatial indices. To see this, we observe that χ is diagonal in its spatial indices, and that the dyadic Green function (the solution to equation (4)) can be written as a product of a Fourier transform F , a diagonal matrix (in Fourier space) and an inverse Fourier transform (see supplementary section S1B) where Π T and Π L are projection matrices that separate the field in a super-position of plane waves with transverse and longitudinal electric-components, respectively. The projection operations are defined as the outer product of the normalized kvectors, Π L = k ⊗ k/ k 2 , and Π T = 1 3 − Π L . The quantity in the rounded brackets can thus be seen to be a diagonal matrix in the two Fourier space indices. As mentioned above, in a pair of recent papers 17,18 , it was shown that this Born series (5) can be made into a convergent series, the only proviso being that the system exhibits solely dissipation (no gain), and that the permittivity is nowhere singular. To achieve this convergence the authors used two ingenious steps. The first is to notice that the complex background permittivity α can be chosen freely, changing the eigenvalues of the operator k 2 0 ↔ Gχ, without affecting the solution. The second is that one can apply a pre-conditioner, Γ, to both sides of the recursion relation (5) as follows: The corresponding so-called 'modified Born series', which is the counterpart of series (5), is then given by which converges when the absolute value of the largest eigenvalue of M is less than 1. This condition has been shown 17,18 to hold when both the preconditioning operator is given as Γ = i αi χ and when the imaginary part of the background permittivity, α i , is larger than the largest value of |∆(x)| = | (x) − α r |, considered over all points, x, in space. The aforementioned work showed that this modified series can be used to very quickly compute the electromagnetic field in large media of moderate index contrast. This result relies on the isotropy of the permittivity, (x) and on the absence of any magnetic effect or chirality. In this paper we generalize this approach to media with any permittivity or permeability, including anisotropic, magnetic, chiral, and bi-anisotropic media. Interestingly we find that proving the convergence of the modified series seems to be much more subtle when these additional optical effects are accounted for. We also discuss the convergence speed and memory requirements of this approach, and present an easy to implement algorithm to simulate electromagnetic wave propagation in an arbitrary medium. The method is implemented as a Python 3 library and available with examples for download at https://uod.box.com/v/macromax. 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 25 , 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.

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 and similarly the susceptibility becomes χ(x) = (x) − α1 3 . Applying this constitutive relation to the Maxwell equations as we did to obtain (2), the vector Helmholtz equation becomes An anisotropic analogue of the modified Born series is given by replacing the scalar preconditioner with a matrix function, Γ(x) ≡ i αi χ(x), in equation (8). The iteration matrix can now be written as a function of the matrix function χ(x): 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 (8) through the requirement max n n ik 2 In Supplementary Section S2 it is shown that this requirement can be satisfied by choosing α i to be larger than ∆ , the largest singular value of ∆ ≡ − α r 1 3 , under the conditions that the eigenvectors of ∆ are orthogonal and that the material is free of gain with non-zero losses. Orthogonal eigenvectors are common in practice, e.g. when there is no point with both anisotropic absorption and birefringence in the material. In what follows we further show that a sufficiently large value of α i can ensure convergence for any material with non-zero losses.
The dyadic Green function, ↔ G, in convergence condition (12) can be written in terms of the identity and an unitary opera- Section S1B). Condition (12) can then be rewritten using the triangle inequality |a + b| ≤ |a| + |b| as max n n|∆ 2 − α 2 i 1 3 |n + | n |χUχ| n | < 2α 2 i (13) where ∆ = χ + iα i 1 3 . To simplify this condition further we apply the Cauchy-Schwarz inequality | a|b | 2 ≤ a|a b|b to remove the unitary operator from the second term, followed by the inequality a 2 + b 2 ≥ 2ab where ∆ = ∆ r + i∆ i , with ∆ r and ∆ i Hermitian matrices which depend respectively on the reactive and dissipative response of the material. As a final step we rewrite and apply the triangle inequality for a second time We are now in a position to give a condition for the convergence of our modified Born series. Assuming that ∆ i is a positive definite operator everywhere, equivalent to assuming a gain-free material that exhibits a dissipative response throughout space, we can restrict α i to be greater than the larger of the following two quantities: If the convergence condition is fulfilled then the inequality (15) will be satisfied, and the modified Born series will converge for the anisotropic medium in question. 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 0 then for vectors |n = |n 0 + η|n ⊥ , A 2 diverges as ∼ 1/η, as η → 0. This is because A 2 is the analogue of a weak value 26 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 27 ). In Supplementary Section S2B 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 (18), because the value of A 2 cannot be readily estimated. Extensive numerical tests (Supp. Sec. S2B) have also indicate the larges singular value ∆ is a tighter bound for α i that is valid for general gainfree 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, ∆ = ∆1 3 , and the basis |n is the position basis, convergence Figure 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 3 ) microrod with a diameter of 20 µm forms a complex scattering pattern instead of a single focus. Although the volume is homogeneous CaCO 3 , 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.
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 3 , 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 28 . Fig. 1a 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.
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 equation (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) 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 ( Fig. 1be). 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 fieldcomponents as the beam exits at the CaCO 3 -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 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 crosspolarizers (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.
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 1 3 , using the Nelder-Mead algorithm 29 . 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 (S1) presented in Supplementary Section S1, 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.
The iteration loop starts on line 6 and repeats until the con- dition 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 .
Algorithm 1 Calculation of the electric field in materials with anisotropic permittivity, .

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 Supplementary Section S1A. The susceptibility, χ, must be redefined for magnetic materials to account for permeability and any electro-magnetic coupling constants. A permeability scale, β, is introduced, and the background permittivity, α, is determined from the modified susceptibility.
The constitutive relations for general linear electromagnetic materials are given by (see e.g. 30 ): Here both the permittivity (x) and permeability, µ(x), are tensors that depend on the position, x. The additional coupling terms ξ(x) and ζ(x) enable arbitrary linear interactions between the electric and magnetic components, which are commonplace in-for instance-metamaterials 31 . The coupling tensors ξ(x) and ζ(x) are also essential to model chiral materials 32 , and non-reciprocal materials such as moving media 33 , and Tellegen media 25 . Substituting (19) into the Maxwell equations and eliminating the magnetic field, H = (iωµ 0 ) −1 µ −1 [∇ × E − ik 0 ζE], leads to the same form of vector Helmholtz equation (10). It is sufficient to generalize the source and susceptibility as follows: where D ≡ k −1 0 ∇× is the differential operator, which calculation is discussed in Supplementary Section S1A. For conciseness, the spatial dependency of the variables is omitted in equations (20) and (21). Note that all these operations can be implemented efficiently using fast-Fourier transforms and diagonalmatrix multiplications. Since vector Helmholtz equation (10) still applies, the same iterative procedure can be used to solve for non-magnetic anisotropic materials.
Although the modified Born series is now formally identical for these general constitutive relations (19), to that for the anisotropic dielectric discussed in section 2, it is important to note that the susceptibility (21) is no longer a diagonal matrix in the position variable. This is because it now contains the differential operators D. Although this affects neither the form of the modified Born series (i.e. as a product of diagonal matrices and Fourier transforms), or 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 D has a magnitude of π/(k 0 ∆r), where ∆r is the grid spacing. Therefore, the required number of iterations increases with the density of the sampling grid when simulating general magnetic materials.
The optimal value of α is determined using a similar method in both the all-dielectric and general cases. Albeit not critical for convergence, the additional parameter β we have introduced in (21) can be used to enable further optimization of the convergence speed.
Algorithm S1 in Supplementary Section S1A generalizes Algorithm 1 to magnetic materials by redefiningσ as a function of two variables, α r and β, and minimizing it in both variables. Using the usual inequalities for matrix norms 34 , the function σ(β, α r ) determines an upper bound for the singular values of the separate terms in the sum (21), and using the largest singular value, σ D = π/(k 0 ∆r), of the discretized differential operator, D. The Nelder-Mead method is used to numerically minimize the productσ(β, α)β, which we found to correlate well with the inverse of the convergence rate. 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 1 3 .
At optical frequencies, it is a common approximation to assume that the permeability, µ, equals µ 0 1 3 in the entire simulation volume. While this is appropriate at such frequencies (where magnetism can in many circumstances be ignored 33,35 ), the modified Born series can-in principle-be applied to any electromagnetic problem, where magnetic effects are important. Indeed, while all our simulations are apparently for optical length scales we must remember that Maxwell's equations are scale invariant so that the same results would be obtained with the spatial frequencies, k, scaled by the same factor as the temporal frequency, ω. In Supplementary Figure S4(a-f), we show that our generalized Born series can simulate impedance matched media (where = µ, which includes transformation media 36 , as well as other interesting wave effects (e.g. 37 ).
Since the constitutive relations (19) include the electricmagnetic coupling tensors, ζ and ξ, one can also use the modified Born series to treat bi-isotropic and more generally bianisotropic 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, Supplementary Figures S4(g and h) show simulations where linearly polarized light is slowly rotated in a chiral solution. As the rotation takes several millimeters to complete, the calculation was performed for a wavelength of 500 nm 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 38,39 . The general algorithm presented here naturally handles negative values for both the permittivity, , and the permeability, µ.

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. We have implemented this numerically, and shown that a large variety of electromagnetic effects can be simulated using this method. The complete is given in Algorithm S1, and implemented as a numerical library that can be found with working examples at https://uod.box.com/v/macromax. 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. Moreover, this opens the door to precalculate approximate solutions using, for instance, the beam propagation method (Supp. Section S1E). 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 (Supp. Section S1C). 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 2 ). The presented numerical method offers a bridge between approximate analytic predictions and experiments. In addition, such largescale 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.

Funding Information
This work is supported by the Philip Leverhulme Prize from the Leverhulme Trust. SARH acknowledges financial support from the Royal Society and TATA (RPG-2016-186).

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.

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 [15][16][17][18][19][20][21][22][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 2 -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
with ⊗ the outer product, while the forward and inverse Fourier transforms are represented by F and 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 × 3dot products for each application of the operator χ on the electric field, E. The source is defined by S ≡ iωµ 0 j/ βk 2 0 . The dyadic Green's function, ↔ G * , is discussed in section S1B in what follows. Note that to prevent issues with numerical precision, in the implementation the factor k −2 0 is moved from the definition of ↔ G to that of S.

S1B The dyadic Green's function
The dyadic Green function, ↔ G, is integral to the calculation of the modified Born series: 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, U.
The vector Helmholtz equation can be separated into a homogeneous, O h , and an inhomogeneous part, O i , as follows The dyadic Green function, ↔ G, is defined as the impulse response solution to the homogeneous part: First, we take the Fourier transform of equation (S5), 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 or equivalently in the operator notation used in the main text (where integrals are subsumed into our '' product) Algorithm S1 A function that implements the general algorithm for both non-magnetic and magnetic materials.

4:
β ← 1/µ 11 5: where σ D = k −1 0 π/ ∆r is proportional to the highest possible spatial frequency 8: β, α r ← argmin β,αr [σ (β, α r ) |β|] ∀β ∈ R >0 , α r ∈ R 9: end if 10: α i ←σ(β, α r ) 11: α ← α r + iα i 12: p ← ∞ the l 2 -norm of the previous update 15: repeat 16: At this point we decompose the identity matrix on the right into two parts 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 asG(k) = g L (k)Π L + g T (k)Π T , finding that equation (S8) separates into two parts from which we can deduce that the Green function is given by which in our operator notation we write as 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 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 α being chosen as a complex number. To find this representation we first note the following identity for any complex number z where z indicates the complex argument of z. Applying identity (S13) to (S12), the Green function operator (S12) is reduced to the following form where α i is the imaginary part of α and the unitary matrix U is given by where k = k is the l 2 -norm of k, while α * is the complex conjugate of α. The operator (S15) is unitary, due to the fact that our Fourier transform operators can be chosen such that F † = F −1 , and thus where we used the fact that Π T Π L = 0, Π T Π T = Π T , and

S2
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 eigen- We note that for numerical purposes it is space and timeefficient 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| 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 × 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, 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, 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 (S1B), 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 S1 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.

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, χ, 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: 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 −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 bandlimited, 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 4,5 .

S1E Robustness to errors
Unlike time-stepped methods such as FDTD, the iteration presented here is robust to numerical errors. The N thcorrection term, ∆E N , is obtained from recursive equation (S1): For the initial field estimate, If we begin with a null field estimate E 0 = 0, this iteration is equivalent to the modified series (S1). However, it should be noted that when E 0 = 0, the iteration corresponds to a different series, yet one that converges under the same conditions and to the same limit: It can be seen that the N th -term of the series differs by M N E 0 from that of series (S1), and that its corresponding residue, 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. 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 allzero 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 twodimensional fast marching 6 , could be leveraged as an initialization of the anisotropic algorithm so to reduce the total runtime.

S2A Geometrical Interpretation
The iteration as given by equation (S1) calculates a series of correction terms of the form M p E 0 . Independently of the initial value E 0 , this series is guaranteed to converge when all eigenvalues of M are less than one in absolute value. Recall that the preconditioned iteration operator, M, is given by definition (S2).
In the previous section it is shown with equation (S14) that the dyadic Green function, ↔ G, can be decomposed in the scaled sum of an identity and a unitary transformation, U: Substitution in equation (S2) allows it to written in terms of just the preconditioner, Γ, and the susceptibility, χ: The operator Γ must be chosen so that the largest eigenvalue of the preconditioned iteration operator, M, is less than 1 in absolute value. In other words, E † max ME max < 1, where E max is the eigenvector of the largest eigenfunction of M. We will now investigate the effects of these operations using a geometrical construction of the complex values that E † max ME max can take and show under what conditions these fall within the unit disk around the origin in the complex plane.  The expression E † max ME max is constructed step-by-step, starting from the permittivity, ≡ χ + α1 3 . The real and imaginary parts of E † max E max can be seen to be: 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 E † max E max in the upper half S4 of the complex plane. The values E † max ( − α r 1 3 ) E max are contained within the same semi-circle, now centered at the origin as shown in Figure S1(b). Although E max is not necessarily also an eigenfunction of χ, it can be seen that when is positive definite, such semi-circle must exist with radius, α i , equal to the largest singular value of − α r 1 3 . The values E † max ( − α1 3 ) E max = E † max χE max can be seen to be contained within the semi-circle centered at −iα i and with its circumference touching the origin, as shown in Figure S1(c). The multiplication by i αi normalizes and rotates the semi-circle 90 • around the origin for the values E † max i αi χ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 + 0i 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 i αik 2 0 . The scaling factor has already been treated together with the susceptibility. To understand the effect of the dyadic Green function, we decompose E max in the eigenfunction basis of i αi χ as E max = ∑ i c i E i , and first consider its effect on the eigenfunctions E i of i αi χ with eigenvalues λ i . Figure S1 αi χ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 E † i ME i as defined in equation (S21) are a subtraction of the identity, a multiplication by the preconditioner Γ, and the re-addition of the identity. Since Γ 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 + 0i. 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 + 0i when E † i i αi χ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 + 0i and make it tangent to the unit disk. It can be seen that this can only be achieved when the operation Γ introduces a phase θ that is equal to that introduced by i αi χ in the previous step. In other words, θ is the complex argument of the previously studied expression E † i i αi χE i . This operation rotates the origin of the black dotted circle onto the real axis. The multiplication by Γ 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 E † i i αi χ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 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 Γ ≡ i αi χ, to ensure that all the eigenvalues are contained within the unit disk.
Similarly one can see that the values E † i 13−U 2 i αi χE i on the circumference of the blue semi-circle are mapped to fall into the unit disk when using the same rotation angle, θ, and scaling factor as E † i i αi χ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 Γ 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 Γ 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 E † i ME i ≤ 1 for all eigenfunctions, E i , of χ; the same is not necessarily true for, E max , an eigenfunction of M that may be a linear combinations of the eigenfunctions E i of χ.
The general expression E † max ME max can be written as the dif-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 i αi χE max = ∑ i |c i | 2 |λ i | 2 = Γ † E max , when Γ is defined as i αi χ. Hence, by defining Γ ≡ i αi χ, it can be seen that equation (S24) cannot be larger than 1, where α i is larger than the largest singular value of ∆ . The series must thus converge when the eigenfunctions of the susceptibility distribution, χ, are orthogonal, a very common situation. Yet, its eigenfunctions will not be orthogonal when the reactive, ∆ r , and dissipative parts, ∆ i , of ∆ = χ + iα i 1 3 = ∆ r + i∆ 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 α i guarantees convergence of the modified Born series where A 1 and A 2 are defined as: It the main text it was noted that without constraining ∆ i to be positive definite, A 2 can be arbitrarily large, when n|∆ i |n is close to zero. We begin this Supplementary 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 a b Figure S3: Numerical check that condition α i > max {A 1 , A 2 } ensures that the numerical radius of M (given by equation (S2)) 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 6 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 7 ). 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 6 points plotted in each panel. In (b) (largest magnitude 0.999997) we used the condition (S32), 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 has an eigenvector |n 0 with eigenvalue λ, and consider |n = |n 0 + η|n ⊥ where η 1. The expression for A 2 inside the maximization is then given by | n|∆∆ i + ∆ i ∆|n | 2 n|∆ i |n = 2λ n 0 |∆|n 0 + η ( n ⊥ |∆ i ∆|n 0 + n 0 |∆∆ i |n ⊥ ) + +ηλ( n ⊥ |∆|n 0 + n 0 |∆|n ⊥ ) + η 2 n ⊥ |∆∆ i + ∆ i ∆|n ⊥ | 2 (λ + η 2 n ⊥ |∆ i |n ⊥ ) If λ = 0 (and the kernel of ∆ i thus contains |n 0 then the above quantity diverges as 1/η as we take η to zero. The quantity A 2 thus becomes infinite. However, if λ is non-zero (but arbitrarily small) then as η → 0, (S28) tends to | n 0 |∆|n 0 | 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 2 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 ∆ i = a 2 (S29) Defining |n = a|n , the quantity A 2 can be written as The numerical radius of a matrix is always less than or equal to its norm 8 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 2 as where · indicates the l 2 -norm. While the spectra of ∆ and e.g. a∆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 |λ max | of ∆. 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) ∆∆ † + ∆ † ∆ , which is bounded by 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 α i . Figure S3b shows a numerical test, where we chose A 2 according to (S32) and repeatedly evaluated the inner product n|M|n for different choices of random complex 40 × 40 matrices ∆ and U. These tests represent 10 6 × 40 = 4 × 10 7 evaluations of the inner product and appear to indicate that α i > ∆ 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, , 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. 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, lowfrequency, 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 mm.