A novel finite element method for simulating surface plasmon polaritons on complex graphene sheets

Abstract

Surface plasmon polaritons (SPPs) are generated on the graphene surface, and provide a window into the nano-optical and electrodynamic response of their host material and its dielectric environment. An accurate simulation of SPPs presents several unique challenges, since SPPs often occur at complex interfaces between materials of different dielectric constants and appropriate boundary conditions at the graphene interfaces are crucial. Here we develop a simplified graphene model and propose a new finite element method accordingly. Stability for the continuous model is established, and extensive numerical results are presented to demonstrate that the new model can capture the SPPs very well for various complex graphene sheets.

Keywords – Maxwell’s equations, finite element time-domain methods, edge elements, graphene, surface plasmon polaritons.

Mathematics Subject Classification (2000): 78M10, 65N30, 65F10, 78-08.

1 Introduction

Graphene is a single layer of carbon atoms arranged in a hexagonal lattice pattern, often described as a ”honeycomb” structure. The 2-D material graphene was first successfully isolated by Novoselov and Geim et al. [33]. In 2010, Geim and Novoselov were awarded Nobel Prizes in Physics for their groundbreaking experiments regarding graphene. Due to graphene’s unique electrical, electromagnetic, and optical characteristics, it has attracted widespread attention, leading to the design of many new systems and equipment with graphene. For example, graphene has played a prominent role in the design of organic light-emitting diodes, solar cells, antennas, and invisibility cloaks (cf. [12, 5, 40]).

In the past two decades, in addition to many great achievements on the physical research on graphene and graphene-based devices, another hot topic has been the numerical simulation of graphene. The famous Kubo formula [13] gives the expression of graphene conductivity, which is a function of many physical parameters such as wavelength, chemical potential, and temperature. With the development of computational electromagnetic in the past, such as the finite difference time-domain (FDTD) method (e.g., [14, 16, 19, 27, 39]) and the finite element method (FEM) (e.g., papers [4, 7, 8, 9, 23, 35, 17, 37, 42], and books [11, 30, 24]), many robust and efficient numerical methods have been proposed to simulate the electromagnetic response of graphene related devices. The numerical approaches for modeling graphene can be classified into two big categories: (1) Treating graphene as a thin plate with finite thickness [34], and converting its surface conductivity into a volumetric conductivity; (2) Taking graphene as a zero-thickness sheet [31]. Due to the easy realization of (1), many published papers and commercial software model graphene as a thin sheet with a finite thickness [6, 43]. However, direct discretization of graphene with a small finite thickness results in extremely fine grids around graphene. This leads to extremely small time steps for time-domain simulations with explicit schemes, which consume enormous amount of memory storage and CPU time.

The interesting physical research on graphene and graphene-based devices has inspired mathematicians to conduct mathematical analysis [1, 15, 21] and modeling of graphene (e.g., [20, 29, 38, 44] and references therein). Here we are interested in simulating the surface plasmon polaritons (SPPs) generated on the graphene surface, since SPPs provide a window into the nano-optical, electrodynamic response of their host material and its dielectric environment [41]. The plasmons occur in the highly sought after terahertz to mid-infrared regime. Note that terahertz waves are used in a variety of applications such as nondestructive analysis, since terahertz waves can be used to analyze the internal structure of objects without damaging them. For example, THz cameras can be used to see what is inside sealed packages. Also, terahertz waves can be found in military applications, e.g., terahertz sensors are used in military communication, detecting biological and explosive warfare agents, and inspecting concealed weapons. However, because of practical difficulties in exciting and detecting the SPP waves in graphene, numerical simulation of wave interactions with graphene materials plays a very important role in designing functional components with graphene. Considering the disadvantage of FDTD methods in handling the complex geometry, which happens quite often in graphene devices, we recently proposed and analyzed some finite element time-domain (FETD) methods for graphene simulation [45, 22, 18, 26]. In [45, 22], we treated the graphene with some thickness (though very thin); while in [26], we successfully developed a new FETD method for simulating surface plasmon polaritons (SPPs) on graphene sheets with zero-thickness. We like to remark that the works of [29, 38] are based on frequency-domain finite element methods. The work of Wilson, Santosa, and Martin [44] is for 2-D time-domain SPP models written as an integro-differential equation in time. In [32], Nicholls et al. proposed the so-called high–order perturbation of surfaces algorithms for simulating the plasmonic response of a perfectly flat sheet of graphene. They solved the Helmholtz equations by using the dispersive Drude model for the surface current.

In this paper, we propose a novel FETD method by treating the graphene sheet as zero-thickness. The new method is based on a reformulated system of governing equations for graphene with electric and magnetic fields as unknowns. Compared to our previous work [26], the new system does not contain the induced current explicitly. Hence, our new method is more efficient in memory and computational cost.

The rest of the paper is organized as follows. In Section 2, we first present and reformulate the time-domain governing equations for modeling the surface plasmon polaritons on graphene sheets. Then we prove the stability for the model. In Section 3, we propose a leapfrog time stepping finite element scheme for solving the graphene model. In Section 4, we present extensive numerical results to demonstrate the effectiveness of our method for simulating the propagation of surface plasmon polaritons on various complex graphene sheets. We conclude the paper in Section 5.

2 The governing equations and stability analysis

In our previous papers [45, 22], we treated the graphene as a homogenized material with an effective permittivity and a small thickness. By ignoring the interband conductivity, we have the $TE_{z}$ model governing equations for simulating surface plasmon propagation on graphene [45, (2.7)-(2.12)]:

$$\epsilon_{0}\partial_{t}\mbox{\boldmath$E$}=\nabla\times H,\quad\mbox{in}~{}~{}\Omega,$$ (2.1)

$$\mu_{0}\partial_{t}H=-\nabla\times\mbox{\boldmath$E$}-K_{s},\quad\mbox{in}~{}~{}\Omega,$$ (2.2)

$$\tau_{0}\partial_{t}\mbox{\boldmath$J$}+\mbox{\boldmath$J$}=\sigma_{0}\mbox{\boldmath$E$},\quad\mbox{on}~{}\Gamma,$$ (2.3)

where we denote the electric field $\mbox{\boldmath$E$}=(E_{x},E_{y})^{\prime}$, magnetic field $H=H_{z}$, $K_{s}$ for an imposed magnetic source function, $\mbox{\boldmath$J$}:=\mbox{\boldmath$J$}_{d}$ (as denoted in [45]) for the induced intraband surface current in graphene, $\epsilon_{0}$ and $\mu_{0}$ for the vacuum permittivity and permeability, the positive constant $\tau_{0}$ for the relaxation time, and the positive constant $\sigma_{0}$ for the graphene surface conductivity. Moreover, we assume that the physical domain $\Omega$ is a bounded Lipschitz polygonal domain in ${\cal R}^{2}$ with boundary $\partial\Omega$, and $\Gamma$ represents the graphene sheet embedded in $\Omega$. Also we define the 2-D curl operators as $\nabla\times H:=(\partial_{y}H,-\partial_{x}H)^{\prime}$ and $\nabla\times\mbox{\boldmath$E$}:=\partial_{x}E_{y}-\partial_{y}E_{x}$.

To complete the model, we assume that (2.1)-(2.3) are subject to the simple perfectly conducting (PEC) boundary condition:

$$\bm{\hat{\nu}}\times\mbox{\boldmath$E$}=0,~{}~{}\mbox{on}~{}\partial\Omega,$$ (2.4)

and the initial conditions

$$\mbox{\boldmath$E$}(\mbox{\boldmath$x$},0)=\mbox{\boldmath$E$}_{0}(\mbox{\boldmath$x$}),~{}H(\mbox{\boldmath$x$},0)=H_{0}(\mbox{\boldmath$x$}),~{}\mbox{\boldmath$J$}(\mbox{\boldmath$x$},0)=\mbox{\boldmath$J$}_{0}(\mbox{\boldmath$x$}),$$ (2.5)

where $\hat{\bm{\nu}}$ denotes the unit outward normal vector on $\partial\Omega$, and $\mbox{\boldmath$E$}_{0},H_{0},\mbox{\boldmath$J$}_{0}$ are some properly given functions.

Recently, we [26] successfully simulated surface plasmon propagation on graphene by treating it as zero-thickness sheet (i.e., appearing as a curve in our 2-D simulations, e.g., Figures 1 and 3 shown later), with the following boundary conditions given on the graphene interface [3, Fig.1]:

$$\hat{n}_{1}\times{\mbox{\boldmath$E$}}_{1}=\hat{n}_{2}\times{\mbox{\boldmath$E$}}_{2},~{}~{}\mbox{on}~{}\Gamma,$$ (2.6)

$$H_{1}-H_{2}=\mbox{\boldmath$J$}\times\bm{\hat{n}},~{}~{}\mbox{on}~{}\Gamma,$$ (2.7)

where $H_{1}$ and $H_{2}$ represent the magnetic fields above and below the interface, respectively, $\bm{\hat{n}}:=(n_{x},n_{y})^{\prime}$ denotes the unit normal vector of the interface pointing upward, and $\hat{n}_{1}$ and $\hat{n}_{2}$ are the unit outward normal vectors from the top and bottom of the interface. Finally, we adopt the 2-D cross product notation $\mbox{\boldmath$J$}\times\bm{\hat{n}}:=J_{x}n_{y}-J_{y}n_{x}$. We like to remark that (2.6)-(2.7) imply that the tangential electric field is continuous across the interface, and the jump of the tangential component of the magnetic field across the interface is equal to the surface current.

In [26], using integration by parts and the interface conditions (2.6)-(2.7), we derived the following weak formulation for simulating graphene as an interface: Find the solution

$$\mbox{\boldmath$E$}\in L^{2}(0,T;H_{0}(curl;\Omega))\cap H^{1}(0,T;(L^{2}(\Omega))^{2}),H\in H^{1}(0,T;L^{2}(\Omega)),\mbox{\boldmath$J$}\in H^{1}(0,T;(L^{2}(\Gamma))^{2}),$$

such that (cf. [26, (2.9)-(2.11)]):

$$\epsilon_{0}(\partial_{t}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})=(H,\nabla\times\mbox{\boldmath${\bf\phi}$})-\langle\mbox{\boldmath$J$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma},$$ (2.8)

$$\mu_{0}(\partial_{t}H,\psi)=-(\nabla\times\mbox{\boldmath$E$},\psi)-(K_{s},\psi),$$ (2.9)

$$\langle\tau_{0}\partial_{t}\mbox{\boldmath$J$},\mbox{\boldmath${\bf\chi}$}\rangle_{\Gamma}+\langle\mbox{\boldmath$J$},\mbox{\boldmath${\bf\chi}$}\rangle_{\Gamma}=\langle\sigma_{0}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\chi}$}\rangle_{\Gamma},$$ (2.10)

hold true for any test functions $\mbox{\boldmath${\bf\phi}$}\in H_{0}(curl;\Omega),\psi\in L^{2}(\Omega)$, and $\mbox{\boldmath${\bf\chi}$}\in(L^{2}(\Gamma))^{2}$. To obtain (2.8), we use integration by parts over $\Omega$ and the boundary condition (2.7). Here we denote $(\cdot,\cdot)$ for the $L^{2}$ inner product over $\Omega$, and $\langle\mbox{\boldmath$J$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}:=\int_{\Gamma}\mbox{\boldmath$J$}\times\bm{\hat{n}}\cdot\mbox{\boldmath${\bf\phi}$}\times\bm{\hat{n}}~{}ds$ for the inner product on $\Gamma$. Finally, we adopt the standard Sobolev space notation

$$H_{0}(curl;\Omega)=\{\mbox{\boldmath$u$}\in(L^{2}(\Omega))^{2}:~{}\nabla\times\mbox{\boldmath$u$}\in L^{2}(\Omega),~{}\bm{\hat{\nu}}\times\mbox{\boldmath$u$}=0~{}\mbox{on}~{}\partial\Omega\}.$$

To derive the governing equations for our new numerical method, we differentiate (2.8) with respect to $t$ and then multiplying the result by $\tau_{0}$. This gives

$$\tau_{0}\epsilon_{0}(\partial_{tt}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})=\tau_{0}(\partial_{t}H,\nabla\times\mbox{\boldmath${\bf\phi}$})-\tau_{0}\langle\partial_{t}\mbox{\boldmath$J$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}.$$ (2.11)

Adding (2.11) and (2.8) together, and using (2.10) with $\mbox{\boldmath${\bf\chi}$}=\mbox{\boldmath${\bf\phi}$}$, we have

$$\tau_{0}\epsilon_{0}(\partial_{tt}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\epsilon_{0}(\partial_{t}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})=\tau_{0}(\partial_{t}H,\nabla\times\mbox{\boldmath${\bf\phi}$})+(H,\nabla\times\mbox{\boldmath${\bf\phi}$})-\sigma_{0}\langle\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}.$$ (2.12)

To develop a more efficient numerical method later, we replace $\partial_{t}H$ in (2.12) by (2.9) with $\psi=\nabla\times\mbox{\boldmath${\bf\phi}$}$ and obtain the following new weak formulation: Find the solution $\mbox{\boldmath$E$}\in L^{2}(0,T;H_{0}(curl;\Omega))\cap H^{2}(0,T;(L^{2}(\Omega))^{2})$, $H\in H^{1}(0,T;L^{2}(\Omega))$, such that

$$\tau_{0}\epsilon_{0}(\partial_{tt}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\epsilon_{0}(\partial_{t}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\frac{\tau_{0}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$},\nabla\times\mbox{\boldmath${\bf\phi}$})$$

$$\quad\quad=(H,\nabla\times\mbox{\boldmath${\bf\phi}$})-\frac{\tau_{0}}{\mu_{0}}(K_{s},\nabla\times\mbox{\boldmath${\bf\phi}$})-\sigma_{0}\langle\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma},$$ (2.13)

$$\mu_{0}(\partial_{t}H,\psi)=-(\nabla\times\mbox{\boldmath$E$},\psi)-(K_{s},\psi),$$ (2.14)

hold true for any test functions $\mbox{\boldmath${\bf\phi}$}\in H_{0}(curl;\Omega)$ and $\psi\in L^{2}(\Omega)$. The initial conditions for the problem (2.13)-(2.14) are as follows:

$$\mbox{\boldmath$E$}(\mbox{\boldmath$x$},0)=\mbox{\boldmath$E$}_{0}(\mbox{\boldmath$x$}),~{}H(\mbox{\boldmath$x$},0)=H_{0}(\mbox{\boldmath$x$}),~{}\partial_{t}\mbox{\boldmath$E$}(\mbox{\boldmath$x$},0)=\epsilon_{0}^{-1}\nabla\times H_{0}(\mbox{\boldmath$x$}).$$ (2.15)

To simplify the notation, we denote the $L^{2}$ norm of $u$ in $\Omega$ as $||u||:=||u||_{L^{2}(\Omega)}$, and the $L^{2}$ norm of $u$ on $\Gamma$ as $||\mbox{\boldmath$u$}||_{\Gamma}:=(\int_{\Gamma}|\mbox{\boldmath$u$}\times\bm{\hat{n}}|^{2}~{}ds)^{1/2}$.

Theorem 2.1.

Denote the energy

$$ENG(t) = \tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}+\mu_{0}\|\tau_{0}^{1/2}\partial_{t}H+\tau_{0}^{-1/2}H\|^{2}+\frac{\mu_{0}}{\tau_{0}}\|H\|^{2}$$ (2.16)

$$\quad+\frac{\epsilon_{0}}{\tau_{0}}\|\mbox{\boldmath$E$}\|^{2}+\frac{1}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}+\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}^{2}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}.$$

Then we have the following continuous stability:

$$ENG(t)\leq\left[ENG(0)+\int_{0}^{t}(\frac{\tau_{0}}{\mu_{0}}\|\partial_{t}K_{s}\|^{2}+\frac{1}{\tau_{0}\mu_{0}}\|K_{s}\|^{2})dt\right]\cdot\exp(C_{*}t),\quad\forall~{}t\in[0,T],$$ (2.17)

where the constant $C_{*}>0$ depends only on the parameter $\tau_{0}$.

Proof. To make our proof easy to follow, we split it into several major parts.

(I)  Choosing $\mbox{\boldmath${\bf\phi}$}=\partial_{t}\mbox{\boldmath$E$}$ in (2.13), we have

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\frac{\tau_{0}}{\mu_{0}}\|\nabla\times\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}\right)+\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}$$

$$\quad\quad=(H,\nabla\times\partial_{t}\mbox{\boldmath$E$})-\frac{\tau_{0}}{\mu_{0}}(K_{s},\nabla\times\partial_{t}\mbox{\boldmath$E$}).$$ (2.18)

Taking the time derivative of (2.2), we obtain

$$\nabla\times\partial_{t}\mbox{\boldmath$E$}=-\mu_{0}\partial_{tt}H-\partial_{t}K_{s}.$$ (2.19)

Replacing $\nabla\times\partial_{t}\mbox{\boldmath$E$}$ in (2.18) by (2.19), we have

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\frac{\tau_{0}}{\mu_{0}}\|\nabla\times\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}\right)+\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}$$

$$\quad\quad=-\mu_{0}(H,\partial_{tt}H)-(H,\partial_{t}K_{s})+\tau_{0}(K_{s},\partial_{tt}H)+\frac{\tau_{0}}{\mu_{0}}(K_{s},\partial_{t}K_{s}).$$ (2.20)

(II)  By choosing $\psi=\tau_{0}\partial_{tt}H$ in (2.14), then replacing $\partial_{tt}H$ by (2.19), we have

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}\right) = -\tau_{0}(\nabla\times\mbox{\boldmath$E$},\partial_{tt}H)-\tau_{0}(K_{s},\partial_{tt}H)$$ (2.21)

$$= \frac{\tau_{0}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$},\nabla\times\partial_{t}\mbox{\boldmath$E$}+\partial_{t}K_{s})-\tau_{0}(K_{s},\partial_{tt}H)$$

$$= \frac{\tau_{0}}{2\mu_{0}}\frac{d}{dt}(||\nabla\times\mbox{\boldmath$E$}||^{2})+\frac{\tau_{0}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$},\partial_{t}K_{s})-\tau_{0}(K_{s},\partial_{tt}H).$$

Adding (2.20) and (2.21) together, and using (2.2), we obtain

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}\right)+\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}$$ (2.22)

$$= -\mu_{0}(H,\partial_{tt}H)-(H,\partial_{t}K_{s})-\tau_{0}(\partial_{t}H,\partial_{t}K_{s})$$

$$= -\mu_{0}\left[\frac{d}{dt}(H,\partial_{t}H)-(\partial_{t}H,\partial_{t}H)\right]-(H,\partial_{t}K_{s})-\tau_{0}(\partial_{t}H,\partial_{t}K_{s}),$$

which can be rewritten as

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}+\mu_{0}\|\tau_{0}^{1/2}\partial_{t}H+\tau_{0}^{-1/2}H\|^{2}-\frac{\mu_{0}}{\tau_{0}}\|H\|^{2}\right)+\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}$$ (2.23)

$$= \mu_{0}||\partial_{t}H||^{2}-(H,\partial_{t}K_{s})-\tau_{0}(\partial_{t}H,\partial_{t}K_{s}).$$

(III) To bound the terms $H$ and $\partial_{t}H$ on the right-hand side (RHS) of (2.23), we first choose $\mbox{\boldmath${\bf\phi}$}=\mbox{\boldmath$E$},\psi=H,\mbox{\boldmath${\bf\chi}$}=\frac{1}{\sigma_{0}}\mbox{\boldmath$J$}$ in (2.8)-(2.10), respectively. Adding the results together, we obtain

$$\frac{1}{2}\frac{d}{dt}\left(\epsilon_{0}\|\mbox{\boldmath$E$}\|^{2}+\mu_{0}\|H\|^{2}+\frac{\tau_{0}}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)+\frac{1}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}=-(K_{s},H).$$ (2.24)

Multiplying (2.24) by $\frac{2}{\tau_{0}}$, we have

$$\frac{1}{2}\frac{d}{dt}\left(\frac{2\epsilon_{0}}{\tau_{0}}\|\mbox{\boldmath$E$}\|^{2}+\frac{2\mu_{0}}{\tau_{0}}\|H\|^{2}+\frac{2}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)+\frac{2}{\tau_{0}\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}=-\frac{2}{\tau_{0}}(K_{s},H).$$ (2.25)

Similarly, taking the time derivative of (2.8)-(2.10), then choosing $\mbox{\boldmath${\bf\phi}$}=\partial_{t}\mbox{\boldmath$E$},\psi=\partial_{t}H,\mbox{\boldmath${\bf\chi}$}=\frac{1}{\sigma_{0}}\partial_{t}\mbox{\boldmath$J$}$ respectively, and adding the results together, we obtain

$$\frac{1}{2}\frac{d}{dt}\left(\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)+\frac{1}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}=-(\partial_{t}K_{s},\partial_{t}H).$$ (2.26)

Multiplying (2.26) by $\tau_{0}$, we have

$$\frac{1}{2}\frac{d}{dt}\left(\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}^{2}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)+\frac{\tau_{0}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}=-\tau_{0}(\partial_{t}K_{s},\partial_{t}H).$$ (2.27)

Now adding (2.23), (2.25), and (2.27) together, we have

$$\frac{1}{2}\frac{d}{dt}\left(2\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}+\mu_{0}\|\tau_{0}^{1/2}\partial_{t}H+\tau_{0}^{-1/2}H\|^{2}+\frac{\mu_{0}}{\tau_{0}}\|H\|^{2}\right.$$ (2.28)

$$\quad\left.+\frac{2\epsilon_{0}}{\tau_{0}}\|\mbox{\boldmath$E$}\|^{2}+\frac{2}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}+\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}^{2}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)$$

$$\quad+\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\frac{2}{\tau_{0}\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}+\frac{\tau_{0}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}$$

$$= \mu_{0}||\partial_{t}H||^{2}-(H,\partial_{t}K_{s})-2\tau_{0}(\partial_{t}H,\partial_{t}K_{s})-\frac{2}{\tau_{0}}(K_{s},H).$$

(IV)  Using the Young’s inequality, we can bound the last three RHS terms of (2.28) as follows:

$$-(H,\partial_{t}K_{s})\leq\frac{\mu_{0}}{2\tau_{0}}\|H\|^{2}+\frac{\tau_{0}}{2\mu_{0}}\|\partial_{t}K_{s}\|^{2},$$ (2.29)

$$-2\tau_{0}(\partial_{t}H,\partial_{t}K_{s})\leq 2\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}}{2\mu_{0}}\|\partial_{t}K_{s}\|^{2},$$ (2.30)

$$-\frac{2}{\tau_{0}}(K_{s},H)\leq\frac{\mu_{0}}{2\tau_{0}}\|H\|^{2}+\frac{1}{2\tau_{0}\mu_{0}}\|K_{s}\|^{2}.$$ (2.31)

Substituting the above estimates (2.29)-(2.31) into (2.28), and dropping the last three non-negative terms on the left hand side of (2.28), we obtain

$$\frac{1}{2}\frac{d}{dt}\left(2\tau_{0}\epsilon_{0}\|\partial_{t}\mbox{\boldmath$E$}\|^{2}+\sigma_{0}\|\mbox{\boldmath$E$}\|^{2}_{\Gamma}+\mu_{0}\|\tau_{0}^{1/2}\partial_{t}H+\tau_{0}^{-1/2}H\|^{2}+\frac{\mu_{0}}{\tau_{0}}\|H\|^{2}\right.$$ (2.32)

$$\quad\left.+\frac{2\epsilon_{0}}{\tau_{0}}\|\mbox{\boldmath$E$}\|^{2}+\frac{2}{\sigma_{0}}\|\mbox{\boldmath$J$}\|^{2}_{\Gamma}+\tau_{0}\mu_{0}\|\partial_{t}H\|^{2}+\frac{\tau_{0}^{2}}{\sigma_{0}}\|\partial_{t}\mbox{\boldmath$J$}\|^{2}_{\Gamma}\right)$$

$$\leq \tau_{0}\mu_{0}(2+\tau_{0}^{-1})||\partial_{t}H||^{2}+\frac{\mu_{0}}{\tau_{0}}\|H\|^{2}+\frac{\tau_{0}}{\mu_{0}}\|\partial_{t}K_{s}\|^{2}+\frac{1}{2\tau_{0}\mu_{0}}\|K_{s}\|^{2}.$$

The proof is complete by applying the Gronwall inequality to (2.32).      

3 The leapfrog finite element scheme and its analysis

To design a finite element method, we partition the physical domain $\Omega$ with $\Gamma$ as an internal boundary by a shape regular triangular mesh $\mathcal{T}_{h}$ with maximum mesh size $h$. Without loss of generality, we consider the following Raviart-Thomas-Nédélec (RTN)’s mixed spaces $U_{h}$ and $\mbox{\boldmath$V$}_{h}$ on triangular elements [24, 30]: For any $r\geq 1$,

$$U_{h}=\{u_{h}\in L^{2}(\Omega):~{}u_{h}|_{K}\in p_{r-1},\forall K\in T_{h}\},$$

$$\mbox{\boldmath$V$}_{h}=\{\mbox{\boldmath$v$}_{h}\in H(curl;\Omega):~{}\mbox{\boldmath$v$}_{h}|_{K}\in(p_{r-1})^{2}\oplus S_{r},\forall K\in T_{h}\},~{}~{}S_{r}=\{\mbox{\boldmath$p$}\in\tilde{p}_{r}^{2},\mbox{\boldmath$x$}\cdot\mbox{\boldmath$p$}=0\},$$

where $p_{r}$ denotes the space of polynomials of degree less than or equal to $r$, and $\tilde{p}_{r}^{2}$ represents the space of homogeneous polynomials of degree $r$.

To handle the PEC boundary condition (2.4), we introduce the subspace

$$\mbox{\boldmath$V$}_{h}^{0}=\{\mbox{\boldmath$v$}_{h}\in\mbox{\boldmath$V$}_{h}:~{}~{}\bm{\hat{\nu}}\times\mbox{\boldmath$v$}_{h}=0\quad\mbox{on}~{}\partial\Omega\}.$$

To construct the fully discrete finite element scheme, we partition the time interval $[0,T]$ uniformly by points $t_{i}=i\tau,i=0,...,N_{t}$, where $\tau=\frac{T}{N_{t}}$ denotes the time step size.

Now we can construct the following leapfrog time stepping scheme: Given proper initial approximations of $\mbox{\boldmath$E$}_{h}^{-1},\mbox{\boldmath$E$}_{h}^{0}\in\mbox{\boldmath$V$}_{h}$, $H_{h}^{\frac{1}{2}}\in U_{h}$, for any $n\geq 0$, find $\mbox{\boldmath$E$}_{h}^{n+1}\in\mbox{\boldmath$V$}_{h}$, $H_{h}^{n+\frac{1}{2}}\in U_{h}$ such that

$$\epsilon_{0}(\delta_{\tau}^{2}\mbox{\boldmath$E$}_{h}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\frac{\epsilon_{0}}{\tau_{0}}(\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\frac{1}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$}^{n}_{h},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})$$

$$\hskip 72.26999pt=\frac{1}{\tau_{0}}(\overline{H}^{n}_{h},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})-\frac{\sigma_{0}}{\tau_{0}}\langle\mbox{\boldmath$E$}_{h}^{n},\mbox{\boldmath${\bf\phi}$}_{h}\rangle_{\Gamma}-\frac{1}{\mu_{0}}(K_{s}^{n},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h}),$$ (3.1)

$$\mu_{0}(\delta_{\tau}H^{n}_{h},\psi_{h})=-(\nabla\times\mbox{\boldmath$E$}_{h}^{n},\psi_{h})-(K_{s}^{n},\psi_{h}),$$ (3.2)

hold true for any test functions $\mbox{\boldmath${\bf\phi}$}_{h}\in\mbox{\boldmath$V$}_{h}^{0},\psi_{h}\in U_{h}$. Here we adopt the following central difference operators and averaging operator in time: For any time sequence function $u^{n}$,

$$\delta_{\tau}u^{n}=\frac{u^{n+\frac{1}{2}}-u^{n-\frac{1}{2}}}{\tau},~{}\delta_{2\tau}u^{n}=\frac{u^{n+1}-u^{n-1}}{2\tau},~{}\delta_{\tau}^{2}u^{n}=\frac{u^{n+1}-2u^{n}+u^{n-1}}{\tau^{2}},~{}\overline{u}^{n}=\frac{u^{n+\frac{1}{2}}+u^{n-\frac{1}{2}}}{2}.$$

Corresponding to the finite element spaces $\mbox{\boldmath$V$}_{h}$ and $U_{h}$, we denote $\Pi_{c}$ and $\Pi_{2}$ for the standard Nédélec interpolation in space $\mbox{\boldmath$V$}_{h}$ and the standard $L^{2}$ projection onto space $U_{h}$, respectively. Furthermore, the following interpolation and projection errors hold true (cf. [24, 30]):

$$||\mbox{\boldmath$u$}-\Pi_{c}\mbox{\boldmath$u$}||_{H(curl;\Omega)}\leq ch^{r}||\mbox{\boldmath$u$}||_{H^{r}(curl;\Omega)},~{}~{}\forall~{}\mbox{\boldmath$u$}\in H^{r}(curl;\Omega),\quad r\geq 1,$$ (3.3)

$$||u-\Pi_{2}u||_{L^{2}(\Omega)}\leq ch^{r}||u||_{H^{r}(\Omega)},~{}~{}\forall~{}u\in H^{r}(\Omega),\quad r\geq 0,$$ (3.4)

where $||u||_{H^{r}(\Omega)}$ denotes the norm for the Sobolev space $H^{r}(\Omega)$, and $||\mbox{\boldmath$u$}||_{H^{r}(curl;\Omega)}:=(||\mbox{\boldmath$u$}||^{2}_{(H^{r}(\Omega))^{2}}+||\nabla\times\mbox{\boldmath$u$}||^{2}_{H^{r}(\Omega)})^{1/2}$ is the norm for the Sobolev space

$$H^{r}(curl;\Omega)=\{\mbox{\boldmath$u$}\in(H^{r}(\Omega))^{2}:~{}\nabla\times\mbox{\boldmath$u$}\in H^{r}(\Omega)\}.$$

The initial conditions (2.15) are discretized as follows:

$$\mbox{\boldmath$E$}_{h}^{0}=\Pi_{c}\mbox{\boldmath$E$}_{0}(\mbox{\boldmath$x$}),$$ (3.5)

$$H_{h}^{-\frac{1}{2}}=\Pi_{2}\left[H(\cdot,0)-\frac{\tau}{2}\partial_{t}H(\cdot,0)\right]=\Pi_{2}\left[H_{0}(\mbox{\boldmath$x$})+\frac{\tau}{2\mu_{0}}(\nabla\times\mbox{\boldmath$E$}_{0}(\mbox{\boldmath$x$})+K_{s}(\mbox{\boldmath$x$},0))\right],$$ (3.6)

$$\frac{\mbox{\boldmath$E$}_{h}^{1}-\mbox{\boldmath$E$}_{h}^{-1}}{2\tau}=\Pi_{c}\left(\epsilon_{0}^{-1}\nabla\times H_{0}(\mbox{\boldmath$x$})\right),$$ (3.7)

where we used the Taylor expansion and the governing equation (2.2).

The implementation of the scheme (3.1)-(3.2) is quite simple. At each time step, we first solve (3.2) for $H_{h}^{n+\frac{1}{2}}$; then solve (3.1) for $\mbox{\boldmath$E$}_{h}^{n+1}$. Of course, at the first time step when $n=0$, we need to use the initial conditions (3.5)-(3.7). It can be seen that the coefficient matrix at each time step is symmetric and positive definite. Hence, the solvability of our scheme (3.1)-(3.2) is guaranteed.

3.1 The stability analysis

To prove the discrete stability for our scheme (3.1)-(3.2), we denote $C_{v}=\frac{1}{\sqrt{\epsilon_{0}\mu_{0}}}$ for the wave propagation speed, and recall the standard inverse estimate:

$$||\nabla\times\mbox{\boldmath${\bf\phi}$}_{h}||\leq C_{\mathrm{in}}h^{-1}||\mbox{\boldmath${\bf\phi}$}_{h}||,\quad\forall~{}\mbox{\boldmath${\bf\phi}$}_{h}\in\mbox{\boldmath$V$}_{h},$$ (3.8)

and the trace inequality:

$$||\mbox{\boldmath${\bf\phi}$}_{h}||_{\Gamma}\leq C_{\mathrm{tr}}h^{-1/2}||\mbox{\boldmath${\bf\phi}$}_{h}||,\quad\forall~{}\mbox{\boldmath${\bf\phi}$}_{h}\in\mbox{\boldmath$V$}_{h},$$ (3.9)

where the positive constants $C_{\mathrm{in}}$ and $C_{\mathrm{tr}}$ are independent of the mesh size $h$.

Theorem 3.1.

Denote the discrete energy for the scheme (3.1)-(3.2):

$$ENG_{m} := \epsilon_{0}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}+\frac{1}{2\mu_{0}}(||\nabla\times\mbox{\boldmath$E$}^{m+1}_{h}||^{2}+||\nabla\times\mbox{\boldmath$E$}^{m}_{h}||^{2})+\mu_{0}||H_{h}^{m+\frac{1}{2}}||^{2}$$ (3.10)

$$+\frac{\sigma_{0}}{2\tau_{0}}(||\mbox{\boldmath$E$}_{h}^{m+1}||^{2}_{\Gamma}+||\mbox{\boldmath$E$}_{h}^{m}||^{2}_{\Gamma})+\frac{\tau}{4\mu_{0}\tau_{0}}(||\nabla\times\mbox{\boldmath$E$}^{m+1}_{h}||^{2}+||\nabla\times\mbox{\boldmath$E$}^{m}_{h}||^{2}).$$

Then under the time step constraint:

$$\tau\leq\min\left(1,\frac{h}{2C_{\mathrm{in}}C_{v}},\frac{h\sqrt{\epsilon_{0}\tau_{0}}}{2C_{\mathrm{tr}}\sqrt{\sigma_{0}}},\frac{h\sqrt{\tau_{0}}}{\sqrt{2}C_{\mathrm{in}}C_{v}},\frac{h\tau_{0}}{C_{\mathrm{in}}C_{v}}\right),$$ (3.11)

we have the following stability: For any $m\in[1,N_{t}-1]$,

$$ENG_{m}\leq C_{*}\left[ENG_{0}+||K_{s}^{0}||^{2}+||K_{s}^{1}||^{2}+\tau\sum_{n=1}^{m}(||K_{s}^{n}||^{2}+||\delta_{2\tau}K_{s}^{n}||^{2})\right],$$ (3.12)

where the positive constant $C_{*}$ is independent of $\tau$ and $h$.

Proof. To make our proof easy to follow, we divite it into several major steps.

(I)  Choosing $\mbox{\boldmath${\bf\phi}$}_{h}=\tau\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}$ in (3.1), and using the following identities

$$\tau(\delta_{\tau}^{2}u^{n},\delta_{2\tau}u^{n})=\frac{1}{2}(||\delta_{\tau}u^{n+\frac{1}{2}}||^{2}-||\delta_{\tau}u^{n-\frac{1}{2}}||^{2}),$$ (3.13)

$$\tau(u^{n},\delta_{2\tau}u^{n})=\frac{1}{4}(||u^{n+1}||^{2}-||u^{n-1}||^{2})-\frac{\tau^{2}}{4}(||\delta_{\tau}u^{n+\frac{1}{2}}||^{2}-||\delta_{\tau}u^{n-\frac{1}{2}}||^{2}),$$ (3.14)

with $u^{n}=\mbox{\boldmath$E$}_{h}^{n}$, we have

$$\frac{\epsilon_{0}}{2}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}-||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}}||^{2})+\frac{\tau\epsilon_{0}}{\tau_{0}}||\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}||^{2}$$

$$+\frac{1}{4\mu_{0}}\left[(||\nabla\times\mbox{\boldmath$E$}^{n+1}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{n-1}_{h}||^{2})-\tau^{2}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}}||^{2})\right]$$

$$+\frac{\sigma_{0}}{4\tau_{0}}\left[(||\mbox{\boldmath$E$}^{n+1}_{h}||^{2}_{\Gamma}-||\mbox{\boldmath$E$}^{n-1}_{h}||^{2}_{\Gamma})-\tau^{2}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}_{\Gamma}-||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}}||^{2}_{\Gamma})\right]$$

$$\hskip 72.26999pt=\frac{\tau}{\tau_{0}}(\overline{H}^{n}_{h},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})-\frac{\tau}{\mu_{0}}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}).$$ (3.15)

Taking $\psi_{h}=\frac{\tau^{2}}{2\mu_{0}\tau_{0}}\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}$ in (3.2), and using (3.14) with $u^{n}=\mbox{\boldmath$E$}_{h}^{n}$, we obtain

$$\frac{\tau}{2\tau_{0}}(H^{n+\frac{1}{2}}_{h}-H^{n-\frac{1}{2}}_{h},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})=-\frac{\tau^{2}}{2\mu_{0}\tau_{0}}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})$$ (3.16)

$$\quad-\frac{\tau}{8\mu_{0}\tau_{0}}\left[(||\nabla\times\mbox{\boldmath$E$}^{n+1}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{n-1}_{h}||^{2})-\tau^{2}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}}||^{2})\right].$$

Taking $\psi_{h}=\tau\overline{H}_{h}^{n}$ in (3.2), we obtain

$$\frac{\mu_{0}}{2}(||H^{n+\frac{1}{2}}_{h}||^{2}-||H^{n-\frac{1}{2}}_{h}||^{2})=-\tau(\nabla\times\mbox{\boldmath$E$}^{n}_{h},\overline{H}_{h}^{n})-\tau(K_{s}^{n},\overline{H}_{h}^{n}).$$ (3.17)

Adding (3.15), (3.16), and (3.17) together, then summing up the result from $n=1$ to any $m\leq N_{t}-1$, we have

$$\frac{\epsilon_{0}}{2}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}-||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2})+\frac{\tau\epsilon_{0}}{\tau_{0}}\sum_{n=1}^{m}||\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}||^{2}$$

$$+\frac{1}{4\mu_{0}}\left[(||\nabla\times\mbox{\boldmath$E$}^{m+1}_{h}||^{2}+||\nabla\times\mbox{\boldmath$E$}^{m}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{1}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{0}_{h}||^{2})\right.$$

$$\hskip 72.26999pt\left.-\tau^{2}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2})\right]$$

$$+\frac{\sigma_{0}}{4\tau_{0}}\left[(||\mbox{\boldmath$E$}^{m+1}_{h}||^{2}_{\Gamma}+||\mbox{\boldmath$E$}^{m}_{h}||^{2}_{\Gamma}-||\mbox{\boldmath$E$}^{1}_{h}||^{2}_{\Gamma}-||\mbox{\boldmath$E$}^{0}_{h}||^{2}_{\Gamma})-\tau^{2}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}_{\Gamma}-||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2}_{\Gamma})\right]$$

$$+\frac{\mu_{0}}{2}(||H^{m+\frac{1}{2}}_{h}||^{2}-||H^{\frac{1}{2}}_{h}||^{2})$$

$$+\frac{\tau}{8\mu_{0}\tau_{0}}\left[(||\nabla\times\mbox{\boldmath$E$}^{m+1}_{h}||^{2}+||\nabla\times\mbox{\boldmath$E$}^{m}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{1}_{h}||^{2}-||\nabla\times\mbox{\boldmath$E$}^{0}_{h}||^{2})\right.$$

$$\hskip 72.26999pt\left.-\tau^{2}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2})\right]$$

$$\quad=\frac{\tau}{\tau_{0}}\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})-\frac{\tau}{\mu_{0}}\sum_{n=1}^{m}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})-\frac{\tau^{2}}{2\mu_{0}\tau_{0}}\sum_{n=1}^{m}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})$$

$$\quad\quad-\tau\sum_{n=1}^{m}(\nabla\times\mbox{\boldmath$E$}^{n}_{h},\overline{H}_{h}^{n})-\tau\sum_{n=1}^{m}(K_{s}^{n},\overline{H}_{h}^{n}).$$ (3.18)

(II)  By the definition of the discrete energy $ENG_{m}$ and dropping the non-negative term $\frac{\tau\epsilon_{0}}{\tau_{0}}\sum_{n=1}^{m}||\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}||^{2}$ on the left hand side of (3.18), we can rewrite (3.18) as follows

$$ENG_{m}\leq ENG_{0}+\frac{\tau^{2}}{2\mu_{0}}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2})$$

$$\quad+\frac{\sigma_{0}\tau^{2}}{2\tau_{0}}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}_{\Gamma}-||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2}_{\Gamma})+\frac{\tau^{3}}{4\mu_{0}\tau_{0}}(||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}-||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{\frac{1}{2}}||^{2})$$

$$\quad+\frac{\tau}{\tau_{0}}\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})-\frac{\tau}{\mu_{0}}\sum_{n=1}^{m}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})-\frac{\tau^{2}}{2\mu_{0}\tau_{0}}\sum_{n=1}^{m}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})$$

$$\quad\quad-\tau\sum_{n=1}^{m}(\nabla\times\mbox{\boldmath$E$}^{n}_{h},\overline{H}_{h}^{n})-\tau\sum_{n=1}^{m}(K_{s}^{n},\overline{H}_{h}^{n})=:ENG_{0}+\sum_{k=1}^{8}Sta_{k}.$$ (3.19)

Now we just need to estimate each $Sta_{k}$ in (3.19). Using the inverse estimate (3.8), it is easy to see that

$$Sta_{1}\leq\frac{\tau^{2}}{2\mu_{0}}||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}\leq\frac{C_{\mathrm{in}}^{2}}{2\mu_{0}}\cdot\tau^{2}h^{-2}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2},$$ (3.20)

and

$$Sta_{3}\leq\frac{\tau}{4\mu_{0}\tau_{0}}\cdot\tau^{2}||\nabla\times\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}\leq\frac{\tau C_{\mathrm{in}}^{2}}{4\mu_{0}\tau_{0}}\cdot\tau^{2}h^{-2}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}.$$ (3.21)

Using the trace inequality (3.8), we have

$$Sta_{2}\leq\frac{\sigma_{0}\tau^{2}}{2\tau_{0}}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}_{\Gamma}\leq\frac{\sigma_{0}C_{\mathrm{tr}}^{2}}{2\tau_{0}}\cdot\tau^{2}h^{-1}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}.$$ (3.22)

Similarly by the inverse estimate (3.8) and the Cauchy-Schwarz inequality, we have

$$Sta_{6} = -\frac{\tau^{2}}{2\mu_{0}\tau_{0}}\sum_{n=1}^{m}(K_{s}^{n},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})\leq\frac{\tau}{2\mu_{0}\tau_{0}}\sum_{n=1}^{m}||K_{s}^{n}||\cdot C_{\mathrm{in}}\tau h^{-1}||\frac{1}{2}(\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}+\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}})||$$ (3.23)

$$\leq \tau\cdot\frac{C_{\mathrm{in}}\tau h^{-1}}{4\mu_{0}\tau_{0}}\sum_{n=1}^{m}\left[||K_{s}^{n}||^{2}+\frac{1}{2}(||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}+||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n-\frac{1}{2}}||^{2})\right].$$

It is also easy to obtain the following estimates:

$$Sta_{7} = -\tau\sum_{n=1}^{m}(\frac{1}{\sqrt{\mu_{0}}}\nabla\times\mbox{\boldmath$E$}^{n}_{h},\sqrt{\mu_{0}}\overline{H}_{h}^{n})$$ (3.24)

$$\leq \frac{\tau}{2}\sum_{n=1}^{m}\left[\frac{1}{\mu_{0}}||\nabla\times\mbox{\boldmath$E$}^{n}_{h}||^{2}+\frac{\mu_{0}}{2}(||H_{h}^{n+\frac{1}{2}}||^{2}+||H_{h}^{n-\frac{1}{2}}||^{2})\right],$$

and

$$Sta_{8}\leq\frac{\tau}{2}\sum_{n=1}^{m}\left[\frac{1}{\mu_{0}}||K_{s}^{n}||^{2}+\frac{\mu_{0}}{2}(||H_{h}^{n+\frac{1}{2}}||^{2}+||H_{h}^{n-\frac{1}{2}}||^{2})\right].$$ (3.25)

The estimates of $Sta_{4}$ and $Sta_{5}$ need a special treatment, since $\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n}$ can not be controlled by $ENG_{m}$. In the next major part, we will carry out the estimates of $Sta_{4}$ and $Sta_{5}$.

(III) Note that

$$\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n})=\frac{1}{2\tau}\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n+1})-\frac{1}{2\tau}\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})$$ (3.26)

$$= \frac{1}{2\tau}\sum_{n=3}^{m+2}(H^{n-\frac{5}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})-\frac{1}{2\tau}\sum_{n=1}^{m}(H^{n-\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})$$

$$= -\frac{1}{2}\sum_{n=3}^{m}(\delta_{\tau}H^{n-2}_{h}+\delta_{\tau}H^{n-1}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})+\frac{1}{2\tau}(H^{m-\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{m+1})$$

$$\quad+\frac{1}{2\tau}(H^{m-\frac{3}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{m})-\frac{1}{2\tau}(H^{\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{0})-\frac{1}{2\tau}(H^{\frac{3}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{1}).$$

Using (3.26) and the inequality $ab\leq\frac{1}{4\delta}a^{2}+\delta b^{2}$, we have

$$Sta_{4} = -\frac{\tau}{2\tau_{0}}\sum_{n=3}^{m}(\delta_{\tau}H^{n-2}_{h}+\delta_{\tau}H^{n-1}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})+\frac{1}{2\tau_{0}}(H^{m-\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{m+1})$$ (3.27)

$$\quad+\frac{1}{2\tau_{0}}(H^{m-\frac{3}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{m})-\frac{1}{2\tau_{0}}(H^{\frac{1}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{0})-\frac{1}{2\tau_{0}}(H^{\frac{3}{2}}_{h},\nabla\times\mbox{\boldmath$E$}_{h}^{1})$$

$$\leq \frac{\tau}{4\tau_{0}}\sum_{n=3}^{m}(\frac{\mu_{0}}{2}||\delta_{\tau}H^{n-2}_{h}||^{2}+\frac{\mu_{0}}{2}||\delta_{\tau}H^{n-1}_{h}||^{2}+\frac{1}{\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{n-1}||^{2})$$

$$\quad+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{m+1}||^{2}+\frac{\mu_{0}}{2\tau_{0}^{2}}||H^{m-\frac{1}{2}}_{h}||^{2}+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{m}||^{2}+\frac{\mu_{0}}{2\tau_{0}^{2}}||H^{m-\frac{3}{2}}_{h}||^{2}$$

$$\quad+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{0}||^{2}+\frac{\mu_{0}}{2\tau_{0}^{2}}||H^{\frac{1}{2}}_{h}||^{2}+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{1}||^{2}+\frac{\mu_{0}}{2\tau_{0}^{2}}||H^{\frac{3}{2}}_{h}||^{2}.$$

By the same technique, we have

$$Sta_{5} = -\frac{\tau}{\mu_{0}}\sum_{n=3}^{m}(\delta_{2\tau}K_{s}^{n-1},\nabla\times\mbox{\boldmath$E$}_{h}^{n-1})+\frac{1}{2\mu_{0}}(K_{s}^{m},\nabla\times\mbox{\boldmath$E$}_{h}^{m+1})$$ (3.28)

$$\quad+\frac{1}{2\mu_{0}}(K_{s}^{m-1},\nabla\times\mbox{\boldmath$E$}_{h}^{m})-\frac{1}{2\mu_{0}}(K_{s}^{1},\nabla\times\mbox{\boldmath$E$}_{h}^{0})-\frac{1}{2\mu_{0}}(K_{s}^{2},\nabla\times\mbox{\boldmath$E$}_{h}^{1})$$

$$\leq \frac{\tau}{2\mu_{0}}\sum_{n=3}^{m}(||\delta_{2\tau}K_{s}^{n-1}||^{2}+||\nabla\times\mbox{\boldmath$E$}_{h}^{n-1}||^{2})$$

$$\quad+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{m+1}||^{2}+\frac{1}{2\mu_{0}}||K_{s}^{m}||^{2}+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{m}||^{2}+\frac{1}{2\mu_{0}}||K_{s}^{m-1}||^{2}$$

$$\quad+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{0}||^{2}+\frac{1}{2\mu_{0}}||K_{s}^{1}||^{2}+\frac{1}{8\mu_{0}}||\nabla\times\mbox{\boldmath$E$}_{h}^{1}||^{2}+\frac{1}{2\mu_{0}}||K_{s}^{2}||^{2}.$$

Substituting all the above estimates of $Sta_{k}$ into (3.19) and collecting like terms, we obtain

$$ENG_{m}\leq ENG_{0}+\left[\frac{C_{\mathrm{in}}^{2}\tau^{2}h^{-2}}{2\mu_{0}}+\frac{\sigma_{0}C_{\mathrm{tr}}^{2}\tau^{2}h^{-2}}{2\tau_{0}}+\frac{\tau C_{\mathrm{in}}^{2}\tau^{2}h^{-2}}{4\mu_{0}\tau_{0}}+\frac{C_{\mathrm{in}}\tau^{2}h^{-1}}{8\mu_{0}\tau_{0}}\right]||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{m+\frac{1}{2}}||^{2}$$

$$\quad+\tau(\frac{C_{\mathrm{in}}\tau h^{-1}}{4\mu_{0}\tau_{0}}+\frac{1}{2\mu_{0}})\sum_{n=1}^{m}||K_{s}^{n}||^{2}+\tau\cdot\frac{C_{\mathrm{in}}\tau h^{-1}}{4\mu_{0}\tau_{0}}\sum_{n=0}^{m-1}||\delta_{\tau}\mbox{\boldmath$E$}_{h}^{n+\frac{1}{2}}||^{2}$$

$$\quad+\tau(\frac{1}{2\mu_{0}}+\frac{1}{4\mu_{0}\tau_{0}}+\frac{1}{2\mu_{0}})\sum_{n=1}^{m}||\nabla\times\mbox{\boldmath$E$}_{h}^{n}||^{2}+\frac{\tau\mu_{0}}{2}||H_{h}^{m+\frac{1}{2}}||^{2}+\tau\mu_{0}\sum_{n=0}^{m-1}||H_{h}^{n+\frac{1}{2}}||^{2}$$

$$\quad+\frac{\tau\mu_{0}}{4\tau_{0}}\sum_{n=1}^{m-1}||\delta_{\tau}H_{h}^{n}||^{2}+\frac{\mu_{0}}{2\tau_{0}^{2}}(||H_{h}^{m-\frac{1}{2}}||^{2}+||H_{h}^{m-\frac{3}{2}}||^{2})+\frac{\mu_{0}}{2\tau_{0}^{2}}(||H_{h}^{\frac{1}{2}}||^{2}+||H_{h}^{\frac{3}{2}}||^{2})$$

$$\quad+\frac{1}{4\mu_{0}}(||\nabla\times\mbox{\boldmath$E$}_{h}^{m+1}||^{2}+||\nabla\times\mbox{\boldmath$E$}_{h}^{m}||^{2})+\frac{1}{4\mu_{0}}(||\nabla\times\mbox{\boldmath$E$}_{h}^{1}||^{2}+||\nabla\times\mbox{\boldmath$E$}_{h}^{0}||^{2})$$

$$\quad+\frac{\tau}{2\mu_{0}}\sum_{n=2}^{m-1}||\delta_{2\tau}K_{s}^{n}||^{2}+\frac{1}{2\mu_{0}}(||K_{s}^{m}||^{2}+||K_{s}^{m-1}||^{2}+||K_{s}^{1}||^{2}+||K_{s}^{2}||^{2}).$$ (3.29)

(IV)  Now all the right hand side terms of (3.29), except $\frac{\tau\mu_{0}}{4\tau_{0}}\sum_{n=1}^{m-1}||\delta_{\tau}H_{h}^{n}||^{2}$ and $\frac{\mu_{0}}{2\tau_{0}^{2}}(||H_{h}^{m-\frac{1}{2}}||^{2}+||H_{h}^{m-\frac{3}{2}}||^{2})$, can be controlled by choosing $\tau$ small enough and using the discrete Gronwall inequality. To bound these two terms, squaring the following identity

$$H_{h}^{m-\frac{1}{2}} = (H_{h}^{m-\frac{1}{2}}-H_{h}^{m-\frac{3}{2}})+(H_{h}^{m-\frac{3}{2}}-H_{h}^{m-\frac{5}{2}})+\cdots+(H_{h}^{\frac{3}{2}}-H_{h}^{\frac{1}{2}})+H_{h}^{\frac{1}{2}}$$ (3.30)

$$= \tau\sum_{n=1}^{m-1}\delta_{\tau}H_{h}^{n}+H_{h}^{\frac{1}{2}},$$

we have

$$||H_{h}^{m-\frac{1}{2}}||^{2} \leq 2(\tau^{2}||\sum_{n=1}^{m-1}\delta_{\tau}H_{h}^{n}||^{2}+||H_{h}^{\frac{1}{2}}||^{2})\leq 2||H_{h}^{\frac{1}{2}}||^{2}+2\tau^{2}(\sum_{n=1}^{m-1}1^{2})(\sum_{n=1}^{m-1}||\delta_{\tau}H_{h}^{n}||^{2})$$ (3.31)

$$\leq 2||H_{h}^{\frac{1}{2}}||^{2}+2T\tau\sum_{n=1}^{m-1}||\delta_{\tau}H_{h}^{n}||^{2}.$$

By the same argument, we have

$$||H_{h}^{m-\frac{3}{2}}||^{2}=||\tau\sum_{n=1}^{m-2}\delta_{\tau}H_{h}^{n}+H_{h}^{\frac{1}{2}}||^{2}\leq 2||H_{h}^{\frac{1}{2}}||^{2}+2T\tau\sum_{n=1}^{m-2}||\delta_{\tau}H_{h}^{n}||^{2}.$$ (3.32)

To bound $\delta_{\tau}H^{n}_{h}$, taking $\psi_{h}=\delta_{\tau}H^{n}_{h}$ in (3.2), we obtain

$$\mu_{0}||\delta_{\tau}H^{n}_{h}||^{2} = -(\nabla\times E_{h}^{n},\delta_{\tau}H^{n}_{h})-(K_{s}^{n},\delta_{\tau}H^{n}_{h})$$ (3.33)

$$\leq (\frac{\mu_{0}}{4}||\delta_{\tau}H^{n}_{h}||^{2}+\frac{1}{\mu_{0}}||\nabla\times E_{h}^{n}||^{2})+(\frac{\mu_{0}}{4}||\delta_{\tau}H^{n}_{h}||^{2}+\frac{1}{\mu_{0}}||K_{s}^{n}||^{2}),$$

which leads to

$$\mu_{0}||\delta_{\tau}H^{n}_{h}||^{2}\leq\frac{2}{\mu_{0}}||\nabla\times E_{h}^{n}||^{2}+\frac{2}{\mu_{0}}||K_{s}^{n}||^{2}.$$ (3.34)

We can use similar techniques to bound the terms $K_{s}^{m}$ and $K_{s}^{m-1}$ in (3.29), even though they are given source functions and we can keep them as they are. When $m$ is even, we have

$$||K_{s}^{m}||^{2}=||2\tau\sum_{n=1}^{m-1}\delta_{2\tau}K_{s}^{n}+K_{s}^{0}||^{2}\leq 2||K_{s}^{0}||^{2}+4T\tau\sum_{n=1}^{m-1}||\delta_{2\tau}K_{s}^{n}||^{2}.$$ (3.35)

When $m$ is odd, we have

$$||K_{s}^{m}||^{2}=||2\tau\sum_{n=2}^{m-1}\delta_{2\tau}K_{s}^{n}+K_{s}^{1}||^{2}\leq 2||K_{s}^{1}||^{2}+4T\tau\sum_{n=2}^{m-1}||\delta_{2\tau}K_{s}^{n}||^{2}.$$ (3.36)

Substituting the estimates of (3.31), (3.32), and (3.34)-(3.36) into (3.29), and choosing $\tau$ small enough, e.g.,

$$\tau\leq 1,~{}~{}\frac{C_{\mathrm{in}}^{2}\tau^{2}h^{-2}}{2\mu_{0}}\leq\frac{\epsilon_{0}}{8}~{}(\mbox{or}~{}\tau\leq\frac{h}{2C_{\mathrm{in}}C_{v}}),~{}~{}\frac{\sigma_{0}C_{\mathrm{tr}}^{2}\tau^{2}h^{-2}}{2\tau_{0}}\leq\frac{\epsilon_{0}}{8}~{}(\mbox{or}~{}\tau\leq\frac{h\sqrt{\epsilon_{0}\tau_{0}}}{2C_{\mathrm{tr}}\sqrt{\sigma_{0}}}),$$

$$\frac{C_{\mathrm{in}}^{2}\tau^{2}h^{-2}}{4\mu_{0}\tau_{0}}\leq\frac{\epsilon_{0}}{8}~{}(\mbox{or}~{}\tau\leq\frac{h\sqrt{\tau_{0}}}{\sqrt{2}C_{\mathrm{in}}C_{v}}),~{}~{}\frac{C_{\mathrm{in}}\tau h^{-1}}{8\mu_{0}\tau_{0}}\leq\frac{\epsilon_{0}}{8}~{}(\mbox{or}~{}\tau\leq\frac{h\tau_{0}}{C_{\mathrm{in}}C_{v}}),$$ (3.37)

which is equivalent to (3.11), we complete the proof by the discrete Gronwall inequality.      

3.2 The error estimate

To carry out the error estimate of our scheme (3.1)-(3.2), we split the solution error for $E$ as follows:

$$\mbox{\boldmath$E$}(\mbox{\boldmath$x$},t_{n})-\mbox{\boldmath$E$}_{h}^{n}=(\mbox{\boldmath$E$}(\mbox{\boldmath$x$},t_{n})-\Pi_{c}\mbox{\boldmath$E$}^{n})-(\mbox{\boldmath$E$}_{h}^{n}-\Pi_{c}\mbox{\boldmath$E$}^{n})=:\mbox{\boldmath$E$}^{n}_{\xi}-\mbox{\boldmath$E$}^{n}_{\eta},$$ (3.38)

where for simplicity we denote $\mbox{\boldmath$E$}^{n}:=\mbox{\boldmath$E$}(\mbox{\boldmath$x$},t_{n})$. Similarly, we can define the solution error for $H$ as follows:

$$H(\mbox{\boldmath$x$},t_{n})-H_{h}^{n}=(H(\mbox{\boldmath$x$},t_{n})-\Pi_{2}H^{n})-(H_{h}^{n}-\Pi_{2}H^{n})=:H^{n}_{\xi}-H^{n}_{\eta}.$$ (3.39)

Integrating (2.13) with respect to $t$ from $t_{n-1}$ to $t_{n+1}$, then divide the result by $2\tau_{0}\tau$, we obtain

$$\epsilon_{0}(\delta_{2\tau}\partial_{t}\mbox{\boldmath$E$}^{n},\mbox{\boldmath${\bf\phi}$})+\frac{\epsilon_{0}}{\tau_{0}}(\delta_{2\tau}\mbox{\boldmath$E$}^{n},\mbox{\boldmath${\bf\phi}$})+\frac{1}{\mu_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}\nabla\times\mbox{\boldmath$E$}~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$})$$

$$= \frac{1}{\tau_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}H~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$})-\frac{1}{\mu_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}K_{s}~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$})-\frac{\sigma_{0}}{\tau_{0}}\langle\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}\mbox{\boldmath$E$}~{}dt,\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}.$$

Integrating (2.14) with respect to $t$ from $t_{n-\frac{1}{2}}$ to $t_{n+\frac{1}{2}}$, then divie the result by $\tau$, we have

$$\mu_{0}(\delta_{\tau}H^{n},\psi)=-(\frac{1}{\tau}\int_{t_{n-\frac{1}{2}}}^{t_{n+\frac{1}{2}}}\nabla\times\mbox{\boldmath$E$}~{}dt,\psi)-(\frac{1}{\tau}\int_{t_{n-\frac{1}{2}}}^{t_{n+\frac{1}{2}}}K_{s}~{}dt,\psi).$$ (3.41)

Now subtracting (3.1) from (3.2) with $\mbox{\boldmath${\bf\phi}$}=\mbox{\boldmath${\bf\phi}$}_{h}$, (3.2) from (3.41) with $\psi=\psi_{h}$, using the error notation we introduced, and after some lengthy but straightforward algebra, we can obtain the error equations:

$$\epsilon_{0}(\delta_{\tau}^{2}\mbox{\boldmath$E$}_{\eta}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\frac{\epsilon_{0}}{\tau_{0}}(\delta_{2\tau}\mbox{\boldmath$E$}_{\eta}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\frac{1}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$}^{n}_{\eta},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})$$ (3.42)

$$= \frac{1}{\tau_{0}}(\overline{H}^{n}_{\eta},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})-\frac{\sigma_{0}}{\tau_{0}}\langle\mbox{\boldmath$E$}_{\eta}^{n},\mbox{\boldmath${\bf\phi}$}_{h}\rangle_{\Gamma}+\epsilon_{0}(\delta_{\tau}^{2}\mbox{\boldmath$E$}_{\xi}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\epsilon_{0}(\delta_{2\tau}\partial_{t}\mbox{\boldmath$E$}^{n}-\delta_{\tau}^{2}\mbox{\boldmath$E$}^{n},\mbox{\boldmath${\bf\phi}$}_{h})$$

$$\quad+\frac{\epsilon_{0}}{\tau_{0}}(\delta_{2\tau}\mbox{\boldmath$E$}_{\xi}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\frac{1}{\mu_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}\nabla\times(\mbox{\boldmath$E$}-\Pi_{c}\mbox{\boldmath$E$}^{n})~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})$$

$$\quad+\frac{1}{\tau_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}(H-\Pi_{2}\overline{H}^{n})~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})+\frac{\sigma_{0}}{\tau_{0}}\langle\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}(\mbox{\boldmath$E$}-\Pi_{c}\mbox{\boldmath$E$}^{n})~{}dt,\mbox{\boldmath${\bf\phi}$}_{h}\rangle_{\Gamma}$$

$$\quad-\frac{1}{\mu_{0}}(\frac{1}{2\tau}\int_{t_{n-1}}^{t_{n+1}}(K_{s}-K_{s}^{n})~{}dt,\nabla\times\mbox{\boldmath${\bf\phi}$}_{h}),$$

and

$$\mu_{0}(\delta_{\tau}H^{n}_{\eta},\psi_{h})+(\nabla\times\mbox{\boldmath$E$}_{\eta}^{n},\psi_{h})=\mu_{0}(\delta_{\tau}H^{n}_{\xi},\psi_{h})$$

$$\quad-(\frac{1}{\tau}\int_{t_{n-\frac{1}{2}}}^{t_{n+\frac{1}{2}}}\nabla\times(\Pi_{c}\mbox{\boldmath$E$}^{n}-\mbox{\boldmath$E$})~{}dt,\psi_{h})-(\frac{1}{\tau}\int_{t_{n-\frac{1}{2}}}^{t_{n+\frac{1}{2}}}(K_{s}^{n}-K_{s})~{}dt,\psi_{h}).$$ (3.43)

Note that the error equation (3.42) has the first five terms in the same form as scheme (3.1), and the rest terms are extra error terms due to the time approximation, and projection or interpolation. Furthermore, the error equation (3.43) has the first two terms in the same form as scheme (3.2), and the rest three terms are error terms due to the time approximation, and projection or interpolation. We want to remark that those extra error terms are $O(\tau^{2}+h^{r})$ by the interpolation error estimate (3.3) and the projection error estimate (3.4). Following the stability proof (due to its technicality, we skip the proof details), we have the following error estimate for our scheme (3.1)-(3.2):

$$||\nabla\times(\mbox{\boldmath$E$}^{m}_{h}-\mbox{\boldmath$E$}^{m})||+||H_{h}^{m+\frac{1}{2}}-H^{m+\frac{1}{2}}||\leq C(\tau^{2}+h^{r}),$$ (3.44)

where $r$ is the degree of our finite element spaces $\mbox{\boldmath$V$}_{h}$ and $U_{h}$.

4 Numerical results

In this section, we present some numerical results solved by our proposed numerical scheme. The first example is developed to test the convergence of our scheme by using FEniCS [28], and the rest are carried out using NGSolve [36] to demonstrate that our reformulated graphene model can still generate surface plasmon polaritons.

4.1 Test of convergence rates

The first example is developed to test the convergence rate of our numerical scheme by a manufactured exact solution given as follows (adapted from our previous work [26]):

$$\mathbf{E}(x,y,t)=\begin{pmatrix}E_{x}\\ E_{y}\\ \end{pmatrix}=\begin{pmatrix}\sin(2\pi x)\sin(2\pi y)\sin(2\pi t)\\ \cos(2\pi x)\cos(2\pi y)\sin(2\pi t)\\ \end{pmatrix},$$

$$H_{1}(x,y,t)=\frac{1}{1+4\pi^{2}}\sin(2\pi x)\sin(2\pi y)\sin(2\pi t),$$

$$H_{2}(x,y,t)=\frac{1}{1+4\pi^{2}}\sin(2\pi x)\sin(2\pi y)(2\pi\cos(2\pi t)-2\pi\exp(-t)),$$

which satisfies the following weak form for the graphene model equations: For any $\mbox{\boldmath${\bf\phi}$}\in H_{0}(curl;\Omega)$ and $\psi\in L^{2}(\Omega)$,

$$\tau_{0}\epsilon_{0}(\partial_{tt}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\epsilon_{0}(\partial_{t}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\frac{\tau_{0}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$},\nabla\times\mbox{\boldmath${\bf\phi}$})=(H_{1},\nabla\times\mbox{\boldmath${\bf\phi}$})-\sigma_{0}\langle\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}$$

$$\hskip 108.405pt+(\mbox{\boldmath$f$}_{1},\mbox{\boldmath${\bf\phi}$})+(\tau_{0}\partial_{t}\mbox{\boldmath$f$}_{1},\mbox{\boldmath${\bf\phi}$})+\frac{\tau_{0}}{\mu_{0}}(f_{2},\nabla\times\mbox{\boldmath${\bf\phi}$}),\quad\mbox{in }\Omega_{1},$$ (4.1)

$$\mu_{0}(\partial_{t}H_{1},\psi)=(-\nabla\times\mbox{\boldmath$E$}+f_{2},\psi),\quad\mbox{in }\Omega_{1},$$ (4.2)

$$\tau_{0}\epsilon_{0}(\partial_{tt}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\epsilon_{0}(\partial_{t}\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$})+\frac{\tau_{0}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$},\nabla\times\mbox{\boldmath${\bf\phi}$})=(H_{2},\nabla\times\mbox{\boldmath${\bf\phi}$})-\sigma_{0}\langle\mbox{\boldmath$E$},\mbox{\boldmath${\bf\phi}$}\rangle_{\Gamma}$$

$$\hskip 108.405pt+(\mbox{\boldmath$f$}_{3},\mbox{\boldmath${\bf\phi}$})+(\tau_{0}\partial_{t}\mbox{\boldmath$f$}_{3},\mbox{\boldmath${\bf\phi}$})+\frac{\tau_{0}}{\mu_{0}}(f_{4},\nabla\times\mbox{\boldmath${\bf\phi}$}),\quad\mbox{in }\Omega_{2},$$ (4.3)

$$\mu_{0}(\partial_{t}H_{2},\psi)=(-\nabla\times\mbox{\boldmath$E$}+f_{4},\psi),\quad\mbox{in }\Omega_{2}.$$ (4.4)

Here the extra source terms $\mbox{\boldmath$f$}_{1},f_{2},\mbox{\boldmath$f$}_{3},f_{4}$ can be derived from the given exact solution $\mbox{\boldmath$E$},\mbox{\boldmath$J$},H_{1}$, and $H_{2}$ as in [26]. Note that the weak form is derived from [26, (4.1)-(4.5)] by following the same steps to get our weak formulation (2.13)-(2.14) with added source terms.

For simplicity, we choose the physical domain $\Omega=(0,1)^{2}$, which is split into two subdomains $\Omega_{1}=(0,1)\times(0.5,1)$ and $\Omega_{2}=(0,1)\times(0,0.5)$ with interface $\Gamma=\{y=0.5,x\in[0,1]\}$. We apply our developed scheme (3.1)-(3.2) to solve (4.1)-(4.4) with physical parameters $\epsilon_{0}=\mu_{0}=\tau_{0}=\sigma_{0}=1$.

Here the added source terms $\mbox{\boldmath$f$}_{1},f_{2}$, and $\mbox{\boldmath$f$}_{3}$ can be calculated from the given exact solution $\mbox{\boldmath$E$},H_{1}$, and $H_{2}$. We use our scheme (3.1)-(3.2) on the the same parameters and physical domain setup as [26]. We solve this example with a fixed small time step size $\tau=1\times 10^{-4}$, and various mesh sizes for $N_{t}=1000$ time steps. The obtained $L^{2}$ errors are presented in Table 1 and Table 2 for the RTN finite element spaces $U_{h}$ and $\mbox{\boldmath$V$}_{h}$ with polynomial degree $r=1,2$, respectively.

Table 1: The obtained errors for $N_{t}=1000,\tau=1\times 10^{-4},r=1$.

h	$\|\mbox{\boldmath$E$}-\mbox{\boldmath$E$}_{h}\|_{L^{2}(\Omega)}$	rate	$\|H-H_{h}\|_{L^{2}(\Omega)}$	rate
1/4	$1.995472\times 10^{-2}$		$4.037259\times 10^{-4}$
1/8	$1.005105\times 10^{-3}$	0.989385	$2.044797\times 10^{-4}$	0.981418
1/16	$5.035177\times 10^{-3}$	0.997231	$1.014174\times 10^{-4}$	1.011653
1/32	$2.518828\times 10^{-3}$	0.999290	$4.906353\times 10^{-5}$	1.047582
1/64	$1.259572\times 10^{-3}$	0.999819	$2.327676\times 10^{-5}$	1.075761
1/128	$6.199001\times 10^{-4}$	1.022826	$1.082733\times 10^{-5}$	1.104213

Table 2: The obtained errors for $N_{t}=1000,\tau=1\times 10^{-4},r=2$.

h	$\|\mbox{\boldmath$E$}-\mbox{\boldmath$E$}_{h}\|_{L^{2}(\Omega)}$	rate	$\|H-H_{h}\|_{L^{2}(\Omega)}$	rate
1/4	$5.028653\times 10^{-2}$		$2.172637\times 10^{-4}$
1/8	$1.290415\times 10^{-3}$	1.962337	$1.035932\times 10^{-4}$	1.068517
1/16	$3.186606\times 10^{-3}$	2.017743	$4.847654\times 10^{-4}$	1.095571
1/32	$7.545130\times 10^{-3}$	2.078403	$1.887592\times 10^{-5}$	1.360740
1/64	$2.112674\times 10^{-3}$	1.836476	$4.771321\times 10^{-6}$	1.984086
1/128	$5.979253\times 10^{-4}$	1.821033	$1.182453\times 10^{-6}$	2.012605

Then we test the convergence rate for $\tau$ by fixing $\tau=\frac{h}{200}$ to satisfy the stability condition. In Tables 3-4 the obtained $L^{2}$ errors for r=1,2 are presented.

Table 3: The obtained errors for $T=0.01,\tau=\frac{h}{200},r=1$.

h	$\|\mbox{\boldmath$E$}-\mbox{\boldmath$E$}_{h}\|_{L^{2}(\Omega)}$	rate	$\|H-H_{h}\|_{L^{2}(\Omega)}$	rate
1/10	$8.441669\times 10^{-3}$		$1.768658\times 10^{-4}$
1/20	$4.125569\times 10^{-3}$	1.032935	$8.359509\times 10^{-5}$	1.081165
1/40	$2.038308\times 10^{-3}$	1.017221	$3.914783\times 10^{-5}$	1.094486
1/80	$1.012953\times 10^{-3}$	1.008805	$1.868672\times 10^{-5}$	1.066919
1/160	$5.069774\times 10^{-4}$	0.998574	$1.035698\times 10^{-5}$	0.851410

Table 4: The obtained errors for $T=0.01,\tau=\frac{h}{200},r=2$.

h	$\|\mbox{\boldmath$E$}-\mbox{\boldmath$E$}_{h}\|_{L^{2}(\Omega)}$	rate	$\|H-H_{h}\|_{L^{2}(\Omega)}$	rate
1/10	$5.648216\times 10^{-3}$		$2.931658\times 10^{-4}$
1/20	$1.368327\times 10^{-3}$	2.045382	$1.216507\times 10^{-4}$	1.268972
1/40	$3.277738\times 10^{-4}$	2.061641	$5.237098\times 10^{-5}$	1.215905
1/80	$7.647986\times 10^{-5}$	2.099549	$1.938811\times 10^{-5}$	1.433595
1/160	$2.140379\times 10^{-5}$	1.837214	$4.807332\times 10^{-6}$	2.011864

4.2 Simulation of surface plasmon polaritons on the graphene sheets

In this section, we provide several numerical examples to illustrate the effectiveness of our graphene model in simulating the propagation of surface plasmon polaritons (SPPs) on graphene sheets.

To effectively demonstrate the SPPs propagating on graphene sheets, we use a perfectly matched layer (PML) to surround the physical domain $\Omega$. Since our current graphene model only involves the electric and magnetic fields, we just adopt the 2-D TEz Berenger’s PML model [2], which can be written as (cf. [25, (12)-(15)]): For any $(\mbox{\boldmath$x$},t)\in\Omega_{\mathrm{pml}}\times(0,T]$,

$$\epsilon_{0}\partial_{t}E_{x}+\sigma_{y}E_{x}=\partial_{y}(H_{zx}+H_{zy}),$$ (4.5)

$$\epsilon_{0}\partial_{t}E_{y}+\sigma_{x}E_{y}=-\partial_{x}(H_{zx}+H_{zy}),$$ (4.6)

$$\mu_{0}\partial_{t}H_{zx}+\frac{\mu_{0}}{\epsilon_{0}}\sigma_{x}H_{zx}=-\partial_{x}E_{y},$$ (4.7)

$$\mu_{0}\partial_{t}H_{zy}+\frac{\mu_{0}}{\epsilon_{0}}\sigma_{y}H_{zy}=\partial_{y}E_{x},$$ (4.8)

where $\sigma_{x}(x)$ and $\sigma_{y}(y)$ are the nonnegative damping functions in the $x$ and $y$ directions, respectively, and $\Omega_{\mathrm{pml}}$ represents the PML region. Here $H_{zx}$ and $H_{zy}$ are the two splitted components of the orginal magnetic field $H_{z}$, i.e., $H_{z}:=H_{zx}+H_{zy}$.

We propose the following finite element scheme for the above PML model in $\Omega_{\mathrm{pml}}$: For any $n\geq 0$, find $\mbox{\boldmath$E$}_{h}^{n+1}\in\mbox{\boldmath$V$}_{h}^{0},H_{zx}^{n+\frac{1}{2}},H_{zy}^{n+\frac{1}{2}}\in U_{h}$ such that

$$\epsilon_{0}(\delta_{2\tau}\mbox{\boldmath$E$}_{h}^{n},\mbox{\boldmath${\bf\phi}$}_{h})+\epsilon_{0}(D_{1}\mbox{\boldmath$E$}_{h}^{n},\mbox{\boldmath${\bf\phi}$}_{h})=(\overline{H}^{n}_{zx,h},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})+(\overline{H}^{n}_{zy,h},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h}),$$ (4.9)

$$\mu_{0}(\delta_{\tau}H^{n}_{zx,h},\psi_{h})+\frac{\mu_{0}}{\epsilon_{0}}(\sigma_{x}\overline{H}_{zx,h}^{n},\psi_{h})=-(\partial_{x}E_{y,h}^{n},\psi_{h}),$$ (4.10)

$$\mu_{0}(\delta_{\tau}H^{n}_{zy,h},\varphi_{h})+\frac{\mu_{0}}{\epsilon_{0}}(\sigma_{y}\overline{H}_{zy,h}^{n},\varphi_{h})=(\partial_{x}E_{x,h}^{n},\varphi_{h}),$$ (4.11)

hold true for any test functions $\mbox{\boldmath${\bf\phi}$}_{h}\in\mbox{\boldmath$V$}_{h}^{0}$, $\psi_{h}$, $\varphi_{h}\in U_{h}$, where $D_{1}=\mathrm{diag}(\sigma_{y},\sigma_{x})$.

To simplify the implementation, we merge the graphene scheme (3.1)-(3.2) and the PML scheme (4.9)-(4.11) together by using subdomain dependent coefficients and rewrite them as follows:

$$\left((\frac{\epsilon_{0}}{\tau^{2}}I+\frac{D_{1}}{2\tau}+\frac{C_{1}\epsilon_{0}}{2\tau\tau_{0}}I)\mbox{\boldmath$E$}^{n+1}_{h},\mbox{\boldmath${\bf\phi}$}_{h}\right)=\left(\frac{2\epsilon_{0}}{\tau^{2}}I\mbox{\boldmath$E$}^{n}_{h},\mbox{\boldmath${\bf\phi}$}_{h}\right)-\left((\frac{\epsilon_{0}}{\tau^{2}}I-\frac{D_{1}}{2\tau}-\frac{C_{1}\epsilon_{0}}{2\tau\tau_{0}}I)\mbox{\boldmath$E$}^{n+1}_{h},\mbox{\boldmath${\bf\phi}$}_{h}\right)$$

$$\hskip 144.54pt-\frac{C_{1}}{\mu_{0}}(\nabla\times\mbox{\boldmath$E$}^{n}_{h},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})-\frac{C_{1}}{\mu_{0}}(K_{sh}^{n},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})-\frac{\sigma_{0}}{\tau_{0}}\langle\mbox{\boldmath$E$}^{n}_{h},\mbox{\boldmath${\bf\phi}$}_{h}\rangle_{\Gamma}$$

$$\hskip 144.54pt+\frac{C_{1}}{2\tau_{0}}(H_{zx,h}^{n+\frac{1}{2}}+H_{zx,h}^{n-\frac{1}{2}},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})+\frac{C_{1}}{2\tau_{0}}(H_{zy,h}^{n+\frac{1}{2}}+H_{zy,h}^{n-\frac{1}{2}},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})$$

$$\hskip 144.54pt+\frac{C_{2}}{\tau}(H_{zx,h}^{n+\frac{1}{2}}-H_{zx,h}^{n-\frac{1}{2}},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h})+\frac{C_{2}}{\tau}(H_{zy,h}^{n+\frac{1}{2}}-H_{zy,h}^{n-\frac{1}{2}},\nabla\times\mbox{\boldmath${\bf\phi}$}_{h}),$$ (4.12)

$$\left((\frac{\mu_{0}}{\tau}+\frac{\mu_{0}\sigma_{x}}{2\epsilon_{0}})H^{n+\frac{1}{2}}_{zx,h},\psi_{h}\right)=\left((\frac{\mu_{0}}{\tau}-\frac{\mu_{0}\sigma_{x}}{2\epsilon_{0}})H^{n-\frac{1}{2}}_{zx,h},\psi_{h}\right)-(\partial_{x}E_{y,h}^{n},\psi_{h})-\frac{1}{2}(K_{s}^{n},\psi_{h}),$$ (4.13)

$$\left((\frac{\mu_{0}}{\tau}+\frac{\mu_{0}\sigma_{y}}{2\epsilon_{0}})H^{n+\frac{1}{2}}_{zy,h},\psi_{h}\right)=\left((\frac{\mu_{0}}{\tau}-\frac{\mu_{0}\sigma_{x}}{2\epsilon_{0}})H^{n-\frac{1}{2}}_{zy,h},\psi_{h}\right)+(\partial_{y}E_{x,h}^{n},\psi_{h})-\frac{1}{2}(K_{s}^{n},\psi_{h}),$$ (4.14)

where we denote the identity matrix $I=\mathrm{diag}(1,1)$, write $H_{h}=H_{zx,h}+H_{zy,h}$, and use the subdomain identity functions

$$C_{1}=\begin{cases}1,\quad\mbox{if}~{}\mbox{\boldmath$x$}\in\Omega\\ 0,\quad\mbox{if}~{}\mbox{\boldmath$x$}\in\Omega_{\mathrm{pml}}\end{cases},C_{2}=\begin{cases}0,\quad\mbox{if}~{}\mbox{\boldmath$x$}\in\Omega\\ 1,\quad\mbox{if}~{}\mbox{\boldmath$x$}\in\Omega_{\mathrm{pml}}\end{cases}.$$ (4.15)

The damping functions $\sigma_{x}$ and $\sigma_{y}$ for the PML are chosen as a fourth-order polynomial:

$$\sigma_{x}(x)=\begin{cases}\sigma_{\max}(\frac{x-b_{x}}{dd})^{4},&\mbox{when}~{}x\geq b_{x},\\ \sigma_{\max}(\frac{x-a_{x}}{dd})^{4},&x\leq a_{x},\\ 0,&\text{elsewhere},\end{cases}$$

where the coefficient $\sigma_{\max}=-\log(\mathrm{err})\cdot 5/(2\cdot\mathrm{dd}\cdot\eta)$ with $\mathrm{err}=10^{-7}$, $\eta=377$ is the impedance of free space, and $dd$ denotes the thickness of the PML in the $x$ direction. The function $\sigma_{y}$ has the same form but varies with respect to the $y$ variable.

In the rest of our simulations, we choose a physical domain $\Omega=[a_{x},b_{x}]~{}\mu\mathrm{m}\times[a_{y},b_{y}]~{}\mu\mathrm{m}$, which is partitioned by a regular unstructured triangular mesh and is surrounded by the Berenger’s PML with thickness $12h_{x}$ and $12h_{y}$ in the $x$ and $y$ directions, respectively. We denote $h_{x}$ and $h_{y}$ for the maximum mesh sizes in the $x$ and $y$ directions. The physical domain for the first three examples is a rectangle domain $\Omega=[-30,30]~{}\mu\mathrm{m}\times[-10,10]~{}\mu\mathrm{m}$ with $h_{x}=0.6~{}\mu\mathrm{m},h_{y}=0.2~{}\mu\mathrm{m}$, and the domain for the last two examples is $\Omega=[-20,20]~{}\mu\mathrm{m}\times[-20,20]~{}\mu\mathrm{m}$ with $h_{x}=h_{y}=0.1~{}\mu\mathrm{m}$. In our simulations, we choose the relaxation time $\tau_{0}=1.2$ picoseconds (ps), the surface conductivity $\sigma_{0}$ given as:

$$\sigma_{0}=-\frac{q^{2}k_{B}T\tau}{\pi\hbar^{2}}\left(\frac{\mu_{c}}{k_{B}T}+2\ln(\exp(-\frac{\mu_{c}}{k_{B}T})+1)\right),$$

where we denote the electron charge $q=1.6022\times e^{-19}\mathrm{C}$, the temperature $T=300~{}\mathrm{K}$, the reduced Plank constant $\hbar=1.0546\times e^{-34}\mathrm{J}\cdot\mathrm{s}$, the Boltzmann constant $k_{B}=1.3806\times e^{-23}\mathrm{J}\cdot\mathrm{K}^{-1}$, and the chemical potential $\mu_{c}=1.5~{}eV$ for Examples 1-4 and $\mu_{c}=0.8~{}eV$ for Example 5.

We use the time step $\tau=8.3\times 10^{-17}~{}s$, and run the simulation for 10000 time steps for Example 3, 20000 time steps for Examples 1, 2, 4, and 100000 time steps for Example 5.

4.2.1 Example 1. Bifurcated graphene sheets

In this example, we firstly reproduce the results in [26, Example 4] with our new scheme. The same setup with exemplary coarse mesh is shown in Figure 1.

Figure 1: A bifurcated graphene sheet buried in vacuum with an exemplary coarse mesh. The structure is constructed by segments $\mathrm{\mathrm{p}}_{1}\mathrm{p}_{2},\mathrm{p}_{2}\mathrm{p}_{3},\mathrm{p}_{3}\mathrm{p}_{4},\mathrm{p}_{2}\mathrm{p}_{5},\mathrm{p}_{5}\mathrm{p}_{6}$. The coordinates for the marked points on the graph are correspondingly $\mathrm{p}_{1}(-30,0)\mu\mathrm{m},\mathrm{p}_{2}(-15,0)\mu\mathrm{m},\mathrm{p}_{3}(0,5)\mu\mathrm{m},\mathrm{p}_{4}(15,5)\mu\mathrm{m},\mathrm{p}_{5}(0,-5)\mu\mathrm{m},\mathrm{p}_{6}(15,-5)\mu\mathrm{m}$.

Snapshots of the magnetic field Hz are presented in Figure 2, which clearly shows that a surface wave propagates along the graphene sheet same as [26, Figure 8]. To show the efficiency of our new algorithm, we made a comparison for the computational time between our new algorithm and the old algorithm proposed in our previous work [26]. The tests are carried out using NGSolve on a Mac mini with an Apple M2 chip. Under the same setup, our new algorithm takes 10,327.50 seconds, while the old algorithm needs 19,114.09 seconds. The new algorithm demonstrates a significant improvement in the computational time.

Figure 2: Snapshots of the contour plots of $H_{z}$ obtained at different time steps: 1000 (top left), 3000 (top right), 4000 (middle left), 5500 (middle right), 6000 (bottom left), and 10000 (bottom right).

Then we carry out a numerical simulation of the SPP propagation on a bifurcated curly graphene sheet. The setup with an exemplary coarse mesh is shown in Figure 3.

Figure 3: A bifurcated graphene sheet buried in vacuum with an exemplary coarse mesh. The structure is constructed by segments $\mathrm{p}_{1}\mathrm{p}_{2},\mathrm{p}_{3}\mathrm{p}_{5},\mathrm{p}_{4}\mathrm{p}_{6}$, and a semicircle centered at $(7,0)\mu\mathrm{m}$ with radius $r=7\mu\mathrm{m}$. The coordinates for the marked points on the graph are correspondingly $\mathrm{p}_{1}(-28,0)\mu\mathrm{m},\mathrm{p}_{2}(0,0)\mu\mathrm{m},\mathrm{p}_{3}(7,7)\mu\mathrm{m},\mathrm{p}_{4}(7,-7)\mu\mathrm{m},\mathrm{p}_{5}(15,7)\mu\mathrm{m},\mathrm{p}_{6}(15,-7)\mu\mathrm{m}$.

A pair of dipole source waves is placed at $(-27\mu\mathrm{m},1\mu\mathrm{m})\text{ and }(-27\mu\mathrm{m},-1\mu\mathrm{m})$ and imposed as $K_{s}=\pm\sin(2\pi f_{0}t)/h$ with $f_{0}=10~{}THz$ and $h=h_{y}$. Snapshots of the magnetic field $H_{z}$ are presented in Figure 4, which clearly show that a surface wave propagates along the graphene sheet.

Figure 4: Snapshots of the contour plots of $H_{z}$ obtained at different time steps: 1000 (top left), 3000 (top right), 5200 (middle left), 8000 (middle right), 10000 (bottom left), and 18800 (bottom right).

4.2.2 Example 2. Four adjacent curved graphene sheets

In this example, we present a numerical simulation of SPPs propagating along four separate but adjacent graphene sheets by our scheme. The simulation setup is shown in Figure 5, where four adjacent curved graphene sheets are embedded inside the physical domain $\Omega$. A pair of dipole source wave is placed at $(-18\mu\mathrm{m},3.5\mu\mathrm{m})\text{ and }(-18\mu\mathrm{m},2.5\mu\mathrm{m})$. The other setup is the same as Example 1.

Figure 5: The setup of four adjacent curved graphene sheets and a sample coarse mesh. The structure is constructed by 4 circular arcs. The coordinates for the marked points on the graph are correspondingly $\mathrm{p}_{1}(-24,0)\mu\mathrm{m},\mathrm{p}_{2}(-18,3)\mu\mathrm{m},\mathrm{p}_{3}(-12,0)\mu\mathrm{m},\mathrm{p}_{4}(-12,-2h_{y})\mu\mathrm{m},\mathrm{p}_{5}(-6,-3-2h_{y})\mu\mathrm{m},\\ \mathrm{p}_{6}(0,-2h_{y})\mu\mathrm{m},\mathrm{p}_{7}(0,0)\mu\mathrm{m},\mathrm{p}_{8}(6,3)\mu\mathrm{m},\mathrm{p}_{9}(12,0)\mu\mathrm{m},\mathrm{p}_{10}(12,-2h_{y})\mu\mathrm{m},\mathrm{p}_{11}(-18,-3-2h_{y})\mu\mathrm{m},\\ \mathrm{p}_{12}(-24,-2h_{y})\mu\mathrm{m}$.

Snapshots of the magnetic field $H_{z}$ are presented in Figure 6, which show that the surface wave propagates along the graphene sheets even though they are separated.

Figure 6: Snapshots of the contour plots of $H_{z}$ at different time steps: 1000 (top left), 3600 (top right), 4400 (middle left), 6000 (middle right), 11000 (bottom left), and 18000 (bottom right).

4.2.3 Example 3. A bulb shaped graphene sheet

In this example, we present a numerical simulation of the SPP propagating along a bulb shaped graphene sheet. The setup with an exemplary coarse mesh is shown in Figure 7. One pair of dipole source wave is placed at $(-15\mu\mathrm{m},2.5\mu\mathrm{m})\text{ and }(-15\mu\mathrm{m},1.5\mu\mathrm{m})$, and another pair is placed at $(-15\mu\mathrm{m},-2.5\mu\mathrm{m})\text{ and }(-15\mu\mathrm{m},-1.5\mu\mathrm{m})$ with the same function given as Example 1. We adopt the same parameters as Example 1.

The numerical results are shown in Figure 8. As we can see, the surface wave propagates along the graphene sheet.

Figure 7: A bulb shaped graphene sheet buried in $\Omega$ with an exemplary coarse mesh. The structure is constructed by two segments $\mathrm{p}_{1}\mathrm{p}_{2},\mathrm{p}_{3}\mathrm{p}_{4}$ and one circular arc. The coordinates for the marked points on the graph are correspondingly $\mathrm{p}_{1}(-15,2)\mu\mathrm{m},\mathrm{p}_{2}(0,2)\mu\mathrm{m},\mathrm{p}_{3}(-15,-2)\mu\mathrm{m},\mathrm{p}_{4}(0,-2)\mu\mathrm{m},\mathrm{p}_{5}(10,0)\mu\mathrm{m}$.

Figure 8: Snapshots of the contour plots of $H_{z}$ obtained at different time steps: 100 (top left), 1800 (top right), 3000 (middle left), 4000 (middle right), 5200 (bottom left), and 10000 (bottom right).

4.2.4 Example 4: Ring resonator graphene interface

In this example, we present a numerical simulation of an optical ring resonator motivated by the paper [10]. The resonator consists of a graphene ring with a radius of $11\mu\mathrm{m}$, centered at the origin (0,0), embeded within the domain $\Omega_{0}$. Two graphene segments, each with $30\mu\mathrm{m}$ in length are positioned at $y=13\mu\mathrm{m}$ and $y=-13\mu\mathrm{m}$, respectively. The setup with a coarse mesh is shown in Figure 8. A pair of dipole source wave is placed at $(-13\mu\mathrm{m},13.5\mu\mathrm{m})\text{ and }(-13\mu\mathrm{m},12.5\mu\mathrm{m})$ with the same parameters as Example 1.

The numerical results are displayed in Figure 10. The input wave propagates through the top graphene segment, reaching the transmission port and exiting to the left of the bottom graphene segment. This behavior mirrors the phenomenon observed in [10, Fig.6(a)], where the input wave enters through the bottom graphene segment and exits to the left of the top graphene segment. The opposite direction in our setup is due to the position of the source wave.

Figure 9: The graphene interface in $\Omega_{2}$ with a coarse mesh.

Figure 10: Snapshots of the contour plots of $H_{z}$ obtained at different time steps: 2000 (top left), 3000 (top middle), 5000 (top right), 7500 (bottom left), 11250 (bottom middle), and 19000 (bottom right).

4.2.5 Example 5: A spiral graphene interface

In this example, we present a numerical simulation of SPP surface wave propagation along a spiral graphene interface inspired by [43] with our scheme. A pair of dipole source wave is placed at $(-16\mu\mathrm{m},-18.5\mu\mathrm{m})\text{ and }(-16\mu\mathrm{m},-17.5\mu\mathrm{m})$.

The obtained numerical magnetic fields $H_{z}$ at various time steps are presented in Figure 12. As we can see, the wave propagates along the graphene interface.

Figure 11: The graphene interface in $\Omega_{2}$ with a coarse mesh. The spiral structure is constructed by 7 semicircles, one quarter circle and one segment. Let $w=4\mu\mathrm{m}$, and the coordinates for the marked points on the graph are correspondingly $p_{1}(-w,0),p_{2}(w,0),p_{3}(-2w,0),p_{4}(2w,0),p_{5}(-3w,0),p_{6}(3w,0),p_{7}(-4w,0),p_{8}(4w,0),p_{9}(0,-4.5w),\\ p_{10}(-4.5w,-4.5w)$.

Figure 12: Snapshots of the contour plots of $H_{z}$ obtained at various time steps: 2000 (top left), 8000 (top middle), 20000 (top right), 32000 (bottom left), 72000 (bottom middle), and 100000 (bottom right).

5 Conclusion

In this paper, we first developed a simplified graphene model by eliminating the surface current variable from the graphene model adopted in our previous work [26]. Then we established the stability for the reformulated PDE model, and proposed a new finite element method for solving it. Extensive numerical simulations were carried out to demonstrate that the reformulated model captures the surface plasmon polaritons very efficiently for various complex graphene sheets. Moreover, our new algorithm significantly improves computational efficiency, as shown by the time comparison results in Example 1, making it a more effective approach for simulating graphene-based plasmonic phenomena.

Acknowledgments. Dr. J. Li is very grateful to Dr. Frederic Marazzato at University of Arizona for some insightful discussions. The authors like to thank the two anonymous referees for their insightful suggestions on improving our paper.

6 Declarations

Funding Zhu’s work is supported by the NSF grant 2136228 “RTG: Program in Computation- and Data-Enabled Science”.

Conflicts of Interests/Competing Interest The authors have no conflicts of interest to declare that are relevant to the content of this article.