Albert Kouchéré Guidzavaï

Hamadjam Abboubakar
email: h.abboubakar@gmail.com

Irépran Damakoa

Nabil Gmati

Mathieu Roche

Tri Nguyen-Huu

Laurent Debreu

Analysis of a mosquito life cycle model

Keywords: Mosquito life cycle, Gonotrophic cycle, Cooperative systems, Sterile Insect Technique (SIT)

The gonotrophic cycle of mosquitoes conditions the frequency of mosquito-human contacts. The knowledge of this important phenomenon in the mosquito life cycle is a fundamental element in the epidemiological analysis of a communicable disease such as mosquito-borne diseases. In this work, we analyse a deterministic model of complete life cycle of mosquitoes which takes into account the principal phases of female mosquitoes gonotrophic cycle, and the Sterile Insect technique combined with the use of insecticide as control measures to fight the proliferation of mosquitoes. We compute the corresponding mosquito reproductive number N * and prove the global asymptotic stability of trivial equilibrium. We prove that the model admits two non-trivial equilibria whenever N * is greater than another threshold, N c , which the total number of sterile mosquitoes depends on. Numerical simulations, using mosquito parameters of the Aedes species, are carried out to illustrate our analytical results and permit to show that the strategy which consists in combining the sterile insect technique with adulticides, when it is well done, effectively combats the proliferation of mosquitoes.

I INTRODUCTION

Many diseases among humans are transmitted by mosquitoes and sandflies. According to the World Health Organization (WHO), over a billion people are infected and over a million die from these vector-borne diseases annually [START_REF] Who | Global brief on vector-borne diseases[END_REF][START_REF] Who | Global vector control response 2017-2030[END_REF]. Although some of these diseases have effective vaccines (yellow fever) or effective curative drugs (malaria and lymphatic filariasis), vector-borne diseases remain real public health problems in developing countries. The effective means of preventing the transmission of most of these vector-borne diseases remain vector control mechanisms. Consequently, a better understanding of the vector life cycle is necessary in order to propose better control mechanisms and thus reduce their proliferation.

The gonotrophic cycle of mosquitoes conditions the frequency of mosquito-human contacts. The knowledge of this important phenomenon in the mosquito life cycle is a fundamental element in the epidemiological analysis of a communicable disease such as mosquito-borne diseases. The frequency of these mosquito-host contacts is indeed dependent on certain factors such as availability of vertebrate hosts, rainfall and temperature.But it depends essentially on the time required for the digestion of the blood meal which is accompanied by the maturation of the ovaries and oviposition and followed by the search for a new blood meal. For Beklemishev [START_REF] Beklemishev | Le cycle gonotrophique, principe de base de la biologie de an. gambiae, Vop Fiziol[END_REF], the gonotrophic cycle would define the interval between the emergence of the first spawning in the case of nulliparous females and two successive ovipositions in the case of parous females (see also [START_REF] Carnevale | Etude du cycle gonotrophique d'anopheles gambiae (diptera, culicidae)(giles, 1902) en zone de forêt dégradée d'afrique centrale[END_REF]). In general, the gonotrophic cycle consists of following distinct phases [START_REF] Arifin | An agentbased model of the population dynamics of anopheles gambiae[END_REF][START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF]:

1. blood meal seeking: search of a vertebrate host by a young female mosquito after its mating with a male mosquito, or an adult female mosquito which has just laid eggs in a breeding site, 2. blood meal digesting: digestion of ingested blood and ovarian maturation, 3. gravid: laying eggs by pregnant females (gravid) in a favourable breeding site, 4. search of another blood meal, . . . Figure 1 presents the principal stages of mosquito life cycle. The time between mating and a Figure 1: Life cycle of female mosquitoes [START_REF] Coon | Mosquitoes rely on their gut microbiota for development[END_REF].

successful blood meal seeking is a key factor in the female mosquito life who increases their mortality rate in the finding of blood to ensure development of eggs. For many mosquitoes species such as Anopheles sp., the search of a blood meal is between 06 pm to 06 am [START_REF] Arifin | An agentbased model of the population dynamics of anopheles gambiae[END_REF][START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF]. Thus, it is in this time interval that some vector control mechanisms such as treated insecticide bed-nets usage and insecticides permit to decrease the number of mosquito-human contacts or increase the mosquito death rate. The technique called sterile insect technique (SIT), was developed to stop the mosquitoes proliferation [START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF]. Indeed, this mosquito control mechanism consists to release around breeding sites, male mosquitoes which was irradiated, in order for them to compete with wild male mosquitoes and mate with female mosquitoes. So, the eggs laying by female mosquitoes who had mated with sterile males will not be able to hatch. Thus, the next generation of mosquito population will decrease according to the number of sterile males mosquitoes whom mated with wild female mosquitoes. This technique has been used to control the proliferation of Aedes albopictus population during the 2005 chikungunya epidemic in Reunion Island [START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF].

There is a long history in the mathematical modelling of life cycle of mosquitoes. After the chikungunya epidemic outbreak in 2005, many mathematical models, which take into account some principal stages of mosquito life cycle, were proposed (see for example [START_REF] Abboubakar | Bifurcation thresholds and optimal control in transmission dynamics of arboviral diseases[END_REF][START_REF] Abboubakar | Backward bifurcation and control in transmission dynamics of arboviral diseases[END_REF][START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF][START_REF] Arifin | An agentbased model of the population dynamics of anopheles gambiae[END_REF][START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF][START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF][START_REF] Moulay | The chikungunya disease: modeling, vector and transmission global dynamics[END_REF][START_REF] Yang | Follow up estimation of aedes aegypti entomological parameters and mathematical modellings[END_REF][START_REF] Yusoff | Stage-structured population dynamics of aedes aegypti[END_REF]). In [START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF][START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF], the authors propose some various models with aquatic development phase of Aedes sp. mosquitoes coupled with the transmission dynamics of chikungunya virus in human populations. To prove the impact of sterile insect techniques on the decrease of Aedes population during the 2005 chikungunya epidemic in Reunion island, Dumont et al. [START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF][START_REF] Dumont | Vector control for the chikungunya disease[END_REF] proposed and studied two compartmental models. They did not modelled the complete mosquito life cycle. In the same idea, Moulay et al. [START_REF] Moulay | The chikungunya disease: modeling, vector and transmission global dynamics[END_REF] (see also Yang et al. [START_REF] Yang | Follow up estimation of aedes aegypti entomological parameters and mathematical modellings[END_REF] and Yusoff et al. [START_REF] Yusoff | Stage-structured population dynamics of aedes aegypti[END_REF]) proposed also a compartmental model in which they modelled eggs and larvae to represent mosquito aquatic phase of Aedes albopictus. They did not take into account the pupae stage and the gonotrophic cycle of mosquitoes. Abboubakar et al. [START_REF] Abboubakar | Bifurcation thresholds and optimal control in transmission dynamics of arboviral diseases[END_REF][START_REF] Abboubakar | Backward bifurcation and control in transmission dynamics of arboviral diseases[END_REF] extended the models of Dumont et al. [START_REF] Dumont | Vector control for the chikungunya disease[END_REF] and Moulay et al. [START_REF] Moulay | The chikungunya disease: modeling, vector and transmission global dynamics[END_REF] by adding the pupa stage in their models. Arifin et al. [START_REF] Arifin | An agentbased model of the population dynamics of anopheles gambiae[END_REF] and Gentile et al. [START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF] proposed two agent-based models to study the dynamic behavior of Anopheles gambiae. Here, we present the corresponding ODE model of Arifin et al. [START_REF] Arifin | An agentbased model of the population dynamics of anopheles gambiae[END_REF] and Gentile et al. [START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF] by the representation of each mosquito life cycle phase in a different compartment, and in which we include the sterile insect technique and the use of insecticide to decrease the population of mosquitoes. We prove the global stability of the trivial equilibrium and show that whenever the mosquito reproductive number N * is greater than one and greater than an other threshold, N c , which depends on the total number of sterile mosquitoes, the model admits two positive equilibria. Numerical simulations are performed to illustrate our analytical results.

The paper is organised as follow. The description of the complete life cycle model of mosquitoes (including gonotrophic cycle, SIT and the use of insecticide) and its mathematical analysis are devoted in Section II. Numerical simulations are performed in Section III.

II THE MODEL AND ITS ANALYSIS

Here, we describe the general mosquito life cycle model.

Aquatic phase

Aquatic phase is divided in three compartments: E for eggs, L for larvae and P for pupae. Eggs are laid with a rate µ b . After 2 to 4 days depending of ambient temperature, eggs become larvae at a rate s. The larvae will in turn become pupae at a rate l. µ E , µ L and µ P denotes the mortality rate of eggs, larvae and pupae, respectively. We assume in this work that, based on the larval population size, the density-dependence affect the survival rates of both larvae and pupae with linear density functions αL and βL, respectively [START_REF] Ai | Mosquito-stage-structured malaria models and their global dynamics[END_REF].

Juvenile stage and mosquito ready to mate

They are two compartments for juvenile mosquitoes. We denote by X the total number of male mosquitoes and by Y the total number of female mosquitoes. Indeed, after the completed aquatic phase, the pupae will become either male mosquitoes X at a rate (1ϑ)η or female mosquitoes Y at a rate ϑη, where ϑ is the sex ratio and η is the pupae development rate. Then comes the phase where a male mosquito will mate with a female mosquito F at a constant rate β 2 . We denote by µ X , µ Y and µ F the mortality rates of male mosquitoes, juvenile female mosquitoes and adult female mosquitoes, respectively.

The gonotrophic cycle

We subdivide the gonotrophic cycle in three compartments: S for total number of female mosquitoes in blood meal seeking stage, D for total number of female mosquitoes in blood meal digesting stage, and G for total number of female mosquitoes in gravid stage. After mating with a male, female mosquito will go to search a vertebrate host for a blood meal rich in protein which will ensure the proper development of their eggs, at a rate ω. Female mosquitoes in blood meal seeking stage die at constant rate µ S . After a successful blood meal, seeking mosquitoes find a place to digest at constant rate γ 1 . Female mosquitoes in blood meal digesting stage die at constant rate µ D , and become gravid G at the constant rate γ 2 . Gravid mosquitoes die at a constant rate µ G . In this stage, gravid mosquitoes find an appropriate breeding site to lay their eggs, and search another vertebrate host to find a new blood meal at a rate ǫ, and then enter in another gonotrophic cycle.

The above assumptions lead to the following nonlinear system of ordinary differential equations

Ė(t) = µ b G(t) -(s + µ E )E(t), (1a) L(t) = sE(t) -(l + µ L + αL(t))L(t), (1b) Ṗ (t) = lL(t) -(η + µ P + β 1 L(t))P (t), (1c) Ẋ(t) = (1 -ϑ)ηP (t) -µ X X(t) (1d) Ẏ (t) = ϑηP (t) -(µ Y + β 2 )Y (t), (1e)
Ḟ (t) = β 2 Y (t) -(µ F + ω)F (t), (1f) Ṡ(t) = ωF (t) + ǫG(t) -(µ S + γ 1 )S(t), (1g) Ḋ(t) = γ 1 S(t) -(µ D + γ 2 )D(t), (1h) Ġ(t) = γ 2 D(t) -(µ G + ǫ)G(t). ( 1i
)
Without loss of generalities, we claim the following results.

Lemma 1:

Every solution ψ(t) = ((E(t), L(t), P (t), X(t), Y (t), F (t), S(t), D(t), G(t))) of ( 1) with its initial conditions in R 9 + is defined and lies in int(R 9 + ) for all t > 0.

Proof. Let f be the right-hand side of (1) and

k 1 = s + µ E , k 2 = l + µ L , k 3 = µ P + η, k 4 = β 2 + µ Y , k 5 = µ F + ω, k 6 = µ S + γ 1 , k 7 = µ D + γ 2 and k 8 = µ G + ε. The Jacobian matrix of f is Df (E, L, P, X, Y, F, S, D, G) =               -k 1 0 0 0 0 0 0 0 µ b s -(k 2 + 2αL) 0 0 0 0 0 0 0 0 l -β 1 P -(k 3 + β 1 L) 0 0 0 0 0 0 0 0 (1 -ϑ)η -µ X 0 0 0 0 0 0 0 ϑη 0 -k 4 0 0 0 0 0 0 0 0 β 2 -k 5 0 0 0 0 0 0 0 0 ω -k 6 0 ε 0 0 0 0 0 0 γ 1 -k 7 0 0 0 0 0 0 0 0 γ 2 -k 8               . (2)
To prove Lemma 1, we proceed by step. Since ( 1) is a smooth system in R 9 , ψ(t) exists on a maximal interval [0, T + ) for some 0 < T + ≤ ∞.

Step 1 We show that ψ(t) > 0 for t ∈ [0, T + ). Using the variation of parameters formula, we have, for t ∈ [0, T + )

E(t) = E(0)e -k 1 t + µ b t 0 G(τ )e -k 1 (t-τ ) dτ, L(t) = L(0)e -t 0 (k 2 +αL(τ ))dτ + s t 0 E(τ )e -t τ (k 2 +αL(η))dη dτ, P (t) = P (0)e -t 0 (k 3 +β 1 L(τ ))dτ + l t 0 L(τ )e -t τ (k 3 +β 1 L(η))dη dτ, X(t) = X(0)e -µ X t + (1 -ϑ)η t 0 P (τ )e -µ X (t-τ ) dτ, Y (t) = Y (0)e -k 4 t + ϑη t 0 P (τ )e -k 4 (t-τ ) dτ, F (t) = F (0)e -k 5 t + β 2 t 0 Y (τ )e -k 5 (t-τ ) dτ, S(t) = S(0)e -k 6 t + ω t 0 F (τ )e -k 6 (t-τ ) dτ + ε t 0 G(τ )e -k 6 (t-τ ) dτ, D(t) = D(0)e -k 7 t + γ 1 t 0 S(τ )e -k 7 (t-τ ) dτ, G(t) = G(0)e -k 8 t + γ 2 t 0 D(τ )e -k 8 (t-τ ) dτ, (3)
Since ψ(t) ∈ R 9 + \ {0}, we have E(0) > 0, L(0) > 0, P (0) > 0, X(0) > 0, Y (0) > 0,F (0) > 0,S(0) > 0, D(0) > 0 or G(0) > 0. If E(0) > 0, it follows from the continuity of E that E(t) > 0, for t ∈ (0, ǫ), for small ǫ > 0. Then it follows from (3) that L(t) > 0 on (0, ǫ], which, in turn, implies that P (t) > 0, X(t) > 0, Y (t) > 0,F (t) > 0, S(t) > 0, D(t) > 0 on (0, ǫ], and then G(t) > 0 on (0, ǫ]. Note that for any t ∈ [ǫ, T + ), as long as E, L, P , X, Y , F , S, D, G are positive on (0, t), all the integrals on the right-hand side of (3) are positive. It then follows from (3) that E(t) > 0, L(t) > 0, P (t) > 0, X(t) > 0, Y (t) > 0, F (t) > 0, S(t) > 0, D(t) > 0 and G(t) > 0, for any t ∈ [ǫ, T + ), as long as E, L, P , X, Y , F , S, D,G, are positive on (0, t). This yields that E, L, P , X 1 , Y , F , S, D and G are positive on (0, T + ). The same assertion holds in the cases that L(0) > 0,

P (0) > 0,X(0) > 0, Y (0) > 0, F (0) > 0, S(0) > 0, D(0) > 0 or G(0) > 0.
Step 2 We next show that T + = +∞. To this aim, we compare ψ(t) with a solution of the variational system of (1) at its trivial equilibrium

ż = Df (0)z (4)
which is a cooperative system in R 9 . Let z(t) be the solution of (4), with z(0

) = ψ(0). Since ψ(t) ∈ R 9 + \ {0}, for t ∈ [ǫ, T + ), it follows, for t ∈ [ǫ, T + ), that ψ(t) = f (ψ(t)) = Df (0)ψ(t)-0, αL 2 (t), β 1 L(t)P (t), 0, 0, 0, 0, 0, 0 T ≤ Df (0)ψ(t).
Then, from Theorem 4, it follows that ψ(t) ≤ z(t), for t ∈ [0, T + ). Since z(t) is defined for all t ≥ 0, it follows that ψ(t) is defined for all t ∈ (0, +∞), i.e. T + = +∞.

The results from Steps 1-2 permit us to conclude the proof of lemma 1.

Define the net reproductive number [START_REF] Moulay | The chikungunya disease: modeling, vector and transmission global dynamics[END_REF], given by

N * = s (s + µ E ) l (l + µ L ) ϑη (η + µ P ) µ b β 2 (β 2 + µ Y ) ω (ω + µ F ) γ 1 γ 2 (k 6 k 7 k 8 -εγ 1 γ 2 ) , (5)
where

k 6 k 7 k 8 -ǫγ 1 γ 2 = γ 1 µ D ǫ + (µ D + γ 2 ) (µ G (µ S + γ 1 ) + µ S ǫ) > 0.
N * measures the average expected number of new adult female offsprings produced by a single female vector during its life time. It can be interpreted as the product of the fraction of eggs that survived and hatched into larvae s/(s + µ E ), the fraction of larvae that survived and progressed into pupae l/(l + µ L ), the fraction of pupae that survived to become juvenile female mosquitoes ϑη/(η + µ P ), the fraction of juvenile female mosquitoes that survived to become adult female mosquitoes µ b β 2 /(β 2 +µ Y ), the fraction of the adult females which survived to give mosquitoes in search of blood meal ω/(ω + µ F ), the fraction of the mosquitoes in search of blood meal which survived to give mosquitoes in gravid (pregnant) state

γ 1 γ 2 /(k 6 k 7 k 8 -ǫγ 1 γ 2 ).
To decrease the proliferation of mosquito population, we use the biological mosquito control called sterile insect technique (SIT) combined with the use of insecticide and treated insecticide bed-nets which permits to increase the mortality rate of mosquitoes in blood meal seeking stage. To this aim, we firstly add one compartment X 2 of sterile male mosquitoes, which would compete with wild male mosquitoes to mate with female mosquitoes. Sterile males are released into the wild in an amount 5 to 10 times greater than that of wild males [START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF]. It only takes one mating of a female mosquito with a sterile male to reduce their fertility since it stores the sperm of the first male which fertilized their to conceive all their eggs throughout their life (about a month). Secondly, we modify the mortality rate of mosquitoes in blood meal seeking stage by adding a constant c m which represents the insecticide killing rate (µ S := µ S + c m ). Thus, the complete life cycle mosquito model with the introduction of sterile male compartment is as follow.

Ė(t) = µ b G(t) -(s + µ E )E(t), (6a)
L(t) = sE(t) -(l + µ L + αL(t))L(t), (6b)
Ṗ (t) = lL(t) -(η + µ P + β 1 L(t))P (t), (6c) Ẋ1 (t) = (1 -ϑ)ηP (t) -µ X 1 X 1 (t), (6d) Ẋ2 (t) = f rΛ 1 -µ X 2 X 2 (t), (6e)
Ẏ (t) = ϑηP (t) -(µ Y + β 2 )Y (t), (6f)
Ḟ (t) = β 2 eX 1 (t) eX 1 (t) + πX 2 (t) Y (t) -(µ F + ω)F (t), (6g)
Ṡ(t) = ωF (t) + ǫG(t) -(µ S + γ 1 + c m )S(t), (6h) Ḋ(t) = γ 1 S(t) -(µ D + γ 2 )D(t), (6i) Ġ(t) = γ 2 D(t) -(µ G + ǫ)G(t). (6j)
In Eq. ( 6), f rΛ 1 represent the efficient release of sterile males with f which denotes the competitivity of the sterile male, r is a parameter that represents the quality of the release, and Λ 1 denotes the rate of release of sterile males [START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF]. It is important to note that we follow Dumont and Tchuenche [START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF] and suppose that the introduction of sterile males only impacts the incoming of immature females in the wild females compartment F .

The trivial equilibrium and its stability analysis

The model system (6) admits a trivial equilibrium Q 1 = (0, 0, 0, 0, X * 2 , 0, 0, 0, 0, 0) whenever N * ≤ 1, where N * is the net reproduction number given by ( 5) with c m = 0, and

X * 2 = f rΛ 1 µ X 2
.

The stability analysis of the trivial equilibrium of model system ( 6) is given by the following result.

Theorem 1:

The trivial equilibrium point Q 1 = (0, 0, 0, 0, X * 2 , 0, 0, 0, 0, 0) of ( 6) is globally asymptotically stable in R 10 + if N * ≤ 1, and unstable if N * > 1.

Proof. The Jacobian matrix of system (6) at the trivial equilibrium Q 1 , with the components being reordered to Q = (E, L, P, X 1 , X 2 , Y, F, S, D, G), has the following form

J (Q 1 ) = A 11 A 12 0 A 22 ,
where

A 11 =         -k 1 0 0 0 0 0 s -k 2 0 0 0 0 0 l -k 3 0 0 0 0 0 (1 -ϑ)η -µ X 1 0 0 0 0 0 0 -µ X 2 0 0 0 ϑη 0 0 -k 4         , A 12 =     0 0 0 µ b 0 0 0 0 0 0 0 0 0 0 0 0     and A 22 =     -k 5 0 0 0 ω -k 6 0 ǫ 0 γ 1 -k 7 0 0 0 γ 2 -k 8     .
The eigenvalues of J (Q 1 ) consist of those of A 11 and A 22 . The eigenvalues of A 11 are given by λ 1 = -k 1 , λ 2 = -k 2 , λ 3 = -k 3 , λ 4 = -µ X 1 , λ 5 = -µ X 2 and λ 6 = -k 4 , which are all negative. It easy to see , using Routh-criterion that A 22 has eigenvalues with negative real part whenever N ≤ 1, where

N = s (s + µ E ) l (l + µ L ) ϑη (η + µ P ) µ b β 2 (β 2 + µ Y ) 1 µ F
. For the global stability, let us consider the following Lyapunov function

L(x, X 2 ) = 9 i=1 a i x i + 1 2 (X 2 -X * 2 ) 2 , (7)
where x = (E, L, P, X 1 , Y, F, S, D, G) and

a 1 = k 4 k 5 k 9 N * , a 2 = k 1 k 4 k 5 k 9 N * /s, a 3 = k 1 k 2 k 4 k 5 k 9 N * /sl, a 4 = N * β 2 ϑµ b ωγ 1 γ 2 , a 5 = N * β 2 ϑµ b ωγ 1 γ 2 , a 6 = N * k 4 ϑµ b ωγ 1 γ 2 , a 7 = N * k 4 k 5 ϑµ b γ 1 γ 2 , a 8 = N * k 4 k 5 k 6 ϑµ b γ 2 and a 9 = N * k 4 k 5 k 6 k 7 ϑµ b .
The Lyapunov derivative of ( 7) is given by

L(x, X 2 ) = 9 i=1 a i ẋi + (X 2 -X * 2 ) Ẋ2 = -(k 1 k 4 k 5 k 9 N * /s)αL 2 -(k 1 k 2 k 4 k 5 k 9 N * /sl)β 1 LP -µ X 2 (X 2 -X * 2 ) 2 + k 1 k 2 k 3 k 4 k 5 k 9 N * sl (N * -1)P -N * β 2 ϑµ b ωγ 1 γ 2 µ X 1 X 1 -N * k 4 β 2 ϑµ b ωγ 1 γ 2 Y + N * k 4 ϑµ b ωγ 1 γ 2 β 2 eX 1 eX 1 + πX 2 Y ≤ -(k 1 k 4 k 5 k 9 N * /s)αL 2 -(k 1 k 2 k 4 k 5 k 9 N * /sl)β 1 LP -µ X 2 (X 2 -X * 2 ) 2 -N * β 2 ϑµ b ωγ 1 γ 2 µ X 1 X 1 + k 1 k 2 k 3 k 4 k 5 k 9 N * sl (N * -1)P (8) It follows that L(x, X 2 ) ≤ 0 for N * ≤ 1. The largest compact invariant set in {(E, L, P, X 1 , X 2 , Y, F, S, D, G) ∈ R + 10 : L(x, X 2 ) = 0} is the singleton {Q 1 }.
It follows from the LaSalle Invariance Principle [12, Chapter 2, Theorem 6.4] that every solution to the equations in [START_REF] Beklemishev | Le cycle gonotrophique, principe de base de la biologie de an. gambiae, Vop Fiziol[END_REF] with initial conditions in R + 10 converge to the trivial equilibrium

Q 1 as t -→ +∞. Thus (x(t), X 2 (t)) -→ (0 R 9 , X *
2 ) when t -→ +∞ for N * ≤ 1. So, From the LaSalle principle we deduce the attractiveness of Q 1 . Since Q 1 is locally asymptotically stable when N * ≤ 1, we deduce that it is not only attractive, but it is also globally asymptotically stable. This ends the proof.

The non-trivial equilibria

In the following, we consider that the net reproductive number N * is greater than unity. We have the following result.

Theorem 2:

Let us define the following thresholds:

ξ = Λ 1 f r µ X 2 µ X 1 π (αk 3 + 2β 1 k 2 ) el(1 -ϑ)ηk 2 and N c = 1 + ξ. (9)
Assume that N * > 1.

1. If N * ≤ N c , model system (6) has no positive equilibrium point.

2. If N * > N c , model system (6) admits two positive equilibrium points

Q +/- 2 = E * +/-, L * +/-, P * +/-, X * 1+/-, X * 2 , Y * +/-, F * +/-, S * +/-, D * +/-, G * +/-where X * 2 = f r∧ 1 µ X 2 , P * = lL * +/- k 3 + β 1 L * +/- , Y * +/-= ϑη k 4 P * +/-, X * 1 +/-= (1 -ϑ)η µ X 1 P * +/-, F * +/-= β 2 eX * 1+/- k 5 (eX * 1+/-+ πX * 2 ) Y * +/-, S * +/-= ωk 7 k 8 k 9 F * +/-, D * +/-= γ 1 k 7 S * +/-, G * +/-= γ 2 k 8 D * +/-, E * +/-= µ b k 1 G * +/-, (10)
and L * +/-are positive solutions of the following equation

P(L * ) = L * ρ 3 L * 3 + ρ 2 L * 2 + ρ 1 L * + ρ 0 = 0, (11)
with

ρ 3 = -αβ 1 eηlϑµ X 2 -Λ 1 αβ 2 1 f µ X 1 πr < 0, ρ 2 = -[eηϑlµ X 2 (αk 3 + β 1 k 2 ) + Λ 1 β 1 f µ X 1 πr(2αk 3 + β 1 k 2 )] < 0, ρ 1 = el(1 -ϑ)ηµ X 2 k 2 k 3 [N * -N c ] and ρ 0 = -Λ 1 f k 2 k 2 3 µ X 1 πr < 0.
Proof. By setting the hand right sides of the equations of ( 6) to zero and expressing the variables in terms of L, we obtain that the non-trivial equilibria come from the resolution of (11) in terms of L. It follows, using Descartes' Rule of Signs, that the polynomial [START_REF] Gentile | Modelling sterile insect technique to control the population of anopheles gambiae[END_REF] has two positive non-trivial roots (L + 1 and L - 2 ) whenever N * > N c , and no positive root otherwise. Thus, we conclude that (6) admits two non-trivial positive equilibria (Q + 2 and Q - 2 ) whenever N * > N c . This ends the proof.

It is important to note that, in the presence of sterile males, the condition N * > 1 is not sufficient to the mosquitoes proliferation (see for example [START_REF] Abboubakar | Bifurcation thresholds and optimal control in transmission dynamics of arboviral diseases[END_REF][START_REF] Abboubakar | Backward bifurcation and control in transmission dynamics of arboviral diseases[END_REF][START_REF] Anguelov | Mathematical modeling of sterile insect technology for control of anopheles mosquito[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF][START_REF] Moulay | The chikungunya disease: modeling, vector and transmission global dynamics[END_REF] ). Indeed, it follows from the item (1) of Theorem 2, that if 1 < N * ≤ N c , there is no non-trivial equilibrium, which implies that there are no mosquitoes. This is illustrated in figure 3.

Note that in absence of sterile males, i.e Λ 1 = 0 or f = r = 0, equation ( 11) becomes

(L * ) 2 [ρ 2 L * 2 + ρ 1 L * + ρ 0 ] = 0, with ρ 2 = -αβ 1 eηlµ X 2 ϑ < 0, ρ 1 = -(αeηk 3 lµ X 2 ϑ + β 1 eηk 2 lµ X 2 ϑ) < 0, and ρ 0 = el(1 -ϑ)ηµ X 2 k 2 k 3 [N * -1]
. So, the system admits only one non-trivial equilibrium whenever N * > 1. Now, consider that N * > 1. Then, if the following condition holds

ρ 1 > 0 ⇐⇒ N c < N * ⇐⇒ Λ 1 < Λ crit 1 = e π µ X 2 f r l(1 -ϑ)η µ X 1 k 2 (αk 3 + 2β 1 k 2 ) (N * -1) , (12)
model system (6) admits two non-trivial equilibria Q + 2 and Q - 2 . Then, we can reformulate theorem 2 as follows.

Theorem 3: Assume that N * > 1. Then there exist Λ crit 1 > 0 such that model ( 6) admits two non-trivial equilibria if 0 < Λ 1 < Λ crit 1 and no non-trivial equilibrium otherwise.

III NUMERICAL SIMULATIONS

In this subsection, we perform some numerical simulations to illustrate results of theorem 3 and 2. To this aim, we use the following parameter values µ b = 250, s = 0.7, µ E = 0.2, µ L = 0.4, µ P = 0.4, l = 0.5, ω = 0.8, ϑ = 0.5, f = 0.7, r = 0.5, α = 0.07, η = 0.5, holds, i.e. when 1 < N * ≤ N c . Indeed, it clear that if there is this condition holds, then there are no non-trivial equilibria, which implies that there are no mosquitoes.

β 1 = 0.08, µ X 1 = 1/14, µ X 2 = 1/8, β 2 = 0.
The effects of varying parameter e to the decrease of the total number larvae are depicted on figure 4. It is clear that if female mosquitoes have a preference to mate with sterile males, the population of larvae will decrease significantly. Now, we numerically study the impact of sterile insect technique combined with the use of conventional insecticide on larvae and mosquitoes in blood meal seeking stage. Because Deltamethrin, for example, is effective only during a couple of hours and taking into account the preservation of the environment, it is not realistic to spray this chemical product continuously. So, we use "pulse control" technique which means that "the control is not continuous in time order is effective only one day every T day" [START_REF] Dumont | Vector control for the chikungunya disease[END_REF]. Also, the continue release of sterile males mosquitoes is not realistic. We thus consider that sterile males will be release after every one week (seven days) [START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF]. Using the above assumptions in numerical simulations, we obtain the following figures on which we clearly see the impact of the use of these two controls technique on larvae and mosquitoes in blood meal seeking stage.

From figure 5, it follows that increases in the number of sterile males permit to decrease in both larvae and female mosquitoes in the blood meal seeking. The same observation can be done if we use insecticide (adulticide). Indeed, in figure 6, it is clear that the periodic use of adulticide permits to decrease both larvae and female mosquitoes in the blood meal seeking. From figure 7, it is clear that the combination of these two controls permits to reduce significantly the total number of larvae and mosquitoes in the blood meal seeking stage.

=33 1 =90 1 =120 1 =164
Figure 5: Effect of release sterile male mosquitoes after every one week on larvae L and blood meal seeking mosquitoes S.

IV CONCLUSION

We analysed a complete life cycle model of mosquitoes, vectors of some diseases in human communities, which take into account the principal stage of the gonotrophic cycle of female mosquitoes. To control the population of eggs, larvae and female mosquitoes in blood meal seeking stage, we combined the sterile insect technique with the use of insecticides to kill female mosquitoes in blood meal seeking stage. To this end, we took into account a compartment of sterile male mosquitoes, and we modified the mortality rate of seeking mosquitoes. Unlike other published works in the literature which consider that sterile and wild mosquitoes have a same probability to mate with a wild female [START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF], we introduce one parameter e to traduce the fact that female mosquitoes, before mating, can make a choice between wild male and sterile male mosquitoes. We prove the global stability of the trivial equilibrium and showed that when the mosquito reproductive N * is greater than a certain threshold N c greater than one, the model admits two positive equilibria. This shows that, according to the quality and the quantity of the release of sterile males, the traditional condition N * > 1 is not sufficient to the proliferation of mosquito populations. So, the sterile insect technique has a great impact in the control of mosquito population.

To illustrate our analytical results, we used parameters of Aedes species and performed numerical simulations. Indeed, we have shown numerically existence of two positive equilibria whenever the mosquito reproduction number N * is greater than one and N c (1 < N c ≤ N * ).

On the other hand, if this condition is violated, then the model with sterile mosquitoes does not have a non-trivial equilibrium. It was clear that if female mosquitoes have a preference to mate with sterile males, the population of larvae will decrease significantly. Also, we used the technique of "pulse control" which contrast with the non-realistic continuous controls and concluded that the combination of sterile insect techniques and insecticides have a great impact on the decrease of total number of larvae and female mosquitoes in blood meal seeking stage, and so, can permit to decrease the total number of females mosquitoes.

A USEFUL DEFINITIONS AND RESULTS.

First of all, we provide a few preliminaries for our proof.

Theorem 4: p. 112 of [START_REF] Walter | Ordinary differential equations[END_REF] Assume that f : U -→ R n is cooperative, where U is open and convex, and that y, z : [t 0 , t 0 + a] -→ U are differentiable. If y(t 0 ) ≤ z(t 0 ), ẏ(t) ≤ f (y(t)), ż(t) = f (z(t))

∀t ∈ [t 0 , t 0 + a], then y ≤ z, on [t 0 , t 0 + a].

Theorem 5: [14, p. 62] Assume that ẏ = f (y) is cooperative in an open convex set U ⊆ R n and that y 0 and y 1 , y 0 < y 1 are the only equilibria in U. If the Jacobian matrix Df (y 0 ) is irreducible and s(Df (y 0 )) := max{Re λ : λ is an eigenvalue of A}, then there exists a unique solution y(t) (up to translation) satisfying ẏ(t) > 0, for all t ∈ R, y(t) → y 0 , as t → -∞, and y(t) → y 1 , as t → +∞.

1 =Figure 2 :

 12 Figure 2: Illustration of theorem 3 for 1 < N * = 61.4689 > N c = 13.144 and Λ 1 = 33 < Λ crit 1 = 164.31766. So, the positive equilibrium Q + 2 is locally stable.

Figure 3 :Figure 4 :

 34 Figure 3: Illustration of item (1) of theorem 3 for 1 < N * = 12.2938 ≤ N c = 15.6194. Parameter values are the same that those use in figure 2. So, model system (1) has no non-trivial equilibrium, which means that there are no mosquitoes where the condition 1 < N * ≤ N c holds.


Gang Pang
email: gangpang@buaa.edu.cn

Yibo Yang
email: ybyang@seas.upenn.edu

Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Shaoqiang Tang

Stability and convergence analysis of artificial boundary conditions for the Schrödinger equation on a rectangular domain

Keywords: Schrödinger equation, artificial boundary condition, stability analysis, convergence analysis AMS Subject Classification: 35J10, 65L20, 65L10, 65L12

published or not. The documents may come

Stability and convergence analysis of artificial boundary conditions for the Schrödinger equation on a rectangular domain

Introduction

Time-dependent Schrödinger equations play a key role in many applications related to physics and engineering [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Hazeltine | The Framework of Plasma Physics[END_REF][START_REF] Schrödinger | An undulatory theory of the mechanics of atoms and molecules[END_REF][START_REF] Senior | Optical Fiber Communications: Principles and Practice[END_REF], including quantum physics, acoustics, electromagnetism and optics. In many cases, the equation is set in the whole space R d (d ≥ 1). Therefore, when one wants to simulate the solution to the initial-value problem, we have to bound the computational domain to get a finite number of spatial degrees of freedom. Doing this, one needs to fix a suitable boundary condition at the fictitious interface of the domain that mimics the behavior of the wave field in the exterior domain. If the physical solution is confined within a finite domain [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF], trivially a homogeneous Dirichlet or Neumann boundary condition can be considered. Nevertheless, this is insufficient in some situations. A classical example is when the wave is outgoing, expanding then to infinity. More advanced boundary conditions must be designed, considering transparent, artificial or absorbing boundary conditions (ABC), according to the properties of the boundary condition. This usually requires a lot of mathematical and numerical analysis to state well-suited boundary conditions that minimize the spurious reflection at the boundary, but also to lead to accurate and stable numerical schemes. For the Schrödinger equation, and over the past two decades, many contributions were done for the oneand two-dimensional cases, for the linear or nonlinear Schrödinger equation, involving possibly some space and time variable potentials, with simple or complicate smooth convex fictitious boundaries. We refer to [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation on circular domains[END_REF][START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation: Fast calculation, approximation, and stability[END_REF][START_REF] Baskakov | Implementation of transparent boundaries for numerical solution of the Schrödinger equation[END_REF][START_REF] Fevens | Absorbing boundary conditions for the Schrödinger equation[END_REF][START_REF] Han | Exact artificial boundary conditions for Schrödinger equation in R 2[END_REF][START_REF] Li | Local artificial boundary conditions for Schrödinger and heat equations by using high-order azimuth derivatives on circular artificial boundary[END_REF] for some examples, the list of references being non exhaustive. Available review papers on artificial boundary conditions are e.g. [START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF][START_REF] Higdon | High-order local non-reflecting boundary conditions: a review[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF] for general PDEs, and [START_REF] Antoine | A Review of artificial boundary conditions for the Schrödinger equation[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] more specifically for Schrödinger-type problems. Many approaches can be developed to build these boundary conditions. Some methods are based on directly working with the continuous PDE, building approximations to the Dirichlet-to-Neumann (DtN) operators, others include constructions of boundary conditions for the semi-discretized equation (in time or space), and finally other techniques are directly applied at the full discrete level. For completeness, let us also mention that Perfectly Matched Layer (PML) [START_REF] Antoine | Perfectly matched layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF][START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF][START_REF]Perfectly matched layers for coupled nonlinear Schrödinger equations with mixed derivatives[END_REF][START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF][START_REF] Nissen | An optimized perfectly matched layer for the Schrödinger Equation[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] can be used for truncating the computational domain.

Here, we consider the time-dependent linear free Schrödinger equation truncated in a rectangular two-dimensional domain. When the fictitious boundary has some corners, like this is the case for the rectangular domain, then the construction of artificial boundary conditions is known to be more complicate since the existence of corners directly impacts the design of the boundary conditions. For example, specific nontrivial compatibility boundary conditions at the corners need to be added to the artificial boundary conditions to get an increased accuracy for time-domain wave equations [START_REF] Bamberger | Second-order absorbing boundary conditions for the wave equation: a solution for the corner problem[END_REF][START_REF] Collino | Conditions absorbantes d'ordre élevé pour les équations de Maxwell dans des domaines rectangulaires[END_REF][START_REF] Hagstrom | A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first-order systems[END_REF][START_REF] Modave | A GPU-accelerated nodal discontinuous Galerkin method with high-order absorbing boundary conditions and corner/edge compatibility[END_REF] and time-harmonic wave scattering [START_REF] Modave | Corner treatments for high-order local absorbing boundary conditions in high-frequency acoustic scattering[END_REF]. In addition, the numerical analysis of such problems remains difficult to conduct. For the PML problem, corners may lead to instabilities in the time-domain [START_REF] Bécache | Remarks on the stability of Cartesian PMLs in corners[END_REF] and require a careful choice of the discretization to stabilize the method. In addition, tuning correctly the PML parameters is not so easy and strongly depends on the configuration, i.e. the PDE under study and its discretization [START_REF] Antoine | Perfectly matched layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF][START_REF]Perfectly matched layers for coupled nonlinear Schrödinger equations with mixed derivatives[END_REF][START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF][START_REF] Nissen | An optimized perfectly matched layer for the Schrödinger Equation[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF].

In the present paper, we focus on the artificial boundary condition method proposed in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] which solves exactly the corner problem for the rectangular shaped domain, and is based on a secondorder spatial semi-discretization of the Schrödinger equation. Therefore, the corresponding boundary condition falls into the framework of semi-discrete artificial boundary conditions. From some numerical simulations, it was shown in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] that, using a Crank-Nicolson scheme, the full approximation is second-order in time. The aim of the present paper is to prove the stability of the fully discrete scheme as well as the second-order error estimate in time, for fixed spatial discretization steps ∆x and ∆y. However, the paper does not address the important and difficult question of designing an accelerated numerical method for evaluating the integral operators arising in the definition of these ABCs. This needs further numerical and computational studies which are out of the scope of the present paper. Since this paper is already technical as well as [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF], the authors strongly advise the reader to first understand the derivation in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] before entering into the details of the present paper, even if the material to understand the paper is self-contained. Let us also remark that the method for constructing ABCs [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] has recently been applied with success in [START_REF] Ji | Artificial boundary conditions for the semi-discretized onedimensional nonlocal Schrödinger equation[END_REF] to design new artificial boundary conditions for the 1D nonlocal Schrödinger equation. This later work, combined with the design of ABCs on rectangular domains developed in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF], should offer the possibility to extend the ABCs to the 2D nonlocal Schrödinger equation (as well as heat equation) and to some peridynamics models. Therefore, the numerical analysis developed in the present paper for the 2D local Schrödinger equation is valuable for understanding how to obtain the convergence and stability results for the more complicate nonlocal case.

The plan of the paper is the following. After the introduction (Section 1), we develop in Section 2 the results related to the exact artificial boundary condition for the discrete spatial rectangular lattice. We also proceed to the full discretization of the initial boundary-value problem. In Section 3, we prove the L 2 -stability of the fully discrete scheme. Section 4 is devoted to the convergence analysis. Finally, we end by a conclusion in Section 5 and two Appendices.

2 Exact boundary conditions for the semi-discrete Schrödinger equation on a rectangular domain

Let us consider the two-dimensional linear Schrödinger equation with unknown ψ := ψ(x, y, t)

ι ψ = - ∂ 2 ψ ∂x 2 + ∂ 2 ψ ∂y 2 , (1)
where (x, y) ∈ R 2 and t > 0 are the space and time variables, respectively. The time derivative is denoted by ψ := ∂ t ψ. In addition, a compactly supported initial data ψ(x, y, t = 0) = ψ 0 (x, y) is added to equation (1). Finally, the complex-valued number ι is such that ι := √ -1. Let us now introduce the uniform spatial steps ∆x and ∆y in the x-and y-directions, respectively. By discretizing (1) with a second-order five points stencil scheme for the Laplacian (assuming here that ∆x = ∆y) and after rescaling the time variable through t = (∆x 2 ) -1 t, we obtain the following semi-discrete system for the exact (e) solution ψ e := (ψ e m,n

) (m,n)∈Z 2 to ι ψe m,n = -(ψ e m+1,n + ψ e m-1,n + ψ e m,n+1 + ψ e m,n-1 -4ψ e m,n ) = -∆ψ e m,n , (2)
with compactly supported initial data ψ e m,n (0) = ψ 0 (m∆x, n∆x) for the uniform lattice (m, n) ∈ Z 2 . In the above equation, we define the discrete Laplacian as: ∆ := ∆ x + ∆ y , with ∆ x,y := ∆ +

x,y ∆ - x,y and we denote by ∆ +

x (∆ - x , respectively) the forward (backward, respectively) finite-difference operator in the x-direction, i.e.

∆ +

x ψ e m,n := ψ e m+1,n -ψ e m,n , ∆ - x ψ e m,n := ψ e m,n -ψ e m-1,n .

Similarly, the operators ∆ + y and ∆ - y are given in the y-direction by ∆ + y ψ e m,n := ψ e m,n+1 -ψ e m,n , ∆ - y ψ e m,n := ψ e m,n -ψ e m,n-1 .

Let us define the Laplace transform h of h as

h(s) := ∞ 0 h(t)e -st dt (3)
and the convolution (f * g) of f and g as

(f * g)(t) := t 0 f (τ )g(t -τ )dτ. (4)
Therefore, the Laplace transform ψe m,n (s) of ψ e m,n is given by ψe m,n (s) := ∞ 0 ψ e m,n (t)e -st dt.

Before presenting the construction of boundary conditions for general compact initial data, we start with a special case, i.e., the so-called point-source problem to compute the discrete Green's functions which is needed for the general case. We assume that the wave function ψ e satisfies Eq. ( 2) and that ψ e m,n (0) = 0, for (m, n) = (0, 0). For a given function ψ e 0,0 (t), from [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] we have, for (m, n) = (0, 0),

ψ e m,n (t) = f e m,n (t) * ψ e 0,0 (t) = t 0 f e m,n ( t)ψ e 0,0 (t -t)d t, (6)
where the Green's function f e m,n (t) comes from the expression

g e m,n (t) = f e m,n (t) * g e 0,0 (t), (7)
with

g e m,n (t) = (ιt) m+n e -4ιt (m + n)! m!n! √ π 2 m+n+1 Γ( m + n + 1 2 )Γ( m + n 2 + 1) × 2 F 3 1 2 (m + n + 1), 1 2 (m + n) + 1; -4t 2 m + 1, n + 1, m + n + 1 . (8)
In the above definition, Γ is the Gamma function and the generalized hypergeometric function 2 F 3 [START_REF] Lebedev | Special Functions and Their Applications[END_REF] is given by

2 F 3 a 1 , a 2 ; z b 1 , b 2 , b 3 = ∞ k=0 (a 1 ) k (a 2 ) k k!(b 1 ) k (b 2 ) k (b 3 ) k z k , (9)
W ((m, n), z) = 1 2 m+n z m+n+1 (m + n)! m!n! × 4 F 3 1 2 (m + n + 1), 1 2 (m + n + 1), 1 2 (m + n) + 1, 1 2 (m + n) + 1; 4 z 2 m + 1, n + 1, m + n + 1 , (10)
and z = (4 -ιs)/2. The function 4 F 3 is the generalized hypergeometric function p F q [START_REF] Lebedev | Special Functions and Their Applications[END_REF], with p = 4 and q = 3. In particular, we have f e 0,0 (t) = δ(t), where δ is the Dirac distribution. Let us now come to the construction of the exact boundary conditions for general initial data. We consider a finite rectangular computational domain and we assume that the initial data is zero outside the selected rectangular box, including its boundary. To simplify the presentation, we suppose that the domain is the square L × L (see Fig. 1). The nodes at the boundary are labelled counterclockwise, with the interior layer set as I k , for 1 ≤ k ≤ 4L -4. More explicitly, we have

I = (1, ) for 1 ≤ ≤ L -1, I +L-1 = ( , L) for 1 ≤ ≤ L -1, I +2L-2 = (L, L + 1 -) for 1 ≤ ≤ L -1 and I +3L-3 = (L + 1 -, 1) for 1 ≤ ≤ L -1.
The outer layer points are denoted by J k for k = 1, • • • , 4L, with the same orientation. Deriving a numerical artificial boundary condition consists in writing ψ e J k (t) as a function of the values of the wave field in the interior layer, completing hence the semi-discrete system (2) to get a problem set in the finite computational region. : Discrete computational lattice (for L = 4): the red bullets represent the outer layer points, the black ones correspond to the interior layer points, the yellow ones are the interior points and the green ones stands for the corner points.

From [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] and for any integers (m, n), the functions ψ e m,n (t) are proved to be analytic with respect to t, and ψe m,n (s) is analytic around ∞. Let us consider that ψ e I k (t), for k = 1, • • • , 4L -4, are given values. Following [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF], we first decompose them as 4L -4 unknown independent sources χ e I k (t) thanks to the previous Green's functions f e m,n following

ψ e I k (t) = 4L-4 =1 f e I k -I (t) * χ e I (t). (11)
For 1 ≤ k ≤ 4L -4 and z = 4 -ιs 2 , it can be proved [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF] that the following relation holds in the Laplace domain ψe

I k (s) = 4L-4 =1 W (I k -I , z) W ((0, 0), z) χe I (s), (12)
where χe I (s) denotes the Laplace transform of the function χ e I (t). From [START_REF] Bécache | Remarks on the stability of Cartesian PMLs in corners[END_REF] and [START_REF] Pang | Eliminating corner effects in square lattice simulation[END_REF], one has the asymptotic estimates [START_REF] Lebedev | Special Functions and Their Applications[END_REF]

W ((m, n), z) = W ((|m|, |n|), z)
W ((m, n), z) W ((0, 0), z) ∼ C |m| |m|+|n| ( 1 2z ) |m|+|n| , |z| → +∞,
where

C
|m| |m|+|n| is the binomial coefficient. Therefore, the matrix (

W (I k -I , z) W ((0, 0), z) ) k, is diagonally
dominant for large values of the parameter s. Consequently, the matrix is invertible and the functions χI (s) exist and can be determined uniquely. Since ψe I k (s) are analytic near the infinity, then the functions χe I (s) also share this property and χ e I (t) are analytical functions. It can also be proved that χ e I (0) = 0. For a point (m, n) outside the inner boundary, we have the representation formula

ψe m,n (s) = 4L-4 =1 W ((m, n) -I , z) W ((0, 0), z) χe I (s). ( 13
)
For any inner boundary point I , a direct computation shows that

(4 -ιs)W ((m, n) -I , z) = W ((m + 1, n) -I , z) + W ((m -1, n) -I , z) +W ((m, n + 1) -I , z) + W ((m, n -1) -I , z) . ( 14
)
Injecting the above relation into [START_REF] Collino | Conditions absorbantes d'ordre élevé pour les équations de Maxwell dans des domaines rectangulaires[END_REF] gives

(4 -ιs) ψe m,n (s) = ψe m+1,n (s) + ψe m-1,n (s) + ψe m,n+1 (s) + ψe m,n-1 (s), (15)
meaning that for any point (m, n) outside the inner boundary, ψ e m,n (t) satisfies [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. It can also be checked that ψ e m,n (0) = 0. Once the sources χ e I (t) have been obtained, an exact boundary condition reads as follows, for 1 ≤ ≤ 4L,

ψ e I (t) = 4L-4 l=1 f e I -I l (t) * χ e I l (t), ψ e J (t) = 4L-4 l=1 f e J -I l (t) * χ e I l (t). (16)
For two sequences u = {u k } and v = {v k }, let us now define the discrete convolution

(u v) k = k l=0 u l v k-l ∆t l := k l=0 u l v k-l ∆t -u 0 v k ∆t 2 -u k v 0 ∆t 2
through the trapezoidal quadrature rule, with ∆t l = ∆t/2 for l = 0, k and ∆t l = ∆t otherwise. For a fixed time step ∆t, (f e m,n ) k , (ψ e m,n ) k and (χ e m,n ) k stand for the quantities f e m,n (k∆t), ψ e m,n (k∆t) and χ e m,n (k∆t), respectively. In addition, the numerical approximation of f e m,n (k∆t), ψ e m,n (k∆t),

χ e m,n (k∆t) is designated by f k m,n , ψ k m,n , χ k m,n
, respectively. In the present paper, the unknowns ψ k m,n

(1 ≤ m, n ≤ L) are computed by applying the Crank-Nicolson scheme [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] to (2)

-ι ψ k m,n -ψ k-1 m,n ∆t = 1 2 ∆(ψ k m,n + ψ k-1 m,n ) = 1 2 (ψ k m-1,n + ψ k m+1,n + ψ k m,n-1 + ψ k m,n+1 -4ψ k m,n ) + 1 2 (ψ k-1 m-1,n + ψ k-1 m+1,n + ψ k-1 m,n-1 + ψ k-1 m,n+1 -4ψ k-1 m,n ), (17)
for 1 ≤ m, n ≤ L. Concerning the boundary condition, ψ k I and ψ k J (1 ≤ ≤ 4L) are discretized from relations [START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF] as

ψ k I = 4L-4 l=1 (f I -I l χ I l ) k , ψ k J = 4L-4 l=1 (f J -I l χ I l ) k . ( 18
)
Based on [START_REF] Bamberger | Second-order absorbing boundary conditions for the wave equation: a solution for the corner problem[END_REF], the terms f k m,n in relation ( 18) can be numerically computed from

(g e m,n ) k = (f m,n g e 0,0 ) k , (19)
where g e m,n (t) can be approximated accurately e.g. by a fourth-order Runge-Kutta scheme. Finally, we can use the values ψ k I to complete the fully discretized system given by ( 17), for 1 ≤ m, n ≤ L, and then prove the stability of the whole system. Remark 1. If the external potential outside the computational domain is a constant V but not zero, then only the kernel function ge m,n (s) needs to be changed to W ((m, n),

4 + V -ιs 2
), and the resulting boundary condition (16) still holds.

3 L 2 -stability analysis of the scheme

In this section, we prove the stability of the scheme ( 17)- [START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF]. Le us first introduce the norm

ψ 2 = ( L m,n=1 |ψ m,n | 2 ) 1/2 (20)
and the Z-transform û of a sequence u := {u k } as

û(z) = ∞ k=0 u k z k . ( 21
)
If

w k = (u v) k , one directly gets ŵ(z) = ∞ k=0 w k z k = ∞ k=0 1 z k k l=0 u l v k-l ∆t l = û(z)v(z)∆t -u 0 ∆t 2 v(z) -v 0 ∆t 2 û(z). ( 22
)
From the relation (g e m,n ) k = (f m,n g e 0,0 ) k and ( 22), one has ĝe m,n (z) = fm,n (z)ĝ e 0,0 (z)∆t -ĝe 0,0 (z)f e m,n (0

) ∆t 2 -f e m,n (z)g e 0,0 (0) ∆t 2 , leading to fm,n (z) = ĝe m,n (z) + ĝe 0,0 (z)f e m,n (0) ∆t 2
(ĝ e 0,0 (z) -

g e 0,0 (0) 2 )∆t , (23)
where ĝe

m,n (z) = ∞ k=0 g e m,n (k∆t) z k . If v k = (u f m,n ) k , with u 0 = 0, we have v(z) = fm,n (z)û(z)∆t -f e m,n (0)û(z) ∆t 2 = ĝe m,n (z) + f e m,n (0)∆t 4
(ĝ e 0,0 (z) -

g e 0,0 (0) 2 ) û(z) = h m,n (z)û(z). ( 24
)
From the definition of the terms g e m,n (t), the quantities h m,n (z) can be uniformly written as

h m,n (z) = ĝe m,n (z) - 1 2 δ 0,0 + i 4 tδ 0,1 ĝe 0,0 (z) - 1 2 . ( 25
)
Here, δ 0,0 = 1 for (m, n) = (0, 0), δ 0,0 = 0 for the other cases. Similarly, we set δ 0,1 = 1 for (m, n) = (0, 1), and δ 0,1 = 0 otherwise. Let us now introduce the index I 1 m , 1 ≤ m ≤ 4L, such that:

I 1 m = (1, m), I 1 L+m = (m, L), I 1 m+2L = (L, L + 1 -m) and I 1 m+3L = (L + 1 -m, 1)
, for 1 ≤ m ≤ L. Therefore, I 1 and J have both 4L elements. Let us set β = e ρ∆t > 1, with ρ > 0. We use this term below to act like a phase change in the spectral space of Z-transformation. Doing this, we have some integrals later (e.g. like for [START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF] or (41)) which are not singular, simplifying the analysis. Taking

ϕ k m,n = 1 β k ψ k m,n in (17), we obtain -ι 2 ∆t (ϕ k m,n -ϕ k-1 m,n ) = ∆(ϕ k m,n + ϕ k-1 m,n ) + (β -1)(∆ + ι 2 ∆t )ϕ k m,n (26)
and

-ι 2 ∆t β(ϕ k m,n -ϕ k-1 m,n ) = β∆(ϕ k m,n + ϕ k-1 m,n ) + (1 -β)(∆ -ι 2 ∆t )ϕ k-1 m,n . (27)
Let us now multiply the conjugate of ( 26) by ϕ k m,n and ( 27) by -β -1 ϕ k-1 m,n (z being the conjugate of a complex number z), and sum over 1 ≤ m, n ≤ L. Then one gets

L n=1 L m=1 (|ϕ k m,n | 2 -|ϕ k-1 m,n | 2 ) ≤ -(1 - 1 β ) L m=1 L n=1 |ϕ k-1 m,n | 2 -(β -1) L m=1 L n=1 |ϕ k n | 2 + ∆t 2β Im L n=1 (ϕ k-1 0,n + βϕ k 0,n ) + x (ϕ k-1 0,n + βϕ k 0,n ) -(ϕ k-1 L+1,n + βϕ k L+1,n ) - x (ϕ k-1 L+1,n + βϕ k L+1,n ) + ∆t 2β Im L m=1 (ϕ k-1 m,0 + βϕ k m,0 ) + y (ϕ k-1 m,0 + βϕ k m,0 ) -(ϕ k-1 m,L+1 + βϕ k m,L+1 ) - y (ϕ k-1 m,L+1 + βϕ k m,L+1
) .

(28) Summing up all these inequalities from k = 1 to K and summing by parts yield

ϕ K 2 2 ≤ ϕ 0 2 2 + ∆t 2β Im 4L m=1 K k=1 (ϕ k-1 Jm + βϕ k Jm )(ϕ k-1 I 1 m + βϕ k I 1 m -ϕ k-1 Jm -βϕ k Jm ) , (29)
where k = K is the index corresponding to the time of computation t K = K∆t.

Let us now determine the sign of the second term appearing in the right-hand side of the above inequality.

Lemma 1. Let us consider the scheme ( 17)- [START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF]. Let ρ be a (relatively large) given positive constant and let us set β = e ρ∆t > 1. Then, there exists a time step ∆t 0 such that, for ∆t ≤ ∆t 0 , the following inequality holds

4L m=1 K k=1 (ϕ k-1 Jm + βϕ k Jm )(ϕ k-1 I 1 m + βϕ k I 1 m -ϕ k-1 Jm -βϕ k Jm ) ≤ 0. ( 30
)
Proof. For h m,n defined by [START_REF] Lebedev | Special Functions and Their Applications[END_REF], we can introduce the three following matrices A, A 1 and B with elements

A ,l (z) = h I 1 -I l (z), A 1 ,l (z) = h I -I l (z) and B ,l (z) = h J -I l (z). ( 31
)
The matrices A and B are of size 4L × (4L -4) while A 1 is a (4L -4) × (4L -4) matrix. Let us now define the vectors

ψI 1 (z) = [ ψI 1 1 , ψI 1 2 , ... ψI 1 4L ] T , ψJ (z) = [ ψJ 1 , ψJ 2 , ... ψJ 4L ] T , χI (z) = [ χI 1 , χI 2 , ... χI 4L-4 ] T , (32)
where T denotes the transposition operation. Therefore, the Z-transform of the boundary conditions can be written in the compact form

ψI 1 (z) = A(z) χI (z), ψI (z) = A 1 (z) χI (z), ψJ (z) = B(z) χI (z). ( 33
)
Following these notations, we obtain

4L l=1 K k=1 Im (ϕ k-1 J l + βϕ k J l )(ϕ k-1 I 1 l + βϕ k I 1 l -ϕ k-1 J l -βϕ k J l ) = 1 2π Im 2π 0 |1 + βe ιφ | 2 ψT J (βe ιφ )( ψI 1 (βe ιφ ) -ψJ (βe ιφ ))dφ = 1 2π Im 2π 0 |1 + βe ιφ | 2 χT I (βe ιφ ) BT (βe ιφ )A(βe ιφ ) χI (βe ιφ )dφ . (34)
To conclude, we only need to prove that

Im VT BT (βe ιφ )A(βe ιφ )V ≤ 0, (35)
for any φ ∈ [0, 2π] and V ∈ C 4L-4 . The analysis is stated by considering three cases for φ.

By definition [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF], we have e (K 0 (x,y)+ρ) t+ιφ e (K 0 (x,y)+ρ) t+ιφ -1 dxdy,

g e m,n (t)
with K 0 (x, y) := 2ι(2 -cos(x) -cos(y)). We set α = ρ 3 and z = βe ιφ . Case 1:

α∆t ≤ φ ≤ π -α∆t or π + α∆t ≤ φ ≤ 2π -α∆t.
Let us first write the following expansions ĝe 0,0 (z) -

1 2 A = A 0 + A 1 ∆t + R 1 ∆t 2 , ḡe 0,0 (z) - 1 2 BT = B1 T ∆t + B2 T ∆t 2 + R 2 ∆t 3 . ( 38
)
The four matrices A 0 , A 1 , B1

T and B2 T are computed explicitly below. The two residual terms R 1 and R 2 then must be estimated to be controlled. Therefore, one gets ĝe 0,0 (z) -

1 2 2 BT A = B1 T A 0 ∆t + (B T 1 A 1 + B T 2 A 0 )∆t 2 +(R 2 A 0 + B T 1 R 1 + B T 2 A 1 )∆t 3 + (B T 2 R 1 + R 2 A 1 )∆t 4 + R 2 R 1 ∆t 5 . (39)
Now, because the coefficients of A and B are given through h m,n (z) by expressions [START_REF] Pang | Eliminating corner effects in square lattice simulation[END_REF] which are related to ĝe m,n (z) by [START_REF] Lebedev | Special Functions and Their Applications[END_REF], we must derive an expansion of ĝe m,n (βe ιφ ) thanks to ∆t. To this end, let us compute the Taylor's expansion

e (K 0 +ρ) t+ιφ e (K 0 +ρ) t+ιφ -1 = - ι 2 cot( φ 2 ) + 1 2 + (K 0 + ρ) ∆t 4 sin 2 ( φ 2 ) +(K 0 + ρ) 2 ι cot( φ 2 ) 8 sin 2 ( φ 2 ) ∆t 2 -(K 0 + ρ) 3 cos(φ) + 2 48 sin 4 ( φ 2 ) ∆t 3 (1 + O( ρ α )). (40)
Plugging ( 40) into (37) needs the computation of the following integrals

1 (2π) 2 2π 0 2π 0 cos(mx) cos(ny)(K 0 (x, y) + ρ)dxdy = (ρ + 4ι)δ 0,0 -ιδ 1,0 -ιδ 0,1 (41) and 1
(2π) 2 2π 0 2π 0 cos(mx) cos(ny)(K 0 (x, y) + ρ) 2 dxdy = (ρ 2 + 8ιρ -20)δ 0,0 -δ 2,0 -δ 0,2 + (8 -2ιρ)δ 0,1 + (8 -2ιρ)δ 1,0 -2δ 1,1 . (42)
The matrices A 0 , A 1 , B 1 , B 2 are given from expressions (105) to (109) in Appendix B. Let us analyze the expansion (39) and its effect on stating the inequality [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] for Case 1.

First, since B1 T A 0 is given by a real-valued diagonal matrix, it does not contribute to [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF]. Next, the diagonal elements of (B

T 1 A 1 + B T 2 A 0 ) are given by - ρ cot 2 ( φ 2 ) 8 sin 2 ( φ 2 ) ι - 3 2 cot 2 ( φ 2 ) sin 2 ( φ 2 )
or -

ρ cot 2 ( φ 2 ) 16 sin 2 ( φ 2 ) ι - 3 4 cot 2 ( φ 2 ) sin 2 ( φ 2 )
, the off-diagonal elements consisting in terms cot 4 (φ/2) 4 and cot 2 (φ/2) sin -2 (φ/2) multiplied by realvalued constants. Therefore, only the imaginary part of the diagonal elements of ( 3 , an analysis shows that its coefficients satisfy cos( φ 2 )

B T 1 A 1 + B T 2 A 0 ) contribute to the inequality (35). Concerning (R 2 A 0 + B T 1 R 1 + B T 2 A 1 )∆t
sin 5 ( φ 2 ) O(ρ 3 )∆t 3 = cot 2 ( φ 2 ) sin 2 ( φ 2 ) O(ρ 3 ) ∆t 3 sin(φ) ≤ cot 2 ( φ 2 ) sin 2 ( φ 2 ) O(ρ 3 ) α ∆t 2 . ( 43
)
Finally, the order of the coefficients of

| (B T 2 R 1 + R T 2 A 1 )∆t 4 | is at most cot 2 ( φ 2 ) sin 2 ( φ 2 ) O(ρ 4 ) α 2 ∆t 2
while the order of the elements of

| R T 2 R 1 ∆t 5 | is cot 2 ( φ 2 ) sin 2 ( φ 2 ) O(ρ 5 ) α 3 ∆t 2 .
From the above discussion, if ρ is large enough and ∆t small enough, we have, for any vector V ∈ C 4L-4 , the following inequalities ĝe 0,0 (z) -

1 2 2 Im VT BT AV ≤ - ρ 16 cot 2 ( φ 2 ) sin 2 ( φ 2 ) ∆t 2 |V| 2 + cot 2 ( φ 2 ) sin 2 ( φ 2 ) ∆t 2 O( ρ 3 α ) + O( ρ 4 α 2 ) + O( ρ 5 α 3 ) |V| 2 ≤ 0. ( 44
)
Case 2: -α∆t ≤ φ ≤ α∆t. Any angle φ can be written as φ = φ ∆t ∆t, where -α ≤ φ ∆t ≤ α. When ∆t → 0, we can write Let us recall that ρ is assumed to be relatively large. From (10), we have [START_REF] Lebedev | Special Functions and Their Applications[END_REF] W

((k, l), z) ∼ 1 2 k+l z k+l+1 C k k+l (1 + O( 1 z ))
when |z| → +∞. Therefore, for I 1 m -I n = (k, ) and ρ large, we obtain

W (I 1 m -I n , ρ + ιφ ∆t 2ι + 2) ∼ 1 2 k+ C k k+ ( 2ι ρ + ι(φ ∆t + 4) ) k+ +1 1 + O 1 ρ + ι(φ ∆t + 4)
.

From (45) and ∆t small enough, we have ĝe

I 1 m -In (z)∆t - 1 2 δ 0,0 ∆t + ι 4 t 2 δ 0,1 = 1 (2π) 2 2π 0 2π 0 cos(kx) cos(ly) ρ + ιφ ∆t + K 0 (x, y) dxdy + O(∆t 2 ) = 1 2ι W ((k, l), ρ + ιφ ∆t 2ι + 2) + O(∆t 2 ). ( 46
)
Let us define I as the (4L -4) × (4L -4) unitary matrix. Then we have ∆t ĝe 0,0 (z) -

1 2 A = 1 ρ + ι(φ ∆t + 4) I 1 (I + O( 1 ρ + ι(φ ∆t + 4) )) + O(∆t 2 ) ( 47
)
and ∆t ĝe 0,0 (z) -

1 2 B = ι (ρ + ι(φ ∆t + 4)) 2 I 1 (I + O( 1 ρ + ι(φ ∆t + 4) )) + O(∆t 2 ). ( 48
)
Since I T 1 I 1 = I 2 , where I 2 is given by (110) (see Appendix B), then one gets ∆t 2 ĝe 0,0 (z) -

1 2 2 BT A = -ιρ + (4 + φ ∆t ) ((φ ∆t + 4) 2 + ρ 2 ) 2 (I 2 + O( 1 ρ + ι(φ ∆t + 4) )) + O(∆t 2 ), (49)
and finally Im VT BT AV ≤ 0. Case 3: π -α∆t ≤ φ ≤ π + α∆t. Let us assume that: φ = π + ∆tφ ∆t . If we take

K 1 = K 0 + ρ + ιφ ∆t and K 2 = ρ + ι(4 + φ ∆t ), then we can write ĝe m,n (e ρ∆t+φι ) = 1 (2π) 2 2π 0 2π 0 cos(mx) cos(ny) ∞ j=0 e -K 0 j t /(e (ρ+ιφ ∆t )∆t ) j (-1) j dxdy = 1 (2π) 2 2π 0 2π 0 cos(mx) cos(ny) 1 1 + e -K 1 ∆t dxdy = 1 2 δ 0,0 + ( K 2 4 )δ 0,0 ∆t - ι 4 δ 0,1 ∆t - ι 4 δ 1,0 ∆t -(K 3 2 δ 0,0 -3K 2 δ 0,2 -3K 2 δ 2,0 + ιδ 3,0 + ιδ 0,3 + ι(9 -3K 2 2 )δ 1,0 +ι(9 -3K 2 2 )δ 0,1 + 3ιδ 1,2 + 3ιδ 2,1 -6K 2 δ 1,1 -12K 2 ) ∆t 3 48 + O(∆t 4 ), (50)
which means that ĝe 0,0 (z) -

1 2 A = ∆t K 2 4 I 1 + O(∆t 2 ) ( 51
)
and ĝe 0,0 (z) -

1 2 B = ∆t 3 ιK 2 2 16 I 1 + O(1) + O(∆t 4 ). (52)
This leads to ĝe 0,0 (z) -

1 2 2 BT A = -∆t 4 (ρ 2 + (4 + φ ∆t ) 2 ) ιK 2 64 I 2 + O( 1 ρ ) + O(∆t 5 ). ( 53
)
In this case, for any small enough time step ∆t, we have : Im VT BT AV ≤ 0.

Collecting the results from the three above cases, we conclude that

Im 4L l=1 K k=1 (ϕ k-1 J l + βϕ k J l )(ϕ k-1 I 1 l + β 1 ϕ k I 1 l -ϕ k-1 J l -β 1 ϕ k J l ) ≤ 0. ( 54
)
By using the previous inequality, we can state the following L 2 -stability result.

Theorem 1. For the scheme ( 17)-( 18), if ρ is a given and large enough positive constant, then there exists ∆t 0 such that, if ∆t ≤ ∆t 0 and

T K = K∆t, ψ K 2 ≤ e ρT K ψ 0 2 . ( 55
)
where r m,n (t, ∆t) is an analytic function such that there exists a constant C m,n,ρ with rm,n (z, ∆t) ĝe 0,0 (z) -

1 2 ≤ C m,n,ρ . (60)
We also have

|h m,n (z)| ≤ C m,n,ρ . (61)
Proof. By using the representation given by (56), we have

g e m,n (k∆t) = (f e m,n * g e 0,0 )(k∆t) = (f e m,n g e 0,0 ) k + r m,n (k∆t, ∆t)∆t 2 ,
which leads to

f e m,n (z) = ĝe m,n (z) + ĝe 0,0 (z)f e m,n (0) ∆t 2
(ĝ e 0,0 (z) -

1 2 )∆t + rm,n (z, ∆t) ĝe 0,0 (z) - 1 2
∆t which is nothing else than (59) since we have the representation [START_REF] Ji | Artificial boundary conditions for the semi-discretized onedimensional nonlocal Schrödinger equation[END_REF] for fm,n (z). Now, let u(t) be a real analytic function such that its Laplace transform ũ(s) is analytic at s = ∞ and can be written as ũ(s) = ∞ k=1 u k s -k . Let us remark here that g m,n (t) and r m,n (t, ∆t) fullfil these properties. Let us consider η and ρ as two large enough real-valued constants, with ρ > η, and z = e ρ∆t+ιφ . Then, we have

û(z) - u(0) 2 = 1 2πι η+ι∞ η-ι∞ ∞ k=0 ũ(s) e sk∆t z k ds - u(0) 2 = 1 4πι η+ι∞ η-ι∞ ∞ k=0 u k s k+1
1 + e -(ρ-s)∆t-ιφ

1 -e -(ρ-s)∆t-ιφ ds

= 1 2 ∞ k=0 1 k! d k ds k 1 + e s-ρ∆t-ιφ 1 -e s-ρ∆t-ιφ s=0 u k ∆t k . ( 62
)
Some computations show that, for any l ≥ 1, we have

∆t l 2 1 l! d l ds l 1 + e s-ρ∆t-ιφ 1 -e s-ρ∆t-ιφ s=0 = ∆t l 2 l m=1 (e -ρ∆t-ιφ ) m (1 -e -ρ∆t-ιφ ) m+1 |α|=l,α i ≥1 1 α! , ( 63
)
setting |α| := m i=1 α i and α! := α 1 !...α m !, for a multi-index α := (α 1 , ..., α m ). In addition, one gets

|α|=l,α i ≥1 1 α! ≤ |α|=l 1 α! = 1 2πι |z|=1 1 z l+1 (1 + z + z 2 2! + ...) m dz = 1 2πι |z|=1 1 z l+1 e mz dz = 1 2π 2π 0 e me ιθ -ι(l-1)θ dθ ≤ e m . ( 64
)
Case 1. π -α∆t ≤ φ ≤ π + α∆t. Let us introduce φ = π + φ ∆t ∆t and take α = ρ 3 as in the previous section. For ∆t small enough and from the inequalities (63)-(64), we have

∆t l 2 1 l! d l ds l 1 + e s-ρ∆t-ιφ 1 -e s-ρ∆t-ιφ s=0 ≤ ∆t l 2 l m=1 (e -ρ∆t ) m |1 + e -ρ∆t-ιφ ∆t ∆t | m+1 e m ≤ C l (e∆t) l . ( 65
)
Now, taking û(z) := rm,n (z, ∆t) in (62) (with z = e ρ∆t+ιφ ), using ( 63)-( 64) and noticing that r m,n (0, ∆t) = 0 (by setting t = 0 in (57)) leads to

|r m,n (z, ∆t)| ≤ ∞ k=1 C k (e∆t) k r m,n (k∆t, ∆t) = O(∆t). (66)
In addition, a direct computation yields ĝe 0,0 (e ρ∆t+ιφ ) -

1 2 = 1 2
1 -e -ρ∆t-ιφ ∆t ∆t 1 + e -ρ∆t-ιφ ∆t ∆t -

2e -ρ∆t-ιφ ∆t ∆t (1 + e -ρ∆t-ιφ ∆t ∆t ) 2 ∆t + O(∆t 2 ) = ρ + ιφ ∆t 4 -1 ∆t + O(∆t) 2 , ( 67
)
giving then the inequality (60). Case 2. -α∆t ≤ φ ≤ α∆t. We again set: φ = ∆tφ ∆t . For small enough ∆t and by using ( 63)-( 64), we obtain

∆t l 2 1 l! d l ds l 1 + e s-ρ∆t-ιφ 1 -e s-ρ∆t-ιφ s=0 ≤ ∆t l 2 l m=1 (e -ρ∆t ) m (1 -e -ρ∆t-ιφ ) m+1 |α|=l,α i ≥1 1 α! ≤ C l 1 |1 -e -ρ∆t-ιφ ∆t ∆t | l m=1 ( e |ρ + φ ∆t | ) m ∆t l-m ≤ C l 1 |1 -e -ρ∆t-ιφ ∆t ∆t | e |ρ + φ ∆t | l , (68)
which means rm,n (e ρ∆t+ιφ , ∆t) = 1 1 -e -ρ∆t-ιφ ṙm,n (0, ∆t)O(

1 ρ + |φ ∆t | ) + O( 1 |ρ + φ ∆t | 2 ) .
We use ṙm,n (z, ∆t) to denote the first-order derivative of r m,n (z, ∆t) with respective to z.

We also have ĝe 0,0 (e ρ∆t+ιφ ) -

1 2 = 1 1 -e -ρ∆t-ιφ (1 + e -ρ∆t-ιφ ) 2 + O( 1 ρ + |φ ∆t | ) , leading to (60). Case 3. α∆t ≤ φ ≤ π -α∆t or π + α∆t ≤ φ ≤ 2π -α∆t. From the inequality (with α = ρ 3 ) ∆t 1 -e -ρ∆t-ιφ = ∆t 1 -e -ιφ 1 + e -ιφ (1 -e -ρ∆t ) 1 -e -ιφ -1 ≤ C 1 α (1 + O( ρ α )) ≤ C α , (69)
we deduce that

∆t l 2 1 l! d l ds l 1 + e s-ρ∆t-ιφ 1 -e s-ρ∆t-ιφ s=0 = ∆t l 2 l m=1 (e -ρ∆t-ιφ ) m (1 -e -ρ∆t-ιφ ) m+1 |α|=l,α i ≥1 1 α! ≤ C 1 |1 -e -ρ∆t-ιφ | l m=1 C m α m ∆t l-m ≤ 1 |1 -e -ρ∆t-ιφ | CC l α l , (70) which means rm,n (e ρ∆t+ιφ , ∆t) = 1 1 -e -ρ∆t-ιφ ṙm,n (0, ∆t)e -ρ∆t-ιφ ∆t 2(1 -e -ρ∆t-ιφ ) + O( 1 α 2 ) .
Furthermore, we also have ĝe 0,0 (e ρ∆t+ιφ ) -

1 2 = 1 1 -e -ρ∆t-ιφ (1 + e -ρ∆t-ιφ ) 2 + O( 1 α ) ,
finally leading to (60).

Concerning the inequality (61), a similar proof can be adapted, starting from the definition of h m,n given by [START_REF] Lebedev | Special Functions and Their Applications[END_REF].

Let now us introduce the following indices

: N 1 m,n = I m -I n and N 2 m,n = J m -I n . We also define I k 1 , I k 2 , ...I k 4L-5 as the ordered indices I 1 , I 2 ... I 4L-4 excluding the k-th one (1 ≤ k ≤ 4L -4). Finally, let us set N 1,k,l m,n = I k m -I l n .
From the boundary condition ( 16), we have by Laplace transform ψe

Im (s) = 4L-4 n=1 f e Im-In (s) χe In (s), 1 ≤ m ≤ 4L -4, ψe Jm (s) = 4L-4 n=1 f e Jm-In (s) χe In (s), 1 ≤ m ≤ 4L. (71)
We define A 1,e as the (4L -4) × (4L -4) matrix such that the (m, n)-th matrix element of

A 1,e is f e N 1 m,n
. In a similar way, B e is a 4L × (4L -4) matrix such that its (m, n)-th element is equal to f e N 2 m,n . As a consequence, system (71) leads to ψe J = Be ( Ã1,e ) -1 ψe

I = Fe ψe I . (72)
(The inversion of the matrix Ã1,e can be obtained from the expression (9).) Let us denote by σ 4L-5 the set of all permutations of the indices from 1 to 4L -5, an element of the set being written as i := (i 1 , i 2 , ..., i 4L-5 ) ∈ σ 4L-5 . Following Cramer's rule, for 1 ≤ l ≤ 4L, one can write that ψe

J l = 1≤m,n≤4L-4 ( Be ) l,m (( Ã1,e ) -1 ) m,n ψe In = 1≤m,n≤4L-4 ψe In 1 det( Ã1,e ) f e N 2 l,m (-1) m+n i∈σ 4L-5 δ i 4L-5 p=1 f e N 1,n,m p,ip , (73)
where δ i is the signature of the permutation for a given index i ∈ σ 4L-5 . Let us now define σ 0 4L-4 as the set σ 4L-4 excluding (1, 2, ..., 4L -4). Multiplying by det( Ã1,e ) on both sides of (73), then inverting by Laplace transform, one has

ψ e J l + j∈σ 0 4L-4 δ j f e h 1 1,j 1 * (... * (f e N 1 4L-4,j 4L-4 * ψ e J l )...) = 1≤m,n≤4L-4 f e N 2 l,m (-1) m+n * ( i∈σ 4L-5 δ i f e N 1,n,m 1,i 1 * (...(f e N 1,n,m 4L-5,i 4L-5 * ψ e In )...)). (74)
From (74) and by using (58), one gets

(ψ e J l ) k + j∈σ 0 4L-4 δ j [f e N 1 1,j 1 (... (f e N 1 4L-4,j 4L-4 ψ e J l )...)] k = 1≤m,n≤4L-4 [f e N 2 l,m (-1) m+n ( i∈σ 4L-5 δ i f e N 1,n,m 1,i 1 (...(f e N 1,n,m 4L-5,i 4L-5 ψ e In )...))] k + H l (k∆t, ∆t)∆t 2 , (75)
where H l (t, ∆t) is an analytic function of t and ∆t which can be expanded thanks to ∆t, the coefficients being some analytic functions of t. From the relation

(u v)(z) = û(z)v(z)∆t -û(z) v 0 2 ∆t -v(z) u 0 2 ∆t,
and applying a Z-transform on both sides of (75) leads to

ψe J l (z) + ψe J l (z) j∈σ 0 4L-4 δ j 4L-4 p=1 f e h 1 p,jp (z)∆t -f e h 1 p,jp (0)∆t/2 = 1≤m,n≤4L-4 ψe In (z) f e h 2 l,m (z)∆t -f e h 2 l,m (0)
∆t/2 (-1) m+n i∈σ 4L-5

δ i 4L-5 p=1 f e h 1,n,m p,ip (z)∆t -f e h 1,n,m p,ip (0)
∆t/2 + Ĥl (z, ∆t)∆t 2 .

(

Now, thanks to (59), one has

ψe J l (z) + ψe J l (z) j∈σ 0 4L-4 δ j 4L-4 p=1 fh 1 p,jp (z)∆t -f e h 1 p,jp (0)∆t/2
= 1≤m,n≤4L-4 ψe In (z) fh 2 l,m (z)∆t -f e h 2 l,m (0)∆t/2 (-1) m+n i∈σ 4L-5 δ i 4L-5 p=1 fh 1,n,m p,ip (z)∆t -fh 1,n,m p,ip (0)∆t/2 +( 1≤m,n≤4L-4 G 1,l,m,n (z, ∆t) ψe In (z))∆t 2 + G 2,l (z, ∆t) ψe J l (z)∆t 2 + Ĥl (z, ∆t)∆t 2 , (77)
where G 1,l,m,n (z, ∆t) and G 2,l (z, ∆t) are some analytic functions equal to linear combinations of products of the functions h p,q rp,q (z, ∆t) ĝ0,0 (z) -g 0,0 (0)/2 .

Next, by using Lemma 2, (58) and (99) (Lemma 4, Appendix A), we can prove the following bounds

|G 1,l,m,n (z, ∆t)| ≤ C 1,l,m,n,ρ , |G 2,l (z, ∆t)| ≤ C 2,l,ρ , (78)
for some well-chosen positive real-valued constants C 1,l,m,n,ρ and C 2,l,ρ which depend on the bounds of ψ e m,n (t) and of its derivatives. Now, by using ( 31), ( 77) is

ψe J l (z) = 1≤m,n≤4L-4 (B) l,m ((A 1 ) -1 ) m,n ψe In (z) + R l (z, ∆t)∆t 2 = 1≤n≤4L-4 (F) l,n ψe In (z) + R l (z, ∆t)∆t 2 , (79)
where F(z) is defined by

ψJ (z) = F(z) ψI (z) := B(z)(A 1 (z)) -1 ψI (z) (80)
and

R l (z, ∆t) = 1 det(A 1 )   1≤m,n≤4L-4 G 1,l,m,n (z, ∆t) ψe In (z) + G 2,l (z, ∆t) ψe J l (z) + Ĥl (z, ∆t)   . (81)
Let us also remark that, if ρ is a given large enough positive constant, then there exists ∆t 0 , such that, for ∆t ≤ ∆t 0 , F being given by (80), for any vector V of size 4L -4 and angle φ, we have

Im VT FT (e ρ∆t+ιφ )I 1 V ≤ 0. (82)
Let us now prove the following error estimate.

Theorem 2. Let us consider the scheme ( 17)- [START_REF] Hagstrom | Radiation boundary conditions for the numerical simulation of waves[END_REF]. Let ρ be a large enough given positive constant, following Theorem 1. Then there exists ∆t 0 , such that, for any ∆t ≤ ∆t 0 and t K = K∆t ≤ T , we have the following estimate

ψ K -ψ e (t K ) 2 ≤ C(L, T, ψ e (0))∆t 2 , ( 83
)
where the positive constant C depends on the length L, the maximal time T and the initial data ψ e (0).

Proof. Let us fix ρ 1 > ρ, where ρ is given in Theorem 1. We define the local error: e k m,n = (ψ e m,n ) kψ k m,n . Then one can see that

- 2ι ∆t (e k m,n -e k-1 m,n ) = ∆(e k-1 m,n + e k m,n ) + r m,n,k-1 2 , ( 84
)
where r m,n,k-1 2

is the local residual. From (79), the Z-transform of the boundary conditions leads to

êJm (z) = 4L-4 n=1 F m,n (z) ψe In (z) + R m (z, ∆t)∆t 2 - 4L-4 n=1 F m,n (z) ψIn (z) = 4L-4 n=1 F m,n (z)ê In (z) + R m (z, ∆t)∆t 2 , (85)
where 1 ≤ m ≤ 4L -4. If we define η k m,n = β -k e k m,n , with β = e ρ 1 ∆t , and after some long calculations, one obtains the following expression for the discrete L 2 -norm of the error

η K 2 2 ≤ 1 -β 2 β 2 K-1 k=0 η k 2 2 + K k=1 [ ∆t 2β L n=1 L m=1 (η k-1 m,n β -k |r m,n,k-1 2 | + η k m,n β -k+1 |r m,n,k-1 2 |)] +Im ∆t 2β 2 K k=1 4L m=1 [(η k-1 Jm + βη k Jm )(η k-1 I 1 m + βη k I 1 m -η k-1 Jm -βη k Jm )]. (86)
Now, based on (85), one gets

Im K k=1 4L m=1 [(η k-1 Jm + βη k Jm )(η k-1 I 1 m + βη k I 1 m -η k-1 Jm -βη k Jm )] = Im 1 2π 2π 0 1 + βe ιφ 2 4L m=1 ēJm (βe ιφ )(ê I 1 m (βe ιφ ) -êJm (βe ιφ )) dφ = 1 2π Im 2π 0 1 + βe ιφ 2 ēT I (βe ιφ ) FT (βe ιφ )(I 1 -F(βe ιφ ))ê I (βe ιφ )dφ + 1 2π Im 2π 0 1 + βe ιφ 2 4L m=1 Rm (βe ιφ , ∆t)ê I 1 m (βe ιφ ) dφ∆t 2 . (87)
The first integral of the above last equality does not contribute to the computation since it is less than zero. This then leads to ∆tIm

K k=1 4L m=1 [(η k-1 Jm + βη k Jm )(η k-1 I 1 m + βη k I 1 m -η k-1 Jm -βη k Jm )] ≤ ∆t 3 1 2π 2π 0 1 + βe ιφ 2 4L m=1 Rm (βe ιφ , ∆t)ê I 1 m (βe ιφ )dφ . (88)
From relations (78), ( 81) and (99) (Lemma 4, Appendix A), we conclude that ∆t 3 1 2π

2π 0 1 + βe ιφ 2 Rl (βe ιφ , ∆t)ê I 1 l (βe ιφ )dφ ≤ β 2 -1 4β 2 K k=0 η k I 1 l 2 + C∆t 5 β 2 (β 2 -1) ( max 1≤m,n≤4L-4 (det F1 (βe ιφ )) -1 G 1,l,m,n (βe ιφ , ∆t) L ∞ [0,2π] + (det F1 (βe ιφ )) -1 G 2,l (βe ιφ , ∆t) L ∞ [0,2π] + ( F1 (βe ιφ )) -1 L ∞ [0,2π] ) 2 ×( max 1≤m≤4L ∆t K k=0 ψ e,k Jm 2 + max 1≤n≤4L ∆t K k=0 ψ e,k In 2 + ∆t K k=0 H k l 2 ) ≤ β 2 -1 4β 2 K k=0 η k I 1 l 2 + ∆t 4 C(L, T, ψ e (0)), (89)
which provides ∆tIm

K k=1 4L m=1 [(η k-1 Jm + βη k Jm )(η k-1 I 1 m + βη k I 1 m -η k-1 Jm -βη k Jm )] ≤ β 2 -1 2β 2 K k=0 η k 2 2 +∆t 4 C(L, T, ψ e (0)). (90)
To get (90) from (89), we use the following estimate for three sequences {a n }, {b n } and {c n }, where their Z-transforms are â(z), b(z) and ĉ(z), respectively,

1 2π 2π 0 |â(e iϕ ) b(e iϕ )ĉ(e iϕ )|dϕ ≤ 1 2π â(e iϕ ) L ∞ [0,2π] b(e iϕ ) L 2 [0,2π] ĉ(e iϕ ) L 2 [0,2π] = 1 2π â(e iϕ ) L ∞ [0,2π] b l 2 c l 2 ≤ ε c 2 l 2 +C(ε)( â(e iϕ ) L ∞ [0,2π] b l 2 ) 2 , ( 91
) setting a 2 l 2 = K k=0 |a k | 2 .
The terms with indices I 1 l belong to the computational domain, 1 ≤ m, n ≤ L, and so they can be controlled by η k 2 . The indices η J l are not in the computational domain. However, the terms which contain an index η J l are negative in (87), so they do not need to be estimated in L 2 -norm [START_REF] Pang | Eliminating corner effects in square lattice simulation[END_REF]. By using (99), we obtain: |r m,n,k+ 1 2 | ≤ C(ψ e (0))∆t 2 . Then, by combining with (86), one gets

η K 2 2 ≤ 1 -β 2 β 2 K-1 k=0 η k 2 2 + β 2 -1 2β 2 K k=0 η k 2 2 +C(L, ψ e (0))∆t 4 + β 2 -1 2β 2 K k=0 η k 2 2 + C∆t 2 (β 2 -1) K k=1 1≤m,n≤L r m,n,k-1 2 2 . (92) Finally, we have η K 2 2 ≤ C(L, T, ψ e (0))∆t 4 , (93)
yielding ψ K -ψ e (K∆t) 2 ≤ C(L, T, ψ e (0))∆t 2 , (94)
ending hence the proof.

Discussion and Conclusion

In this paper, we proved the L 2 -stability of the Crank-Nicolson discretization of the semi-discrete Schrödinger equation with the artificial boundary condition derived in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF]. In addition, the full scheme is proved to be second-order in time. These results confirm the numerical simulations presented in [START_REF] Ji | Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains[END_REF]. Nevertheless, to fully analyze the scheme, some points remain to prospect. For example, the presented convergence analysis of the scheme does not include yet the spatial convergence analysis which will be studied in a forthcoming paper. Furthermore, fast evaluation schemes for integral operators should be used in conjunction with the ABC method to have an efficient scheme. This needs to be further prospected from the point of view of the design of the algorithm but also concerning the effect of these additional approximations on the numerical analysis of the full scheme. Finally, the extension to the nonlinear case can be obtained thanks to time-splitting schemes [START_REF] Antoine | Perfectly matched layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF][START_REF] Li | Stability and error analysis for a second-order fast approximation of the one-dimensional Schrödinger equation under absorbing boundary conditions[END_REF] and the methodology should also yield the construction of ABCs for 2D nonlocal Schrödinger equations in rectangular domains, combined with the recent results in [START_REF] Ji | Artificial boundary conditions for the semi-discretized onedimensional nonlocal Schrödinger equation[END_REF].

A Proof of Lemma 4

Since the initial condition is compactly supported, then it is square integrable, i.e. ψ e (0) * ,2 < ∞, where the L 2 -norm on the infinite lattice is defined by

ψ e * ,2 = ( (m,n)∈Z 2 |ψ e m,n | 2 ) 1/2 .
Lemma 3. The behavior of the exact solution ψ e m,n of (2) at infinity is such that

|ψ e m,n (t)| ∼ O 1 |m| + |n| for |m| + |n| → ∞. Proof. Setting λe I k (s) = χe I k (s) ge 0,0 (s)
, for 1 ≤ k ≤ 4L -4, and from ( 16), the boundary conditions can be written as

ψ e I k (t) = 4L-4 =1 g e I k -I (t) * λ e I (t), 1 ≤ k ≤ 4L -4. (95)
Then, for any term ψ e m,n (t) outside the computational domain, one has from ( 13)

ψ e m,n (t) = 4L-4 =1 g e (m,n)-I * λ e I (t). (96)
By integration by parts for (36) and for fixed t, we have

|g e m,n (t)| ∼ O 1 max(|m|, |n|) ∼ O 1 |m| + |n| (97)
and Let us now prove the following Lemma.

|ψ e m,n (t)| = | 4L-4 =1 g e (m,n)-I * λ e I (t)| ∼ O 1 |m| + |n| for |m| + |n| → ∞,
Lemma 4. We have the following L 2 -norm (mass) conservation property ∀t > 0, ψ e (t) * ,2 = ψ e (0) * ,2 .

Furthermore, the following inequality holds ψ e,(k) (t) * ,2 = ψ e,(k) (0) * ,2 ≤ 10 k ψ e (0) * ,2 , (99

)
and ψ e : R + 0 → 2 (Z 2 ) is an analytic function.

Proof. By multiplying (2) by ψ e m,n (t), taking the complex conjugate of (2) and multiplying it by -ψ e m,n (t), then summing up the two equations for (m, n) such that -N ≤ m, n ≤ N , and summing by parts the resulting equation, one gets We conclude that ψ e,(k+1) (0) * ,2 ≤ 10 ψ e,(k) (0) * ,2

and ψ e,(k) (t) * ,2 = ψ e,(k) (0) * ,2 ≤ 10 k ψ e (0) * ,2 .

(102)

The following inequality also holds: |ψ Finally, let us remark that the results above extend to the case of the Schrödinger equation with a compactly supported potential term given by V , i.e.

ι ψ = - ∂ 2 ψ ∂x 2 + ∂ 2 ψ ∂y 2 + V (x, y, t)ψ, (103)
where each function V (m∆x, n∆x, t) is analytic with respect to t.

B Matrices expressions

The matrix I 1 can be written as

I 1 =               1 I L-2 1 2 I L-2 1 2 I L-2 1 2 I L-2 1               , (104)
where I L-2 is the identity matrix of order (L -2), 1 2 = [1 1] T , and all other entries are zero. In addition, some computations show that

A 0 = - ι cot( φ 2 ) 2 I 1 , B 1 = - ι cot 2 ( φ 2 ) 4 I 1 (105)
and

A 1 =                      a 1 b 1 b 1 C 1 b 1 a 1 b 1 b 1 a 1 b 1 C 1 b 1 a 1 b 1 b 1 a 1 b 1 C 1 b 1 a 1 b 1 b 1 a 1 b 1 C 1 a 1 b 1 b 1                      , (106)
with

a 1 = - ρ + 4ι 4 sin 2 ( φ 2 ) , b 1 = - ι cot 2 ( φ 2 ) 4 ,
and the (L -2) × L matrix C 1 is given by

C 1 =           b 1 a 1 b 1 b 1 a 1 b 1 b 1 a 1 b 1 . . . b 1 a 1 b 1           . ( 107
)
We also can prove that the matrix B 2 is such that .

B 2 =                      a 2 2b 2 b 2 C 2 2b 2 a 2 b 2 b 2 a 2 2b 2 C 2 2b 2 a 2 b 2 b 2 a 2 2b 2 C 2 2b 2 a 2 b 2 b 2 a 2 2b 2 C 2 a 2 b 2 2b 2                      , (108)
The (L -2) × L matrix C 2 has the following form

C 2 =           2b 2 a 2 2b 2 2b 2 a 2 2b 2 2b 2 a 2 2b 2 . . . 2b 2 a 2 2b 2           . ( 109
)
It is finally easy to get I T 1 I 1 = I 2 . The matrix I 2 can be written as

I 2 =             2 I L-2 2 I L-2 2 I L-2 2 I L-2             . (110)

  with (a) k = a(a + 1)...(a + k -1), for k ≥ 1 and (a) 0 = 1. The function 2 F 3 is the inverse Laplace transform of the function W ((m, n),

  Figure1: Discrete computational lattice (for L = 4): the red bullets represent the outer layer points, the black ones correspond to the interior layer points, the yellow ones are the interior points and the green ones stands for the corner points.

  ) cos(ny)e -i(4-2 cos(x)-2 cos(y))t dxdy (36) and, by using the definition of Z-transform, we obtain ĝe m,n (βe ιφ ) ) cos(ny) ∞ j=0 e -K 0 (x,y)j t /(e (ρ t+ιφ) ) j dxdy

  0 (x,y)j t /(e (ρ+ιφ ∆t )∆t ) j ∆tdxdy -) cos(ny)ρ + ιφ ∆t + K 0 (x, y) dxdy + O(∆t 2 ) = 1 2ι W ((m, n), ρ + ιφ ∆t 2ι + 2) + O(∆t 2 )(45) uniformly for -α∆t ≤ φ ≤ α∆t, with an error O(∆t 2 ).

  (t)| → 0 for |m| + |n| → ∞. (98) If we define ψ e,(k) m,n (t) as the k-th order time derivative of ψ e m,n (t), then, we can also prove that |ψ e,(k) m,n (t)| ∼ O 1 |m| + |n| when |m| + |n| → ∞.

  m,n (t)| 2 -|ψ e m,n (0)| 2 ) n (0) + ψ e -N,n (t)) + x (ψ e -N,n (0) + ψ e -N,n (t)) -(ψ e N,n (0) + ψ e N,n (t)) - x (ψ e N,n (0) + ψ e N,n(t)) N (0) + ψ e m,-N (t)) + y (ψ e m,-N (0) + ψ e m,-N (t)) -(ψ e m,N (0) + ψ e m,N (t)) - y (ψ e m,N (0) + ψ e m,N (t)) .(100)When N → ∞, from Lemma 3, the left side of (100) behaves like 8N × O(1/N 2 ) = O(8/N ), and thus we obtain ∀t > 0, ψ e (t) * ,2 = ψ e (0) * ,2 .In addition, taking the k-th derivative of the semi-discrete equation (2), we have ιψ e,(k+1)

  (0)| ≤ 10 k ψ e (0) * ,2 , meaning that each function ψ e m,n (t) is analytic.

with a 2

 2 = -(2ρ + 8ι) cot(

m,n | 2 ≤ 5|ψ e,(k) m,n-1 | 2 + 5|ψ e,(k) m,n+1 | 2 + 5|ψ e,(k) m+1,n | 2 + 5|ψ e,(k) m-1,n | 2 + 80|ψ e,(k) m,n | 2 .Summing up over all the indices m and n, we have m,n |ψ e,(k+1) m,n | 2 ≤ 100 m,n |ψ e,(k) m,n | 2 .

Acknowledgment

G. PANG and S. TANG are supported partially by NSFC under contract numbers 11521202, 11832001, 11890681 and 11988102. X. ANTOINE thanks the support of the CNRS LIASFMA (Laboratoire International Associé Sino-Français en Mathématiques Appliquées), University of Lorraine and of the French National Research Agency project NABUCO, grant ANR-17-CE40-0025. Parts of the work were done while X. ANTOINE was a visiting Professor at Peking University (Department of Mechanics and Engineering Science).

Proof. From ϕ K 2 2 ≤ ϕ 0 2 2 , we have: e -ρK∆t ψ K 2 ≤ ψ 0 2 which gives (55).

Remark 2. While the Crank-Nicolson scheme is usually unconditionally stable in • 2 -norm, by the current approach we only proved the discretization (17)-( 18) to be conditionally stable in Theorem 1.

Further explorations are expected on unconditional stability.

Error analysis of the scheme

For two analytic functions u(t) and v(t), and by using iterated integration by parts, we have the following equalities

for k ≥ 0, where

has a series expansion thanks to ∆t, with coefficients which are analytic functions with respect to t.

For any analytic functions u 1 (t),..., u n (t), by induction, one has