Energy Stable Nodal Discontinuous Galerkin Methods for Nonlinear Maxwell’s Equations in Multi-dimensions

Abstract

In this work, we extend the energy stable discontinuous Galerkin (DG) schemes proposed in Bokil et al. (J Comput Phys 350:420–452, 2017), for the time domain Maxwell’s equations augmented with a class of nonlinear constitutive polarization laws, to higher dimensions. The nontrivial discrete temporal treatment of the nonlinearity in the ordinary differential equations that encode the Kerr and Raman effects (Bokil et al. 2017), is first generalized to higher spatial dimensions. To further improve the computational efficiency in dealing with the nonlinearity, we apply nodal DG methods in space. Energy stability is proved for the semi-discrete in time and in space schemes as well as for the fully-discrete schemes. Under some assumptions on the strength of nonlinearity, error estimates are established for the semi-discrete in space methods, and, in particular, optimal accuracy is achieved for the methods on Cartesian meshes with $Q^k$-type elements and alternating fluxes. Attention is paid to the role of the nodal form of the DG discretizations in the analysis. We numerically validate the accuracy, energy stability, and computational efficiency of the proposed schemes using manufactured solutions. We further illustrate the performance of the methods through physically relevant experiments involving spatial soliton propagation and airhole scattering in realistic glasses.

Appendices

Appendix A: Proof of Theorem 3

This section will be devoted to the proof of Theorem 3. We will start with some preparatory results on projection operators and their approximation properties, as well as inverse inequalities.

1.1 A.1 Preliminaries

We first consider the commonly used $L^2$ projection $\pi _{xh}$ and Gauss-Radau projections $\pi _{xh}^{\pm }$ onto $U_{xh}^k = \left\{ u\in L^2([x_a, x_b]): u|_{I_{i}}\in P^k(I_i), \;\forall i\right\}$ in the x direction on $[x_a, x_b]$. Here $P^k(I)$ consists of polynomials of degree up to k on an interval I.

• 1. $L^2$ projection $\pi _{xh}$: $L^2([x_a,x_b])\rightarrow U_{xh}^k$, satisfying, $\forall i$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}(\pi _{xh}w)(x)v(x) dx = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}w(x)v(x) dx, \quad \forall v\in P^{k}(I_i); \end{aligned}$$

• 2. Gauss-Radau projection $\pi _{xh}^{+}$: $H^1([x_a,x_b])\rightarrow U_{xh}^k$, satisfying, $\forall i$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}(\pi _{xh}^{+}w)(x)v(x) dx = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}w(x)v(x)dx, \quad \forall v\in P^{k-1}(I_i), \end{aligned}$$

  and

  $$\begin{aligned} \Big (\pi _{xh}^{+}w\Big )(x_{i-\frac{1}{2}}^{+}) = w(x_{i-\frac{1}{2}}^{+}); \end{aligned}$$

• 3. Gauss-Radau projection $\pi _{xh}^{-}$: $H^1([x_a,x_b]) \rightarrow U_{xh}^k$, satisfying, $\forall i$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}(\pi _{xh}^{-}w)(x)v(x) dx = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}w(x)v(x) dx, \quad \forall v\in P^{k-1}(I_i), \end{aligned}$$

  and

  $$\begin{aligned} \Big (\pi _{xh}^{-}w\Big )(x_{i+\frac{1}{2}}^{-}) = w(x_{i+\frac{1}{2}}^{-}). \end{aligned}$$

The one-dimensional $L^2$ projection $\pi _{yh}$ and Gauss-Radau projections $\pi _{yh}^{\pm }$ can be defined similarly in the y direction onto $U_{yh}^k = \left\{ u\in L^2([y_a, y_b]): u|_{J_{j}}\in P^k(J_j), \;\forall j\right\}$. One can further define two-dimensional projection operators as tensor products of one-dimensional ones as follows [9, 20].

• 1. $\Pi _{h}^{\pm ,0}=\pi _{xh}^{\pm }\otimes \pi _{yh}$: $H^{2}(\Omega )\rightarrow V_h^k$, satisfying, $\forall v\in Q^{k}(K_{ij})$ and $\forall i,j$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left( \Pi _{h}^{\pm ,0}w\right) (x,y) \frac{\partial v(x,y)}{\partial x} dxdy = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x,y) \frac{\partial v(x,y)}{\partial x} dxdy \end{aligned}$$

  and

  $$\begin{aligned} \int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\Big (\Pi _{h}^{\pm ,0}w\Big )(x_{i\mp \frac{1}{2}}^{\pm },y)v(x_{i\mp \frac{1}{2}}^{\pm },y)dy = \int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x_{i\mp \frac{1}{2}}^{\pm },y)v(x_{i\mp \frac{1}{2}}^{\pm },y)dy; \end{aligned}$$

• 2. $\Pi _{h}^{0,\pm }=\pi _{xh}\otimes \pi _{yh}^{\pm }$: $H^{2}(\Omega )\rightarrow V_h^k$, satisfying, $\forall v\in Q^{k}(K_{ij})$ and $\forall i,j$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left( \Pi _{h}^{0,\pm }w\right) (x,y) \frac{\partial v(x,y)}{\partial y} dxdy = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x,y) \frac{\partial v(x,y)}{\partial y} dxdy \end{aligned}$$

  and

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Big (\Pi _{h}^{0,\pm }w\Big )(x,y_{j\mp \frac{1}{2}}^{\pm })v(x,y_{j\mp \frac{1}{2}}^{\pm })dx = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}w(x,y_{j\mp \frac{1}{2}}^{\pm })v(x,y_{j\mp \frac{1}{2}}^{\pm })dx; \end{aligned}$$

• 3. $\Pi _{h}^{\pm ,\pm }=\pi _{xh}^{\pm }\otimes \pi _{yh}^{\pm }$: $H^{2}(\Omega )\rightarrow V_h^k$, satisfying, $\forall v\in Q^{k-1}(K_{ij})$ and $\forall i,j$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left( \Pi _{h}^{\pm ,\pm }w\right) (x,y) v(x,y) dxdy= & {} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x,y) v(x,y) dxdy,\\ \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\Big (\Pi _{h}^{\pm ,\pm }w\Big )(x,y_{j\mp \frac{1}{2}}^{\pm })v(x,y_{j\mp \frac{1}{2}}^{\pm })dx= & {} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}w(x,y_{j\mp \frac{1}{2}}^{\pm })v(x,y_{j\mp \frac{1}{2}}^{\pm })dx,\\ \int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\Big (\Pi _{h}^{\pm ,\pm }w\Big )(x_{i\mp \frac{1}{2}}^{\pm },y)v(x_{i\mp \frac{1}{2}}^{\pm },y)dy= & {} \int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x_{i\mp \frac{1}{2}}^{\pm },y)v(x_{i\mp \frac{1}{2}}^{\pm },y)dy \end{aligned}$$

  and

  $$\begin{aligned} \Big (\Pi _{h}^{\pm ,\pm }w\Big )(x_{i\mp \frac{1}{2}}^{\pm },y_{j\mp \frac{1}{2}}^{\pm }) = w(x_{i\mp \frac{1}{2}}^{\pm },y_{j\mp \frac{1}{2}}^{\pm }); \end{aligned}$$

• 4.) $\Pi _{h}^{0,0}=\pi _{xh}\otimes \pi _{yh}$: $L^2(\Omega )\rightarrow V_h^k$, satisfying, $\forall v\in Q^{k}(K_{ij})$ and $\forall i,j$,

  $$\begin{aligned} \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\left( \Pi _{h}^{0,0}w\right) (x,y) v(x,y) dxdy = \int _{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\int _{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}w(x,y) v(x,y) dxdy. \end{aligned}$$

The following approximation results are standard for the analysis of DG methods, and can be established following classical arguments [8].

Lemma A1

(Approximation properties) Let ${\mathcal {P}}_h$ be either the interpolation operator ${\mathcal {I}}_h$, or any of the projection operators $\Pi _{h}^{\pm , 0}, \Pi _{h}^{0,\pm }, \Pi _{h}^{\pm ,\pm }, \Pi _{h}^{0,0}$. There exist constants $C_\star$ and $C_k$, such that $\forall w\in H^{k+1}(\Omega )$, there hold

$$\begin{aligned} \Vert w-{\mathcal {P}}_hw\Vert _{L^2(K)}^2 + h \Vert w-{\mathcal {P}}_hw\Vert _{L^2(\partial K)}^2 \le C_\star h^{2k+2}\Vert w\Vert _{H^{k+1}(K)}^2,~\forall K\in {\mathcal {T}}_h \end{aligned}$$

(A1)

and

$$\begin{aligned} \Vert {\mathcal {P}}_hw\Vert _\infty \le C_k\Vert w\Vert _\infty . \end{aligned}$$

(A2)

As a direct consequence of (A2), there holds

$$\begin{aligned} \Vert w-{\mathcal {P}}_hw\Vert _\infty \le C_k \Vert w\Vert _\infty . \end{aligned}$$

(A3)

Associated with terms in our methods involving the interpolation operator to deal with nonlinearity more efficiently, it is convenient to work with a discrete norm $\Vert \cdot \Vert _h$ on $V_h^k$, which is equivalent to the standard $L^2$ norm.

Lemma A2

(Norm equivalence [10]) Define

$$\begin{aligned} \Vert v\Vert _h = \left( \sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\sum _{m,n=0}^{k}\left| v(x_{im},y_{jn})\right| ^2\Delta x_i \Delta y_j\right) ^{1/2}, \end{aligned}$$

then $\Vert \cdot \Vert _h$ is a norm on $V_h^k$. In addition, there exist constants $C_\star , \widetilde{C_\star }>0$, such that

$$\begin{aligned} \widetilde{C_\star }\Vert v\Vert _h\le \Vert v\Vert \le C_\star \Vert v\Vert _h, \quad \forall v\in V_h^k. \end{aligned}$$

Below are some immediate results of Lemma A2.

Lemma A3

There exists a constant $C_\star$, such that

$$\begin{aligned} \left| \int _{\Omega }{\mathcal {I}}_h (fg) d\Omega \right| \le C_\star \Vert f\Vert \Vert g\Vert , \quad \forall f,g\in W_h(\Omega ). \end{aligned}$$

(A4)

Proof

By the definition of ${\mathcal {I}}_h$, and the exactness of $(k+1)$-point Gaussian-Legendre quadrature rule, the boundedness of $\{{\widehat{\omega }}_m\}_{m=0}^k$, Cauchy-Schwartz inequality, and norm equivalence in Lemma A2, we have

$$\begin{aligned} \left| \int _{\Omega }{\mathcal {I}}_h (fg) d\Omega \right|&=\left| \sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\frac{\Delta x_i}{2}\frac{\Delta y_j}{2}\sum _{m,n=0}^{k}{\widehat{\omega }}_m {\widehat{\omega }}_n f(x_{im},y_{jn})g(x_{im},y_{jn}) \right| \\&\le C_k \Vert f\Vert _h \Vert g\Vert _h \le C_\star \Vert f\Vert \Vert g\Vert . \end{aligned}$$

$\square$

Finally, our analysis will need some inverse inequality [8].

Lemma A4

(Inverse inequality) There exists a constant $C_\star$, such that

$$\begin{aligned} h^{2}\Vert \nabla v\Vert _{L^2(K)}^2 + h \Vert v\Vert _{L^2(\partial K)}^2 \le C_\star \Vert v\Vert _{L^{2}(K)}^2,\quad \forall v\in V_h^{k}, \; \forall K\in {\mathcal {T}}_h. \end{aligned}$$

(A5)

1.2 A.2 Proof of Theorem 3

Proof

The proof will proceed similarly as that in [3]. The DG methods in [3] are of the modal form and for the 1D version of model (1). One difference in our analysis here comes from the use of the interpolation operator ${\mathcal {I}}_h$ in the nodal DG setting to deal with nonlinear terms. Being two-dimensional also adds complexity to the use of projection operators in the proof.

For any component u of the solution (such as $H_z, E_x, E_y$ etc), we first decompose its error as $u-u_h = (u-{\mathcal {P}}^u_{h}u) - (u_h-{\mathcal {P}}^u_{h}u)=\eta _u -\xi _u$, where the projection ${\mathcal {P}}^u_h$ is taken to be $\Pi _h^{0,0}$ except for ${\mathcal {P}}_h^{E_x}$, ${\mathcal {P}}_h^{E_y}$, ${\mathcal {P}}_h^{H_z}$ when the numerical fluxes are alternating. When alternating fluxes are used, without loss of generality, we only consider one of the four possible cases, particularly, we consider Alternating I, given by (28) with $\dagger =+, \ddagger =-, \natural =+, \sharp =-$, and take

$$\begin{aligned} {\mathcal {P}}_h^{E_x} = \Pi _{h}^{0,+},~{\mathcal {P}}_h^{E_y} = \Pi _{h}^{+,0},~{\mathcal {P}}_h^{H_z} = \Pi _{h}^{-,-}. \end{aligned}$$

We initialize the method by setting $\xi _u=u_h-{\mathcal {P}}^u_{h}u=0$ at $t=0$ for all solution components.

${\mathbf{Step \,\,\mathbf{1}: }}$ With the consistency of numerical fluxes, from (25), we can obtain the error equations:

$$\begin{aligned}&\mu _0(\partial _{t}\eta _{H_{z}},\phi ) + {\mathcal {B}}_h^{\varvec{E}}(\eta _{E_x},\eta _{E_y},\phi ) =\mu _0(\partial _{t}\xi _{H_{z}},\phi ) + {\mathcal {B}}_h^{\varvec{E}}(\xi _{E_x},\xi _{E_y},\phi ),\quad \forall \phi \in V_h^k, \end{aligned}$$

(A6a)

$$\begin{aligned}&\Big (\partial _{t}(D_{x}-D_{xh}),\phi \Big ) + {\mathcal {B}}_{xh}^{H}(\eta _{H_z},\phi ) = {\mathcal {B}}_{xh}^{H}(\xi _{H_z},\phi ),\quad \forall \phi \in V_h^k, \end{aligned}$$

(A6b)

$$\begin{aligned}&\Big (\partial _{t}(D_{y}-D_{yh}),\phi \Big ) + {\mathcal {B}}_{yh}^{H}(\eta _{H_z},\phi ) = {\mathcal {B}}_{yh}^{H}(\xi _{H_z},\phi ),\quad \forall \phi \in V_h^k, \end{aligned}$$

(A6c)

$$\begin{aligned}&\partial _{t}\varvec{\eta }_{P}- \varvec{\eta }_{J}=\partial _{t}\varvec{\xi }_{P}- \varvec{\xi }_{J}, \end{aligned}$$

(A6d)

$$\begin{aligned}&\partial _{t}\varvec{\eta }_{J} + \gamma \varvec{\eta }_{J} + \omega _{0}^2\varvec{\eta }_{P} - \omega _{p}^2 \varvec{\eta }_{E} = \partial _{t}\varvec{\xi }_{J} + \gamma \varvec{\xi }_{J} + \omega _{0}^2\varvec{\xi }_{P} - \omega _{p}^2 \varvec{\xi }_{E}, \end{aligned}$$

(A6e)

$$\begin{aligned}&\partial _{t}\eta _{Q}-\eta _{\sigma } = \partial _{t}\xi _{Q}-\xi _{\sigma }, \nonumber \\&\partial _{t}\eta _{\sigma } + \gamma _v\eta _{\sigma } + \omega _v^2\eta _{Q} - \partial _{t}\xi _{\sigma } - \gamma _v\xi _{\sigma } - \omega _v^2\xi _{Q} = \omega _v^2\left( |\varvec{E}|^2-{\mathcal {I}}_h\left( |\varvec{E}_h|^2\right) \right) \end{aligned}$$

(A6f)

$$\begin{aligned}&\quad = \omega _v^2\left( |\varvec{E}|^2-{\mathcal {I}}_h\left( |\varvec{E}|^2\right) \right) + \omega _v^2{\mathcal {I}}_h\left( |\varvec{E}|^2 - |\varvec{E}_h|^2\right) , \end{aligned}$$

(A6g)

along with

$$\begin{aligned} \varvec{D}-\varvec{D}_{h} =&\,\, \epsilon _0\epsilon _{\infty }(\varvec{\eta _{E}}-\varvec{\xi _{E}}) + \epsilon _0(\varvec{\eta _{P}}-\varvec{\xi _{P}}) + \epsilon _{0}a(1-\theta )\left( |\varvec{E}|^2\varvec{E}-{\mathcal {I}}_h\left( |\varvec{E}_{h}|^2\varvec{E_{h}}\right) \right) \nonumber \\&+ \epsilon _{0}a\theta (Q\varvec{E}-{\mathcal {I}}_h(Q_{h}\varvec{E}_{h})). \end{aligned}$$

(A7)

Note that ${\mathcal {P}}_h^{\varvec{P}}={\mathcal {P}}_h^{\varvec{J}}=\Pi _h^{0,0}$, then the equation (1d) implies

$$\begin{aligned} \partial _{t}\varvec{\eta }_{P}- \varvec{\eta }_{J}=\varvec{0}. \end{aligned}$$

(A8)

Therefore (A6d) further gives

$$\begin{aligned} \partial _{t}\varvec{\xi }_{P}- \varvec{\xi }_{J}=\varvec{0}. \end{aligned}$$

(A9)

Similarly, (A6f) gives

$$\begin{aligned} \partial _{t}\xi _{Q}-\xi _{\sigma }=0. \end{aligned}$$

(A10)

Recall that the choice of ${\mathcal {P}}_h^{\varvec{E}}$ depends on the numerical fluxes, then the equation (1e) implies $\partial _{t}\varvec{\eta }_{J} + \gamma \varvec{\eta }_{J} + \omega _{0}^2\varvec{\eta }_{P} - \omega _{p}^2 (\varvec{E} - \Pi _h^{0,0} \varvec{E})=0,$ this together with (A6e) leads to

$$\begin{aligned} \partial _{t}\varvec{\xi }_{J} + \gamma \varvec{\xi }_{J} + \omega _{0}^2\varvec{\xi }_{P} - \omega _{p}^2 \varvec{\xi }_{E}=\omega _{p}^2 \varvec{\mho _E}, \end{aligned}$$

(A11)

where

$$\begin{aligned} \varvec{\mho _E}=(\varvec{E} - \Pi _h^{0,0} \varvec{E})-\varvec{\eta }_{E}=\left\{ \begin{array}{ll} \varvec{0}&{}\text{ central } \text{ flux } \text{ for } \varvec{E}\\ \Pi _h^{\varvec{E}} \varvec{E} - \Pi _h^{0,0} \varvec{E} &{}\text{ alternating } \text{ flux } \text{ for } \varvec{E} \end{array} \right. . \end{aligned}$$

(A12)

${\mathbf{Step \,\,\mathbf{2}: }}$ Take $\phi = \xi _{H_z}$ in (A6a), $\phi = \xi _{E_x}$ in (A6b), and $\phi = \xi _{E_y}$ in (A6c), sum them up and use (30), we have,

$$\begin{aligned}&\frac{\mu _0}{2}\frac{d}{dt}\Vert \xi _{H_z}\Vert ^2 - \mu _0(\partial _{t}\eta _{H_{z}},\xi _{H_z}) - {\mathcal {B}}_h^{\varvec{E}}(\eta _{E_x},\eta _{E_y},\xi _{H_z}) - {\mathcal {B}}_{xh}^{H}(\eta _{H_z},\xi _{E_x}) - {\mathcal {B}}_{yh}^{H}(\eta _{H_z},\xi _{E_y}) \nonumber \\&\quad =\Big (\partial _{t}(\varvec{D}-\varvec{D}_{h}),\varvec{\xi _{E}}\Big ). \end{aligned}$$

(A13)

Differentiate (A7) with respect to time t, and use that $\partial _{t}$ and ${\mathcal {I}}_h$ commute, we have

$$\begin{aligned} \Big (\partial _{t}(\varvec{D}-\varvec{D}_{h}),\varvec{\xi _{E}}\Big )&= \epsilon _0\epsilon _{\infty }\Big ( \partial _{t}(\varvec{\eta _{E}}-\varvec{\xi _{E}}),\varvec{\xi _{E}}\Big ) + \epsilon _0\Big (\partial _{t}(\varvec{\eta _{P}}-\varvec{\xi _{P}}),\varvec{\xi _{E}}\Big )\nonumber \\&\quad \quad + \epsilon _{0}a(1-\theta )\Big (\partial _{t}(|\varvec{E}|^2\varvec{E})-{\mathcal {I}}_h\left( \partial _{t}(|\varvec{E}|^2\varvec{E})\right) ,\varvec{\xi _{E}} \Big ) \nonumber \\&\quad \quad + \epsilon _{0}a(1-\theta )\Big ({\mathcal {I}}_h\partial _{t}\Big (|\varvec{E}|^2\varvec{E}-|\varvec{E}_{h}|^2\varvec{E_{h}}\Big ),\varvec{\xi _{E}}\Big )\nonumber \\&\quad \quad + \epsilon _{0}a\theta \Big (\partial _{t}(Q\varvec{E})-{\mathcal {I}}_h\left( \partial _{t}(Q\varvec{E})\right) ,\varvec{\xi _{E}}\Big ) + \epsilon _{0}a\theta \Big ({\mathcal {I}}_h\partial _{t}\Big (Q\varvec{E}-Q_{h}\varvec{E}_{h}\Big ),\varvec{\xi _{E}}\Big ).\nonumber \\ \end{aligned}$$

(A14)

In next few steps, we will work with the terms on the right side of (A14) based on the error equations.

$\mathbf{Step \,\,\mathbf{2}.1: \,\,\mathbf{the} \,\,\mathbf{fourth} \,\,\mathbf{term} \,\,\mathbf{in} \,\, (A14)}$. Using $\varvec{E}_h=\varvec{E}-\varvec{\eta _{E}}+\varvec{\xi _{E}}$ and some direct manipulation, one can get

$$\begin{aligned}&\partial _{t}\Big (|\varvec{E}|^2\varvec{E} - |\varvec{E}_{h}|^2\varvec{E}_{h}\Big )\varvec{\cdot }\varvec{\xi _{E}}\nonumber \\&\quad =\partial _{t}\Big (|\varvec{E}|^2\varvec{E}\Big )\varvec{\cdot }\varvec{\xi _{E}} - 2(\varvec{E}_{h}\varvec{\cdot }\partial _{t}\varvec{E}_h)(\varvec{E}_h\varvec{\cdot }\varvec{\xi _{E}}) - |\varvec{E}_h|^2\partial _{t}\varvec{E}_h\varvec{\cdot }\varvec{\xi _{E}}\nonumber \\&\quad =- \frac{3}{4}\partial _{t}\Big (|\varvec{\xi _{E}}|^4\Big ) - 2\partial _{t}\Big (|\varvec{\xi _{E}}|^2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big ) \nonumber \\&\qquad - \partial _{t}\Big (((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})^2\Big )- \frac{1}{2}\partial _{t}\Big (|\varvec{E}-\varvec{\eta _{E}}|^2|\varvec{\xi _{E}|^2}\Big )\nonumber \\&\qquad - \frac{1}{2}\partial _{t}\Big (|\varvec{E}-\varvec{\eta _{E}})|^2\Big )|\varvec{\xi _{E}}|^2 - 2\Big ((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big )\Big (\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big ) \nonumber \\&\qquad - |\varvec{\xi _{E}}|^2\Big (\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big )\nonumber \\&\qquad - |\varvec{E}-\varvec{\eta _{E}}|^2\Big (\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big ) - \partial _{t}\Big (|\varvec{E}-\varvec{\eta _{E}}|^2\Big )\Big ((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big )\nonumber \\&\qquad + \partial _{t}\Big (|\varvec{E}|^2\varvec{E}\Big )\varvec{\cdot }\varvec{\xi _{E}}. \end{aligned}$$

(A15)

And the last three terms in (A15) can be further written as

$$\begin{aligned}&-\partial _{t}\Big (|\varvec{E}-\varvec{\eta _{E}}|^2(\varvec{E}-\varvec{\eta _{E}}) - |\varvec{E}|^2\varvec{E}\Big )\varvec{\cdot }\varvec{\xi _{E}}\nonumber \\&\quad = \partial _{t}\Big (|\varvec{\eta _{E}}|^2\varvec{\eta _{E}} - |\varvec{\eta _{E}}|^2\varvec{E} - 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{\eta _{E}} + 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{E} + |\varvec{E}|^2\varvec{\eta _{E}}\Big )\varvec{\cdot }\varvec{\xi _{E}}. \end{aligned}$$

(A16)

With (A16) and (A15), and Lemma 1, we now have

$$\begin{aligned}&\Big ({\mathcal {I}}_h\partial _{t}(|\varvec{E}|^2\varvec{E}-|\varvec{E}_{h}|^2\varvec{E_{h}}),\varvec{\xi _{E}}\Big )=\int _{\Omega }{\mathcal {I}}_h\Big (\partial _{t}(|\varvec{E}|^2\varvec{E} - |\varvec{E}_{h}|^2\varvec{E}_{h})\varvec{\cdot }\varvec{\xi _{E}}\Big )d\Omega \nonumber \\&\quad =- \frac{3}{4}\frac{d}{dt}\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^4\right) d\Omega - 2\frac{d}{dt}\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega \nonumber \\&\qquad - \frac{d}{dt}\int _{\Omega } {\mathcal {I}}_h((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})^2 d\Omega \nonumber \\&\qquad - \frac{1}{2}\frac{d}{dt}\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{E}-\varvec{\eta _{E}}|^2|\varvec{\xi _{E}}|^2\right) d\Omega - \frac{1}{2}\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(|\varvec{E}-\varvec{\eta _{E}})|^2)|\varvec{\xi _{E}}|^2\right) d\Omega \nonumber \\&\qquad - 2\int _{\Omega }{\mathcal {I}}_h\left( ((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})(\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})\right) d\Omega \nonumber \\&\qquad - \int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2(\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})\right) d\Omega \nonumber \\&\qquad +\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(|\varvec{\eta _{E}}|^2\varvec{\eta _{E}} - |\varvec{\eta _{E}}|^2\varvec{E} - 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{\eta _{E}} \right. \nonumber \\&\qquad \left. +\, 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{E} + |\varvec{E}|^2\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega . \end{aligned}$$

(A17)

$\mathbf{Step \,\,\mathbf{2}.2: \,\,\mathbf{the} \,\,\mathbf{sixth} \,\,\mathbf{term} \,\,\mathbf{in} \,\, (A14)}$. We consider

$$\begin{aligned} \partial _{t}(Q\varvec{E}-Q_{h}\varvec{E}_{h})\varvec{\cdot }\varvec{\xi _{E}}&=\partial _{t}\Big (Q\varvec{E} - (Q-\eta _Q+\xi _Q)(\varvec{E}-\varvec{\eta _{E}}+\varvec{\xi _{E}})\Big )\varvec{\cdot }\varvec{\xi _{E}}\nonumber \\&=\partial _{t}\Big (Q\varvec{\eta _{E}} + \eta _Q(\varvec{E}-\varvec{\eta _{E}})\Big )\varvec{\cdot }\varvec{\xi _{E}} - \frac{1}{2}\partial _{t}(Q-\eta _Q)|\varvec{\xi _{E}}|^2 \nonumber \\&\quad - \frac{1}{2}\partial _{t}\Big ((Q-\eta _Q)|\varvec{\xi _{E}}|^2\Big ) - \partial _{t}\xi _Q(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\nonumber \\&\quad - \xi _Q\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}} - \frac{1}{2}\partial _{t}\xi _Q|\varvec{\xi _{E}}|^2 - \frac{1}{2}\partial _{t}\Big (\xi _Q|\varvec{\xi _{E}}|^2\Big ). \end{aligned}$$

(A18)

Observe that

$$\begin{aligned} \Big (|\varvec{E}|^2-|\varvec{E}_h|^2\Big )\xi _{\sigma } =&\Big (|\varvec{E}|^2-|\varvec{E}-\varvec{\eta _{E}}+\varvec{\xi _{E}}|^2\Big )\xi _{\sigma }\nonumber \\ =&\Big (|\varvec{E}|^2-|\varvec{E}-\varvec{\eta _{E}}|^2 - 2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}} - |\varvec{\xi _{E}}|^2\Big )\xi _{\sigma }\nonumber \\ =&\Big (2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}} - |\varvec{\eta _{E}}|^2 - 2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\Big )\xi _{\sigma } - \partial _{t}\xi _{Q}|\varvec{\xi _{E}}|^2, \end{aligned}$$

(A19)

where we have used (A10) in the last equality. With (A19), (A10) and (A6g), and Lemma 1, we have

$$\begin{aligned}&\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}\xi _{Q}|\varvec{\xi _{E}}|^2\right) d\Omega \nonumber \\&\quad = - \int _{\Omega }{\mathcal {I}}_h\left( |\varvec{E}|^2-|\varvec{E}_h|^2\right) \xi _{\sigma }d\Omega \nonumber \\&\qquad + \int _{\Omega }{\mathcal {I}}_h\Big ((2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2-2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})\xi _{\sigma }\Big )d\Omega \nonumber \\&\quad = \int _{\Omega }\left( |\varvec{E}|^2-{\mathcal {I}}_h\left( |\varvec{E}|^2\right) \right) \xi _{\sigma } d\Omega \nonumber \\&\qquad + \int _{\Omega }{\mathcal {I}}_h\Big ((2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2-2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})\xi _{\sigma }\Big )d\Omega \nonumber \\&\qquad - \frac{1}{\omega _v^2}\int _{\Omega }\left( \partial _{t}\eta _{\sigma }+\gamma _v\eta _{\sigma } + \omega _v^2\eta _Q\right) \xi _{\sigma }d\Omega \nonumber \\&\qquad + \frac{1}{2\omega _v^2}\frac{d}{dt}\int _{\Omega }\xi _{\sigma }^2d\Omega + \frac{\gamma _v}{\omega _v^2}\int _{\Omega }\xi _{\sigma }^2d\Omega + \frac{1}{2}\frac{d}{dt}\int _{\Omega }\xi _{Q}^2d\Omega . \end{aligned}$$

(A20)

Combine (A18) and (A20), and use Lemma 1, we have

$$\begin{aligned}&\Big ({\mathcal {I}}_h\partial _{t}\Big (Q\varvec{E}-Q_{h}\varvec{E}_{h}\Big ),\varvec{\xi _{E}}\Big )=\int _{\Omega }{\mathcal {I}}_h\Big (\partial _{t}(Q\varvec{E}-Q_{h}\varvec{E}_{h})\varvec{\cdot }\varvec{\xi _{E}}\Big )d\Omega \nonumber \\&\quad = - \frac{d}{dt}\int _{\Omega }\left( \frac{1}{2}{\mathcal {I}}_h\left( (Q-\eta _Q)|\varvec{\xi _{E}}|^2\right) + \frac{1}{2}{\mathcal {I}}_h\left( \xi _Q|\varvec{\xi _{E}}|^2\right) + \frac{1}{4\omega _v^2}\xi _{\sigma }^2 + \frac{1}{4}\xi _{Q}^2\right) d\Omega \nonumber \\&\qquad - \frac{\gamma _v}{2\omega _v^2}\int _{\Omega }\xi _{\sigma }^2d\Omega \nonumber \\&\qquad + \int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(Q\varvec{\eta _{E}} + \eta _Q(\varvec{E}-\varvec{\eta _{E}}))\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega - \frac{1}{2}\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(Q-\eta _Q)|\varvec{\xi _{E}}|^2\right) d\Omega \nonumber \\&\qquad - \int _{\Omega }{\mathcal {I}}_h\left( \xi _Q\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega - \frac{1}{2}\int _{\Omega } {\mathcal {I}}_h\Big ((2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2)\xi _{\sigma }\Big ) d\Omega \nonumber \\&\qquad + \frac{1}{2\omega _v^2}\int _{\Omega }(\partial _{t}\eta _{\sigma }+\gamma _v\eta _{\sigma }+\omega _v^2\eta _Q)\xi _{\sigma }d\Omega - \frac{1}{2}\int _{\Omega }\left( |\varvec{E}|^2-{\mathcal {I}}_h\left( |\varvec{E}|^2\right) \right) \xi _{\sigma } d\Omega . \end{aligned}$$

(A21)

Here we have dropped the term $\int _\Omega {\mathcal {I}}_h\Big ((\partial _{t}\xi _Q-\xi _\sigma )\varvec{(E-\eta _E)\cdot \xi _E}\Big )d\Omega$ due to (A10).

$\mathbf{Step \,\,\mathbf{2}.3: \,\,\mathbf{the} \,\,\mathbf{first} \,\,\mathbf{and} \,\,\mathbf{second} \,\,\mathbf{terms} \,\,\mathbf{in} \,\, (A14)}$. Use (A11), (A9), and (A8), we have

$$\begin{aligned}&(\partial _{t}(\varvec{\eta _{P}}-\varvec{\xi _{P}}), \varvec{\xi _{E}})\nonumber \\&\quad = (\partial _{t}\varvec{\eta _{P}}, \varvec{\xi _{E}}) -\frac{1}{\omega _p^2}(\partial _{t}\varvec{\xi _{P}}, \partial _{t}\varvec{\xi _{J}}+\gamma \varvec{\xi _{J}}+\omega _0^2\varvec{\xi _{P}}-\omega _p^2\varvec{\mho _E})\nonumber \\&\quad = (\varvec{\eta _{J}}, \varvec{\xi _{E}}) + (\varvec{\mho _E}, \varvec{\xi _{J}}) - \frac{1}{2\omega _p^2}\frac{d}{dt}\int _\Omega |\varvec{\xi _J}|^2 d\Omega \nonumber \\&\qquad -\frac{\gamma }{\omega _p^2}\int _\Omega |\varvec{\xi _J}|^2 d\Omega - \frac{\omega _0^2}{2\omega _p^2}\frac{d}{dt}\int _\Omega |\varvec{\xi _{P}}|^2 d\Omega . \end{aligned}$$

(A22)

In addition

$$\begin{aligned} (\partial _{t}(\varvec{\eta _{E}}-\varvec{\xi _{E}}), \varvec{\xi _{E}})=(\partial _{t}\varvec{\eta _{E}}, \varvec{\xi _{E}})-\frac{1}{2}\frac{d}{dt}\int _\Omega |\varvec{\xi _E}|^2 d\Omega . \end{aligned}$$

(A23)

Now we can gather (A17), (A21)–(A23) for (A13)–(A14), and come up with the following identity,

$$\begin{aligned} \frac{d}{dt}\widehat{{\mathcal {E}}_h}+ \frac{\epsilon _0\gamma }{\omega _p^2}\int _{\Omega }|\varvec{\xi _{J}}|^2d\Omega +\frac{\epsilon _0a\theta \gamma _v}{2\omega _v^2}\int _{\Omega }\xi _{\sigma }^2d\Omega = \sum _{i=0}^4 \Lambda _i, \end{aligned}$$

(A24)

where

$$\begin{aligned} \widehat{{\mathcal {E}}_h}&= \int _{\Omega }\Big (\frac{\mu _0}{2}\xi _{H_z}^2 + \frac{\epsilon _0\epsilon _{\infty }}{2}|\varvec{\xi _{E}}|^2 + \frac{\epsilon _0}{2\omega _p^2}|\varvec{\xi _{J}}|^2 + \frac{\epsilon _0\omega _0^2}{2\omega _p^2}|\varvec{\xi _{P}}|^2 + \frac{\epsilon _0a\theta }{4\omega _v^2}\xi _{\sigma }^2\Big )d\Omega \\&\quad + \int _{\Omega }{\mathcal {I}}_h\left( \frac{\epsilon _0a\theta }{2}\xi _{Q}|\varvec{\xi _{E}}|^2 + \frac{3\epsilon _0 a(1-\theta )}{4}|\varvec{\xi _E}|^4 + \frac{\epsilon _0 a\theta }{4}\xi _Q^2\right) d\Omega \\&\quad + \frac{\epsilon _0a\theta }{2}\int _{\Omega }{\mathcal {I}}_h\left( (Q-\eta _Q)|\varvec{\xi _{E}}|^2\right) d\Omega \\&\quad + \epsilon _0a(1-\theta )\int _{\Omega }{\mathcal {I}}_h\left( 2|\varvec{\xi _{E}}|^2(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}} + \left( (\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) ^2 + \frac{1}{2}|\varvec{E}-\varvec{\eta _{E}}|^2|\varvec{\xi _{E}}|^2\right) d\Omega ,\\ \Lambda _0&= \epsilon _{0}a(1-\theta )\Big (\partial _{t}(|\varvec{E}|^2\varvec{E})-{\mathcal {I}}_h\left( \partial _{t}(|\varvec{E}|^2\varvec{E})\right) ,\varvec{\xi _{E}} \Big )+ \epsilon _{0}a\theta \Big (\partial _{t}(Q\varvec{E})-{\mathcal {I}}_h\left( \partial _{t}(Q\varvec{E})\right) ,\varvec{\xi _{E}}\Big )\\&\quad - \frac{\epsilon _{0}a\theta }{2}\Big (|\varvec{E}|^2-{\mathcal {I}}_h(|\varvec{E}|^2), \xi _{\sigma }\Big ), \quad \Lambda _1 = \mu _0(\partial _{t}\eta _{H_{z}},\xi _{H_z}),\\ \Lambda _2&= {\mathcal {B}}_h^{\varvec{E}}(\eta _{E_x},\eta _{E_y},\xi _{H_z}) + {\mathcal {B}}_{xh}^{H}(\eta _{H_z},\xi _{E_x}) + {\mathcal {B}}_{yh}^{H}(\eta _{H_z},\xi _{E_y}),\\ \Lambda _3&= \frac{\epsilon _0a\theta }{2\omega _v^2}(\partial _{t}\eta _{\sigma }+\gamma _v\eta _{\sigma } +\omega _v^2\eta _Q,\xi _{\sigma }) - \frac{\epsilon _0a\theta }{2}({\mathcal {I}}_h (2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2),\xi _{\sigma })\\ \Lambda _4&= \epsilon _0\epsilon _\infty (\partial _{t}\varvec{\eta _{E}},\varvec{\xi _{E}}) + \epsilon _0(\varvec{\eta _{J}},\varvec{\xi _{E}}) + \epsilon _0 (\varvec{\mho _E}, \varvec{\xi _{J}}) - \epsilon _0a\theta \int _{\Omega }{\mathcal {I}}_h\left( \xi _Q\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega \\&\quad - 2\epsilon _0a(1-\theta )\int _{\Omega }{\mathcal {I}}_h\Big (((\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})(\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}})\Big )d\Omega \\&\quad - \frac{\epsilon _0a(1-\theta )}{2}\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}\left( |\varvec{E}-\varvec{\eta _{E}}|^2\right) |\varvec{\xi _{E}}|^2\right) d\Omega \\&\quad - \epsilon _0a(1-\theta )\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega \\&\quad + \epsilon _0 a\theta \int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(Q\varvec{\eta _{E}}+\eta _Q(\varvec{E}-\varvec{\eta _{E}}))\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega \\&\quad - \frac{\epsilon _0a\theta }{2}\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(Q-\eta _Q)|\varvec{\xi _{E}}|^2\right) d\Omega \\&\quad + \epsilon _0a(1-\theta )\int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(|\varvec{\eta _{E}}|^2\varvec{\eta _{E}} - |\varvec{\eta _{E}}|^2\varvec{E} - 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{\eta _{E}} \right. \\&\quad \left. + 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{E} + |\varvec{E}|^2\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega . \end{aligned}$$

$\mathbf{Step \,\,\mathbf{3}: }$ Next we will estimate terms on both sides of (A24).

$\mathbf{Step \,\,\mathbf{3}.1: \,\,\mathbf{to} \,\,\mathbf{estimate} \,\,\mathbf{the} \,\,\mathbf{left} \,\,\mathbf{hand} \,\,\mathbf{side} of \,\, (A24)}$. Let $\rho _{err}, \kappa _{err}\in (0,1)$ be two arbitrary constants, then by Cauchy-Schwartz inequality, Lemmas 1-2, we have

$$\begin{aligned} \widehat{{\mathcal {E}}_h} \ge&\int _{\Omega }\left( \frac{\mu _0}{2}\xi _{H_z}^2 + \frac{\epsilon _0\epsilon _{\infty }\kappa _{err}}{2}|\varvec{\xi _{E}}|^2 + \frac{\epsilon _0}{2\omega _p^2}|\varvec{\xi _{J}}|^2 + \frac{\epsilon _0\omega _0^2}{2\omega _p^2}|\varvec{\xi _{P}}|^2 + \frac{\epsilon _0a\theta }{4\omega _v^2}\xi _{\sigma }^2\right) d\Omega \nonumber \\&+ \int _{\Omega }\left( \frac{\epsilon _0 a(1-\theta )\rho _{err}}{12}{\mathcal {I}}_h\left( |\varvec{\xi _E}|^4\right) + \frac{\epsilon _0 a\theta \rho _{err}}{4}\xi _Q^2\right) d\Omega \nonumber \\&+ \int _{\Omega }{\mathcal {I}}_h\left( \frac{\epsilon _0a\theta }{2}\xi _{Q}|\varvec{\xi _{E}}|^2 + \frac{\epsilon _0 a(1-\theta )(1-\rho _{err})}{12}|\varvec{\xi _E}|^4 + \frac{\epsilon _0 a\theta (1-\rho _{err})}{4}\xi _Q^2\right) d\Omega \nonumber \\&+ \frac{\epsilon _0}{2}\int _{\Omega }{\mathcal {I}}_h\Big (\left( \epsilon _{\infty }(1-\kappa _{err})+a\theta (Q-\eta _Q)\right) |\varvec{\xi _{E}}|^2\Big ) d\Omega \nonumber \\&+ \frac{3\epsilon _0a(1-\theta )}{2}\int _{\Omega }{\mathcal {I}}_h\left( \left( (\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}+\frac{2}{3}|\varvec{\xi _{E}}|\right) ^2\right) d\Omega \ge {\mathcal {E}}_h^{(mod)}, \end{aligned}$$

(A25)

where

$$\begin{aligned} {\mathcal {E}}_h^{(mod)}=&\int _{\Omega }\left( \frac{\mu _0}{2}\xi _{H_z}^2 + \frac{\epsilon _0\epsilon _{\infty }\kappa _{err}}{2}|\varvec{\xi _{E}}|^2 + \frac{\epsilon _0}{2\omega _p^2}|\varvec{\xi _{J}}|^2 + \frac{\epsilon _0\omega _0^2}{2\omega _p^2}|\varvec{\xi _{P}}|^2\right) d\Omega \\&+ \int _{\Omega }\left( \frac{\epsilon _0a\theta }{4\omega _v^2}\xi _{\sigma }^2 + \frac{\epsilon _0 a(1-\theta )\rho _{err}}{12}{\mathcal {I}}_h\left( |\varvec{\xi _E}|^4\right) + \frac{\epsilon _0 a\theta \rho _{err}}{4}\xi _Q^2\right) d\Omega . \end{aligned}$$

Here (24) in Lemma 2, as well as Conditions 1-2 in (39)–(40) have been used. In fact, by Condition 1, we have

$$\begin{aligned} \frac{\epsilon _0a\theta }{2}\xi _{Q}|\varvec{\xi _{E}}|^2 + \frac{\epsilon _0 a(1-\theta )(1-\rho _{err})}{12}|\varvec{\xi _E}|^4 + \frac{\epsilon _0 a\theta (1-\rho _{err})}{4}\xi _Q^2\ge 0, \end{aligned}$$

while with Condition 2, there holds

$$\begin{aligned} \epsilon _{\infty }(1-\kappa _{err})+a\theta (Q-\eta _Q)\ge & {} \epsilon _{\infty }(1-\kappa _{err})-a\theta \Vert \Pi _{h}^{0,0}Q\Vert _{\infty }\\\ge & {} \epsilon _{\infty }(1-\kappa _{err})-a\theta C_{k}\Vert Q\Vert _{\infty }\ge 0. \end{aligned}$$

$\mathbf{Step \,\,\mathbf{3}.2: \,\,\mathbf{to} \,\,\mathbf{estimate} \,\, \Lambda _0 \,\,\mathbf{and} \,\, \Lambda _1 \,\,\mathbf{in} \,\,(A24)}$. For $\Lambda _0$, we apply Cauchy-Schwartz inequality and the approximation properties in Lemma A1, and obtain

$$\begin{aligned} |\Lambda _0| \le CC(\kappa _{err}) h^{2k+2} + \frac{\epsilon _0\epsilon _{\infty }\kappa _{err}}{12}\Vert \varvec{\xi _E}\Vert ^2 + \frac{\epsilon _{0}a\theta }{8\omega _v^2}\Vert \xi _{\sigma }\Vert ^2. \end{aligned}$$

(A26)

For $\Lambda _1$, with that $\partial _{t}$ and ${\mathcal {P}}_h^{H_z}$ commute, and Lemma A1, we have

$$\begin{aligned} |\Lambda _1|\le \mu _0\Vert \partial _{t}H_{z} - {\mathcal {P}}_h^{H_z}\partial _{t}H_{z}\Vert ^2 +\frac{\mu _0}{4}\Vert \xi _{H_{z}}\Vert ^2 \le Ch^{2k+2} +\frac{\mu _0}{4}\Vert \xi _{H_z}\Vert ^2. \end{aligned}$$

(A27)

$\mathbf{Step \,\,\mathbf{3}.3: \,\,\mathbf{to} \,\,\mathbf{estimate} \,\, \Lambda _2 \,\,\mathbf{in} \,\, (A24)}$. For $\Lambda _2$, when the central fluxes (27) are used, we have

$$\begin{aligned} |\Lambda _2|&= \left| \sum _{i=1}^{N_x}\int _{y_{a}}^{y_{b}}\left( \widehat{\eta _{E_{y}}}(y)[\xi _{H_z}] + \widetilde{\eta _{H_{z}}}(y)[\xi _{E_y}]\right) _{x_{i+\frac{1}{2}}}dy\right. \nonumber \\&\quad \left. - \sum _{j=1}^{N_y}\int _{x_{a}}^{x_{b}} \left( \widehat{\widehat{\eta _{E_{x}}}}(x)[\xi _{H_z}] + \widetilde{\widetilde{\eta _{H_{z}}}}(x)[\xi _{E_x}]\right) _{y_{j+\frac{1}{2}}}dx\right| \nonumber \\&\le CC(\kappa _{err})h^{2k} + \frac{\mu _0}{4}\Vert \xi _{H_z}\Vert ^2 + \frac{\epsilon _0\epsilon _{\infty }\kappa _{err}}{12}\Vert \varvec{\xi _{E}}\Vert ^2. \end{aligned}$$

(A28)

The inverse inequality in (A5) is applied. When the alternating fluxes (28) are used with $\dagger =+, \ddagger =-, \natural =+, \sharp =-$, there holds

$$\begin{aligned} |\Lambda _2|&= \left| (\eta _{H_{z}},\partial _x\xi _{E_y}) + \sum _{i=1}^{N_x}\int _{y_{a}}^{y_{b}} \widetilde{\eta _{H_{z}}}(y)[\xi _{E_y}]_{x_{i+\frac{1}{2}}}dy - (\eta _{H_{z}},\partial _y\xi _{E_x}) \right. \nonumber \\&\quad \left. - \sum _{j=1}^{N_y}\int _{x_{a}}^{x_{b}} \widetilde{\widetilde{\eta _{H_{z}}}}(x)[\xi _{E_x}]_{y_{j+\frac{1}{2}}}dx\right| \nonumber \\&\le CC(\kappa _{err})h^{2k+2} + \frac{\epsilon _0\epsilon _{\infty }\kappa _{err}}{12}\Vert \varvec{\xi _{E}}\Vert ^2. \end{aligned}$$

(A29)

Here for the case of the alternating fluxes, we have used the super-convergence result established in [20] (see its Lemma 3.4) to gain the extra $h^2$ factor, and this is important for us to obtain the optimal error estimates.

$\mathbf{Step \,\,\mathbf{3}.4: \,\,\mathbf{to} \,\,\mathbf{estimate} \,\, \Lambda _3 \,\,\mathbf{in} \,\, (A24)}$. For $\Lambda _3$, we first bound the second term based on Lemma 1 and Lemma A3,

$$\begin{aligned}&\Big |\big ({\mathcal {I}}_h\left( 2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2\right) ,\xi _{\sigma }\big )\Big |=\Big |\int _\Omega {\mathcal {I}}_h\Big (\left( 2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2\right) \xi _{\sigma }\Big )d\Omega \Big | \nonumber \\&\quad \le C_\star \left\| 2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2\right\| \Vert \xi _{\sigma }\Vert \le C_\star \omega _v^2\left\| 2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2\right\| ^2 +\frac{1}{8\omega _v^2}\Vert \xi _{\sigma }\Vert ^2. \end{aligned}$$

(A30)

Using the regularity assumption on $\varvec{E}$ and Lemma A1, we get

$$\begin{aligned} \Vert \varvec{E}\varvec{\cdot }\varvec{\eta _E}\Vert \le \Vert \varvec{E}\Vert _{\infty }\Vert \varvec{\eta _E}\Vert ,\quad \Vert |\varvec{\eta _E}|^2\Vert \le \Vert \varvec{\eta _E}\Vert _{\infty }\Vert \varvec{\eta _E}\Vert \le C_k \Vert \varvec{E}\Vert _{\infty }\Vert \varvec{\eta _E}\Vert , \end{aligned}$$

(A31)

and therefore

$$\begin{aligned}&\Big |\big ({\mathcal {I}}_h\left( 2\varvec{E}\varvec{\cdot }\varvec{\eta _{E}}-|\varvec{\eta _{E}}|^2\right) ,\xi _{\sigma }\big )\Big |\le C\Vert \varvec{\eta _E}\Vert ^2 +\frac{1}{8\omega _v^2}\Vert \xi _{\sigma }\Vert ^2. \end{aligned}$$

(A32)

Apply Young’s inequality and the approximation properties in Lemma A1, we now bound $\Lambda _3$ as

$$\begin{aligned} |\Lambda _3| \le&\frac{\epsilon _0a\theta \gamma _v}{2\omega _v^2}\left( \frac{1}{4}\Vert \eta _{\sigma }\Vert ^2+\Vert \xi _{\sigma }\Vert ^2 \right) +\frac{\epsilon _0a\theta }{2\omega _v^2}\left( 2\Vert \partial _{t}\eta _\sigma +\omega _v^2\eta _{Q}\Vert ^2+ C\Vert \varvec{\eta _E}\Vert ^2 + \frac{1}{4}\Vert \xi _{\sigma }\Vert ^2\right) \nonumber \\ \le&Ch^{2k+2} + \frac{\epsilon _0a\theta \gamma _v}{2\omega _v^2}\Vert \xi _{\sigma }\Vert ^2 + \frac{\epsilon _0a\theta }{8\omega _v^2}\Vert \xi _{\sigma }\Vert ^2. \end{aligned}$$

(A33)

$\mathbf{Step \,\,\mathbf{3}.5: \,\,\mathbf{to} \,\,\mathbf{estimate} \,\, \Lambda _4 \,\,\mathbf{in} \,\, (A24)}$.

For $\Lambda _4$, using Cauchy Schwartz inequality, Young’s inequality, Lemma 1, and Lemma A3, we have

$$\begin{aligned} |\Lambda _4|&\le \epsilon _0\epsilon _\infty \left( \frac{1}{4\alpha _1}\Vert \partial _{t}\varvec{\eta _{E}}\Vert ^2+\alpha _1\Vert \varvec{\xi _{E}}\Vert ^2\right) + \epsilon _0\left( \frac{1}{4\alpha _2}\Vert \varvec{\eta _{J}}\Vert ^2+\alpha _2\Vert \varvec{\xi _{E}}\Vert ^2\right) \\&\quad +\frac{\epsilon _0 \omega _p^2}{2}\Vert \varvec{\mho _E}\Vert ^2+ \frac{\epsilon _0}{2\omega _p^2} \Vert \varvec{\xi _{J}}\Vert ^2\\&\quad + 3\epsilon _0a(1-\theta )\Vert \varvec{E}-\varvec{\eta _{E}}\Vert _\infty \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}})\Vert _\infty \Vert \varvec{\xi _{E}}\Vert ^2 + \epsilon _0a(1-\theta )\\&\quad \Big (\frac{1}{\alpha _3} \int _\Omega {\mathcal {I}}_h(|\varvec{\xi _{E}}|^4) d\Omega + \frac{\alpha _3}{4} \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}}) \Vert _\infty ^2 \Vert \varvec{\xi _{E}}\Vert ^2 \Big )\\&\quad + \epsilon _0a\theta C_\star \Big (\frac{1}{4\alpha _4}\left\| \partial _{t}\Big (Q\varvec{\eta _{E}}+\eta _Q(\varvec{E}-\varvec{\eta _{E}})\Big )\right\| ^2+\alpha _4\Vert \varvec{\xi _{E}}\Vert ^2\Big )\\&\quad + \frac{\epsilon _0a\theta }{2}\Vert \partial _{t}(Q-\eta _Q)\Vert _{\infty }\Vert \varvec{\xi _{E}}\Vert ^2 \\&\quad + \epsilon _0a(1-\theta )C_\star \left( \frac{1}{4\alpha _5}\left\| \partial _{t}\Big (|\varvec{\eta _{E}}|^2\varvec{\eta _{E}} - |\varvec{\eta _{E}}|^2\varvec{E} \right. \right. \\&\quad \left. \left. - 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{\eta _{E}} + 2(\varvec{E}\varvec{\cdot }\varvec{\eta _{E}})\varvec{E} + |\varvec{E}|^2\varvec{\eta _{E}}\Big )\right\| ^2 + \alpha _5\Vert \varvec{\xi _{E}}\Vert ^2\right) \\&\quad + \frac{\epsilon _0a\theta }{2} \Big (\frac{1}{4\alpha _6}\Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}})\Vert _\infty ^2\Vert \varvec{\xi _{E}}\Vert ^2 + \alpha _6\Vert \xi _Q\Vert ^2\Big ), \quad \forall \alpha _i>0, i=1,\cdots , 6. \end{aligned}$$

We only choose two terms from $\Lambda _4$ to illustrate how the estimation above is obtained,

$$\begin{aligned}&\left| \int _{\Omega }{\mathcal {I}}_h\left( \partial _{t}(Q-\eta _Q)|\varvec{\xi _{E}}|^2\right) d\Omega \right| \\&\quad \le \left| \sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\frac{\Delta x_i}{2}\frac{\Delta y_j}{2}\sum _{m,n=0}^{k}{\widehat{\omega }}_m {\widehat{\omega }}_n \Big (\partial _{t}(Q-\eta _Q)|\varvec{\xi _{E}}|^2\Big )(x_{im},y_{jn})\right| \\&\quad \le \Vert \partial _{t}(Q-\eta _Q)\Vert _\infty \left| \sum _{i=1}^{N_x}\sum _{j=1}^{N_y}\frac{\Delta x_i}{2}\frac{\Delta y_j}{2}\sum _{m,n=0}^{k}{\widehat{\omega }}_m {\widehat{\omega }}_n |\varvec{\xi _{E}}|^2(x_{im},y_{jn})\right| \\&\quad = \Vert \partial _{t}(Q-\eta _Q)\Vert _\infty \Vert \varvec{\xi _{E}}\Vert ^2, \end{aligned}$$

and

$$\begin{aligned}&\left| \int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2\partial _{t}(\varvec{E}-\varvec{\eta _{E}})\varvec{\cdot }\varvec{\xi _{E}}\right) d\Omega \right| \quad (\text{ by } \text{ Lemmas } \text{1-2 })\\&\quad \le \int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2| \partial _{t}(\varvec{E}-\varvec{\eta _{E}})| |\varvec{\xi _{E}}|\right) d\Omega = \int _{\Omega }{\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2\right) {\mathcal {I}}_h\left( |\partial _{t}(\varvec{E}-\varvec{\eta _{E}})|\; |\varvec{\xi _{E}}|\right) d\Omega \\&\quad \le \Vert {\mathcal {I}}_h\left( |\varvec{\xi _{E}}|^2\right) \Vert \Vert {\mathcal {I}}_h\left( |\partial _{t}(\varvec{E}-\varvec{\eta _{E}})|\; |\varvec{\xi _{E}}|\right) \Vert \\&\quad \le \left( \int _\Omega {\mathcal {I}}_h(|\varvec{\xi _{E}}|^4) d\Omega \right) ^{1/2} \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}}) \Vert _\infty \Vert \varvec{\xi _{E}}\Vert \\&\quad \le \frac{1}{\alpha _3} \int _\Omega {\mathcal {I}}_h(|\varvec{\xi _{E}}|^4) d\Omega + \frac{\alpha _3}{4} \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}}) \Vert _\infty ^2 \Vert \varvec{\xi _{E}}\Vert ^2. \end{aligned}$$

We now specify $\alpha _i,i=1,\cdots ,6$ as follows,

$$\begin{aligned} \epsilon _0\epsilon _\infty \alpha _1 = \epsilon _0\alpha _2 = \epsilon _0a\theta C_\star \alpha _4 = \epsilon _0a(1-\theta )C_\star \alpha _5 = \frac{1}{4}\frac{\epsilon _0\epsilon _\infty \kappa _{err}}{12},\quad \alpha _3 = \frac{12}{\rho _{err}},\quad \alpha _6 = \frac{\rho _{err}}{2}. \end{aligned}$$

Moreover we restrict the strength of nonlinearity such that

$$\begin{aligned}&3\epsilon _0a(1-\theta )\Vert \varvec{E}-\varvec{\eta _{E}}\Vert _\infty \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}})\Vert _\infty + \left( \epsilon _0a(1-\theta )\frac{\alpha _3}{4} + \frac{\epsilon _0a\theta }{8\alpha _6}\right) \Vert \partial _{t}(\varvec{E}-\varvec{\eta _{E}})\Vert _{\infty }^2 \\&\quad + \frac{\epsilon _0a\theta }{2}\Vert \partial _{t}(Q-\eta _Q)\Vert _{\infty } \le \frac{\epsilon _0\epsilon _\infty \kappa _{err}}{4}. \end{aligned}$$

This can be guaranteed if

$$\begin{aligned} 3a(1-\theta ) C_k^2\Vert \varvec{E}\Vert _\infty \Vert \partial _{t}\varvec{E}\Vert _\infty + (12-11\theta )\frac{a C_k^2}{4\rho _{err}} \Vert \partial _{t}\varvec{E}\Vert _{\infty }^2 + \frac{a\theta }{2} C_k\Vert \partial _{t}Q\Vert _{\infty } \le \frac{\epsilon _\infty \kappa _{err}}{4}, \end{aligned}$$

holds, which is also Condition 3 in (41). Based on what we have by far, with the approximation properties from Lemma A1, we arrive at

$$\begin{aligned} |\Lambda _4|\le & {} CC(\kappa _{err},\rho _{err})h^{2k+2} + \frac{\epsilon _0\epsilon _\infty \kappa _{err}}{3}\Vert \varvec{\xi _E}\Vert ^2 + \frac{\epsilon _0 a(1-\theta )\rho _{err}}{12} \int _\Omega {\mathcal {I}}_h(|\varvec{\xi _{E}}|^4) d\Omega \nonumber \\&+ \frac{\epsilon _0a\theta \rho _{err}}{4}\Vert \xi _Q\Vert ^2+\frac{\epsilon _0}{2\omega _p^2} \Vert \varvec{\xi _{J}}\Vert ^2. \end{aligned}$$

(A34)

The combination of (A24)–(A34) gives

$$\begin{aligned} \frac{d\widehat{{\mathcal {E}}_h}}{dt}\le \widehat{{\mathcal {E}}_h} + CC(\kappa _{err},\rho _{err})h^{2r}, \end{aligned}$$

where

$$\begin{aligned} r=\left\{ \begin{array}{ll} k,\quad &{}\text{ for } \text{ central } \text{ numerical } \text{ fluxes, }\\ k+1, \quad &{}\text{ for } \text{ alternating } \text{ numerical } \text{ fluxes. } \end{array} \right. \end{aligned}$$

Finally, using Gronwall inequality and $\widehat{{\mathcal {E}}_h}(0)=0$, we get ${\mathcal {E}}_h^{(mod)}\le \widehat{{\mathcal {E}}_h}\le CC(\kappa _{err},\rho _{err})h^{r}$. We further apply a triangle inequality and the approximation results in Lemma A1, and conclude

$$\begin{aligned} \Vert u-u_h\Vert \le \Vert \xi _u\Vert + \Vert \eta _u\Vert \le CC(\kappa _{err},\rho _{err})h^{r},\quad u=H_z,\varvec{E},\varvec{P},\varvec{J},\sigma ,Q. \end{aligned}$$

$\square$

Appendix B: Proof of Theorem 4

Proof

Apply two time steps to (42a) and (42j), sum them up after taking $\phi =H_{zh}^{n+1/2}$, one gets

$$\begin{aligned}&\mu _0\left( H_{zh}^{n+3/2}-H_{zh}^{n+1/2},H_{zh}^{n-1/2}\right) + \Delta t {\mathcal {B}}_{h}^{\varvec{E}}(E_{xh}^{n+1}+E_{xh}^{n},E_{yh}^{n+1}+E_{yh}^{n},H_{zh}^{n+1/2}) = 0. \end{aligned}$$

(B1)

Take $\phi =E_{xh}^{n+1}+E_{xh}^{n}$ in (42b), $\phi =E_{yh}^{n+1}+E_{yh}^{n}$ in (42c), sum them up, we have

$$\begin{aligned}&\displaystyle \left( \varvec{D}_{h}^{n+1}\!-\!\varvec{D}_{h}^{n} , \varvec{E}_{h}^{n+1} + \varvec{E}_{h}^{n}\right) + \Delta t{\mathcal {B}}_{xh}^{H}(H_{zh}^{n+1/2},E_{xh}^{n+1}+E_{xh}^{n})+ \Delta t{\mathcal {B}}_{yh}^{H}(H_{zh}^{n+1/2},E_{yh}^{n+1} + E_{yh}^{n}) = 0.\nonumber \\ \end{aligned}$$

(B2)

Add (B2) to (B1) and using the identity (30), we obtain

$$\begin{aligned}&\mu _0\int _{\Omega }H_{zh}^{n+3/2}H_{zh}^{n+1/2}d\Omega - \mu _0\int _{\Omega }H_{zh}^{n+1/2}H_{zh}^{n-1/2}d\Omega \nonumber \\&\quad + \int _{\Omega }\left( \varvec{D}_{h}^{n+1}-\varvec{D}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega =0. \end{aligned}$$

(B3)

By (42d), we have

$$\begin{aligned}&\int _{\Omega }\left( \varvec{D}_{h}^{n+1}-\varvec{D}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\epsilon _0\epsilon _\infty \int _{\Omega }|\varvec{E}_{h}^{n+1}|^2d\Omega - \epsilon _0\epsilon _\infty \int _{\Omega }|\varvec{E}_{h}^{n}|^2d\Omega + \epsilon _0\int _{\Omega }\left( \varvec{P}_{h}^{n+1}-\varvec{P}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\qquad + \epsilon _0a(1-\theta )\int _{\Omega }\left( \varvec{Y}_{h}^{n+1}-\varvec{Y}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\qquad + \epsilon _0a\theta \int _{\Omega }\left( {\mathcal {I}}_h\left( Q_{h}^{n+1}\varvec{E}_{h}^{n+1}\right) -{\mathcal {I}}_h\left( Q_{h}^{n}\varvec{E}_{h}^{n}\right) \right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega . \end{aligned}$$

(B4)

For the nonlinear Kerr term, by (42e) and Lemma 1, we get

$$\begin{aligned}&\int _{\Omega }\left( \varvec{Y}_{h}^{n+1}-\varvec{Y}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\int _{\Omega }{\mathcal {I}}_h\left( \left( |\varvec{E}^{n+1}|^2+|\varvec{E}^{n}|^2-\varvec{E}^{n+1}\varvec{\cdot }\varvec{E}^{n}\right) \left( \varvec{E}^{n+1}-\varvec{E}^{n}\right) \right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\qquad +\frac{1}{2}\int _{\Omega }{\mathcal {I}}_h\left( (\varvec{E}^{n+1} +\varvec{E}^{n})\varvec{\cdot }(\varvec{E}^{n+1} - \varvec{E}^{n})(\varvec{E}^{n+1}+\varvec{E}^{n})\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\frac{1}{2}\int _{\Omega }{\mathcal {I}}_h\Big (\left( 2|\varvec{E}^{n+1}|^2+2|\varvec{E}^{n}|^2-2\varvec{E}^{n+1}\varvec{\cdot }\varvec{E}^{n}+|\varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}|^2 \right) \left( |\varvec{E}^{n+1}|^2-|\varvec{E}^{n}|^2\right) \Big )d\Omega \nonumber \\&\quad =\frac{3}{2}\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{E}^{n+1}|^4\right) d\Omega - \frac{3}{2}\int _{\Omega }{\mathcal {I}}_h\left( |\varvec{E}^{n}|^4\right) d\Omega . \end{aligned}$$

(B5)

For the Lorentz term, using (42f) and (42g), we have

$$\begin{aligned}&\int _{\Omega }\left( \varvec{P}_{h}^{n+1}-\varvec{P}_{h}^{n}\right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\frac{2}{\omega _p^2}\int _{\Omega }\left( \varvec{P}_{h}^{n+1}-\varvec{P}_{h}^{n}\right) \varvec{\cdot }\left( \frac{\varvec{J}_{h}^{n+1}-\varvec{J}_{h}^{n}}{\Delta t} + \gamma \frac{\varvec{J}_{h}^{n+1}+\varvec{J}_{h}^{n}}{2} + \omega _{0}^2\frac{\varvec{P}_{h}^{n+1}+\varvec{P}_{h}^{n}}{2}\right) d\Omega \nonumber \\&\quad =\frac{1}{\omega _p^2}\int _{\Omega }|\varvec{J}_{h}^{n+1}|^2d\Omega - \frac{1}{\omega _p^2}\int _{\Omega }|\varvec{J}_{h}^{n}|^2d\Omega + \frac{\gamma \Delta t}{2\omega _p^2} \int _{\Omega }|\varvec{J}_{h}^{n+1}+\varvec{J}_{h}^{n}|^2d\Omega \nonumber \\&\qquad + \frac{\omega _0^2}{\omega _p^2}\int _{\Omega }|\varvec{P}_{h}^{n+1}|^2d\Omega - \frac{\omega _0^2}{\omega _p^2}\int _{\Omega }|\varvec{P}_{h}^{n}|^2d\Omega . \end{aligned}$$

(B6)

For the nonlinear Raman term, we have

$$\begin{aligned}&\int _{\Omega }\left( {\mathcal {I}}_h\left( Q_{h}^{n+1}\varvec{E}_{h}^{n+1}\right) -{\mathcal {I}}_h\left( Q_{h}^{n}\varvec{E}_{h}^{n}\right) \right) \varvec{\cdot }\left( \varvec{E}_{h}^{n+1}+\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\int _{\Omega }{\mathcal {I}}_h\left( Q_{h}^{n+1}|\varvec{E}_{h}^{n+1}|^2\right) d\Omega - \int _{\Omega }{\mathcal {I}}_h\left( Q_{h}^{n}|\varvec{E}_{h}^{n}|^2\right) d\Omega \nonumber \\&\qquad + \int _{\Omega }(Q_{h}^{n+1}-Q_{h}^{n}){\mathcal {I}}_h\left( \varvec{E}_{h}^{n+1}\varvec{\cdot }\varvec{E}_{h}^{n}\right) d\Omega , \end{aligned}$$

(B7)

where the last term can be further proceeded based on (42h) and (42i),

$$\begin{aligned}&\int _{\Omega }(Q_{h}^{n+1}-Q_{h}^{n}){\mathcal {I}}_h\left( \varvec{E}_{h}^{n+1}\varvec{\cdot }\varvec{E}_{h}^{n}\right) d\Omega \nonumber \\&\quad =\frac{1}{\omega _v^2}\int _{\Omega }(Q_{h}^{n+1}-Q_{h}^{n})\left( \frac{\sigma _{h}^{n+1}-\sigma _{h}^{n}}{\Delta t} + \gamma _v \frac{\sigma _{h}^{n+1}+\sigma _{h}^{n}}{2} + \omega _{v}^2\frac{Q_{h}^{n+1}+Q_{h}^{n}}{2}\right) d\Omega \nonumber \\&\quad = \frac{1}{2\omega _v^2}\int _{\Omega }(\sigma _{h}^{n+1})^2d\Omega - \frac{1}{2\omega _v^2}\int _{\Omega }(\sigma _{h}^{n})^2d\Omega \nonumber \\&\qquad + \frac{\gamma _v\Delta t}{4\omega _v^2}\int _{\Omega }(\sigma _{h}^{n+1}+\sigma _{h}^{n})^2d\Omega + \frac{1}{2}\int _{\Omega }(Q_{h}^{n+1})^2d\Omega - \frac{1}{2}\int _{\Omega }(Q_{h}^{n})^2d\Omega . \end{aligned}$$

(B8)

Combine (B3)–(B8), we reach the energy relation (43) with the discrete energy ${\mathcal {E}}_h^n$ defined in (44).

The next step is to derive the time step condition under which the discrete energy ${\mathcal {E}}_h^n$ is nonnegative. Based on (42a) and (42j), we have

$$\begin{aligned}&\mu _0\left( H_{zh}^{n+1/2}-H_{zh}^{n-1/2},\phi \right) - \Delta t\sum _{i=1}^{N_x}\int _{y_{a}}^{y_{b}} \widehat{E_{yh}^{n}}[\phi ]_{x_{i+\frac{1}{2}}}dy + \Delta t\sum _{j=1}^{N_y}\int _{x_{a}}^{x_{b}} \widehat{\widehat{E_{xh}^{n}}}[\phi ]_{y_{j+\frac{1}{2}}}dx\nonumber \\&\quad - \Delta t(E_{yh}^{n},\partial _x\phi ) + \Delta t(E_{xh}^{n},\partial _y\phi ) = 0. \end{aligned}$$

(B9)

Take $\phi =H_{zh}^{n+1/2}$, using the inverse inequality (A5), we obtain

$$\begin{aligned}&\mu _0\int _{\Omega }\left( H_{zh}^{n+1/2}\right) ^2d\Omega - \mu _0\int _{\Omega }H_{zh}^{n+1/2}H_{zh}^{n-1/2}d\Omega \nonumber \\&\quad = \Delta t\sum _{i=1}^{N_x}\int _{y_{a}}^{y_{b}} \widehat{E_{yh}^{n}}\left[ H_{zh}^{n+1/2}\right] _{x_{i+\frac{1}{2}}}dy - \Delta t\sum _{j=1}^{N_y}\int _{x_{a}}^{x_{b}} \widehat{\widehat{E_{xh}^{n}}}\left[ H_{zh}^{n+1/2}\right] _{y_{j+\frac{1}{2}}}dx\nonumber \\&\qquad + \Delta t(E_{yh}^{n},\partial _xH_{zh}^{n+1/2}) - \Delta t(E_{xh}^{n},\partial _yH_{zh}^{n+1/2}) \nonumber \\&\quad \le \frac{\Delta t C_\star }{h}\left( \int _{\Omega }\left| \varvec{E}_{h}^{n}\right| ^2d\Omega + \int _{\Omega }\left( H_{zh}^{n+1/2}\right) ^2d\Omega \right) , \end{aligned}$$

(B10)

hence,

$$\begin{aligned} \frac{\mu _0}{2}\int _{\Omega }H_{zh}^{n+1/2}H_{zh}^{n-1/2}d\Omega \ge \left( \frac{\mu _0}{2}-\frac{\Delta t C_\star }{2h}\right) \int _{\Omega }\left( H_{zh}^{n+1/2}\right) ^2d\Omega - \frac{\Delta t C_\star }{2h}\int _{\Omega }\left| \varvec{E}_{h}^{n}\right| ^2d\Omega . \end{aligned}$$

This implies, if we restrict the time step $\Delta t$ such that,

$$\begin{aligned} \frac{\mu _0}{2}-\frac{\Delta t C_\star }{2h}\ge 0,\quad \frac{\epsilon _0\epsilon _\infty }{2}-\frac{\Delta t C_\star }{2h}\ge 0, \end{aligned}$$

then by Lemma 1 and with the condition $\theta \in [0,\frac{3}{4}]$, we have

$$\begin{aligned} {\mathcal {E}}_h^n \ge \int _{\Omega } \frac{\epsilon _0a\theta }{4}{\mathcal {I}}_h\left( \left( |\varvec{E}_h^n|^2 + \left( Q_h^n\right) ^2\right) ^2\right) + \frac{\epsilon _0a(3-4\theta )}{4}{\mathcal {I}}_h\left( |\varvec{E}_h^n|^4 \right) d\Omega \ge 0. \end{aligned}$$

$\square$