A Crank-Nicolson Finite Element Method and the Optimal Error Estimates for the modified Time-dependent Maxwell-Schr\"{o}dinger Equations

TL;DR

This paper introduces a Crank-Nicolson finite element method for the time-dependent Maxwell-Schrödinger equations, providing optimal error estimates and validating them through numerical tests, advancing numerical analysis for quantum-electromagnetic interactions.

Contribution

The paper develops a new finite element scheme with optimal error bounds for Maxwell-Schrödinger equations without time-step restrictions.

Findings

Optimal energy-norm error estimates derived

Numerical tests confirm theoretical error bounds

Method applicable to quantum-electromagnetic interaction problems

Abstract

In this paper we consider the initial-boundary value problem for the time-dependent Maxwell-Schr\"{o}dinger equations, which arises in the interaction between the matter and the electromagnetic field for the semiconductor quantum devices. A Crank-Nicolson finite element method for solving the problem is presented. The optimal energy-norm error estimates for the numerical algorithm without any time-step restrictions are derived. Numerical tests are then carried out to confirm the theoretical results.

Tables

Table 1. Table 1: H 1 superscript 𝐻 1 H^{1} error of linear FEM with h = 1 M ℎ 1 𝑀 h=\frac{1}{M} and Δ ​ t = h 1 2 Δ 𝑡 superscript ℎ 1 2 \Delta t=h^{\frac{1}{2}} .

		${\Vert {𝐀}_h^k-𝐀(\cdot ,t_k)\Vert }_{{𝐇}^1}$
t	M=25	M=50	M=100	Order
1.0	9.7908e-01	4.8501e-01	2.1082e-01	1.11
2.0	7.6414e-01	3.7807e-01	1.7681e-01	1.06
3.0	6.3094e-01	3.1006e-01	1.5308e-01	1.02
4.0	7.2739e-01	3.5204e-01	1.7705e-01	1.02
		${\Vert {\Psi }_h^k-\Psi (\cdot ,t_k)\Vert }_{{\mathscr{H}}^1}$
t	M=25	M=50	M=100	Order
1.0	6.8289e-01	3.3004e-01	1.5046e-01	1.09
2.0	8.0035e-01	3.6227e-01	1.6032e-01	1.16
3.0	4.1192e-01	1.9485e-01	1.0286e-01	1.00
4.0	2.3418e-01	1.1430e-01	5.6022e-02	1.03

Table 2. Table 2: H 1 superscript 𝐻 1 H^{1} error of quadratic FEM with h = Δ ​ t = 1 M . ℎ Δ 𝑡 1 𝑀 h=\Delta t=\frac{1}{M}.

		${\Vert {𝐀}_h^k-𝐀(\cdot ,t_k)\Vert }_{{𝐇}^1}$
t	M=25	M=50	M=100	Order
1.0	3.3770e-02	8.4984e-03	2.2743e-03	1.95
2.0	2.2786e-02	5.7068e-03	1.4966e-03	1.96
3.0	3.4016e-02	8.9115e-03	2.4360e-03	1.90
4.0	3.9787e-02	9.0740e-03	2.3467e-03	2.04
		${\Vert {\Psi }_h^k-\Psi (\cdot ,t_k)\Vert }_{{\mathscr{H}}^1}$
t	M=25	M=50	M=100	Order
1.0	2.8221e-02	6.8364e-03	1.8024e-03	1.98
2.0	4.5738e-02	1.1608e-02	2.6719e-03	2.05
3.0	3.2712e-02	8.5698e-03	2.2101e-03	1.94
4.0	2.1868e-02	5.1495e-03	1.3250e-03	2.02

Equations

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle\mathrm{i}\hbar\frac{\partial\Psi(\mathbf{x},t)}{\partial t}=\left\{\frac{1}{2m}\left[\mathrm{i}\hbar\nabla+q\mathbf{A}(\mathbf{x},t)\right]^{2}+q\Phi(\mathbf{x},t)+V_{0}\right\}\Psi(\mathbf{x},t),}\\[5.69054pt] {\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle-\frac{\partial}{\partial t}\nabla\cdot\big{(}\epsilon\mathbf{A}(\mathbf{x},t)\big{)}-\nabla\cdot\big{(}\epsilon\nabla\Phi(\mathbf{x},t)\big{)}=q|\Psi(\mathbf{x},t)|^{2},\,\,(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\epsilon\frac{\partial^{2}\mathbf{A}(\mathbf{x},t)}{\partial t^{2}}+\nabla\times\big{(}{\mu}^{-1}\nabla\times\mathbf{A}(\mathbf{x},t)\big{)}+\epsilon\frac{\partial(\nabla\Phi(\mathbf{x},t))}{\partial t}=\mathbf{J}_{q}(\mathbf{x},t),}\\[5.69054pt] {\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\mathbf{J}_{q}=-\frac{\mathrm{i}q\hbar}{2m}\big{(}\Psi^{\ast}\nabla{\Psi}-\Psi\nabla{\Psi}^{\ast}\big{)}-\frac{|q|^{2}}{m}|\Psi|^{2}\mathbf{A},}\\[5.69054pt] {\displaystyle\Psi,\Phi,\mathbf{A}\,\,\mathrm{subject\ to\ the\ appropriate\ initial\ and\ boundary\ conditions},}\end{array}\right.$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}\frac{\partial\Psi}{\partial t}+\frac{1}{2}\left(\mathrm{i}\nabla+\mathbf{A}\right)^{2}\Psi+V_{0}\Psi=0,\,\,(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\frac{\partial^{2}\mathbf{A}}{\partial t^{2}}+\nabla\times(\nabla\times\mathbf{A})+\frac{\mathrm{i}}{2}\big{(}\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}\big{)}+|\Psi|^{2}\mathbf{A}=0,\,\,(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] \end{array}\right.$$

$$-\frac{\partial }{\partial t}\nabla \cdot \mathbf{A}(\mathbf{x},t)=|\Psi (\mathbf{x},t){|}^2,(\mathbf{x},t)\in \Omega \times (0,T),$$

$$-\nabla \cdot \mathbf{A}(\mathbf{x},t)+\nabla \cdot \mathbf{A}(\mathbf{x},0)-{\int }_0^t\rho (\mathbf{x},\tau )d\tau =0,$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}\frac{\partial\Psi}{\partial t}+\frac{1}{2}\left(\mathrm{i}\nabla+\mathbf{A}\right)^{2}\Psi+V_{0}\Psi=0,\,\,(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\frac{\partial^{2}\mathbf{A}}{\partial t^{2}}+\nabla\times(\nabla\times\mathbf{A})-\gamma\nabla(\nabla\cdot\mathbf{A})+\frac{\mathrm{i}}{2}\big{(}\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}\big{)}+|\Psi|^{2}\mathbf{A}}\\[5.69054pt] {\displaystyle\quad\,+\,\gamma\nabla(\nabla\cdot\mathbf{A}(\mathbf{x},0))-\gamma\nabla\int_{0}^{t}\rho(\mathbf{x},\tau)d\tau=0,\,\,(\mathbf{x},t)\in\Omega\times(0,T).}\\[5.69054pt] \end{array}\right.$$

$$\rho (\mathbf{x},t)\approx \rho (\mathbf{x},0)+t\frac{\partial \rho }{\partial t}(\mathbf{x},0),\forall x\in \Omega$$

$${\int }_0^t\rho (\mathbf{x},\tau )d\tau \approx t\rho (\mathbf{x},0)+\frac{1}{2}{t}^2\frac{\partial \rho }{\partial t}(\mathbf{x},0),\forall x\in \Omega .$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\int_{0}^{t}\rho(\mathbf{x},\tau)d\tau=\int_{0}^{t_{1}}\rho(\mathbf{x},\tau)d\tau+\int_{t_{1}}^{t}\rho(\mathbf{x},\tau)d\tau}\\[5.69054pt] {\displaystyle\quad\quad\quad\approx\int_{0}^{t_{1}}\rho(\mathbf{x},\tau)d\tau+(t-t_{1})\rho(\mathbf{x},t_{1})+\frac{1}{2}(t-t_{1})^{2}\frac{\partial\rho}{\partial t}(\mathbf{x},t_{1}),\quad\forall x\in\Omega}.\end{array}$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}\frac{\partial\Psi}{\partial t}+\frac{1}{2}\left(\mathrm{i}\nabla+\mathbf{A}\right)^{2}\Psi+V_{0}\Psi=0,\,\,(\mathbf{x},t)\in\Omega\times(t_{j-1},t_{j}],}\\[5.69054pt] {\displaystyle\frac{\partial^{2}\mathbf{A}}{\partial t^{2}}+\nabla\times(\nabla\times\mathbf{A})-\gamma\nabla(\nabla\cdot\mathbf{A})+\frac{\mathrm{i}}{2}\big{(}\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}\big{)}}\\[5.69054pt] {\displaystyle\quad\quad+\,\,|\Psi|^{2}\mathbf{A}=\mathbf{f}_{j}(\mathbf{x},t),\quad\,\,(\mathbf{x},t)\in\Omega\times(t_{j-1},t_{j}],}\\[5.69054pt] \end{array}\right.$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}\frac{\partial\Psi}{\partial t}+\frac{1}{2}\left(\mathrm{i}\nabla+\mathbf{A}\right)^{2}\Psi+V_{0}\Psi=0,\,\,(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\frac{\partial^{2}\mathbf{A}}{\partial t^{2}}+\nabla\times(\nabla\times\mathbf{A})-\gamma\nabla(\nabla\cdot\mathbf{A})+\frac{\mathrm{i}}{2}\big{(}\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}\big{)}}\\[5.69054pt] {\displaystyle\quad\quad+\,\,|\Psi|^{2}\mathbf{A}=\mathbf{g}(\mathbf{x},t),\quad\,\,(\mathbf{x},t)\in\Omega\times(0,T).}\\[5.69054pt] \end{array}\right.$$

$$\Psi (\mathbf{x},t)=0,\mathbf{A}(\mathbf{x},t)\times \mathbf{n}=0,\nabla \cdot \mathbf{A}(\mathbf{x},t)=0,(\mathbf{x},t)\in \partial \Omega \times (0,T),$$

$$\Psi (\mathbf{x},0)={\Psi }_0(\mathbf{x}),\mathbf{A}(\mathbf{x},0)={\mathbf{A}}_0(\mathbf{x}),{\mathbf{A}}_t(\mathbf{x},0)={\mathbf{A}}_1(\mathbf{x}),$$

$${\mathbf{H}}_t^1(\Omega )=\{\mathbf{A}|\mathbf{A}\in {\mathbf{H}}^1(\Omega ),\mathbf{A}\times \mathbf{n}=0\mathrm{on}\partial \Omega \}$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{H}}_t^1(\Omega )}:={\left[\Vert \nabla \cdot \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )}^2+\Vert \nabla \times \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )}^2\right]}^{\frac{1}{2}},$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}(\frac{\partial\Psi}{\partial t},\varphi)+\frac{1}{2}(\left(\mathrm{i}\nabla+\mathbf{A}\right)\Psi,\left(\mathrm{i}\nabla+\mathbf{A}\right)\varphi)+(V_{0}\Psi,\varphi)=0,\quad\forall\varphi\in\mathcal{H}_{0}^{1}(\Omega)}\\[5.69054pt] {\displaystyle(\frac{\partial^{2}\mathbf{A}}{\partial t^{2}},\mathbf{v})+(\nabla\times\mathbf{A},\nabla\times\mathbf{v})+\gamma(\nabla\cdot\mathbf{A},\nabla\cdot\mathbf{v})+(\frac{\mathrm{i}}{2}\big{(}\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}\big{)},\mathbf{v})}\\[5.69054pt] {\displaystyle\quad\quad\quad+\,(|\Psi|^{2}\mathbf{A},\mathbf{v})=(\mathbf{g},\mathbf{v}),\qquad\qquad\qquad\qquad\qquad\forall\mathbf{v}\in\mathbf{H}^{1}_{t}(\Omega),}\end{array}\right.$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\partial U^{k}=(U^{k}-U^{k-1})/\Delta t,\quad\partial^{2}U^{k}=(\partial U^{k}-\partial U^{k-1})/\Delta t,}\\[5.69054pt] {\displaystyle\overline{U}^{k}=(U^{k}+U^{k-1})/2,\quad\widetilde{U}^{k}=(U^{k}+U^{k-2})/2,}\\[5.69054pt] \end{array}$$

$$Y_h^r=\{u_h\in C(\overline{\Omega }):u_h{|}_K\in P_r(K),\forall K\in {\mathcal{T}}_h\}.$$

$$V^{r}_{h}=Y^{r}_{h}\cap H_{0}^{1}(\Omega),\quad\mathcal{V}^{r}_{h}=V^{r}_{h}\oplus{\rm i}V^{r}_{h},\quad\mathbf{V}^{r}_{h}=\big{(}Y^{r}_{h}\big{)}^{3}\cap\mathbf{H}^{1}_{\rm t}(\Omega).$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\|u-I_{h}u\|_{\mathcal{W}^{s}_{p}}\leq Ch^{m-s}\|u\|_{\mathcal{W}^{m}_{p}}\quad\forall\;u\in\mathcal{W}^{m}_{p}(\Omega),}\\[5.69054pt] {\displaystyle\|\mathbf{v}-{\bm{\pi}}_{h}\mathbf{v}\|_{\mathbf{W}^{s}_{p}}\leq Ch^{m-s}\|\mathbf{v}\|_{\mathbf{W}^{m}_{p}}\quad\forall\;\mathbf{v}\in\mathbf{W}^{m}_{p}(\Omega).}\end{array}$$

$$\mathbf{A}(\cdot ,-\Delta t)=\mathbf{A}(\cdot ,0)-\Delta t\frac{\partial \mathbf{A}}{\partial t}(\cdot ,0)={\mathbf{A}}_0-\Delta t{\mathbf{A}}_1,$$

$${\Psi }_h^0=I_h{\Psi }_0,{\mathbf{A}}_h^0={\mathbf{\pi }}_h{\mathbf{A}}_0,{\mathbf{A}}_h^0-{\mathbf{A}}_h^{-1}=\Delta t{\mathbf{\pi }}_h{\mathbf{A}}_1,$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}(\partial\Psi_{h}^{k},\varphi)+\frac{1}{2}\left((\mathrm{i}\nabla+\overline{\mathbf{A}}_{h}^{k})\overline{\Psi}_{h}^{k},(\mathrm{i}\nabla+\overline{\mathbf{A}}_{h}^{k})\varphi\right)+(V_{0}\overline{\Psi}_{h}^{k},\varphi)=0,\quad\forall\varphi\in\mathcal{V}_{h}^{r}}\\[5.69054pt] {\displaystyle(\partial^{2}\mathbf{A}_{h}^{k},\mathbf{v})+\left(\frac{\mathrm{i}}{2}\big{(}(\Psi_{h}^{k-1})^{\ast}\nabla{\Psi_{h}^{k-1}}-\Psi_{h}^{k-1}\nabla{(\Psi_{h}^{k-1})}^{\ast}\big{)},\mathbf{v}\right)+(\nabla\times\widetilde{\mathbf{A}}_{h}^{k},\nabla\times\mathbf{v})}\\[5.69054pt] {\displaystyle\quad+\gamma(\nabla\cdot\widetilde{\mathbf{A}}_{h}^{k},\nabla\cdot\mathbf{v})+\big{(}|\Psi_{h}^{k-1}|^{2}\frac{\overline{\mathbf{A}}_{h}^{k}+\overline{\mathbf{A}}_{h}^{k-1}}{2},\mathbf{v}\big{)}=(\mathbf{g}^{k-1},\mathbf{v}),\quad\forall\mathbf{v}\in\mathbf{V}^{r}_{h},}\end{array}\right.$$

$$\begin{array}[]{@{}l@{}}{\displaystyle B(\mathbf{A};\Psi,\varphi)=\left((\mathrm{i}\nabla+\mathbf{A})\Psi,(\mathrm{i}\nabla+\mathbf{A})\varphi\right),}\\[5.69054pt] {\displaystyle D(\mathbf{A},\mathbf{v})=\gamma(\nabla\cdot\mathbf{A},\nabla\cdot\mathbf{v})+(\nabla\times\mathbf{A},\nabla\times\mathbf{v}),}\\[5.69054pt] {\displaystyle f(\Psi,\varphi)=\frac{\mathrm{i}}{2}(\varphi^{\ast}\nabla\Psi-\Psi\nabla\varphi^{\ast}).}\end{array}$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}(\frac{\partial\Psi}{\partial t},\varphi)+\frac{1}{2}B(\mathbf{A};\Psi,\varphi)+(V_{0}\Psi,\varphi)=0,\quad\forall\varphi\in\mathcal{H}_{0}^{1}(\Omega),}\\[5.69054pt] {\displaystyle(\frac{\partial^{2}\mathbf{A}}{\partial t^{2}},\mathbf{v})+D(\mathbf{A},\mathbf{v})+(f(\Psi,\Psi),\mathbf{v})+(|\Psi|^{2}\mathbf{A},\mathbf{v})=(\mathbf{g},\mathbf{v}),\quad\forall\mathbf{v}\in\mathbf{H}^{1}_{t}(\Omega),}\end{array}\right.$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}(\partial\Psi_{h}^{k},\varphi)+\frac{1}{2}B(\overline{\mathbf{A}}_{h}^{k};\overline{\Psi}_{h}^{k},\varphi)+(V_{0}\overline{\Psi}_{h}^{k},\varphi)=0,\quad\forall\varphi\in\mathcal{V}_{h}^{r}}\\[5.69054pt] {\displaystyle(\partial^{2}\mathbf{A}_{h}^{k},\mathbf{v})+D(\widetilde{\mathbf{A}}_{h}^{k},\mathbf{v})+\left(f(\Psi_{h}^{k-1},\Psi_{h}^{k-1}),\mathbf{v}\right)}\\[5.69054pt] {\displaystyle\quad+\big{(}|\Psi_{h}^{k-1}|^{2}\frac{(\overline{\mathbf{A}}_{h}^{k}+\overline{\mathbf{A}}_{h}^{k-1})}{2},\mathbf{v}\big{)}=(\mathbf{g}^{k-1},\mathbf{v}),\quad\forall\mathbf{v}\in\mathbf{V}^{r}_{h}.}\end{array}\right.$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\Psi,\Psi_{t},\Psi_{tt}\in{L}^{\infty}(0,T;\mathcal{H}^{r+1}(\Omega)),\quad\Psi_{ttt}\in{L}^{\infty}(0,T;\mathcal{H}^{1}(\Omega)),}\\[5.69054pt] {\displaystyle\Psi_{tttt}\in L^{2}(0,T;\mathcal{L}^{2}(\Omega)),\quad\mathbf{A},\mathbf{A}_{t},\mathbf{A}_{tt}\in{L}^{\infty}(0,T;\mathbf{H}^{r+1}(\Omega)),}\\[5.69054pt] {\displaystyle\mathbf{A}_{ttt}\in{L}^{\infty}(0,T;\mathbf{H}^{1}(\Omega)),\mathbf{A}_{tttt}\in L^{2}(0,T;\mathbf{L}^{2}(\Omega)).}\end{array}$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\Psi_{0}\in\mathcal{H}^{r+1}(\Omega)\cap\mathcal{H}_{0}^{1}(\Omega),\,\,\,\mathbf{A}_{0},\mathbf{A}_{1}\in\mathbf{H}^{r+1}(\Omega)\cap\mathbf{H}_{t}^{1}(\Omega),\,\,\,\mathbf{g}\in C(0,T;\mathbf{L}^{2}(\Omega)).}\end{array}$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\max_{1\leq k\leq M}\big{\{}\|\Psi_{h}^{k}-\Psi^{k}\|_{\mathcal{L}^{2}(\Omega)}^{2}+\|\nabla(\Psi_{h}^{k}-\Psi^{k})\|_{\mathbf{L}^{2}(\Omega)}^{2}+\|\mathbf{A}_{h}^{k}-\mathbf{A}^{k}\|_{\mathbf{L}^{2}(\Omega)}^{2}}\\[5.69054pt] {\displaystyle\quad+\|\nabla\cdot(\mathbf{A}_{h}^{k}-\mathbf{A}^{k})\|_{{L}^{2}(\Omega)}^{2}+\|\nabla\times(\mathbf{A}_{h}^{k}-\mathbf{A}^{k})\|_{\mathbf{L}^{2}(\Omega)}^{2}\big{\}}\leq C\big{\{}h^{2r}+(\Delta t)^{4}\big{\}},}\end{array}$$

$$\Vert u{\Vert }_{L^p}\le C\Vert u{\Vert }_{H^1},\Vert \mathbf{v}{\Vert }_{{\mathbf{L}}^p}\le C\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}^1},1\le p\le 6(d=2,3),$$

$$\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}^1}\le C(\Vert \nabla \times \mathbf{v}{\Vert }_{{\mathbf{L}}^2}+\Vert \nabla \cdot \mathbf{v}{\Vert }_{L^2}),\mathbf{v}\in {\mathbf{H}}_t^1(\Omega ),$$

Full text

1 Introduction

When the characteristic size of the semiconductor device reaches the wavelength of an electron, the quantum effects become important even dominant and can not be neglected. The accurate electromagnetic theory for the case is quantum electrodynamics (QED), i.e. the second quantization for the matter and quantization for the electromagnetic field. However, so far it is extremely difficult even impossible to employ QED to analyze the interaction between the matter and the electromagnetic field for some complex systems. The semiclassical (or semi-quantum) electromagnetic models are widely used in the semiconductor quantum devices. The basic idea is that we use the Maxwell’s equations for the electromagnetic field while we use the Schrödinger equation of the non-relativistic quantum mechanics for the matter (see [Ez, Sch]). The Maxwell-Schrödinger coupled system (M-S) is written as follows:

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$ is a bounded Lipschitz polyhedral convex domain, $\Psi^{\ast}$ denotes the complex conjugate of $\Psi$, $\epsilon$ and $\mu$ respectively denote the electric permittivity and the magnetic permeability of the material and $V_{0}$ is the constant potential energy.

It is well-known that the solutions of the Maxwell-Schrödinger equations (1) are not unique. In fact, for any function $\chi:\Omega\times(0,T)\rightarrow\mathbb{R}$, if $(\Psi,\Phi,\mathbf{A})$ is a solution of (1), then $(\exp(i\chi)\Psi,\Phi-\partial_{t}\chi,\mathbf{A}+\nabla\chi)$ is also a solution of (1). It is often assumed that the further equations can be adjoined to the Maxwell-Schrödinger equations by means of a gauge transformation. In this paper we consider the M-S system (1) under the temporal gauge (also called Weyl gauge), i.e. $\Phi=0$.

In this paper we employ the atomic units, i.e. $\hbar=m=q=1$. For simplicity, we also assume that $\epsilon=\mu=1$ without loss of generality. Hence, $\Psi$ and $\mathbf{A}$ satisfy the following Maxwell-Schrödinger equations :

[TABLE]

Here we omit the initial and boundary conditions for $\Psi$ and $\mathbf{A}$ temporarily.

Under the temporal gauge, the second equation in (1) involving the divergence of $\mathbf{A}$ can be rewritten as

[TABLE]

which can be derived from (2) if the solutions of (2) are sufficiently smooth and the initial datas are consistent.

Integrating with respect to $t$ on the both sides of (3), we have

[TABLE]

where $\rho(\mathbf{x},t)=|\Psi(\mathbf{x},t)|^{2}$.

For the purpose of theoretical analysis, we take the gradient of (4), multiply it by a parameter $\gamma>0$ and add it to the second equation of (2), to obtain

[TABLE]

The parameter $\gamma$ is referred to as the penalty factor. The choice of $\gamma$ depends on how much emphasis one places on the equality (3). In this paper, we keep $\gamma$ fixed. To avoid the difficulty for integro-differential equations, assuming that the change of the density function $\rho(\mathbf{x},t)$ is smooth with respect to $t$ for all $\mathbf{x}\in\Omega$, we give an approximation of $\int_{0}^{t}\rho(\mathbf{x},\tau)d\tau$ as follows.

First denoting by $t_{j}=jT/M,j=0,1,\cdots,M$, we divide the time interval $[0,T]$ into $M$ subintervals $[0,t_{1}],(t_{1},t_{2}],\cdots,(t_{M-1},T]$. For $t\in[0,t_{1}]$, $\rho(\mathbf{x},t)$ is approximated by Taylor expansion and the initial conditions:

[TABLE]

and

[TABLE]

The computation of $\frac{\displaystyle\partial\rho}{\displaystyle\partial t}(\mathbf{x},0)$ involves $\frac{\displaystyle\partial\Psi}{\displaystyle\partial t}(\mathbf{x},0)$, the time derivative of initial wave function. Here we assume the initial conditions are consistent and so we can obtain $\frac{\displaystyle\partial\Psi}{\displaystyle\partial t}(\mathbf{x},0)$ from Schrödinger’s equation. Given an approximation of $\int_{0}^{t}\rho(\mathbf{x},\tau)d\tau$ in $[0,t_{1}]$, we can solve the coupled differential equations (5) in the subinterval $[0,t_{1}]$ and integrate the density function to obtain $\int_{0}^{t_{1}}\rho(\mathbf{x},\tau)d\tau$. Then for $t\in(t_{1},t_{2}]$, $\int_{0}^{t}\rho(\mathbf{x},\tau)d\tau$ can be calculated as follows:

[TABLE]

Now we solve the Maxwell-Schrödinger equations (5) in the subinterval $(t_{1},t_{2}]$. Repeating the above procedure, we can solve the Maxwell-Schrödinger equations (5) in the subinterval $(t_{2},t_{3}],\cdots,(t_{M-1},T]$ successively. Therefore, we decompose the original system in $[0,T]$ into M system in $[0,t_{1}],(t_{1},t_{2}],\cdots,(t_{M-1},T]$, respectively. For $t\in(t_{j-1},t_{j}],j=1,2,\cdots,M$, the Maxwell-Schrödinger equations (5) can be rewritten as follows:

[TABLE]

where $\mathbf{f}_{j}(\mathbf{x},t)$ is the known function.

** Remark 1.1****.**

We can get the modified Maxwell-Schrödinger equations (6) under the assumption that the change of the density function $\rho(\mathbf{x},t)$ is smooth with respect to $t$ for all $\mathbf{x}\in\Omega$. If the initial wave function $\Psi(\mathbf{x},0)$ is the eigenfunction of the stationary Schrödinger equation and the incoming electromagnetic field is weak and can be considered as a small perturbation to the quantum system, this assumption is reasonable. The choice of the number $M$ of subintervals depends on the initial wave function, the incoming electromagnetic field and $T$.

In this paper, we consider the following modified Maxwell-Schrödinger equations:

[TABLE]

The boundary conditions are

[TABLE]

and the initial conditions are

[TABLE]

where $\mathbf{A}_{t}$ denotes the derivative of $\mathbf{A}$ with respect to the time $t$, $\mathbf{n}=(n_{1},n_{2},n_{3})$ is the outward unit normal to the boundary $\partial\Omega$. We assume that $\nabla\cdot\mathbf{A}_{0}=\nabla\cdot\mathbf{A}_{1}=0$ on $\partial\Omega$.

** Remark 1.2****.**

The boundary condition $\Psi(x)=0$ on $\partial\Omega$ implies that the particle is confined in a whole domain $\Omega$. The boundary condition $\mathbf{A}(\mathbf{x},t)\times\mathbf{n}=0$ on $\partial\Omega$ is referred to as the perfect conductive boundary condition. The boundary condition $\nabla\cdot\mathbf{A}(\mathbf{x},t)=0$ on $\partial\Omega$ can be deduced from the boundary condition of $\,\Psi$ and (4) if the initial conditions $\mathbf{A}_{0}$ and $\mathbf{A}_{1}$ satisfy $\nabla\cdot\mathbf{A}_{0}=\nabla\cdot\mathbf{A}_{1}=0$ on $\partial\Omega$. As for the determination of the boundary conditions for the vector potential $\mathbf{A}$, we refer to [5].

** Remark 1.3****.**

The existence and uniqueness of the solution for the time-dependent Maxwell-Schrödinger equations (2) have been investigated in [11, 13, 19, 20, 21, 26, 31]. However, the results of the well-posedness of the problem were obtained only for the Cauchy problem in $\mathbf{R}^{d}$, $d\geq 1$ instead of the initial-boundary value problem. To the best of our knowledge, the existence and uniqueness of the solution for the Maxwell-Schrödinger equations in a bounded domain seem to be open. For the modified equations (7)-(9), we will investigate the existence of solutions in another paper.

Many authors have discussed the numerical methods for the time-dependent Maxwell-Schrödinger equations. We recall some important studies about the problem. Sui and his collaborators [28] used the finite-difference time-domain (FDTD) method to solve the Maxwell-Schrödinger equations and to simulate a simple electron tunneling problem. Pierantoni, Mencarelli and Rozzi [23] applied the transmission line matrix method(TLM) to solve the Maxwell’s equations and employed the FDTD method to solve the Schrödinger equation, and did the simulation for a carbon nanotube between two metallic electrodes. Ahmed and Li [1] used the FDTD method for the Maxwell-Schrödinger system to simulate plasmonics nanodevices. The numerical studies listed above all include a step where they extract the vector potential $\mathbf{A}$ and the scalar potential $\Phi$ from the electric field $\mathbf{E}$ and the magnetic field $\mathbf{H}$ after solving the Maxwell’s equations involving $\mathbf{E}$ and $\mathbf{H}$. Recently, Ryu [24] employed directly the FDTD scheme to discretize the Maxwell-Schrödinger equations (1) under the Lorentz gauge and to simulate a single electron in an artificial atom excited by an incoming electromagnetic field. Other related studies on this topic have been reported in [22, 25, 30] and the references therein.

There are few results on the finite element method (FEM) of the Maxwell-Schrödinger equations and the convergence analysis. In this paper we will present a Crank-Nicolson finite element method for solving the problem (7)-(9), i.e. the finite element method in space and the Crank-Nicolson scheme in time. Then we will derive the optimal error estimates for the proposed method. Roughly speaking, compared with explicit algorithms such as the FDTD method, our method is more stable and suffers from less restriction in the time step-size since we use the Crank-Nicolson scheme in the time direction. Moreover, our method is more appropriate to deal with materials with discontinuous electromagnetic coefficients than the FDTD method. our work is motivated by [18] in which Mu and Huang proposed an alternating Crank-Nicolson method for the time-dependent Ginzburg-Landau equations. The optimal error estimates were derived under the time step restrictive conditions $\Delta t\leq O(h^{\frac{11}{12}})$ for the two-dimension model and $\Delta t\leq O(h^{2})$ for the three-dimension model, where $h$ and $\Delta t$ are the spatial mesh size and the time step, respectively. The related convergence results associated with the time-dependent Ginzburg-Landau equations can be also given in [3, 4, 6, 8, 17]. It should be emphasized that although the time-dependent Ginzburg-Landau model is somehow formally similar to the time-dependent Maxwell-Schrödinger system, there exists the essential difference between them. The former is classified as a parabolic system and the latter is a hyperbolic system. The main key point in our work is how to avoid using the finite element inverse estimates when dealing with the nonlinear terms. The new ideas are to derive the energy-norm error estimates for the Schrödinger’s equation, and to employ some tricks to eliminate the nonlinear terms both in the Schrödinger’s equation and in Maxwell’s equations, respectively.

The remainder of this paper is organized as follows. In section 2, a Crank-Nicolson scheme with the Galerkin finite element approximation for the modified Maxwell-Schrödinger equations (7)-(9) is developed. In section 3, the stability estimates are given. The optimal error estimates for the numerical solution without any restriction on time step are derived in section 4. Finally, the numerical testes are then carried out to confirm the theoretical results.

Throughout this paper, we denote by $C$ a generic positive constant independent of the mesh size and the time step without distinction.

2 A Crank-Nicolson Galerkin finite element scheme

In this section, we present a numerical scheme for the modified Maxwell-Schrödinger equations (7)-(9) using Galerkin finite element method in space and the Crank-Nicolson scheme in time. To start with, here and afterwards, we assume that $\Omega$ is a bounded Lipschitz polygonal convex domain in $\mathbb{R}^{2}$ (or a bounded Lipschitz polyhedron convex domain in $\mathbb{R}^{3}$).

We introduce the following notation. Let $W^{s,p}(\Omega)$ denote the conventional Sobolev spaces of the real-valued functions. As usual, $W^{s,2}(\Omega)$ and $W^{s,2}_{0}(\Omega)$ are denoted by $H^{s}(\Omega)$ and $H^{s}_{0}(\Omega)$ respectively. We use $\mathcal{W}^{s,p}(\Omega)=\{u+\mathrm{i}v\,|\,u,v\in W^{s,p}(\Omega)\}$ and $\mathcal{H}^{s}(\Omega)=\{u+\mathrm{i}v\,|\,u,v\in H^{s}(\Omega)\}$ with calligraphic letters for Sobolev spaces of the complex-valued functions, respectively. Furthermore, let $\mathbf{W}^{s,p}(\Omega)=[W^{s,p}(\Omega)]^{d}$ and $\mathbf{H}^{s}(\Omega)=[H^{s}(\Omega)]^{d}$ with bold faced letters be Sobolev spaces of the vector-valued functions with $d$ components ($d$=2, 3). $L^{2}$ inner-products in $H^{s}(\Omega)$, $\mathcal{H}^{s}(\Omega)$ and $\mathbf{H}^{s}(\Omega)$ are denoted by $(\cdot,\cdot)$ without ambiguity.

In particular, we introduce the following subspace of $\mathbf{H}^{1}(\Omega)$:

[TABLE]

The semi-norm on $\mathbf{H}^{1}_{t}(\Omega)$ is defined by

[TABLE]

which is equivalent to the standard $\mathbf{H}^{1}(\Omega)$-norm $\|\mathbf{u}\|_{\mathbf{H}^{1}(\Omega)}$, see [12].

To take into account the time-independence, for a time $T>0$ fixed, let $L^{p}(0,T;X)$ be the Bochner space defined in [27] for $p\in[1,\infty]$ and a Banach space $X$.

The weak formulation of the Maxwell-Schrödinger system (7)- (9) can be specified as follows: given $\mathbf{g}\in{L}^{2}(0,T;\mathbf{L}^{2}(\Omega))$, find $(\Psi,\mathbf{A})\in\mathcal{H}_{0}^{1}(\Omega)\times\mathbf{H}_{t}^{1}(\Omega)$ such that $\forall t\in(0,T)$,

[TABLE]

with the initial conditions $\Psi_{0}\in\mathcal{H}^{1}_{0}(\Omega)$, $\mathbf{A}_{0}\in\mathbf{H}^{1}_{t}(\Omega)$ and $\mathbf{A}_{t}(\cdot,0)\in\mathbf{L}^{2}(\Omega)$.

Let M be a positive integer and let $\Delta t=T/M$ be the time step. For any k=1,2,$\cdots,M$, we introduce the following notation:

[TABLE]

for any given sequence $\{U^{k}\}_{0}^{M}$ and denote $u^{k}=u(\cdot,t^{k})$ for any given functions $u\in C(0,T;\,X)$ with a Banach space $X$.

Let $\mathcal{T}_{h}=\{K\}$ be a regular partition of $\Omega$ into triangles in $\mathbb{R}^{2}$ or tetrahedrons in $\mathbb{R}^{3}$ without loss of generality, where the mesh size $h=\mathrm{max}_{K\in\mathcal{T}_{h}}\{diam(K)\}$. For any $K\in\mathcal{T}_{h}$, we denote by $P_{r}(K)$ the spaces of polynomials of degree $r\;(r\geq 1)$ defined on $K$. We now define the standard Lagrange finite element space

[TABLE]

We have the following finite element subspaces of $H_{0}^{1}(\Omega)$, $\mathcal{H}_{0}^{1}(\Omega)$ and $\mathbf{H}^{1}_{\rm t}(\Omega)$

[TABLE]

We shall approximate the wave function $\Psi$ and the vector potential $\mathbf{A}$ by the functions in $\mathcal{V}_{h}^{r}$ and $\mathbf{V}_{h}^{r}$, respectively. Let $I_{h}$ and ${\bm{\pi}}_{h}$ be the conventional pointwise interpolation operators on $\mathcal{V}_{h}^{r}$ and $\mathbf{V}_{h}^{r}$, respectively. For $0\leq s\leq m\leq r+1$, $m\geq 2$, $2\leq p\leq\infty$, standard finite element theory gives that [2]:

[TABLE]

For convenience, assume that the function $\mathbf{A}$ is defined in the interval $[-\Delta t,T]$ in terms of the time variable $t$. We can compute $\mathbf{A}(\cdot,-\Delta t)$ by

[TABLE]

which leads to an approximation to $\mathbf{A}^{-1}$ with second order accuracy.

A Crank-Nicolson Galerkin finite element approximation to the Maxwell-Schrödinger system (10) is formulated as follows:

[TABLE]

and find $(\Psi_{h}^{k},\mathbf{A}_{h}^{k})\in\mathcal{V}_{h}^{r}\times\mathbf{V}^{r}_{h}$ such that for $k=1,2,\cdots,M$,

[TABLE]

where $\overline{\mathbf{A}}_{h}^{k}$, $\overline{\Psi}_{h}^{k}$ and $\widetilde{\mathbf{A}}_{h}^{k}$ have been defined in (11).

For convenience, we define the following bilinear forms:

[TABLE]

Then the variational forms of the modified Maxwell-Schrödinger equations and the discrete system can be written as follows:

[TABLE]

and for $k=1,2,\cdots,M$,

[TABLE]

Note that after discretization in time and space, the Maxwell equation and Schrödinger equation in the discrete system (20) are decoupled. At each time step, we only need to solve the two discrete linear equations alternately.

In this paper we assume that the modified Maxwell-Schrödinger equations (19) has one and only one weak solution $(\Psi,\mathbf{A})$ and the following regularity conditions are satisfied:

[TABLE]

For the initial conditions $(\Psi_{0},\mathbf{A}_{0},\mathbf{A}_{1})$ and the right hand function $\mathbf{g}(\mathbf{x},t)$, we assume that

[TABLE]

We now give the main convergence result in this paper as follows:

Theorem 1**.**

Suppose that $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$ is a bounded Lipschitz polyhedral convex domain. Let $(\Psi,\mathbf{A})$ be the unique solution of the modified Maxwell-Schrödinger equations (19), and let $(\Psi_{h}^{k},\mathbf{A}_{h}^{k})$ be the numerical solution of the full discrete scheme (20) associated with (19). Under the assumptions (21) and (22), we have the following error estimates

[TABLE]

where $\Psi^{k}=\Psi(\cdot,t^{k})$, $\mathbf{A}^{k}=\mathbf{A}(\cdot,t^{k})$, and $C$ is a constant independent of $h$, $\Delta t$.

The proof of Theorem 1 will be given in section 4.

3 Stability estimates

In this section we derive some stability estimates for the numerical solutions of the full discrete system (20), which play an important role in the error estimates in the next section.

For convenience, we list the following imbedding inequalities and interpolation inequalities in Sobolev spaces (see, e.g., [14] and [12]), and use them in the sequel:

[TABLE]

[TABLE]

[TABLE]

where $\|u\|_{L^{p}}=\|u\|_{L^{p}(\Omega)}$, $\|u\|_{{H}^{1}}=\|u\|_{{H}^{1}(\Omega)}$, $\|\mathbf{v}\|_{\mathbf{L}^{2}}=\|\mathbf{v}\|_{\mathbf{L}^{2}(\Omega)}$ and $\|\mathbf{v}\|_{\mathbf{H}^{1}}=\|\mathbf{v}\|_{\mathbf{H}^{1}(\Omega)}$.

We first give the definition of the discrete energy functional of (20) as follows:

[TABLE]

where \overline{D(\mathbf{A}_{h}^{k},\mathbf{A}_{h}^{k})}=\frac{\displaystyle 1}{\displaystyle 2}\Big{[}D(\mathbf{A}_{h}^{k},\mathbf{A}_{h}^{k})+D(\mathbf{A}_{h}^{k-1},\mathbf{A}_{h}^{k-1})\Big{]}.

Lemma 2**.**

For the solution of (20), for $k=1,2\cdots,M$, we have

[TABLE]

where C is a constant independent of k, $h$ and $\Delta t$.

Proof.

Choosing $\varphi=\overline{\Psi}_{h}^{k}$ in $(\ref{eq:2-9})_{1}$ and taking the imaginary part, we can complete the proof of $(\ref{eq:3-9})_{1}$. Let us turn to the proof of $(\ref{eq:3-9})_{2}$. It is obvious that

[TABLE]

By a direct computation, we get

[TABLE]

and consequently

[TABLE]

We thus have

[TABLE]

It is not difficult to check that

[TABLE]

We choose $\varphi=\partial\Psi_{h}^{k}$ in $(\ref{eq:2-9})_{1}$ and take the real part. Combining (30) and (31) gives

[TABLE]

Taking $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2\Delta t}(\mathbf{A}_{h}^{k}-\mathbf{A}_{h}^{k-2})=\frac{\displaystyle 1}{\displaystyle\Delta t}(\overline{\mathbf{A}}_{h}^{k}-\overline{\mathbf{A}}_{h}^{k-1})=\frac{\displaystyle 1}{\displaystyle 2}(\partial\mathbf{A}_{h}^{k}+\partial\mathbf{A}_{h}^{k-1})$ in $(\ref{eq:2-9})_{2}$, and combining with (32), we get

[TABLE]

It follows that

[TABLE]

Multiply (33) by $\Delta t$, sum $k=1,2,\cdots,M$, to discover

[TABLE]

Now $(\ref{eq:3-9})_{2}$ follows from the discrete Gronwall’s inequality and thus we complete the proof of Lemma 2.   ∎

** Remark 3.1****.**

Lemma 2 shows that the numerical scheme presented in this paper for the modified Maxwell-Schrödinger equations (19) is stable in some senses.

Theorem 3**.**

The solution of the full discrete system (20) fulfills the following estimates

[TABLE]

and

[TABLE]

where $C$ is a constant independent of $h$, $\Delta t$.

Proof.

(34) is the direct result of Lemma  2. Next we give the proof of (35). Since the semi-norm in $\mathbf{H}_{t}^{1}(\Omega)$ is equivalent to $\mathbf{H}^{1}$-norm, from (34) we get

[TABLE]

Then Sobolev’s imbedding theorem implies that

[TABLE]

with $\,1\leq p\leq 6$ for $d=3$ and $1\leq p<\infty$ for $d=2$.

Using Young’s inequality and the interpolation inequalities (26), we further prove

[TABLE]

Hence we have

[TABLE]

Consequently, we obtain

[TABLE]

Combining (36), (37) and (38), we complete the proof of (35).       ∎

4 The error estimates

In this section, we will give the proof of Theorem 1. Let $(I_{h}\Psi,{\bm{\pi}}_{h}\mathbf{A})$ denote the interpolation functions of $(\Psi,\mathbf{A})$ in $\mathcal{V}_{h}^{r}\times\mathbf{V}_{h}^{r}$. Set $e_{\Psi}=I_{h}\Psi-\Psi$, $e_{\mathbf{A}}={\bm{\pi}}_{h}\mathbf{A}-\mathbf{A}$. By applying the interpolation error estimates (14) and the regularity assumptions (21), we have

[TABLE]

where $C$ is a constant independent of $h$.

For convenience, we give the following identities, which will be used frequently in the sequel.

[TABLE]

Let $\theta_{\Psi}^{k}=\Psi_{h}^{k}-I_{h}\Psi^{k}$, $\theta_{\mathbf{A}}^{k}=\mathbf{A}_{h}^{k}-{\bm{\pi}}_{h}\mathbf{A}^{k}$. By using the error estimates of the interpolation operators, we only need to estimate $\theta_{\Psi}^{k}$ and $\theta_{\mathbf{A}}^{k}$. Subtracting (19) from (20), we get the following equations for $\theta_{\Psi}^{k}$ and $\theta_{\mathbf{A}}^{k}$:

[TABLE]

and

[TABLE]

where $\overline{\mathbf{A}}_{h}^{k}$, $\overline{\Psi}_{h}^{k}$, $\widetilde{\theta^{k}_{\mathbf{A}}}$ and $\widetilde{{\bm{\pi}}_{h}\mathbf{A}^{k}}$ are similarly given in (11).

Now we briefly describe the outline of the proof of (23). First, we take $\varphi=\overline{\theta}_{\Psi}^{k}$ in (41) and obtain the estimate of $\|\theta_{\Psi}^{M}\|_{\mathcal{L}^{2}}$. Second, we choose $\varphi=\theta_{\Psi}^{k}-\theta_{\Psi}^{k-1}$ in (41) and derive the energy-norm estimate for $\theta^{M}_{\Psi}$. Finally, let $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2\Delta t}(\theta_{\mathbf{A}}^{k}-\theta_{\mathbf{A}}^{k-2})$ in (42) and acquire the estimate involving $\theta_{\mathbf{A}}^{k}$. Combining the above three estimates, we will complete the proof of (23).

4.1 Estimates for (41)

To begin with, choosing $\varphi=\overline{\theta}_{\Psi}^{k}$, $k=1,2,\cdots,M$ as the test function in (41), we get

[TABLE]

where

[TABLE]

Using the error estimates (39) for the interpolation operator $I_{h}$ and the regularity of $\Psi$ in (21), it is easy to see that

[TABLE]

We observe that

[TABLE]

and

[TABLE]

It follows from (39), (45) and (46) that

[TABLE]

Notice that

[TABLE]

By using (24)-(26) and (29), we prove

[TABLE]

Taking the imaginary part of (43), we have

[TABLE]

and therefore

[TABLE]

Here we have used the fact $\|\theta^{0}_{\Psi}\|^{2}_{\mathcal{L}^{2}}\leq Ch^{2r+2}$.

To proceed further, we take $\varphi=\Delta t\partial{\theta_{\Psi}^{k}}=\theta_{\Psi}^{k}-\theta_{\Psi}^{k-1}$, $k=1,2,\cdots,M$ in (41), to find

[TABLE]

where

[TABLE]

By virtue of (40), we get

[TABLE]

It follows from (39) and (53) that

[TABLE]

To estimate the term $J_{2}^{(k)}$, we rewrite it as

[TABLE]

By applying a standard argument, we find that

[TABLE]

We recall (29) and rewrite $J_{3}^{(k)}$ as follows:

[TABLE]

By employing (39), (40), the regularity assumption (21) and Young’s inequality, we can prove the following estimate of $\sum_{k=1}^{M}{J_{3}^{(k)}}$ .

[TABLE]

The proof is standard but tedious. Due to space limitations, we omit it here.

In order to estimate $\sum\limits_{k=1}^{M}J_{4}^{(k)}$, we rewrite $J_{4}^{(k)}$ in the following form:

[TABLE]

By applying (35) and (40), we deduce

[TABLE]

In order to estimate $|\sum\limits_{k=1}^{M}J_{4}^{(k),3}|$, we rewrite it as follows.

[TABLE]

Note that

[TABLE]

By applying the Young’s inequality and (35), we can estimate the first two terms on the right side of (61) as follows.

[TABLE]

From (34) and (35), we further deduce

[TABLE]

where we have used the fact:

[TABLE]

Hence we get the following estimate:

[TABLE]

Employing (39) and integrating by parts, we discover

[TABLE]

By using the Young’s inequality, we can estimate the first three terms on the right side of (66) as follows:

[TABLE]

Using (39), the last term on the right side of (66) can be estimated by

[TABLE]

Hence we get

[TABLE]

Reasoning as before, we can estimate $Q_{3}$ as follows:

[TABLE]

Combining (65), (69) and (70) implies

[TABLE]

It follows from (58), (59) and (71) that

[TABLE]

Now take the real part of (52), and we get

[TABLE]

Similarly to (30), we have

[TABLE]

Substituting (74) into (73) and summing over $k=1,2,\cdots,M$, we get

[TABLE]

Combining (54), (55), (57) and (72) implies

[TABLE]

By employing Theorem 3, we discover

[TABLE]

Setting

[TABLE]

and applying (76) and (77), we obtain

[TABLE]

Arguing as in the proof of Theorem 3, we discover

[TABLE]

and thus

[TABLE]

Substituting (81) into (79) , we obtain

[TABLE]

Multiplying (51) with $(C+1)$ and adding to (82), we end up with

[TABLE]

4.2 Estimates for (42)

Setting

[TABLE]

we rewrite (42) as follows:

[TABLE]

We first estimate $\sum\limits_{k=1}^{M}K_{1}^{(k)}$, $\sum\limits_{k=1}^{M}K_{3}^{(k)}$ and $\sum\limits_{k=1}^{M}K_{4}^{(k)}$. Under the regularity assumption of $\mathbf{A}$ in (21), we have

[TABLE]

By applying the regularity assumption (21), the interpolation error estimates (39) and Theorem 3, it is easy to deduce

[TABLE]

In order to estimate $\sum\limits_{k=1}^{M}K_{4}^{(k)}$, we rewrite $K_{4}^{(k)}$ in the following form:

[TABLE]

We observe that

[TABLE]

We obtain from (39) and (89) that

[TABLE]

Similarly, from (39) and (89), we deduce

[TABLE]

Combining (88), (90) and (91) gives

[TABLE]

From (86), (87) and (92), we obtain

[TABLE]

Now taking $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2\Delta t}(\theta_{\mathbf{A}}^{k}-\theta_{\mathbf{A}}^{k-2})=\frac{\displaystyle 1}{\displaystyle 2}(\partial\theta_{\mathbf{A}}^{k}+\partial\theta_{\mathbf{A}}^{k-1})$ in (85), we find

[TABLE]

Note that

[TABLE]

Thus we have

[TABLE]

Here we have used the definition of $J_{5}^{(k)}$ in (78).

Multiplying (94) by $\Delta t$, and using (93) and (95), we obtain

[TABLE]

Since $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2}(\partial\theta_{\mathbf{A}}^{k}+\partial\theta_{\mathbf{A}}^{k-1})$, by applying (40) and the Young’s inequality, we get

[TABLE]

Combining (83), (96) and (97) implies

[TABLE]

By applying the discrete Gronwall’s inequality, we have

[TABLE]

Combine (99) with the interpolation error estimates (39) and we can complete the proof of Theorem 1.

5 Numerical tests

To validate the developed algorithm and to confirm the theoretical analysis reported in this paper, we present numerical simulations for the following case studies.

** Example 5.1****.**

We consider the Maxwell-Schrödinger system (2), where the initial-boundary conditions are as follows:

[TABLE]

Here we take $\Omega=(0,1)^{3}$, $V_{0}=0$, $T=0.5$, and the time step $\Delta t=0.0025$. Note that the initial wave function $\Psi(\mathbf{x},0)$ is the eigenfunction of the stationary Schrödinger’s equation.

The numerical results are displayed in Fig. 5.1.

** Remark 5.1****.**

The numerical results illustrated in Fig. 5.1 clearly show that the change of $\rho(\mathbf{x},t)$ is smooth with respect to $t$ and the assumption on which the modified Maxwell-Schrödinger equations are based is valid in this case.

** Example 5.2****.**

We consider the modified Maxwell–Schrödinger’s equations as follows:

[TABLE]

where the initial-boundary conditions are given in (8)-(9). We take $\Omega=(0,1)^{3}$, $T=4$, $\gamma=1$ and $V_{0}=5$. The exact solution $(\Psi,\mathbf{A})$ of (100) is given by

[TABLE]

[TABLE]

The functions $f$ and $\mathbf{g}$ in (100) are chosen correspondingly to the exact solution $(\Psi,\mathbf{A})$.

A uniform tetrahedral partition is generated with $M+1$ nodes in each direction and $6M^{3}$ elements in total. We solve the system(100) by the proposed Crank-Nicolson Galerkin finite element scheme (20) with linear elements and quadratic elements, respectively. To confirm our error analysis, we take $\Delta t=h^{\frac{1}{2}}$ for the linear element method and $\Delta t=h$ for the quadratic element method respectively. Numerical results for the linear element method and the quadratic element method at time $t=1.0,2.0,3.0,4.0$ are listed in Tables 1 and 2, respectively.

** Remark 5.2****.**

Numerical results in Tables 1 and 2 are in good agreement with the theoretical analysis, see Theorem 1.

** Example 5.3****.**

We consider the following modified Maxwell–Schrödinger’s equations

[TABLE]

with

[TABLE]

In this example we take $\Omega=(0,1)^{3}$, $T=10.0$ , $\gamma=1$, $V_{0}=0$. Using the mesh in Example 5.2 with M = 50, we solve the system (101) by the proposed Crank-Nicolson Galerkin finite element scheme (20) with linear elements. The time step $\Delta t=0.0025$.

In Fig. 5.2 we display the numerical results of $\rho=|\Psi|^{2}$ on the line $x_{1}=x_{2}=x_{3}$ and on the intersection $x_{1}=0.5$ at time $t=2.0,4.0,6.0,8.0,10.0$.

In Fig. 5.3 we plot the numerical results of $\rho=|\Psi|^{2}$ on the intersections $x_{1}=0.4,\,x_{2}=0.4\,$ and $x_{3}=0.4$ at time $t=10.0$, respectively.

Althouth the Maxwell-Schrödinger system and the time-dependent Ginzburg-Landau equations are somehow formally similar, they describe the different physical phenomenons. The time-dependent Ginzburg-Landau equations describe the vortex dynamics of superconductor [8, 9, 10] while the Maxwell-Schrödinger equations describe the wave packet dynamics of an electron. As can be seen in Fig. 5.2-5.3, the wave packet of the electron is located at the center of the computational domain at first and the external and its self-induced electromagnetic fields cause the motion of the wave packet. Unlike the time-dependent Ginzburg-Landau equations, no stable state is observed in our computation.

6 Conclusions

We have presented the optimal $H^{1}$ error estimates of a Crank-Nicolson Galerkin finite element method for the modified Maxwell-Schrödinger equations , which are derived from the original equations under some assumptions. The techniques used in this paper may also be applied to other nonlinear PDEs, such as the Ginzburg-Landau equations. The original Maxwell-Schrödinger system is challenging and difficult to perform numerical computation and theoretical analysis. Our work can serve as an elementary attempt for the numerical analysis of this system. We will study the original system in a further work using the mixed finite element method.