Error estimates of the Crank-Nicolson Galerkin method for the time-dependent Maxwell-Schr\"{o}dinger equations under the Lorentz gauge

TL;DR

This paper develops and analyzes an alternating Crank-Nicolson finite element method for solving the time-dependent Maxwell-Schrödinger equations under the Lorentz gauge, providing optimal error estimates and confirming results through numerical experiments.

Contribution

It introduces a new numerical scheme with proven optimal error bounds for Maxwell-Schrödinger equations under the Lorentz gauge.

Findings

Optimal error estimates are established for the proposed method.

Numerical experiments confirm the theoretical convergence rates.

Abstract

In this paper we study the numerical method and the convergence for solving the time-dependent Maxwell-Schr\"{o}dinger equations under the Lorentz gauge. An alternating Crank-Nicolson finite element method for solving the problem is presented and the optimal error estimate for the numerical algorithm is obtained by a mathematical inductive method. Numerical examples 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	4.9855e-01	2.3894e-01	1.1887e-01	1.03
2.0	6.7234e-01	3.2057e-01	1.7574e-01	0.97
3.0	4.6375e-01	2.2119e-01	1.0373e-01	1.08
4.0	5.8205e-01	3.0487e-01	1.4287e-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	3.0718e-01	1.4419e-01	7.4104e-02	1.02
2.0	4.3713e-01	2.1289e-01	1.1430e-01	0.97
3.0	3.0004e-01	1.4202e-01	6.9620e-02	1.05
4.0	2.2543e-01	1.2316e-01	6.0316e-02	0.95
		${\Vert {\varphi }_h^k-\varphi (\cdot ,t_k)\Vert }_{H^1}$
t	M=25	M=50	M=100	Order
1.0	1.9412e-01	8.6435e-02	4.3191e-02	1.07
2.0	1.2473e-01	7.0233e-02	3.3148e-02	0.96
3.0	1.1031e-01	6.1394e-02	2.9217e-02	0.96
4.0	8.7695e-02	4.1380e-02	2.1815e-02	1.00

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	6.5062e-02	1.6679e-02	4.1805e-03	1.98
2.0	5.5496e-02	1.5018e-02	3.7016e-03	1.95
3.0	4.1903e-02	1.3291e-02	2.4825e-03	2.04
4.0	3.4765e-02	8.4562e-03	2.0952e-03	2.03
		${\Vert {\psi }_h^k-\psi (\cdot ,t_k)\Vert }_{{\mathscr{H}}^1}$
t	M=25	M=50	M=100	Order
1.0	1.1930e-02	3.2391e-03	8.0967e-04	1.94
2.0	1.0237e-02	2.9786e-03	5.9668e-04	2.05
3.0	2.2342e-02	4.8413e-03	1.3683e-03	2.01
4.0	1.3203e-02	3.1681e-03	7.8548e-04	2.04
		${\Vert {\varphi }_h^k-\varphi (\cdot ,t_k)\Vert }_{H^1}$
t	M=25	M=50	M=100	Order
1.0	3.1070e-02	7.8114e-03	1.8648e-03	2.03
2.0	2.7189e-02	8.2275e-03	1.8301e-03	1.95
3.0	3.1162e-02	6.9836e-03	1.7682e-03	2.07
4.0	2.5343e-02	6.4257e-03	1.6146e-03	1.99

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.$$

$$\nabla \cdot \mathbf{A}+\frac{\partial \varphi }{\partial t}=0.$$

$$\nabla \cdot \mathbf{A}=0.$$

$$\varphi =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+\phi\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})-\nabla(\nabla\cdot\mathbf{A})+\frac{\mathrm{i}}{2}\big{(}\psi^{*}\nabla{\psi}-\psi\nabla{\psi}^{*}\big{)}}\\[5.69054pt] {\displaystyle\qquad\qquad+|\psi|^{2}\mathbf{A}=0,\,\,\quad(\mathbf{x},t)\in\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\frac{\partial^{2}\phi}{\partial t^{2}}-\Delta\phi=|\psi|^{2},\,\,(\mathbf{x},t)\in\Omega\times(0,T)}.\end{array}\right.$$

$$\psi (\mathbf{x},t)=0,\varphi (\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),$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\psi(\mathbf{x},0)=\psi_{0}(\mathbf{x}),\quad\phi(\mathbf{x},0)=\phi_{0}(\mathbf{x}),\quad\phi_{t}(\mathbf{x},0)=\phi_{1}(\mathbf{x}),}\\[5.69054pt] {\displaystyle\qquad\qquad\mathbf{A}(\mathbf{x},0)=\mathbf{A}_{0}(\mathbf{x}),\quad\mathbf{A}_{t}(\mathbf{x},0)=\mathbf{A}_{1}(\mathbf{x}),}\end{array}$$

$$\nabla \cdot {\mathbf{A}}_0+{\varphi }_1=0,\nabla \cdot {\mathbf{A}}_1+\Delta {\varphi }_0+|{\psi }_0{|}^2=0.$$

$${\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 }_{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)+(\phi\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})+(\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})=0,\qquad\qquad\qquad\qquad\qquad\forall\mathbf{v}\in\mathbf{H}^{1}_{t}(\Omega),}\\[5.69054pt] {\displaystyle(\frac{\partial^{2}\phi}{\partial t^{2}},\eta)+(\nabla\phi,\nabla\eta)=(|\psi|^{2},\eta),\quad\quad\forall\eta\in H_{0}^{1}(\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}$$

$$\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,$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\psi_{h}^{0}=R_{h}\psi_{0},\quad\mathbf{A}_{h}^{0}={\pi}_{h}\mathbf{A}_{0},\quad\mathbf{A}_{h}^{-1}=\mathbf{A}_{h}^{0}-\Delta t{\pi}_{h}\mathbf{A}_{1},}\\[5.69054pt] {\displaystyle\qquad\qquad\phi_{h}^{0}=I_{h}\phi_{0},\quad\phi_{h}^{-1}=\phi_{h}^{0}-\Delta tI_{h}\psi_{1},}\end{array}$$

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\mathrm{i}({\partial}\psi_{h}^{k},\varphi)+\frac{1}{2}\left((\mathrm{i}\nabla+\overline{\mathbf{A}}^{k}_{h})\overline{\psi}_{h}^{k},(\mathrm{i}\nabla+\overline{\mathbf{A}}_{h}^{k})\varphi\right)+\left((V_{0}+\overline{\phi}_{h}^{k})\overline{\psi}_{h}^{k},\varphi\right)=0,\,\,\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+(\nabla\cdot\widetilde{\mathbf{A}}_{h}^{k},\nabla\cdot\mathbf{v})+\big{(}|\psi_{h}^{k-1}|^{2}\widetilde{\mathbf{A}}_{h}^{k},\mathbf{v}\big{)}=0,\quad\forall\mathbf{v}\in\mathbf{V}^{r}_{h},}\\[5.69054pt] {\displaystyle(\partial^{2}\phi^{k}_{h},\eta)+(\nabla\widetilde{\phi}_{h}^{k},\nabla\eta)=(|\psi_{h}^{k-1}|^{2},\eta),\quad\forall\eta\in V_{h}^{r}.}\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})=(\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)+(\phi\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})=0,\quad\forall\mathbf{v}\in\mathbf{H}^{1}_{t}(\Omega),}\\[5.69054pt] {\displaystyle(\frac{\partial^{2}\phi}{\partial t^{2}},\eta)+(\nabla\phi,\nabla\eta)=(|\psi|^{2},\eta),\quad\quad\forall\eta\in H_{0}^{1}(\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)+(\overline{\phi}_{h}^{k}\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)+\big{(}|\psi_{h}^{k-1}|^{2}\widetilde{\mathbf{A}}_{h}^{k},\mathbf{v}\big{)}=0,\quad\forall\mathbf{v}\in\mathbf{V}^{r}_{h},}\\[5.69054pt] {\displaystyle(\partial^{2}\phi^{k}_{h},\eta)+(\nabla\widetilde{\phi}_{h}^{k},\nabla\eta)=(|\psi_{h}^{k-1}|^{2},\eta),\quad\forall\eta\in V_{h}^{r}.}\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));}\\[5.69054pt] {\displaystyle\phi,\phi_{t},\phi_{tt}\in{L}^{\infty}(0,T;{H}^{r+1}(\Omega)),\phi_{ttt}\in{L}^{\infty}(0,T;{H}^{1}(\Omega)),}\\[5.69054pt] {\displaystyle\phi_{tttt}\in L^{2}(0,T;{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);\,\,\,\phi_{0},\phi_{1}\in H^{r+1}(\Omega)\cap H_{0}^{1}(\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\qquad\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}}\\[5.69054pt] {\displaystyle\qquad+\|\phi_{h}^{k}-\phi^{k}\|^{2}_{L^{2}(\Omega)}+\|\nabla(\phi_{h}^{k}-\phi^{k})\|^{2}_{\mathbf{L}^{2}(\Omega)}\big{\}}\leq C\big{\{}h^{2r}+(\Delta t)^{4}\big{\}},\,\,r\geq 1,}\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 ),$$

$$\Vert u{\Vert }_{L^3}^2\le \Vert u{\Vert }_{L^2}\Vert u{\Vert }_{L^6},$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\sum_{k=1}^{M}{(a_{k}-a_{k-1})b_{k}}=a_{M}b_{M}-a_{0}b_{1}-\sum_{k=1}^{M-1}{a_{k}(b_{k+1}-b_{k})},}\\[5.69054pt] {\displaystyle\sum_{k=1}^{M}{(a_{k}-a_{k-1})b_{k}}=a_{M}b_{M}-a_{0}b_{0}-\sum_{k=1}^{M}{a_{k-1}(b_{k}-b_{k-1})}.}\end{array}$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\|e_{\psi}\|_{\mathcal{L}^{2}}+h\|e_{\psi}\|_{\mathcal{H}^{1}}\leq Ch^{r+1},\quad\|e_{\mathbf{A}}\|_{\mathbf{L}^{2}}+h\|e_{\mathbf{A}}\|_{\mathbf{H}^{1}}\leq Ch^{r+1},}\\[5.69054pt] {\displaystyle\|e_{\phi}\|_{{L}^{2}}+h\|e_{\phi}\|_{{H}^{1}}\leq Ch^{r+1},\,\,\,\|R_{h}\psi\|_{\mathcal{L}^{\infty}}+\|{\pi}_{h}\mathbf{A}\|_{\mathbf{L}^{\infty}}+\|I_{h}\phi\|_{{L}^{\infty}}+\|\nabla R_{h}\psi\|_{\mathbf{L}^{3}}\leq C,}\end{array}$$

$$\Vert {\psi }_h^k{\Vert }_{{\mathcal{L}}^2}^2={\Vert {\psi }_h^0\Vert }_{{\mathcal{L}}^2}^2,k=1,2\cdots ,M.$$

$$\begin{array}[]{@{}l@{}}{\displaystyle B(\mathbf{A};\psi,\varphi)-B(\widetilde{\mathbf{A}};\psi,\varphi)=\left((\mathbf{A}+\widetilde{\mathbf{A}})\psi\varphi^{*},\mathbf{A}-\widetilde{\mathbf{A}}\right)+2(f(\psi,\varphi),\mathbf{A}-\widetilde{\mathbf{A}}),}\\[5.69054pt] {\displaystyle\mathrm{Re}\left[B\left(\overline{\mathbf{A}}^{k};\overline{\psi}^{k},\partial{\psi}^{k}\right)\right]=-\left(\frac{1}{2}(\overline{\mathbf{A}}^{k}+\overline{\mathbf{A}}^{k-1})|{\psi}^{k-1}|^{2},\frac{1}{2}(\partial\mathbf{A}^{k}+\partial\mathbf{A}^{k-1})\right)}\\[5.69054pt] {\displaystyle\quad-\left(f({\psi}^{k-1},{\psi}^{k-1}),\frac{1}{2}(\partial\mathbf{A}^{k}+\partial\mathbf{A}^{k-1})\right)+\frac{1}{2}\partial B(\overline{\mathbf{A}}^{k};{\psi}^{k},{\psi}^{k}).}\end{array}$$

$$\Vert {\theta }_{\psi }^k{\Vert }_{{\mathcal{H}}^1}^2+\Vert \partial {\theta }_{\mathbf{A}}^k{\Vert }_{{\mathbf{L}}^2}^2+\Vert {\theta }_{\mathbf{A}}^k{\Vert }_{{\mathbf{H}}^1}^2+\Vert \partial {\theta }_{\varphi }^k{\Vert }_{L^2}^2+\Vert {\theta }_{\varphi }^k{\Vert }_{H^1}^2\le C_{\ast }\left\{h^{2r}+(\Delta t{)}^4\right\},$$

$$\Vert {\theta }_{\psi }^0{\Vert }_{{\mathcal{H}}^1}^2+\Vert \partial {\theta }_{\mathbf{A}}^0{\Vert }_{{\mathbf{L}}^2}^2+\Vert {\theta }_{\mathbf{A}}^0{\Vert }_{{\mathbf{H}}^1}^2+\Vert \partial {\theta }_{\varphi }^0{\Vert }_{L^2}^2+\Vert {\theta }_{\varphi }^0{\Vert }_{H^1}^2\le C_0\left\{h^{2r}+(\Delta t{)}^4\right\}.$$

Full text

1 Introduction

Light-matter interaction at nanoscale is a central topic in the study of optical properties of nanophotonic systems, for example, metallic nanostructures and quantum dots. In view of practical numerical simulation, a semiclassic model is often used for modelling light-matter interaction. The basic idea is to use the classical Maxwell’s equations for the electromagnetic field and the Schrödinger equation for the matter. In this paper, we study the following Maxwell-Schrödinger coupled system, which describes the interaction between an electron and its self-consistent generated and external electromagnetic fields.

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$ is a bounded Lipschitz polygonal convex domain in $\mathbb{R}^{2}$ (or a bounded Lipschitz polyhedron convex domain in $\mathbb{R}^{3}$), $\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) lack uniqueness. 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). To obtain mathematically well-posed equations, some extra constraint, commonly known as gauge choice, is often enforced on the solutions of the Maxwell-Schrödinger equations. The most common gauges are listed below.

(i) The Lorentz gauge

[TABLE]

(ii) The Coulomb gauge

[TABLE]

(iii) The temporal gauge

[TABLE]

For simplicity, we employ the atomic units, i.e. $\hbar=m=q=1$, and assume that $\epsilon=\mu=1$ without loss of generality.

In this paper, we consider the time-dependent Maxwell-Schrödinger equations under the Lorentz gauge as follows:

[TABLE]

The boundary conditions are

[TABLE]

and the initial conditions are

[TABLE]

where $\phi_{t}$ and $\mathbf{A}_{t}$ denote the derivative of $\phi$ and $\mathbf{A}$ with respect to the time $t$ respectively, $\mathbf{n}=(n_{1},n_{2},n_{3})$ is the outward unit normal to the boundary $\partial\Omega$.

The gauge choice and the equations (2) impose the following constraints on the initial datas:

[TABLE]

Remark 1.1**.**

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$ and $\phi(\mathbf{x},t)=0$ on $\partial\Omega$ are direct results of the perfect conductive boundary (PEC). The boundary condition $\nabla\cdot\mathbf{A}=0$ on $\partial\Omega$ can be deduced from the boundary condition of $\,\phi$ and the gauge choice. As for the determination of the boundary conditions for the vector potential $\mathbf{A}$ and the scalar potential $\phi$, we refer to Chew (2014).

The local and global well-posedness of solutions on all of $\mathbb{R}^{3}$ for the time-dependent Maxwell-Schrödinger equations (1) have been studied in, for example, Ginibre & Velo (2003), Guo et al. (1995), Nakamitsu & Tsutsumi (1986), Nakamura & Wada (2005), Nakamura & Wada (2007), Shimomura (2003) and Wada (2012). 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.

There are a number of results for the numerical methods for the coupled Maxwell-Schrödinger equations. We recall some interesting studies. Sui et al. (2007) proposed the finite-difference time-domain (FDTD) method to solve the Maxwell-Schrödinger equations in the simulation of electron tunneling problem. Pierantoni et al. (2008) studied a carbon nanotube between two metallic electrodes by solving the Maxwell-Schrödinger equations with the transmission line matrix (TLM)/FDTD hybrid method. Ahmed & Li (2012) used the FDTD method for the Maxwell-Schrödinger system to simulate plasmonics nanodevices. However, in the numerical studies listed above, the EM fields are all described by the Maxwell’s equations involving electric fields $\mathbf{E}$ and magnetic fields $\mathbf{H}$, instead of the $\mathbf{A}$-$\phi$ formulation. Recently, Ryu (2015) applied the FDTD scheme to discretize the Maxwell-Schrödinger equations (1) under the Lorentz gauge and to simulate the interaction between a single electron in an artificial atom and an incoming electromagnetic field. For more results on this topic, we refer to Ohnuki et al. (2013), Sato & Yabana (2014), Turati & Hao (2012) and the references therein.

There are few results on the finite element method (FEM) and its convergence analysis of the Maxwell-Schrödinger equations (1). In this paper we will present an alternating Crank-Nicolson finite element method for solving the problem (2)-(4), 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. Our work is motivated by Gao et al. (2014) in which Gao and his collaborators proposed a linearized Crank-Nicolson Galerkin method for the time-dependent Ginzburg-Landau equations and derived an optimal error estimate via a mathematical inductive method under the assumption that $h$ and $\Delta t$ are sufficiently small. Here $h$ and $\Delta t$ are the spatial mesh size and the time step, respectively. Compared to the time-dependent Ginzburg-Landau model, the error analysis of numerical schemes for the time-dependent Maxwell-Schrödinger system is much more difficult. The main difficulties and tricky parts in this paper are the estimates of the current term $\mathbf{J}_{q}$ and the error analysis for the wave function $\psi$. In particular we derive the energy-norm error estimates for the Schrödinger’s equation.

The remainder of this paper is organized as follows. In section 2, we introduce some notation and propose a decoupled alternating Crank-Nicolson scheme with the Galerkin finite element approximation for the Maxwell-Schrödinger equations (2)-(4). The proof of the main theorem (see Theorem 2.1) in this paper will be given in section  3. Finally, some numerical tests are carried out to validate the theoretical results in this paper.

Throughout the paper the Einstein summation convention on repeated indices is adopted. By $C$ we shall denote a positive constant independent of the mesh size $h$ and the time step $\Delta t$ without distinction.

2 An alternating Crank-Nicolson Galerkin finite element scheme

In this section, we present a numerical scheme for the Maxwell-Schrödinger equations (2)-(4) using Galerkin finite element method in space and a decoupled alternating 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 Girault & Raviart (1986).

The weak formulation of the Maxwell-Schrödinger system (2)-(4) can be specified as follows: find $(\psi,\mathbf{A},\phi)\in\mathcal{H}_{0}^{1}(\Omega)\times\mathbf{H}_{t}^{1}(\Omega)\times H_{0}^{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)$, $\phi_{0}\in H_{0}^{1}(\Omega)$, $\mathbf{A}_{t}(\cdot,0)\in\mathbf{L}^{2}(\Omega)$ and $\phi_{t}(\cdot,0)\in 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}=\{e\}$ 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}_{e\in\mathcal{T}_{h}}\{diam(e)\}$. For a given partition $\mathcal{T}_{h}$, let $\mathcal{V}_{h}^{r}$ , $\mathbf{V}_{h}^{r}$ and $V_{h}^{r}$ denote the corresponding $r$-th order finite element subspaces of $\mathcal{H}_{0}^{1}(\Omega)$ , $\mathbf{H}^{1}_{t}(\Omega)$ and $H_{0}^{1}(\Omega)$, respectively. Let $R_{h}$, ${\pi}_{h}$ and $I_{h}$ be the conventional point-wise interpolation operators on $\mathcal{V}_{h}^{r}$, $\mathbf{V}_{h}^{r}$ and $V_{h}^{r}$, respectively.

For convenience, assume that the function $\mathbf{A}$ and $\phi$ 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. $\phi^{-1}$ can be approximated in a similar way.

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

[TABLE]

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

[TABLE]

Note that $\overline{\mathbf{A}}_{h}^{k}$, $\overline{\phi}_{h}^{k}$, $\overline{\psi}_{h}^{k}$, $\widetilde{\mathbf{A}}_{h}^{k}$ and $\widetilde{\phi}_{h}^{k}$ are defined in (7), and $(\psi_{h}^{k-1})^{\ast}$ denotes the complex conjugate of $\psi_{h}^{k-1}$. For convenience, we define the following bilinear forms:

[TABLE]

Then the variational form of the Maxwell-Schrödinger system (6) and its discrete system (7) can be reformulated as follows:

[TABLE]

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

[TABLE]

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

[TABLE]

For the initial conditions $(\psi_{0},\mathbf{A}_{0},\mathbf{A}_{1},\phi_{0},\phi_{1})$ , we assume that

[TABLE]

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

Theorem 2.1**.**

Suppose that the Maxwell-Schrödinger coupled system (12) has a unique solution $(\psi,\mathbf{A},\phi)$ satisfying (14) and (15). Let $(\psi_{h}^{k},\mathbf{A}_{h}^{k},\phi_{h}^{k})$ be the fully discrete numerical solution of $(\psi,\mathbf{A},\phi)$ defined in (13). Then there exist two positive constants $h_{0}>0$ and $\Delta t_{0}>0$, such that when $h<h_{0}$, $\Delta t<\Delta t_{0}$, we have the following error estimates:

[TABLE]

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

3 The proof of Theorem 2.1

In this section, we will give the proof of Theorem 2.1.

3.1 Preliminaries

For convenience, we list some imbedding inequalities and interpolation inequalities in Sobolev spaces (see, e.g., Ladyzhenskaya et al. (1968) and Girault & Raviart (1986)), 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)}$.

The following identities will be used frequently in this paper.

[TABLE]

Let $(R_{h}\psi,{\pi}_{h}\mathbf{A},I_{h}\phi)$ denote the interpolation function of $(\psi,\mathbf{A},\phi)$ in $\mathcal{V}_{h}^{r}\times\mathbf{V}_{h}^{r}\times V_{h}^{r}$. Set $e_{\psi}=R_{h}\psi-\psi$, $e_{\mathbf{A}}={\pi}_{h}\mathbf{A}-\mathbf{A}$, $e_{\phi}=I_{h}\phi-\phi$. By applying standard finite element theory and the regualrity conditions in (14), we have

[TABLE]

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

The following lemmas will be useful in the proof of Theorem 2.1.

Lemma 3.1**.**

For the solution of (13), we have

[TABLE]

Lemma 3.2**.**

For $k=1,2\cdots,M$, the following identities hold for the bilinear functional $B(\mathbf{A};\psi,\varphi)$ defined in (11):

[TABLE]

Lemma 3.1 can be proved by choosing $\varphi=\overline{\psi}_{h}^{k}$ in $(\ref{eq:2-9})_{1}$ and taking the imaginary part. A direct calculation gives (23) in Lemma 3.2.

Let $\theta_{\psi}^{k}=\psi_{h}^{k}-R_{h}\psi^{k}$, $\theta_{\mathbf{A}}^{k}=\mathbf{A}_{h}^{k}-{\pi}_{h}\mathbf{A}^{k}$, $\theta_{\phi}^{k}=\phi_{h}^{k}-I_{h}\phi^{k}$. By using the error estimates of the interpolation operators (21), we only need to estimate $\theta_{\psi}^{k}$ , $\theta_{\mathbf{A}}^{k}$ and $\theta_{\phi}^{k}$. We will prove the following estimate:

[TABLE]

where $k=0,1,\cdots,M$.

By using the regularity assumption of the initial conditions (15) and the error estimates of the interpolation operators (21), we get

[TABLE]

We use the mathematical inductive method to show (24). By (25), if we require $C_{\ast}\geq C_{0}$, then (24) holds for $k=0$. We assume that (24) holds for $0\leq k\leq m-1$. In the rest of this section, we will find $C_{\ast}$ such that (24) holds for $k\leq m$, where $C_{\ast}$ is independent of $k\,$, $h\,$, $\Delta t$.

Subtracting (12) from (13), we obtain the following equations for $\theta_{\mathbf{A}}^{k}$, $\theta_{\phi}^{k}$ and $\theta_{\psi}^{k}$ :

[TABLE]

[TABLE]

[TABLE]

The key steps of the proof of (24) are now briefly described. In order to find $C_{\ast}$ and to show that (24) holds for $k=m$, we first take $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2\Delta t}(\theta_{\mathbf{A}}^{k}-\theta_{\mathbf{A}}^{k-2})$ in (26) and $\eta=\frac{\displaystyle 1}{\displaystyle 2\Delta t}(\theta_{\phi}^{k}-\theta_{\phi}^{k-2})$ in (27) and obtain the estimates of $\theta_{\mathbf{A}}^{m}$ and $\theta_{\phi}^{m}$. Then we choose $\varphi=\overline{\theta}_{\psi}^{k}$ in (28) and give the estimate of $\|\theta_{\psi}^{m}\|_{\mathcal{L}^{2}}$. Finally, we take $\varphi=\partial\theta_{\psi}^{k}$ in (28) and make use of the above estimates of $\theta_{\mathbf{A}}^{m}$ and $\theta_{\phi}^{m}$ to derive the energy-norm estimate for $\theta^{m}_{\psi}$. Using the above estimates, we can complete the proof of (24).

3.2 Estimates for (26)

If we set

[TABLE]

then we rewrite (26) as follows:

[TABLE]

We take $\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})=\overline{\partial\theta}_{\mathbf{A}}^{k}$ in (30) and get

[TABLE]

which leads to

[TABLE]

Here we have used the fact that

[TABLE]

Now we estimate $\sum\limits_{k=1}^{m}I_{j}^{k}(\overline{\partial\theta}_{\mathbf{A}}^{k})$, $j=1,2,3,4$. Under the regularity assumption of $\mathbf{A}$ in (14), we can prove

[TABLE]

We rewrite the term $\Delta t\sum_{k=1}^{m}I_{2}^{k}(\overline{\partial\theta}_{\mathbf{A}}^{k})$ as follows:

[TABLE]

Applying (20), the regularity assumption and Young’s inequality, we get

[TABLE]

and

[TABLE]

We thus have

[TABLE]

We observe that

[TABLE]

The first four terms in (38) can be estimated by a standard argument, i.e.

[TABLE]

We notice that

[TABLE]

and obtain

[TABLE]

By using the assumption of the induction, we have

[TABLE]

If we choose some sufficiently small $h$ and $\Delta t$ such that $CC^{\frac{1}{2}}_{\ast}\left\{(\Delta t)^{2}+h^{r}\right\}\leq 1$, then we get

[TABLE]

Combining (38)-(42) implies

[TABLE]

To estimate $\sum\limits_{k=1}^{m}I_{4}^{k}(\mathbf{v})$, we rewrite $I_{4}^{k}(\overline{\partial\theta}_{\mathbf{A}}^{k})$ as follows:

[TABLE]

Observing

[TABLE]

and applying (21), we get

[TABLE]

Using (21), we similarly prove

[TABLE]

Hence we have

[TABLE]

Substituting (33), (37), (43) and (48) into (32), we get

[TABLE]

We estimate the last term on the right side of (49). It follows from the definition of bilinear functional $f(\varphi,\psi)$ in (11) that

[TABLE]

For $\Delta t\leq h^{\frac{1}{2}}$, by the assumption of the induction and the inverse inequalities, we have

[TABLE]

We choose a sufficiently small $h>0$ such that $CC_{\ast}^{\frac{1}{2}}h^{\frac{1}{2}}\leq 1$ and obtain

[TABLE]

Consequently

[TABLE]

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

As for $h^{\frac{1}{2}}\leq\Delta t$, by the assumption of the induction, we discover

[TABLE]

Now choose a sufficiently small $\Delta t>0$ to find $2C_{\ast}^{\frac{1}{2}}\Delta t\leq 1$, in which case we have

[TABLE]

It follows that

[TABLE]

Combining (51), (52) and (49) implies

[TABLE]

We choose a sufficiently small $\Delta t>0$ such that $C\Delta t\leq\frac{1}{8}$ and find

[TABLE]

which leads to

[TABLE]

by the assumption of the induction.

Applying the assumption of the induction again and (55), for $k=1,2,\cdots,m$, we deduce

[TABLE]

when $\|\partial\theta_{\mathbf{A}}^{k}\|_{\mathbf{L}^{2}}+D^{\frac{1}{2}}({\theta}_{\mathbf{A}}^{k},{\theta}_{\mathbf{A}}^{k})\leq 1.$

3.3 Estimates for (27)

Setting

[TABLE]

we rewrite (27) as follows:

[TABLE]

We take $\eta=\frac{\displaystyle 1}{\displaystyle 2}(\partial\theta_{\phi}^{k}+\partial\theta_{\phi}^{k-1})=\overline{\partial\theta}_{\phi}^{k}$ in (58) and obtain

[TABLE]

Multiply (59) by $\Delta t$ and sum $k=1,2,\cdots,m$ to discover

[TABLE]

Applying the regularity assumption of $\phi$ in (14), we deduce

[TABLE]

$\sum_{k=1}^{m}J_{2}^{k}(\overline{\partial\theta}_{\phi}^{k})$ can be bounded by a argument similar to the estimate of $\sum_{k=1}^{m}I_{2}^{k}(\overline{\partial\theta}_{\mathbf{A}}^{k})$ in (35) and (36).

[TABLE]

Observing

[TABLE]

and using (40), we get

[TABLE]

Here we have used $\|\theta_{\psi}^{k-1}\|_{\mathcal{L}^{4}}\leq C\|\theta_{\psi}^{k-1}\|_{\mathcal{H}^{1}}\leq 1.$

Substituting (61), (62) and (63) into (60) yields

[TABLE]

Similarly to (56), for $k=1,2,\cdots,m$, we can prove

[TABLE]

where $C$ is a constant independent of $h$, $\Delta t$ and $C_{\ast}$.

3.4 Estimates for (28)

We rewrite (28) as follows:

[TABLE]

where

[TABLE]

Taking $\varphi=\overline{\theta}_{\psi}^{k}$ in (66), and observing the imaginary part of the above equation, we have

[TABLE]

It is obvious that

[TABLE]

Using the error estimates (21) for the interpolation operator $R_{h}$ and the regularity of $\psi$ in (14), we give the following estimate

[TABLE]

Note that

[TABLE]

and

[TABLE]

Hence we obtain

[TABLE]

To estimate $Q_{4}^{k}(\overline{\theta}_{\psi}^{k})$, we rewrite it as follows:

[TABLE]

It follows from Lemma 3.1, the regularity assumption (14), the properties of the interpolation operators and (72) that

[TABLE]

We observe that

[TABLE]

It follows from Lemma 3.2 that

[TABLE]

Using (56), we prove

[TABLE]

and thus

[TABLE]

From (67), summing over $k=1,2,\cdots,m$ and combining (68), (71), (73) and (77), we have

[TABLE]

To proceed further, we take $\varphi=\partial{\theta_{\psi}^{k}}=\frac{\displaystyle 1}{\displaystyle\Delta t}(\theta_{\psi}^{k}-\theta_{\psi}^{k-1})$, $k=1,2,\cdots,m$ in (26), to find

[TABLE]

We take the real part of (79) and use (23) to get

[TABLE]

which leads to

[TABLE]

It follows from (56) that

[TABLE]

Combining (51) and (52) gives

[TABLE]

Substituting (82) and (83) into (81), we have

[TABLE]

We now proceed to estimate \sum_{k=1}^{m}\mathrm{Re}\big{[}Q_{j}^{k}(\partial{\theta_{\psi}^{k}})\big{]}, $j=1,\cdots,5$. By virtue of (20), we get

[TABLE]

From (21) and (85), we deduce

[TABLE]

We observe that

[TABLE]

It is not difficult to prove

[TABLE]

We estimate $|\sum\limits_{k=1}^{m}\mathrm{Re}[J_{2}^{k,2}]|$ by

[TABLE]

Hence we have

[TABLE]

From the definition of the bilinear functional $B(\mathbf{A};\psi,\varphi)$ in (11), we rewrite $\Delta tQ_{3}^{k}(\partial{\theta_{\psi}^{k}})$ as follows:

[TABLE]

Analogous to the estimate of $\sum_{k=1}^{m}Q_{1}^{k}(\partial{\theta_{\psi}^{k}})$, we use (20), (21), the regularity assumption (14) and Young’s inequality to prove

[TABLE]

Due to space limitations, we omit the proof of (91).

We observe that

[TABLE]

and have

[TABLE]

The first two terms can be estimated as follows.

[TABLE]

By applying (65), the last term on the right hand side of (92) can be bounded as follows:

[TABLE]

Combining (92)-(94) gives

[TABLE]

$\Delta tQ_{5}^{k}(\partial{\theta_{\psi}^{k}})$ can be decomposed as follows:

[TABLE]

Following the lines of the proof of (95) and using (23), we prove

[TABLE]

To estimate $|\sum\limits_{k=1}^{m}R_{5}^{k,3}|$, we rewrite it as follows:

[TABLE]

Note that

[TABLE]

By applying Young’s inequality and (56), we can estimate the first two terms on the right side of (99) by

[TABLE]

Observing

[TABLE]

and using (56), we get

[TABLE]

From (99)-(102), we thus have

[TABLE]

From (98), using (20) and integrating by parts, we get

[TABLE]

Using (21) and Young’s inequality, we can estimate the first three terms on the right side of (104):

[TABLE]

The last term on the right side of (104) can be estimated by

[TABLE]

We thus get

[TABLE]

We similarly prove

[TABLE]

Adding (103), (107) and (108) together, we find

[TABLE]

It follows from (96), (97) and (109) that

[TABLE]

Substituting (86), (89), (91), (95) and (110) into (84), we get

[TABLE]

Since

[TABLE]

we have

[TABLE]

It follows from (56) and the interpolation inequality (19) that

[TABLE]

We thus obtain

[TABLE]

Multiplying (78) with $(C+1)$ and adding to (114), we get

[TABLE]

Adding (55) and (64) to (115), we end up with

[TABLE]

Now by applying the discrete Gronwall’s inequality and choosing a sufficiently small $\Delta t$ such that $C\Delta t\leq\frac{\displaystyle 1}{\displaystyle 2}$, we conclude

[TABLE]

If we take $C_{\ast}\geq C\exp\left(2TC\right)$, then (24) holds for $k=m$. By virtue of the mathematical induction, we complete the proof of (24). Furthermore, using the triangle inequality and the equality $\theta_{\mathbf{A}}^{M}=\Delta t\sum\limits_{k=1}^{M}{\partial\theta_{\mathbf{A}}^{k}}+\theta_{\mathbf{A}}^{0}$, we can complete the proof of Theorem 2.1.

4 Numerical testes

In this section, we present numerical experiments to illustrate the error estimates.

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

[TABLE]

where the initial-boundary conditions are given in (3)-(4).

Let $\Omega=(0,1)^{3}$, $T=4$ and $V_{0}=5$. The exact solution $(\psi,\mathbf{A},\phi)$ of (118) is defined by

[TABLE]

[TABLE]

[TABLE]

The functions $f(\mathbf{x},t)$, $\mathbf{g}(\mathbf{x},t)$ and $l(\mathbf{x},t)$ at the right hand side of (118) are chosen correspondingly to the exact solution $(\psi,\mathbf{A},\phi)$.

We partition a whole domain $\Omega$ into quasi-uniform tetrahedrons with $M+1$ nodes in each direction and $6M^{3}$ elements in total. The system (118) is solved by the proposed Crank-Nicolson Galerkin finite element scheme (13) with linear elements and quadratic elements, respectively. To confirm the convergence rate of the proposed method, we take $\Delta t=h^{\frac{1}{2}}$ for the linear element method and $\Delta t=h$ for the quadratic element method respectively. The numerical results for the linear element method and the quadratic element method at time $t=1.0,2.0,3.0,4.0$ are displayed in Tables 1 and 2, respectively. It clearly shows that they are in good agreement with the error estimates presented in Theorem 2.1.