Mathematical and numerical analysis of the time-dependent Maxwell--Schr\"{o}dinger Equations in the Coulomb gauge

TL;DR

This paper establishes the global existence of weak solutions for the time-dependent Maxwell--Schrödinger equations in the Coulomb gauge, proposes an energy-conserving finite element scheme, and confirms its accuracy through numerical tests.

Contribution

It introduces a novel energy-conserving finite element method for the Maxwell--Schrödinger system and provides rigorous error analysis without time-step restrictions.

Findings

Proved global existence of weak solutions.

Developed an energy-conserving finite element scheme.

Achieved optimal error estimates and validated with numerical results.

Abstract

In this paper, we consider the initial-boundary value problem for the time-dependent Maxwell--Schr\"{o}dinger equations in the Coulomb gauge. We first prove the global existence of weak solutions to the equations. Next we propose an energy-conserving fully discrete finite element scheme for the system and prove the existence and uniqueness of solutions to the discrete system. The optimal error estimates for the numerical scheme without any time-step restrictions are then derived. Numerical results are provided to support our theoretical analysis.

Tables

Table 1. Table 1: H 1 superscript 𝐻 1 H^{1} error of linear FEM with h = 1 N ℎ 1 𝑁 h=\frac{1}{N} and τ = h 1 2 𝜏 superscript ℎ 1 2 \tau=h^{\frac{1}{2}} .

	${\Vert {\Psi }_h^M-\Psi (\cdot ,2)\Vert }_{{\mathscr{H}}^1}$	${\Vert {𝐀}_h^M-𝐀(\cdot ,2)\Vert }_{{𝐇}^1}$	${\Vert {\varphi }_h^M-\varphi (\cdot ,2)\Vert }_{H^1}$
N=25	3.3481e-01	2.8513e-01	4.0193e-02
N=50	1.5649e-01	1.1960e-01	1.6210e-02
N=100	6.9587e-02	4.5826e-02	7.0549e-03
order	1.13	1.33	1.25

Table 2. Table 2: H 1 superscript 𝐻 1 H^{1} error of quadratic FEM with h = 1 N ℎ 1 𝑁 h=\frac{1}{N} and τ = h 𝜏 ℎ \tau=h .

	${\Vert {\Psi }_h^M-\Psi (\cdot ,2)\Vert }_{{\mathscr{H}}^1}$	${\Vert {𝐀}_h^M-𝐀(\cdot ,2)\Vert }_{{𝐇}^1}$	${\Vert {\varphi }_h^M-\varphi (\cdot ,2)\Vert }_{H^1}$
N=25	2.3605e-02	1.5911e-02	6.8377e-03
N=50	6.1967e-03	4.1329e-03	1.6405e-03
N=100	1.5204e-03	9.4441e-04	3.9781e-04
order	1.97	2.04	2.05

Equations

$$\begin{array}[]{lll}{\displaystyle\mathrm{i}\hbar\frac{\partial\Psi}{\partial t}=\Big{\{}\frac{1}{2m}\left[\mathrm{i}\hbar\nabla+q\mathbf{A}\right]^{2}+q\phi+V\Big{\}}\Psi\quad{\rm in}\,\,\,\Omega_{T},}\\[5.69054pt] \end{array}$$

$$\mathbf{E}=-\nabla \varphi -\frac{\partial \mathbf{A}}{\partial t},\mathbf{B}=\nabla \times \mathbf{A},$$

$$\begin{array}[]{lll}{\displaystyle\nabla\times\mathbf{E}+\frac{\partial\mathbf{B}}{\partial t}=0,\;\quad\quad\nabla\cdot\mathbf{B}=0,}\\[5.69054pt] {\displaystyle\frac{1}{{\mu}}\nabla\times\mathbf{B}-\epsilon\frac{\partial\mathbf{E}}{\partial t}=\mathbf{J},\quad\nabla\cdot(\epsilon\mathbf{E})=\rho.}\end{array}$$

$$\rho=q|\Psi|^{2},\qquad\mathbf{J}=-\frac{\mathrm{i}q\hbar}{2m}\big{(}\Psi^{\ast}\nabla{\Psi}-\Psi\nabla{\Psi}^{\ast}\big{)}-\frac{|q|^{2}}{m}|\Psi|^{2}\mathbf{A},$$

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

$$\begin{array}[]{@{}l@{}}{\displaystyle\mathcal{E}(t)=\int_{\Omega}\big{(}\frac{1}{2}\big{|}\left(i\nabla+q\mathbf{A}\right)\Psi(t,\mathbf{x})\big{|}^{2}+V|\Psi(t,\mathbf{x})|^{2}}\\[5.69054pt] {\displaystyle\qquad\qquad+\frac{\epsilon}{2}|\mathbf{E}(t,\mathbf{x})|^{2}+\frac{1}{2\mu}|\mathbf{B}(t,\mathbf{x})|^{2}\big{)}\mathrm{d}\mathbf{x}.}\end{array}$$

$$\Psi ⟶{\Psi }^{\prime }=e^{\mathrm{i}\chi }\Psi ,\mathbf{A}⟶{\mathbf{A}}^{\prime }=\mathbf{A}+\nabla \chi ,\varphi ⟶{\varphi }^{\prime }=\varphi -\frac{\partial \chi }{\partial t},$$

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

$$\left\{\begin{array}[]{lll}{\displaystyle\Psi(\mathbf{x},t)=0,\quad\mathbf{A}(\mathbf{x},t)\times\mathbf{n}=0,\quad\phi(\mathbf{x},t)=0\quad\mathrm{on}\;\partial\Omega\times(0,T),}\\[5.69054pt] {\displaystyle\Psi(\mathbf{x},0)=\Psi_{0}(\mathbf{x}),\quad\mathbf{A}(\mathbf{x},0)=\mathbf{A}_{0}(\mathbf{x}),\quad\mathbf{A}_{t}(\mathbf{x},0)=\mathbf{A}_{1}(\mathbf{x})\quad\rm{in}\;\Omega,}\end{array}\right.$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\mathcal{E}(t)=\int_{\Omega}\left(\frac{1}{2}\big{|}\left[i\nabla+q\mathbf{A}\right]\Psi\big{|}^{2}+V|\Psi|^{2}+\frac{1}{2}\big{|}\frac{\partial\mathbf{A}}{\partial t}\big{|}^{2}+\frac{1}{2}|\nabla\times\mathbf{A}|^{2}+\frac{1}{2}|\nabla\phi|^{2}\right)\;\mathrm{d}\mathbf{x}.}\end{array}$$

$$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{H}_{t}^{1}(\Omega)=\{\mathbf{v}\in\mathbf{H}^{1}(\Omega)\,|\,\,\mathbf{v}\times\mathbf{n}=0\,\,\,{\rm on}\,\,\partial\Omega\},}\\[5.69054pt] {\displaystyle\mathbf{H}^{1}_{t,0}(\Omega)=\{\mathbf{v}\in\mathbf{H}_{t}^{1}(\Omega)\,|\,\,\nabla\cdot\mathbf{v}=0\},}\\[5.69054pt] {\displaystyle\mathbf{L}^{2}_{0}(\Omega)=\{\mathbf{v}\in\mathbf{L}^{2}(\Omega)\,|\,\,\nabla\cdot\mathbf{v}=0\;{\rm weakly}\,\}}.\end{array}$$

$$|\mathbf{u}{|}_{{\mathbf{H}}_t^1}:=\Vert \nabla \cdot \mathbf{u}{\Vert }_{L^2(\Omega )}+\Vert \nabla \times \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )},|\mathbf{u}{|}_{{\mathbf{H}}_{t,0}^1}:=\Vert \nabla \times \mathbf{u}{\Vert }_{{\mathbf{L}}^2(\Omega )},$$

$$\displaystyle\Psi\in C\big{(}[0,T];\mathcal{L}^{2}(\Omega)\big{)}\cap L^{\infty}\big{(}0,T;\mathcal{H}_{0}^{1}(\Omega)\big{)},\;\frac{\partial\Psi}{\partial t}\in L^{\infty}\big{(}0,T;\mathcal{H}^{-1}(\Omega)\big{)},$$

$$\displaystyle\phi\in C\big{(}[0,T];{L}^{2}(\Omega)\big{)}\cap L^{\infty}\big{(}0,T;H_{0}^{1}(\Omega)\big{)},\;\frac{\partial\phi}{\partial t}\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)},$$

$$\displaystyle\mathbf{A}\in C\big{(}[0,T];\,\mathbf{L}^{2}(\Omega)\big{)}\cap L^{\infty}\big{(}0,T;\,\mathbf{H}^{1}_{t,0}(\Omega)\big{)},$$

$$\displaystyle\frac{\partial\mathbf{A}}{\partial t}\in C\big{(}[0,T];\,\big{(}\mathbf{H}^{1}_{t,0}(\Omega)\big{)}^{\prime}\,\big{)}\cap L^{\infty}\big{(}0,T;\,\mathbf{L}^{2}(\Omega)\big{)},$$

$$\int_{0}^{T}\bigg{[}\mathrm{i}\big{(}\Psi,\,\frac{\partial\widetilde{\Psi}}{\partial t}\big{)}+\frac{1}{2}\big{(}\left(\mathrm{i}\nabla+\mathbf{A}\right)\Psi,\,\left(\mathrm{i}\nabla+\mathbf{A}\right)\widetilde{\Psi}\big{)}+\big{(}V\Psi,\widetilde{\Psi}\big{)}+\big{(}\phi\Psi,\widetilde{\Psi}\big{)}\bigg{]}\mathrm{d}t=0,$$

$$\int_{0}^{T}\bigg{[}\big{(}\mathbf{A},\frac{\partial^{2}\widetilde{\mathbf{A}}}{\partial t^{2}})+\big{(}\nabla\times\mathbf{A},\nabla\times\widetilde{\mathbf{A}}\big{)}+\big{(}\frac{\mathrm{i}}{2}(\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}),\widetilde{\mathbf{A}}\big{)}+\big{(}|\Psi|^{2}\mathbf{A},\widetilde{\mathbf{A}}\big{)}\bigg{]}\mathrm{d}t=0,$$

$$\displaystyle\int_{0}^{T}\left[\big{(}\nabla\phi,\,\nabla\widetilde{\phi}\big{)}-\big{(}|\Psi|^{2},\,\widetilde{\phi}\big{)}\right]\mathrm{d}t=0,$$

$$\displaystyle\mathbf{A}\in C\big{(}[0,T];\,\mathbf{L}^{2}(\Omega)\big{)}\cap L^{\infty}\big{(}0,T;\,\mathbf{H}_{t}^{1}(\Omega)\big{)},$$

$$\displaystyle\frac{\partial\mathbf{A}}{\partial t}\in C\big{(}[0,T];\,(\mathbf{H}_{t}^{1}(\Omega))^{\prime}\big{)}\cap L^{\infty}\big{(}0,T;\,\mathbf{L}^{2}(\Omega)\big{)},$$

$$\begin{array}[]{lll}{\displaystyle\int_{0}^{T}\bigg{[}\big{(}\mathbf{A},\frac{\partial^{2}\widetilde{\mathbf{A}}}{\partial t^{2}})+\big{(}\nabla\times\mathbf{A},\nabla\times\widetilde{\mathbf{A}}\big{)}+\big{(}\frac{\mathrm{i}}{2}(\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}),\,\widetilde{\mathbf{A}}\big{)}}\\[5.69054pt] {\displaystyle\qquad\quad-\big{(}\frac{\partial\phi}{\partial t},\,\nabla\cdot\widetilde{\mathbf{A}}\big{)}+\big{(}|\Psi|^{2}\mathbf{A},\,\widetilde{\mathbf{A}}\big{)}\bigg{]}\mathrm{d}t=0,\quad\forall\widetilde{\mathbf{A}}\in C_{0}^{2}\big{(}(0,T);\,\mathbf{H}_{t}^{1}(\Omega)\big{)}}.\end{array}$$

$$\begin{array}[]{lll}{\displaystyle\int_{0}^{T}\Big{[}\big{(}\frac{\partial\rho}{\partial t},\,\varphi\big{)}-\big{(}\mathbf{J},\,\nabla\varphi\big{)}\Big{]}\eta(t)\mathrm{d}t=0,}\\[5.69054pt] {\displaystyle\displaystyle\int_{0}^{T}\Big{[}\big{(}\frac{\partial\phi}{\partial t},\,\Delta\varphi\big{)}+\big{(}\frac{\partial\rho}{\partial t},\,\varphi\big{)}\Big{]}\eta(t)\mathrm{d}t=0,}\end{array}$$

$$\rho =|\Psi {|}^2,\mathbf{J}=-\frac{\mathrm{i}}{2}({\Psi }^{\ast }\nabla \Psi -\Psi \nabla {\Psi }^{\ast })-|\Psi {|}^2\mathbf{A}.$$

$$\big{(}\frac{\partial\phi}{\partial t},\,\Delta\varphi\big{)}+\big{(}\mathbf{J},\,\nabla\varphi\big{)}=0,\quad\forall\varphi\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega).$$

$$\mathbf{v}={\mathbf{v}}_0+\nabla \phi ,{\mathbf{v}}_0\in {\mathbf{H}}_{t,0}^1(\Omega ),\phi \in H_0^1(\Omega )\cap H^2(\Omega ),$$

$$\Vert {\mathbf{v}}_0{\Vert }_{{\mathbf{H}}^1}\le C\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}^1},\Vert \phi {\Vert }_{H^2}\le C\Vert \mathbf{v}{\Vert }_{{\mathbf{H}}^1}.$$

$$\big{(}\frac{\partial\mathbf{A}}{\partial t},\,\nabla\varphi)=0,\;\forall\varphi\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega).$$

$$\frac{\partial\mathbf{A}}{\partial t}\in C\big{(}[0,T];\,(\mathbf{H}_{t}^{1}(\Omega))^{\prime}\big{)}$$

$$\begin{array}[]{lll}{\displaystyle\int_{0}^{T}\bigg{[}\big{(}\mathbf{A},\frac{\partial^{2}\widetilde{\mathbf{A}}}{\partial t^{2}})+\big{(}\nabla\times\mathbf{A},\nabla\times\widetilde{\mathbf{A}}\big{)}+\big{(}\frac{\mathrm{i}}{2}(\Psi^{*}\nabla{\Psi}-\Psi\nabla{\Psi}^{*}),\,\widetilde{\mathbf{A}}\big{)}+\big{(}|\Psi|^{2}\mathbf{A},\,\widetilde{\mathbf{A}}\big{)}}\\[5.69054pt] {\displaystyle\qquad\quad-\big{(}\frac{\partial\phi}{\partial t},\,\nabla\cdot\widetilde{\mathbf{A}}\big{)}\bigg{]}\mathrm{d}t=0,\;\forall\widetilde{\mathbf{A}}\in C_{0}^{2}\big{(}(0,T);\,\nabla\big{(}H_{0}^{1}(\Omega)\cap H^{2}(\Omega)\big{)}\big{)}.}\end{array}$$

$$\int_{0}^{T}\big{(}\nabla\cdot\mathbf{A},\,\varphi)\frac{{\rm d}^{2}\eta}{{\rm d}t^{2}}\,\mathrm{d}t=0,$$

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

$$\left\{\begin{array}[]{@{}l@{}}{\displaystyle-\Delta u_{k}=\lambda_{k}u_{k}\quad in\;\Omega,}\\[5.69054pt] {\displaystyle u_{k}=0\quad on\;\partial\Omega,}\end{array}\right.$$

Full text

1 Introduction

In this paper, we consider one of the fundamental equations of nonrelativitic quantum mechanics, the Maxwell–Schrödinger (M-S) system, which describes the time-evolution of an electron within its self-consistent generated and external electromagnetic fields. In this system, the Schrödinger’s equation can be written as follows:

[TABLE]

where $\Omega_{T}=\Omega\times(0,T)$. Here $\Psi$, $m$, and $q$ denote the wave function, the mass, and the charge of the electron, respectively. $V$ is the time-independent potential energy and is assumed to be bounded in this paper. The magnetic potential $\mathbf{A}$ and the electric potential $\phi$ are obtained by solving the following equations:

[TABLE]

where the electric fields $\mathbf{E}$ and the magnetic fields $\mathbf{B}$ satisfy the Maxwell’s equations:

[TABLE]

Here $\epsilon$ and $\mu$ denote the electric permittivity and the magnetic permeability of the material, respectively. The charge density $\rho$ and the current density $\mathbf{J}$ are defined as follows

[TABLE]

where $\Psi^{\ast}$ denotes the conjugate of $\Psi$.

Substituting (2) into (3) and combining (1) and (4), we have the following M-S system

[TABLE]

We assume that $\Omega$ is a bounded domain in $\mathbb{R}^{d}$ ($d=2,\,3$). The total energy of the system, at time $t$, are defined as follows

[TABLE]

For a smooth solution $(\Psi,\mathbf{A},\phi)$ satisfying certain appropriate boundary conditions, the energy is a conserved quantity.

It is well known that the solutions of the above M-S system are not uniquely determined. In fact, the M-S system is invariant under the gauge transformation:

[TABLE]

for any sufficiently smooth function $\chi:\Omega\times(0,T)\rightarrow\mathbb{R}$. That is, if $(\Psi,\mathbf{A},\phi)$ satisfies M-S system, then so does $(\Psi^{\prime},\mathbf{A}^{\prime},\phi^{\prime})$.

In view of the gauge freedom, to obtain mathematically well-posed equations, we can impose some extra constraint, commonly known as gauge choice, on the solutions of the M-S system. In this paper, we study the M-S system in the Coulomb gauge, i.e. $\nabla\cdot\mathbf{A}=0$.

By employing the atomic units, i.e. $\hbar=m=q=1$, and assuming that $\epsilon=\mu=1$, the M-S system in the Coulomb gauge (M-S-C) can be reformulated as follow:

[TABLE]

In this paper, the M-S-C system (8) is considered in conjunction with the following initial boundary conditions:

[TABLE]

with $\nabla\cdot\mathbf{A}_{0}=\nabla\cdot\mathbf{A}_{1}=0$.

In the M-S-C system, the energy $\mathcal{E}(t)$ takes the following form

[TABLE]

Remark 1

The boundary condition $\Psi(\mathbf{x},t)=0$ implies that the particle is confined in the domain $\Omega$. The boundary conditions $\mathbf{A}(\mathbf{x},t)\times\mathbf{n}=0$ and $\phi(\mathbf{x},t)=0$ denote the perfect conductive boundary condition (PEC). We refer readers to Weng for the determination of the boundary conditions for the vector potential $\mathbf{A}$ and the scalar potential $\phi$ in different electromagnetic environment.

The existence and uniqueness of (smooth) solutions to the time-dependent M-S equations (5) in all of $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$ have been studied in Guo ; Nak ; Nak-1 ; Nak-2 . However, these results don’t hold for bounded domains because some important tools used in these work can’t apply to bounded domains. For example, Strichartz estimates and many tools from Fourier analysis. In this paper, we will prove the existence of weak solutions to the M-S-C system in a bounded smooth domain by Galerkin’s method and compactness arguments. To the best of our knowledge, this is the first result on the existence of weak solutions to the inital-boundary problem of the M-S-C system in a bounded smooth domain.

In recent years, with the development of nanotechnology, there has been considerable interest in developing physical models and numerical methods to simulate light-matter interaction at the nanoscale. Due to the natural coupling of the electromagnetic fields and quantum effects, the Maxwell–Schrödinger model is widely used in simulating self-induced transparency in quantum dot systems karni , laser-molecule interaction Lor-1 , carrier dynamics in nanodevices Pi and molecular nanopolaritonics Lop . However, in these existing methods, the Maxwell’s equations of field type (3), instead of the potential type in (5), are usually coupled to the Schrödinger’s equation through the dipole approximation or by extracting the vector potential $\mathbf{A}$ and the scalar potential $\phi$ from the electric field $\mathbf{E}$ and the magnetic field $\mathbf{H}$ Ah-1 ; Oh ; Sui ; Tur . In part because there exists robust numerical algorithms for the Maxwell’s equations (3), for example, the time domain finite difference (FDTD) method, the transmission line matrix (TLM) method, etc. Recently, RYU Ryu used the FDTD scheme to discretize the Maxwell–Schrödinger equations (5) directly in the Lorentz gauge to simulate a single electron in an artificial atom excited by an incoming electromagnetic field. But so far, there are rather limited studies on the numerical algorithms of the M-S system (5) as well as their convergence analysis.

In this paper we will present a fully discrete finite element method for solving the problem (8)-(9) and show that it is equivalent to a fully discrete Crank–Nicolson scheme based on mixed finite element method. We will show that our scheme maintains the conservation properties of the original system. Compared with the commonly used method which couples the Maxwell’s equations of field type with the Schrödinger’s equation and solves the system by the FDTD method, our method keeps the total charge and energy of the discrete system conserved and may suffer from less restriction in the time step-size since we use the Crank–Nicolson scheme in the time direction. In this paper we establish the optimal error estimates for the proposed method without any restrictions on the spatial mesh step $h$ and the time step $\tau$. In general it is very difficult to derive error estimates without any restrictions on the spatial mesh step and the time step for the highly complicated, nonlinear equations since the inverse inequalities are usually used to bound the nonlinear terms. In this paper we avoid using the inverse inequalities due to two aspects. On the one hand, we deduce some stability estimates of the approximate solutions by using the conservation properties of our scheme. More importantly, we take advantage of the special structures of the system and make some difficult nonlinear terms in the Schrödinger’s equation and the Maxwell’s equations respectively cancel out. To the best of our knowledge, this is the first theoretical analysis on the numerical algorithms for the M-S-C system (8).

The rest of this paper is organized as follows. In section 2 we introduce some notation and prove the existence of weak solutions to the M-S-C system (8)-(9). In section 3, we present two fully discrete finite element schemes for the M-S-C system and show that they are equivalent. Section 4 is devoted to the proof of energy-conserving property of the discrete system and some stability estimates of the approximate solutions. In section 5, we prove the existence and uniqueness of solutions to the discrete system. The optimal error estimates without any restrictions on the time step are derived in section 6. We provide some numerical experiments in section 7 to confirm our theoretical analysis.

2 Global existence of weak solutions to the M-S-C system

In this section, we study the existence of weak solutions to the M-S-C system (8) together with the initial-boundary conditions (9) in a bounded smooth domain. For simplicity, We introduce some notation below.

For any nonnegative integer $s$, we denote $W^{s,p}(\Omega)$ as the conventional Sobolev spaces of the real-valued functions defined in $\Omega$ and $W^{s,p}_{0}(\Omega)$ as the subspace of $W^{s,p}(\Omega)$ consisting of functions whose traces are zero on $\partial\Omega$. As usual, we denote $H^{s}(\Omega)\,=\,W^{s,2}(\Omega)$, $H^{s}_{0}(\Omega)\,=\,W^{s,2}_{0}(\Omega)$, and $L^{p}(\Omega)\,=\,W^{0,p}(\Omega)$, respectively. We use $\mathcal{H}^{s}(\Omega)=\{u+\mathrm{i}v\,|\,u,v\in H^{s}(\Omega)\}$ and $\mathcal{L}^{p}(\Omega)=\{u+\mathrm{i}v\,|\,u,v\in L^{p}(\Omega)\}$ with calligraphic letters for Sobolev spaces and Lebesgue spaces of the complex-valued functions, respectively. Furthermore, let $\mathbf{H}^{s}(\Omega)=[H^{s}(\Omega)]^{d}$ and $\mathbf{L}^{p}(\Omega)=[L^{p}(\Omega)]^{d}$ with bold faced letters be Sobolev spaces and Lebesgue spaces of the vector-valued functions with $d$ components ($d$=2, 3). The dual spaces of $\mathcal{H}_{0}^{s}(\Omega)$, ${H}_{0}^{s}(\Omega)$, and $\mathbf{H}_{0}^{s}(\Omega)$ are denoted by $\mathcal{H}^{-s}(\Omega)$, ${H}^{-s}(\Omega)$, and $\mathbf{H}^{-s}(\Omega)$, respectively. $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 consider the following subspaces of $\mathbf{H}^{1}(\Omega)$ and $\mathbf{L}^{2}(\Omega)$:

[TABLE]

The semi-norms on $\mathbf{H}_{t}^{1}(\Omega)$ and $\mathbf{H}^{1}_{t,0}(\Omega)$ are defined by

[TABLE]

both of which are equivalent to the standard $\mathbf{H}^{1}(\Omega)$-norm $\|\mathbf{u}\|_{\mathbf{H}^{1}(\Omega)}\,$Gir .

To take into account the time dependence, for any Banach space $W$ and integer $s\geq 1$, we define function spaces C\big{(}[0,T],\;W\big{)}, C\big{(}(0,T),\;W\big{)}, and C^{s}_{0}\big{(}(0,T),\;W\big{)} consisting of $W$-valued functions in $C[0,T]$, $C(0,T)$, and $C^{s}_{0}(0,T)$, respectively.

We now give two definitions of weak solution to the M-S-C system (8) together with the initial-boundary conditions (9).

Definition 1 (Weak solution I)

$(\Psi,\,\mathbf{A},\,\phi)$ is a weak solution of type I to (8)-(9), if

[TABLE]

with the initial condition $\Psi(\cdot,0)\in\mathcal{H}^{1}_{0}(\Omega)$, $\mathbf{A}(\cdot,0)\in\mathbf{H}^{1}_{t,0}(\Omega)$, $\mathbf{A}_{t}(\cdot,0)\in\mathbf{L}_{0}^{2}(\Omega)$, and the variational equations

[TABLE]

[TABLE]

[TABLE]

hold for all \widetilde{\Psi}\in C_{0}^{1}\big{(}(0,T);\,\mathcal{H}_{0}^{1}(\Omega)\big{)}, \widetilde{\mathbf{A}}\in C_{0}^{2}\big{(}(0,T);\,\mathbf{H}^{1}_{t,0}(\Omega)\big{)} and \widetilde{\phi}\in C\big{(}(0,T);\\ \,H_{0}^{1}(\Omega)\big{)}.

Definition 2 (Weak solution II)

$(\Psi,\,\mathbf{A},\,\phi)$ is a weak solution of type II to (8)-(9), if (11a), (11b), (12) and (14) in Definition 1 are satisfied and

[TABLE]

[TABLE]

The following theorem shows that the above two definitons of weak solutions are equivalent.

Theorem 2.1

The weak solutions to the M-S-C system defined in Definiton 1 and Definition 2 are equivalent.

Proof

It suffices to show that the vector potential $\mathbf{A}$ in Definiton 1 and Definiton 2 are consistent.

For any $\varphi(\mathbf{x})\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega),\;\eta(t)\in C_{0}^{1}(0,\,T)$, by choosing $\widetilde{\Psi}=\eta\Psi\varphi$ in (12), $\widetilde{\phi}=\frac{{\rm d}\eta}{{\rm d}t}\varphi$ in (14), and taking the imaginary part of (12), we have

[TABLE]

where

[TABLE]

Since $\eta(t)$ is arbitrary, from (17) we see that

[TABLE]

For any $\mathbf{v}\in\mathbf{H}_{t}^{1}(\Omega)$, we have the Helmholtz decomposition Gir :

[TABLE]

We first prove that the vector potential $\mathbf{A}$ given by Definition 1 satisfies (15a), (15b), and (16) in Definition 2. Obviously $\mathbf{A}$ satisfies (15a). Thanks to (11c) and (11d), we deduce that

[TABLE]

Then

[TABLE]

follows from (11d), (20)-(22) and thus $\mathbf{A}$ satisfies (15b). Since $\nabla\cdot\mathbf{A}=0$, by using (19), we have

[TABLE]

By applying (13), (24), and the Helmholtz docomposition (20), we find that $\mathbf{A}$ satisfies (16).

Next we assume that $\mathbf{A}$ is the vector potential given by Definition 2. It is easy to see that $\mathbf{A}$ satisfies (11d) and (14). For any $\eta(t)\in C_{0}^{2}(0,T),\;\varphi\in H_{0}^{1}(\Omega)\cap H^{2}(\Omega)$, take $\widetilde{\mathbf{A}}=\eta\nabla\varphi$ in (16) and employ (19) to find

[TABLE]

which implies that

[TABLE]

Consequently $\mathbf{A}$ satisfies (11c) and we complete the proof of this theorem.

Remark 2

Theorem 2.1 shows that weak solutions II satisfy $\nabla\cdot\mathbf{A}=0$ implicity.

Next we use the Galerkin method and compactness arguments to prove the existence of weak solutions to (8)-(9).

We first introduce two lemmas to construct finite dimensional subspaces of $\mathcal{H}_{0}^{1}(\Omega)$, $\mathbf{H}^{1}_{t,0}(\Omega)$, and $H^{1}_{0}(\Omega)$.

Lemma 1

Suppose that $\Omega\subset\mathbb{R}^{d}(d=2,3)$ is a bounded smooth domain. Then there exists a sequence $\{u_{k}\}_{k=1}^{\infty}$ being an orthogonal basis of $H_{0}^{1}(\Omega)$ as well as an orthonormal basis of $L^{2}(\Omega)$. Here $u_{k}\in H_{0}^{1}(\Omega)$ is an eigenfunction corresponding to $\lambda_{k}$:

[TABLE]

for $k=1,2\cdots.$

The proof of Lemma 1 is given in evans . It is worth pointing out that the conclusion is also true for complex-valued functions, i.e. there exists a sequence $\{\varphi_{k}\}_{k=1}^{\infty}$ being an orthogonal basis of $\mathcal{H}_{0}^{1}(\Omega)$ as well as an orthonormal basis of $\mathcal{L}^{2}(\Omega)$.

Lemma 2

Suppose that $\Omega\subset\mathbb{R}^{d}(d=2,3)$ is a bounded smooth domain. Then there exists an orthonormal basis $\{\mathbf{v}_{k}\}_{k=1}^{\infty}$ of $\mathbf{L}^{2}(\Omega)$, where $\mathbf{v}_{k}\in\mathbf{H}^{1}_{t,0}(\Omega)$ is an eigenfunction corresponding to $\mu_{k}$:

[TABLE]

for $k=1,2\cdots.$ Furthermore, $\{\mathbf{v}_{k}\}_{k=1}^{\infty}$ is an orthogonal basis of $\mathbf{H}^{1}_{t,0}(\Omega)$.

Proof

Let L:\mathbf{H}^{1}_{t,0}(\Omega)\rightarrow\big{(}\mathbf{H}^{1}_{t,0}(\Omega)\big{)}^{\prime} be defined by

[TABLE]

By the Lax-Milgram theorem, L^{-1}:\big{(}\mathbf{H}^{1}_{t,0}(\Omega)\big{)}^{\prime}\rightarrow\mathbf{H}^{1}_{t,0}(\Omega) exists and is bounded. Since $\mathbf{H}^{1}_{t,0}(\Omega)$ is compactly embedded into $\mathbf{L}^{2}(\Omega)$, $S=L^{-1}$ is a bounded, linear, compact operator mapping $\mathbf{L}^{2}(\Omega)$ into itself. It is easy to show that $S$ is self-adjoint. Then by the Hilbert-Schmidt theorm, there exists a countable orthonormal basis $\{\mathbf{v}_{k}\}_{k=1}^{\infty}$ of $\mathbf{L}^{2}(\Omega)$ consisting of eigenfunctions of $S$. The proof of $\{\mathbf{v}_{k}\}_{k=1}^{\infty}$ being an orthogonal basis of $\mathbf{H}^{1}_{t,0}(\Omega)$ is straightfoward.

Let $\mathcal{Z}_{n}$, $\mathbf{Z}_{n}$, and $Z_{n}$ be n-dimensional subspaces of $\mathcal{H}_{0}^{1}(\Omega)$, $\mathbf{H}^{1}_{t,0}(\Omega)$, and $H^{1}_{0}(\Omega)$, respectively,

[TABLE]

where $\{\varphi_{k}\}_{k=1}^{\infty}$, $\{\mathbf{v}_{k}\}_{k=1}^{\infty}$, and $\{u_{k}\}_{k=1}^{\infty}$ are given in Lemma 1 and 2.

For each $n$, we can construct the Galerkin approximate solutions of weak solutions to the M-S-C system in the sense of Definition 1 as follows.

Find $\Psi_{n}(t)=\sum_{j=1}^{n}{h_{jn}(t)\varphi_{j}\in\mathcal{Z}_{n}}$, $\mathbf{A}_{n}(t)=\sum_{j=1}^{n}{g_{jn}(t)\mathbf{v}_{j}\in\mathbf{Z}_{n}}$, and $\phi_{n}(t)=\sum_{j=1}^{n}{f_{jn}(t)u_{j}\in Z_{n}}$ such that

[TABLE]

for any $t\in(0,T)$, $j=1,2,\cdots,n$. Here $\mathcal{P}_{n}$ and $\mathbf{P}_{n}$ denote the orthogonal projection onto $\mathcal{Z}_{n}$ and $\mathbf{Z}_{n}$, respectively.

Using the local existence and uniqueness theory on ODEs, we can show that the nonlinear differential system (27) has a unique local solution defined on some interval $[0,T_{n}]$. Next we derive some $a\,\,priori$ estimates to extend the local solution to a global solution defined on $[0,T]$. In this paper, the following lemma will be used frequently.

Lemma 3

Let $2<p<6$. Suppose that $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$ is a bounded Lipschitz domain. Then for each $\epsilon>0$, there exists some constant $C_{\epsilon}$ depending on $\epsilon$ (and on $p$ and $\Omega$) such that

[TABLE]

Lemma 1 can be proved by applying Sobolev’s embedding thorems, Poincar$\acute{\rm e}$’s inequality, and the following lemma in temma .

Lemma 4

Let $W_{0}$, $W$, and $W_{1}$ be three Banach spaces such that $W_{0}\subset W\subset W_{1}$, the injection of $W$ into $W_{1}$ being continuous, and the injection of $W_{0}$ into $W$ is compact. Then for each $\epsilon>0$, there exists some constant $C_{\epsilon}$ depending on $\epsilon$ (and on the spaces $W_{0}$, $W$, $W_{1}$) such that

[TABLE]

We define the energy $\mathcal{E}_{n}(t)$ of the approximate system (27) as follows.

[TABLE]

Lemma 5

For any $t\in(0,T]$, if the solution $(\Psi_{n}(t),\mathbf{A}_{n}(t),\phi_{n}(t))$ of the approximate system (27) exists, we have the conservation of total charge and energy

[TABLE]

Proof

$(\ref{eq:6-3-0})_{1}$ can be proved by multiplying the first equation of (27) by $h_{jn}$, summing $j=1,2,\cdots,n$, and taking the imaginary part.

To prove $(\ref{eq:6-3-0})_{2}$, we first multiply the first equation of (27) by $\frac{d}{dt}h_{jn}$, sum $j=1,2,\cdots,n$, and take the real part. We discover

[TABLE]

Multiplying the second equation of (27) by $\frac{d}{dt}g_{jn}$ , summation gives

[TABLE]

Differentiating both sides of the third equation of (27) with respect to $t$, multiplying it by $f_{jn}$ and summing $j=1,2,\cdots,n$, we obtain

[TABLE]

Adding (30), (31), and (32) together completes the proof of Lemma 2.

We now establish some estimates of the solution $(\Psi_{n}(t),\mathbf{A}_{n}(t),\phi_{n}(t))$.

Theorem 2.2

For any $t\in(0,T]$, if the solution $(\Psi_{n}(t),\mathbf{A}_{n}(t),\phi_{n}(t))$ exists, then it satisfies the estimates

[TABLE]

[TABLE]

where $C$ is independent of $t$ and $n$.

Proof

By the definition of initial data $(\Psi_{n}(0),\mathbf{A}_{n}(0),\phi_{n}(0))$, it is easy to show that $\mathcal{E}_{n}(0)\leq C$. Thus by applying (29), we have

[TABLE]

Since the semi-norm in $\mathbf{H}^{1}_{t,0}(\Omega)$ is equivalent to $\mathbf{H}^{1}$-norm, 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 Lemma 1, we further prove

[TABLE]

From (35), (38) and Lemma 2, we deduce

[TABLE]

Consequently, we obtain

[TABLE]

By applying Poincar$\acute{\rm e}$’s inequality and (35), we see that

[TABLE]

Therefore, (33) is proved by combining (36), (39), and (40).

To estimate $\|\frac{\partial\Psi_{n}}{\partial t}\|_{\mathcal{H}^{-1}}$, we first fix $\omega\in\mathcal{H}_{0}^{1}(\Omega)$ with $\|\omega\|_{\mathcal{H}_{0}^{1}}\leq 1$. Note that $\{\varphi_{n}\}_{n=1}^{\infty}$ is an orthogonal basis of $\mathcal{H}_{0}^{1}(\Omega)$ as well as an orthonormal basis of $\mathcal{L}^{2}(\Omega)$. Thus we can write $\omega={\omega}^{1}+{\omega}^{2}$, where ${\omega}^{1}\in\mathcal{Z}_{n}={\rm span}\{\varphi_{k}\}_{k=1}^{n}$ and $({\omega}^{2},\varphi_{k})=0$ for $k=1,\cdots,n.$ It is clear that $\|{\omega}^{1}\|_{\mathcal{H}_{0}^{1}}\leq 1$. Then the first equation of (27) implies that

[TABLE]

Thus by applying (33), we obtain

[TABLE]

which implies that

[TABLE]

Similarly, we can prove

[TABLE]

It remains to show $\|\frac{\partial\phi_{n}}{\partial t}\|_{L^{2}}\leq C.$ In order to estimate $\|\frac{\partial\phi_{n}}{\partial t}\|_{L^{2}}$, we fix $f\in L^{2}(\Omega)$ and find $\psi_{n}\in Z_{n}={\rm span}\{u_{k}\}_{k=1}^{n}$ such that

[TABLE]

Then by differentiating both sides of the third equation of (27) with respect to $t$ and using (44), we have

[TABLE]

Thus from (42) and (39), we deduce

[TABLE]

Next we will prove $\|\psi_{n}\|_{H^{2}}\leq C\|f\|_{L^{2}}$. Note that each basis function $u_{j}\in Z_{n}$ satisfies $-\Delta u_{j}=\lambda_{j}u_{j}$. By (44), we have

[TABLE]

It follows that

[TABLE]

which implies that

[TABLE]

It is easy to show $\|\nabla\psi_{n}\|_{{L}^{2}}\leq C\|f\|_{L^{2}}$. Hence we get

[TABLE]

Combining (46) and (50), we find $\|\frac{\partial\phi_{n}}{\partial t}\|_{L^{2}}\leq C.$ Thus we complete the proof of Theorem 2.2.

Using the above energy estimates, we have

Corollary 1

Given $T>0$, the nonlinear differential system (27) has a unique global solution $(\Psi_{n},\mathbf{A}_{n},\phi_{n})$, which satisfies (33) and (34).

We now quote a compactness lemma which can be found in Simon .

Lemma 6

Let $B_{0}$, $B$, $B_{1}$ be three Banach spaces such that $B_{0}\subset B\subset B_{1}$ with continuous embedding and the embedding $B_{0}\subset B$ is compact. Suppose $F$ is a bounded set in $L^{\infty}(a,b;B_{0})$ such that $\frac{\partial F}{\partial t}=\{\frac{\partial f}{\partial t}\,|\,f\in F\}$ is bounded in $L^{p}(a,b;B_{1})$ for some $p>1$. Then $F$ is relatively compact in $C([0,T];B)$.

From Lemma 6 and Theorem 2.2, we deduce that there exists $\Psi\in C([0,T];\\ \mathcal{L}^{p}(\Omega))\cap L^{\infty}((0,T);\mathcal{H}_{0}^{1}(\Omega))$, $\mathbf{A}\in C([0,T];\mathbf{L}^{p}(\Omega))\cap L^{\infty}((0,T);\mathbf{H}^{1}_{t,0}(\Omega))$, $\phi\in C([0,T];{L}^{p}(\Omega))\cap L^{\infty}((0,T);H_{0}^{1}(\Omega))$ and a subsequence $\{(\Psi_{n_{k}},\mathbf{A}_{n_{k}},\phi_{n_{k}})\}$ such that as $n_{k}\rightarrow\infty$

[TABLE]

Here $1<p<\infty$ for $d=2$ and $1<p<6$ for $d=3$.

Furthermore, we have the following convergence properties for time derivatives of $\{(\Psi_{n_{k}},\mathbf{A}_{n_{k}},\phi_{n_{k}})\}$.

[TABLE]

Passing to the limits $n_{k}\rightarrow\infty$ in our Galerkin appriximations, we obtain

Theorem 2.3

Given $T>0$, there exists a weak solution $(\Psi,\mathbf{A},\phi)$ to the M-S-C system (8)-(9) in the sense of Definition 1, which satisfies the conservation of the total energy:

[TABLE]

where $\mathcal{E}(t)$ is given in (10).

Here we omit the proof of the weak limit satisfies (12)-(14) since the technique is standard.

Remark 3

By making a slight modification of the proof in this section, Theorem 2.3 holds for the M-S-C system with bounded coefficients:

[TABLE]

where $0<\alpha_{0}\leq m(\mathbf{x})\leq\alpha_{1}$, $0<\beta_{0}\leq\mu(\mathbf{x})\leq\beta_{1}$.

In particular, Theorem 2.3 holds for the M-S-C system with rapidly oscillating discontinuous coefficients arising from the modeling of a heterogeneous structure with a periodic microstructure, i.e. $m(\mathbf{x})=a(\frac{\mathbf{x}}{\varepsilon})$, $\mu(\mathbf{x})=b(\frac{\mathbf{x}}{\varepsilon})$, where $a({\bm{\xi}})$, $b({\bm{\xi}})$ are 1-periodic in ${\bm{\xi}}$ and $\varepsilon>0$ is a small parameter. Furthermore, if $\mathcal{E}(0)$, the initial energy, is independent of $\varepsilon$, then $\mathcal{E}(t)$ is also independent of $\varepsilon$.

3 Fully discrete finite element scheme

In this section, we consider the fully discretization of the M-S-C system (8)-(9) by the Galerkin finite element method in space together with the Crank-Nicolson scheme in time. In the following of the paper, we assume that $\Omega$ is a bounded Lipschitz polyhedron convex domain in $\mathbb{R}^{3}$.

Let us first triangulate the space domain $\Omega$ and assume that $\mathcal{T}_{h}=\{K\}$ is a regular partition of $\Omega$ into tetrahedrons of maximal diameter $h$. Without loss of generality, we assume that $0<h<1$. We denote by $P_{r}(K)$ the space of polynomials of degree $r$ defined on the element $K$. In the rest of this paper, we assume that $r\leq 2$ unless otherwise specified. For a given partition $\mathcal{T}_{h}$, we define the classical 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]

Let $\mathcal{I}_{h}$, $I_{h}$, and $\mathbf{I}_{h}$ be the commonly used Lagrange interpolation on $\mathcal{X}_{h}^{r}$, $X_{h}^{r}$, and $\mathbf{X}_{h}$, respectively. For $r=1,2$, $1\leq s\leq r+1$, we have the following interpolation error estimates Bre :

[TABLE]

We approximate the scalar potential $\phi$ and the wave function $\Psi$ in $X_{h}^{r}$ and $\mathcal{X}_{h}^{r}$ respectively, and find the approximate solution of the vector potential $\mathbf{A}$ in a subspace of $\mathbf{X}_{h}$:

[TABLE]

It is important to note that $\mathbf{X}_{0h}\nsubseteq\mathbf{H}^{1}_{t,0}(\Omega)$ since for each $\mathbf{v}_{h}\in\mathbf{X}_{0h}$, we only have $\rho_{h}(\nabla\cdot\mathbf{v}_{h})=0$, where $\rho_{h}$ is the orthogonal projection of $L^{2}(\Omega)$ onto $X_{h}^{1}.$

We now claim that there exists an interpolation operator ${\bm{i}}_{h}:\mathbf{H}^{1}_{t,0}(\Omega)\rightarrow\mathbf{X}_{0h}$, such that for every $\mathbf{v}\in\mathbf{H}^{1}_{t,0}(\Omega)\cap\mathbf{H}^{r+1}(\Omega)$,

[TABLE]

By the mixed finite element theory Brezzi ; Gir , we can ensure (59) by applying (57c) and the following discrete inf-sup condition: there exists a positive constant $\beta$, independent of $h$, such that

[TABLE]

For \widetilde{\mathbf{X}}_{h}=\big{(}Y^{2}_{h}\big{)}^{3}\cap\mathbf{H}^{1}_{0}(\Omega), $\widetilde{X}_{h}=Y^{1}_{h}$, the following discrete inf-sup condition for Hood-Taylor element is proved in Brezzi by Verfürth’s trick:

[TABLE]

The technique used in the proof of (61) can be applied directly to prove (60) by virtue of the fact that $\widetilde{\mathbf{X}}_{h}\subset\mathbf{X}_{h}$, $X_{h}^{1}\subset\widetilde{X}_{h}$ and the following continuous inf-sup condition:

[TABLE]

For more details, see Brezzi . Thus (59) is verified.

Let ${\bm{\pi}}_{h}:\mathbf{H}_{\rm t}^{1}(\Omega)\rightarrow\mathbf{X}_{0h}$ be a Ritz projection as follows: $\forall\,\mathbf{A}\in\mathbf{H}_{\rm t}^{1}(\Omega)$, find ${\bm{\pi}}_{h}\mathbf{A}\in\mathbf{X}_{0h}$ such that

[TABLE]

Owing to (59), we have the following error estimate of ${\bm{\pi}}_{h}$.

[TABLE]

To define our fully discrete scheme, we divide the time interval $(0,T)$ into $M$ uniform subintervals using the nodal points

[TABLE]

with $t^{k}=k\tau$ and $\tau=T/M$. We denote $u^{k}=u(\cdot,t^{k})$ for any given functions u\in C\big{(}(0,T);\,W\big{)} with a Banach space $W$. For a given sequence $\{u^{k}\}_{k=0}^{M}$, we introduce the following notation:

[TABLE]

For convenience, Let us assume that $\mathbf{A}^{-1}$ is defined by

[TABLE]

which is an approximation of $\mathbf{A}(\cdot,-\tau)$ with second order accuracy.

Using the above notation, we can formulate our first fully discrete finite element scheme for the M-S-C system as follows.

Scheme (I).    For $k=0,1,\cdots,M$, find $(\Psi_{h}^{k},\mathbf{A}_{h}^{k},\phi^{k}_{h})\in\mathcal{X}_{h}^{r}\times\mathbf{X}_{0h}\times X_{h}^{r}$ such that

[TABLE]

and for any $\varphi\in\mathcal{X}_{h}^{r}$, $\mathbf{v}\in\mathbf{X}_{0h}$, $u\in X_{h}^{r}$, the following equations hold:

[TABLE]

Remark 4

In Section 2 we prove that weak solutions in Definition 1 and Definition 2 are equivalent. However, from the standpoint of finite element numerical computation, we choose to approximate weak solutions of the M-S-C system in the sense of Definition 2 instead of Definition 1 since it is very difficult to construct finite element subspaces of $\mathbf{H}^{1}_{t,0}(\Omega)$. We know that weak solutions of type II in Definition 2 imply that $\mathbf{A}$ is divergence-free although we only require $\mathbf{A}$ in $\mathbf{H}^{1}_{\rm t}(\Omega)$. In the discrete level, if we approximate $\mathbf{A}$ in $\mathbf{X}_{h}$, it is difficult to degisn time integration schemes to ensure a discrete analogue of $\nabla\cdot\mathbf{A}=0$, i.e. $P_{h}(\nabla\cdot\mathbf{A}_{h})=0$, where $P_{h}$ denotes the orthogonal projection of $L^{2}(\Omega)$ onto some finite element space. Thus we approximate $\mathbf{A}$ in $\mathbf{X}_{0h}$ to enforce the projection of $\nabla\cdot\mathbf{A}_{h}$ onto $X_{h}^{1}$ vanishes. Moreover, for the purpose of theoretical analysis, we add an extra term \big{(}\nabla\cdot\widetilde{\mathbf{A}}_{h}^{k},\,\nabla\cdot\mathbf{v}\big{)} to the discrete system (68). It turns out that this term is indispensable to the proof of the error estimates and the existence of solutions to the discrete system (68).

Apart from introducing the subspace $\mathbf{X}_{0h}$ of $\mathbf{X}_{h}$, we can also introduce a Lagrangian multiplier $p_{h}^{k}$ to relax the divergence-free constraint of $\mathbf{A}_{h}^{k}$ at each time step. We now give another fully discrete scheme based on the mixed finite element method as follows.

[TABLE]

and for $k=1,2,\cdots,M$, find $(\Psi_{h}^{k},\,\mathbf{A}_{h}^{k},\,p_{h}^{k},\,\phi^{k}_{h})\in\mathcal{X}_{h}^{r}\times\mathbf{X}_{h}\times X_{h}^{1}\times X_{h}^{r}$ satisfying

[TABLE]

At each time step, the equation

[TABLE]

in scheme (I) and

[TABLE]

in scheme (II) are decoupled from the other two equations, respectively. Due to the discrete inf-sup condition (60) and the coercivity of the bilinear functional $a^{k}_{h}(\cdot,\,\cdot)$ in $\mathbf{X}_{h}$, where

[TABLE]

we know that there exists a unique solution to (71) and (72), respectively. It is easy to see that $\mathbf{A}_{h}^{k}$ in (72) satisfies (71) and thus the two above equations admit the same solution $\mathbf{A}_{h}^{k}$. Consequently, scheme (I) and scheme (II) are mathematically equivalent. However, scheme (I) is easier to perform theoretical analysis while scheme (II) is easier to carry out numerical computation.

At each time step, we first solve (71) or (72) and obtain $\mathbf{A}_{h}^{k}$. Then we substitute it into (70) and solve the nonlinear subsystem concerning $\Psi_{h}^{k}$ and $\phi_{h}^{k}$. The existence and uniqueness of solutions to this subsystem is proved in Section 5. In practical computations, we can apply the Picard simple iterative method or the Newton iterative method to solve the nonlinear subsystem.

For convenience, we define the following bilinear forms:

[TABLE]

Then (68) in scheme (I) can be rewritten as follows: for $k=1,2,\cdots,M$,

[TABLE]

In this paper we assume that the M-S-C system (8)-(9) has one and only one weak solution $(\Psi,\mathbf{A},\phi)$ in the sense of Definition 2 and the following regularity conditions are satisfied:

[TABLE]

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

[TABLE]

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

Theorem 3.1

Suppose that $\Omega\subset\mathbb{R}^{3}$ is a bounded Lipschitz polyhedral convex domain. Let $(\Psi,\mathbf{A},\phi)$ be the unique solution to the M-S-C system (8)-(9) , and let $(\Psi_{h}^{k},\mathbf{A}_{h}^{k},\phi_{h}^{k})$ be the numerical solution to the discrete system (67)-(68). Under the assumptions (76) and (77), 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})$, $r\leq 2$, and $C$ is a constant independent of $h$ and $\tau$.

4 Stability estimates

In this section we first show the discrete system (67)-(68) maintains the conservation of the total charge and energy. Then we deduce some stability estimates of the discrete solutions, which will be used to derive the error estimates in next section.

First we define the energy of the discrete system (67)-(68) as follows:

[TABLE]

Lemma 7 and Theorem 4.1 in the following are the discrete analogues of Lemma 5 and Theorem 2.2, respectively.

Lemma 7

For $k=1,2\cdots,M$, the solution $(\Psi^{k}_{h},\mathbf{A}_{h}^{k},\phi_{h}^{k})$ of the discrete system (67)-(68) satisfies

[TABLE]

Proof

The proof the this lemma is very similar to its continuous counterpart. For $(\ref{eq:3.4})_{1}$ we can simply choose $\varphi=\overline{\Psi}_{h}^{k}$ in $(\ref{eq:2-8})_{1}$ and take its imaginary part.

To prove $(\ref{eq:3.4})_{2}$, we first notice that

[TABLE]

We also have the following identities by direct calculations

[TABLE]

from which we deduce

[TABLE]

Thus we get

[TABLE]

Also we have

[TABLE]

By choosing $\varphi=\partial_{\tau}\Psi_{h}^{k}$ in $(\ref{eq:2-8})_{1}$, taking the real part of the equation and combining with (84) and (85), we get

[TABLE]

Next by taking $\mathbf{v}=\frac{1}{2\tau}(\mathbf{A}_{h}^{k}-\mathbf{A}_{h}^{k-2})$ and adding it to (86), we obtain

[TABLE]

Finally it is easy to deduce the following equation from the last equation of $(\ref{eq:2-8})$:

[TABLE]

Take $u=\overline{\phi}_{h}^{k}$ in (88), insert it into (87) and we complete the proof of $(\ref{eq:3.4})_{2}$.

Theorem 4.1

The solution of the discrete system (74) fulfills the following estimate

[TABLE]

where $C$ is independent of $h$, $\tau$ and $k$.

The proof of this theorem is very similar to its continuous counterpart, i.e. Theorem 2.2, and thus we omit the proof.

5 Solvability of the discrete system

In this section we consider the existence and uniqueness of the solutions to the discrete system (67)-(68). To prove it, we first introduce a useful lemma in Browder as follows.

Lemma 8

Let \big{(}H,\langle\cdot,\cdot\rangle\big{)} be a finite-dimensional inner product space, $\|\cdot\|_{H}$ be the associated norm, and $g:H\longrightarrow H$ be continuous. Assume that

[TABLE]

Then there exists a $z_{0}\in H$ such that $g(z_{0})=0$ and $\|z_{0}\|_{H}\leq\alpha$.

Theorem 5.1

For any \big{(}\Psi_{h}^{k-1},\mathbf{A}_{h}^{k-2},\mathbf{A}_{h}^{k-1},\phi_{h}^{k-1}\big{)} satisfies (89), there exists a solution \big{(}\Psi_{h}^{k},\mathbf{A}_{h}^{k},\phi_{h}^{k}\big{)} to the discrete system (68). Furthermore, if the time step $\tau$ is sufficiently small, the solution is unique.

Proof

As noted in Section 3, to solve the discrete system (67)-(68), we need to solve (71) and the following subsystem alternately.

[TABLE]

Since we have proved the solvability of (71) in Section 3, we only need to consider the solvability of (90), which can be rewritten as follows:

[TABLE]

For a given $h$, assume that $\{u_{1},u_{2},\cdots,u_{N}\}$ is a basis of $X_{h}^{r}$. Note that $\mathcal{X}_{h}^{r}=X_{h}^{r}\oplus{\rm i}X_{h}^{r}$. Then $\overline{\phi}_{h}^{k}$, $\overline{\Psi}_{h}^{k}$, and $\Psi_{h}^{k-1}$ can be written as follows

[TABLE]

where

[TABLE]

Denote by

[TABLE]

where

[TABLE]

Using the above notation, we can write (91) in the form of matrix:

[TABLE]

or a more compact form

[TABLE]

Now we define a finite dimensional space \big{(}H,\langle\cdot,\cdot\rangle\big{)} and a mapping $g:H\longrightarrow H$ as following:

[TABLE]

where $(\vec{e}_{2})^{\ast}$ denotes the conjugate transpose of $\vec{e}_{2}$ and $W$, $S$ and $Q$ are defined in (94)-(95). Obviously $g$ is continuous. Moreover,

[TABLE]

which implies that

[TABLE]

Thus the existence of $\vec{b}$ and $\vec{a}=W^{-1}\vec{r}(\vec{b})$ for (97) follows from Lemma 8. Combining with the existence of $\mathbf{A}_{h}^{k}$, we obtain the existence of \big{(}\Psi_{h}^{k},\mathbf{A}_{h}^{k},\phi_{h}^{k}\big{)} to the discrete system (68).

Now we study the uniqueness of the solutions to (90). Let $(\Psi_{h}^{k},\phi_{h}^{k})$ and $(\hat{\Psi}_{h}^{k},\hat{\phi}_{h}^{k})$ be two solutions of (90). Set $\eta=\Psi_{h}^{k}-\hat{\Psi}_{h}^{k},\;\psi=\phi_{h}^{k}-\hat{\phi}_{h}^{k}.$ They satisfy

[TABLE]

By choosing $\varphi=\eta$ in the first equation of (101) and taking its imaginary part, we obtain

[TABLE]

Take $u=\psi$ in the second equation of (101) and we get

[TABLE]

By substituting (103) into (102) and taking $\tau$ sufficiently small, we find $\eta=0$ and consequently obtain the uniqueness of the solutions.

6 The error estimates

We now turn to the proof of Theorem 3.1. Set $e_{\Psi}=\mathcal{I}_{h}\Psi-\Psi$, $e_{\mathbf{A}}={\bm{\pi}}_{h}\mathbf{A}-\mathbf{A}$ and $e_{\phi}=I_{h}\phi-\phi$, where $\mathcal{I}_{h}$, ${\bm{\pi}}_{h}$ and $I_{h}$ are defined in Section 3. From the regularity assumption and the interpolation error estimates (57a), (57b) and (64), we deduce

[TABLE]

By using standard finite element theory Bre and the regularity assumption (76), we also have

[TABLE]

In the rest of this paper, we need the following discrete integration by parts formulas.

[TABLE]

To simplify the notation , we denote by $\theta_{\Psi}^{k}=\Psi_{h}^{k}-\mathcal{I}_{h}\Psi^{k}$, $\theta_{\mathbf{A}}^{k}=\mathbf{A}_{h}^{k}-{\bm{\pi}}_{h}\mathbf{A}^{k}$ and $\theta_{\phi}^{k}=\phi_{h}^{k}-I_{h}\phi^{k}$. In view of the interpolation error estimates (104), we only need to prove that for $k=1,2,\cdots,M$, there holds

[TABLE]

By assuming that

[TABLE]

the weak solution $(\Psi,\mathbf{A},\phi)$ of the M-S-C system in sense of Definition 2 satisfies

[TABLE]

Subtracting the discrete sysytem (75) from (109) and noting that $\nabla\cdot\mathbf{A}^{k-1}=0$, we have

[TABLE]

[TABLE]

[TABLE]

In the following of the section, we analyze the above three error equations term by term, respectively.

6.1 Estimates for (110)

First, by taking $\varphi=\overline{\theta}_{\Psi}^{k}$ in (110), the imaginary part of the equation implies

[TABLE]

Now we are going to estimate the terms $V_{i}^{k}(\overline{\theta}_{\Psi}^{k}),\,i=1,\cdots,5$ one by one. By the error estimates for the interpolation operator $\mathcal{I}_{h}$ and the regularity of $\Psi$ in (76), we see that

[TABLE]

Thanks to

[TABLE]

and

[TABLE]

we get

[TABLE]

In order to estimate $V_{4}^{k}(\overline{\theta}_{\Psi}^{k})$, we first decompose it as follows.

[TABLE]

By using Theorem 4.1, the regularity assumption, and the properties of the interpolation operators, we obtain from (118) that

[TABLE]

Notice that

[TABLE]

By applying (82) and Theorem 4.1, it is easy to see that

[TABLE]

Now multiplying (113) by $\tau$, summing over $k=1,2,\cdots,m$, and applying the above estimates , we have

[TABLE]

Next, we take $\varphi=\partial_{\tau}{\theta_{\Psi}^{k}}$ in (110), which gives

[TABLE]

From the real part of (123) and (82), we obtain

[TABLE]

which yields

[TABLE]

From Theorem 4.1, we deduce

[TABLE]

Denoting by

[TABLE]

we have

[TABLE]

Now let us estimate \sum_{k=1}^{m}\mathrm{Re}\big{[}V_{j}^{k}(\partial_{\tau}{\theta_{\Psi}^{k}})\big{]}, $j=1,\cdots,5$ term by term. In light of (106), we get

[TABLE]

By employing the regularity assumption and the error estimates of interpolation operators, we see that

[TABLE]

The second term can be rewritten as

[TABLE]

Arguing as before, we obtain

[TABLE]

By the definition of the bilinear functional $B(\mathbf{A};\psi,\varphi)$ in (74), we can rewrite $\tau V_{3}^{k}(\partial_{\tau}{\theta_{\Psi}^{k}})$ as follows:

[TABLE]

By employing (106), (57a), the regularity assumption (76), and the Young’s inequality, we can prove the following estimate

[TABLE]

by some standard but tedious arguments which are analogous to the estimate of $\sum_{k=1}^{m}V_{1}^{k}(\partial_{\tau}{\theta_{\Psi}^{k}})$. Due to space limitations, we omit the proof here.

To estimate the term $\tau V_{4}^{k}(\partial_{\tau}{\theta_{\Psi}^{k}})$, we rewrite it by

[TABLE]

Arguing as before, we can obtain

[TABLE]

The real part of the last two terms on the right hand side of (134) can be decomposed as follows:

[TABLE]

Combining (134)-(136) and setting

[TABLE]

we obtain

[TABLE]

We now focus on the analysis of $\tau V_{5}^{k}(\partial_{\tau}{\theta_{\Psi}^{k}})$, which can be rewritten as

[TABLE]

By applying (82) and (106) and arguing as before, we deduce

[TABLE]

In order to estimate $|\sum\limits_{k=1}^{m}V_{5}^{k,3}|$, we rewrite it as follows.

[TABLE]

We decompose the term $T_{1}$ as follows.

[TABLE]

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

[TABLE]

Since

[TABLE]

we deduce

[TABLE]

by applying Theorem 4.1.

Hence we get the estimate of $T_{1}$ as follows.

[TABLE]

By virtue of (106) and integrating by parts, we discover

[TABLE]

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

[TABLE]

The last term at the right side of (147) satisfies the following estimate.

[TABLE]

Hence we get

[TABLE]

Reasoning as before, we can estimate $T_{3}$ as follows.

[TABLE]

Combining (146), (150), and (151) implies

[TABLE]

and thus

[TABLE]

Now substituting (130), (131), (133), (138), and (153) into (128), we have

[TABLE]

Arguing as in the proof of Theorem 2.2, we have

[TABLE]

Consequently, by inserting (155) into (154), we find

[TABLE]

By substituting (122) into (156), we end up with

[TABLE]

6.2 Estimates for (111)

Taking $\mathbf{v}=\frac{\displaystyle 1}{\displaystyle 2\tau}(\theta_{\mathbf{A}}^{k}-\theta_{\mathbf{A}}^{k-2})$ in (111), we see that

[TABLE]

which leads to

[TABLE]

Now we estimate $\sum\limits_{k=1}^{m}U_{i}^{k}(\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k})$, $i=1,2,3,4,5$. Under the regularity assumption of $\mathbf{A}$ in (76), we have

[TABLE]

Applying (106), the regularity assumption, and the Young’s inequality, we can bound $\tau\sum_{k=1}^{m}U_{2}^{k}(\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k})$ as follows.

[TABLE]

Since $\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k}\in\mathbf{X}_{0h}$, we have

[TABLE]

from which we deduce

[TABLE]

By applying (106) and reasoning as before, we have

[TABLE]

By applying Theorem 4.1 and the regularity assumption, the terms $\tau\sum_{k=1}^{m}\\ U_{4}^{k}(\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k})$ can be estimated by a standard argument.

[TABLE]

To estimate $\sum\limits_{k=1}^{m}U_{5}^{k}(\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k})$, we first rewrite $U_{5}^{k}(\overline{\partial_{\tau}\theta}_{\mathbf{A}}^{k})$ as follows.

[TABLE]

A simple calculation shows that

[TABLE]

and consequently, we have

[TABLE]

by applying (104) and (105).

Similarly, by employing (105), we deduce

[TABLE]

Therefore, we have

[TABLE]

Substituting (160), (161), (164), (165), and (170) into (159) and recalling the definition of $J_{1}^{k}$ in (127), we obtain

[TABLE]

6.3 Estimates for (112)

We can deduce the estimate of $\nabla{\theta}^{k}_{\phi}$ by a standard argument.

[TABLE]

From (112) we know that

[TABLE]

By taking $u=\overline{\theta}_{\phi}^{k}$ in (173) and recalling the definition of $J_{2}^{k}$ in (137), we find

[TABLE]

which implies that

[TABLE]

Employing the error estimates of interpolation operators and the regularity assumption, we obtain

[TABLE]

It follows that

[TABLE]

By combining (157), (171), and (177), we finally obtain

[TABLE]

which yields the desired estimate (107) by using the discrete Gronwall’s inequality.

7 numerical experiments

In this section, we present two numerical examples to confirm our theoretical analysis.

Example 1

To verify the conservation of total charge and energy of our scheme, we consider the M-S-C system (8)-(9) with the initial datas:

[TABLE]

In this example we take $\Omega=(0,1)^{3}$, $V(\mathbf{x})=5.0$, $T=4.0$, and the time step $\tau=0.01$.

The domain is partitioned into uniform tetrahedrals with $51$ nodes in each direction and $6\times 50^{3}$ elements in total. We solve the M-S-C system by the scheme (70) using mixed finite element method with $r=2$.

The evolutions of the total charge and energy of the discrete system are displayed in Fig.1, which clearly show that our algorithm almost exactly keeps the conservation of the total charge and energy of the system.

Example 2

We consider the following M-S-C system:

[TABLE]

with the exact solution $(\Psi,\mathbf{A},\phi)$

[TABLE]

where $V(\mathbf{x})=\frac{1}{2}(x_{1}^{2}+x_{2}^{2}+x_{3}^{2})$ and the right-hand side functions $g(\mathbf{x},t)$, $\mathbf{f}(\mathbf{x},t)$ and $h(\mathbf{x},t)$ are determined by the exact solution.

In this example, we set $\Omega=(0,1)^{3}$ and $T=2$. We take a uniform tetrahedral partition with $N+1$ nodes in each direction and $6N^{3}$ elements in total as in the example 1. The M-S-C system (179) are solved by the proposed scheme (70) with $r=1,2$, which are denoted by linear element method and quadratic element method, respectively. To test the convergence order of our scheme, we pick $\tau=h^{\frac{1}{2}}$ for the linear element method and $\tau=h$ for the quadratic element method respectively. We present numerical results for the linear element method and the quadratic element method at the final time $T=2.0$ in Tables 1 and 2, respectively. We can clearly see that the convergence rate of the quadratic element method agrees with our theoretical analysis while the linear element method has better convergence order than our theoretical analysis, which is partially because we use a quadratic element approximation of the vector potential $\mathbf{A}$.