Fourth-Order Accurate Compact Scheme for First-Order Maxwell’s Equations

We construct a compact fourth-order scheme, in space and time, for the time-dependent Maxwell’s equations given as a ﬁrst-order system on a staggered (Yee) grid. At each time step, we update the ﬁelds by solving positive deﬁnite second-order elliptic equations. We develop compatible boundary conditions for these elliptic equations while maintaining a compact stencil. The proposed scheme is compared computationally with a non-compact scheme and with a convolutional dispersion relation preserving (DRP) scheme


Introduction
Time-dependent Maxwell's equations have been established over 125 years as the theoretical representation of electromagnetic phenomena [9,10].However, the accuracy requirement of high-frequency simulations remains a challenge due to the pollution effect.This was investigated first for the convection equation [8] and later for the Helmholtz equation [1,5].
We consider a uniform discretization of = [0, 1] 3 , x = y = z = h = L N where L is the length of the domain in a direction and N is the number of grid points in that direction.The wavenumber is given by ω = f v where f is the frequency of the wave and v is its velocity.The pollution effect states that, for a p-order time accurate scheme, the quantity ω p+1 h p t should remain constant for the error to remain constant as the wave number varies.A similar phenomenon applies also to spatial errors [1,5].This effect motivates the need for higher-order schemes, i.e., larger p.Using straightforward central differences, higher-order accuracy requires a larger stencil.This has two disadvantages.Firstly, the larger the stencil, the more work may be needed to invert a matrix with a larger bandwidth and more non-zero entries.Even more serious are the difficulties near boundaries.A large stencil requires some modification near the boundaries where all the points needed in the stencil are not available.This raises questions about the efficiency and stability of such schemes.In addition, higher-order compact methods can achieve the same error while using fewer grid nodes, making them more efficient in terms of both storage and CPU time than non-compact or low-order methods (see e.g., [3] for the wave equation, [13] for Helmholtz equation, [17,18] for Maxwell's equations, and references therein).The Yee scheme (also known as finite-difference time-domain, or FDTD) was introduced by Yee [16] and remains the common numerical method used for electromagnetic simulations [15].Although it is only a second-order method, it is still preferred for many applications because it preserves important structural features of Maxwell's equations that other methods fail to capture.One of the main novelties presented by Yee was the staggered-grid approach, where the electric field E and magnetic field H do not live at the same discrete space or time locations but at separate nodes on a staggered lattice.From the perspective of differential forms in space-time, it becomes clear that the staggered-grid approach is more faithful to the structure of Maxwell's equations that are dual to one another and hence E, H naturally live on two staggered, dual meshes [14].The question of staggered versus non-staggered high-order schemes has been studied in [6] (see also [15, pp. 63-109]), where it has been shown that, for a given order of accuracy, a staggered scheme is more accurate and efficient than a non-staggered scheme.
In this paper, we present a new compact fourth-order accurate scheme, in both space and time, for Maxwell's equations in three dimensions using an equation-based method on a staggered grid.At each time step, we update the solution by solving uncoupled (positive definite) second-order elliptic equations using the conjugate gradient method.This procedure involves a non-trivial treatment at the boundaries on which boundary conditions are not given explicitly but have to be deduced by the equation itself while maintaining the compact stencil.While the development of the scheme is done in 3D, the simulations are only two-dimensional.

Maxwell's Equations
Let E, H, D, B, and j be (real-valued) vector fields, defined for (t, x) ∈ [0, ∞) × , which take values in R 3 .E and H are the electric and magnetic field strength, respectively, D is the displacement, B is the magnetic flux density.In addition ρ and j are the densities of the extraneous electric charges and currents, c is the speed of light in vacuum, and ε and μ are the electric permittivity and magnetic permeability, respectively.The Maxwell equations on , in first-order differential vector form, are given by where ∇× denotes the curl operator and ∇• denotes the divergence operator with respect to x.The first pair of equations (1) represents the Ampére law (law of magnetic circulation) and Faraday law (law of electromagnetic induction).The second pair of equations (1) represent

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.the Gauss laws of electricity and magnetism, respectively.The third pair of equations ( 1) is the constitutive relations, where we are assuming that the permittivity ε and permeability μ are positive scalar constants which corresponds to a homogeneous and isotropic dielectric medium that fills the region .The constitutive relations in (1) are the same as they would have been in the case of static fields.Hence, the equations (1) are considered only for the range of frequencies below the threshold values where the dispersion of material characteristics starts to manifest itself.For many actual materials, these threshold frequencies are high and this limitation presents little to no loss of generality.For further detail on the derivation and interpretation of Maxwell's equations (1) we refer the reader to [10,Chapter IX].
By taking the divergence of the Ampére law and substituting the Gauss law of electricity, we arrive at the continuity equation for the charges and currents: This equation represents the physical property of the conservation of the electric charge.
At the same time, as it is derived from Maxwell's equations, we conclude that it provides a necessary solvability condition for the system (1).Also, by taking the divergence of the Faraday law we obtain ∂ ∂t (∇ • B) = 0, which implies that as long as the magnetic field is solenoidal at the initial moment, it will remain solenoidal at all subsequent times as well.At t = 0, we assume that smooth initial conditions are given for D and B. In addition, we assume that = [0, 1] 3 and n × E = 0 on ∂ (see [11,Section 8]).
By differentiating the Ampére law in (1) with respect to time t, substituting ∂ ∂t ∇ × B = μ ∂ ∂t ∇ × H from the Faraday law, and employing the identity ∇ × ∇ × E = − E + ∇(∇ • E) along with the Gauss law of electricity and the constitutive relations, we obtain the (inhomogeneous) wave equation for the electric field: ( A similar argument yields the wave equation for the magnetic field: From Eqs. ( 2) and (3) we conclude, in particular, that the propagation speed of electromagnetic waves in a homogeneous isotropic dielectric is equal to c √ εμ .For convenience, we introduce a simplified setup with simplified notation and rescale Maxwell's equations accordingly.Namely, from here on we assume that the propagation speed in the material is equal to one so that c = √ εμ.Let us also denote J = 4π c j, absorb the coefficient of 4π into ρ: ρ → 4πρ, and define Z = μ ε .Then, the Maxwell's system of Eq. ( 1) is recast as Equation ( 4) are to be supplemented with the initial conditions at t = 0: and boundary conditions n × E = 0 on ∂ .Moreover, the second order wave Eqs. ( 2) and (3) reduce to and respectively.
To determine the number of imposed boundary conditions, we need to examine the eigenvalues for the reduced one-dimensional equations at that boundary.A positive eigenvalue indicates a variable that enters from outside the domain and so needs to be specified.A negative eigenvalue indicates the variable is determined from the inside and so cannot be imposed and instead is determined by the interior equations.A zero eigenvalue is more ambiguous.
Remark 1 An imposed boundary condition must contain only derivatives of a lower order than the differential equation.Thus, for a first-order system, all boundary conditions must be of Dirichlet type.A Neumann condition is not considered an imposed boundary condition.
To derive the boundary conditions for E and H, we analyze the characteristics of the homogeneous counterpart to system (1), ρ = 0 and J = 0, in the direction normal to the boundary (cf.[7,Section 9.1]).The two divergence equations can be considered as constraints on the initial data and do not affect the characteristics.Let w be the vector (E x , E y , E z , H x , H y , H z ) t .Then, considering only the x space derivatives (without loss of generality) we arrive at The eigenvalues of A are c = ± 1 μ (twice) and 0 (twice).Thus, two boundary conditions need to be imposed, two are determined from the interior and two are not clear.We shall see that we impose E = 0 (two conditions) and H ⊥ = 0 (one condition).
The structure of the paper is as follows: In Sect. 2 we define all the necessary quantities with a detailed description of the grid placement.We then describe the type of the given boundary conditions, In Sect. 3 we present our scheme in space and time and then discuss. in more detail.the Neumann boundary conditions.Section 4 is devoted to the stability analysis of the scheme.Finally, Sect. 5 gives numerical confirmation of the theoretical results.

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.

Preliminaries
We consider a uniform discretization of = [0, 1] 3 , x = y = z = h and h t is the time-step.We introduce the following notations: denotes the Laplacian with respect to x.
-For any operator T : R M → R N , T := sup To discretize the equations, we introduce a staggered mesh in both space and time as in the Yee scheme [16].With this arrangement, all space derivatives are spread over a single mesh width, and the central time and space derivatives are centered at the same point, similar to that of the Yee scheme.E is evaluated at time nh t while H and J are evaluated at time (n+ 1 2 )h t .For the spatial discretization, we define the following meshes: For each s = x, y, z we define the interior meshes: By convention, x 0 = y 0 = z 0 = 0 and x N = y N = z N = 1.We denote 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Fig. 1 The grid points of H y h projected on the (x, y) plane in the case N = 5 See Fig. 1 for an illustration of the sets H y h projected on the x, y plane.The discretized functions E h , H h are then defined as: With this arrangement, the boundary condition n × E = 0 on ∂ E h implies that In short, E x = 0 when y = 0, 1 and z = 0, 1; E y = 0 when x = 0, 1 and z = 0, 1;E z = 0 when x = 0, 1 and y = 0, 1.

Boundary Conditions for Time Derivatives
We will now derive the boundary conditions for the time derivatives of the fields E and H rather than the fields themselves.We do so since time marching with our compact scheme derived in Sect. 3 involves solving elliptic equations for the time derivatives of the fields (see Eqs. (12a) and (12b)).By ( 9), we have the following Dirichlet-type boundary conditions: 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.
The first three rows follow from (10).The other six rows are an implication of (11) Moreover, differentiating the Gauss law of electricity ∇ • E = ρ ε with respect to time and using (9) we see that, for the remaining variables, the Neumann-type conditions hold: The Faraday law ∂ t H = − 1 Z ∇ × E yields the following Dirichlet-type boundary conditions: In Table 1, we provide the values of the second-order derivatives at the boundaries of that we need for subsequent analysis and that can be derived from Eqs. ( 10) and ( 11) by differentiation.
The Ampère law Substituting the values from Table 1 into these equations, we obtain the following six Neumann-type boundary conditions: 3 The Scheme

Equation Based Derivation
We now extend the second-order accurate Yee scheme to a fourth-order accurate scheme, in both space and time, while maintaining the compact Yee stencil.Using Maxwell's equations, the idea is to use a Taylor expansion in time to the next order and then replace the resulting third-order time derivatives with space derivatives.The fourth-order Taylor expansion applied to (4) and combined with Eqs. ( 5) and ( 6) yields: Similarly, we have for H: Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Thus, at each time step, we obtain several uncoupled modified Helmholtz equations which are positive definite.Note: the equations might be coupled at the boundary for more complicated situations. where Equation ( 12) are supplemented by the boundary conditions for the functions δ t E and δ t H, as discussed in Sect. 2.

Modified Helmholtz Equation
In a Cartesian coordinate system, each of the Eqs.(12a) and (12b) gets split into three scalar modified-Helmholtz equations: We use the fourth-order accurate compact finite-difference scheme [13] for solving elliptic equations of the type (14).Consider an equally-spaced mesh of dimension N x × N y × N z and size h in each direction.We denote by D x x , D yy , and D zz the standard second-order central-difference operators and define Then, the fourth-order accurate scheme for ( 14) is given by [13]: If F is known at the grid nodes to fourth-order accuracy, then ( 14) can be simplified further to Letting κ 2 = 24 h 2 t in ( 12) and ( 15) gives rise to six elliptic equations of type ( 14) with φ and F specified in Table 3. Once the discretization in space has been applied (Sect.3.4), Eq. ( 12) lead to six positive definite linear systems to be solved independently by conjugate gradients.Note: it is difficult to use multigrid since the grid dimension for various equations is not the same and hence not all equations can have 2 n grid nodes in every direction.

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Table 3 The expressions for φ and F in (14) (P(J) is given by ( 13))

Neumann Boundary Conditions
We reemphasize that Neumann conditions for a first-order system are not imposed boundary conditions.They are rather an implication of the PDEs in the interior.Let = [0, 1] 3 .We define ghost nodes outside the numerical grid so that the Neumann conditions are satisfied with fourth-order accuracy (see Fig. 2).
We propose two methods for the implementation of the Neumann boundary conditions.Method 1 (equation-based) Consider, without loss of generality, Eq. ( 14) with the Neumann boundary condition ∂ x φ(0, y, z) = g(y, z).Differentiating (14) with respect to x (cf.[13]), we have: 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Next, we show how to approximate the Neumann-type boundary conditions to fourth-order accuracy for ∂ t E using Method 1 and for ∂ t H using Method 2.
Proof φ x (0, y, z) equals to ∂ t E x at x = 0, which has already been specified.Moreover, we have Differentiating (19) twice with respect to y and once with respect to t, we arrive at If x = 0, then ∂ t yyy E x = ∂ tzyy E x = 0, and hence if x = 0, then ∂ xyy φ = ∂ tyy (ρ/ε).Differentiating (19) twice with respect to z and once with respect to t, we get If x = 0, then ∂ tzzz E x = ∂ tyzz E x = 0, and hence at x = 0 we have This quantity, in turn, can be derived via the Ampère law And the latter is known through the Neumann boundary condition for The second assertion follows immediately by the same argument.
By repeating the previous proof for the remaining faces of the cube = [0, 1] 3 and field components E y and E z , we obtain the following Corollary 1 Neumann boundary conditions for Eq.(12a) admit a compact fourth-order accurate discretization using Method 1.
Next, we will show how to use Method 2 to build a compact discretization of the Neumann boundary conditions for ∂ t H. Lemma 2 Let φ(x, y, z) = ∂ t H y be evaluated at some fixed t 0 .Then, the term φ x x x at x = 0, 1 can be approximated to fourth-order accuracy.

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.(2024) 100:31 Proof Assume without loss of generality that Z = 1.The equation )) and boundary conditions ∂ t E z = 0 at x = 0, 1 imply that at x = 0, 1 Then, we use ∂ t ∇ • H = 0 and replace ∂ x H x with −∂ ty H y − ∂ tz H z .This yields: To complete the proof, we have to approximate ∂ tyyx H y and ∂ txyz H z .This is done as follows.
We use (20) again to obtain At x = 0, 1, the first equation of (11) implies that ∂ t yyy H x = 0.
To approximate ∂ txyz H z , we use the y-component of the Ampère law (4a): Therefore, at x = 0, 1 we have: Since ∂ t H x = 0 at x = 0, 1, we finally derive: This completes the proof.
By repeating the previous procedure for H x and H z , we obtain the following corollaries.

Corollary 2 Neumann boundary conditions for Eq. (12b) admit a compact fourth-order accurate discretization using Method 2.
Corollary 3 Assume that ρ = 0 and J = 0.Then, all Neumann boundary conditions obtained for ∂ t E and ∂ t H are homogeneous and all the derivatives in (17), ( 18) vanish.In particular, 2 ) is a fourth-order accurate approximation to the homogeneous Neumann condition at the boundary point x (see also Sect.5).
We emphasize that the right-hand side, F, needs to be approximated with fourth-order accuracy.This is done using a fourth-order Padé approximation for the curl operator (see details in "Appendix A.1").Then, we solve Eqs. 1, 2, and 3 from Table 3 using the scheme (15), and subsequently solve Eqs. 4, 5, and 6 from Table 3 using the scheme (16).

Definition 4
For any s = x, y, z and G = H s , E s , let

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.
(2024) 100:31 Page 13 of 25 31 We omit the superscript G whenever it does not lead to any ambiguity.
Hereafter, r = h t h will denote the CFL number.As shown in "Appendix A.2", if r < 3 + √ 21, then P G 2 and P G 1 are symmetric positive definite matrices and the following estimates hold: In particular,

The Numerical Scheme
Let ⎞ ⎠ be the block matrix representing a spatial Padé fourth-order finite-difference approximation of the curl operator.See "Appendix A.1" for the full definition of the matrix curl h .
Let E n h and E n h be given on E h and let H n+ 1 2 h be given on H h .To advance in time, we take two steps.

2). Extend H n+3/2 h
to H h using the appropriate boundary condition.

Stability Analysis
We assume, with no loss of generality, that ρ = 0, J = 0, and Z = 1.The scheme can therefore be written compactly: We approximate in Eqs.(23a) and (23b) using the standard second-order difference Laplacian h (same as in ( 15)).Then, we recast (23a), (23b) as: where the operators P 1 and P 2 are introduced in Definition 4. Next, assume the solution is in the form of a plane wave: 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.Let Q = −P −1 1 P 2 curl h and let λ denote an eigenvalue of Q.We substitute the plane wave solution into (24) and derive Therefore, In particular, |h t λ| ≤ 2 implies that |σ | ≤ 1.Then, using the definition of the matrix A given by Eq. ( 35) and Remark 2 (see "Appendix A.1"), we can derive the stability condition in the following form: (see "Appendix A.2" for detail).Therefore, the scheme is stable provided that r ≤ min The inequality gives rise to a sufficient stability condition In the following section, we will examine a two-dimensional reduction of Maxwell's equation.A modification of Remark 2 readily implies the stability condition (cf. [18]).The numerical results summarized in

Numerical Simulations Data
For the computations, we consider the scaled two-dimensional TM system without any current or charges.Thus, the equations are reduced to 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Since H is given at time moments n + 1/2 we need to specify its initial condition at the time t = t 2 .By Taylor series, H t (0) is given by (4b) and H tt (0) is given by ( 6).Differentiating (6) yields: where, again, H t is given by (4b).Similarly, where H tt is given by ( 6).At the nodes on or next to a boundary, we have either Dirichlet or Neumann conditions as given in Table 2.For (28), we have assumed J = 0.For the computations, we further assume We consider the analytical solutions We use the mean absolute error where the subscripts "numer."and "true" refer to the numerical and analytical solution, respectively.The ghost nodes are easily defined using Corollary 3. The ghost points for H x (see Fig. 2) are defined similarly.

Comparable Schemes: C4, NC, AI h
We compare our proposed scheme, denoted by C4, with the following two schemes, which are second-order accurate in time.These schemes are obtained according to the Yee updating rules: ), Content courtesy of Springer Nature, terms of use apply.Rights reserved.
where D x and D y are finite difference operators that approximate the first derivative on different stencils.For the derivatives in the x direction, we consider a general stencil whereas for the derivatives in the y direction we use the transposed stencil S t .In particular, For grid nodes near the boundary, we use the standard fourth-order accurate one-sided finitedifference approximation of the first derivatives.By a Taylor expansion, we obtain second-order accuracy provided that and fourth-order accuracy if the additional constraints hold: We define the following stencils: The non-compact fourth-order accurate scheme, denoted by NC, exploits the stencil K 4 with a = 0. Its order of accuracy is O(h 4 ) + O(h 2 t ).Our last comparable scheme is a data-driven scheme, AI h .We consider this approach since it has been recently shown to reduce the numerical dispersion for the wave equation [12].The general framework for a data-driven scheme is described as follows.Let {(X n , Y n )} n be n-tuples of data-points such that for any n, X n ∈ R m and Y n ∈ R l .We define a network N a as a function from R m → R l which depends on parameters a.We define a loss function of the form and obtain the optimal parameters a and the corresponding optimal network N a .The minimization process can be done using variations of the gradient descent method.The data-driven scheme AI h uses the stencil K 2 , where the free parameters are obtained by a minimization process over the given training data.This scheme is of order O(h 2 ) + O(h 2 t ).The details of the minimization process and the selected training data are given in "Appendix A.3".

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved., Z = 1

Observations
In Table 4, we verify the fourth-order accuracy of our scheme C4 and provide the rates of grid convergence for the other two (NC, AI h ) schemes as well.
In Table 5, we examine the effect of the CFL number on the mean error and verify the results of our stability analysis (Sect.4).For the points-per-wavelength (PPW) ratio of 64 (k x = k y = 2), the mean error is smaller for scheme C4, as expected.Moreover, the NC scheme is more accurate than AI since it is of a higher spatial order, For the PPW ratio of approximately 6.4, the AI scheme shows a similar error to that of C4 even though it is of only second order.This is because the AI scheme was trained on coarse grids with a small number of points per wavelength.
In Fig. 3, we examine the mean error as a function of the wave number.As expected, as the PPW ratio decreases, the C4 scheme no longer outperforms the other schemes.This is consistent with the Taylor approximation where the local truncation error depends on higherorder derivatives that increase for shorter wavelengths (larger wavenumbers k 2 x + k 2 y ).

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.In Figs. 4 and 5, we examine how the error behaves as a function of the time step.As expected, scheme C4 is more accurate when the wave number is low ( k 2 x + k 2 y = 2) and the local truncation errors are relatively small.However, we see that for k x = k y = 11 the non-compact scheme is more accurate for r = 1 , the scheme NC is more accurate than AI while the situation is the opposite for r = 5 6 √ 2 . This is not a surprise since the dispersion effects may grow as the time step increases [2].
For higher wave numbers ( k 2 x + k 2 y = 29 and k 2 x + k 2 y = 51), the fourth-order accuracy of C4 and the spatial fourth-order accuracy of NC lose their advantage due to the pollution effect, and as the wave number increases the AI scheme demonstrates a similar accuracy even though it is only second-order accurate in both time and space.

Conclusions
We have constructed a compact implicit scheme for the 3D Maxwell's equations.The scheme is fourth-order accurate in both space and time and can be useful for long-time simulations, including scenarios where the PPW ratio is not very high.The scheme is compact and main-

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.
tains its accuracy near the boundaries using equation-based approximations.The temporal grid is staggered, which allows us to solve 3 scalar uncoupled elliptic equations at each half-time step.That can be done in parallel using the conjugate gradient method.The elliptic equations are strictly positive definite, which yields rapid convergence of the iterative scheme.With a fixed CFL number, approximately three iterations of the conjugate gradient are required for convergence regardless of the grid size.This is likely explained by the fact that the amount of positivity in Eq. ( 12) increases as the grid length decreases because the quantity 24 h 2 t becomes larger as h decrease and the CFL number remains fixed.Overall, the method is efficient in both CPU and memory requirements.Although we have tested numerically only 2D examples, the implementation in 3D is straightforward.
Since the method we have presented is a finite-difference scheme, the pollution effect cannot be avoided completely.For higher wavenumbers, lower-order methods obtained by minimizing a certain loss function over the stencil parameters (such as data-driven) can outperform higher-order schemes.

A.1 Operator Details
Let ⊗ denote the Kronecker product of matrices and let I , J , K be index sets.Let Next, we consider a fourth-order accurate approximation for the first derivatives: Assume that f (x) is known at N points x 0 , x 1 , .., x N −1 .Then, one estimates f (x + h/2) at the N − 1 points x 1 2 , .., x N −3/2 to fourth order as:

123
Content courtesy of Springer Nature, terms of use apply.Rights reserved.
We define the following operators from I Finally, we define a fourth-order approximation of the curl operator as a matrix that operates on the vectorized tensor of dimension and by [18], A −1 ∼ 6 5 .
We recall the standard second-order finite difference matrices for the second derivative.We define the operators D x x , D yy , D zz from I × J × K to I × J × K using the square matrix

A.2 Operator Estimates
Let r = h t h be the CFL number.Consider the finite-difference operators and 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.
Let σ (•) denote the spectrum of a given operator.The inclusion The operator where and −4 ≤ S ≤ 12. Hence, As a result, r < 3 + √ 21 implies that −P 2 , P 1 are positive definite and

A.3 The Data-Driven Scheme AI h
We follow the approach of [12].We wish to find the optimal parameters a, b, d for evaluation of the first derivative using the stencil in Yee updating rules (see Sect. 5.2).The network will be a function that takes a solution E n , H n+1/2 and returns E n+1 , H n+3/2 using Yee updating rules.We generate a data set of analytical solutions using (30a) as follows.We fix h = 1 16 , CFL number r and the corresponding h t .We also fix the final time T .The number of spatial grid points is then denoted by N and the number of times steps is denoted by N t .For any (k x , k y ) in {12, 13, 14, 15} 2 we generate analytical solutions E z , H x , H y defined by (30a).For any 0 ≤ n ≤ N t we let 123 Content courtesy of Springer Nature, terms of use apply.Rights reserved.1.

Terms and Conditions
Springer Nature journal content, brought to you courtesy of Springer Nature Customer Service Center GmbH ("Springer Nature").Springer Nature supports a reasonable amount of sharing of research papers by authors, subscribers and authorised users ("Users"), for small-scale personal, non-commercial use provided that all copyright, trade and service marks and other proprietary notices are maintained.By accessing, sharing, receiving or otherwise using the Springer Nature journal content you agree to these terms of use ("Terms").For these purposes, Springer Nature considers academic use (by researchers and students) to be non-commercial.These Terms are supplementary and will apply in addition to any applicable website terms and conditions, a relevant site licence or a personal subscription.These Terms will prevail over any conflict or ambiguity with regards to the relevant terms, a site licence or a personal subscription (to the extent of the conflict or ambiguity only).For Creative Commons-licensed articles, the terms of the Creative Commons license used will apply.We collect and use personal data to provide access to the Springer Nature journal content.We may also use these personal data internally within ResearchGate and Springer Nature and as agreed share it, in an anonymised way, for purposes of tracking, analysis and reporting.We will not otherwise disclose your personal data outside the ResearchGate or the Springer Nature group of companies unless we have your permission as detailed in the Privacy Policy.While Users may use the Springer Nature journal content for small scale, personal non-commercial use, it is important to note that Users may not: use such content for the purpose of providing other users with access on a regular or large scale basis or as a means to circumvent access control; use such content where to do so would be considered a criminal or statutory offence in any jurisdiction, or gives rise to civil liability, or is otherwise unlawful; falsely or misleadingly imply or suggest endorsement, approval , sponsorship, or association unless explicitly agreed to by Springer Nature in writing; use bots or other automated methods to access the content or redirect messages override any security feature or exclusionary protocol; or share the content in order to create substitute for Springer Nature products or services or a systematic database of Springer Nature journal content.
In line with the restriction against commercial use, Springer Nature does not permit the creation of a product or service that creates revenue, royalties, rent or income from our content or its inclusion as part of a paid for service or for other commercial gain.Springer Nature journal content cannot be used for inter-library loans and librarians may not upload Springer Nature journal content on a large scale into their, or any other, institutional repository.These terms of use are reviewed regularly and may be amended at any time.Springer Nature is not obligated to publish any information or content on this website and may remove it or features or functionality at our sole discretion, at any time with or without notice.Springer Nature may revoke this licence to you at any time and remove access to any copies of the Springer Nature journal content which have been saved.To the fullest extent permitted by law, Springer Nature makes no warranties, representations or guarantees to Users, either express or implied with respect to the Springer nature journal content and all parties disclaim and waive any implied warranties or warranties imposed by law, including merchantability or fitness for any particular purpose.Please note that these rights do not automatically extend to content, data or other material published by Springer Nature that may be licensed from third parties.If you would like to use or distribute our Springer Nature journal content to a wider audience or on a regular basis or in any other manner not expressly permitted by these Terms, please contact Springer Nature at onlineservice@springernature.com

Fig. 2
Fig. 2 Left: The grid nodes for H x projected onto the (x, y) plane in the case N x = N y = 5.The unit square is bounded by the dashed line.Solid circles show the mesh points of H x h .Hollow circles denote the ghost points induced by the Neumann boundary conditions at y = 0, 1. Right: The grid nodes for H y projected onto the (x, y) plane in the case N x = N y = 5.The unit square is bounded by the dashed frame.Solid stars show the mesh points of H y h .Hollow stars denote ghost points induced by the Neumann boundary conditions at x = 0, 1

Table 1
Second-order derivatives in time and mixed space-time derivatives at the boundaries of

Table 2
Boundary conditions for the time derivatives of the Cartesian components of E and H With these arrangements, the derivatives ∂ t E and ∂ t H are supplied with well-defined boundary conditions, as summarized in Table 2.In short, on ∂ we have ∂ t E = ∂ t H ⊥ = 0, while ∂ t E ⊥ and ∂ t H obey a Neumann condition.As discussed earlier, Dirichlet conditions are true boundary conditions while Neumann conditions are derived from the interior equations.

Table 4
Grid convergence for the three schemes we have chosen.N is the grid dimension

Table 5
The