Weizhu Bao

Émi Carles
email: remi.carles@math.cnrs.fr

AND Su § Chunmei

Qinglin Tang
email: qinglintang@scu.edu.cn

ERROR ESTIMATES OF A REGULARIZED FINITE DIFFERENCE METHOD FOR THE LOGARITHMIC SCHR ÖDINGER EQUATION *

Keywords: Logarithmic Schrödinger equation, logarithmic nonlinearity, regularized logarithmic Schrödinger equation, semi-implicit finite difference method, error estimates, convergence rate AMS subject classifications. 35Q40, 35Q55, 65M15, 81Q05

We present a regularized finite difference method for the logarithmic Schrödinger equation (LogSE) and establish its error bound. Due to the blow-up of the logarithmic nonlinearity, i.e. ln ρ → -∞ when ρ → 0 + with ρ = |u| 2 being the density and u being the complex-valued wave function or order parameter, there are significant difficulties in designing numerical methods and establishing their error bounds for the LogSE. In order to suppress the round-off error and to avoid blow-up, a regularized logarithmic Schrödinger equation (RLogSE) is proposed with a small regularization parameter 0 < ε ≪ 1 and linear convergence is established between the solutions of RLogSE and LogSE in term of ε. Then a semi-implicit finite difference method is presented for discretizing the RLogSE and error estimates are established in terms of the mesh size h and time step τ as well as the small regularization parameter ε. Finally numerical results are reported to confirm our error bounds.

1. Introduction. We consider the logarithmic Schrödinger equation (LogSE) which arises in a model of nonlinear wave mechanics (cf. [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF]), (1.1) i∂ t u(x, t) + ∆u(x, t) = λu(x, t) ln(|u(x, t)| 2 ), x ∈ Ω, t > 0, u(x, 0) = u 0 (x), x ∈ Ω, where t is time, x ∈ R d (d = 1, 2, 3) is the spatial coordinate, λ ∈ R\{0} measures the force of the nonlinear interaction, u := u(x, t) ∈ C is the dimensionless wave function or order parameter and Ω = R d or Ω ⊂ R d is a bounded domain with homogeneous Dirichlet or periodic boundary condition 1 fixed on the boundary. It admits applications to quantum mechanics [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF][START_REF] Gaussons | Solitons of the logarithmic Schrödinger equation[END_REF], quantum optics [START_REF] Buljan | Incoherent white light solitons in logarithmically saturable noninstantaneous nonlinear media[END_REF][START_REF] Krolikowski | Unified model for partially coherent solitons in logaritmically nonlinear media[END_REF], nuclear physics [START_REF] Hefter | Application of the nonlinear Schrödinger equation with a logarithmic inhomogeneous term to nuclear physics[END_REF], transport and diffusion phenomena [START_REF] Hansson | Propagation of partially coherent solitons in saturable logarithmic media: A comparative analysis[END_REF][START_REF] Martino | Logarithmic Schrödinger-like equation as a model for magma transport[END_REF], open quantum systems [START_REF] Hernandez | General properties of Gausson-conserving descriptions of quantal damped motion[END_REF][START_REF] Yasue | Quantum mechanics of nonconservative systems[END_REF], effective quantum gravity [START_REF] Zloshchastiev | Logarithmic nonlinearity in theories of quantum gravity: Origin of time and observational consequences[END_REF], theory of superfluidity and Bose-Einstein condensation [START_REF] Avdeenkov | Quantum bose liquids with logarithmic nonlinearity: Self-sustainability and emergence of spatial extent[END_REF].

The logarithmic Schrödinger equation enjoys three conservation laws, mass, momentum and energy [START_REF]Semilinear Schrödinger equations[END_REF][START_REF] Cazenave | Équations d'évolution avec non linéarité logarithmique[END_REF], like in the case of the nonlinear Schrödinger equation with a power-like nonlinearity (e.g. cubic):

M (t) : = u(•, t) 2 L 2 (Ω) = Ω |u(x, t)| 2 dx ≡ Ω |u 0 (x)| 2 dx = M (0), P (t)
: = Im Ω u(x, t)∇u(x, t)dx ≡ Im Ω u 0 (x)∇u 0 (x)dx = P (0), t ≥ 0,

E(t) : = Ω |∇u(x, t)| 2 dx + λF (|u(x, t)| 2 ) dx ≡ Ω |∇u 0 (x)| 2 + λF (|u 0 (x)| 2 ) dx = E(0), (1.2)
where Im f and f denote the imaginary part and complex conjugate of f , respectively, and

(1.3) F (ρ) = ρ 0 ln(s)ds = ρ ln ρ -ρ, ρ ≥ 0.
On a mathematical level, the logarithmic nonlinearity possesses several features that make it quite different from more standard nonlinear Schrödinger equations. First, the nonlinearity is not locally Lipschitz continuous because of the behavior of the logarithm function at the origin. Note that in view of numerical simulation, this singularity of the "nonlinear potential" λ ln(|u(x, t)| 2 ) makes the choice of a discretization quite delicate. The second aspect is that whichever the sign of λ, the nonlinear potential energy in E has no definite sign. In fact, whether the nonlinearity is repulsive/attractive (or defocusing/focusing) depends on both λ and the value of the density ρ := ρ(x, t) = |u(x, t)| 2 . When λ > 0, then the nonlinearity λρ ln ρ is repulsive when ρ > 1; and respectively, it is attractive when 0 < ρ < 1. On the other hand, when λ < 0, then the nonlinearity λρ ln ρ is attractive when ρ > 1; and respectively, it is repulsive when 0 < ρ < 1. Therefore, solving the Cauchy problem for (1.1) is not a trivial issue, and constructing solutions which are defined for all time requires some work; see [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF][START_REF] Cazenave | Équations d'évolution avec non linéarité logarithmique[END_REF][START_REF] Guerrero | Global H 1 solvability of the 3D logarithmic Schrödinger equation[END_REF]. Essentially, the outcome is that if u 0 belongs to (a subset of) H 1 (Ω), (1.1) has a unique, global solution, regardless of the space dimension d (see also Theorem 2.2 below).

Next, the large time behavior reveals new phenomena. A first remark suggests that nonlinear effects are weak. Indeed, unlike what happens in the case of a homogeneous nonlinearity (classically of the form λ|u| p u), replacing u with ku (k ∈ C \ {0}) in (1.1) has only little effect, since we have i∂ t (ku) + ∆(ku) = λku ln |ku| 2 -λ(ln |k| 2 )ku .

The scaling factor thus corresponds to a purely time-dependent gauge transform: ku(x, t)e -itλ ln |k| 2 solves (1.1) (with initial datum ku 0 ). In particular, the size of the initial datum does not influence the dynamics of the solution. In spite of this property which is reminiscent of linear equations, nonlinear effects are stronger in (1.1) than in, say, cubic Schrödinger equations in several respects. For Ω = R d , it was established in [START_REF] Cazenave | Stable solutions of the logarithmic Schrödinger equation[END_REF] that in the case λ < 0, no solution is dispersive (not even for small data, in view of the above remark), while if λ > 0, the results from [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] show that every solution disperses, at a faster rate than for the linear equation.

In view of the gauge invariance of the nonlinearity, for Ω = R d , (1.1) enjoys the standard Galilean invariance: if u(x, t) solves (1.1), then, for any v ∈ R d , so does u(x -2vt, t)e iv•x-i|v| 2 t .

A remarkable feature of (1.1) is that it possesses a large set of explicit solutions. In the case Ω = R d : if u 0 is Gaussian, u(•, t) is Gaussian for all time, and solving (1.1) amounts to solving ordinary differential equations [START_REF] Bia Lynicki-Birula | Nonlinear wave mechanics[END_REF]. For simplicity of notation, we take the one-dimensional case as an example. If the initial data in (1.1) with Ω = R is taken as

u 0 (x) = b 0 e -a 0 2 x 2 +ivx , x ∈ R,
where a 0 , b 0 ∈ C and v ∈ R are given constants satisfying α 0 := Re a 0 > 0 with Re f denoting the real part of f , then the solution of (1.1) is given by [START_REF] Ardila | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF][START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] (

1.4) u(x, t) = b 0 r(t) e i(vx-v 2 t)+Y (x-2vt,t) , x ∈ R, t ≥ 0, with (1.5) Y (x, t) = -iφ(t) -α 0 x 2 2r(t) 2 + i ṙ(t) r(t) x 2 4 , x ∈ R, t ≥ 0,
where φ := φ(t) ∈ R and r := r(t) > 0 solve the ODEs [START_REF] Ardila | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF][START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF]

φ = α 0 r 2 + λ ln |b 0 | 2 -λ ln r, φ(0) = 0, r = 4α 2 0 r 3 + 4λα 0 r , r(0) = 1, ṙ(0) = -2 Im a 0 .
(1.6)

In the case λ < 0, the function r is (time) periodic (in agreement with the absence of dispersive effects). In particular, if a 0 = -λ > 0, it follows from (1.6) that r(t) ≡ 1 and φ(t) = φ 0 t with φ 0 = λ ln(|b 0 | 2 ) -1 , which generates the uniformly moving Gausson as [START_REF] Ardila | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF][START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] (

1.7) u(x, t) = b 0 e λ 2 (x-2vt) 2 +i(vx-(φ0+v 2 )t) , x ∈ R, t ≥ 0.
As a very special case with b 0 = e 1/2 and v = 0 such that φ 0 = 0, one can get the static Gausson as

(1.8) u(x, t) = e 1/2 e λ|x| 2 /2 , x ∈ R, t ≥ 0.
This special solution is orbitally stable [START_REF] Cazenave | Stable solutions of the logarithmic Schrödinger equation[END_REF][START_REF] Cazenave | Orbital stability of standing waves for some nonlinear Schrödinger equations[END_REF]. On the other hand, in the case λ > 0, it is proven in [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] that for general initial data (not necessarily Gaussian), there exists a universal dynamics. For extensions to higher dimensions, we refer to [START_REF] Ardila | Orbital stability of Gausson solutions to logarithmic Schrödinger equations[END_REF][START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] and references therein. Therefore, (1.1) possesses several specific features, which make it quite different from the nonlinear Schrödinger equation. Different numerical methods have been proposed and analyzed for the nonlinear Schrödinger equation with smooth nonlinearity (e.g. cubic nonlinearity) in the literature, such as the finite difference methods [START_REF] Bao | Uniform error estimates of finite difference methods for the nonlinear Schrödinger equation with wave operator[END_REF][START_REF]Optimal error estimates of finite difference methods for the Gross-Pitaevskii equation with angular momentum rotation[END_REF], finite element methods [START_REF] Akrivis | On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schrödinger equation[END_REF][START_REF] Karakashian | A space-time finite element method for the nonlinear Schrödinger equation: the continuous Galerkin method[END_REF] and the time-splitting pseudospectral methods [START_REF] Bao | Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation[END_REF][START_REF] Taha | Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrödinger equation[END_REF]. However, they cannot be applied to the LogSE (1.1) directly due to the blow-up of the logarithmic nonlinearity, i.e. ln ρ → -∞ when ρ → 0 + . The main aim of this paper is to present a regularized finite difference method for the LogSE (1.1) by introducing a proper regularized logarithmic Schrödinger equation (RLogSE) and then discretizing the RLogSE via a semi-implicit finite difference method. Error estimates will be established between the solutions of LogSE and RLogSE as well as their numerical approximations.

The rest of the paper is organized as follows. In Section 2, we propose a regularized version of (1.1) with a small regularization parameter 0 < ε ≪ 1, and analyze its properties, as well as the convergence of its solution to the solution of (1.1). In Section 3, we introduce a semi-implicit finite difference method for discretizing the regularized logarithmic Schrödinger equation, and prove an error estimate, in which the dependence of the constants with respect to the regularization parameter ε is tracked very explicitly. Finally, numerical results are provided in Section 4 to confirm our error bounds and to demonstrate the efficiency and accuracy of the proposed numerical method.

Throughout the paper, we use H m (Ω) and • H m (Ω) to denote the standard Sobolev spaces and their norms, respectively. In particular, the norm and inner product of L 2 (Ω) = H 0 (Ω) are denoted by • L 2 (Ω) and (•, •), respectively. Moreover, we adopt A B to mean that there exists a generic constant C > 0 independent of the regularization parameter ε, the time step τ and the mesh size h such that A ≤ C B, and c means the constant C depends on c.

2.

A regularized logarithmic Schrödinger equation. It turns out that a direct simulation of the solution of (1.1) is very delicate, due to the singularity of the logarithm at the origin, as discussed in the introduction. Instead of working directly with (1.1), we shall consider the following regularized logarithmic Schrödinger equation (RLogSE) with a samll regularized parameter 0 < ε ≪ 1 as

(2.1) i∂ t u ε (x, t) + ∆u ε (x, t) = λu ε (x, t) ln (ε + |u ε (x, t)|) 2 , x ∈ Ω, t > 0, u ε (x, 0) = u 0 (x), x ∈ Ω.
2.1. Conserved quantities. For the RLogSE (2.1), it can be similarly deduced that the mass, momentum, and energy are conserved.

Proposition 2.1. The mass, momentum, and 'regularized' energy are formally conserved for the RLogSE (2.1):

M ε (t) := Ω |u ε (x, t)| 2 dx ≡ Ω |u 0 (x)| 2 dx = M (0), P ε (t) := Im Ω u ε (x, t)∇u ε (x, t)dx ≡ Im Ω u 0 (x)∇u 0 (x)dx = P (0), t ≥ 0, E ε (t) := Ω |∇u ε (x, t)| 2 + λF ε (|u ε (x, t)| 2 ) (x, t)dx ≡ Ω |∇u 0 (x)| 2 + λF ε (|u 0 (x)| 2 ) dx = E ε (0), (2.2)
where

(2.3) F ε (ρ) = ρ 0 ln(ε + √ s) 2 ds = ρ ln (ε + √ ρ) 2 -ρ + 2ε √ ρ -ε 2 ln (1 + √ ρ/ε) 2 , ρ ≥ 0.
Proof. The conservation for mass and momentum is standard, and relies on the fact that the right hand side of (2.1) involves u ε multiplied by a real number. For the energy E ε (t), we compute

d dt E ε (t) = 2 Re Ω ∇u ε • ∇∂ t u ε + λu ε ∂ t u ε ln(ε + |u ε |) 2 -λu ε ∂ t u ε (x, t)dx + 2λ Ω ∂ t |u ε | ε + |u ε | 2 -ε 2 ε + |u ε | (x, t)dx = 2 Re Ω ∂ t u ε -∆u ε + λu ε ln(ε + |u ε |) 2 (x, t)dx = 2 Re Ω i|∂ t u ε | 2 (x, t)dx = 0, t ≥ 0,
which completes the proof. Note however that since the above 'regularized' energy involves L 1 -norm of u ε for any ε > 0, E ε is obviously well-defined for u 0 ∈ H 1 (Ω) when Ω has finite measure, but not when Ω = R d . This aspect is discussed more into details in Subsections 2.3.3 and 2.4.

The Cauchy problem.

For α > 0 and Ω = R d , denote by L 2 α the weighted

L 2 space L 2 α := {v ∈ L 2 (R d ), x -→ x α v(x) ∈ L 2 (R d )},
where x := 1 + |x| 2 , with norm

v L 2 α := x α v(x) L 2 (R d ) .
In the case where Ω is bounded, we simply set L 2 α = L 2 (Ω). Regarding the Cauchy problems (1.1) and (2.1), we have the following result.

0 ∈ H 1 0 (Ω) ∩ L 2 α , for some 0 < α ≤ 1. • There exists a unique, global solution u ∈ L ∞ loc (R; H 1 0 (Ω) ∩ L 2 α ) to (1.1), and a unique, global solution u ε ∈ L ∞ loc (R; H 1 0 (Ω) ∩ L 2 α ) to (2.1). • If in addition u 0 ∈ H 2 (Ω), then u, u ε ∈ L ∞ loc (R; H 2 (Ω)). • In the case Ω = R d , if in addition u 0 ∈ H 2 ∩L 2 2 , then u, u ε ∈ L ∞ loc (R; H 2 ∩L 2 2 )
. Proof. This result can be proved by using more or less directly the arguments invoked in [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF]. First, for fixed ε > 0, the nonlinearity in (2.1) is locally Lipschitz, and grows more slowly than any power for large |u ε |. Therefore, the standard Cauchy theory for nonlinear Schrödinger equations applies (see in particular [START_REF]Semilinear Schrödinger equations[END_REF]Corollary 3.3.11 and Theorem 3.4.1]), and so if

u 0 ∈ H 1 0 (Ω), then (2.1) has a unique solution u ε ∈ L ∞ loc (R; H 1 0 (Ω)).
Higher Sobolev regularity is propagated, with controls depending on ε in general.

A solution u of (1.1) can be obtained by compactness arguments, by letting ε → 0 in (2.1), provided that we have suitable bounds independent of ε > 0. We have

i∂ t ∇u ε + ∆∇u ε = 2λ ln (ε + |u ε |) ∇u ε + 2λ u ε ε + |u ε | ∇|u ε |.
The standard energy estimate (multiply the above equation by ∇u ε , integrate over Ω and take the imaginary part) yields, when Ω = R d or when periodic boundary conditions are considered, 1 2

d dt ∇u ε 2 L 2 (Ω) ≤ 2|λ| Ω |u ε | ε + |u ε | |∇|u ε || |∇u ε |dx ≤ 2|λ ∇u ε 2 L 2 (Ω) .
Gronwall lemma yields a bound for u ε in L ∞ (0, T ; H 1 (Ω)), uniformly in ε > 0, for any given T > 0. Indeed, the above estimate uses the property

Im Ω ∇u ε • ∆∇u ε dx = 0,
which needs not be true when Ω is bounded and u ε satisfies homogeneous Dirichlet boundary conditions. In that case, we use the conservation of the energy E ε (Proposition 2.1), and write

∇u ε (t) 2 L 2 (Ω) ≤ E ε (u 0 ) + 2|λ| Ω |u ε (x, t)| 2 |ln (ε + |u ε (x, t)|)| dx + 2ε|λ| u ε (t) L 1 (Ω) + 2|λ|ε 2 Ω |ln (1 + |u ε (x, t)|/ε)| dx 1 + ε|Ω| 1/2 u ε (t) L 2 (Ω) + Ω |u ε (x, t)| 2 |ln (ε + |u ε (x, t)|)| dx 1 + Ω |u ε (x, t)| 2 |ln (ε + |u ε (x, t)|)| dx, t ≥ 0,
where we have used Cauchy-Schwarz inequality and the conservation of the mass M ε (t). Writing, for 0 < η ≪ 1,

Ω |u ε | 2 |ln (ε + |u ε |)| dx ε+|u ε |>1 |u ε | 2 (ε + |u ε |) η dx + ε+|u ε |<1 |u ε | 2 (ε + |u ε |) -η dx u ε L 2 (Ω) + u ε 2+η L 2+η (Ω) + u ε 2-η L 2-η (Ω) 1 + ∇u ε dη/2 L 2 (Ω) ,
where we have used the interpolation inequality (see e.g. [START_REF] Nirenberg | On elliptic partial differential equations[END_REF])

u L p (Ω) u 1-α L 2 (Ω) ∇u α L 2 (Ω) + u L 2 (Ω) , p = 2d d -2α , 0 ≤ α < 1,
we obtain again that u ε is bounded in L ∞ (0, T ; H 1 (Ω)), uniformly in ε > 0, for any given T > 0.

In the case where Ω is bounded, compactness arguments show that u ε converges to a solution u to (1.1); see [START_REF]Semilinear Schrödinger equations[END_REF][START_REF] Cazenave | Équations d'évolution avec non linéarité logarithmique[END_REF]. When Ω = R d , compactness in space is provided by multiplying (2.1) with x 2α u ε and integrating in space:

d dt u ε 2 L 2 α = 4α Im x • ∇u ε x 2-2α u ε (t) dx x 2α-1 u ε L 2 (Ω) ∇u ε L 2 (Ω) ,
where we have used Cauchy-Schwarz inequality. Recalling that 0 < α ≤ 1,

x 2α-1 u ε L 2 (Ω) ≤ x α u ε L 2 (Ω) = u ε L 2 α ,
and we obtain a bound for u ε in L ∞ (0, T ; H 1 (Ω) ∩ L 2 α ) which is uniform in ε. Uniqueness of such a solution for (1.1) follows from the arguments of [START_REF] Cazenave | Équations d'évolution avec non linéarité logarithmique[END_REF], involving a specific algebraic inequality, generalized in Lemma 2.4 below. Note that at this stage, we know that u ε converges to u by compactness arguments, so we have no convergence estimate. Such estimates are established in Subsection 2.3.

To prove the propagation of the H 2 regularity, we note that differentiating twice the nonlinearity in (2.1) makes it unrealistic to expect direct bounds which are uniform in ε. To overcome this difficulty, the argument proposed in [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] relies on Kato's idea: instead of differentiating the equation twice in space, differentiate it once in time, and use the equation to infer H 2 regularity. This yields the second part of the theorem.

To establish the last part of the theorem, we prove that u ∈ L ∞ loc (R; L 2 2 ) and the same approach applies to u ε . It follows from (1.1) that

d dt u(t) 2 L 2 2 = -2 Im R d x 4 u(x, t)∆u(x, t)dx = 8 Im R d x 2 u(x, t) x • ∇u(x, t)dx ≤ 8 u(t) L 2 2 x • ∇u(t) L 2 (R d ) . (2.4)
By Cauchy-Schwarz inequality and integration by parts, we have

x • ∇u(t) 2 L 2 (R d ) ≤ d j=1 d k=1 R d x 2 j ∂u(x, t) ∂x k ∂u(x, t) ∂x k dx = -2 R d u(x, t) x • ∇u(x, t)dx - R d |x| 2 u(x, t)∆u(x, t)dx ≤ 1 2 x • ∇u(t) 2 L 2 (R d ) + 2 u(t) 2 L 2 (R d ) + 1 2 u(t) 2 L 2 2 + 1 2 ∆u(t) 2 L 2 (R d ) ,
which yields directly that

x • ∇u(t) L 2 (R d ) ≤ 2 u(t) L 2 (R d ) + u(t) L 2 2 + ∆u(t) L 2 (R d )
. This together with (2.4) gives that

d dt u(t) L 2 2 ≤ 4 x • ∇u(t) L 2 (R d ) ≤ 4 u(t) L 2 2 + 8 u(t) L 2 (R d ) + 4 ∆u(t) L 2 (R d ) .
Since we already know that u ∈ L ∞ loc (R; H 2 (R d )), Gronwall lemma completes the proof.

Remark 2.1. We emphasize that if u 0 ∈ H k (R d ), k ≥ 3, we cannot guarantee in general that this higher regularity is propagated in (1.1), due to the singularities stemming from the logarithm. Still, this property is fulfilled in the case where u 0 is Gaussian, since then u remains Gaussian for all time. However, our numerical tests, in the case where the initial datum is chosen as the dark soliton of the cubic Schrödinger equation multiplied by a Gaussian, suggest that even the H 3 regularity is not propagated in general.

2.3.

Convergence of the regularized model. In this subsection, we show the approximation property of the regularized model (2.1) to (1.1).

2.3.1.

A general estimate. We prove: Lemma 2.3. Suppose the equation is set on Ω, where Ω = R d , or Ω ⊂ R d is a bounded domain with homogeneous Dirichlet or periodic boundary condition, then we have the general estimate:

(2.5) d dt u ε (t) -u(t) 2 L 2 (Ω) ≤ 4|λ| u ε (t) -u(t) 2 L 2 (Ω) + ε u ε (t) -u(t) L 1 (Ω) .
Before giving the proof of Lemma 2.3, we introduce the following lemma, which is a variant of [12, Lemma 9.

|Im ((f ε (z 1 ) -f ε (z 2 )) (z 1 -z 2 ))| ≤ |z 1 -z 2 | 2 , z 1 , z 2 ∈ C. Proof. Notice that Im [(f ε (z 1 ) -f ε (z 2 )) (z 1 -z 2 )] = 1 2 [ln(ε + |z 1 |) -ln(ε + |z 2 |)] Im(z 1 z 2 -z 1 z 2 ). Supposing, for example, 0 < |z 2 | ≤ |z 1 |, we can obtain that |ln(ε + |z 1 |) -ln(ε + |z 2 |)| = ln 1 + |z 1 | -|z 2 | ε + |z 2 | ≤ |z 1 | -|z 2 | ε + |z 2 | ≤ |z 1 -z 2 | |z 2 | , and
|Im(z 1 z 2 -z 1 z 2 )| = |z 2 (z 1 -z 2 ) + z 2 (z 2 -z 1 ))| ≤ 2|z 2 | |z 1 -z 2 |.
Otherwise the result follows by exchanging z 1 and z 2 .

Proof. (Proof of Lemma 2.3) Subtracting (1.1) from (2.1), we see that the error function e ε := u ε -u satisfies

i∂ t e ε + ∆e ε = λ u ε ln(ε + |u ε |) 2 -u ln(|u| 2 ) .
Multiplying the error equation by e ε (t), integrating in space and taking the imaginary parts, we can get by using Lemma 2.4 that 1 2

d dt e ε (t) 2 L 2 (Ω) = 2λ Im Ω [u ε ln(ε + |u ε |) -u ln(|u|)] (u ε -u)(x, t)dx ≤ 2|λ| e ε (t) 2 L 2 (Ω) + 2|λ| Ω e ε u [ln(ε + |u|) -ln(|u|)] (x, t)dx ≤ 2|λ| e ε (t) 2 L 2 (Ω) + 2ε|λ| e ε (t) L 1 (Ω) ,
where we have used the general estimate 0 ≤ ln(1 + |x|) ≤ |x|.

Convergence for bounded domain.

If Ω has finite measure, then we can have the following convergence behavior.

Proposition 2.5. Assume that Ω has finite measure, and let u 0 ∈ H 2 (Ω). For any T > 0, we have

(2.6) u ε -u L ∞ (0,T ;L 2 (Ω)) ≤ C 1 ε, u ε -u L ∞ (0,T ;H 1 (Ω)) ≤ C 2 ε 1/2 ,
where C 1 depends on |λ|, T , |Ω| and C 2 depends on |λ|, T , |Ω| and u 0 H 2 (Ω) . Proof. Note that e ε (t

) L 1 (Ω) ≤ |Ω| 1/2 e ε (t) L 2 (Ω) , then it follows from (2.5) that d dt e ε (t) L 2 (Ω) ≤ 2|λ| e ε (t) L 2 (Ω) + 2ε|λ||Ω| 1/2 .
Applying Gronwall's inequality, we immediately get that

e ε (t) L 2 (Ω) ≤ e ε (0) L 2 (Ω) + ε|Ω| 1/2 e 2|λ|t = ε|Ω| 1/2 e 2|λ|t .
The convergence rate in H 1 follows from the property u ε , u ∈ L ∞ loc (R; H 2 (Ω)) and the Gagliardo-Nirenberg inequality [START_REF] Leoni | A first course in Sobolev spaces[END_REF],

∇v L 2 (Ω) v 1/2 L 2 (Ω) ∆v 1/2 L 2 (Ω) ,
which completes the proof. Remark 2.2. The weaker rate in the H 1 estimate is due to the fact that Lemma 2.3 is not easily adapted to H 1 estimates, because of the presence of the logarithm. Differentiating (1.1) and (2.1) makes it hard to obtain the analogue in Lemma 2.3. This is why we bypass this difficulty by invoking boundedness in H 2 and interpolating with the error bound at the

L 2 level. If we have u ε , u ∈ L ∞ loc (R; H k (Ω)) for k > 2, then the convergence rate in H 1 (Ω) can be improved as e ε L ∞ (0,T ;H 1 (Ω)) ε k-1 k ,
by using the inequality (see e.g. [START_REF] Nirenberg | On elliptic partial differential equations[END_REF]):

v H 1 (Ω) v 1-1/k L 2 (Ω) v 1/k H k (Ω) .
2.3.3. Convergence for the whole space. In order to prove the convergence rate of the regularized model (2.1) to (1.1) for the whole space, we need the following lemma.

Lemma 2.6.

For d = 1, 2, 3, if v ∈ L 2 (R d ) ∩ L 2 2 , then we have (2.7) v L 1 (R d ) ≤ C v 1-d/4 L 2 (R d ) v d/4 L 2 2 ,
where C > 0 depends on d.

Proof. Applying the Cauchy-Schwarz inequality, we can get for fixed r > 0,

v L 1 (R d ) = |x|≤r |v(x)|dx + |x|≥r |x| 2 |v(x)| |x| 2 dx r d/2 |x|≤r |v(x)| 2 dx 1 2 + |x|≥r |x| 4 |v(x)| 2 dx 1 2 |x|≥r 1 |x| 4 dx 1 2 r d/2 v L 2 (R d ) + r d/2-2 v L 2 2 .
Then (2.7) can be obtained by setting

r = v L 2 2 / v L 2 (R d ) 1/2 . Proposition 2.7. Assume that Ω = R d , 1 ≤ d ≤ 3, and let u 0 ∈ H 2 (R d ) ∩ L 2 2 .
For any T > 0, we have

u ε -u L ∞ (0,T ;L 2 (R d )) ≤ C 1 ε 4 4+d , u ε -u L ∞ (0,T ;H 1 (R d ))) ≤ C 2 ε 2 4+d
,

where C 1 depends on d, |λ|, T , u 0 L 2 2 and C 2 depends on additional u 0 H 2 (R d ) . Proof. Applying (2.7) and the Young's inequality, we deduce that

ε e ε (t) L 1 (R d ) ≤ εC d e ε (t) 1-d/4 L 2 (R d ) e ε (t) d/4 L 2 2 ≤ C d e ε (t) 2 L 2 (R d ) + ε 8 4+d e ε (t) 2d 4+d L 2 2
, which together with (2.5) gives that

d dt e ε (t) 2 L 2 (R d ) ≤ 4|λ|(1 + C d ) e ε (t) 2 L 2 (R d ) + 4C d |λ|ε 8 4+d e ε (t) 2d 4+d L 2 2 .
Gronwall lemma yields

e ε (t) L 2 (R d ) ≤ ε 4 4+d e ε (t) d 4+d L 2 2 e tC d,|λ| .

The proposition follows by recalling that

u ε , u ∈ L ∞ loc (R; H 2 (R d ) ∩ L 2 2 ). Remark 2.3. If we have u ε , u ∈ L ∞ loc (R; L 2 m ) for m > 2,
then by applying the inequality

ε v L 1 (R d ) ε v 1-d 2m L 2 (R d ) v d 2m L 2 m v 2 L 2 (R d ) + ε 4m 2m+d v 2d 2m+d L 2 m ,
which can be proved like above, the convergence rate can be improved as

u ε -u L ∞ (0,T ;L 2 (R d )) ε 2m 2m+d . Remark 2.4. If in addition u ε , u ∈ L ∞ loc (R; H s (R d )) for s > 2, then the conver- gence rate in H 1 (R d ) can be improved as e ε L ∞ (0,T ;H 1 (R d ))) ≤ Cε 2m 2m+d s-1 s ,
by using the Gagliardo-Nirenberg inequality:

∇v L 2 (R d ) ≤ C v 1-1/s L 2 (R d ) ∇ s v 1/s L 2 (R d ) .
The previous two remarks apply typically in the case of Gaussian initial data.

Convergence of the energy.

In this subsection we will show the convergence of the energy

E ε (u 0 ) → E(u 0 ). Proposition 2.8. For u 0 ∈ H 1 (Ω) ∩ L 1 (Ω), the energy E ε (u 0 ) converges to E(u 0 ) with |E ε (u 0 ) -E(u 0 )| ≤ 4 ε|λ| u 0 L 1 (Ω) .
Proof. It can be deduced from the definition that

|E ε (u 0 ) -E(u 0 )| = 2|λ| ε u 0 L 1 (Ω) + Ω |u 0 (x)| 2 [ln(ε + |u 0 (x)| -ln(|u 0 (x)|)] dx -ε 2 Ω ln (1 + |u 0 (x)|/ε) dx ≤ 4 ε|λ| u 0 L 1 (Ω) , which completes the proof. Remark 2.5. If Ω is bounded, then H 1 (Ω) ⊆ L 1 (Ω). If Ω = R d , then Lemma 2.6 (and its natural generalizations) shows that H 1 (R) ∩ L 2 1 ⊆ L 1 (R), and if d = 2, 3, H 1 (R d ) ∩ L 2 2 ⊆ L 1 (R d
). Remark 2.6. This regularization is reminiscent of the one considered in [START_REF] Carles | Universal dynamics for the defocusing logarithmic Schrödinger equation[END_REF] in order to prove (by compactness arguments) that (1.1) has a solution,

(2.8) i∂ t u ε (x, t) + ∆u ε (x, t) = λu ε (x, t) ln ε + |u ε (x, t)| 2 , x ∈ Ω, t > 0.
With that regularization, it is easy to adapt the error estimates established above for (2.1). Essentially, ε must be replaced by √ ε (in Lemma 2.3, and hence in its corollaries).

3. A regularized semi-implicit finite difference method. In this section, we study the approximation properties of a finite difference method for solving the regularized model (2.1). For simplicity of notation, we set λ = 1 and only present the numerical method for the RLogSE (2.1) in 1D, as extensions to higher dimensions are straightforward. When d = 1, we truncate the RLogSE on a bounded computational interval Ω = (a, b) with homogeneous Dirichlet boundary condition (here |a| and b are chosen large enough such that the truncation error is negligible):

(3.1) i∂ t u ε (x, t) + ∂ xx u ε (x, t) = u ε (x, t) ln(ε + |u ε (x, t)|) 2 , x ∈ Ω, t > 0, u ε (x, 0) = u 0 (x), x ∈ Ω; u ε (a, t) = u ε (b, t) = 0, t ≥ 0,
3.1. A finite difference scheme and main results on error bounds. Choose a mesh size h := ∆x = (b -a)/M with M being a positive integer and a time step τ := ∆t > 0 and denote the grid points and time steps as

x j := a + jh, j = 0, 1, • • • , M ; t k := kτ, k = 0, 1, 2, . . .

Define the index sets

T M = {j | j = 1, 2, • • • , M -1}, T 0 M = {j | j = 0, 1, • • • , M }.
Let u ε,k j be the approximation of u ε (x j , t k ), and denote u ε,k = (u ε,k 0 , u ε,k 1 , . . . , u ε,k M ) T ∈ C M+1 as the numerical solution vector at t = t k . Define the standard finite difference operators

δ c t u k j = u k+1 j -u k-1 j 2τ , δ + x u k j = u k j+1 -u k j h , δ 2 x u k j = u k j+1 -2u k j + u k j-1 h 2 .
Denote

X M = v = (v 0 , v 1 , . . . , v M ) T | v 0 = v M = 0 ⊆ C M+1 ,
equipped with inner products and norms defined as (recall that

u 0 = v 0 = u M = v M = 0 by Dirichlet boundary condition) (u, v) = h M-1 j=1 u j v j , u, v = h M-1 j=0 u j v j , u ∞ = sup j∈T 0 M |u j |; u 2 = (u, u), |u| 2 H 1 = δ + x u, δ + x u , u 2 H 1 = u 2 + |u| 2 H 1 . (3.2)
Then we have for u, v ∈ X M ,

(3.3) (-δ 2 x u, v) = δ + x u, δ + x v = (u, -δ 2 x v).
Consider a semi-implicit finite difference (SIFD) discretization of (3.1) as following

(3.4) iδ c t u ε,k j = - 1 2 δ 2 x (u ε,k+1 j + u ε,k-1 j ) + u ε,k j ln(ε + |u ε,k j |) 2 , j ∈ T M , k ≥ 1.
The boundary and initial conditions are discretized as

(3.5) u ε,k 0 = u ε,k M = 0, k ≥ 0; u ε,0 j = u 0 (x j ), j ∈ T 0 M .
In addition, the first step u ε,1 j can be obtained via the Taylor expansion as

(3.6) u ε,1 j = u ε,0 j + τ u 1 (x j ), j ∈ T 0 M ,
where

u 1 (x) := ∂ t u ε (x, 0) = i u ′′ 0 (x) -u 0 (x) ln(ε + |u 0 (x)|) 2 , a ≤ x ≤ b.
Let 0 < T < T max with T max the maximum existence time of the solution u ε to the problem (3.1) for a fixed 0 ≤ ε ≪ 1. By using the standard von Neumann analysis, we can show that the discretization (3.4) is conditionally stable under the stability condition

(3.7) 0 < τ ≤ 1 2 max{| ln ε|, ln(ε + max j∈TM |u ε,k j |)} , 0 ≤ k ≤ T τ .
Define the error functions e ε,k ∈ X M as

(3.8) e ε,k j = u ε (x j , t k ) -u ε,k j , j ∈ T 0 M , 0 ≤ k ≤ T τ ,
where u ε is the solution of (3.1). Then we have the following error estimates for (3.4) with (3.5) and (3.6). Theorem 3.1 (Main result). Assume that the solution u ε is smooth enough over

Ω T := Ω × [0, T ], i.e. (A) u ε ∈ C [0, T ]; H 5 (Ω) ∩ C 2 [0, T ]; H 4 (Ω) ∩ C 3 [0, T ]; H 2 (Ω) ,
and there exist ε 0 > 0 and C 0 > 0 independent of ε such that

u ε L ∞ (0,T ;H 5 (Ω)) + ∂ 2 t u ε L ∞ (0,T ;H 4 (Ω)) + ∂ 3 t u ε L ∞ (0,T ;H 2 (Ω)) ≤ C 0 ,
uniformly in 0 ≤ ε ≤ ε 0 . Then there exist h 0 > 0 and τ 0 > 0 sufficiently small with h 0 ∼ √ εe -CT | ln(ε)| 2 and τ 0 ∼ √ εe -CT | ln(ε)| 2 such that, when 0 < h ≤ h 0 and 0 < τ ≤ τ 0 satisfying the stability condition (3.7), we have the following error estimates

e ε,k ≤ C 1 (ε, T )(h 2 + τ 2 ), 0 ≤ k ≤ T τ , e ε,k H 1 ≤ C 2 (ε, T )(h 2 + τ 2 ), u ε,k ∞ ≤ Λ + 1, (3.9)
where

Λ = u ε L ∞ (ΩT ) , C 1 (ε, T ) ∼ e CT | ln(ε)| 2 , C 2 (ε, T ) ∼ 1 ε e CT | ln(ε)| 2 and C depends on C 0 .
The error bounds in this Theorem show not only the quadratical convergence in terms of the mesh size h and time step τ but also how the explicit dependence on the regularization parameter ε. Here we remark that the Assumption (A) is valid at least in the case of taking Gaussian as the initial datum.

Define the error functions e ε,k ∈ X M as

(3.10) e ε,k j = u(x j , t k ) -u ε,k j , j ∈ T 0 M , 0 ≤ k ≤ T τ ,
where u ε is the solution of the LogSE (1.1) with Ω = (a, b). Combining Proposition 2.5 and Theorem 3.1, we immediately obtain (see an illustration in the following diagram):

u ε,k O(h 2 +τ 2 ) / / O(ε)+O(h 2 +τ 2 ) * * U U U U U U U U U U U u ε (•, t k ) O(ε) u(•, t k ) Corollary 3.2.
Under the assumptions of Proposition 2.5 and Theorem 3.1, we have the following error estimates

e ε,k ≤ C 1 ε + C 1 (ε, T )(h 2 + τ 2 ), e ε,k H 1 ≤ C 2 ε 1/2 + C 2 (ε, T )(h 2 + τ 2 ), 0 ≤ k ≤ T τ , (3.11)
where C 1 and C 2 are presented as in Proposition 2.5, and C 1 (ε, T ) and C 2 (ε, T ) are given in Theorem 3.1.

Error estimates. Define the local truncation error

ξ ε,k j ∈ X M for k ≥ 1 as ξ ε,k j = iδ c t u ε (x j , t k ) + 1 2 δ 2 x u ε (x j , t k+1 ) + δ 2 x u ε (x j , t k-1 ) -u ε (x j , t k ) ln(ε + |u ε (x j , t k )|) 2 , j ∈ T M , 1 ≤ k < T τ , (3.12)
then we have the following bounds for the local truncation error. Lemma 3.3 (Local truncation error). Under Assumption (A), we have

ξ ε,k H 1 h 2 + τ 2 , 1 ≤ k < T τ .
Proof. By Taylor expansion, we have

(3.13) ξ ε,k j = iτ 2 4 α ε,k j + τ 2 2 β ε,k j + h 2 12 γ ε,k j ,
where

α ε,k j = 1 -1 (1 -|s|) 2 ∂ 3 t u ε (x j , t k + sτ )ds, β ε,k j = 1 -1 (1 -|s|)∂ 2 t u ε xx (x j , t k + sτ )ds, γ ε,k j = 1 -1 (1 -|s|) 3 ∂ 4 x u ε (x j + sh, t k+1 ) + ∂ 4 x u ε (x j + sh, t k-1 ) ds.
By the Cauchy-Schwarz inequality, we can get that

α ε,k 2 = h M-1 j=1 |α ε,k j | 2 ≤ h 1 -1 (1 -|s|) 4 ds M-1 j=1 1 -1 ∂ 3 t u ε (x j , t k + sτ ) 2 ds = 2 5 1 -1 ∂ 3 t u ε (•, t k + sτ ) 2 L 2 (Ω) ds - 1 -1 M-1 j=0 xj+1 xj (|∂ 3 t u ε (x, t k + sτ )| 2 -|∂ 3 t u ε (x j , t k + sτ )| 2 )dxds = 2 5 1 -1 ∂ 3 t u ε (•, t k + sτ ) 2 L 2 (Ω) ds - 1 -1 M-1 j=0 xj+1 xj ω xj ∂ x |∂ 3 t u ε (x ′ , t k + sτ )| 2 dx ′ dωds ≤ 2 5 1 -1 ∂ 3 t u ε (•, t k + sτ ) 2 L 2 (Ω) + 2h ∂ 3 t u ε x (•, t k + sτ ) L 2 (Ω) ∂ 3 t u ε (•, t k + sτ ) L 2 (Ω) ds ≤ max 0≤t≤T ∂ 3 t u ε L 2 (Ω) + h ∂ 3 t u ε x L 2 (Ω) 2 ,
which yields that when h ≤ 1,

α ε,k ≤ ∂ 3 t u ε L ∞ (0,T ;H 1 (Ω)) .
Applying the similar approach, it can be established that

β ε,k ≤ 2 ∂ 2 t u ε L ∞ (0,T ;H 3 (Ω)) .
On the other hand, we can obtain that

γ ε,k 2 ≤ h 1 -1 (1 -|s|) 6 ds M-1 j=1 1 -1 ∂ 4 x u ε (x j + sh, t k+1 ) + ∂ 4 x u ε (x j + sh, t k-1 ) 2 ds ≤ 4h 7 M-1 j=1 1 -1 ∂ 4 x u ε (x j + sh, t k+1 ) 2 + ∂ 4 x u ε (x j + sh, t k-1 ) 2 ds ≤ 8 7 ∂ 4 x u ε (•, t k-1 ) 2 L 2 (Ω) + ∂ 4 x u ε (•, t k+1 ) 2 L 2 (Ω) ≤ 4 u ε 2 L ∞ (0,T ;H 4 (Ω)) , which implies that γ ε,k ≤ 2 u ε L ∞ (0,T ;H 4 (Ω))
. Hence by Assumption (A), we get

ξ ε,k τ 2 ∂ 3 t u ε L ∞ (0,T ;H 1 (Ω)) + ∂ 2 t u ε L ∞ (0,T ;H 3 (Ω)) + h 2 u ε L ∞ (0,T ;H 4 (Ω)) C0 τ 2 + h 2 .
Applying δ + x to ξ ε,k and using the same approach, we can get that

|ξ ε,k | H 1 τ 2 ∂ 3 t u ε L ∞ (0,T ;H 2 (Ω)) + ∂ 2 t u ε L ∞ (0,T ;H 4 (Ω)) + h 2 u ε L ∞ (0,T ;H 5 (Ω)) C0 τ 2 + h 2 ,
which completes the proof.

For the first step, we have the following estimates. Lemma 3.4 (Error bounds for k = 1). Under Assumption (A), the first step errors of the discretization (3.6) satisfy e ε,0 = 0, e ε,1

H 1 τ 2 .
Proof. By the definition of u ε,1 j in (3.6), we have

e ε,1 j = τ 2 1 0 (1 -s)u ε tt (x j , sτ )ds, which implies that e ε,1 τ 2 ∂ 2 t u ε L ∞ (0,T ;H 1 (Ω)) τ 2 , |e ε,1 | H 1 τ 2 ∂ 2 t u ε L ∞ (0,T ;H 2 (Ω)) τ 2 ,
and the proof is completed.

Proof. [Proof of Theorem 3.1] We prove (3.9) by induction. It follows from Lemma 3.4 that (3.9) is true for k = 0, 1.

Assume (3.9) is valid for k ≤ n ≤ T τ -1. Next we need to show that (3.9) still holds for k = n + 1. Subtracting (3.4) from (3.12), we get the error equations

(3.14) iδ c t e ε,m j = - 1 2 (δ 2 x e ε,m+1 j +δ 2 x e ε,m-1 j )+r ε,m j +ξ ε,m j , j ∈ T M , 1 ≤ m ≤ T τ -1,
where r ε,m ∈ X M represents the difference between the logarithmic nonlinearity

(3.15) r ε,m j = u ε (x j , t m ) ln(ε + |u ε (x j , t m )|) 2 -u ε,m j ln(ε + |u ε,m j |) 2 , 1 ≤ m ≤ T τ -1.
Multiplying both sides of (3.14) by 2τ (e ε,m+1 j + e ε,m-1 j ), summing together for j ∈ T M and taking the imaginary parts, we obtain for 1 ≤ m < T /τ ,

(3.16) e ε,m+1 2 -e ε,m-1 2 = 2τ Im(r ε,m + ξ ε,m , e ε,m+1 + e ε,m-1 ) ≤ 2τ r ε,m 2 + ξ ε,m 2 + e ε,m+1 2 + e ε,m-1 2 . Summing (3.16) for m = 1, 2, • • • , n (n ≤ T τ -1), we obtain e ε,n+1 2 + e ε,n 2 ≤ e ε,0 2 + e ε,1 2 + 2τ e ε,n+1 2 + 2τ n-1 m=0 ( e ε,m 2 + e ε,m+1 2 ) + 2τ n m=1 r ε,m 2 + ξ ε,m 2 . (3.17) For m ≤ n, when |u ε,m j | ≤ |u ε (x j , t m )|, we write r ε,m j as |r ε,m j | = e ε,m j ln(ε + |u ε (x j , t m )|) 2 + 2u ε,m j ln ε + |u ε (x j , t m )| ε + |u ε,m j | ≤ 2 max{ln(ε -1 ), | ln(ε + Λ)|}|e ε,m j | + 2|u ε,m j | ln 1 + |u ε (x j , t m )| -|u ε,m j | ε + |u ε,m j | ≤ 2|e ε,m j |(1 + max{ln(ε -1 ), | ln(ε + Λ)|}). On the other hand, when |u ε (x j , t m )| ≤ |u ε,m j |, we write r ε,m j as |r ε,m j | = e ε,m j ln(ε + |u ε,m j |) 2 + 2u ε (x j , t m ) ln ε + |u ε (x j , t m )| ε + |u ε,m j | ≤ 2 max{ln(ε -1 ), | ln(ε + 1 + Λ)|}|e ε,m j | + 2|u ε (x j , t m )| ln 1 + |u ε,m j | -|u ε (x j , t m )| ε + |u ε (x j , t m )| ≤ 2|e ε,m j |(1 + max{ln(ε -1 ), | ln(ε + 1 + Λ)|})
, where we use the assumption that u ε,m ∞ ≤ Λ + 1 for m ≤ n. Thus it follows that

r ε,m 2 | ln(ε)| 2 e ε,m 2 ,
when ε is sufficiently small. Thus when τ ≤ 1 2 , by using Lemmas 3.3, 3.4 and (3.17), we have

e ε,n+1 2 + e ε,n 2 e ε,0 2 + e ε,1 2 + τ n-1 m=0 ( e ε,m 2 + e ε,m+1 2 ) + τ n m=1 r ε,m 2 + ξ ε,m 2 (h 2 + τ 2 ) 2 + τ | ln(ε)| 2 n-1 m=0
( e ε,m 2 + e ε,m+1 2 ).

We emphasize here that the implicit multiplicative constant in this inequality depends only on C 0 , but not on n. Applying the discrete Gronwall inequality, we can conclude that

e ε,n+1 2 e CT | ln(ε)| 2 (h 2 + τ 2 ) 2 ,
for some C depending on C 0 , which gives the error bound for e ε,k with k = n + 1 in (3.9) immediately.

To estimate |e ε,n+1 | H 1 , multiplying both sides of (3.14) by 2(e ε,m+1 j -e ε,m-1 j

) for m ≤ n, summing together for j ∈ T M and taking the real parts, we obtain

|e ε,m+1 | 2 H 1 -|e ε,m-1 | 2 H 1 = -2 Re r ε,m + ξ ε,m , e ε,m+1 -e ε,m-1 = 2τ Im r ε,m + ξ ε,m , -δ 2 x (e ε,m+1 + e ε,m-1 ) = 2τ Im δ + x (r ε,m + ξ ε,m ), δ + x (e ε,m+1 + e ε,m-1 ) ≤ 2τ |r ε,m | 2 H 1 + |ξ ε,m | 2 H 1 + |e ε,m+1 | 2 H 1 + |e ε,m-1 | 2 H 1 . (3.18)
To give the bound for δ +

x r ε,m , for simplicity of notation, denote

u ε,m j,θ = θu ε (x j+1 , t m ) + (1 -θ)u ε (x j , t m ), v ε,m j,θ = θv ε,m j+1 + (1 -θ)v ε,m j ,
for j ∈ T M and θ ∈ [0, 1]. Then we have

δ + x r ε,m j = 2δ + x (u ε (x j , t m ) ln(ε + |u ε (x j , t m )|)) -2δ + x (u ε,m j ln(ε + |u ε,m j |)) = 2 h 1 0 [u ε,m j,θ ln(ε + |u ε,m j,θ |)] ′ (θ)dθ - 2 h 1 0 [v ε,m j,θ ln(ε + |v ε,m j,θ |)] ′ (θ)dθ = I 1 + I 2 + I 3 ,
where

I 1 := 2δ + x u ε (x j , t m ) 1 0 ln(ε + |u ε,m j,θ |)dθ -2δ + x u ε,m j 1 0 ln(ε + |v ε,m j,θ |)dθ, I 2 := δ + x u ε (x j , t m ) 1 0 |u ε,m j,θ | ε + |u ε,m j,θ | dθ -δ + x u ε,m j 1 0 |v ε,m j,θ | ε + |v ε,m j,θ | dθ, I 3 := δ + x u ε (x j , t m ) 1 0 (u ε,m j,θ ) 2 |u ε,m j,θ |(ε + |u ε,m j,θ |) dθ -δ + x u ε,m j 1 0 (v ε,m j,θ ) 2 |v ε,m j,θ |(ε + |v ε,m j,θ |)
dθ.

Then we estimate I 1 , I 2 and I 3 , separately. Similar as before, we have

|I 1 | ≤ 2|δ + x u ε (x j , t m )| 1 0 ln ε + |u ε,m j,θ | ε + |v ε,m j,θ | dθ + 2 δ + x e ε,m j 1 0 ln(ε + |v ε,m j,θ |) dθ = 2|δ + x u ε (x j , t m )| 1 0 ln 1 + |u ε,m j,θ | -|v ε,m j,θ | ε + min{|u ε,m j,θ |, |v ε,m j,θ |} dθ + 2 δ + x e ε,m j 1 0 ln(ε + |v ε,m j,θ |) dθ ≤ 2 ε |δ + x u ε (x j , t m )| |e ε,m j | + |e ε,m j+1 | + 2 δ + x e ε,m j max{ln(ε -1 ), | ln(ε + 1 + Λ)|} 1 ε |e ε,m j | + |e ε,m j+1 | + ln(ε -1 ) δ + x e ε,m j , and
|I 2 | = δ + x u ε (x j , t m ) 1 0 |u ε,m j,θ | ε + |u ε,m j,θ | - |v ε,m j,θ | ε + |v ε,m j,θ | dθ + δ + x e ε,m j 1 0 |v ε,m j,θ | ε + |v ε,m j | dθ ≤ |δ + x e ε,m j | + |δ + x u ε (x j , t m )| 1 0 ε|u ε,m j,θ -v ε,m j,θ | (ε + |u ε,m j,θ |)(ε + |v ε,m j,θ |) dθ ≤ |δ + x e ε,m j | + |δ + x u ε (x j , t m )| ε 1 0 |u ε,m j,θ -v ε,m j,θ |dθ |δ + x e ε,m j | + 1 ε |e ε,m j | + |e ε,m j+1 | .
In view of the inequality that

(u ε,m j,θ ) 2 |u ε,m j,θ |(ε + |u ε,m j,θ |) - (v ε,m j,θ ) 2 |v ε,m j,θ |(ε + |v ε,m j,θ |) = (u ε,m j,θ ) 2 -u ε,m j,θ v ε,m j,θ |u ε,m j,θ |(ε + |u ε,m j,θ |) + u ε,m j,θ v ε,m j,θ |u ε,m j,θ |(ε + |u ε,m j,θ |) - (v ε,m j,θ ) 2 |v ε,m j,θ |(ε + |v ε,m j,θ |) ≤ |u ε,m j,θ -v ε,m j,θ | ε + u ε,m j,θ (v ε,m j,θ ) 2 (u ε,m j,θ -v ε,m j,θ ) + εv ε,m j,θ (u ε,m j,θ |v ε,m j,θ | -|u ε,m j,θ |v ε,m j,θ ) |u ε,m j,θ ||v ε,m j,θ |(ε + |u ε,m j,θ |)(ε + |v ε,m j,θ |) ≤ 4|u ε,m j,θ -v ε,m j,θ | ε ,
we can obtain that

I 3 |δ + x e ε,m j | + 1 ε |e ε,m j | + |e ε,m j+1 | .
Thus we can conclude that

|δ + x r ε,m j | 1 ε |e ε,m j | + |e ε,m j+1 | + ln(ε -1 ) δ + x e ε,m j . Summing (3.18) for m = 1, 2, • • • , n (n ≤ T τ -1), we obtain |e ε,n+1 | 2 H 1 + |e ε,n | 2 H 1 ≤ |e ε,0 | 2 H 1 + |e ε,1 | 2 H 1 + τ n m=1 |r ε,m | 2 H 1 + |ξ ε,m | 2 H 1 + τ |e ε,n+1 | 2 H 1 + τ n-1 m=0 (|e ε,m | 2 H 1 + |e ε,m+1 | 2 H 1 ).
Thus when τ ≤ 1/2, by using Lemmas 3.3 and 3.4, we have

|e ε,n+1 | 2 H 1 + |e ε,n | 2 H 1 |e ε,0 | 2 H 1 + |e ε,1 | 2 H 1 + τ n m=1 1 ε 2 |e ε,m | 2 H 1 + |ξ ε,m | 2 H 1 + τ | ln(ε)| 2 n-1 m=0 |e ε,m | 2 H 1 + |e ε,m+1 | 2 H 1 e CT | ln(ε)| 2 ε 2 (h 2 + τ 2 ) 2 + τ | ln(ε)| 2 n-1 m=0 (|e ε,m | 2 H 1 + |e ε,m+1 | 2 H 1 ).
Applying the discrete Gronwall's inequality, we can get that

|e ε,n+1 | 2 H 1 e CT | ln(ε)| 2 (h 2 + τ 2 ) 2 /ε 2 ,
which establishes the error estimate for e ε,k H 1 for k = n + 1. Finally the boundedness for the solution u ε,k can be obtained by the triangle inequality

u ε,k ∞ ≤ u ε (•, t k ) L ∞ (Ω) + e ε,k ∞ ,
and the inverse Sobolev inequality [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF] e ε,k

∞ e ε,k H 1 ,
which completes the proof of Theorem 3.1. Case II: A general initial data, i.e. u 0 in (1.1) is chosen as

(4.2) u 0 (x) = tanh(x)e -x 2 , x ∈ R,
which is the multiplication of a dark soliton of the cubic nonlinear Schrödinger equation and a Gaussian. Notice that in this case, the logarithmic term ln |u 0 | 2 is singular at x = 0.

The RLogSE (2.1) is solved numerically by the SIFD (3.4) on domains Ω = [-12 , 12] and Ω = [-16, 16] for Case I and II, respectively. To quantify the numerical errors, we introduce the following error functions:

e ε (t k ) := u(•, t k ) -u ε (•, t k ), e ε (t k ) := u ε (•, t k ) -u ε,k , e ε (t k ) := u(•, t k ) -u ε,k , e ε E := |E(u) -E ε (u ε )|. (4.3)
Here u and u ε are the exact solutions of the LogSE (1.1) and RLogSE (2.1), respectively, while u ε,k is the numerical solution of the RLogSE (2.1) obtained by the SIFD (3.4). The 'exact' solution u ε is obtained numerically by the SIFD (3.4) with a very small time step, e.g. τ = 0.01/2 9 and a very fine mesh size, e.g. h = 1/2 15 . Similarly, the 'exact' solution u in Case II is obtained numerically by the SIFD (3.4) with a very small time step and a very fine mesh size as well as a very small regularization parameter ε, e.g. ε = 10 -14 . The energy is obtained by the trapezoidal rule for approximating the integrals in the energy (1.2) and (2.2). 4.1), which confirms the error bounds in Corollary 3.2.

Conclusion.

In order to overcome the singularity of the log-nonlinearity in the logarithmic Schrödinger equation (LogSE), we proposed a regularized logarithmic Schrödinger equation (RLogSE) with a regularization parameter 0 < ε ≪ 1 and established linear convergence between RLogSE and LogSE in terms of the small regularization parameter. Then we presented a semi-implicit finite difference method

Theorem 2 . 2 .

 22 Let λ ∈ R and ε > 0. Consider (1.1) and (2.1) on Ω = R d , or bounded Ω with homogeneous Dirichlet or periodic boundary condition. Consider an initial datum u

  3.5], established initially in [13, Lemme 1.1.1]. Lemma 2.4. Let ε ≥ 0 and denote f ε (z) = z ln(ε + |z|), then we have

4 .

 4 Numerical results. In this section, we test the convergence rate of the regularized model (2.1) and the SIFD(3.4). To this end, we take d = 1, Ω = R and λ = -1 in the LogSE (1.1) and consider two different initial data: Case I: A Gaussian initial data, i.e. u 0 in (1.1) is chosen as(4.1) u 0 (x) = 4 -λ/πe ivx+ λ 2 x 2 , x ∈ R, with v = 1.In this case, the LogSE (1.1) admits the moving Gausson solution (1.7) with v = 1 and b 0 = 4 -λ/π as the exact solution.

4. 1 .

 1 Convergence rate of the regularized model. Here we consider the error between the solutions of the RLogSE (2.1) and the LogSE (1.1). Fig. 4.1 shows e ε , e ε H 1 , e ε ∞ (the definition of the norms is given in (3.2)) at time t = 0.5 for Cases I & II, while Fig. 4.2 depicts e ε E (0.5) for Cases I & II and time evolution of e ε (t) with different ε for Case I. For comparison, similar to Fig. 4.1, Fig. 4.3 displays the convergent results from (2.8) to (1.1). From Figs. 4.1, 4.2 & 4.3 and additional numerical results not shown here for brevity, we can draw the following conclusions: (i) The solution of the RLogSE (2.1) converges linearly to that of the LogSE (1.1) in terms of the regularization parameter ε in both L 2 -norm and L ∞ -norm, and respectively, the convergence rate becomes O( √ ε) in H 1 -norm for Case II. (ii) The regularized energy E ε (u ε ) converges linearly to the energy E(u) in terms of ε. (iii) The constant C in (2.6) may grow linearly with time T and it is independent of ε. (iv) The solution of (2.8) converges at O( √ ε) to that of (1.1) in both L 2 -norm and L ∞ -norm, and respectively, the convergence rate becomes O(ε 1/4 ) in H 1 -norm for Case II. Thus (2.1) is much more accurate than (2.8) for the regularization of the LogSE (1.1). (v) The numerical results agree and confirm our analytical results in Section 2.

2 || e ε || H 1 Fig. 4 . 1 :Fig. 4 . 2 :

 214142 Fig. 4.1: Convergence of the RLogSE (2.1) to the LogSE (1.1), i.e. the error e ε (0.5) in different norms vs the regularization parameter ε for Case I (left) and Case II (right).

4. 2 . 1 || e ε || H 1 Fig. 4 . 3 :

 21143 Fig. 4.3: Convergence of the RLogSE (2.8) to the LogSE (1.1), i.e. the error e ε (0.5) in different norms vs the regularization parameter ε for Case I (left) and Case II (right).

2 Fig. 4 . 4 :

 244 Fig. 4.4: Convergence of the SIFD (3.4) to the RLogSE (2.1), i.e. errors e ε (0.5) vs τ (with h = 75τ /64) under different ε for Case I initial data.

* This work was partially supported by the Ministry of Education of Singapore grant R-146-000-223-112 (MOE2015-T2-2-146) (W. Bao).

for discretizing RLogSE and proved second-order convergence rates in terms of mesh size h and time step τ . Finally, we established error bounds of the semi-implicit finite difference method to LogSE, which depend explicitly on the mesh size h and time step τ as well as the small regularization parameter ε. Our numerical results confirmed our error bounds and demonstrated that they are sharp.


Jean-Michel Portal
email: jean-michel.portal@univ-amu.fr

Marc Bocquet

Santhosh Onkaraiah

Mathieu Moreau

Hassen Aziza

Damien Deleruyelle

Kholdoun Torki

Elisa Vianello

Alexandre Levisse

Bastien Giraud

Olivier Thomas

Design and Simulation of a 128 kb Embedded Nonvolatile Memory Based on a Hybrid RRAM (HfO_2 )/28 nm FDSOI CMOS Technology

Keywords: Embedded non-volatile memory, memory architecture, resistive switching memory, RRAM, I

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

two strong markets: the high performance, high density, high capacity memories for computation in servers, and the low energy, low cost embedded-memories devoted to autonomous connected nodes. The growth of these two markets is strongly related to each other and is at the basis of Internet of Things market. High performance computing applications for servers require high quantities of NVM [START_REF] Lynch | Big data: How do your data grow?[END_REF] in order to store and process tremendous amount of data generated by the autonomous connected nodes that are spread in the environment [START_REF] Gubbia | Internet of Things (IoT): A vision, architectural elements, and future directions[END_REF]. To this aim, high performance Volatile and Non-Volatile memories are used. In classic computing, the first elements of memory hierarchy are the volatile memories: SRAM and DRAM. The DRAM memory is used as a buffer between NVM and SRAM cache memories. At the end of the hierarchy are NVM memories HDD [START_REF] Shiroishi | Future Options for HDD Storage[END_REF] and NAND Flash technologies [START_REF] Helm | A 128Gb MLC NAND-Flash device using 16nm planar cell[END_REF]. The main constraints of these memories are (i) overall power consumption (from the DRAM refreshment to the NVM programming operations), (ii) communication cost between processing units and memories, and (iii) production costs per bit. Embedded memories for autonomous connected objects are more focused in reducing energy consumption and the overall chip cost. These device operations are based on classical memory hierarchy. Von Neumann or Harvard architectures are considered. SRAM memories are used to store the variables while NVM is used to store the instruction code. Embedded NVM are based on NOR flash technologies [START_REF] Ogura | A 90nm Floating Gate "B4-Flash" Memory Technology-Breakthrough of the Gate Length Limitation on NOR Flash Memory[END_REF] [6] [START_REF] Fastow | A 45nm NOR Flash Technology with Self-Aligned Contacts and 0.024µm2 Cell Size for Multi-level Applications[END_REF] down to 40nm commercial CMOS technologies. For both server applications (standalone NVM) and autonomous connected nodes (embedded NVM), scaling down the technology nodes is a real struggle. NAND flash technology complexity (air gap [START_REF] Seo | Highly reliable M1X MLC NAND flash memory cell with novel active air-gap and p+ poly process integration technologies[END_REF], vertical stacking [START_REF] Park | Three-Dimensional 128 Gb MLC Vertical nand Flash Memory With 24-WL Stacked Layers and 50 MB/s High-Speed Programming[END_REF], multiple patterning [START_REF] Pikus | Advanced multi-patterning and hybrid lithography techniques[END_REF]) and production costs are becoming unaffordable. On the embedded side, NOR flash is facing extremely complex structure [START_REF] Do | Scaling of split-gate flash memory and its adoption in modern embedded non-volatile applications[END_REF], and the co-integration with advanced CMOS nodes leads to reliability issues due to the high voltages needed for NVM programming operations. Moreover, scaling down floating gate technologies leads to

Design and Simulation of a 128kb Embedded

Non-Volatile Memory based on a Hybrid RRAM (HfO 2 ) / 28nm FDSOI CMOS Technology

Jean-Michel Portal, Marc Bocquet, Santhosh Onkaraiah, Mathieu Moreau, Hassen Aziza, Damien Deleruyelle, Kholdoun Torki, Elisa Vianello, Alexandre Levisse, Bastien Giraud, Olivier Thomas, Fabien Clermidy reduced endurance and retention [START_REF] Grupp | The bleak future of NAND flash memory[END_REF]. In order to continue the scaling of NVM technologies, for both Embedded and Standalone applications, emerging resistive switching technologies (i.e. RRAM) are extensively investigated [START_REF] Clermidy | Resistive memories: Which applications?[END_REF] because of their Back End of Line (BEoL) compatible structure and their high scalability [START_REF] Lee | Scaling trends and challenges of advanced memory technology[END_REF]. Among all the RRAM technologies, the Oxide-Based RRAM (OxRAM), because of their simple stack structure, fast switching, and common material appear as a promising solution for Flash memories replacement [START_REF] Benoist | Advanced CMOS Resistive RAM Solution as Embedded Non-Volatile Memory[END_REF]. Several metal oxides can be used to obtain the OxRAM behavior such as AlOx, NiOx, TaOx, TiOx or HfOx [16] [17] [18] [START_REF] Kwon | Atomic structure of conducting nanofilaments in TiO2 resistive switching memory[END_REF]. Integrating RRAMs as Flash memories replacement can be a solution either for Standalone Memories or for Embedded Memories. In the context of Standalone Memories, the supporting CMOS technology is totally dedicated to the memory integration and, as a consequence, the associated CMOS technology does not suffer from the memory requirements. In the context of Embedded Memories, CMOS technology is a limitation due to integration severe constraints. For example, in NOR Flash integration at advanced CMOS nodes, High Voltage Thick Oxide Transistors are integrated at the cost of additional masks, higher thermal budget and higher occupied area [START_REF] Tanzawa | High-voltage transistor scaling circuit techniques for high-density negative-gate channel-erasing NOR flash memories[END_REF].

The co-integration of RRAM at sub 30nm CMOS technologies leads to reliability issues on the Thin Oxide Transistors and is usually solved by using the IO transistors (thicker gate oxide process with no impact on the thin oxide transistor performances) [START_REF] Fackenthal | 19.7 A 16Gb ReRAM with 200MB/s write and 1GB/s read in 27nm technology[END_REF]. The direct impact is a higher maximum voltage, but also a larger footprint (up to 15 times for an equivalent drive). Several Embedded RRAM Memories are presented considering thick gate oxide transistors, in [START_REF] Sheu | A 4Mb embedded SLC resistive-RAM macro with 7.2ns read-write random-access time and 160ns MLC-access capability[END_REF] a 4Mb HfO2-based RRAM memory in 180nm CMOS technology with write verify technic is demonstrated, in [START_REF] Chang | A High-Speed 7.2-ns Read-Write Random Access 4-Mb Embedded Resistive RAM (ReRAM) Macro Using Process-Variation-Tolerant Current-Mode Read Schemes[END_REF] a 4Mb embedded RRAM memory is demonstrated with process variation control circuits in 130nm CMOS technology. In [START_REF] Chang | A Low-Voltage Bulk-Drain-Driven Read Scheme for Sub-0.5 V 4 Mb 65 nm Logic-Process Compatible Embedded Resistive RAM (ReRAM) Macro[END_REF] a 4Mb macro using a unipolar RRAM is embedded in 65nm CMOS technology. A 1Mb RRAM memory embedded 28nm bulk CMOS is demonstrated in [START_REF] Chang | Embedded 1Mb ReRAM in 28nm CMOS with 0.27-to-1V Read Using Swing-Sample-and-Couple Sense Amplifier and Self-Boost-Write-Termination Scheme[END_REF] with write assist circuits. However, in these papers the reliability issue on MOS transistors during programming operations and during the forming step is not considered. Others studies using thick oxide transistors and Conductive Bridge Random Access Memories (CBRAM) [START_REF] Fackenthal | 19.7 A 16Gb ReRAM with 200MB/s write and 1GB/s read in 27nm technology[END_REF] or Phase Change Memories (PCM) [START_REF] Borghi | r. a. t. a. 1. w. throughput[END_REF] are reported.

In [START_REF] Shen | High-K metal gate contact RRAM (CRRAM) in pure 28nm CMOS logic process[END_REF], a testchip with 28nm thin-gate oxide bitcells was presented and the MOS transistors reliability is shown for a 3V programming voltage. This paper presents the architecture of an embedded 128kb memory cut based on a hybrid RRAM (HfO 2 ) and thin-gate oxide 28nm Fully Depleted Silicon On Insulator (FDSOI) CMOS technology suitable for NOR Flash replacement. Validation of the proposed architecture is performed through post-layout simulation by using RRAM compact model calibrated on CEA-Leti RRAM samples [START_REF] Vianello | Resistive Memories for Ultra-Low-Power embedded computing design[END_REF] and implemented under transistor-level simulator Eldo [START_REF]Mentor Graphics Website[END_REF] and CMOS 28nm FDSOI design kit from STMicroelectronics.

The rest of the paper is organized as follows: In the next section, the RRAM technology is briefly introduced together with the compact model calibrated on state of the art devices. In section III, the bit-cell and the memory array organization are described with the different modes of operation (FORMING/SET, RESET and READ). Section IV is dedicated to the full memory macro-cell description, including peripheral circuits, addressing hierarchy and scheduler. The validation of the macro-cell through simulation is presented in Section V. Finally, Section VI gives some concluding remarks and highlights on the proposed embedded memory.

II. RRAM OVERVIEW: TECHNOLOGY AND COMPACT MODEL

RRAM based on HfO 2 are studied as part of the bit-cell of the proposed memory architecture. The RRAM stack is composed of a 5 nm thick HfO 2 resistive switching layer embedded in-between a TiN/Ti Top Electrode (TE) and a TiN Bottom Electrode (BE). The resistive switching layers are deposited by Atomic Layer Deposition (ALD), whereas the metallic electrodes are deposited by Physical Vapor Deposition (PVD) [START_REF] Vianello | Resistive Memories for Ultra-Low-Power embedded computing design[END_REF].

RRAM modeling used in this study is based on the work presented in [START_REF] Bocquet | Robust compact model for bipolar oxidebased resistive switching memories[END_REF]. This approach relies on electric fieldinduced creation/destruction of oxygen vacancies within the switching layer. The model enables continuous accounting for FORMING, SET and RESET operations into a single master equation in which the resistance is controlled by the radius of a conductive filament (namely r CF ).

After calibration, the model satisfactorily matches quasistatic and dynamic experimental data measured on actual HfO 2 -based memory elements. Moreover, to account for the variability of RRAM technology, two corners cases were simulated. They include the two extreme behaviors observed experimentally: one favoring the SET mechanism and slowing the RESET, and the other being the exact opposite.

The Fig. 1 show the quasi-static behavior of RRAM devices and the good modeling correlation. The corners encompass the full range of features that ensure to take into account the worst case of FORMING, SET and RESET. The Fig. 2.a describes the dependence of the switching time according to the amplitude of the programming pulse voltage. It is important to note that rapid operations -low consumptions -need higher CMOS standard voltage. Furthermore, the Fig. 2.b. underlines the interest of high voltage applied during the RESET to ensure a high resistance value. These two central behaviors are perfectly captured by the model implemented and extreme behaviors are included within in the corner simulations.

III. BIT-CELL AND MEMORY ARRAY ORGANIZATION

In this section, the bit-cell structure is detailed together with memory array organization. Based on the array arrangement, biasing conditions for selected and unselected cells in the array are discussed for the different modes of operation.

A. Bit-cell structure

The bit-cell is based on a 2T1R structure [START_REF] Jovanovic | Design considerations for reliable OxRAM-based non-volatile flipflops in 28nm FD-SOI technology[END_REF], as described in Fig. 3.a and exhibits one NMOS and one PMOS to access the Top Electrode (TE) of the RRAM, whereas the Bottom Electrode (BE) is connected to the Reset Line (RL). In addition to the classical Bit-Line (BL) and Word-Line (WL), two others lines are present namely Reset-Line (RL) and Set-Line (SL). As depicted Fig. 3.b, to avoid any CMOS core process modification, the RRAM elements are introduced in the Metal Insulator Metal (MIM) stack between the first and second 8x metallization level in place of the decoupling capacitance. On one hand, this solution is fully compatible with the standard CMOS 28nm FDSOI process flow but on the other hand this integration scheme is achieved at the cost of a larger RRAM foot-print on the upper metallization level.

The additional lines RL and SL are used to perform FORMING/SET and RESET/READ operations on two different paths, as illustrated Fig. 4. Doing so FORMING/SET voltages are applied on the top electrode of the RRAM through the PMOS (Fig. 4.a), whereas ground (gnd) is applied on the top-electrode through a NMOS and RESET voltage is applied on the RL during a RESET operation (Fig. 4.b). With this scheme, degradation of the voltage level for the different operations is reduced, since positive voltages are applied through PMOS transistor whereas ground voltages are applied through NMOS. Moreover, the compliance current, during FORMING / SET operations, is defined by the sizing of the PMOS transistor and by the V DDM biasing. Given that the current flowing through the cell during the RESET operation must be above the FORMING / SET current, the NMOS transistor has to be larger than the PMOS for an equivalent biasing.

B. Memory array organization

The arrangement of the memory array is given in Fig. 5.a with the schematic view and Fig. 5.b with the corresponding layout view. The RL and WL lines are shared per row whereas BL and SL are shared per column. To access a cell for a given operation, biasing of all the four access lines is mandatory, whereas inhibition voltages must be applied on unselected cells. Depending, on the operation unselected cells are inhibited by turning OFF the access transistor (NMOS and PMOS) with V GS below the threshold voltage or by having a voltage difference on the RRAM close to zero volt. Fig. 6 summarizes the biasing of the different access lines for four cases: selected bit-cell, unselected bit-cell on the same row, unselected bit-cell on the same column, unselected bit-cell on the rest of the array without any common lines with the selected cell. Three potentials are used, V DD equal to V DDM in the array, ground (gnd) and a high voltage (HV). HV takes a value close to 2.5×V DDM for the FORMING operations and a value close to 2×V DDM for the SET and RESET operations. Fig. 6 depicts a single cell access, but this memory array organization enables to perform all the programming operations by selecting an entire or partial row or a full or partial array. This feature offers the capability to have the best compromise between speed and consumption.

In the proposed memory macro-cell, the array is defined by the bank size. A bank is composed of 1024 cells organized on 32 rows per 32 columns, in other words, 32 words of 32 bits.

To program a word, a RESET operation is first performed on the bits to reset and followed by a SET operation on the remaining bits of the word. Similarly to NOR Flash memory, two programming phases are used to write a word. But, in our case, both operations are selective contrary to Flash memory, where erase operation is applied to all bits in the word followed by a selective write operation. In this section the full macro-cell architecture is detailed, starting with the peripheral blocks to address the bank array, up to the full macro-cell hierarchy including scheduler finite state machine to generate timing and internal signals. It is worth to note that the macro-cell communication bus is fully compatible with AMBA 3 AHB lite protocol [START_REF]AMBA 3 AHB-Lite Protocol Specification Documentation v1.0[END_REF].

A. Peripheral block description

The macro-cell architecture is massively multi-bank, thus all peripheral circuits to program and read the content of the bit-cells are introduced at bank level. In this subsection, levelshifters used for programming operations as well as sense amplifier used for reading operation are detailed. Banks are powered with two power supply V DD =V DDana and a higher voltage for FORMING/SET/RESET operations named HV. During programming operation, level shifters are used to drive HV on RL and SL access lines, while classic buffers are used to drive gnd/V DDana on BL and WL access lines. Table I gives a summary of possible biasing of the different access lines. The architecture of level shifters acting on SL and RL are represented in Fig. 7 and Fig. 8 respectively with their schematic and layout views. Level shifter structures are designed with cascade MOS since the voltage difference between any MOS transistor nodes has to be below or equal to VDD to ensure reliability and avoid gate-oxide breakdown.

Only bit-cells on a common row are activated at a time during programming operation, thus since SL level shifters are shared per column, they only have to drive single cell FORMING/SET compliance current, limiting the sizing of their output stage.

On the contrary, RL level shifters are shared per row, thus during a RESET operation, they may drive up to the 32 bitcells of the addressed word. Thus, the sizing of the RL level shifter output stage is able to drive 32 times the current of a single cell. It is to be noted on the layout view of the RL level shifter (Fig. 8.b), that a single output buffer is multiplied 32 times to form the complete output stage. Similar to the SL level shifter, the input of the RL level shifter is driven by logic gate biased in the V DD domain. In order to sense the value of each bit-cells, sense amplifiers are added on top of the corresponding columns. Applying a current I READ through the resistive element and comparing the resulting voltage V READ with a voltage reference V REF gives the logic value of the bit-cell. I READ is generated thanks to a Wilson current mirror structure. An operational full swing amplifier is used to discriminate between V READ and V REF .

Since, RRAM technology is still in development, the reference voltage V REF is provided externally for characterization purpose. Indeed, by trimming V REF value, knowing I READ , it is possible to extract the resistance value of all the bit-cells in the array and thus extract Low Resistance State (LRS) and High Resistance State (HRS) distributions on chip. Doing so at the end of the characterization procedure, V REF value is set between the voltage distribution corresponding to LRS and HRS distributions. Moreover, an external pad gives a direct access to the bit-cell content to extract RRAM resistance value in order to verify the value extracted from V REF trimming. This pad can be connected to the sense amplifier of the first column (BL 0 ) of each bank. A bank is composed of 32 by 32 bit-cells to avoid voltage drop on the access line. Indeed, since the whole circuit is designed with GO1 devices the voltage budget is limited, especially during the FORMING steps. Thus, all peripheral circuits, i.e. level shifters, sensing circuitry are implemented at bank level, as depicted Fig. 10 with the layout view. Moreover, addressing signal gates all control signals at bank level. Thus, unselected bank are completely inactive, with no internal signal change, preventing any disturb or extra powerconsumption.

B. Bank organization and addressing hierarchy

For the selected bank an acknowledgement signal is generated to properly apply all control signal and to avoid any temporal drift between signals due to hierarchy routing. The top-analog circuit is divided into four sectors of eight pages and each page contains four banks. Control signals, generated by the digital scheduler, are enabled at each level of the hierarchy. Moreover, buffers are inserted on all input signals (data, address and control) and tri-state buffer are inserted on the data-out signals, to prevent delay issue.

C. Scheduler description

The scheduler is a Finite State Machine (FSM) compatible with AMBA 3 AHB lite protocol, which generates all necessary timing and internal signal to drive the top analog circuit. The timing can be trimmed since a programmable timer gives the time reference. Indeed, the targeted bus clock is in the range of a few ten's of MHz, in other words a few hundreds of 'ns' period, while FORMING/RESET/SET operations are in a range of 'ms' to 'µs' depending on the targeted HV voltage and RRAM technology variability. Thus, to enhance yield, timing has to be programmable in order to fit the voltage/duration dependency of the RRAM technology for the outlier cells of the memory array.

The scheduler states are Ready, Read, Write. The Write state is decomposed in a second FSM, since a write operation can be of two kinds FORMING or RESET/SET. Moreover, a write/verify procedure is embedded in the scheduler. After a write, which is a RESET operation followed by a SET operation, a read is performed and data are compared to data in, in case of mismatch, up to 10 cycles of RESET/SET operations can be applied. After 10 cycles, if Write operation still fails the HRESP signal indicates an error. This write/verify procedure has been implemented to tackle cycleto-cycle variability of the RRAM technology.

The full layout of the 128kb Non-Volatile Memory based on a Hybrid RRAM (HfO2) / 28nm FDSOI CMOS Technology is given Fig. 11 and its main simulated features are summarized Table II • READ sequentially each addressed word of the bank to retrieve the inverse checkerboard pattern. Fig. 12 shows the FORMING operation of a full bank, all the words are sequentially formed, it is worth to note that the FORMING process duration for a word is 100µs, thus for a full bank, it represents 3.2ms. The Fig. 12 highlights the voltage, current and Conductive Filament (CF) radius variation for the first bit-cell of the first word (Cell0,0) and for the last bit-cell of the last word (Cell31,31). For all the bitcell, the compliance current during FORMING is set to 62µA, representing an overall current of nearly 2mA for a full word. The evolution of the CF radius clearly shows the FORMING operation. The biasing conditions extracted from the simulation are:

• All the BL and WL remains to VDDana=1.2V during the FORMING process, • SL is set to HVFORMING=2.75V, instead of 2.8V due to voltage drop in the SL level-shifter, • RL is set to 5mV, instead of ground due to voltage drop in the RL level-shifter for the selected word, whereas it is set to HV FORMING for unselected words. Thus the selected RRAM are biased to 2.75V during FORMING. One can notice a very low leakage on the already formed cell of 0.76µA (see I CELL (cell 0,0 ) during the FORMING of the word 31 in Fig. 12). Fig. 13 shows simulation results for the overall READ/PROGRAM process on a bank (READ to verify FORMING, PROGRAM checkerboard pattern, READ to verify checkerboard pattern, PROGRAM inverse checkerboard pattern, READ to verify inverse checkerboard pattern). Functionality is validated, since the programmed patterns are successfully read. Fig. 14 and Fig. 15 highlight respectively a PROGRAM operation and a READ operation. The PROGRAM operation is divided into two steps, in a first step bit-cells corresponding to DATA_IN='1' are RESET, whereas in a second step bitcells corresponding to DATA_IN='0' are SET. On Fig. 14 DATA_IN 0 ='1', thus the cell 0,0 is RESET, with a current of 87µA for a voltage on the RRAM of 2.32-0.11 = 2.21V. In the SET step, the cell remains unchanged. Each steps, RESET and SET, takes 80ns for a global PROGRAM time of 170ns per word. Global current consumption is in a range of 1.5mA (SET all the cell) to 2.9mA (RESET all the cell). It is important to notice that the RESET current, as shown in Fig. 14, is not constant during the RESET phase. Fig. 15 exhibits the two steps of a READ operation, in the first step the current source is enabled to pre-charge the bit-line for 40ns, during the second step the sense amplifier is activated for 2ns to differentiate between HRS and LRS state. The current through the cell during the READ operation is below 2µA. The READ operation duration together with the voltage on the RRAM (worst case 0.8V for a HRS RRAM) allow to avoid any disturb on the cell. This assumption is validated since there is no CF radius change during the READ process presented Fig. 15.

The same simulations are performed for the corner cases obtained from the worst-case measurements on the RRAM cell as depicted in Section II together with the SS corner of the CMOS core process. The To guarantee the functionality of the macro-cell in the worstcase scenarios, timing has to be adapted, due to the voltage/timing characteristic of the RRAM. Doing so for nominal voltage configuration, timings and energy per cell are strongly degraded versus typical simulation as represented Fig. 16 (Typical (V 1 ) versus Corners (V 1 )). To recover typical results at nominal voltages for worst-case scenarios, voltages have to be increased as represented Fig. 16 by V 2 and V 3 voltage sets. However, the maximal stress voltage defined as gate/source or gate/drain voltage ramps up to 1.8V for Moreover, it is interesting to note that the best efficiency is achieved with higher voltages but this track of optimization must be carefully used considering reliability issue of the standard logic MOS devices. Since January 2015 he is the project leader of Silicon Impulse an IC competence center helping companies to design innovative products based on the latest low-power semiconductor technologies. He is author or co-author of 75 articles in international refereed journals and conferences and 25 patents.

Figure 1 :

 1 Figure1: Experimental I(V) characteristics for Electroforming, Set, and Reset measured on a large number of memory elements reflecting the device-todevice variability presented in[START_REF] Vianello | Resistive Memories for Ultra-Low-Power embedded computing design[END_REF] and simulation results including corners definition.

Figure 2 :

 2 Figure 2: Experimental (a) switching time for Forming, Set and Reset operations as a function of voltage and (b) RHRS as a function of stop voltage during Reset operation presented in [28] and corresponding simulation results.

Figure 3 :

 3 Figure 3: (a) Schematic view of the bit-cell with 2T1R and four access lines (BL, WL, SL and RL) (b) Layout of the bit-cell, with a RRAM element introduces in the Metal Insulator Metal (MIM) stack to avoid CMOS process modification at the cost of large cell area.

Figure 4 :

 4 Figure 4: (a) FORMING/SET are performed through the PMOS, SL is set to the positive operation voltage, RL is grounded whereas NMOS transistor is off (b) RESET is performed through NMOS, RL is set to the positive operation voltage, BL is grounded whereas PMOS transistor is off. For READ operation, BL is set to the read voltage and RL is grounded, here also the PMOS transistor is off. To reduce voltage degradation, two different paths are used through PMOS and NMOS, in order to fit the bipolar behavior of the RRAM.

Figure 5 :

 5 Figure 5: (a) Schematic view of a two by two bit-cell array (b) Layout view of a two by two bit-cell array. SL and BL lines are shared in column; RL and WL are shared in row.

Figure 6 :

 6 Figure 6: (a) Array biasing condition for a FORMING/SET operation performed on a single cell (b) Array biasing condition for a RESET operation performed on a single cell (c) Array biasing condition for a READ operation performed on a single cell. Biasing condition of the selected cell can be extended to any cells on a row to perform parallel operations. Biasing conditions can also be applied on the full array to perform global operation.

  MACRO-CELL ARCHITECTURE OVERVIEW

Figure 7 :

 7 Figure 7: (a) Schematic view of the level shifter used to drive SL access line (b) Layout view of the level shifter used to drive SL access line. The SL level shifter output swing is defined between VDDM and HV. The level shifter input is controlled with standard voltage logic (VDD domain).

Figure 8 :

 8 Figure 8: (a) Schematic view of the level shifter used to drive RL access line (b) Layout view of the level shifter used to drive RL access line. The RL level shifter output swing is defined between gnd and HV. The level shifter input is controlled with standard voltage logic (VDD domain).

Figure 9 :

 9 Figure 9: (a) Schematic view of the sensing circuitry connected to a bit-line including read current generation, output direct access and operational amplifier (b) Layout view of the 32 sensing circuitry associated to the 32 bitline of a single bank.

Fig. 9 .

 9 Fig.9.a illustrates the sense amplifier architecture, including CDMA standing for Current Direct Memory Access to provide I OUT on the first column of each bank from an external PAD for characterization purpose, doing so internal current mirror is disconnected. Finally, it is worth noting that internal current mirror as well as operational amplifier can be disconnected from V DD when no sensing operation is required. The layout of the sense amplifier together with current mirror is given Fig.9.b.

Figure 10 :

 10 Figure 10: Layout view of a full bank including a 32 by 32 bit-cells with surrounding peripherals circuits, i.e. 32 SL level shifters below the array, 32 BL drivers and 2 sensing circuitry on top of the array and 32 RL level shifter and 32 WL drivers on the right of the array.

Figure 11 :

 11 Figure 11: Layout view of the full macro-cell with top analog (right) and scheduler (left).

  .

Figure 12 :

 12 Figure 12: Simulation results (corner Typical) of the FORMING process for a full bank. Current through RRAM of the cell0,0 (Word 0, Bit 0) and of the cell31,31 (Word 31, Bit 31), CF radius evolution during forming, as well as access-line voltage to the cell0,0 and cell31,31 are plotted.

Figure 13 :

 13 Figure 13: Simulation results (corner Typical) of the full PROGRAM/READ process for a full bank. CF radius evolution of the cell0,0 RRAM (Word 0, Bit 0) and of the cell31,31 RRAM (Word 31, Bit 31), during PROGRAM&READ, as well as access-line voltage to the cell0,0 and cell31,31 are plotted. Finally, the similarity of the DATA_IN and DATA_OUT values shows the success of PROGRAM and READ operations.

Figure 14 :

 14 Figure 14: Simulation results (corner Typical) of a PROGRAM operation on cell0,0. Current and CF radius of the cell0,0 RRAM (Word 0, Bit 0) and of the cell31,31 RRAM (Word 31, Bit 31), during PROGRAM, as well as access-line voltage to the cell0,0 and cell31,31 are plotted. The program operation is composed of two steps, RESET/SET.

Fig. 16 .

 16 summarizes the time and energy per cell for different operations, comparing typical (CMOS & RRAM TT) and worst-case results using the corners case (CMOS SS & RRAM corner Slow FORMING/SET and Slow-corner RESET) for the nominal voltages, defined as V1 (HV FORMING =2.8V, HV SET/RESET =2.4V) and for two others voltages V2 (HV FORMING =3.0V, HV SET/RESET =2.6V) and V3 (HV FORMING =3.2V, HV SET/RESET =2.8V).

Figure 15 :

 15 Figure 15: Simulation results (corner Typical) of a READ operation for a full word. CF radius evolution of the cell0,0 RRAM (Word 0, Bit 0) and of the cell31,31 RRAM (Word 31, Bit 31), during PROGRAM&READ, as well as access-line voltage to the cell0,0 and cell31,31 are plotted. Pre-charge and sense phase are detailed.

Figure 16 :

 16 Figure 16: Operation duration and energy per cell for typical and corners cases, for 3 different voltage sets.

FORMING and 1 .

 1 4V for SET/RESET with set V 2 and respectively 2V and 1.6V with set V 3 . These voltages remain acceptable for the FDSOI 28nm standard logic CMOS technology. VI. CONCLUSION In this paper, a full 128kb Embedded Non-Volatile Memory based on a Hybrid RRAM (HfO 2 ) & 28nm FDSOI CMOS Technology is presented. The key points of the architecture are the use of standard logic MOS exclusively, avoiding any high voltage MOS usage, program/verify procedure to mitigate cycle to cycle (C2C) variability and direct bit-cell read access for characterization purpose. This architecture has been fully validated through an extensive set of postlayout simulations at different voltage levels and using typical and the most pessimistic corners (MOS SS and RRAM worst FORMING/SET and RESET corners). The memory is functional at all corners for the different set of voltages owing to the time/voltage dependencies of RRAM.

  25 papers in international conferences and journals. He is the main author of a book chapter and the main inventor or coinventor of 15 patents. Fabien Clermidy obtained his Ph.D in microelectronic from INPG, Grenoble in 1999 and his supervisor degree in 2011. He is a pioneer in designing Network-on-Chip based multicore. He was the leader of the second generation of Networkon-Chip based multicore dedicated to 3GPP-LTE. At this period, his team elaborated one of the first 3D multi-core prototypes embedding a WIDE-IO DRAM memory called WIOMING. He is currently managing the digital circuit laboratory implied in the development of new architectures using emerging technologies such as 3D TSV, 3D monolithic integration and emerging memories. He has published 2 books, more than 75 journal and conferences papers and is author or co-author of[START_REF] Benoist | Advanced CMOS Resistive RAM Solution as Embedded Non-Volatile Memory[END_REF] patents. Alexandre Levisse received his B.S. (Electrical Engineering) degree in 2012 and his M.S. (Electrical Engineering) degree in 2014, both from Aix-Marseille University, France. He is currently PhD student in CEA-Leti (Grenoble, France) and IM2NP (Aix-Marseille University, France). His research interests include emerging resistive memories with emphasis on circuit design, architecture and crossbar architecture. Elisa Vianello received the Ph.D. degree in microelectronics from the University of Udine, Udine, Italy, and the Polytechnic Institute of Grenoble, Grenoble, France, in 2009. She has been a Scientist with CEA-LETI, Grenoble, since 2011. Olivier Thomas received the M.S. Electrical Engineering degree in 2001 and the Ph.D. degree in microelectronics in 2004. He joined the CEA-LETI Laboratory in the Center for Innovation in Micro & Nanaotechnology (MINATEC), Grenoble, France in 2005. He was first involved in the development of low-power and low-leakage design solutions for digital wireless applications in 65nm Partially-Depleted SOI technology in collaboration with STMicroelectronics. From 2006 to 2010, he was in charge of low power SRAM and Digital design projects in Thin Film SOI technologies. His work was focused on efficient and simple multiple-VT design solutions. From 2010 to 2012, he was a visiting researcher at Berkeley Wireless Research Center (BWRC) of University of California at Berkeley. He worked on methodologies to characterize on large-scale static/dynamic SRAM performances. Back to CEA-LETI, from 2012 to 2014, he launched and led a advanced memory design group at LETI.

TABLE I VOLTAGE

 I SWING ON THE BIT-CELL ACCESS LINES

	Access lines	Voltage swing
	WL	gnd to VDDana
	BL	gnd to VDDana
	SL	VDDana to HV
	RL	gnd to HV

  scaled Flash memories. In 2005 he joined the Memory group of Im2np (Institut Matériaux Microélectronique Nanosciences de Provence) and became associate professor at Aix-Marseille Université. His research topics include nanoscale electrical characterization by scanning probe microscopy and physical modeling of emerging memory devices such as RRAM. In 2016 he became Professor at the Institut des Nanotechnologies de Lyon (INSA de Lyon) where he currently works on plastic electronics. Hassen Aziza received his B.S. and M.S. degrees in Electrical Engineering, both from University of Marseille, France. He received his Ph.D. degree in 2002 from the University of Marseille, France. Hassen Aziza is currently associate professor at Aix-Marseille University-IM2NP laboratory (Institute of materials, microelectronics and nanosciences of Provence). His research fields cover design, test and reliability of conventional non-volatile memories (Flash & EEPROM) as well as emerging memories (Resistive RAM). He is (co)author of more than 90 papers in international conferences and journals and is (co)inventor of 4 patents. Marc Bocquet received the Ph.D. degree in micro and nanoelectronics from University Grenoble, Grenoble, France, in 2009. He became an Associate Professor with the University of Marseille -Polytech'Marseille, in 2010, and he is a member of the Memories Team, IM2NP. He has conducted several studies on understanding the physical mechanisms in dielectrics to link the physical and chemical properties to the electricalperformance/reliability of memory devices: flash and resistive memory. Kholdoun Torki received the Ph.D. degree in Microelectronics from the Institut National Polytechnique de Grenoble, France, in 1990. He joined CMP as senior engineer in 1990, later on joining CNRS/CMP in 1994. He is currently Technical Director at CMP since 2002. His research interest includes Deep Submicron design methodologies, Non Volatile Memory CMOS co-integration, and 3D-IC integration. He authored and co-authored more than 100 scientific papers, coauthored 2 patents, designed more than 30 ASIC circuits, and participated or coordinated 15 European and National projects. He is member of the Board of Directors at iRoC Technologies. Bastien Giraud received the Ph.D. degree in 2008 from Telecom ParisTech France. The PhD thesis focused on SRAM design in Double Gate FDSOI. In 2009, he was postdoctoral researcher at UC Berkeley working on low power circuits and SRAM variability. From 2010, he works at CEA/Leti as a circuit designer specialized in memory and low power circuit in advanced technologies. His research interests include resilient memory with assist technics, energy efficiency, specific design technics and non-volatile memories. His current research are focused on SRAM ULV and robust, smart CAM, logic in memory, crossbar, using advanced CMOS technologies and non-volatile RRAM technologies such as CBRAM, OxRAM and PCRAM. He has published more than

Jean-Michel

Damien Deleruyelle received the Ph.D. degree in micro and nanoelectronics from Aix-Marseille University in 2004 after a thesis carried out at CEA-Leti (Grenoble, France) on ultra-