Crank-Nicolson Galerkin method for Maxwell-Schrödinger Equations

# Error estimates of the Crank-Nicolson Galerkin method for the time-dependent Maxwell-Schrödinger equations under the Lorentz gauge

## Abstract

In this paper we study the numerical method and the convergence for solving the time-dependent Maxwell-Schrö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. error estimates; Crank-Nicolson; Galerkin method; Maxwell-Schrödinger.

\jno

drnxxx

\shortauthorlist

C. P. Ma et al.

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

 Missing or unrecognized delimiter for \big (1)

where , is a bounded Lipschitz polygonal convex domain in (or a bounded Lipschitz polyhedron convex domain in ), denotes the complex conjugate of , and respectively denote the electric permittivity and the magnetic permeability of the material and 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 , if is a solution of (1), then 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

 ∇⋅A+∂ϕ∂t=0.

(ii) The Coulomb gauge

 ∇⋅A=0.

(iii) The temporal gauge

 ϕ=0.

For simplicity, we employ the atomic units, i.e. , and assume that without loss of generality.

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

 Missing or unrecognized delimiter for \big (2)

The boundary conditions are

 ψ(x,t)=0,ϕ(x,t)=0,A(x,t)×n=0,∇⋅A(x,t)=0,(x,t)∈∂Ω×(0,T), (3)

and the initial conditions are

 ψ(x,0)=ψ0(x),ϕ(x,0)=ϕ0(x),ϕt(x,0)=ϕ1(x),A(x,0)=A0(x),At(x,0)=A1(x), (4)

where and denote the derivative of and with respect to the time respectively, is the outward unit normal to the boundary .

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

 ∇⋅A0+ϕ1=0,∇⋅A1+Δϕ0+|ψ0|2=0. (5)
{remark}

The boundary condition on implies that the particle is confined in a whole domain . The boundary condition and on are direct results of the perfect conductive boundary (PEC). The boundary condition on can be deduced from the boundary condition of and the gauge choice. As for the determination of the boundary conditions for the vector potential and the scalar potential , we refer to [Chew(2014)Chew].

The local and global well-posedness of solutions on all of for the time-dependent Maxwell-Schrödinger equations (1) have been studied in, for example, [Ginibre & Velo(2003)Ginibre & Velo], [Guo et al.(1995)Guo, Nakamitsu & Strauss], [Nakamitsu & Tsutsumi(1986)Nakamitsu & Tsutsumi], [Nakamura & Wada(2005)Nakamura & Wada], [Nakamura & Wada(2007)Nakamura & Wada], [Shimomura(2003)Shimomura] 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)Sui, Yang, Yun & Wang] 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)Pierantoni] 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)Ahmed & Li] 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 and magnetic fields , instead of the - formulation. Recently, [Ryu(2015)Ryu] 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)Ohnuki], [Sato & Yabana(2014)Sato & Yabana], [Turati & Hao(2012)Turati & Hao] 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)Gao, Li & Sun] 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 and are sufficiently small. Here and 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 and the error analysis for the wave function . 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) 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 we shall denote a positive constant independent of the mesh size and the time step 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 is a bounded Lipschitz polygonal convex domain in (or a bounded Lipschitz polyhedron convex domain in ). We introduce the following notation. Let denote the conventional Sobolev spaces of the real-valued functions. As usual, and are denoted by and respectively. We use and with calligraphic letters for Sobolev spaces of the complex-valued functions, respectively. Furthermore, let and with bold faced letters be Sobolev spaces of the vector-valued functions with components (=2, 3). inner-products in , and are denoted by without ambiguity.

In particular, we introduce the following subspace of :

 H1t(Ω)={A|A∈H1(Ω),A×n=0on∂Ω}

The semi-norm on is defined by

 ∥u∥H1t(Ω):=[∥∇⋅u∥2L2(Ω)+∥∇×u∥2L2(Ω)]12,

which is equivalent to the standard -norm , see [Girault & Raviart(1986)Girault & Raviart].

The weak formulation of the Maxwell-Schrödinger system (2)-(4) can be specified as follows: find such that ,

 Unknown environment '% (6)

with the initial conditions , , , and .

Let M be a positive integer and let be the time step. For any k=1,2,, we introduce the following notation:

 ∂Uk=(Uk−Uk−1)/Δt,∂2Uk=(∂Uk−∂Uk−1)/Δt,¯¯¯¯Uk=(Uk+Uk−1)/2,˜Uk=(Uk+Uk−2)/2, (7)

for any given sequence and denote for any given functions with a Banach space .

Let be a regular partition of into triangles in or tetrahedrons in without loss of generality, where the mesh size . For a given partition , let , and denote the corresponding -th order finite element subspaces of , and , respectively. Let , and be the conventional point-wise interpolation operators on , and , respectively.

For convenience, assume that the function and is defined in the interval in terms of the time variable . We can compute by

 A(⋅,−Δt)=A(⋅,0)−Δt∂A∂t(⋅,0)=A0−ΔtA1, (8)

which leads to an approximation to with second order accuracy. 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:

 ψ0h=Rhψ0,A0h=πhA0,A−1h=A0h−ΔtπhA1,ϕ0h=Ihϕ0,ϕ−1h=ϕ0h−ΔtIhψ1, (9)

and find such that for ,

 Unknown environment '% (10)

Note that , , , and are defined in (7), and denotes the complex conjugate of . For convenience, we define the following bilinear forms:

 B(A;ψ,φ)=((i∇+A)ψ,(i∇+A)φ),D(A,v)=(∇⋅A,∇⋅v)+(∇×A,∇×v),f(ψ,φ)=i2(φ∗∇ψ−ψ∇φ∗). (11)

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

 Unknown environment '% (12)

and for ,

 Extra open brace or missing close brace (13)

In this paper we assume that the Maxwell-Schrödinger equations (12) has one and only one weak solution and the following regularity conditions are satisfied:

 ψ,ψt,ψtt∈L∞(0,T;Hr+1(Ω)),ψttt∈L∞(0,T;H1(Ω)),ψtttt∈L2(0,T;L2(Ω));A,At,Att∈L∞(0,T;Hr+1(Ω)),Attt∈L∞(0,T;H1(Ω)),Atttt∈L2(0,T;L2(Ω));ϕ,ϕt,ϕtt∈L∞(0,T;Hr+1(Ω)),ϕttt∈L∞(0,T;H1(Ω)),ϕtttt∈L2(0,T;L2(Ω)). (14)

For the initial conditions , we assume that

 ψ0∈Hr+1(Ω)∩H10(Ω);A0,A1∈Hr+1(Ω)∩H1t(Ω);ϕ0,ϕ1∈Hr+1(Ω)∩H10(Ω). (15)

We now give the main convergence result in this paper as follows: {theorem} Suppose that the Maxwell-Schrödinger coupled system (12) has a unique solution satisfying (14) and (15). Let be the fully discrete numerical solution of defined in (13). Then there exist two positive constants and , such that when , , we have the following error estimates:

 max1≤k≤M{∥ψkh−ψk∥2L2(Ω)+∥∇(ψkh−ψk)∥2L2(Ω)+∥Akh−Ak∥2L2(Ω)+∥∇⋅(Akh−Ak)∥2L2(Ω)+∥∇×(Akh−Ak)∥2L2(Ω)+∥ϕkh−ϕk∥2L2(Ω)+∥∇(ϕkh−ϕk)∥2L2(Ω)}≤C{h2r+(Δt)4},r≥1, (16)

where , , , and is a constant independent of , .

## 3 The proof of Theorem 2

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

### 3.1 Preliminaries

For convenience, we list some imbedding inequalities and interpolation inequalities in Sobolev spaces (see, e.g., [Ladyzhenskaya et al.(1968)Ladyzhenskaya, Solonnikov & Uralceva] and [Girault & Raviart(1986)Girault & Raviart]), and use them in the sequel.

 ∥u∥Lp≤C∥u∥H1,∥v∥Lp≤C∥v∥H1,1≤p≤6(d=2,3), (17)
 ∥v∥H1≤C(∥∇×v∥L2+∥∇⋅v∥L2),v∈H1t(Ω), (18)
 ∥u∥2L3≤∥u∥L2∥u∥L6, (19)

where , , and .

The following identities will be used frequently in this paper.

 M∑k=1(ak−ak−1)bk=aMbM−a0b1−M−1∑k=1ak(bk+1−bk),M∑k=1(ak−ak−1)bk=aMbM−a0b0−M∑k=1ak−1(bk−bk−1). (20)

Let denote the interpolation function of in . Set , , . By applying standard finite element theory and the regualrity conditions in (14), we have

 ∥eψ∥L2+h∥eψ∥H1≤Chr+1,∥eA∥L2+h∥eA∥H1≤Chr+1,∥eϕ∥L2+h∥eϕ∥H1≤Chr+1,∥Rhψ∥L∞+∥πhA∥L∞+∥Ihϕ∥L∞+∥∇Rhψ∥L3≤C, (21)

where is a constant independent of .

The following lemmas will be useful in the proof of Theorem 2. {lemma} For the solution of (13), we have

 ∥ψkh∥2L2=∥ψ0h∥2L2,k=1,2⋯,M. (22)
{lemma}

For , the following identities hold for the bilinear functional defined in (11):

 Extra open brace or missing close brace (23)

Lemma 3.1 can be proved by choosing in and taking the imaginary part. A direct calculation gives (23) in Lemma 3.1.

Let , , . By using the error estimates of the interpolation operators (21), we only need to estimate , and . We will prove the following estimate:

 ∥θkψ∥2H1+∥∂θkA∥2L2+∥θkA∥2H1+∥∂θkϕ∥2L2+∥θkϕ∥2H1≤C∗{h2r+(Δt)4}, (24)

where .

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

 ∥θ0ψ∥2H1+∥∂θ0A∥2L2+∥θ0A∥2H1+∥∂θ0ϕ∥2L2+∥θ0ϕ∥2H1≤C0{h2r+(Δt)4}. (25)

We use the mathematical inductive method to show (24). By (25), if we require , then (24) holds for . We assume that (24) holds for . In the rest of this section, we will find such that (24) holds for , where is independent of , , .

Subtracting (12) from (13), we obtain the following equations for , and :

 (∂2θkA,v)+D(˜θkA,v)=(∂2Ak−1∂t2−∂2πhAk,v)+D(Ak−1−˜πhAk,v)+(|ψk−1|2Ak−1−|ψk−1h|2˜Akh,v)+(f(ψk−1,ψk−1)−f(ψk−1h,ψk−1h),v),∀v∈Vrh, (26)
 (∂2θkϕ,η)+(∇˜θkϕ,∇η)=(∂2ϕk−1∂t2−∂2Ihϕk,η)+(∇(ϕk−1−˜Ihϕk),∇η)+(|ψk−1h|2−|ψk−1|2,η),∀η∈Vrh, (27)
 −2i(∂θkψ,φ)+B(¯¯¯¯¯Akh;¯¯¯θkψ,φ)=2i⎛⎝∂Rhψk−∂ψk−12∂t,φ⎞⎠+2V0(ψk−12−¯¯¯¯ψkh,φ)+B(Ak−12;(ψk−12−Rh¯¯¯¯ψk),φ)+2(ϕk−12ψk−12−¯¯¯ϕkh¯¯¯¯ψkh,φ)+(B(Ak−12;Rh¯¯¯¯ψk,φ)−B(¯¯¯¯¯Akh;Rh¯¯¯¯ψk,φ)),∀φ∈Vrh. (28)

The key steps of the proof of (24) are now briefly described. In order to find and to show that (24) holds for , we first take in (26) and in (27) and obtain the estimates of and . Then we choose in (28) and give the estimate of . Finally, we take in (28) and make use of the above estimates of and to derive the energy-norm estimate for . Using the above estimates, we can complete the proof of (24).

### 3.2 Estimates for (26)

If we set

 Ik1(v)=(∂2Ak−1∂t2−∂2πhAk,v),Ik2(v)=D(Ak−1−˜πhAk,v),Ik3(v)=(|ψk−1|2Ak−1−|ψk−1h|2˜Akh,v),Ik4(v)=(f(ψk−1,ψk−1)−f(ψk−1h,ψk−1h),v),v∈Vrh, (29)

then we rewrite (26) as follows:

 (∂2θkA,v)+D(˜θkA,v)=Ik1(v)+Ik2(v)+Ik3(v)+Ik4(v). (30)

We take in (30) and get

 (∂2θkA,12(∂θkA+∂θk−1A))+D(˜θkA,12(∂θkA+∂θk−1A))=12Δt[∥∂θkA∥2L2−∥∂θk−1A∥2L2]+14Δt[D(θkA,θkA)−D(θk−2A,θk−2A)]=Ik1(¯¯¯¯¯¯∂θkA)+Ik2(¯¯¯¯¯¯∂θkA)+Ik3(¯¯¯¯¯¯∂θkA)+Ik4(¯¯¯¯¯¯∂θkA), (31)

which leads to

 12∥∂θmA∥2L2+14D(θmA,θmA)+14D(θm−1A,θm−1A)=12∥∂θ0A∥2L2+14D(θ0A,θ0A)+14D(θ−1A,θ−1A)+Δtm∑k=1[Ik1(¯¯¯¯¯¯∂θkA)+Ik2(¯¯¯¯¯¯∂θkA)+Ik3(¯¯¯¯¯¯∂θkA)+Ik4(¯¯¯¯¯¯∂θkA)]≤Ch2r+Δtm∑k=1[Ik1(¯¯¯¯¯¯∂θkA)+Ik2(¯¯¯¯¯¯∂θkA)+Ik3(¯¯¯¯¯¯∂θkA)+Ik4(¯¯¯¯¯¯∂θkA)]. (32)

Here we have used the fact that

 ∥∂θ0A∥2L2+D(θ0A,θ0A)+D(θ−1A,θ−1A)≤Ch2r.

Now we estimate , . Under the regularity assumption of in (14), we can prove

 Missing or unrecognized delimiter for \right (33)

We rewrite the term as follows:

 Δtm∑k=1Ik2(¯¯¯¯¯¯∂θkA)=12Δtm∑k=1D(Ak−1−12(Ak+Ak−2),(∂θkA+∂θk−1A))+14Δtm∑k=1D((Ak+Ak−2)−πh(Ak+Ak−2),(∂θkA+∂θk−1A)). (34)

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

 12Δtm∑k=1D(Ak−1−12(Ak+Ak−2),(∂θkA+∂θk−1A))≤C{h2r+(Δt)4}+164D(θmA,θmA)+164D(θm−1A,θm−1A)+CΔtm∑k=0D(θkA,θkA) (35)

and

 14Δtm∑k=1D((Ak+Ak−2)−πh(Ak+Ak−2),(∂θkA+∂θk−1A))≤C{h2r+(Δt)4}+164D(θmA,θmA)+164D(θm−1A,θm−1A)+CΔtm∑k=0D(θkA,θkA). (36)

We thus have

 Missing \left or extra \right (37)

We observe that

 Ik3(¯¯¯¯¯¯∂θkA)=(|ψk−1|2(Ak−1−πh˜Ak),¯¯¯¯¯¯∂θkA)+(|ψk−1|2(πh˜Ak−˜Akh),¯¯¯¯¯¯∂θkA)+((|ψk−1|2−|Rhψk−1|2)πh˜Ak,¯¯¯¯¯¯∂θkA)+((|ψk−1|2−|Rhψk−1|2)˜θkA,¯¯¯¯¯¯∂θkA) +((|Rhψk−1|2−|ψk−1h|2)πh˜Ak,¯¯¯¯¯¯∂θkA)+((|