Back and Forth Error Compensation and Correction Method for Linear Hyperbolic Systems with Application to the Maxwell's equations

We study the Back and Forth Error Compensation and Correction (BFECC) method for linear hyperbolic PDE systems. The BFECC method has been applied to schemes for advection equations to improve their stability and order of accuracy. Similar results are established in this paper for schemes for linear hyperbolic PDE systems with constant coefficients. We apply the BFECC method to central difference scheme and Lax-Friedrichs scheme for the Maxwell's equations and obtain second order accurate schemes with larger CFL number than the classical Yee scheme. The method is further applied to schemes on non-orthogonal unstructured grids. The new BFECC schemes for the Maxwell's equations operate on a single non-staggered grid and are simple to implement on unstructured grids. Numerical examples are given to demonstrate the effectiveness of the new schemes.


Introduction
The goal of this paper is to study finite difference schemes for the Maxwell's equations that are based on the back and forth error compensation and correction (BFECC) method [6]. Extensive studies have been done on finite difference time domain (FDTD) schemes for the Maxwell's equation [23]. Compared with other methods, for example finite element schemes, FDTD methods are very efficient, easy to implement, and are able to model behaviors over all frequencies simultaneously [23]. The classical Yee scheme [25] is originally designed for uniform orthogonal grids. For non-uniform orthogonal grids, Yee scheme is known to be second order globally (though the local truncation error is first order) [18,19]. It can also be generalized for irregular nonorthogonal grids, such as the Nonorthogonal FDTD scheme [20], the Generalized Yee scheme [8] and the Overlapping Yee scheme [16]. These schemes require generation of (nonorthogonal or unstructured) staggered grids for E and H and the formulation and implementation on the unstructured staggered grids can be complicated. In this paper, we propose a simple finite difference scheme based on the BFECC method that requires very few modifications when changing from uniform non-staggered grids to unstructured non-staggered grids.
Back and Forth Error Compensation and Correction (BFECC) method is introduced in [6,7] to obtain a higher order scheme based on a lower order scheme for advection equations. Given a scheme for advection equations, the idea of BFECC method is to improve its accuracy by estimating using forward and then backward advections and correcting its leading order error. Suppose L is a r-th order linear scheme for scalar linear advection equations, where r is an odd integer, the in general the BFECC scheme based on L is (r + 1)-th order accurate, and is stable as long as scheme L has an amplification factor no more than 2, thus has a larger CFL number than L [6,7]. In this paper, we extend the BFECC method to linear hyperbolic systems, and show that similar accuracy and stability improvement can be achieved.
The BFECC method has been applied to level set interface computation and fluid simulations [6,7,11,12,13]. A two-step unconditionally stable MacCormack scheme and its generalization are developed in [22] for fluid simulations. The property that BFECC stabilizes even an unstable scheme (with its amplification factor no more than 2) is very helpful for systems because one doesn't have to compute the local characteristic information for constructing a low diffusion stable scheme. With the new extension to linear hyperbolic systems, we propose BFECC schemes for the Maxwell's equations which are second order accurate, easy to implement, and have larger CFL numbers than that of the classic Yee scheme [25]. Given the accuracy improving ability of the BFECC method, we propose to use a simple first order scheme that is based on the least square local linear approximation as the underlying scheme for BFECC on unstructured grids which is very easy to implement after being stabilized. Numerical examples show that the scheme remains to be second order on non-orthogonal grids.
The rest of the paper is organized as follows. In Section 2, we discuss the BFECC method for linear hyperbolic PDE systems and prove the stability and accuracy theorems. In Section 3, we apply BFECC to the Maxwell's equations. On uniform orthogonal grids, we use central difference and Lax-Friedrichs schemes as the underlying schemes for the BFECC method. Order of accuracy and CFL numbers for the corresponding schemes are discussed. On unstructured grids, we present a first order scheme based on the least square local linear approximation and use it as the underlying scheme. The divergence of the magnetic field and the perfectly matched layer [2] implementation are also discussed. Numerical examples are presented in Section 4. We conclude the paper in Section 5. A detailed error analysis of the BFECC applied to the central difference scheme is presented in the appendix.

BFECC method for homogeneous linear hyperbolic PDE systems with constant coefficients
In this section, we discuss the BFECC method for homogeneous linear hyperbolic PDE systems with constant coefficients. Denote u(x, t) the vector of unknown functions, where x = (x 1 , x 2 , ..., x d ) T ∈ R d and t ∈ R are the spatial and temporal variables. Consider a homogeneous linear hyperbolic PDE system with constant coefficients in the following form: where A i , i = 1, 2, ..., d are real constant matrices, and any linear combination d i=1 α i A i is diagonalizable with real eigenvalues. When all the coefficient matrices A i are symmetric, we say it is a symmetric linear hyperbolic system. We solve this system numerically with a finite difference scheme. For simplicity of discussion, we assume a uniform orthogonal grid is used and discuss the scheme in the whole space. Denote the mesh sizes ∆x = (∆x 1 , ∆x 2 , ..., ∆x d ), and ∆t n = t n+1 − t n (we omit subscript n when ∆t n is the same for all n). Denote the numerical solution where j = (j 1 , j 2 , ..., j n ) is the multi-index vector. Denote U n = U n j : ∀j the collection of numerical solution at all grid points at the time t n .
Suppose L is a numerical scheme for this system, i.e.
In this paper, all the schemes we discussed are linear schemes, i.e. L is a linear operator. We define L * the backward update step from t n+1 to t n by applying L to the time-reversed system: By applying the Back and Forth Error Compensation and Correction (BFECC) steps [6,7], we obtain a new scheme L BF ECC which updates the solution in three steps: (1) Solve forward. U n+1 = L U n + e (1) , where e (1) = 1 2 U n −Ũ n . U n andŨ n should have been the same if there were no numerical error. Therefore e (1) provides an estimate of the value lost during the forward step, which is then compensated to U n before performing the final forward step. In general, for linear advection equations, BFECC can improve the order of accuracy by one for odd order schemes and also improve stabilities of the schemes (see [6,7]). We establish similar results for systems of equations in the following theorems with the help of techniques in [26,7].
In the following discussion, we consider system (2.1) in d i=1 [0, 1] with periodic boundary conditions. And we assume the numerical scheme L is a linear scheme. Let ∆x j = 1 Nj for j = 1, 2, ..., d. The numerical solutions are then defined at any time on be the set for the dual indices of the finite Fourier series. Expand U n as a finite Fourier series where j ∈ D N and x j = (j 1 ∆x 1 , j 2 ∆x 2 , ..., j d ∆x d ).
Since scheme L is a linear scheme, the coefficients of the Fourier series get updated as Remark Note that scheme L is l 2 stable if the spectral radius ρ(Q L (k)) < 1 for all k ∈ F N or Q L (k) is diagonalizable and ρ(Q L (k)) ≤ 1 for all k ∈ F N .
Denote Q L * (k) the Fourier symbol matrix of L * . Then Fourier symbol matrix Q B for the BFECC scheme based on L is 2.1. Stability. In general, BFECC method improves the stability of an underlying scheme L for the scalar hyperbolic equation u t + v · ∇u = 0 [6,7]. It increases the CFL numbers of conditionally stable schemes (for example, the upwind scheme) and makes unstable schemes (for example, the central difference scheme) conditionally stable. We generalize this property of BFECC method for linear hyperbolic systems with constant coefficients. The result is summarized in the following theorem.
Theorem 2.1. Let L be a linear scheme for system 2.1. Suppose Q L and Q L * satisfies the following conditions Proof. We first show that λ j (Q B ) = 1 + 1 2 (1 − |λ j (Q L )| 2 ) λ j (Q L ) under the assumptions in the theorem, where λ j (Q L ) and λ j (Q B ) are eigenvalues of Q L and Q B , respectively, j = 1, 2, ..., d Let X = Re(Q L ) and Y = Im(Q L ). SinceQ L Q L = Q LQL , we have Since X and Y are diagonalizable with real eigenvalues and they commute, there is a basis set of real eigenvectors {v j } j=1,2,...,n that diagonalizes X and Y simultaneously. Then v i 's are also eigenvectors of Q L andQ L , and the corresponding eigenvalues are complex conjugate of each other, i.e. λ j (Q L ) =λ j (Q L ) for j = 1, 2, ..., d.
By the assumption Q L * =Q L , we get Q B = Q L I + 1 2 (I −Q L Q L ) , and thus therefore the conclusion of the theorem follows.

Remarks
1. Under the assumption of the theorem, Fourier symbol matrix Q B has a complete (real) eigenvector basis, so |ρ(Q B )| ≤ 1 implies l 2 stability. 2. Condition 1 follows the same assumption in the BFECC method for advection equations [7], condition 2 requires that the scheme treats backward temporal direction the same as forward temporal direction, and condition 3 usually follows from the diagonalizability of coefficient matrix of the system.
In particular, these assumptions on Q L and Q L * are satisfied for several classical schemes. For example, consider the following one dimensional hyberbolic system where A is diagonalizable with real eigenvalues. Let L represent the central difference scheme for this system, i.e.
Let λ = ∆t/∆x, then Here h = ∆x. We will continue to denote the spatial mesh size by h when there is no ambiguity. Let M represents the Lax-Friedrichs scheme for this system, then It is easy to see both schemes satisfy the assumptions of the theorem. 3. A easier-to-check (but more restrictive) alternative for condition 3 in the theorem is to require Q L being complex symmetric. This implies X and Y are real symmetric matrices, so they are diagonalizable with real eigenvalues. We will show in Section 3 that this condition is satisfied for the central difference scheme and Lax-Friedrichs scheme for the Maxwell's equations.

2.2.
Accuracy. In general, BFECC method improves the accuracy of odd order schemes for advection equations [6,7]. We extend this result to linear hyperbolic PDE systems with constant coefficients. Expand the solution into Fourier series and plug in system (2.1) to obtain where P (ik) is a matrix with entries that are homogeneous linear polynomials in ik with real coefficients, k = (k 1 , k 2 , ..., k d ) T . Therefore Assume ∆x 1 = ∆x 2 = ... = ∆x d = h, and fix ∆t/h during the mesh refinement. We first quote a theorem of Lax [15], Theorem 2.2. For the linear hyperbolic PDE system (2.1) with constant coefficients, a scheme L is r-th order accurate if and only if its Fourier symbol matrix Q L satisfies Here the O(|kh| r+1 ) term is a matrix whose entries are O(|kh| r+1 ) terms as h → 0. The "only if" part is stated in theorem 2.1 of Lax's paper [15] for linear hyperbolic systems with variable coefficients. When the coefficients are constant, Lax's argument can also be used to show that the "if" part is also true.
We have the following theorem, which is an extension of theorem 4 in [7] to homogeneous linear hyperbolic systems with constant coefficients. Theorem 2.3. Suppose Q L * (k) =Q L (k) for any k ∈ Z d and scheme L is r-th order accurate for system 2.1 with constant coefficient matrices, where r is an odd integer, then the BFECC scheme L BF ECC based on L is (r + 1)-th order accurate.
Proof. Since L is r-th order accurate, by the Theorem-2.2 [15], we have where Q r+1 (ikh) is a matrix with entries that are homogeneous degree r + 1 polynomials in ik with real coefficients.
By the assumption, ThenQ Therefore L BF ECC is a (r + 1)-th order accurate scheme.

2.3.
Alternative view of BFECC method for hyperbolic PDE systems. In some cases, we can view the BFECC method for systems as applying the BFECC method for advection equations to the Riemann invariants. Consider a one dimensional hyperbolic PDE system with constant coefficients where Λ is a diagonal matrix with eigenvalues of A as entries, define w = V −1 u, then the system is equivalent to Suppose now we have a r-th order scheme L for system (2.3), with r being odd, Note that this scheme updates each component W i independently from other components. Then it gives a r-th order scheme M for system-2.2, By theorem 4 in [7], applying BFECC to L produces an (r + 1)-th order scheme: Applying BFECC to M gives us therefore it is an (r + 1)-th order scheme for system (2.2) following the results for scalar equations for L B However, not all schemes for system (2.2) come from schemes for system (2.3) that update components of w independently. Also, it is numerically more costly to decouple the system, especially in multi dimensions. In these cases, theorem 2.1 and theorem 2.3 provide the stability and accuracy improvement results.

BFECC schemes for the Maxwell's equations
In this section, we discuss the BFECC schemes for the Maxwell's equations. We show that BFECC turns the central difference scheme and Lax-Friedrichs scheme into stable second order accurate schemes with larger CFL numbers than that of the Yee scheme on uniform rectangular grids. On non-orthogonal and unstructured grid, we discuss schemes based on least square linear approximation.
Consider the dimensionless Maxwell's equations in a medium with zero conductivity [23] where ǫ r and µ r are the relative permittivity and permeability, respectively. We assume they are constant in the following discussion.
To simplify the discussion for schemes, we use this Maxwell's equations in this section and refer to E ′ and H ′ as E and H.
For simplicity, denote E = E z , H = H y and we have: The central difference scheme on a uniform rectangular grid for the above system is: where λ = ∆t/∆x, E n j and H n j denote the numerical solutions E n j ≈ E(j∆x, t n ) and H n j ≈ H(j∆x, t n ). With periodic boundary conditions, E n j and H n j can be expanded uniquely as finite Fourier series: where k ∈ F N is the dual index, C n k and D n k are the Fourier coefficients for E and H, respectively.
Plug the finite Fourier series into the central difference scheme, we get where Q L is the Fourier symbol matrix. Since the spectral radius of Q L is greater than 1 for most k ∈ F N , the central difference scheme is a first order scheme that is unstable and cannot be directly used to solve the Maxwell's equations. Applying BFECC method to the central difference scheme stabilizes it and also improves the order of accuracy to second order. Solving Maxwell's equations in the backward temporal direction is equivalent to changing λ to −λ in the scheme, therefore Q L * = Q L . An easy calculation shows that Q L * Q L = Q L Q L * . The real and imaginary part of Q L are both diagonalizable with real eigenvalues. Therefore the conditions of theorem 2.1 and 2.3 are satisfied. We see that BFECC based on the central difference scheme is 2nd order accurate and l 2 stable if and only if ρ(Q L ) ≤ 2. Since the eigenvalues of Q L are 1 ± iλ sin(2πkh), the stability condition reduces to max k∈FN 1 + λ 2 sin 2 (2πkh) ≤ √ 3. Therefore BFECC based on the central difference scheme is a 2nd order accurate scheme and is stable if ∆t/∆x ≤ √ 3. An explicit calculation of the Fourier symbol matrix can be found in appendix A, which verifies that it is 2nd order accurate and stable if ∆t/∆x ≤ √ 3. Remark. In Section 4, we apply schemes discussed in this section to Maxwell's equations with variable permittivities. The schemes discussed in this section can be simply adapted to the case with variable permittivities. For example, for the following system where ǫ i and µ i are the permittivity permeability respectively at grid point x i . Other first order underlying schemes discussed in this paper can be similarly adapted to the variable coefficient case.
The central difference scheme is where λ x = ∆t/∆x and λ y = ∆t/∆y.
Expand H x , H y and E z into Fourier series: where (k, l) ∈ F N are dual indices and C n k,l , D n k,l and E n k,l are Fourier coefficients for H x , H y and E z , respectively.
Plug into the central difference scheme L, we get and Y = Im(Q L ). Similar to the one dimensional case, solving the equation backward in time amounts to switching the signs of λ x and λ y in the scheme. Therefore we have Q L * = I − iY = Q L , and Q L * Q L = Q L Q L * = I + Y 2 . I and Y are both symmetric real matrices, so they are diagonalizable with real eigenvalues. the conditions for theorem 2.1 and 2.3 are satisfied, and therefore BFECC based the central difference scheme is a 2nd order accurate scheme and is stable if ρ(Q L ) ≤ 2. The eigenvalues of Q L are It is satisfied if 2 ∆x is sufficient for stability, which implies a CFL factor An explicit calculation of the Fourier symbol matrix for the BFECC scheme is shown in appendix A.
3.3. BFECC based on the central difference scheme -three dimensional case. Similar to the one and two dimensional cases, we can also check the conditions of theorem 2.1 and theorem 2.3, and find that BFECC based on the central difference scheme is second order accurate and l 2 stable if Note that this still implies a CFL factor equal to one in three dimensions if ∆x = ∆y = ∆z. We summarize the results in the following theorem.
3.4. BFECC based on the Lax-Friedrichs scheme. We study BFECC based on the Lax-Friedrichs scheme M for the Maxwell's equations. In one dimension, the scheme is Write the one dimensional Maxwell's equations as Then the Fourier symbol matrix of the Lax-Friedrichs scheme is Q M = cos(kh)I + iλ sin(kh)A, wherek = 2πk is the angular wave number. It satisfies the conditions in theorem 2.1 and theorem 2.3, so BFECC based on the Lax-Friedrichs scheme is second order accurate and is stable if and only if |ρ(Q M )| ≤ 2, i.e., This is true if λ 2 ≤ 4. For two dimensional Maxwell's equations (3.4), write the equations as where λ x = ∆t/∆x, λ y = ∆t/∆y and λ z = ∆t/∆z. The stability and accuracy results are summerized as follows: Theorem 3.2. BFECC based on the Lax-Friedrichs scheme for Maxwell's equations in free space on uniform rectangular grid is 2nd order accurate. It is stable in the (1/∆x) 2 +(1/∆y) 2 and ∆t ≤ 7 2 min(∆x, ∆y); or (3) in three-dimensional case, ∆t ≤ 2 (1/∆x) 2 + (1/∆y) 2 + (1/∆z) 2 and ∆t ≤ √ 3 min(∆x, ∆y, ∆z).

3.5.
BFECC based on interpolation of the central difference and the Lax-Friedrichs schemes. The Lax-Friedrichs schems is more diffusive than the central difference scheme as the underlying scheme for BFECC. However, when there are discontinuities in the coefficients of the equations, the latter scheme may generate some numerical artifacts in the vicinities of the discontinuities. An interpolation between the two schemes could combine the strengths of both schemes. Let θ ∈ [0, 1]. A θ-scheme L θ is formally L θ = (1 − θ)L + θM, where L is the central difference scheme and M is the Lax-Friedrichs scheme for Maxwell's equations.
Using aforementioned notations, for one dimensional Maxwell's equations, the scheme is Its Fourier symbol matrix is which satisfies all the conditions in Theorem 2.1 and 2.3, and BFECC based on the θ-scheme is second order accurate. Note that Since f (x) = x 2 is convex, we have where Q L and Q M are the Fourier symbol matrices for the central difference and Lax-Friedrichs schemes, respectively. And the CFL number of the θ-scheme is between √ 3 and 2. Similarly, for the two dimensional Maxwell's equations (3.4), the θ-scheme is where U n i,j , A 1 and A 2 are defined as in Section 3.2. And its Fourier symbol matrix where q θ = 1 − θ + θ cos(kxhx)+cos(kyhy) 2 . The spectral radius ρ(Q θ ) satisfies In the inequality, we again use the convexity of f (x) = x 2 and the special form of q θ . Therefore, the constant in the CFL condition similar to (3.5) would be between √ 3 and 2. The analysis for three dimensional Maxwell's equations is similar, and the result is summarized as follows. where c θ ∈ [ √ 3, 2] depends only on θ.
3.6. Least square local linear approximation for non-rectangular grids. A special case of L θ is based on the linear least square fitting, which can also be used on irregular grids conveniently. In order to adapt to non-orthogonal grids, we consider some simple first order underlying schemes based on linear least squares. Least squares method significantly improves the robustness of polynomial approximation in multi dimensions. In WENO-type schemes for solving nonlinear conservation laws on unstructured meshes, least squares (high degree) polynomial fitting has been used, see for example [1,9].
To design an explicit scheme for the Maxwell's equations, we need approximations for spatial derivatives such as ∂Ez ∂x and ∂Hx ∂y at the time t n to update field variables E and H. A natural approach is to locally fit a linear function for each component of a field variable using the function values at a grid point and its neighbors, and then use the spatial derivatives of the linear function as approximations.
Consider for example the approximation of H x and its derivatives at a grid point (x i , y j ). Denote this point (x 0 , y 0 ). Suppose its neighboring grid points are (x 1 , y 1 ), (x 2 , y 2 ), ..., (x K , y K ), where K ≥ 2, and denote (H x ) i = H x (x i , y i ) for i = 0, 1, ..., K. A linear functionĤ x (x, y) =â +b(x − x 0 ) +ĉ(y − y 0 ) can be determined to fit the numerical values of H x at (x j , y j ), j = 0, 1, ..., K, by using least squares fitting. This is a local procedure, and has to be done at every point at which the scheme is evaluated.
We denote the approximated spatial derivatives at (x 0 , y 0 ) by ∂Ĥx ∂x and ∂Ĥx ∂y , and the approximated function value at (x 0 , y 0 ) byĤ x (x 0 , y 0 ) or Ĥ Similarly let ∂Êz ∂x , ∂Êz ∂y , ∂Ĥy ∂x , ∂Ĥy ∂y be the least square approximation of E z and H y 's partial derivatives at (x 0 , y 0 ). An explicit scheme similar to the central difference scheme is (for Maxwell's equations in two dimensions, (3.4)): where (E z ) n i,j denotes the numerical solution (E z ) n i,j ≈ E z (x i , y j , t n ), similarly for (H x ) n i,j and (H y ) n i,j . The set of grid points near (x i , y j ) used for least squares fitting in this paper are (x i , y j ), (x i±1 , y j ) and (x i , y j±1 ). When the grid is a uniform rectangular grid, then the above least square approximation for spatial derivatives is the central difference approximation if the same set of neighboring points are used, and (3.6) is just the central difference scheme. We refer to (3.6) as the least square central difference scheme.
One could use the least square approximated field values as well as the least square approximated derivatives in the scheme, i.e.
Here the subscript (i, j) and superscript n indicate that the approximation is done in a neighborhood of grid point (x i , y j ) using field values at time level t n . Note are weighted averages of field values at (i, j) and its neighbors, therefore this scheme is similar to the θ-scheme on uniform rectangular grids. When the grid is a uniform rectangular grid (possibly with ∆x = ∆y), scheme (3.7) reduces to the θ-scheme with θ = 0.8. We refer to this scheme as the least square θ-scheme. Both schemes are first order accurate, because the least square gradient approximation are first order accurate, and the least square field value approximation is second order accurate. Function approximation by least squares fitting have been well studied (see e.g. [4,24]). For completeness, we give a short discussion on the accuracy of the least squares fitting. Without loss of generality, we can assume (x 0 , y 0 ) = (0, 0). In a neighborhood of (0, 0) with radius O(h), rewrite function u(x, y) as u(x, y) = a + bx + cy + f (x, y) = l(x, y) + f (x, y), and function values at these grid points are collected in vector U = L + F , where L = l(x 0 , y 0 ), ..., l(x K , y K ) T and F = f (x 0 , y 0 ), ..., f (x K , y K ) T . Then we havê where || · || denotes the l 2 norm. In the above, we use the fact that A(A T A) −1 A T is an orthogonal projection. Suppose A is a (K + 1) × 3 matrix of full rank, so its smallest singular value σ 3 (A) > 0. Suppose σ 3 (A) ≥ Dh for some constant D > 0, then we have So the problem reduces to a geometric condition σ 3 (A) ≥ Dh for some D > 0 for the selected neighboring grid points. It can be easily verified that the rectangular mesh and the hexagonal mesh both satisfy this condition. For example, a rectangular grid of size h has σ 3 (A) = 2h, and a uniform hexagonal grid with edge length h has σ 3 (A) = √ 3h (using a grid point and its 6 adjacent grid points in the least squares fitting).
Next, to show the least square field value approximation is second order accurate, we notice that (x 0 , y 0 ) = (0, 0 and the first component of A(θ − θ) iŝ a +bx 0 +ĉy 0 − (a + bx 0 + cy 0 ) =â − a =û 0 − u(x 0 , y 0 ), whereû 0 is the least square field value approximation. Therefore Therefore the least square field value approximation is second order accurate. The order of accuracy result is summarized in the following theorem. Similar to the central difference scheme, the least square central difference scheme is usually numerically unstable. We can apply the BFECC method to improve the stability and accuracy. The least square θ-scheme is conditionally stable, and applying BFECC also improves its stability and accuracy. On a uniform rectangular grid, BFECC based on the least square central difference and least square θ-scheme are second order accurate and stable with CFL number √ 3 and a CFL number between √ 3 and 2, respectively. On non-uniform or non-orthogonal grids, our current analysis is not sufficient to prove the stability and order of accuracy. Numerical examples in Section 4 show that BFECC based on the least square θ-scheme is conditionally stable and second order accurate. We omit the examples for BFECC based on the least square central difference scheme, which is also second order in our experiments with smooth solutions (not reported here) but is likely to have numerical artifacts at places where the coefficients of the equations have jump discontinuities.
Remark. As will be discussed in Section 3.8, on a uniform rectangular grid, the central difference scheme and the BFECC scheme based on it preserve the divergence free property of the magnetic field. On a non-rectangular grid, the least square schemes and the corresponding BFECC schemes don't have this property. The flexibility of least square gradient approximation allows an option to reduce the divergence error. We can add a penalty term λ ∂Ĥx ∂x + ∂Ĥy ∂y 2 to the minimization functional of the least squares method, where λ ≥ 0 is a parameter. The Gauss's law for the electric field can similarly be incorporated into the least squares. We will study them in the future.
3.7. Point shifted algorithm for grid generation. It is often necessary to model curved material interfaces in computational eletromagnetics. The simplest treatment with a staircased approximation for the curved boundary can lead to large errors [5,23]. Local subcell methods [23] model curved interfaces/boundaries by modifying the update rule near them. In these cells, the integral form of the Maxwell's equations are usually used to update the field, e.g., the contour path method [10].
Using BFECC based on the least square central difference scheme (3.6) or BFECC based on the least square θ-scheme (3.7), we can locally deform the grid near a curved interface to conform with the interface, and avoid switching to the integral form of the Maxwell's equations in these deformed cells. In this section, we describe a simple point shifted algorithm [17] for shifting nearby grid points to the interface. It is used for numerical examples of scattering in Section 4.
Given a uniform rectangular grid in two dimensions, denote the grid points G rec = {(x i , y j ) : x i = i∆x, y j = j∆y, i = 0, 1, ..., N x , j = 0, 1, ..., N y }. Let C be a closed curve, e.g., the boundary of a scattering object. The point shifted algorithm shifts nearby grid points to the interface for distances less than half of the grid size so that the topological structure of the grid remains unchanged. And the new grid point set G C = {(x i ,ỹ j ) : i = 0, 1, ..., N x , j = 0, 1, ..., N y } conforms with curve C. It does so by finding the intersections of the grid lines and C, and shifts the nearest grid points to the intersection points. Find the nearest point (x i * , y j * ) in G rec to (x k ,ŷ k ), when there is a tie, break the tie arbitrarily. Set (x i * ,ỹ j * ) = (x k ,ŷ k ). end 4. Return G C .

Remark.
A optional smoothing step can be added after the point shift to make the grid deformation more smooth. Denote the uniform rectangular grid points x i,j and the point shifted grid pointx i,j , where i = 0, 1, ..., N x and j = 0, 1, ..., N y . First compute the point shift deformation d i,j =x i,j − x i,j . Second, copy d i,j tõ d i,j , and for every (i, j) such that d i,j = 0 (i.e. unshifted points), set This has the effect of smoothing out the point shift deformation. Third, assign new locations to the shifted grid points for i = 0, 1, ..., N x and j = 0, 1, ..., N y . Note the shifted grid points that lie on the curve C are unaffected by this smoothing step, only their neighbors get shifted in the smoothing step. This step can be repeated multiple times to smooth out the deformation to points that are further away from the curve C. Smoothing helps reduce grid deformation near the interface, and can be helpful when complicated interfaces are involved. Figure 1 shows examples of non-rectangular grids after applying the point shifted algorithm. The subfigure (a) is a uniform rectangular grid shifted to conform a circle without smoothing, the subfigure (b) is the same grid shifted to conform a circle, with a smoothing step, and the subfigure (c) is a uniform rectangular grid shifted to conform a more complicated curve, without smoothing. Grid (a) and (c) are use in the scattering numerical examples in Section 4. We didn't use the smoothing step since the material interfaces in our numerical examples are simple and solutions on grids without smoothing already has expected order of accuracy. Note that the topologies of these grids have not been changed by the algorithm, making the implementation almost as simple as on a uniform rectangular grid.  3.8. Divergence of the magnetic field. The magnetic field satisfies the divergence free condition in the Maxwell's equations as long as it does so initially. We show that the central difference scheme conserves the numerical divergence of the magnetic field when the grid is a uniform rectangular grid. Therefore, BFECC based on the central difference scheme also conserves the numerical divergence of the magnetic field.
The numerical divergence of the magnetic field at the time t n is: 2∆y .
Using the central difference scheme to update H x and H y , we get: ∆t.
Therefore (∇ · H) n+1 i,j = (∇ · H) n i,j . Similar arguments show that for the Lax-Friedrichs scheme and the θ-scheme, is a convex combination of ∇ · H at (x i , y j ) and its neighboring grid points (these two schemes only conserve i,j (∇ · H) i,j ). In particular, if (∇ · H) i,j = 0 for all i and j initially, then this property holds for all subsequent t n .
For irregular grids, the divergence free property is no longer guaranteed. But divergence penalty terms can be added to the minimization functional in the least square gradient approximation to reduce the divergence error, as discussed in Section 3.6.
3.9. Perfectly Matched Layer. Perfectly matched layers are commonly used as the absorbing boundary condition for problems in unbounded domains [2]. We consider combining the unsplit convolutional perfectly matched layer [14] with the BFECC method. Here we adapt the implementation in [21] and discuss it in the two dimensional case. The three dimensional case will be similar.
In the lossless domain, consider With the unsplit convolutional perfectly match layers, the equations in the perfectly matched layers are: where ζ w (t) = −σ w e −σwt u(t), w = x, y, u(t) is the unit step function, and σ x , σ y are chosen conductivity parameters in the perfectly matched layers (PMLs). For PMLs adjacent to a boundary perpendicular to the x-axis, we choose σ x > 0 and σ y = 0; and for PMLs adjacent to a boundary perpendicular to the y-axis, we choose σ y > 0 and σ x = 0. To implement BFECC in the perfectly matched layers, we first denote To update the field variables in the PMLs, we separate the terms dependent on the time t n from others in the convolution integrals (in order to stabilize them later by BFECC), approximate them using least squares, and obtain the least square central difference scheme: HereĤ x ,Ĥ y andÊ z are corresponding linear approximation functions obtained by the least squares fitting.
To apply BFECC to this scheme, we combine all the terms on the right hand side that involve spatial derivatives. For example, the equation for E z becomes ∆t is treated as a source term. In the first two steps of the BFECC method, we ignore this source term. It is only there in the third step of BFECC. We can see that this requires very little modification to the scheme used in the computational domain.
Similarly, we can also use BFECC based on the least square θ-scheme in the PMLs.

Numerical examples
4.1. 1D periodic solution. We consider the following periodic initial condition for the 1D Maxwell's equations The solution that satisfies the given initial condition is We solve the system with BFECC based on the central difference scheme from t = 0 to t = 0.6 with ∆t/∆x = 0.38, 0.98 and 1.7, and compare the numerical solutions with the exact solution.
The order of accuracy result is summarized in Table-1. The results confirm that BFECC based on the central difference scheme is second order accurate. Also note that the scheme is stable for ∆t = 1.7∆x, for which the classical Yee scheme becomes unstable.
The exact solution is We solve the system with BFECC based on the least square θ-scheme from t = 0 to t = 2.5 with ∆t/∆x = 0.25, and compare the solutions with the exact solutions. The problem is solved in four grids: (a) uniform rectangular grid; (b) nonrectangular grid obtained by a smooth perturbation from (a); (c) non-rectangular grid with a global circular grid deformation; and (d) non-rectangular grid with grid points shifted to a circular interface. The grids are shown in Figure-2 and the order of accuracy is shown in Table 2. We see the numerical orders of accuracy are all above 2, showing the effectiveness of the BFECC method on non-orthogonal grids.
(c) Non-rectangular with global circular deformation.
(d) Non-rectangular with local deformation.  4.3. Scattering by a dielectric cylinder. In this example, we solve the 2D Maxwell's equations in TM z mode with BFECC based on the least square θ-scheme for the scattering problem by a dielectric cylinder.
The incident wave is a z-polarized plane wave travelling in the x direction, i.e. (E z ) inc = sin(ω(x − t)), (H x ) inc = 0 and (H y ) inc = − sin(ω(x − t)), where ω = 2π/0.6 is the angular frequency. The computational domain is [0, 1] × [0, 1]. A dielectric cylinder with ǫ = 2.25 and µ = 1 and radius 0.24 is placed in the center of the computation domain. The surrounding medium has ǫ = 1 and µ = 1. Perfectly match layers are used as absorbing boundaries, and the total-field/scattered-field formulation [23] is used to introduce plane waves into the computational domain.
Two grids are used in computation: (a) a uniform rectangular grid is used and the material interface is approximated by stair-casing; and (b) a point shifted grid in which intersection points of the uniform rectangular grid and the material interface are computed and the closest rectangular grid points are moved to the intersection points, and see Figure 1 (a). We use a simple treatment for the material interface: if a grid point falls inside the dielectric cylinder, ǫ = 2.25 and µ = 1 are used during the update of E and H, otherwise, ǫ = 1 and µ = 1 are used. Other interface treatments will be studied in the future.
BFECC based on the least square θ-scheme is used instead of BFECC based on the least square central difference scheme is used. The larger numerical dissipation is helpful when there is material discontinuity. When BFECC based on the least square central difference scheme is used, there are small spurious oscillations presented in the numerical solution due to the material discontinuity.
Since the CFL condition for BFECC based on the least square θ-scheme only requires ∆t ≤ √ 3 √ (1/∆x) 2 +(1/∆y) 2 , here we take ∆t = ∆x = ∆y. Smaller ∆t values have also been experimented, giving similar results as presented here.
The numerical solution on the point-shifted grid at t = 3.8 is shown in Figure-3 and is compared with the analytic Mie solution [3] in Figure 4. BFECC based on the least square θ-scheme scheme is able to generate smooth solutions without any spurious oscillation. t = 3.8 is chosen since the solution seems to reach the steady state at this time. The scheme has also been tested for several thousands time steps (up to t = 12) and the solution remains stable.
The grid refinement analysis for numerical solutions on uniform rectangular grids and point shifted grids is shown in Table 3. Here the numerical solution on a 320 × 320 grid is taken as the accurate solution, and all errors (in l 2 ) are computed with respect to this numerical solution. We can see that the BFECC scheme essentially achieves second order accuracy.    4.4. Scattering by a dielectric object of complicated shape. In this example, BFECC bases on the least square θ-scheme is used to solve a scattering problem by a dielectric object of more complicated shape. The material setup is the same as in the previous example. Notice that the object has sharp corners and cavities inside, as shown in Figure 5. The two grids used for computations are (a) a uniform rectangular grid with staircasing approximation for material interface and (b) a point shifted grid, see Figure 1 (b). The contour plot of E z at t = 3.6 is shown in Figure 5. Taking the 320 × 320 numerical solution as the reference, the numerical errors (in l 2 ) are shown in Table 4. Again we see that the BFECC scheme is stable and has second order accuracy.

Conclusion
We study the Back and Forth Error Compensation and Correction (BFECC) Method for linear hyperbolic PDE systems and establish the stability and accuracy properties of BFECC for homogeneous linear hyperbolic systems with constant coefficients. The method is then applied to the Maxwell's equations. On uniform orthogonal grids, BFECC based on the central difference and BFECC based on the Lax-Friedrichs schemes are proved to be second order accurate and have larger CFL number than the classical Yee scheme. On non-orthogonal or unstructured grids, the BFECC method is applied to a first order scheme based on least square gradient approximation. Numerical examples demonstrate the effectiveness of the BFECC schemes for Maxwell's equations. In particular BFECC based on the least square central difference scheme or the least square θ-scheme is easy to implement on non-orthogonal grids, has larger CFL numbers and second order accuracy in the numerical examples we have tested. We plan to test the BFECC schemes on unstructured grids with adaptive refinement for more complicated application problems in the future.

Acknowledgement
The authors thank Jinjie Liu for helpful discussions on the Yee scheme.
Therefore ∆t/∆x = λ < √ 3 ensures l 2 stability for the BFECC scheme. Accuracy Write and plug in the Maxwell's equations, we get: where matrix G is defined by the last equality. Calculate the matrix exponential, we get: While with BFECC based on the central difference scheme, we have: Note λh = ∆t, we see: By the Theorem-2.2, we see BFECC based on the central difference scheme is a second order accurate scheme.

  
It is then easy to see the scheme is l 2 stable if and only if max s x k ,s y l |λ 2,3 | ≤ 1, which is not true since |λ 2,3 | 2 = 1 + λ 2 x (s x k ) 2 + λ y (s y l ) 2 > 1 for s x k = 0 or s y l = 0. Therefore the central difference scheme is unconditionally unstable.