# Low-regularity integrators for nonlinear Dirac equations

In this work, we consider the numerical integration of the nonlinear Dirac equation and the Dirac-Poisson system (NDEs) under rough initial data. We propose a ultra low-regularity integrator (ULI) for solving the NDEs which enables optimal first-order time convergence in H^r for solutions in H^r, i.e., without requiring any additional regularity on the solution. In contrast to classical methods, ULI overcomes the numerical loss of derivatives and is therefore more efficient and accurate for approximating low regular solutions. Convergence theorems and the extension of ULI to second order are established. Numerical experiments confirm the theoretical results and underline the favourable error behaviour of the new method at low regularity compared to classical integration schemes.

## Authors

• 17 publications
• 152 publications
• 10 publications
08/17/2020

### Embedded exponential-type low-regularity integrators for KdV equation under rough data

In this paper, we introduce a novel class of embedded exponential-type l...
07/31/2020

### Numerical integrators for continuous disordered nonlinear Schrödinger equation

In this paper, we consider the numerical solution of the continuous diso...
10/16/2019

### Optimal convergence of a second order low-regularity integrator for the KdV equation

In this paper, we establish the optimal convergence result of a second o...
03/29/2022

### A second-order low-regularity correction of Lie splitting for the semilinear Klein–Gordon equation

The numerical approximation of the semilinear Klein–Gordon equation in t...
09/02/2021

### A second order low-regularity integrator for the nonlinear Schrödinger equation

In this paper, we analyse a new exponential-type integrator for the nonl...
07/30/2019

### On the lack of interior regularity of the p-Poisson problem with p>2

In this note we are concerned with interior regularity properties of the...
12/27/2021

### Implicit regularity and linear convergence rates for the generalized trust-region subproblem

In this paper we develop efficient first-order algorithms for the genera...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Numerical integrators for solving the semi-linear dispersive equation

 ∂tu=Lu+N(u),t>0,

usually require smoothness of the solution , i.e., the boundedness of in some Sobolev space, where

denotes a skew-adjoint linear differential operator and

denotes some nonlinear function [52]. In particular, traditional methods, such as explicit and implicit Runge-Kutta methods, splitting schemes and exponential integrators [30], can only reach their optimal convergence rate in () for solutions in for some depending on the order of spatial differentiation in and . This additional regularity requirement on the solution (i.e., ) is indeed introduced by the numerical approximation and necessary for (optimal) convergence. Recently, for certain problems a new class of integrators could be constructed which allow optimal convergence without any loss of numerical derivatives in the time discretization (i.e., ), see for instance [41] for the one-dimensional quadratic nonlinear Schrödinger equation. In this work, we identify the nonlinear Dirac equations, i.e., a class of nonlinear dispersive equations with cubic nonlinearities, for which such low-regularity integrators can be designed. As important models in particle physics and relativistic quantum mechanics, the Dirac-type equations have been widely considered in studies of two-dimensional materials [11, 23], quantum field theory [46, 51] and Bose-Einstein condensates [26]. We shall consider in this paper the following one-dimensional nonlinear Dirac equation (NDE) [12, 35] as the model problem:

 i∂tΦ=−iα∂xΦ+βΦ+VΦ+F(Φ),t>0, x∈R, (1.1a) Φ(0,x)=Φ0(x),x∈R, (1.1b)

where is the unknown spinorfield, is the initial data, denotes the electric potential which is either given as an external real-valued function or a time-dependent function determined by a self-consistent Poisson equation [11]

 −∂xxV=|Φ|2,t≥0, x∈Rwith∫RVdx≡0,t≥0,

denotes the cubic Thirring-type nonlinearity with a given parameter denoting the strength of the nonlinear interaction [46], where denotes the complex conjugate transpose of , and

 α=(0110),β=(100−1),

are the Pauli matrices.

Equation (1.1) occurs as a core part of many related models, such as the Gross-Neveu system [33], the Dirac-Klein-Gordon system [16, 53], the Maxwell-Dirac system [18, 32] and the Chern-Simons-Dirac system [10]. Due to these applications, the numerical solutions of the NDE (1.1) are widely interested in computational physics [4, 11, 24, 32, 35, 40, 45], especially the bound states, the solitary waves and the interaction dynamics of solitons [1, 29, 55]. These existing works in computational physics or applied mathematics mainly assume that the initial data of NDE is smooth enough for theoretical and/or numerical studies. While in general physical situations or real applications with unavoidable noise, the initial input in (1.1b) for the NDE may not be a smooth function or a function with high regularity. For such initial data, mathematical analysis of NDE has been intensively carried out in the literature. The global existence of the NDE (1.1) was firstly established in [19] for initial data, i.e., . Later, many efforts [5, 7, 15, 37, 38, 48] were made to obtain the local and/or global well-posedness of (1.1) by allowing initial data with lower regularities. Among them, Selberg and Tesfahun proved the local well-posedness for initial data in almost critical space with any [48], and Candy extended their result to the critical case for any [15]. We refer interested readers to [14]

for a systematical review and some refined Strichartz estimates. Corresponding studies have also been made in high space dimensions case

[8, 20, 44] and in coupling case with Klein-Gordon equation or Maxwell equations et al. [9, 10, 16, 18, 21, 33, 43, 47].

However, the numerical methods proposed for solving the NDE in the literature so far all lead to numerical loss of regularity. More precisely, for the numerical solutions to have a -th () order of convergence in time in the Sobolev space , the finite difference time domain methods in [11, 28, 45] and the Runge-Kutta discontinuous Galerkin methods in [49] both require initial data at least in , and the splitting spectral methods in [4, 24, 32, 40] and the exponential integrator spectral methods in [2, 3, 54] require initial data in . We refer to [55] for a detailed numerical comparison of these methods with smooth initial data. For initial data with less regularity than the required, these numerical schemes will fail to reach their optimal convergence rates, suffering from sever order reduction and hence become less efficient and accurate (see also Section 5 for a numerical illustration of this phenomenon). We shall summarize and emphasize the limitations of these methods in the paper by writing down their rigorous error estimates with critical regularity requirements.

To enable efficient numerical integrations of NDE with rough initial data, we shall present a class of integrators in the spirit of low-regularity methods [31, 41, 42], for solving (1.1). The methods are derived under the framework of the nested Picard iterative integrators (NPI) [13], while we modify the evolution operator, instead of using the one with “free Dirac operator” [6, 13]

. This allows for a decomposition in eigenspaces which ensures exact integrations of terms in the Fourier frequency space in an explicit and efficient way.

The main novelty of the paper is to introduce a new class of ultra low-regularity integrator (ULI) integrators for Dirac equations. This new class will offer optimal first-order time convergence in for any -data (for ), i.e., no auxiliary smoothness is needed from the solution at all if one does not take the spatial discretization into account. Therefore, we say that the scheme is a ultra low-regularity integrator. In the literature, such ULI methods could only be designed so far for the quadratic Schrödinger equation [41]. On the other hand, for the cubic nonlinear Schrödinger equations or the KdV equation, some addition regularity requirement has to be imposed [31, 34, 41, 42].

Based on the first-order ULI method, we build another scheme that offers optimal second-order convergence in for -data (). In particular, second-order convergence in time holds in for solutions in . Up to our knowledge no scheme which allows for such a generous convergence has been proposed in literature so far and compared to the classical Strang splitting scheme, one space derivative is saved by our new second-order method. Note that the requirement that the solution is in can be seen naturally as the spatial convergence requires some smoothness of the solution. We focus on first- and second-order integrators. Our new framework can, however, be (in principle) generalized to arbitrary high-order low-regularity integrators. Rigorous convergence results are established and numerical results are presented to underline the performance of the proposed methods in comparison to classical schemes.

The rest of this paper is organized as follows. In Section 2, we review some standard numerical methods for solving the NDE (1.1). Thereby we pay particular attention on the required regularities. In Section 3, we present the first-order ULI scheme and its convergence result. Extensions of ULI to the second order are made in Section 4. Numerical results are presented in Section 5 and some conclusions are drawn in Section 6.

## 2. Convergence of standard methods

In this section, we shall review some of the popular integrators including finite difference methods, exponential integrators and splitting methods for NDE (1.1). For simplicity of presentation, we will assume that is a given external smooth function in this section and we will focus on the time integrations. As for spatial discretization, finite difference methods [11], discontinuous Galerkin methods [49] and spectral methods [12, 35] can all be applied. To simplify numerical implementations, we truncate the whole space problem (1.1) onto the torus :

 i∂tΦ=−iα∂xΦ+βΦ+VΦ+F(Φ),t>0, x∈T, (2.1a) Φ(0,x)=Φ0(x),x∈T, (2.1b)

and impose periodic boundary conditions so that Fourier spectral/pseudospectral method can be easily applied. To simplify the notations, we shall omit the space variable in the following, i.e., we will write . We denote by the time step, as the time grids, and define as the numerical solution. We write

for the standard Sobolev norm, and as for a vector field

on , we set .

For each of the reviewed standard methods in the following, we write down their convergence theorems to address their critical regularity requirements for convergence. The proofs of the theorems are given in Appendix A.

### 2.1. Finite difference methods

As the most traditional numerical discretization, finite difference methods have been widely applied for solving the NDE [3, 28, 45, 49] and coupled systems, such as the Dirac-Poisson and Klein-Gordon-Dirac systems [11, 53]. Here, we present two semi-implicit finite difference integrators which are free from CFL conditions.

A first-order semi-implicit finite difference integrator (FD1) for the NDE (2.1) reads:

 iΦn+1−Φnτ=−iα∂xΦn+1+βΦn+VΦn+F(Φn),n≥0, (2.2a) Φ0=Φ0, (2.2b)

and a second-order semi-implicit finite difference integrator (FD2) reads:

 iΦn+1−Φn−12τ=−iα∂xΦn+1+Φn−12+βΦn+VΦn+F(Φn),n≥0, (2.3a) Φ0=Φ0,Φ1=Φ0−iτ[−iα∂xΦ0+βΦ0+VΦ0+F(Φ0)], (2.3b)

where in the second-order scheme is obtained by a first-order Taylor expansion of equation (2.1a) at .

###### Theorem 2.1.

(Convergence of finite difference methods) Let denote the numerical solution of the FD1 scheme (2.2) for solving the NDE (2.1). Let and where , for some . Then there exist constants depending on , , and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ.

Further under assumption where , for from the FD2 scheme (2.3), there exist constants depending on , , , and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ2.

### 2.2. Classical exponential integrators

Exponential integrators [30] have been intensively developed and analyzed for various evolution equations, particularly for equations involving a stiff linear term. They have been recently considered to solve NDEs in the nonrelativistic limit regime in [2, 3, 54]. In the following we recall their construction in case of the Dirac equation.

By writing the NDE (2.1) for under the Duhamel’s formula, we have

 Φ(tn+1)=e−iτTΦ(tn)−i∫τ0e−i(τ−s)TG(Φ(tn+s))ds,n≥0, (2.4)

with

 T:=−iα∂x+β,G(Φ(tn+s))=VΦ(tn+s)+F(Φ(tn+s)).

Here is known as the free Dirac operator [6]. By approximating in the integrant with , the first-order classical Gautschi-type exponential integrator (EI1) [25] reads

 Φn+1=e−iτTΦn−iτφ1(−iτT)Gn, n≥0, (2.5)

with and

 φ1(z)=ez−1z. (2.6)

On the other hand, by extrapolation as in [3], the second-order Gautschi-type exponential integrator (EI2) reads

 Φn+1=e−iτTΦn−iτφ1(−iτT)Gn−iτφ2(−iτT)(Gn−Gn−1),n≥1, (2.7)

with .

###### Theorem 2.2.

(Convergence of exponential integrators) Let denote the numerical solution of the EI1 scheme (2.5) for solving the NDE (2.1). Let and where , for some . Then there exist constants depending on , and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ.

Further under assumption where , then for from the EI2 scheme (2.7), there exist constants depending on , , and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ2.

### 2.3. Splitting methods

The time splitting or operator splitting method can be considered as one of the most popular numerical techniques for solving evolutions equations of first order in time [39]. They have been proposed in the literature for the NDEs [3, 24, 40] and for Maxwell-Dirac system [4, 32]. The idea of the time splitting methods is to split (2.1) into the two systems

 Φkin(t):i∂tΦ=−iα∂xΦ+βΦ,t>0, x∈T,

and

 Φpon(t):i∂tΦ=VΦ+F(Φ),t>0, x∈T.

Both of the resulting flows can be exactly integrated in time, i.e.,

 Φkin(t): Φ(t)=e−itTΦ(0),Φpon(t): Φ(t)=e−it(V⋅Id+λ(|ϕ1|2−|ϕ2|2)β)Φ(0),

where is the identity matrix. The two flows are then combined as for the Lie splitting scheme:

 Φn−=e−iτ(V⋅Id+λ(|ϕn1|2−|ϕn2|2)β)Φn,n≥0, (2.8a) Φn+1=e−iτTΦn−, (2.8b)

or combined as for the Strang splitting scheme:

 Φn−=e−iτ2(V⋅Id+λ(|ϕn1|2−|ϕn2|2)β)Φn,n≥0, (2.9a) Φn+=e−iτTΦn−, (2.9b) Φn+1=e−iτ2(V⋅Id+λ(|ϕn+,1|2−|ϕn+,2|2)β)Φn+, (2.9c)

with .

###### Theorem 2.3.

(Convergence of splitting methods) Let denote the numerical solution of the Lie splitting scheme (2.8) for solving the NDE (2.1). Let and for some . Then there exist constants depending on and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ.

Further under assumption , then for from the Strang splitting scheme (2.9), there exist constants depending on and , such that for all and , we have

 ∥Φ(tn)−Φn∥r≤Cτ2.

## 3. First-order ultra low-regularity integrators

In this section, we derive the ultra low-regularity integrators of first-order accuracy for solving the NDEs and present their convergence theorems. We shall consider in a sequel the case of an external electrical potential and the case of Dirac-Poisson system.

### 3.1. NDE with external field

We firstly consider the NDE (2.1) with a given function , under the same notations as in Section 2. Similar to exponential integrators in Section 2.2, our method is also constructed based on the integral form of the NDE. However, we apply Duhamel’s formula for (2.1) for in an alternative way as follows

 Φ(tn+s)=e−sα∂xΦ(tn)−i∫s0e−(s−ρ)α∂x[βΦ(tn+ρ)+VΦ(tn+ρ)+F(Φ(tn+ρ))]dρ. (3.1)

Here, we use the evolution operator instead of in (2.4), which is crucial to the success of the method. It is known that the ‘free Dirac operator’ is diagonalizable in Fourier space and it can be decomposed as [6]

 T=√Id−∂xx ΠT+−√Id−∂xx ΠT−,

with the projectors defined as

 ΠT±=12[Id±(Id−∂xx)−1/2T].

As can be seen that the spatial differentiation operator is involved nonlinearly in the above decomposition. This works well in the numerical studies for the nonrelativistic limit regime [12, 13], but here it makes things difficult in designing low-regularity methods. In contrast, the operator can be decomposed as

 α∂x=∂xΠ+−∂xΠ− (3.2)

with

 Π+=12(1111),Π−=12(1−1−11),

where is linearly involved. It can be verified that , and

 (α∂x)k=(∂x)kΠ++(−∂x)kΠ−,k∈N,

which implies the following relation

 esα∂x=es∂xΠ++e−s∂xΠ−. (3.3)

Therefore, a nested Picard iteration based on (3.1) can be effectively computed in analogous to [13] as follows.

Letting for in (3.1), we find

 Φ(tn+ρ)=e−ρα∂xΦ(tn)+O(τ). (3.4)

Here (and in the following) denotes a remainder of order in time which does not require any additional regularity of the solution, i.e.,

 Φ(tn+ρ)=e−ρα∂xΦ(tn)+O(τ) if ∥Φ(tn+ρ)−e−ρα∂xΦ(tn)∥r≤cτ

with depending on and .

Plugging (3.4) into (3.1) and setting , we get

 Φ(tn+1)=e−τα∂xΦ(tn)−ie−τα∂x[I1(tn)+I2(tn)+I3(tn)]+O(τ2), (3.5)

where

 I1(tn):=∫τ0eρα∂xβe−ρα∂xΦ(tn)dρ,I2(tn):=∫τ0eρα∂xVe−ρα∂xΦ(tn)dρ,I3(tn):=∫τ0eρα∂xF(e−ρα∂xΦ(tn))dρ. (3.6)

For , by using the relation (3.3), we have

 eρα∂xβe−ρα∂xΦ(tn)=(eρ∂xΠ++e−ρ∂xΠ−)β(e−ρ∂xΠ++eρ∂xΠ−)Φ(tn)=(e2ρ∂xβΠ−+e−2ρ∂xβΠ+)Φ(tn),

since and . Therefore, we have the exact integration:

 I1(tn)=τ[φ1(2τ∂x)βΠ−+φ1(−2τ∂x)βΠ+]Φ(tn), (3.7)

where is defined in (2.6).

For , we can analogously write the integrand as

 eρα∂xVe−ρα∂xΦ(tn)=(eρ∂xVe−ρ∂xΠ++e−ρ∂xVeρ∂xΠ−)Φ(tn).

What is interesting is that the integration of the above function can be performed in an exact and explicit way in the physical space. For some general function , we have

 =∑l∈Z∑l1,l2∈Zl1+l2=leilx∫τ0ei(l−l2)ρdρˆVl1ˆψl2=∑l∈Z∑l1,l2∈Zl1+l2=leilxτφ1(iτl1)ˆVl1ˆψl2 =τ(φ1(τ∂x)V)ψ,

where denotes the Fourier coefficients of , and similarly

 ∫τ0e−ρ∂xVeρ∂xψdρ=τ(φ1(−τ∂x)V)ψ.

Hence, we have the exact integration for as

 I2(tn)=τ[(φ1(τ∂x)V)Π++(φ1(−τ∂x)V)Π−]Φ(tn). (3.8)

This exact integration is expressed explicitly in the physical space so that in practice it can be obtained efficiently by means of fast Fourier transform, which is crucial for the success of our method.

For the integration of the nonlinear term in , firstly we have

 g:= (e−ρα∂xΦ(tn))∗β(e−ρα∂xΦ(tn))= [e−ρ∂x(Φ(tn))∗eρ∂xβΠ−+eρ∂x(Φ(tn))∗e−ρ∂xβΠ+]Φ(tn)= 12[(eρ∂xϕ−(tn))(e−ρ∂x¯¯¯¯¯¯¯ϕ+(tn))+(eρ∂x¯¯¯¯¯¯¯ϕ−(tn))(e−ρ∂xϕ+(tn))],

where

 ϕ±(tn):=ϕ1(tn)±ϕ2(tn),Φ(tn)=(ϕ1(tn),ϕ2(tn))T.

 I3(tn)=∫τ0eρα∂xgβe−ρα∂xΦ(tn)dρ=∫τ0[eρ∂xgeρ∂xβΠ−+e−ρ∂xge−ρ∂xβΠ+]Φ(tn)dρ,

which as we shall see can also be done explicitly in the physical space. Let us consider three general complex-valued scalar functions , and we find

 ∫τ0eρ∂x[(eρ∂xψ)(e−ρ∂xϕ)eρ∂xφ]dρ=∫τ0∑l∈Zeilxeilρ∑l1,l2,l3∈Zl1+l2+l3=leil1ρˆψl1e−il2ρˆϕl2eil3ρˆφl3dρ =∑l∈Z∑l1,l2,l3∈Zl1+l2+l3=leilx∫τ0e2i(l1+l3)ρdρˆψl1ˆϕl2ˆφl3=τϕφ1(2τ∂x)(ψφ).

Similarly,

 ∫τ0e−ρ∂x[(eρ∂xψ)(e−ρ∂xϕ)e−ρ∂xφ]dρ=τψφ1(−2τ∂x)(ϕφ).

Hence, we have for :

 I3(tn)=τλ2[¯¯¯¯¯¯¯ϕ+(tn)φ1(2τ∂x)(ϕ−(tn)βΠ−Φ(tn))+ϕ+(tn)φ1(2τ∂x)(¯¯¯¯¯¯¯ϕ−(tn)βΠ−Φ(tn))+ϕ−(tn)φ1(−2τ∂x)(¯¯¯¯¯¯¯ϕ+(tn)βΠ+Φ(tn))+¯¯¯¯¯¯¯ϕ−(tn)φ1(−2τ∂x)(ϕ+(tn)βΠ+Φ(tn))]. (3.9)

In summary of (3.5)-(3.9), the detailed scheme of the first-order ultra low-regularity integrator (ULI) for integrating NDE (2.1) with given reads: denote for , let and then

 Φn+1= e−τα∂xΦn−ie−τα∂x(In1+In2+In3)=:Θext(Φn),n≥0, (3.10)

where

 In1=τ[φ1(2τ∂x)βΠ−+φ2(−2τ∂x)βΠ+]Φn,In2=τ[(φ1(τ∂x)V)Π++(φ1(−τ∂x)V)Π−]Φn, In3=τλ4(¯¯¯¯¯¯¯ϕn+φ1(2τ∂x)(ϕn−)2+ϕn+φ1(2τ∂x)|ϕn−|2+ϕn−φ1(−2τ∂x)|ϕn+|2+¯¯¯¯¯¯¯ϕn−φ1(−2τ∂x)(ϕn+)2¯¯¯¯¯¯¯ϕn+φ1(2τ∂x)(ϕn−)2+ϕn+φ1(2τ∂x)|ϕn−|2−ϕn−φ1(−2τ∂x)|ϕn+|2−¯¯¯¯¯¯¯ϕn−φ1(−2τ∂x)(ϕn+)2), (3.11)

with .

The proposed ULI scheme, i.e., (3.10) with (3.11), is fully explicit. In practical computations, the spatial discretization of ULI could easily be done by the Fourier pseudospectral method [50], where the computational cost at each time level is with the number of the total Fourier modes. Thus, the ULI scheme is of similar computational costs as the standard methods in Section 2.

### 3.2. NDE with consistent field

Next, we present the integration for the NDE coupled with the Poisson equation:

 i∂tΦ=−iα∂xΦ+βΦ+VΦ+F(Φ),t>0, x∈T, (3.12a) −∂xxV=|Φ|2,t≥0, x∈T,∫TVdx=0,t≥0, (3.12b) Φ(0)=Φ0,x∈T, (3.12c)

where periodic boundary conditions are imposed for both and .

The Duhamel’s formula for (3.12) reads as

 Φ(tn+s)=e−sα∂xΦ(tn)−i∫s0e−(s−ρ)α∂x[βΦ(tn+ρ)+V(tn+ρ)Φ(tn+ρ)+F(Φ(tn+ρ))]dρ. (3.13)

Let in (3.13), and by adopting as before, the approximation goes the same as (3.5) but with (re-)defined as:

 I2(tn):=−∫τ0eρα∂x(∂−1xx|e−ρα∂xΦ(tn)|2)e−ρα∂xΦ(tn)dρ. (3.14)

Here, we define for some general function ,

 ∂−1xxψ(x):=−∑l≠01l2ˆψleilx,

where denotes the natural  inverse operator of .

The integrations of