Convergence of an Explicit Iterative Leap-frog Discontinuous Galerkin Method for Time-domain Maxwell's Equations in Anisotropic Materials

We propose an explicit iterative leap-frog discontinuous Galerkin method for time-domain Maxwell's equations in anisotropic materials and derive its convergence properties. The a priori error estimates are illustrated by numerical means in some experiments. Motivated by a real application which encompasses modeling electromagnetic wave's propagation through the eye's structures, we simulate our model in a 2D domain aiming to represent a simple example of light scattering in the outer nuclear layer of the retina.


Introduction
The human retina is a complex structure in the eye that is responsible for the sense of vision. It is part of the central nervous system and it is composed by several layers, namely the outer nuclear layer that comprises the cells bodies of light sensitive photoreceptors cells, rods and cones (see Fig. 1) [12]. For many diseases that affect the eye, the diagnosis is not straightforward. The sensitivity of this structure makes medical analysis particularly complicated. Most of the diagnoses are made either by direct observation, with the possible injection of dyes, to enhance certain parts of the organ, or by numbing the eye and directly measuring its inner pressure or thickness. There are a number of eye-related pathologies that can be identified by the detailed analysis of the retinal layers [17].
Optical Coherence Tomography (OCT) is an increasingly popular noninvasive technique that has been successfully used as a diagnostic tool in ophthalmology in the past decades. This method allows the assessment of the human retina in vivo and has been shown to provide functional information. By analysing data acquired through OCT, several retinal pathologies, such as diabetic retinopathy, or macular edema, can be detected in their early stages, before noticeable morphologic alterations on the retina [17]. As OCT standard techniques only provide structural information [15], it is necessary to expand OCT data analysis to account for both structural and functional information. OCT provides also the possibility of evaluating different elements in measuring the retinal nerve fiber layer (RNFL), namely the tendency of RNFL thinning in glaucoma and other diseases that involve optic nerve atrophy. Waveguides with induced anisotropy may worth to be considered for modeling biological waveguides [10]. Maxwell's equations are a fundamental set of partial differential equations which describe electromagnetic wave interactions with materials. The electromagnetic fields in space are classically described by two field vectors, E and H, respectively electric field and magnetic field. Here we shall consider the time domain Maxwell's equations in the transverse electric (TE) mode, as in [8], where the only non-vanishing components of the electromagnetic fields are E x , E y and H z . Using the following notation for the vector and scalar curl operators and assuming no conductivity effects, the equations in the non-dimensional form are where E = (E x , E y ) represents the electric field components and H = (H z ) represents the magnetic field component. These equations are set and solved on a bounded polygonal domain Ω ⊂ R 2 . The electric permittivity of the medium, , and the magnetic permeability of the medium, µ, are varying in space, being µ a scalar function and an anisotropic tensor We assume that electric permittivity tensor is symmetric and uniformly positive definite for almost every (x, y) ∈ Ω, and that it is uniformly bounded with a strictly positive lower bound, i.e., there are constants > 0 and > 0 such that, for almost every (x, y) ∈ Ω, |ξ| 2 ≤ ξ T (x, y)ξ ≤ |ξ| 2 , ∀ξ ∈ R 2 . We also assume that there are constants µ > 0 and µ > 0 such that, for almost every (x, y) ∈ Ω, µ ≤ µ(x, y) ≤ µ.
Equations (1) must be complemented by initial conditions E(x, y, 0) = E 0 (x, y) and H(x, y, 0) = H 0 (x, y), (x, y) ∈ Ω, and by proper boundary conditions. Motivated by our application of interest, here we consider absorbing boundary conditions which mimic an open space by absorbing the incident radiation in the truncated computational domain. The first order Silver-Müller absorbing boundary conditions (SM-ABC) are defined by where n = (n x , n y ) T is the unit outward normal vector to the boundary and c is the speed with which a wave travels along the direction of the unit normal, defined, using the effective permittivity ef f = det( )/(n T n) (see [8]), by c = 1/ √ µ ef f .
The attention to the development of high-order accurate methods for solving timedomain Maxwell's equations in complex geometries brings to the use of discontinuous Galerkin (DG) methods [7]. The one-step explicit time integration methods, like leap-frog schemes, are computationally efficient per update cycle and easy to implement. The leap-frog DG method in anisotropic materials which is discussed in [4] leads to a locally implicit method for the case of SM-ABC. In [1] a fully explicit in time leap-frog DG method is investigated in the same framework. The error estimates derived therein show that the method is only of first order convergent in time when SM-ABC are considered.
In the present work we propose an iterative predictor-corrector method based on the explicit method investigated in [1], resulting a fully explicit method that is second order convergent in time for the SM-ABC case. In the Section 2 we prove that the explicit iterative method converges to a second order in time implicit method and we deduce the a priori error estimates for the fully discrete scheme. In Section 3 we illustrate the theoretical results with some numerical examples and, in Section 4, we apply the numerical method to a computational model that aims to simulate the light scattering through the outer nuclear layer of the retina.
This work was developed in the framework of a more general project that aims to develop a computational model to simulate the electromagnetic wave's propagation through the eye's structures in order to create a virtual optical coherence tomography scan [14].
2 An explicit iterative leap-frog discontinuous Galerkin method 2.1 Numerical scheme Assume that the computational domain Ω is a bounded polygonal set that is partitioned into K triangular elements T k such that Ω = ∪ k T k . For simplicity, we consider that the resulting mesh T h is conforming. The finite element space is then taken to be V N = {v ∈ L 2 (Ω) 3 : v| T k ∈ P N (T k ) 3 }, where P N (T k ) denotes the space of polynomials of degree less than or equal to N on T k . On each element T k , the solution fields E x (x, y, t), E y (x, y, t), H z (x, y, t) are approximated by the piecewise polynomial functions E xk (x, y, t),Ê yk (x, y, t),Ĥ zk (x, y, t). The approximate fields are allowed to be discontinuous across element boundaries. In this way, we introduce the notation for the jumps of the field values across the interfaces of the elements, [Ê] =Ê − −Ê + and [Ĥ] =Ĥ − −Ĥ + , where the superscript " + " denotes the neighboring element and the superscript " − " refers to the local cell. Furthermore we introduce, respectively, the cell-impedances and cell-conductances Z ± = µ ± c ± and Y ± = (Z ± ) −1 . At the outer cell boundaries we set Z + = Z − . The coupling between elements is introduced via the numerical flux as in [1].
To define the fully discrete scheme, we divide the time interval into M subintervals by the points 0 = t 0 < t 1 < · · · < t M = T , where t m = m∆t, ∆t is the time step size and T + ∆t/2 ≤ T f . The unknowns related to the electric field are approximated at integer time-stations t m and are denoted byÊ m k =Ê k (., t m ). The unknowns related to the magnetic field are approximated at half-integer timestations t m+1/2 = (m + 1 2 )∆t and are denoted byĤ m+1/2 k =Ĥ k (., t m+1/2 ). With the above setting , we can now formulate the iterative leap-frog DG method. The process starts with an approximation to the initial data which we denote by The (n + 1) th inner iteration of the iterative scheme, for n = 0, 1, 2, . . ., is: find where (·, ·) T k and (·, ·) ∂T k denote the classical L 2 (T k ) and L 2 (∂T k ) inner-products andÊ [m+1/2,n] andĤ [m+1,n] are the average approximationŝ The parameter α ∈ {0, 1} can be used to control dissipation. Taking α = 0 yields a non dissipative central flux while α = 1 corresponds to the classic upwind flux.
The boundary conditions are discretised as in [1,5], this is, for both upwind and central fluxes, consider α = 1 for the numerical flux at the outer boundary and The current time step m + 1 is terminated when the stopping criterion is satisfied for some pre-defined small constant tol. Then the correspondent numerical solution is denoted by (Ê m+1 ). If we only perform one iteration (n = 0) we obtain the explicit method considered in [1]. If we perform two iterations (n = 0 and n = 1) we obtain a predictor-corrector type method.

Convergence result
We will show that, under a suitable stability condition, the solution of the iterative predictor-corrector scheme (4)- (6) converges to the solution of the underlying implicit method. The implicit method is defined as follows: given an initial where we consider the average approximationsẼ [m+1 We note that the the numerical solutions are defined implicitly, since the upwind fluxes involve the unknownsẼ m+1 Let h k be the diameter of the triangle T k ∈ T h , and h be the maximum element diameter, that is, We assume that the mesh is regular in the sense that there exists a constant τ > 0 such that for all T k ∈ T h , h k τ k ≤ τ , where τ k denotes the maximum diameter of a ball inscribed in T k . It may be proved (see [13]) that, for any u ∈ P N (T k ), the following trace inequality where f k is an edge of T k and C τ a positive constant independent of h k and N but dependent on the shape-regularity τ .
Let us now define the difference between two successive numeric values of the electromagnetic fields by for n = 0, 1, 2, · · · . The following theorem gives upper bounds for δ nÊ m+1 x k , δ nÊ m+1 y k and δ nĤ m+3/2 z k .

Theorem 1
The the solution of the iterative predictor-corrector scheme (4)- (6) converges to the solution of the method (7)-(9) provided that the stability condition of the underlying explicit method (i.e.,(4)-(6) taking only the iteration n = 0) is satisfied, that is (see [1]) with where C τ satisfies the trace inequality (11), C inv is a positive constant independent of h k and N , and Z k and Y k denote respectively the cell-impedance Z and the cellconductance Y inside the triangle T k ∈ T h .
Proof: The stability condition (12) ensures that δ 0Ê m+1 L 2 (Ω) and δ 0Ĥ Let us denote by F int the set of internal edges and F ext the set of edges that belong to the boundary δΩ. Let v k be the set of indices of the neighbouring elements of T k . For each i ∈ v k , we consider the internal edge f ik = T i ∩ T k , and we denote by n ik the unit normal oriented from T i towards T k . For each boundary edge f k = T k ∩ δΩ, n k is taken to be the unitary outer normal vector to f k .
Taking the difference of (4)-(6) between two successive iterations, n + 1 and n, and replacing u k , v k and w k by, respectively, δ nÊ m+1 x k , δ nÊ m+1 y k and δ nĤ m+3/2 z k and summing over all elements T k ∈ T h , we obtain Consequently, considering (11), we obtain Taking the following condition into account (that results from (12)) we conclude the proof. The next theorem establishes that the implicit method is second order convergent in time and arbitrary high order in space and so, with the previous result, we may conclude that the same occurs for the iterative explicit scheme.
Theorem 2 Let us consider the implicit leap-frog DG method (7)-(9) complemented with the discrete boundary conditions defined in Section 2.1 and suppose that the solution of the Maxwell's equations (1) complemented by (3) has the following regularity: where C inv and C τ are the positive constants defined in the previous theorem, then holds, where C is a generic positive constant independent of ∆t and the mesh size h. Proof: Follows the steps of the proof of Theorem 4.2 in [1].

Numerical results
To illustrate the theoretical results of the previous section, we consider the model problem (1)  To illustrate the order of convergence in space, we fix ∆t = 10 −5 , except when N = 4 where we consider ∆t = 10 −6 . In Fig. 2 we plot the discrete L 2 -error of theẼ x component of electric field depending on the maximum element diameter of each mesh, for different degrees for the polynomial approximation, for both central and upwind fluxes. The vertical and horizontal axis are scaled logarithmically. The numerical convergence rate is approximated by the slope of the linear regression line. For central flux, the numerical convergence rate is close to the value estimated in Theorem 2, O(h N ), and for upwind flux we observe higher order of convergence, up to O(h N +1 ) in some cases. Similar results were obtained forẼ y andH z .
To visualize the convergence order in time, the polynomial degree and the number of elements have been set to N = 8 and K = 800, respectively. The results plotted in Fig. 3 illustrate the first order convergency in time for the explicit leap-frog method

Modeling scattered electromagnetic wave's propagation through eye's structures
This work is part of a research project which aims to develop a cellular model of the human retina able to simulate different retinal/cellular conditions and how these changes are translated to an Optical Coherence Tomography scan [14]. Simulating the full complexity of the retina, in particular the variation of the size and shape of each structure, distance between them and the respective refractive indexes, requires a rigorous approach that can be achieved by solving Maxwell's equations. As the interest is to acquire the backscattered light intensity, we start this section by the scattered field formulation. Then we build up a two dimensional model which tries to represent a single nucleus of the outer nuclear layer (ONL) of the retina. The performance of our method is examined by simulating the light scattering in this 2D domain. The evolution of the scattering field intensity in time is obtained using the predictor-corrector DG method.

Optical Coherence Tomography (OCT)
Optical Coherence Tomography (OCT) its an imaging technology that produces high-resolution cross-sectional images of the internal microstructure of living tissue, widely used in ophthalmological exams. This technology's working principle is analogous to ultrasound, but it uses light instead of sound to locate subtle differences in the tissue being analysed. Discontinuities in the refractive index of the tissue give rise to light scattering, with some light backscattered to the detector. Factors such as the shape and size of the scatterer, wavelength of the incident light and refractive index differences have an impact on the amount of backscattered light. During a scan, the OCT machine directs a light beam into the retina and extracts, through interferometry, the backscattered light intensity of retinal structures and their depth location in an A-scan (see Fig. 4). By transversely moving the light beam, several A-scans can be collected into a cross-sectional image -a Bscan. Usually, several cross-sectional images are acquired by probing an azimuthal direction and combined into a volume.

The scattered field formulation
We can exploit the linearity of the Maxwell's equations in order to separate the electromagnetic fields (E, H) into incident fields (E i , H i ) and scattered components (E s , H s ), i.e., E = E s + E i and H = H s + H i .
Assuming that the incident field is also a solution of the Maxwell's equations we obtain in the same way as in [18], the scattered field formulation, with the source terms where i and µ i represent, respectively, the relative permittivity and permeability of the medium in which the incident field propagates in the absence of scatterers (in the background medium). Additionally, using this formulation it is straightforward to specify an incident wave using an analytic formula. The intensity of the light that hits the OCT detector defines the output signal. Hence, we are interested in computing the scattered field intensity,

Light scattering in the outer nuclear layer
We use our numerical model to simulate light scattering in the ONL. This layer has a special relevance among the retina's layers as it consistently presents the characteristics of diabetic macular edema [2,3]. The ONL is mostly populated by the cells bodies of light sensitive photoreceptor cells (rods and cons). Thus, we postulate that the main contribution to light scattering in this layer comes from the nucleus [16], as it is the biggest organelle in the soma and presents a high refractive index difference to the surrounding medium. As such, the ONL can be modelled as a population of spherical nuclei in an homogenous medium. As a proof of concept we present a simple simulation in a two dimensional square domain which contains circles that aims to represent, respectively, a single nucleus and three nuclei in the ONL. The permittivity inside the circles and in the background domain has different values. Let us consider equations (14)- (16), in Ω = (−1, 1) 2 , complemented with SM-ABC and null initial conditions. The absorbing boundary conditions are chosen for the model as they avoid undesirable reflections that invade the computational domain. In the first experiment we will consider the case of just one circle: C = {(x, y) ∈ Ω : x 2 + y 2 < 0.25}. In the second example we will consider the case of three circles: C 1 = {(x, y) ∈ Ω : x 2 + (y − 0.5) 2 < 0.01}; C 2 = {(x, y) ∈ Ω : x 2 + y 2 < 0.01}; C 3 = {(x, y) ∈ Ω : x 2 + (y + 0.5) 2 < 0.01}. In the experiments the relative permittivity and permeability and magnetic permeability are considered as constants, i = 1 and µ i = µ = 1. The electric permittivity is considered as a diagonal matrix with xx (x, y) = yy (x, y) = 1.2 for (x, y) inside the circles and xx (x, y) = yy (x, y) = 1 otherwise. For the incident wave we consider the planar wave E i y (x, t) = cos(10(x − t)). For the simulation we used with predictor-corrector DG method defined in Section 2, considering α = 0 (central flux) and the approximation polynomial degree N = 4. The time step was chosen to be ∆t = 0.002 and the final simulation time is T = 0.8. The meshes are illustrated in Fig. 5. The evolution in time of the scattered field  intensity (20) is plotted in Fig. 6. These results show that the scatterers are clearly identified. With this model, we can simulate more complex cellular structures only by changing the electric permittivity tensor .

Conclusions
We presented an iterative explicit leap-frog DG method for time dependent Maxwell's equations in anisotropic media, considering SM-ABC. The numerical scheme is fully explicit and converges to a second order in time implicit method. The results of a set of numerical experiments support the theoretical results. Moreover we developed a 2D model which simulates the light scattering by a single nucleus in the outer nuclear layer of the retina. This work was elaborated in the framework of a more general project with a real application (see [3,14]).