Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

Florian M Éhats
email: florian.mehats@univ-rennes1.fr

Xiaofei Zhao

UNIFORMLY ACCURATE NUMERICAL SCHEMES FOR THE NONLINEAR DIRAC EQUATION IN THE NONRELATIVISTIC LIMIT REGIME

Keywords: nonlinear Dirac equation, nonrelativistic limit, highly oscillatory equations, uniform accuracy, two-scale formulation AMS subject classifications. 35Q55, 65M12, 74Q10

We apply the two-scale formulation approach to propose uniformly accurate (UA) schemes for solving the nonlinear Dirac equation in the nonrelativistic limit regime. The nonlinear Dirac equation involves two small scales ε and ε 2 with ε → 0 in the nonrelativistic limit regime. The small parameter causes high oscillations in time which brings severe numerical burden for classical numerical methods. We transform our original problem as a two-scale formulation and present a general strategy to tackle a class of highly oscillatory problems involving the two small scales ε and ε 2 . Suitable initial data for the two-scale formulation is derived to bound the time derivatives of the augmented solution. Numerical schemes with uniform (with respect to ε ∈ (0,1]) spectral accuracy in space and uniform first order or second order accuracy in time are proposed. Numerical experiments are done to confirm the UA property.

1. Introduction The nonlinear Dirac equation has been widely considered in many physical and mathematical areas, such as the electron self-interacting [START_REF] Gross | Dynamical sysmmetry breaking in asymptotically free field theories[END_REF], the gravity theory [START_REF] Kibble | Lorentz invariance and the gravitational field[END_REF], and recent studies in material graphene and Bose-Einstein condensates [START_REF] Ablowitz | Nonlinear waves in shallow honeycomb lattices[END_REF][START_REF] Fefferman | Honeycomb lattice potentials and Dirac points[END_REF][START_REF] Haddad | The nonlinear Dirac equation in Bose-Einstein condensates: superfluid fluctuations and emergent theories from relativistic linear stability equations[END_REF]. In this paper, we consider the nonlinear reduced Dirac equation

i∂ t Φ ε = - i ε α∂ x Φ ε + 1 ε 2 βΦ ε + [V e (t) + V m (t)α]Φ ε + λ(βΦ ε ,Φ ε )βΦ ε , (1.1)
for t > 0, x ∈ R, where the unknown Φ ε = (φ ε 1 ,φ ε 2 ) T = Φ ε (t,x) is a bi-dimensional complexvalued column vector and is interpreted as the wave function subject to the initial condition Φ ε (0,x) = Φ 0 (x),

x ∈ R.

λ ∈ R denotes the coupling constant and ε ∈ (0,1] is a dimensionless parameter which is inversely proportional to the speed of light. V e (t) = V e (t,x) and V m (t) = V m (t,x) are two given real-valued scalar functions representing the electrical potential and magnetic potential [START_REF] Bechouche | Semi-nonrelativistic limits of the Dirac equation with external time-dependent electromagnetic field[END_REF][START_REF] Huang | A time-splitting spectral scheme for the Maxwell-Dirac system[END_REF][START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF], respectively, and the matrices α and β are known as the Pauli matrices, i.e.

α = 0 1 1 0 , β = 1 0 0 -1 .
In particular, we have (βΦ

ε ,Φ ε ) = |φ ε 1 | 2 -|φ ε 2 | 2 .
The nonlinear Dirac equation (1.1) has been widely considered in the literature [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bournaveas | Theory and numerical approximations for a nonlinear 1 + 1 Dirac system[END_REF][START_REF] Wang | An efficient adaptive mesh redistribution method for a nonlinear Dirac equation[END_REF][START_REF] Pecher | Local well-posedness for the nonlinear Dirac equation in two space dimensions, Commun[END_REF] as a reduced mathematical model. It was originally derived in [START_REF] Soler | Classical, stable, nonlinear spinor field with positive rest energy[END_REF] as a four-component equations in 3D for describing the spinor field with nonlinear coupling. A more complicated cubic nonlinearity could be considered as in [START_REF] Thirring | A soluble relativistic field theory[END_REF][START_REF] Xu | Numerical methods for nonlinear Dirac equation[END_REF]. We refer to [START_REF] Bechouche | Semi-nonrelativistic limits of the Dirac equation with external time-dependent electromagnetic field[END_REF][START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF] for the nondimensionalization and the model reduction to the two-component problem (1.1).

The model (1.1) conserves the mass M (t) := R |Φ ε | 2 dx and the energy E(t), provided that the two external potentials V e ,V m are independent of time,

E(t) := R 1 ε 2 (Φ ε ,βΦ ε ) - i ε (Φ ε ,α∂ x Φ ε ) + V e |Φ ε | 2 + V m (Φ ε ,αΦ ε ) + λ(βΦ ε ,Φ ε )|Φ ε | 2 dx.
With a fixed ε > 0, the well-posedness of the Cauchy problem (1.1) has been studied and established. We refer to [START_REF] Bournaveas | Theory and numerical approximations for a nonlinear 1 + 1 Dirac system[END_REF][START_REF] Huh | Global strong solution to the Thirring model in critical space[END_REF][START_REF] Najman | The nonrelativistic limit of the nonlinear Dirac equation[END_REF][START_REF] Pecher | Local well-posedness for the nonlinear Dirac equation in two space dimensions, Commun[END_REF][START_REF] Th | Solutions of nonlinear Dirc equations[END_REF] and the references therein for detailed analytical results. Various numerical methods including finite difference time domain (FDTD) methods and operator splitting methods have also been considered in [START_REF] Shao | Interaction of solitary waves with a phase shift in a nonlinear Dirac model[END_REF][START_REF] Xu | Numerical methods for nonlinear Dirac equation[END_REF][START_REF] Wang | An efficient adaptive mesh redistribution method for a nonlinear Dirac equation[END_REF][START_REF] Bournaveas | Theory and numerical approximations for a nonlinear 1 + 1 Dirac system[END_REF][START_REF] Cao | Time-Compact Scheme for the One-Dimensional Dirac Equation[END_REF] for solving the nonlinear Dirac equation in the classical regime, i.e. ε ≈ 1. When ε → 0, which corresponds to the speed of light going to infinity and is known as the nonrelativistic limit, the nonlinear Dirac equation (1.1) has been shown to converge in the energy space to the nonlinear Schrödinger equations [START_REF] Bechouche | Semi-nonrelativistic limits of the Dirac equation with external time-dependent electromagnetic field[END_REF][START_REF] Najman | The nonrelativistic limit of the nonlinear Dirac equation[END_REF]. Both the analytical results and the numerical studies [START_REF] Bechouche | Semi-nonrelativistic limits of the Dirac equation with external time-dependent electromagnetic field[END_REF][START_REF] Najman | The nonrelativistic limit of the nonlinear Dirac equation[END_REF][START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | A uniformly accurate (UA) multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime[END_REF] show that the solution Φ ε of (1.1) propagates waves with wavelength O(ε 2 ) in time t and wavelength O(1) in space x, as ε 1. The small temporal wavelength makes the solution highly oscillatory in time, and as a consequence it causes severe numerical burden for classical numerical discretizations. Results in [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | Error estimates of numerical methods for the nonlinear Dirac equation in the nonrelativistic limit regime[END_REF] show that, in order to capture the correct solution in the limit regime, one needs to use time step ∆t = O(ε 2 ) and mesh size ∆x = O( √ ε) for FDTD, and ∆t = O(ε 2 ), ∆x = O(1) for exponential integrators or time splitting spectral methods. Thus designing numerical methods that allow the use of step size independent of ε becomes a very needed and challenging issue. Recently, in study of the nonlinear Klein-Gordon equation in nonrelativistic limit regime, schemes with uniformly accurate (UA) properties have been proposed in [START_REF] Bao | A uniformly accurate multiscale time integrator pseudospectral method for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | Uniformly accurate multiscale time integrators for highly oscillatory second order differential equations[END_REF][START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF]. The UA schemes can use ∆t = O(1), ∆x = O(1) for corrected approximations when ε 1 and are significantly better than classical methods in the limit regime and asymptotic preserving (AP) methods in the intermediate regime, as shown in [START_REF] Bao | A numerical study and comparison of numerical methods for the nonlinear Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF]. Among the UA schemes, there are two approaches so far. One is based on the multiscale expansion of the solution [START_REF] Bao | A uniformly accurate multiscale time integrator pseudospectral method for the Klein-Gordon equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | Uniformly accurate multiscale time integrators for highly oscillatory second order differential equations[END_REF], and the other uses the two-scale formulation approach [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]. The first approach strongly relies on the pre-knowledge of the expansion and the polynomial type nonlinearity. In fact, a multiscale approach has already been considered for the linear Dirac in [START_REF] Bao | A uniformly accurate (UA) multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime[END_REF], leading to a first order UA scheme. However, this proposed scheme and its corresponding numerical analysis become rather complicated in the nonlinear case. In contrast, the two-scale formulation approach provides a general strategy to design UA schemes for a class of oscillatory problems [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]. Nevertheless, the nonlinear Dirac equation (1.1) which contains two scales ε and ε 2 , does not belong to the type of problems in [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF], which means that the strategy in [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF] cannot be directly applied. In fact, problems involving two scales are usual in collisional kinetic equations in the diffusion limit, but the presence of these two scales (ε and ε 2 ) also appears in the context of oscillatory kinetic equations, as mentioned in [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Frénod | Long time behaviour of an exponential integrator for a Vlasov-Poisson system with strong magnetic field[END_REF]. This extra scaling in the equations causes more difficulties for achieving UA property and a suitable two-scale formulation has to be found.

In this work, we are going to propose a two-scale formulation and design UA schemes for solving the nonlinear Dirac equation in the nonrelativistic limit regime. A suitable formulation is presented on which we construct a second order in time UA scheme.

Extensive numerical experiments are done to show the UA property. Extension of this strategy to the case of oscillatory kinetic equations with strong magnetic field is the subject of a forthcoming paper.

The rest of the paper is organized as follows. In Section 2, we introduce the twoscale formulation and construct the suitable initial data for this augmented problem. In Section 3, we present the numerical schemes based on this two-scale formulation. Numerical results are reported in Section 4 and some conclusions are finally drawn in Section 5.

Two scale formulation

To filter out the main oscillation in (1.1), let us introduce the filtered data

u ε = e it/ε 2 0 0 e -it/ε 2 Φ ε . (2.1)
It solves the equation

∂ t u ε = - 1 ε A(t/ε 2 )∂ x u ε + F (t,t/ε 2 ,u ε ), t > 0, x ∈ R, (2.2)
with

A(τ ) = 0 e 2iτ e -2iτ 0 , and
F (t,τ,u ε ) = -i[V e (t) + V m (t)A(τ )]u ε -iλ(βu ε ,u ε )βu ε . (2.3)
In this paper, we assume that V e and V m are C 1 functions of time. For future convenience, we denote the nonlinearity without the magnetic potential by

F e (t,u ε ) := -i[V e (t)u ε + λ(βu ε ,u ε )βu ε ].
Now, following [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF], we separate the fast time variable τ = t/ε 2 from the slow time variable t in the solution u ε (t,x) of (2.2) and consider the two-scale formulation for the augmented unknown

U ε = U ε (t,τ,x)∈ C 2 : ∂ t U ε + 1 ε 2 ∂ τ U ε = - 1 ε A(τ )∂ x U ε + F (t,τ,U ε ), t > 0, τ ∈ T, x ∈ R. (2.4)
Here U ε (t,τ,x) is 2π-periodic in τ and T := R/(2πZ) denotes the torus. With initial data U ε (0,τ,x) satisfying

U ε (0,0,x) = u ε (0,x) = Φ 0 (x), x ∈ R,
the solution of the two-scale problem (2.4) enables to recover the solution of the filtered equation (2.2) by setting

U ε t, t ε 2 ,x = u ε (t,x), t ≥ 0, x ∈ R.
In the two-scale formulation (2.4), τ is considered as an additional independent variable, which is periodic. The advantage one can get from this formulation is that now the initial data U ε (0,τ,x) is only prescribed at a single point τ = 0, so there is some freedom to choose the initial data. By choosing a suitable initial data, we will succeed to bound the time derivatives of U ε uniformly with respect to ε, which allows to design uniformly accurate numerical schemes.

A toy model

In order to find out a suitable initial data for problems like (2.4), we propose to analyze the following simple toy model:

∂ t u + 1 ε 2 ∂ τ u + i a(τ ) ε u = 0, with u(0,τ ) = u in (τ ). (2.5)
Here, u(t,τ ) is a scalar unknown, t ≥ 0 and τ ∈ T. Moreover, we assume (as for the augmented problem (2.4)) that only u in (0) = u 0 is prescribed, that a is a real-valued smooth periodic function and that we have the property 2π 0 a(τ )dτ = 0.

Our aim is to derive expressions for u in (τ ) in order to ensure conditions (i) u in (0) = u 0 , and (ii) that time-derivative of u, up to some order p ≥ 1, are uniformly bounded with respect to ε, a property which enables to construct easily uniformly accurate numerical schemes for (2.5).

The exact solution of (2.5) is given by

u(t,τ ) = e -iεb(τ ) e iεb(τ -t/ε 2 ) u in (τ -t/ε 2 ), where b(τ ) = τ 0 a(s)ds. Note that b is periodic, since the average of a vanishes. It is clear that ∂ k u dt k (t,τ ) = e -iεb(τ ) (-1) k ε 2k d k ds k e iεb(s) u in (s) s=τ -t/ε 2 .
Therefore, in order to get uniformly bounded derivatives for k = 0,...,p, we must have

d p ds p e iεb(s) u in (s) = O(ε 2p ),
which is satisfied as soon as

u in (τ ) = Q p-1 (τ )e -iεb(τ ) + O(ε 2p ).
Here Q p-1 is any polynomial of degree ≤ p -1. Due to the smoothness and periodicity of u in , we have necessarily that Q is a constant polynomial, thus

u in (τ ) = Ce -iεb(τ ) + O(ε 2p ).
Therefore, the following initial data is suitable to ensure (i) and (ii):

u in (τ ) = u 0 + u 0 2p-1 k=1 (-ib(τ )) k k! ε k .
In particular, one observes that one needs to choose the initial data as an expansion in powers ε up to the order ε 2p-1 in order to ensure the boundedness of the time derivatives of u up to the order p. This crucial property will guide our analysis in the next subsection. Of course, for (2.4), we do not have an exact solution, but we will obtain iteratively the expansion in powers of ε of the suitable initial data by using Chapman-Enskog techniques.

Suitable initial data for the augmented problem (2.4)

In order to define the initial data U ε (0,τ ) where here and after we omit the space variable x for simplicity, let us perform formally the Chapman-Enskog expansion of U ε [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF]. To this aim, we introduce the operators L and Π for a periodic function h(τ ) : T → C 2 as

Lh := ∂ τ h, Πh := 1 2π 2π 0 h(τ )dτ,
and when Πh = 0, the operator L is invertible with

(L -1 h)(τ ) = (I -Π) τ 0 h(θ)dθ.
We perform the Chapman-Enskog expansion by setting

U ε (t,τ ) = U (t) + h(t,τ ), with U = ΠU ε .
We observe that ΠA = 0 and write the micro-macro formulation of (2.4) as

∂ t U = - 1 ε Π(A(τ )∂ x h) + Π(F (t,τ,U + h)), (2.6)
∂ t h = - 1 ε 2 Lh - 1 ε A(τ )∂ x U - 1 ε (I -Π)(A(τ )∂ x h) + (I -Π)F (t,τ,U + h). (2.7)
Based on this formulation, it appears that, when ε → 0, U ε (t,τ ) converges to U (t) satisfying

∂ t U = C∂ 2 x U + Π(F (t,τ,U )) = F e (t,U ), (2.8)
where the matrix C is given below by (2.12). In particular, when ε is small, the solution Φ ε of the nonlinear Dirac equation (1.1) will formally be close to Φ defined by Φ = e -it/ε 2 0 0 e it/ε 2 U.

The above limit model (2.8) is a system of coupled nonlinear Schrödinger equations.

When there is no magnetic potential, i.e. V m = 0, the convergence of the model (2.2) to (2.8) as ε → 0 has been proved rigorously in some energy space in [START_REF] Najman | The nonrelativistic limit of the nonlinear Dirac equation[END_REF].

Uniform boundedness of time derivatives up to order 2

In this subsection, we derive the expression of U ε (0,τ ) such that the time derivatives of U ε up to order 2 are uniformly bounded. Following the previous subsection on the toy model, we need an expansion up to the order ε 3 . In the following formal calculations, we thus assume that derivatives of h until order 2 are bounded.

By applying the inverse of L to (2.7), we get

h = -εL -1 A∂ x U -εL -1 (I -Π)A∂ x h + ε 2 L -1 (I -Π)F (t,τ,U + h) -ε 2 L -1 ∂ t h. (2.9)
Then we further have

∂ x h = -εL -1 A∂ 2 x U -εL -1 (I -Π)A∂ 2 x h + ε 2 L -1 (I -Π)∂ x F (t,τ,U + h) -ε 2 L -1 ∂ tx h, (2.10)
∂ t h = -εL -1 A∂ tx U -εL -1 (I -Π)A∂ tx h + ε 2 L -1 (I -Π) ∂ t F (t,τ,U + h) + ∂ u F (t,τ,U + h)(∂ t U + ∂ t h) -ε 2 L -1 ∂ 2 t h. (2.11)
Hence, we can expand the unknowns in powers of ε and find

h = -εL -1 A∂ x U + ε 2 L -1 (I -Π)AL -1 A∂ 2 x U + ε 2 L -1 (I -Π)F (t,τ,U ) + O(ε 3 ).
Therefore,

U ε (0,τ ) =U -εL -1 A∂ x U + ε 2 L -1 (I -Π)AL -1 A∂ 2 x U + ε 2 L -1 (I -Π)F (t,τ,U ) + O(ε 3 ). Direct computations yield B(τ ) := L -1 A(τ ) = - i 2 0 e 2iτ -e -2iτ 0 and C ≡ C(τ ) := A(τ )B(τ ) = i 2 1 0 0 -1 (2.12) so (I -Π)C(τ ) = 0 and U ε (0,τ ) = U (0) -εB(τ )∂ x U (0) + ε 2 L -1 (I -Π)F (0,τ,U ) + O(ε 3 ). Now we use that U ε (0,0) = Φ 0 to get U (0) = Φ 0 + εB(0)∂ x U (0) -ε 2 f 0 (0) + O(ε 3 ) = Φ 0 + εB(0)∂ x Φ 0 + ε 2 B(0) 2 ∂ 2 x Φ 0 -ε 2 f 0 (0) + O(ε 3 ), (2.13)
where

f 0 (τ ) := L -1 (I -Π)F (0,τ,Φ 0 ) = -iV m (0)B(τ )Φ 0 ,
and then

U ε (0,τ ) =Φ 0 -ε(B(τ ) -B(0))∂ x Φ 0 -ε 2 (B(τ ) -B(0))B(0)∂ 2 x Φ 0 + ε 2 (f 0 (τ ) -f 0 (0)) + O(ε 3 ).
We thus get an expression U ε 2 for the suitable initial data up to

O(ε 3 ), i.e. U ε (0,τ ) = U ε 2 (τ ) + O(ε 3
), as follows:

U ε 2 (τ ) :=Φ 0 + iε 2 0 e 2iτ -1 1 -e -2iτ 0 ∂ x Φ 0 + ε 2 4 1 -e 2iτ 0 0 1 -e -2iτ ∂ 2 x Φ 0 - ε 2 2 V m (0) 0 e 2iτ -1 1 -e -2iτ 0 Φ 0 .
For future use, we denote

U 2 (τ ) = U ε 2 (τ ) -ε 2 (f 0 (τ ) -f 0 (0)
), and also define

U ε 1 (τ ) := Φ 0 + iε 2 0 e 2iτ -1 1 -e -2iτ 0 , which gives U ε (0,τ ) = U ε 1 (τ ) + O(ε 2 ).
Now we push the asymptotic expansion to the next order. Inserting (2.10) and (2.11) into (2.9) and letting t = 0, we get

h(0,τ ) = -εB∂ x U -εL -1 (I -Π)A -εB∂ 2 x U + ε 2 L -1 (I -Π)C∂ 3 x U +ε 2 L -1 (I -Π)∂ x F (0,τ,U ) + ε 2 L -1 (I -Π)F (0,τ,U + h) + ε 3 L -1 B∂ tx U + O(ε 4 ) = -εB∂ x U + ε 2 L -1 (I -Π)F (0,τ,U ε 1 ) + ε 3 L -1 B∂ tx U -ε 3 L -1 (I -Π)AL -1 (I -Π)∂ x F (0,τ,U ) + O(ε 4 ).
(2.14)

From (2.6) and (2.10), we have

∂ t U = C∂ 2 x U + ΠF (t,τ,U ) + O(ε). (2.15)
Noting F e (t,U ) = ΠF (t,τ,U ), we find

h(0,τ ) = -εB∂ x U + ε 2 L -1 (I -Π)F (0,τ,U ε 1 ) + ε 3 L -1 B C∂ 3 x U + ∂ x F e (0,U ) -ε 3 L -1 (I -Π)AL -1 (I -Π)∂ x F (0,τ,U ) + O(ε 4 ).
Using (2.13), we then get

h(0,τ ) = -εB(τ ) ∂ x Φ 0 + εB(0)∂ 2 x Φ 0 + ε 2 B 2 (0)∂ 3 x Φ 0 -ε 2 ∂ x f 0 (0) + ε 2 f 1 (τ ) + ε 3 L -1 B C∂ 3 x Φ 0 + ∂ x F e (0,Φ 0 ) -ε 3 L -1 (I -Π)A∂ x f 0 + O(ε 4 ), (2.16)
where

f 1 (τ ) := L -1 (I -Π)F (0,τ,U ε 1 ).
Based on

U ε (0,τ ) = U (0) + h(0,τ ) = Φ 0 + h(0,τ ) -h(0,0),
and remarking that L -1 B = -A/4 and L -1 (I -Π)A∂ x f 0 ≡ 0, we can update the expansion of U (0) as

U (0) =Φ 0 + εB(0)∂ x Φ 0 + ε 2 B 2 (0)∂ 2 x Φ 0 -ε 2 f 1 (0) + ε 3 B 3 (0)∂ 3 x Φ 0 -ε 3 B(0)∂ x f 0 (0) + ε 3 4 A(0) C∂ 3 x Φ 0 + ∂ x F e (0,Φ 0 ) + O(ε 4 ), (2.17)
and then get an expression U ε 3 of the prepared initial data up to O(ε 4 ), i.e. U ε (0,τ ) = U ε 3 (τ ) + O(ε 4 ), as follows:

U ε 3 (τ ) =U 2 (τ ) + ε 2 [f 1 (τ ) -f 1 (0)] -ε 3 [B(τ ) -B(0)]B 2 (0)∂ 3 x Φ 0 (2.18) - ε 3 4 (A(τ ) -A(0)) C∂ 3 x Φ 0 + ∂ x F e (0,Φ 0 ) + ε 3 (B(τ ) -B(0))∂ x f 0 (0) =U 2 (τ ) + ε 2 [f 1 (τ ) -f 1 (0)] + iε 3 4 0 e 2iτ -1 1 -e -2iτ 0 ∂ 3 x Φ 0 + iε 3 4 1 -e 2iτ 0 0 1 -e -2iτ ∂ x (V m (0)Φ 0 ) + ε 3 4 0 1 -e 2iτ 1 -e -2iτ 0 ∂ x F e (0,Φ 0 ). (2.19)
2.2.2. Uniform boundedness of time derivatives up to order 3 Repeating the procedure above, we are going to get the next order expansion. Here we assume that derivatives of h until order 3 are bounded. From (2.11), we have

∂ t h = -εB∂ tx U + ε 2 L -1 (I -Π)[∂ t F (t,τ,U ) + ∂ u F (t,τ,U )∂ t U ] + O(ε 3 ).
Plugging the above expansion and (2.14) into (2.9) and letting t = 0, we get

h(0,τ ) = -εB∂ x U + εL -1 (I -Π)A εB∂ 2 x U -ε 2 f 1 (τ ) -ε 3 L -1 B∂ txx U + ε 2 L -1 (I -Π)F (0,τ,U ε 2 ) + ε 3 L -1 B∂ tx U -ε 4 f t (τ ) + O(ε 5 ) = -εB∂ x U + ε 2 L -1 (I -Π)F (0,τ,U ε 2 ) -ε 3 L -1 (I -Π)Af 1 - ε 3 4 A∂ tx U -ε 4 f t (τ ) + O(ε 5 ), (2.20)
where

f t (τ ) := L -2 (I -Π) ∂ t F (0,τ,Φ 0 ) + ∂ u F (0,τ,Φ 0 )(C∂ 2 x Φ 0 + F e (0,Φ 0 )) = i 4 A(τ ) ∂ t V m (0)Φ 0 + V m (0)(C∂ 2 x Φ 0 + F e (0,Φ 0 )) .
From (2.6), we find

∂ t U (0) = C∂ 2 x U + ΠF (0,τ,U ε 1 ) -εΠA∂ x f 1 + O(ε 2 ).
Then (2.20) becomes

h(0,τ ) = -εB(τ )∂ x U + ε 2 L -1 (I -Π)F (0,τ,U ε 2 ) -ε 3 L -1 (I -Π)Af 1 - ε 3 4 A(τ ) C∂ 3 x U + Π∂ x F (0,τ,U ε 1 ) -εΠA∂ 2 x f 1 -ε 4 f t (τ ) + O(ε 5 ). (2.21)
Combining (2.21) with (2.17) and noting A(0)C = 4B 3 (0), we get

h(0,τ ) = -εB(τ ) ∂ x Φ 0 + εB(0)∂ 2 x Φ 0 + ε 2 B 2 (0)∂ 3 x Φ 0 -ε 2 ∂ x f 1 (0) + 2ε 3 B 3 (0)∂ 4 x Φ 0 -ε 3 B(0)∂ 2 x f 0 (0) + ε 3 4 A(0)∂ 2 x F e (0,Φ 0 ) + ε 2 L -1 (I -Π)F (0,τ,U ε 2 ) -ε 3 L -1 (I -Π)Af 1 - ε 3 4 A(τ )C ∂ 3 x Φ 0 + εB(0)∂ 4 x Φ 0 - ε 3 4 A(τ ) Π∂ x F (0,τ,U ε 1 ) -εΠA∂ 2 x f 1 -ε 4 f t (τ ) + O(ε 5 ). (2.22)
Then we have

U (0) =Φ 0 + εB(0)∂ x Φ 0 + ε 2 B 2 (0)∂ 2 x Φ 0 -ε 2 f ε 2 (0) + 2ε 3 B 3 (0)∂ 3 x Φ 0 -ε 3 B(0)∂ x f 1 (0) + ε 3 4 A(0)Π∂ x F (0,τ,U ε 1 ) + ε 3 g ε (0) + 3ε 4 B 4 (0)∂ 4 x Φ 0 -ε 4 B 2 (0)∂ 2 x f 0 (0) + ε 4 4 B(0)A(0)∂ 2 x F e (0,Φ 0 ) - ε 4 4 A(0)ΠA∂ 2 x f 1 + ε 4 f t (0) + O(ε 5 ), (2.23)
where

f 2 (τ ) := L -1 (I -Π)F (0,τ,U ε 2 ), g 1 (τ ) := L -1 (I -Π)Af 1 .
Thus, we get the prepared initial data U ε (0,τ ) = U ε 4 (τ ) + O(ε 5 ), with

U ε 4 (τ ) = U 2 (τ ) + ε 2 [f 2 (τ ) -f 2 (0)] + iε 3 4 0 e 2iτ -1 1 -e -2iτ 0 ∂ 3 x Φ 0 (2.24) + ε 3 4 0 1 -e 2iτ 1 -e -2iτ 0 Π∂ x F (0,τ,U ε 1 ) -ε 3 [g 1 (τ ) -g 1 (0)] - iε 3 2 0 e 2iτ -1 1 -e -2iτ 0 ∂ x f 1 (0) - 3ε 4 16 e 2iτ -1 0 0 e -2iτ -1 ∂ 4 x Φ 0 - ε 4 8 0 e 2iτ -1 1 -e -2iτ 0 ∂ 2 x (V m (0)Φ 0 ) + ε 4 4 0 e 2iτ -1 e -2iτ -1 0 ΠA∂ 2 x f 1 - iε 4 8 1 -e 2iτ 0 0 e -2iτ -1 ∂ 2 x F e (0,Φ 0 ) -ε 4 (f t (τ ) -f t (0)).
To get the next order expansion, we recall that we assume

∂ t h,∂ 2 t h,∂ 3 t h = O(1) as ε → 0, which indicates that ∂ t h = -εB∂ tx U + ε 2 L -1 (I -Π) d dt F (t,τ,U (t) + h(t)) + ε 3 L -1 B∂ ttx U -ε 3 L -1 (I -Π)AL -1 (I -Π) d dt ∂ x F (t,τ,U (t)) + O(ε 4 ). Using d dt F (t,τ,U (t) + h(t)) = ∂ t F (t,τ,U + h) + ∂ u F (t,τ,U + h)(∂ t U + ∂ t h) = ∂ t F (t,τ,U -εB∂ x U ) + ∂ u F (t,τ,U -εB∂ x U )(∂ t U -εB∂ tx U ) + O(ε 2 )
and denoting

H 0 (τ ) :=L -1 (I -Π)[∂ t F (0,τ,Φ 0 ) + ∂ u F (0,τ,Φ 0 )C∂ 2 x Φ 0 + ∂ u F (0,τ,Φ 0 )F e (0,Φ 0 )], H 1 (τ ) :=∂ t F (0,τ,U ε 1 ) + ∂ u F (0,τ,U ε 1 )C(∂ 2 x Φ 0 + εB(0)∂ 3 x Φ 0 ) + ∂ u F (0,τ,U ε 1 )[ΠF (0,τ,U ε 1 ) -εΠA∂ x f 1 ] -ε∂ u F (0,τ,U ε 1 )B(τ ) C∂ 3 x Φ 0 + ∂ x F e (0,Φ 0 ) , we get ∂ t h(0,τ ) = -εB(τ )∂ tx U + ε 2 L -1 (I -Π)H 1 + ε 3 L -1 B∂ ttx U -ε 3 L -1 (I -Π)A∂ x H 0 + O(ε 4 ). ( 2

.25)

A detailed computation gives

H 0 (τ ) = -iB(τ ) ∂ t V m (0)Φ 0 + V m (0)(C∂ 2 x Φ 0 + F e (0,Φ 0 )) .
Plugging (2.16) into (2.6) at t = 0, we have

∂ t U (0) =C∂ 2 x Φ 0 + ΠF (0,τ,U ε 2 ) + εCB(0)∂ 3 x Φ 0 + ε 2 2 C∂ 4 x Φ 0 -εΠA∂ x f 1 -ε 2 C∂ 2 x f 0 (0) + ε 2 4 ∂ 2
x F e (0,Φ 0 ) + O(ε 3 ), and

∂ 2 t U (0) =C 2 ∂ 4 x Φ 0 + C∂ 2 x F e (0,Φ 0 ) + Z e + O(ε),
where

Z e = ∂ u F e (0,Φ 0 ) C∂ 2 x Φ 0 + F e (0,Φ 0 ) + ∂ t F e (0,Φ 0 ).
Combining the above two identities with (2.25), and inserting them together with (2.23) and (2.22) into (2.9), we finally obtain the prepared initial data U ε (0,τ ) = U ε 5 (τ ) + O(ε 6 ), with

U ε 5 (τ ) = U 2 (τ ) + ε 2 (f 3 (τ ) -f 3 (0)) + iε 3 4 0 e 2iτ -1 1 -e -2iτ 0 ∂ 3 x Φ 0 (2.26) + ε 3 4 0 1 -e 2iτ 1 -e -2iτ 0 Π∂ x F (0,τ,U ε 2 ) -ε 3 (g 2 (τ ) -g 2 (0)) - iε 3 2 0 e 2iτ -1 1 -e -2iτ 0 ∂ x f 2 (0) + ε 4 4 e 2iτ -1 0 0 e -2iτ -1 ∂ 2 x f 1 (0) - iε 4 8 1 -e 2iτ 0 0 e -2iτ -1 Π∂ 2 x F (0,τ,U ε 1 ) - 3ε 4 16 e 2iτ -1 0 0 e -2iτ -1 ∂ 4 x Φ 0 + iε 4 2 0 e 2iτ -1 1 -e -2iτ 0 ∂ x g 1 (0) + ε 4 (v(τ ) -v(0)) -ε 4 (w 1 (τ ) -w 1 (0)) + ε 4 4 0 e 2iτ -1 e -2iτ -1 0 ΠA∂ 2 x f 1 + 3iε 5 16 0 e 2iτ -1 1 -e -2iτ 0 ∂ 5 x Φ 0 + ε 5 (w 0 (τ ) -w 0 (0)) - 3ε 5 16 0 e 2iτ -1 e -2iτ -1 0 ∂ 3 x F e (0,Φ 0 ) + iε 5 8
1 -e 2iτ 0 0 e -2iτ -1 ΠA∂ 3

x f 1 -

ε 5 8 e 2iτ -1 0 0 1 -e -2iτ ∂ x Z m - iε 5 8 e 2iτ -1 0 0 e -2iτ -1 ∂ 3 x (V m (0)Φ 0 ) - iε 5 8 0 e 2iτ -1 1 -e -2iτ 0 ∂ x Z e where f 3 (τ ) := L -1 (I -Π)F (0,U ε 3 ), g 2 (τ ) = L -1 (I -Π)A∂ x f 2 , v(τ ) := L -1 (I -Π)AL -1 (I -Π)A∂ x f 1 , Z m = ∂ t V m (0)Φ 0 + V m (0)(C∂ 2 x Φ 0 + F e (0,Φ 0 )), w 0 (τ ) = L -2 (I -Π)A∂ x H 0 , w 1 (τ ) := L -2 (I -Π)H 1 .
The augmented problem (2.4), endowed with one of the initial data U ε 2k-1 (k ∈ {1,2,3}), provides a solution which has uniformly bounded time derivatives up to the order k, respectively. The proof of this property, which will not be done here, can be derived in an analogous way as in [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]. The detailed proof would be lengthy and tedious, so we omit it for brevity.

Numerical methods

From the analytical results in [START_REF] Bechouche | Semi-nonrelativistic limits of the Dirac equation with external time-dependent electromagnetic field[END_REF][START_REF] Najman | The nonrelativistic limit of the nonlinear Dirac equation[END_REF], the solution of the Dirac equation (1.1) decays very fast at infinity for localized initial data. Thus, to give the numerical discretization, we truncate the Dirac equation (1.1) in the whole space to a bounded interval (a,b) which is large enough, and impose periodic boundary conditions as done in [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | A uniformly accurate (UA) multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Xu | Numerical methods for nonlinear Dirac equation[END_REF][START_REF] Wang | An efficient adaptive mesh redistribution method for a nonlinear Dirac equation[END_REF], i.e.

i∂ t Φ ε = - i ε α∂ x Φ ε + 1 ε 2 βΦ ε + [V e (t) + V m (t)α]Φ ε + λ(βΦ ε ,Φ ε )βΦ ε , x ∈ (a,b), Φ ε (t,a) = Φ ε (t,b), ∂ x Φ ε (t,a) = ∂ x Φ ε (t,b), t ≥ 0; Φ ε (0,x) = Φ 0 (x).
(3.1)

Uniformly accurate schemes

Based on the two-scale formulation, we give the numerical discretization to solve the nonlinear Dirac equation. Due to the truncation (3.1), the corresponding two-scale problem (2.4) reads

∂ t U ε + 1 ε 2 ∂ τ U ε = - 1 ε A(τ )∂ x U ε + F (t,τ,U ε ), t > 0, τ ∈ T, x ∈ (a,b), (3.2a) U ε (t,τ,x) = U ε (t,τ + 2π,x), t ≥ 0, τ ∈ T, x ∈ (a,b), (3.2b) U ε (t,τ,a) = U ε (t,τ,b), t ≥ 0, τ ∈ T, (3.2c)
with the constructed initial data

U ε (0,τ,x) = U ε,0 (τ,x), τ ∈ T, x ∈ (a,b). (3.3)
Let ∆t > 0 be the time step, and denote t n = n∆t for n = 0,1,... Let us choose mesh sizes ∆x = (b -a)/N and ∆τ = 2π/N τ with N and N τ two positive even integers. We also denote the grid points in x and τ by x j := a + j∆x, j = 0,1,...,N ; τ j = j∆τ, j = 0,1,...,N τ .

For first order time discretization, inspired by the scheme in [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], we propose the following semi-implicit Euler method:

U n+1 (τ,x) -U n (τ,x) ∆t + 1 ε 2 ∂ τ U n+1 (τ,x) = - 1 ε A(τ )∂ x U n+1 (τ,x) + F n (τ,x), (3.4)
where we denote U n (τ,x) ≈ U ε (t n ,τ,x) and F n (τ,x) ≈ F (t n ,τ,U ε (t n ,τ,x)), for n = 0,1,...

3.1.1.

Uniform accuracy of the semi-implicit scheme (3.4) We recall that the initial data U ε 3 (τ ) given by (2.19) has been obtained by Chapman-Enskog expansions, in order to ensure the uniform boundedness of U ε , ∂ t U ε and ∂ 2 t U ε with respect to ε. Although we omit the rigorous proof of this property, we assume that this holds and use it to prove the uniform accuracy of the semi-discrete scheme (3.4) initialized with U 0 = U ε 3 . In fact, we have the following result. Theorem 3.1. Assume that the solution U ε (t,τ,x) of (2.4) has uniformly bounded (with respect to ε) derivatives for t ∈ [0,T ], up to order two. Let U n (τ,x) be the solution of the numerical scheme (3.4), initialized with U 0 = U ε,0 , and ∆t > 0. Then, for n∆t ≤ T , we have

U ε (n∆t) -U n L ∞ τ L 2
x ≤ C∆t, where C is a constant independent of n, ∆t and ε.

Proof. Taking the Fourier transform of (3.4) with respect to the variable x, we get

U n+1 -U n ∆t + 1 ε 2 ∂ τ U n+1 = - 1 ε A(τ )iµ U n+1 + F n ,
where µ is the Fourier variable associated to x. In terms of the following quantity

Q(V ) = V + ∆t ε A(τ )iµV + ∆t ε 2 ∂ τ V, the scheme reads Q U n+1 = U n + ∆t F n . (3.5)
Step 1: invertibility property of Q. We claim that Q is invertible on L ∞ τ L 2 µ and that

Q -1 (W ) L ∞ τ L 2 µ ≤ W L ∞ τ L 2 µ . (3.6)
Indeed, consider the equation Q(V ) = W and set

Z(τ ) := e -iτ 0 0 e iτ V.
We get the equation

Z + iΩZ + ∆t ε 2 ∂ τ Z = e -iτ 0 0 e iτ W, with Ω = ∆t/ε 2 µ∆t/ε µ∆t/ε -∆t/ε 2 .
Then,

d dτ e ε 2 ∆t (I+iΩ)τ Z = ∆t ε 2 e ε 2 ∆t (I+iΩ)τ e -iτ 0 0 e iτ W.
We integrate this expression between τ and τ + 2π and use the periodicity of Z(τ ) to obtain e

τ +2π ∆t ε 2 e i τ +2π ∆t ε 2 Ω -e τ ∆t ε 2 e i τ ∆t ε 2 Ω Z(τ ) = ∆t ε 2 τ +2π τ e s ∆t ε 2 e i s ∆t ε 2 Ω e -is 0 0 e is W (s)ds.
Using that Ω is real-valued and symmetric, we get, with the notation

Z = |Z 1 | 2 + |Z 2 | 2 for all Z ∈ C 2 , that e τ +2π ∆t ε 2 e i τ +2π ∆t ε 2 Ω Z(τ ) -e τ ∆t ε 2 e i τ ∆t ε 2 Ω Z(τ ) ≥ e τ +2π ∆t ε 2 e i τ +2π ∆t ε 2 Ω Z(τ ) -e τ ∆t ε 2 e i τ ∆t ε 2 Ω Z(τ ) = e τ +2π ∆t ε 2 -e τ ∆t ε 2 Z(τ ) ,
and that

∆t ε 2 τ +2π τ e s ∆t ε 2 e i s ∆t ε 2 Ω e -is 0 0 e is W (s)ds L 2 µ ≤ ∆t ε 2 τ +2π τ e s ∆t ε 2 e -is 0 0 e is W (s) L 2 µ ≤ e τ +2π ∆t ε 2 -e τ ∆t ε 2 W L ∞ τ L 2 µ .
These two inequalities give Z(τ

) L 2 µ ≤ W L ∞ τ L 2
µ , which proves the claim.

Step 2: Taylor expansion of the exact solution. The exact solution of (2.4) satisfies

U ε (t n ) -U ε (t n+1 ) + ∆t∂ t U ε (t n+1 ) = R n with R n = t n+1 t n (s -t n )∂ 2 t U ε (s)ds which yields Q( U ε (t n+1 )) = U ε (t n ) + ∆t F (t n+1 ,τ, U ε (t n+1 )) -R n .
Here we have denoted G( U ) = G(U ).

We now introduce the error E n = U ε (t n ) -U n and subtracting (3.5) from this equation, we get

Q(E n+1 ) = E n + ∆t F (t n+1 ,τ, U ε (t n+1 )) -F (t n ,τ, U n ) -R n .
(3.7)

Step 3: error estimate. We observe from the expression (2.3) of F and the assumptions on V e and V m that F is locally Lipschitz continuous with respect to t and

U ε . Let us fix M > 0 such that U ε L ∞ t L ∞ τ L 2 x ≤ M for t ∈ [0,T ].
Let N be the number of time discretization points and ∆t = T /N . Let n 0 be the largest integer n ≤ N such that

U n L ∞ τ L 2
x ≤ 2M . Our aim is to prove that n 0 = N and to estimate the error. Let us proceed by contradiction and assume that n 0 ≤ N -1. Then we have

F (t n+1 ,τ, U ε (t n+1 )) -F (t n ,τ, U n ) L ∞ τ L 2 µ ≤ C ∆t + E n L ∞ τ L 2 µ ,
where C, here and after, is a generic constant which only depends on M and T . Applying (3.6) to (3.7), we deduce that

E n+1 L ∞ τ L 2 µ ≤ (1 + C∆t) E n L ∞ τ L 2 µ + C∆t 2 + R n L ∞ τ L 2
µ . Now, in order to estimate R n , we apply the assumption made in Theorem 3.1 saying that ∂ 2 t U is uniformly bounded with respect to ε.

Thus R n L ∞ τ L 2 µ ≤ C∆t 2 . Therefore E n+1 L ∞ τ L 2 µ ≤ (1 + C∆t) E n L ∞ τ L 2 µ + C∆t 2 , for all n ≤ n 0 , which implies E n+1 L ∞ τ L 2 µ ≤ C∆t.
In particular, if n 0 ∆t < T and if we choose C∆t < M then we deduce from this estimate that

U n0+1 L ∞ τ L 2 µ ≤ U ε (t n0+1 ) L ∞ τ L 2 µ + E n0+1 L ∞ τ L 2 µ ≤ 2M. This contradicts the fact that n 0 is the largest integer n ≤ N such that U n L ∞ τ L 2
x ≤ 2M . We conclude that n 0 = N . It is also clear from the above estimates that we have

E n L ∞ τ L 2
µ ≤ C∆t, for all n ≤ N , which concludes the proof.

to ∆x and ∆τ are shown in Fig. 4.5. The behaviour of the spatial errors of the UA1 method and of the other two numerical examples are similar, so the results are omitted here for brevity.

The error between the solution of the nonlinear Dirac equation (1.1) and the solution of the limit model (2.8), i.e.

φ ε 1 (t) -φ 1 (t) L ∞ + φ ε 2 (t) -φ 2 (t) L ∞ ,
are shown in Fig. 4.6 at t = 0.5.

Based on the numerical results, we can have the following observation: (i) The UA1 scheme with initial data U ε 3 shows uniformly first order temporal accuracy and the UA2 scheme with U ε 5 shows uniformly second order temporal accuracy in all the three cases. In space discretization, the UA schemes have uniformly spectral accuracy in both x and the artificial τ .

(ii) The convergence rate from the nonlinear Dirac equation to the limit model is of order O(ε), as ε → 0.

Conclusion

We proposed some uniformly accurately (UA) schemes for solving the nonlinear Dirac equation in the nonrelativistic limit regimes. Our approach is based on a suitable two-scale formulation which offers a general strategy for constructing UA schemes for a class of highly oscillatory problems involving two small scales. We derive correct initial data for this augmented formulation, using a Chapman-Enskog expansion. This allows us to construct a UA scheme with second order accuracy in time and spectral accuracy in space. Numerical tests were done to show the UA property. Our approach can also be applied to solve the oscillatory kinetic equations with diffusion scaling. This the subject of a work in progress. results of using U ε 0 (first row); results of using U ε 1 (second row) results of using U ε 2 (third row); results of using U ε 3 (last row).

Fig. 4 . 1 :

 41 Fig.4.1: Temporal error of the UA1 method in Example I with respect to ∆t and ε: results of using U ε 0 (first row); results of using U ε 1 (second row) results of using U ε 2 (third row); results of using U ε 3 (last row).

Fig. 4 . 2 :

 42 Fig.4.2: Temporal error of the UA2 method in Example I with respect to ∆t and ε: results of using U ε 2 (first row); results of using U ε 3 (second row) results of using U ε 4 (third row); results of using U ε 5 (last row).

Fig. 4 . 3 :

 43 Fig.4.3: Temporal error of the UA2 method in Example II with respect to ∆t and ε: results of using U ε 2 (first row); results of using U ε 3 (second row) results of using U ε 4 (third results of using U ε 5 (last row).

Fig. 4 . 4 :

 44 Fig. 4.4: Temporal error of the UA2 method with U ε 5 in Example III with respect to ∆t and ε: results of the nonlinear Dirac equation with magnetic potential case.

Fig. 4 . 5 :

 45 Fig. 4.5: Spatial error of the UA2 method in Example III with respect to N (= (ba)/∆x) and N τ .

Fig. 4 . 6 :

 46 Fig. 4.6: The maximum error in solution between the nonlinear Dirac equation and the limit model (2.8) in Example III.

Acknowledgements. This work was supported by the French ANR project MOONRISE ANR-14-CE23-0007-01. M. Lemou is supported by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03.

3.1.2. The fully discretized scheme In the x-direction, we apply the Fourier transform and obtain for l = -N/2,...,N/2 -1,

where µ l = 2πl b-a and

U n l (τ )e iµ l (x-a) , F n (τ,x) =

F n l (τ )e iµ l (x-a) , with

F n (τ,x j )e iµ l (xj -a) .

In terms of the two components U ε = (u 1 ,u 2 ) T and F = (f 1 ,f 2 ) T , the numerical scheme (3.8) reads

2

Then we discretize the τ -direction as follows

Introducing the Fourier pseudo-differential matrix D τ [START_REF] Shen | Spectral Methods: algorithms, analysis and applications[END_REF]:

le il(τn-τm) , n,m = 1,...,N τ , we derive the full discretization of (3.9) as

2

2

for l = -N 2 ,..., N 2 -1, where here and after e ±2iτ are interpreted as diagonal matrices. Solving (3.10) and choosing (u 0 1 ,u 0 2 ) T = U ε,0 , we get the detailed scheme of the first order method (UA1) for n = 0,1,..., and l = -N/2,...,N/2 -1,

2

where

and where Id denotes the identity matrix. The inverse of the above matrices A τ ,B l τ can be computed numerically once for all, so the UA1 method (3.11) is explicit in practise. Also noticing that B l τ = -B -l τ , we only need to store and compute the inverse of B l τ for l = -N/2,...,-1. The computational cost of the UA1 scheme is O(N N 2 τ log(N )) per time step.

3.1.3. The uniformly accurate second order scheme Again inspired by [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], our second order method starts with a prediction step from the first order method

followed by the correction step

where F n+1/2 (τ,x) := F (t n+1/2 ,τ,U n+1/2 (τ,x)). Then by applying the Fourier pseudospectral discretization similarly as before, we end up with the following detailed second order scheme (UA2). Choosing (u 0 1 ,u 0 2 ) T = U ε,0 , then for n = 0,1,..., and l = -N/2,..., N/2 -1,

2

2

with

where

Similarly, we remark that the inverse of matrices A + τ ,A for l = -N/2,...,-1. The computational cost of the UA2 scheme is also O(N N 2 τ log(N )) per time step. Thanks to the semi-implicit approximations in time, the time step ∆t of UA1 and UA2 methods are free from any CFL-type conditions on ∆τ , ∆x or stability condition on ε. After obtaining U n (τ,x), we take τ = t n /ε 2 and consider the inverse of (2.1) to get the approximation of Φ ε (t n ,x), i.e.

Numerical results

For comparison purpose, we define the initial data

We shall test the error of the proposed UA1 scheme (3.11) and the UA2 scheme (3.12) under those choices of initial data for (3.3) by performing the following three numerical experiments. The computational domain is chosen as (a,b) = (-8,8). The 'exact' solution is obtained numerically by the UA2 scheme with very small step sizes, e.g. ∆t = 10 -6 ,∆x = 1/64,N τ = 64. We shall focus on the time discretization error.

Example I: (nonlinear without magnetic potential) We take the potential and initial data in (1.1) as

, and choose λ = 0.5. We solve the Dirac equation (3.2) till t = 0.5, and show the error of the numerical solution (φ n 1 ,φ n 2 ) with n = t/∆t as

The time discretization errors of the UA1 and UA2 methods under different ε are shown in Fig. 4.1 and Fig. 4.2, respectively, where we take the spatial mesh size of the numerical method small enough, e.g. ∆x = 1/64,N τ = 32.

Example II: (linear with magnetic potential) We use the same setup in Example I but without nonlinearity and with a non-zero magnetic potential, i.e.

The temporal errors of the UA2 method at t = 0.5 under different ε are shown in Fig.

4.3.

Example III: (nonlinear with magnetic potential) We keep both the nonlinearity and magnetic potential, i.e.

V m = (x + 1) 2 1 + x 2 , λ = 0.5.

The time discretization error of the UA2 method at t = 0.5, with initial data U ε 5 , is shown in Fig. 4.4. The corresponding spatial errors of the UA2 method with respect


D Prot
email: dprot@talend.com

O Bellenguez-Morineau

P M|pmtn

Open Qm|pmtn

= p|Cmax Open

Keywords: Scheduling, Precedence constraints, Complexity 1 Introduction

    How the structure of precedence constraints may change the complexity class of scheduling problems D Prot, Odile Bellenguez-Morineau To cite this version: D Prot, Odile Bellenguez-Morineau. How the structure of precedence constraints may change the complexity class of scheduling problems. 2017.

particular directed acyclic graph called intree. The fact that the precedence graph takes a particular form may transform the complexity of the problem. Most of the time, adding precedence constraints makes problem harder, since the empty graph is included in many graph classes. Yet, it may make the problem easier if it is not the case, for example for the class of graphs of bounded width. For this reason, the idea of this survey is to consider the complexity results in scheduling theory according to the structure of precedence constraints.

We assume that the reader is conversant with the theory of N P-completeness, otherwise the book by [START_REF] Garey | Computers and intractability: a guide to the theory of NPcompleteness[END_REF] is a good entry point. We will discuss, in the conclusion of this paper, the parameterized complexity. The reader can refer to the book of [START_REF] Downey | Parameterized complexity[END_REF] if needed. [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF] offer a large set of complexity results for scheduling problems with precedence constraints. [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling: a survey[END_REF] and [START_REF] Lawler | Sequencing and scheduling: Algorithms and complexity[END_REF] propose two surveys of complexity results for scheduling problems, but they are not necessarily focused on precedence constraints. For complexity results with the most classical precedence constraints (i.e., chains, trees, series-parallel and arbitrary precedence constraints), the reader can refer to the book by [START_REF] Brucker | Scheduling algorithms[END_REF], and/or to the websites [START_REF] Dürr | The scheduling zoo[END_REF] and [START_REF] Brucker | Complexity results for scheduling problems[END_REF]. Note that this website has been recently updated in order to include most of the results discussed in this article. [START_REF] Möhring | Computationally tractable classes of ordered sets[END_REF] proposes a very interesting survey dedicated to specific partial ordered sets, and studies their structure. Some applications to scheduling are also presented. Since this survey was conducted, many results arised in scheduling theory for specific precedence graphs, we hence believe that a new survey would be beneficial for the scheduling community. We restrict our field on purpose to complexity results and do not talk about approximations results [START_REF] Williamson | The design of approximation algorithms[END_REF]), despite the fact that many results arised in this field recently in scheduling theory, such as [START_REF] Svensson | Hardness of precedence constrained scheduling on identical machines[END_REF] and [START_REF] Levey | A Lasserre-based (1+ )-approximation for P m|p j = 1, prec|Cmax[END_REF]. We believe that it is important to limit the scope of the survey, in order to be as complete as possible in a given field.

The paper is organized as follows: In Section 2, we introduce all the specific types of precedence constraints that will be studied in this paper, while scheduling notations are recalled in Section 3.

Section 4 is dedicated to single-machine scheduling problems, and Sections 5 and 6 are respectively dedicated to non-preemptive and preemptive parallel machines scheduling problems. Each of these three sections is based on the following structure: we first give the polynomial results, then the N P-hard cases and last the most interesting open problems. Finally in Section 7 we give some concluding remarks.

Special types of precedence constraints

In this section we introduce, for the sake of completeness, the special types of precedence constraints that have already been studied in scheduling theory. Precedence constraints between jobs are easily modelled by a directed acyclic graph and we will use it as long as possible. Nevertheless,precedence relations may also be treated as partially ordered set (poset) in the terminology of order theory, and some definitions at the end of this section are much easier to understand from this point of view.

First, let us recall some basic graph theory definitions. Let G = (X, A) be a directed acyclic graph, where X denotes the set of vertices and A, the set of arcs.

A directed acyclic graph is a collection of chains if each vertex has at most one successor and at most one predecessor. An inforest (resp. outforest) is a directed acyclic graph where each vertex has at most one successor (resp. predecessor). An intree (resp. outtree) is a connected inforest (resp. outforest).

We call forest (resp. tree) a graph that is either an inforest or an outforest (resp. either an intree or an outtree). An opposing forest is a collection of inforests and outforests.

For each vertex x ∈ X, we can compute its level h(x) which corresponds to the longest directed path starting from x in G. The height of a DAG, denoted h(G), corresponds to the number of levels -1 in this graph, as illustrated with Figure 1. Directed acyclic graphs of bounded height correspond to directed acyclic graphs where the height is bounded by a constant.

Definition 1 (Level order graph) A directed acyclic graph is a level order graph if each vertex of a given level l is a predecessor of all the vertices of level l -1 (see Figure 2). [1978].

Definition 2 (sp-graph) A graph consisting of a single vertex is a sp-graph. Given two sp-graphs

G 1 = (X 1 , A 1 ) and G 2 = (X 2 , A 2 ), the graph G = (X 1 ∪ X 2 , A 1 ∪ A 2 ) is a sp-graph (this is called parallel composition). Given two sp-graphs G 1 = (X 1 , A 1 ) and G 2 = (X 2 , A 2 ), the graph G = (X 1 ∪ X 2 , A 1 ∪ A 2 ∪ X 1 × X 2 ) is a sp-graph (this is called series composition).
An example is given in Figure 3. Note that sp-graph can also be defined by the forbidden subgraph of Fig. 4 The forbidden subgraph for sp-graphs

A divide-and-conquer-graph (DC-graph) is a special sp-graph built using symmetries. It can be used to model divide-and-conquer algorithms (for example binary search, merge sort, number multiplication...)

and is formally defined as follows:

Definition 3 (DC-graph) A single vertex is a DC-graph; given two vertices s and t and k DC-graphs

(X 1 , A 1 ), . . . , (X k , A k ), the graph (∪ n k=1 X k ∪ {s} ∪ {t}, ∪ n k=1 (A k ∪ (s × X k ) ∪ (X k × t))) is a DC-graph.
This definition is illustrated with Figure 5. Definition 4 (Interval order graph) The interval order for a collection of intervals on the real line is the partial order corresponding to their left-to-right precedence relation, so that the interval order graph is the Hasse diagram of an interval order. such that for any two vertices x and y, (x, y) ∈ A ⇐⇒ ex ≤ sy.

An example is given in Figure 6. [START_REF] Papadimitriou | Scheduling interval-ordered tasks[END_REF] show that interval order graphs can also be defined by the forbidden induced subgraph presented in Figure 7. Larger classes of graphs were defined by forbidden subgraphs, such as quasi-interval order graphs and over-interval order graphs (respectively in [START_REF] Moukrim | Optimal scheduling on parallel machines for a new order class[END_REF] and [START_REF] Chardon | The coffman-graham algorithm optimally solves uet task systems with overinterval orders[END_REF]), the quasi-interval order graphs being strictly included in over-interval order graphs. The corresponding forbidden subgraphs are drawn in Figures 8 and9. To ease the reading, the following definitions will be given within the order theory paradigm. We will hence talk about a partial order set P = (X, P ) rather than a directed acyclic graph G = (X, A) to describe the precedence graph, and a partial order P corresponds to the precedence constraints.

Definition 5 (Antichain) Given a partial order set P = (X, P ), an antichain is a subset S of X such that any two elements of S are incomparable.

Definition 6 (Width) Given a poset P = (X, P ), the width of a poset is the size of a maximum antichain.

By extension, for a given directed acyclic graph G = (X, A) we define the width of the graph to be the width of the corresponding poset, denoted by w(G).

The Am-order (first introduced in [START_REF] Moukrim | A relation between multiprocessor scheduling and linear programming[END_REF]) contains the over-interval order for any integer m ≥ 2 and is defined in the following way:

Definition 7 (Am-order) Let P = (X, P ) be a poset. For any two antichains A and B of size at most m, let us define the four sets :

max(A, B) = {x ∈ A ∪ B|∃y ∈ A ∪ B, y P x}, min(A, B) = {x ∈ A ∪ B|∃y ∈ A ∪ V, x P y}, max(A, B) = (A ∩ B) ∪ max(A, B), and min(A, B) = (A ∩ B) ∪ min(A, B).
P is an Am order if and only if there do not exist two antichains A and B of size m at most, such

that |max(A, B)| ≥ m + 1 or |min(A, B)| ≥ m + 1.
We call Am-order graph a directed acyclic graph for which the set of arcs corresponds to an Am-order.

Definition 8 (Linear extension) Given a partial order P over a set X, a linear extension of P over X is a total order respecting P .

Definition 9 (Dimension) The dimension of a poset P = (X, P ) is the mimimum number t of linear extensions 1 , . . . , t such that x P y ⇐⇒ ∀l ∈ 1..t, x l y. In other words, if x||y (x and y are incomparable in P ), then there are at least two linear extensions, one with x y and another one with y x.

An interesting point is that series-parallel graphs are strictly included in directed acyclic graphs of dimension 2.

The fractional dimension of a poset is extending the notion of dimension (see [START_REF] Brightwell | Fractional dimension of partial orders[END_REF]).

Definition 10 (Fractional dimension) For any integer k, t(k) denotes the minimum number of linear extensions such that for any two incomparable elements x and y, there are at least k extensions with x y and k with y x. The fractional dimension is the limit {t(k)/k} as k tends to +∞.

The diagram in Figure 10 In this paper, we will use the standard α|β|γ notation introduced in [START_REF] Graham | Optimization and approximation in deterministic sequencing and scheduling: a survey[END_REF], and updated in [START_REF] Brucker | Scheduling algorithms[END_REF]. We define below the different notations used all along the paper.

The α-field is used for the machine environment. α = 1 corresponds to a single machine problem; if α = P m, there is a fixed number m of identical parallel machines. If this number is arbitrary, it is noted

α = P . Similarly, if α = Qm (resp. α = Q),
it corresponds to a fixed (resp. arbitrary) number of uniform parallel machines, i.e., each machine i has a speed s i and the processing time of a job j on machine i is equal to p j /s i .

The field β ⊂ {β 1 , β 2 , β 3 , β 4 } describes the jobs characteristics, the possible entries that we will deal with are the following ones:

β 1 ∈ {pmtn, •} : pmtn means that preemption is allowed, i.e., a job may be interrupted and finished

later. If β 1 = •, preemption is forbidden.
β 2 describes the precedence constraints. If β 2 = •, there is no precedence constraint, whereas β 2 = prec means that the precedence graph is a general directed acyclic graph. This field can take many values according to the structure of the directed acyclic graph, as presented in the previous section.

For sake of completeness, all the acronyms are recalled hereinafter: chains, intree, outtree, opp.f orest It appears there that classical literature often abuses of terms intree and outtree to handle in fact inf orest and outf orest, that can be mixed in opp.f orest for some problems.

(
β 3 ∈ {r j , •} : if β 3 = r j , each job j has a given release date. If β 3 = •, the release date is 0 for each job.

β 4 represents the processing time of a job j. If β 4 = •, there is one processing time p j for each job j. If p j = p, all the jobs have the same processing time p. When p j = 1, we use the acronym UET which stands for Unit Execution Time. In some cases, the processing time may increase or decrease with either the position of the job or the starting time of the job. If the job is in position r on a machine, the processing time will be noted p

[r] j . If the processing time is depending of the starting time t of the job, it will be written p

(t) j .
The γ-field is related to the objective function of the problem. Let C j be the completion time of job j. The makespan is defined by Cmax = max j C j and the total completion time by j C j . It is clear that, in general, makespan and total completion time are not equivalent. Nevertheless, for some problems, we can show that there exists an ideal schedule, which both minimizes makespan Cmax and total completion time C j . Due-date related objectives are also studied; if d j is the due date of job j, then the lateness of job j is defined by

L j = C j -d j .
The tardiness is T j = max(0, C j -d j ) and the unit penalty U j is equal to one if C j > d j and to zero otherwise. We then can define the maximum lateness Lmax = max L j , the total tardiness T j , the total number of late jobs U j . A weight w j may also exist for each job j, leading to the corresponding objective functions: the total weighted completion time w j C j , the total weighted tardiness w j T j and the total weighted number of late jobs w j U j . All the functions presented so far are regular functions, i.e., they are non-decreasing with C j .

Single machine problems

For single machine problems, most of the interesting results for this survey are related to the total weighted completion time w j C j . It is mainly due to the fact that 1|chains, p j = 1| U j (see [START_REF] Karel | Complexity results for scheduling chains on a single machine[END_REF]) and 1|chains, p j = 1| T j (see [START_REF] Leung | Minimizing total tardiness on a single machine with precedence constraints[END_REF]) are already N P-hard, while 1|prec, r j |Cmax is solvable in polynomial time. Nevertheless, as we will see in Section 4.3, there

are some interesting open problems for other criteria when considering preemption.

Polynomial cases

Lawler [1978] uses Sidney's theory [START_REF] Sidney | Decomposition algorithms for single-machine sequencing with precedence relations and deferral costs[END_REF]) to derive a polynomial-time algorithm to solve problem 1|sp -graph| w i C i . The most important results are based on the concept of a module in a precedence graph G = (N, A): A non-empty subset M ⊂ N is a module if, for each job j ∈ N -M , exactly one of the three conditions holds: 1. j must precede every job in M , 2. j must follow every job in M , 3. j is not constrained to any job in M . Using this concept leads to a very powerful theorem stating that there exists an optimal sequence consistent with any optimal sequence of any module.

A large improvement on this problem has been done recently, by using order theory and by proving that the problem is a special case of the vertex cover problem. More precisely, in Correa and Schulz [2005] the authors conjecture that 1|prec| w j C j is a special case of vertex cover and prove that, under this conjecture, the problem 1|prec| w j C j is polynomial if the precedence graph is of dimension 2. In [START_REF] Ambühl | Single machine precedence constrained scheduling is a vertex cover problem[END_REF], the authors prove the conjecture and hence the result provided in [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF] is considerably extended (since series-parallel graphs are strictly included in DAGs of dimension 2).

Using the same methodology as in [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF], [START_REF] Wang | Single-machine scheduling with deteriorating jobs under a series-parallel graph constraint[END_REF] extend these results to the case where jobs are deteriorating, i.e., processing times are an increasing function of their starting time, and they show that 1|sp -graph, p

(t) j = p j (1 + at)| w j C j and 1|sp -graph, p (t)
j = p j + α j t|Cmax can be solved in a polynomial time, where p j , a, and α j are positive constants and t is the starting time of the job. The same framework is used in [START_REF] Wang | Single-machine scheduling with precedence constraints and positiondependent processing times[END_REF] for position-dependent processing times, and the authors show that 1|sp -graph, p

[r] j = a j + b j r|Cmax and 1|sp -graph, p

[r] j = a j -b j r|Cmax, where r is the position of the job, a j and b j positive constants, are polynomially solvable. [START_REF] Gordon | Single machine scheduling models with deterioration and learning: handling precedence constraints via priority generation[END_REF] propose a more general framework than the one introduced in Wang et al. [2008], that was first presented in [START_REF] Monma | Sequencing with series-parallel precedence constraints[END_REF]. They extend it to different models with deterioration and learning, that were aimed at minimizing either the total weighted completion time or the makespan, with series-parallel precedence constraints.

Minimum NP-hard cases

The problem 1|prec| w j C j is known to be N P-hard (see [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF]). In order to identify for which precedence graph the problem may be polynomially solvable, we present here the minimal (with respect to precedence constraints) N P-hard cases: the problem is still N P-hard even if the precedence graph is:

of indegree at most 2 (see [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF]). Note that this result also stands with equal weights (i.e., w j = 1)

of bounded height (proof is straightforward by using the transformation proposed in [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF]).

an interval order graph (see [START_REF] Ambühl | On the approximability of single-machine scheduling with precedence constraints[END_REF]). This is a strong difference with parallel machine results for which interval order graphs provide often polynomial algorithms, as we will see in the dedicated section.

of fractional dimension (see Definition 10) greater or equal to 3 (see [START_REF] Ambühl | On the approximability of single-machine scheduling with precedence constraints[END_REF]).

Open problems

The problem 1|prec| w i C i has been widely studied, and the boundary between polynomial and N Phard cases is globally well-defined. However, there are still boundaries to defined, as illustrated by the following examples. First, if the fractional dimension f dim of the precedence graph lies in interval ]2, 3[, the problem is open (it is polynomially solvable if f dim ≤ 2 since the fractional dimension of a DAG is less than or equal to the dimension of this DAG and the problem is solvable in polynomial time if the precedence graph is of dimension 2). This is an interesting question but the gap is rather limited. A wider open question is when we consider the problem with equal weights, i.e., 1|prec| C i . The problem is still N P-hard, even if the precedence graph is of indegree at most 2 (see [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF]). Nevertheless, we may hope that the problem becomes polynomial for precedence graphs with dimension larger than two. It is also possible that the problem is solvable in polynomial time if the precedence graph is an interval order graph, and even for larger classes like quasi-interval order graphs and over-interval order graphs.

For the criteria related to tardiness T j and unit penalty U j , it is surprising to see that preemptive problems with precedence constraints have not been studied yet. More precisely, [START_REF] Tian | An O(n 2 ) algorithm for scheduling equal-length preemptive jobs on a single machine to minimize total tardiness[END_REF] have shown that 1|pmtn, r j , p j = p| T j is solvable in polynomial time, but it is still open whether adding precedence constraints leads this problem to be N P-hard or not; the problem remains also open when considering the more general criterion w j T j . The same outline appears when considering unit penalty: 1|r j ; pmtn; p j = p| w j U j is solvable in polynomial time (with an algorithm in O(n 10 ) by [START_REF] Baptiste | Polynomial time algorithms for minimizing the weighted number of late jobs on a single machine with equal processing times[END_REF], and in O(n 4 ) by Baptiste et al. [2004b]); nevertheless, nothing has been shown when adding precedence constraints to the problem. The first research avenue is to study this problem including Chains as a first step and determine whether the problem is polynomial or not.

Parallel machines without preemption

When considering non-preemptive scheduling problems, whatever the structure of precedence graphs, we will mainly focus on problems with equal processing times, since problems P 2||Cmax and P 2|chains| C j are already NP-hard (see [START_REF] Karel Lenstra | Complexity of machine scheduling problems[END_REF] and [START_REF] Du | Scheduling chain-structured tasks to minimize makespan and mean flow time[END_REF]).

Polynomial cases

Makespan criterion and arbitrary number of machines

Three seminal works on parallel machines with precedence constraints are the approaches of [START_REF] Hu | Parallel sequencing and assembly line problems[END_REF], [START_REF] Papadimitriou | Scheduling interval-ordered tasks[END_REF] and [START_REF] Möhring | Computationally tractable classes of ordered sets[END_REF] in which the authors are respectively dealing with trees, interval order graphs and graphs of bounded width.

In [START_REF] Hu | Parallel sequencing and assembly line problems[END_REF], the author proves that problem P |tree, p j = p|Cmax is polynomially solvable by a list scheduling algorithm where the highest priority is given to the job with the highest level (this strategy is called HLF, for Highest Level First). It is unlikely to find a precedence graph that strictly includes trees for which the problem is solvable in polynomial time, since it was proven in [START_REF] Garey | Scheduling opposing forests[END_REF] that scheduling opposing forest is N P-hard. The hardness of the latter problem is mainly due to the arbitrary number of machines, since it is solvable in polynomial time for any fixed number of machines.

In [START_REF] Papadimitriou | Scheduling interval-ordered tasks[END_REF], the authors prove that problem P |io, p j = p|Cmax is polynomially solvable by a list scheduling algorithm where the highest priority is given to the job with the largest number of successors. This result has been improved twice; first [START_REF] Moukrim | Optimal scheduling on parallel machines for a new order class[END_REF] shows that the same algorithm gives an optimal solution if interval order graphs are replaced by quasi-interval order graphs, that properly contains the former. [START_REF] Chardon | The coffman-graham algorithm optimally solves uet task systems with overinterval orders[END_REF] show that the same result does not stand for over-interval order graphs, but the Coffman-Graham algorithm (see [START_REF] Coffman | Optimal scheduling for two-processor systems[END_REF]) can be applied to solve problem P |oio, p j = p|Cmax to optimality.

In [START_REF] Möhring | Computationally tractable classes of ordered sets[END_REF], the author studies the problem with bounded width, equal processing times and the makespan criterion. He shows that the problem can be solved in polynomial time by using dynamic programming on the digraph of order ideals. [START_REF] Middendorf | Transversal graphs for partially ordered sets: Sequencing, merging and scheduling problems[END_REF] proposed their own approach to solve the problem P |w(G) ≤ k, r j , p j = 1|f for any regular function f . More precisely, their algorithm consists in searching a shortest path in the related transversal graph.

When adding release dates, the problem is already N P-hard for intrees, yet it is solvable in polynomial time for outtrees (see Brucker et al. [1977], and [START_REF] Monma | Linear-time algorithms for scheduling on parallel processors[END_REF] for a linear algorithm for the latter problem). Note that there is a strong relationship between scheduling with release dates and Cmax criterion and scheduling with Lmax, by simply looking at the schedule in the reverse way, and reversing the precedence constraints. Hence problems P |intree, p j = p|Lmax is polynomially solvable and P |outtree, p j = p|Lmax is N P-hard. [START_REF] Kubiak | Optimality of HLF for scheduling divide-and-conquer UET task graphs on identical parallel processors[END_REF] recently opened new perspectives, since they show that P |DC -graph, p j = p|Cmax is solvable in polynomial time. They more precisely proved that the Highest Level First strategy (used in [START_REF] Hu | Parallel sequencing and assembly line problems[END_REF]) solves the problem to optimality when the precedence graph is a divide-and-conquer graph (see Definition 3).

Makespan criterion and fixed number of machines

Let us now focus on the case where the number of machines is fixed. Recall that the problem P m|prec, p j =

p|Cmax is still open (this problem is known as [OPEN8] in the book by [START_REF] Garey | Computers and intractability: a guide to the theory of NPcompleteness[END_REF]) for m ≥ 3, and it was solved for m = 2 in Coffman Jr and [START_REF] Coffman | Optimal scheduling for two-processor systems[END_REF]. For opposing forests, [START_REF] Garey | Scheduling opposing forests[END_REF] propose an optimal polynomial algorithm (of complexity O(n m 2 +2m-5 log n)) that consists in a divide and conquer approach, that uses the HLF algorithm as a subroutine. A new algorithm with complexity O(n 2m-2 log n) has been proposed by [START_REF] Dolev | Profile scheduling of opposing forests and level orders[END_REF], who also show that the problem is polynomially solvable for level order graphs (that are strictly included in series-parallel graphs). [START_REF] Dolev | Scheduling precedence graphs of bounded height[END_REF] solve the case where the precedence graph is of bounded height by proposing an algorithm of time complexity O(n h(m-1)+1 ). Recently, [START_REF] Aho | On a parallel machine scheduling problem with precedence constraints[END_REF] show that P m|prec, p j = p|Cmax is solvable in polynomial time when the precedence graph is of bounded height and the maximum degree is bounded. This result is in fact a special case of the one proposed in Dolev

and Warmuth [1984].

Other criteria and/or machine environment

The other well-studied criterion for parallel machine environment is the total completion time C j .

One of the reasons for that is that, for some problems, it is equivalent to solve the total completion time and the makespan since they admit an ideal schedule. For example, ideal schedules exist when considering two machines, arbitrary precedence constraints and equal processing times, hence problem P 2|prec, p j = p| C j is polynomially solvable with CG-algorithm (see [START_REF] Coffman | Optimal scheduling for two-processor systems[END_REF]).

When adding release dates, [START_REF] Baptiste | On preemption redundancy in scheduling unit processing time jobs on two parallel machines[END_REF] show that P 2|prec, r j , p j = 1| C j is solvable in polynomial time by reducing it to a shortest path problem.

For an arbitrary number of machines, if the precedence graph forms an outtree and the processing times are UET, the same result holds and hence P |outtree, p j = p| C j is solvable in polynomial time.

Note that problem P 3|intree, p j = p| C j is not ideal, see [START_REF] Huo | Minimizing mean flow time for uet tasks[END_REF] for a counterexam-ple. Nevertheless, for any fixed number of machines, problem P m|intree, p j = p| C j is solvable in polynomial time (see Baptiste et al. [2004a]). Adding release dates maintains the same property: the algorithm proposed in [START_REF] Brucker | A polynomial algorithm for p-pj= 1, rj, outtreecj[END_REF] solves simultaneously problems P |outtree, r j , p j = 1|Cmax and P |outtree, r j , p j = 1| C j . An improvement of this algorithm has been proposed in [START_REF] Huo | Minimizing total completion time for uet tasks with release time and outtree precedence constraints[END_REF].

For interval order graphs, [START_REF] Möhring | Computationally tractable classes of ordered sets[END_REF] notices that P |io, p j = p| C j is solvable in polynomial time since the proof of the algorithm for P |io, p j = p|Cmax (in [START_REF] Papadimitriou | Scheduling interval-ordered tasks[END_REF])

only uses swaps between tasks, and this property is verified by the makespan and the total completion time. It has been noticed recently that the same result holds for overinterval orders (that properly contain interval orders), the problem admits also an ideal solution, so P |oio, p j = p| C j is solvable in polynomial time (see [START_REF] Wang | Ordonnancement parallèle non-préeemptif avec contraintes de précédence[END_REF]).

When considering uniform parallel machines, only few results are available; problem Q2|chains, p j = p|Cmax is solvable in polynomial time (see [START_REF] Brucker | Scheduling identical jobs with chain precedence constraints on two uniform machines[END_REF]). If one processor is going b times faster than the other (with b an integer), the problem Q2|tree, p j = p|Cmax is also polynomially solvable (see [START_REF] Kubiak | Optimal scheduling of unit-time tasks on two uniform processors under tree-like precedence constraints[END_REF]); the problem is also ideal, and hence Q2|tree, p j = p| C j is also solvable in polynomial time.

Minimum NP-hard cases

The interesting results for this survey are the N P-hardness of P |prec, p j = p|Cmax and P |prec, p j = p| C j (see [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF]). Note that the proof for the two problems also holds when the precedence graph is of bounded height, and that the problem P |opp.f orest, p j = p|Cmax is also N Phard (see [START_REF] Garey | Scheduling opposing forests[END_REF]). When adding release dates, the corresponding problem is already N P-hard for intrees, for both the makespan and the total completion time (see Brucker et al. [1977]). Finally, what may not be obvious is that even if P ||Cmax is strongly N P-hard, the problem P |w(G) ≤ K|Cmax is solvable in pseudo-polynomial time [START_REF] Middendorf | Transversal graphs for partially ordered sets: Sequencing, merging and scheduling problems[END_REF]). This results is due to the fact that the class of empty graphs is not included in bounded-width graphs class.

Open problems

For an arbitrary number of machines and the makespan criterion, the boundary between polynomially solvable and N P-hard problems seems very sharp, we believe that the efforts should not be focused on these problems. When the number of machines is fixed, this boundary is much larger. Surprisingly, to the best of our knowledge, no other structures of precedence graphs than the ones introduced in previous section have been studied for problem Pm|prec, p j = p|Cmax. In our opinion, it could be a good opportunity to work on more general precedence graphs on this problem, to be able to arbitrate if Pm|prec, p j = p|Cmax is solvable in polynomial time or N P-hard. The most natural extension in our opinion is to consider series-parallel graphs, since it is a generalization of opposing forests, level order graphs and DC-graphs, for which the problem is polynomially solvable.

For the total completion time criterion, the two most intriguing problems are P |intree, p j = p| C j and P |outtree, r j , p j = p| C j . For the former problem, the interest lies in the fact that there exists an ideal schedule for outtree precedences, but not for intree precedences. Nevertheless, we do believe that this problem admits an optimal polynomial algorithm. For the latter problem, an algorithm exists for p j = 1 (i.e., release dates are multiple of the processing time, see [START_REF] Brucker | A polynomial algorithm for p-pj= 1, rj, outtreecj[END_REF]), and hence the gap to p j = p seems small.

For a fixed number of uniform parallel machines and the makespan criterion, the set of open problems is wide, since the only polynomial algorithm is for Q2|chains, p j = p|Cmax [START_REF] Brucker | Scheduling identical jobs with chain precedence constraints on two uniform machines[END_REF]), and the problem Qm|prec, r j , p j = p|Lmax is still open. The same behavior occurs for the total completion time: [START_REF] Dessouky | Scheduling identical jobs on uniform parallel machines[END_REF] proved that Qm|r j , p j = p| C j is solvable in polynomial time, and problem Qm|prec, r j , p j = p| C j is still open. We hence believe that this set of problems deserves a deeper study. A first approach may consist in trying to adapt the algorithms available for identical parallel machines.

6 Parallel machines with preemption Timkovsky shows very strong links between preemption and chains, including the fact that a large set of scheduling problems with preemption can be reduced to problems without preemption, with UET tasks, and where each job is replaced by a chain of jobs (see Theorem 3.5 in Timkovsky [2003]). This interesting result can be applied in many cases, but does not hold for the total completion time criterion. Moreover, according to the structure of the precedence graph, the resulting graph may not have the same structure.

For example, an intree where each job is replaced by a chain of jobs remains an intree, whereas it does not hold for interval order graphs (two independent jobs will be replaced by two chains of parallel jobs, which is not an interval order graph). That is why we will examine more precisely what happens in this section.

Polynomial cases

Makespan criterion

Since P |tree, p j = p|Cmax is polynomially solvable (see [START_REF] Hu | Parallel sequencing and assembly line problems[END_REF]), by Timkovsky's result, so is P |pmtn, tree|Cmax.

The first polynomial algorithm for this problem is proposed in [START_REF] Muntz | Preemptive scheduling of real-time tasks on multiprocessor systems[END_REF] with a running time of O(n 2 ). An algorithm in O(n log m) was then proposed in Gonzalez and Johnson [1980].

Note that the latter algorithm produces at most O(n) preemptions whereas the former may obtain a schedule with O(nm) preemptions. [START_REF] Lawler | Preemptive scheduling of. precedence-constrained jobs on parallel machines[END_REF] studies the case with release dates and outtree, and

shows that it can be solved in O(n 2 ) with a dynamic priority list algorithm (i.e., priorities may change according to what has already been scheduled).

Timkovsky's result can not be applied to precedence graphs such as interval order graphs. Yet, it was proven that the problem is also solvable in polynomial time for this precedence structure; first, [START_REF] Sauer | Preemptive scheduling of interval orders is polynomial[END_REF] show it for a fixed number of machines P m|pmtn, io|Cmax by using a linear programming approach based on the set of jobs scheduled at each instant. Later [START_REF] Djellab | Scheduling preemptive jobs with precedence constraints on parallel machines[END_REF] proposes another linear program that solves the problem for an arbitrary number of machines, i.e., P |pmtn, io|Cmax. For a fixed number m of machines, [START_REF] Moukrim | Optimal preemptive scheduling on a fixed number of identical parallel machines[END_REF] extend the result of [START_REF] Sauer | Preemptive scheduling of interval orders is polynomial[END_REF], by proposing an linear programming approach for Am-order graphs (which properly contain interval order graphs).

Other criteria and/or machine environment [START_REF] Du | Scheduling chain-structured tasks to minimize makespan and mean flow time[END_REF] explain that P 2|pmtn, chains| C j is strongly N P-hard by showing that preemption is useless for this problem, that is why we will only focus on the UET case for the total completion time criterion. [START_REF] Chardon | The coffman-graham algorithm optimally solves uet task systems with overinterval orders[END_REF] prove that P 2|pmtn, prec, p j = p| C j is solvable in polynomial time by adapting the Coffman-Graham algorithm. Their algorithm finds an ideal schedule. By proving that preemption is redundant, [START_REF] Baptiste | On preemption redundancy in scheduling unit processing time jobs on two parallel machines[END_REF] prove that P 2|pmtn, outtree, r j , p j = 1| C j is solvable in polynomial time. When the precedence graph is an intree, Coffman Jr et al. [2012] prove that the problem is not ideal, and [START_REF] Chen | Normal-form preemption sequences for an open problem in scheduling theory[END_REF] provide a deep analysis of the structure of preemption. Using the same methodology than [START_REF] Baptiste | On preemption redundancy in scheduling unit processing time jobs on two parallel machines[END_REF], [START_REF] Brucker | A polynomial algorithm for p-pj= 1, rj, outtreecj[END_REF] show that the problem is solvable in polynomial time with an outtree and an arbitrary number of machines.

Moreover, they provide a O(n 2 ) algorithm, that admits a O(n log n) implementation according to [START_REF] Huo | Minimizing total completion time for uet tasks with release time and outtree precedence constraints[END_REF]. [START_REF] Lushchakova | Two machine preemptive scheduling problem with release dates, equal processing times and precedence constraints[END_REF] slightly improves the result of [START_REF] Baptiste | On preemption redundancy in scheduling unit processing time jobs on two parallel machines[END_REF] and proposes an algorithm of complexity O(n 2 ) for the problem P 2|pmtn, outtree, r j , p j = p| C j .

Minimum NP-hard cases

For the makespan criterion, [START_REF] Ullman | Complexity of sequencing problems[END_REF] shows that P |pmtn, prec, p j = p|Cmax is N P-hard. For the total completion time criterion, if we carefully look at the N P-hardness proof of P |prec, p j = p| C j in [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF], we can see that preemption is useless for the instance constructed in the reduction and hence the preemptive version is still N P-hard for precedence graphs of bounded height:

Theorem 1 P |pmtn, prec, p j = p| C j is N P-hard even if the precedence graph is of bounded height.

Proof. We just need to modify slightly the proof in the reduction from Clique of [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF].

Clique : G = (V, E) is an undirected graph and k an integer. Does G have a clique on k vertices?

Let us recall this reduction; we denote by v (resp. e) the number of vertices (resp. edges) of G. We also define following parameters: l = k(k-1) 2

, k = v -k and l = e -l (we use the notations of the original article). We construct an instance of P |pmtn, prec, p j = p| C j with m = max{k, l + k , l } + 1 machines and n = 3m jobs:

for each vertex i ∈ V , there is a job J i .

for each edge [i, j] ∈ E, there is a job J

[i,j] . -dummy jobs J h,t with h ∈ D t , t ∈ {1, 2, 3}, D 1 = {1, . . . , m -k}, D 2 = {1, . . . , m -l -k }, D 3 = {1, . . . , m -l }.
There are precedence constraints between any job J g,t and J h,t+1 , for any g ∈ D t , h ∈ D t+1 , t = 1, 2.

Moreover, for any edge [i, j] ∈ E, there is a precedence between J i and J [i,j] . Finally, the question is whether there is a schedule such that C j ≤ 6m.

Clearly, if Clique has a solution, so is the scheduling problem:

dummy jobs J h,t are scheduled during the time-interval [t -1, t],

the k jobs corresponding to the vertices of the clique are scheduled during the first period,

the l jobs corresponding to the edges of the clique, and the k jobs corresponding to the remaining vertices are scheduled at the second period, all the remaining jobs are scheduled at the third period.

This solution has a total completion time of exactly 6m. Note that it does not use preemption, and is such that Cmax = 3.

Conversly, let us prove that if there is a schedule such that C j ≤ 6m, then there is a clique of size k. First, we can easily prove that if there is a schedule such that C j ≤ 6m, then there is a schedule such that Cmax ≤ 3. Indeed, if there exists a job j * such that C j * > 3, then in the best case (i.e., if there is no precedence constraint), we have

C j ≥ 1 * m + 2 * m + 3 * (m -1) + C j * > 6m.
We hence know that there is no idle time and that the dummy jobs J h,t are scheduled during interval [t -1, t]. To conclude, we just need to see that if there is no clique of size k then, whatever the schedule on interval [0, 1] without idle time, the number of eligible jobs at time 1 is strictly less than k + l, which implies an idle time and hence no schedule such that Cmax ≤ 3.

Open problems

Preemptive parallel machine scheduling problems did not receive as much attention as their non-preemptive counterpart, hence the set of open problems is wider.

For the makespan objective and an arbitrary number of machines, when the precedence graph is an interval order graph, it is known to be solvable in polynomial time, as in the UET non-preemptive case.

Since for the UET non-preemptive problem, new classes strictly including interval order graphs (namely quasi-interval order graphs and over-interval order graphs) have led to polynomial algorithms, the same question arises for P |pmtn, qio|Cmax and P |pmtn, oio|Cmax.

When the number of machines is fixed, there is a wide set of open problems, since Qm|pmtn, prec, r j |Lmax

is the maximal open problem. A good challenge may be for example to try to fix the complexity of P m|pmtn, prec, p j = p|Cmax, or at least to try to find new precedence graphs for which the problem is polynomial, by taking advantage of the fact that tasks are UET (even if it may not always be helpful with preemption).

For the total completion time criterion, P |pmtn, outtree, p j = p| C j is maximal polynomially solvable, and we just show that P |pmtn, prec, p j = p| C j is N P-hard even if the precedence graph is of bounded height. It could be interesting to consider other structures of precedence graph to derive polynomial algorithms. In a similar way, Baptiste et al. [2004a] show that P |pmtn, p j = p| T j is solvable in polynomial time, and adding precedence constraints makes the problem N P-hard, but there is no other result available in the literature, it hence would be interesting to search for new polynomial cases by testing different precedence graphs.

For the problem P 2|pmtn, prec, r j , p j = 1| C j , the gap is even wider: it is polynomial if the precedence graph is an outtree (see [START_REF] Baptiste | On preemption redundancy in scheduling unit processing time jobs on two parallel machines[END_REF]) but all the other cases are open, from P 2|pmtn, intree, r j , p j = 1| C j to P 2|pmtn, prec, r j , p j = 1| C j .

Conclusion

In this paper, we surveyed the complexity results for scheduling problems with precedence constraints, and we have seen that single machine scheduling problems have been much more studied than others.

This looks quite normal since the single machine problem is the scheduling problem that is the closest to order theory. Nevertheless, we show that there still are a few open problems for the single machine case.

We believe that the most interesting problems for which the complexity is open lie in the parallel machine case; more precisely, we do conjecture that P m|sp -graph, p j = p|Cmax is solvable in polynomial time; this result will be a large breakthrough since series-parallel graphs are most of the time studied for single machine problems.

Another approach to understand the complexity of scheduling problems is to deal with the parameterized complexity (see [START_REF] Downey | Parameterized complexity[END_REF]). There are only very few results on parameterized complexity of scheduling problems. One can cite [START_REF] Michael | On the parametric complexity of schedules to minimize tardy tasks[END_REF] who show that, if the precedence graph is of bounded width w (it is equal to the size of a maximum antichain), then problem 1|prec| T j ≤ k is FPT when parameterized by (w, k). The most recent result on the subject is that P ||Cmax is FPT for parameter pmax = max p j (see [START_REF] Mnich | Scheduling and fixed-parameter tractability[END_REF]). A good graph measure is a powerful tool for parameterized complexity. For general (undirected) graphs, the creation of the treewidth (see [START_REF] Robertson | Graph minors. ii. algorithmic aspects of tree-width[END_REF]) helped to discover many results in graph theory, including the Courcelle's theorem [START_REF] Courcelle | The monadic second-order logic of graphs. i. recognizable sets of finite graphs[END_REF]). For directed graphs, and more specifically DAGs, none of the existing measures (see [START_REF] Ganian | Digraph width measures in parameterized algorithmics[END_REF]) is satisfactory. In our opinion, a major breakthrough will be achieved when one will be able to find a good measure on directed acyclic graphs.

Fig. 2 A

 2 Fig.1The height of a DAG

  Figure 4.

  Fig. 5 A DC-graph

Fig. 9

 9 Fig.8The three forbidden subgraphs for quasi-interval order graphs

Fig. 10

 10 Fig. 10 Hasse diagram for different classes of directed acyclic graphs

  opposing forest), io (interval order graph), qio (quasi-interval order graph), oio (over-interval order graph), sp -graph (series-parallel graph), DC -graph (divide-and-conquer graph), lo (level order graph), h(G) ≤ k (directed acyclic graph with height bounded by k), dim ≤ k (directed acyclic graph with dimension bounded by k), f dim ≤ k (directed acyclic graph with fractional dimension bounded by k), w(G) ≤ k (directed acyclic graph with width bounded by k).

Appendix A: List of results

For an easier reading of all the complexity results that are reviewed in this survey, we proposed a synthesis in the following tables. In each table, we write the polynomial cases, some open cases (the ones that seem the most promising in our opinion) and the N P-hard problems.

Problem

Complexity Reference 1|sp -graph| w j C j P Lawler [1978] 1|dim ≤ 2| w j C j P Ambühl and Mastrolilli [2009] 1|sp -graph, p [START_REF] Wang | Single-machine scheduling with deteriorating jobs under a series-parallel graph constraint[END_REF] 1|sp -graph, p j (t) = p j + α j t|Cmax P [START_REF] Wang | Single-machine scheduling with deteriorating jobs under a series-parallel graph constraint[END_REF] 1|sp -graph, p

[r] j = a j + b j r|Cmax P Wang and [START_REF] Wang | Single-machine scheduling with precedence constraints and positiondependent processing times[END_REF] 1|sp -graph, p

[r] j = a j -b j r|Cmax P Wang and [START_REF] Wang | Single-machine scheduling with precedence constraints and positiondependent processing times[END_REF] 1|sp -graph, p

[r] j = p j r|Cmax P Gordon et al. [2008] 1|sp -graph, p

[r] [START_REF] Gordon | Single machine scheduling models with deterioration and learning: handling precedence constraints via priority generation[END_REF] 1|sp -graph, p [START_REF] Gordon | Single machine scheduling models with deterioration and learning: handling precedence constraints via priority generation[END_REF] 1|sp -graph, p [START_REF] Gordon | Single machine scheduling models with deterioration and learning: handling precedence constraints via priority generation[END_REF] 1|sp -graph, p [START_REF] Gordon | Single machine scheduling models with deterioration and learning: handling precedence constraints via priority generation[END_REF] 1|pmtn, r j , p j = p| T j P Tian et al. [2006] 1|pmtn, r j , p j = p| w j U j P Baptiste [1999] 1|2 < f dim < 3| w j C j Open 1|io| C j Open 1|pmtn, chains, r j , p j = p| T j Open 1|pmtn, prec, r j , p j = p| w j T j Open 1|pmtn, chains, r j , p j = p| U j Open 1|pmtn, prec, r j , p j = p| w j U j Open 1|prec, p j = p| w j C j N P-hard [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF] 1|prec| C j N P-hard [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF] 1|indegree ≤ 2| C j N P-hard [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF] 1|h(G) ≤ k| w j C j N P-hard [START_REF] Lawler | Sequencing jobs to minimize total weighted completion time subject to precedence constraints[END_REF] 1|io| w j C j N P-hard Ambühl et al. [2011] 1|f dim ≥ 3| w j C j N P-hard Ambühl et al. [2011] 1|chains, p j = 1| U j N P-hard [START_REF] Karel | Complexity results for scheduling chains on a single machine[END_REF] 1|chains, p j = 1| T j N P-hard [START_REF] Leung | Minimizing total tardiness on a single machine with precedence constraints[END_REF] P m|sp -graph, p j = p|Cmax Open P m|prec, p j = p|Cmax Open P |intree, p j = p| C j Open P |outtree, r j , p j = p| C j Open Qm|prec, r j , p j = p|Lmax Open Qm|prec, r j , p j = p| C j Open P 2||Cmax N P-hard [START_REF] Karel Lenstra | Complexity of machine scheduling problems[END_REF] P 2|chains| C j N P-hard [START_REF] Du | Scheduling chain-structured tasks to minimize makespan and mean flow time[END_REF] P |opp.f orest, p j = p|Cmax N P-hard [START_REF] Garey | Scheduling opposing forests[END_REF] P |h(G) ≤ k, p j = p|Cmax N P-hard [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF] P |intree, r j , p j = p|Cmax N P-hard Brucker et al. [1977] P |h(G) ≤ k, p j = p| C j N P-hard [START_REF] Karel | Complexity of scheduling under precedence constraints[END_REF] P |intree, r j , p j = 1| C j N P-hard Lenstra P 2|pmtn, chains, p j = 1| U j N P-hard Baptiste et al. [2004a] Table 3 Complexity results for parallel machine problems with preemption