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


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.




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.


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

(ii) The Coulomb gauge

(iii) The temporal gauge

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:


The boundary conditions are


and the initial conditions are


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:


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 :

The semi-norm on is defined by

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 ,


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:


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


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:


and find such that for ,


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


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


and for ,


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:


For the initial conditions , we assume that


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:


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.


where , , and .

The following identities will be used frequently in this paper.


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


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


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


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:


where .

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


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 :


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


then we rewrite (26) as follows:


We take in (30) and get


which leads to


Here we have used the fact that

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


We rewrite the term as follows:


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




We thus have


We observe that