N Ben Abdallah

A Faraj

An improved transient algorithm for the simulation of a resonant tunneling diode *

Keywords: Nonlinear Schrödinger equation, numerical scheme, resonant tunneling diode, resonant states, time-dependent. Subject classifications: 35Q40, 35Q55, 65M06, 65Z05, 81-08, 81V99

A fast algorithm is proposed for the simulation of the transient evolution of a resonant tunneling diode by a multiscale approach. The problem is modeled by the time-dependent Schrödinger-Poisson system. By decomposing the wave function into a non resonant part and a resonant part, the fast algorithm is designed by combining standard finite difference method for the Schrödinger-Poisson equation with proper time-dependent and/or nonlinear transmission boundary condition. In addition, with a suitable interpolation of the non resonant one, an accurate and fast algorithm is presented for the computation of the resonant part via the projection method. With this two scale decomposition, the new numerical method can save the computational time significantly.

Introduction

The Resonant Tunneling Diode (RTD) is a heterostructure made of a superposition of thin layers of semiconductor of different types used to create double or multiple potential barriers. The possibility of resonant tunneling through multiple barriers and the resulting negative differential resistance were firstly discussed in [START_REF] Iogansen | The possibility of resonance transmission of electrons in crystals through a system of barriers[END_REF][START_REF] Tsu | Tunneling in a finite superlattice[END_REF]. Such physical effects make the RTD useful in the conception of signal generators, detectors and mixers, multi-valued logic switches, low-power amplifiers, local oscillators, frequency locking circuits, and also as generators of multiple high frequency harmonics, see [START_REF] Mennemann | Transient Schrödinger-Poisson simulations of a highfrequency resonant tunneling diode oscillator[END_REF][START_REF] Mizuta | The Physics and Applications of Resonant Tunneling Diodes[END_REF]. In the mathematical and numerical models for the RTD, it is natural to incorporate the quantum effects described above. In particular, due to the small size of the device, the kinetic models for semi-conductors [START_REF] Greengard | A boundary-value problem for the stationary Vlasov-Poisson equations: The plane diode[END_REF][START_REF] Poupaud | Boundary value problems for the stationary Vlasov-Maxwell system[END_REF] have to be replaced by their natural extension to the quantum level: the Wigner equation or equivalently a system of infinitely many Schrödinger equations. In the last models, the electron interaction can be represented thanks to a Poisson equation, which leads to the Wigner-Poisson [START_REF] Brezzi | The three-dimensional Wigner-Poisson problem: Existence, uniqueness and approximation[END_REF][START_REF] Lions | Sur les mesures de Wigner[END_REF] or Schrödinger-Poisson system [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF][START_REF] Ben Abdallah | A mathematical model for the transient evolution of a resonant tunneling diode[END_REF][START_REF] Nier | The dynamics of some quantum open systems with short-range nonlinearities[END_REF]. Both approaches are able to recover the negative differential resistance in the I-V characteristics, however, it was shown in [START_REF] Jona-Lasinio | On Schrödinger Equations with Concentrated Nonlinearities[END_REF][START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF] that the second approach could lead to hysteresis phenomena in agreement with physical observations. For the Schrödinger-Poisson system, which will be adopted in the present work, the nonlinear coupling effects take place in the middle of the structure. This active region is out of equilibrium and connected to the exterior through access zones which allow the injection of electrons. Therefore, a first important issue is to find suitable transparent boundary conditions for the Schrödinger equation at the boundary of the active region. Such boundary conditions have been derived and analyzed in [START_REF] Antoine | Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schrödinger equation[END_REF][START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF][START_REF] Ben Abdallah | A mathematical model for the transient evolution of a resonant tunneling diode[END_REF][START_REF] Frensley | Boundary conditions for open quantum systems driven far from equilibrium[END_REF] and, at the discrete level, in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF][START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation on circular domains[END_REF][START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF][START_REF] Lent | The quantum transmitting boundary method[END_REF][START_REF] Schulte | Discrete transparent boundary conditions for the Schrödinger equation -a compact higher order scheme[END_REF]. More recently, the discrete transparent boundary conditions approach for the simulation of quantum devices was compared to Perfectly Matched Layers in [START_REF] Mennemann | Perfectly Matched Layers versus discrete transparent boundary conditions in quantum device simulations[END_REF]. Other theoretical studies were devoted to the Schrödinger-Poisson system, dealing with the wellposedness of the nonlinear problem [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF][START_REF] Ben Abdallah | A mathematical model for the transient evolution of a resonant tunneling diode[END_REF][START_REF] Nier | The dynamics of some quantum open systems with short-range nonlinearities[END_REF] and the derivation of simplified models for the RTD [START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells I[END_REF][START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells II[END_REF][START_REF] Jona-Lasinio | On Schrödinger Equations with Concentrated Nonlinearities[END_REF][START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF]. Such simplified models were used, in the stationary case, in [START_REF] Bonnaillie-Noël | Computing the steady states for an asymptotic model of quantum transport in resonant heterostructures[END_REF] to perform fast simulations, which agreement with the numerical solution of the complete Schrödinger-Poisson system was verified in [START_REF] Bonnaillie-Noël | Simulation of resonant tunneling heterostructures: numerical comparison of a complete Schrödinger-Poisson system and a reduced nonlinear model[END_REF]. The simulation of the RTD is a challenging problem due to the presence of quantum resonances which create stiff variations of the wave function with respect to the frequency variable and increase dramatically the numerical complexity. In this context, the reduction of the simulation time provided by the simplified models described above is of particular interest. However, fast simulations with a higher accuracy were obtained in the stationary case in [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] by working on the complete Schrödinger-Poisson system. In addition to a new WKB finite element method (see also [START_REF] Ben Abdallah | An accelerated algorithm for 2D simulations of the quantum ballistic transport in nanoscale MOSFETs[END_REF][START_REF] Negulescu | Numerical analysis of a multiscale finite element scheme for the resolution of the stationary Schrödinger equation[END_REF]), the reduction of the simulation time in this approach is provided by an adapted treatment of the transmission peaks using resonant modes. In the time-dependent case, the supplementary numerical cost imposed by resonances is even more important and its reduction even more challenging. Indeed, the numerical complexity of the transparent boundary conditions of each time-dependent Schrödinger equation is of order one with respect to the number of time iterations, and the last is necessarily big since the period of the high energy free plane waves is much smaller than the scaling time of the system. Although advances (based on the adiabatic theorem) has been made in [START_REF] Faraj | Adiabatic evolution of 1D shape resonances: an artificial interface conditions approach[END_REF][START_REF] Faraj | An explicit model for the adiabatic evolution of quantum observables driven by 1D shape resonances[END_REF] in order to find simplified models for the time-dependent simulations, an algorithm which extends the algorithm of [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] is still missing for the resolution of the transient Schrödinger-Poisson system in the presence of resonances. This is the aim of the present work. In our RTD model, the classical observables, like electronic and current densities, are given by integrals on the frequency variable involving the wave function. To compute correctly these integrals, a uniform frequency mesh has to be very thin to be able to capture the resonance peaks. The high number of Schrödinger equations to be solved, and therefore the high numerical cost, of this standard approach was reduced in the stationary regime in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF] using adaptive refinement. The last method does however not extend to the transient case since the resonant frequencies do move as time varies. Lately, inspired by an idea from [START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF], the algorithm in [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] was proposed for the discretization of the stationary system. The method consists in decomposing the wave function in a resonant part living inside the double barrier and a non resonant part which is mostly localized outside. The non resonant part is smooth with respect to the frequency variable whereas the resonant one has sharp peaks at resonant frequencies. The latter is computed by a projection method after a precomputation of resonant states. Due to the sharpness of the resonance peaks, this process can not succeed without an accurate value of resonance energies and widths. Therefore, we first present an improvement of the method used in [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] to compute these values. The problem of finding quantum resonances has interested many authors, e.g. [START_REF] Alves | Numerical simulations on resonance poles -the trapping case of two plane cracks[END_REF][START_REF] Bindel | Symmetry of Bound and Antibound States in the Semiclassical Limit[END_REF][START_REF] Kim | Cartesian PML approximation to resonances in open systems in IR 2[END_REF][START_REF] Lenoir | Variational formulations for the determination of resonant states in scattering problems[END_REF][START_REF] Riss | Calculation of resonance energies and widths using the complex absorbing potential method[END_REF]]. In the one dimensional case, it can be written as a differential equation with homogeneous transparent boundary conditions selecting outgoing functions, see [START_REF] Faraj | Méthodes asymptotiques et numériques pour le transport quantique résonnant[END_REF][START_REF] Patel | Développement de modèles macroscopiques pour des systèmes quantiques nonlinéaires hors-équilibre[END_REF]. Starting from this approach, we write the problem of finding resonances at the discrete level. Since the spectral parameter is involved nonlinearly in the boundary conditions, we obtain a nonlinear eigenvalue problem for which we derive a Newton-like method which extends to holomorphic matrices the method proposed in [START_REF] Peters | Inverse Iteration, Ill-Conditioned Equations and Newton's Method[END_REF]. Additional references concerning the resolution of nonlinear eigenvalue problems can be found in [START_REF] Guillaume | Nonlinear Eigenproblems[END_REF][START_REF] Szyld | Local convergence of Newton-like methods for degenerate eigenvalues of nonlinear eigenproblems[END_REF]. In the present article, the Fermi-Dirac distribution involved in the densities operates like a frequency cutoff and only the first mode is important in the resonant part of the wave function. The multi-mode approximation follows the same line and finds applications in the simulation of the stationary two-dimensional RTD, see [START_REF] Ben Abdallah | A Generalized Stationary Algorithm for Resonant Tunneling: Multi-Mode Approximation and High Dimension[END_REF].

The key point in the algorithm of [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] for the reduction of the simulation time is to use two scales for the frequency mesh. For the non-resonant part, the frequency mesh can be chosen coarse enough thus reducing the number of Schrödinger equations to be solved. For the resonant part, a single spatial problem (the computation of the resonant mode) is solved for all the frequency points. Therefore a refined frequency mesh can be used for the projection on the resonant mode, capturing this way the resonant peak without increasing significantly the numerical cost. Due to the scale difference between the two parts of the wave function, the projection method requires an interpolation of the non resonant part. In the present work, we show that the time-dependent algorithm is accurate only if the interpolation takes into account the time oscillations coming from the Schrödinger equation. In particular, we verify that the constant interpolation, which is used in the stationary case, is no more adapted in the time-dependent case. For stability reasons, we solve the time-dependent Schrödinger equation with a semi-implicit scheme. In our nonlinear framework, such a scheme requires an adapted extrapolation of the potential at half-time step, see [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF], leading to an Adams-Bashforth-Crank-Nicolson finite difference method. We show that such an extrapolation allows bigger time steps than the extrapolation introduced in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF] for our Schrödinger-Poisson system and also considered in [START_REF] Besse | A Relaxation Scheme for the Nonlinear Schrödinger Equation[END_REF] for a different kind of nonlinearity. In the stationary case, the nonlinearity of the Schrödinger-Poisson system is dealt with by using the self-consistent scheme in [START_REF] Gummel | A self-consistent iterative scheme for one-dimensional steady state transistor calculations[END_REF], see also [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF]. More numerical methods for the time-dependent nonlinear Schrödinger equation can be found in [START_REF] Bao | Numerical Study of Time-Splitting Spectral Discretizations of Nonlinear Schrödinger Equations in the Semiclassical Regimes[END_REF][START_REF] Jin | The Gaussian Beam Methods for Schrödinger-Poisson Equations[END_REF]. Together with the Adams-Bashforth-Crank-Nicolson scheme, we use the discrete boundary conditions in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF][START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF] for the time-dependent Schrödinger equation since they have the advantage to be verified exactly by the whole space discrete solution. We note that these boundary conditions can be used to write several kind of finite difference schemes, including splitting methods, to solve time-dependent nonlinear Schrödinger equations, see [START_REF] Zisowsky | Discrete artificial boundary conditions for nonlinear Schrödinger equations[END_REF]. Applying additional results from [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF][START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF] to our RTD model, we verify numerically that, compared to the results in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF], the transparent boundary conditions for the Schrödinger equation can be fitted to improve the accuracy or simplified to reduce the computational time. In order to be able to compare the time-dependent solution to the stationary one in a meaningful manner, the discretization of the space variable must be the same for the two regimes. Therefore, we will employ a standard finite difference method along with the stationary boundary conditions in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF] to solve the stationary Schrödinger equation and concentrate on the frequency variable to handle the resonance peaks. Similarly, in our algorithm, the resonant part is computed using finite difference in order to be consistent with the non resonant one. For that purpose, we write the discrete boundary conditions for the resonant mode as an extension to complex valued energies of the homogeneous stationary boundary conditions. These boundary conditions are verified exactly by the resonant mode thanks to the following strategy, which is similar to the strategy in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF][START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF]: first we consider the discretization of the differential equation on the whole space, and then we derive the boundary conditions for the difference scheme directly on a purely discrete level. The paper is organized as follows. In section 2, are presented the stationary and time-dependent Schrödinger-Poisson systems. The stationary case is dealt with in section 3 where are given the reference method, the adaptive refinement method and two methods based on the projection on the resonant mode. Our method to compute resonances and resonant modes is also introduced. In section 4, the time-dependent reference method is recalled and our improved transient algorithm is presented. The numerical results are given in section 5. The different methods are compared in the stationary and in the time-dependent case. In the time-dependent case, the comparison shows the improvement provided by our new algorithm and the existence of two resonant peaks is verified. In the appendix are given the finite difference methods with discrete transparent boundary conditions to solve the stationary and time-dependent Schrödinger equation. For the time-dependent scheme, a stability result is proven in the situation of a continuous influx in the device. The discrete problem for the computation of resonances is also obtained and the properties of its solution are presented.

The model

It consists in an infinite number of Schrödinger equations coupled to the Poisson equation. The Schrödinger equation involves the time-dependent Hamiltonian

H(t) = - 2 2m ∂ 2 x + U (t, x) + V (t, x),
where is the reduced Planck constant, m is the effective mass of the electron and x ∈ IR is the position variable. The domain occupied by the device is the interval [0, L], L > 0. The selfconsistent potential V (t, x) is due to Coulomb interaction and depends nonlinearly on the wave functions. The external potential U (t, x) describes the double barrier and the applied bias, see Figure 1, and is given as if a data of the problem. Given the points

0 < a 1 < a 2 < a 3 < b 3 < b 2 < b 1 < L ,
the external potential writes

U (t, x) = v 0 1 [a2,a3[∪]b3,b2] (x) -B(t) x -a 1 b 1 -a 1 1 [a1,b1[ (x) + 1 [b1,+∞[ (x) ,
where v 0 ≥ 0 and B(t) ≥ 0 are scalars representing respectively the height of the barrier and the amplitude of the applied bias. The points a 1 and b 1 are the extremities of the diode. The interval [a 3 , b 3 ] is the quantum well. The problem is the following infinite system of Schrödinger equations

X X X X X X X 6 v 0 ? -B(t) 0 L a 1 b 1 a 2 b 2 a 3 b 3 Figure 1: External potential U (t, x). i ∂ t Ψ k (t, x) = - 2 2m ∂ 2 x Ψ k (t, x) + (U (t, x) + V (t, x)) Ψ k (t, x) , x ∈ IR Ψ k (0, x) = Φ k (x) , (2.1)
where the initial condition Φ k (x) is a generalized eigenfunction of the initial Hamiltonian H(0) corresponding to the frequency k ∈ IR. The functions Φ k (x) are defined more in details at the end of this section. The self-consistent potential satisfies the Poisson equation:

-∂ 2 x V (t, x) = q 2 ε (n(t, x) -n D (x)) , x ∈ (0, L) V (t, 0) = V (t, L) = 0 (2.2)
and is extended by 0 outside (0, L). In (2.2), q is the elementary charge of the electron, ε is the dielectric constant, n D (x) is the doping equal to

n D (x) = n 1 D 1 [0,a1[ (x) + 1 ]b1,L] (x) + n 2 D 1 [a1,b1] (x) , with n 1 D > n 2 D ≥ 0.
The electron density n(t, x) is given by:

n(t, x) = I R g(k)|Ψ k (t, x)| 2 dk . (2.3)
In our model, the injection profile g(k) is the one dimensional Fermi-Dirac distribution:

g(k) = mk B T 2π 2 2 ln 1 + exp E F - 2 k 2 2m k B T , (2.4)
where k B is the Boltzmann constant, T is the temperature of the semiconductor and E F is the Fermi level.

In the stationary regime, the external potential is the potential U I (x) = U (0, x) corresponding to the initial bias B I = B(0). The time-independent problem is the following infinite system of equations

- 2 2m ∂ 2 x Φ k (x) + (U I (x) + V I (x)) Φ k (x) = E k Φ k (x) , x ∈ IR , (2.5)
with scattering conditions

Φ k (x) = e ikx + r(k)e -ikx , x < 0 Φ k (x) = t(k)e i √ k 2 +2mBI / 2 x , x > L for k ≥ 0 (2.6) and Φ k (x) = t(k)e -i √ k 2 -2mBI / 2 x , x < 0 Φ k (x) = e ikx + r(k)e -ikx , x > L for k < 0 , (2.7)
where

E k = 2 k 2 2m , k ≥ 0 2 k 2 2m -B I , k < 0 . (2.8)
As in the time-dependent case, the self-consistent potential V I (x) satisfies the Poisson equation:

-∂ 2 x V I (x) = q 2 ε (n I (x) -n D (x)) , x ∈ (0, L) V I (0) = V I (L) = 0 (2.9)
and is extended by 0 outside (0, L). The electron density n I (x) is given by:

n I (x) = I R g(k)|Φ k (x)| 2 dk .
(2.10)

3 The stationary algorithm revisited

Recalling the standard algorithm

To solve the stationary nonlinear Schrödinger-Poisson system (2.5)(2.9)(2.10), we will use the same algorithm as in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF]. It is based on a Gummel iteration, see [START_REF] Gummel | A self-consistent iterative scheme for one-dimensional steady state transistor calculations[END_REF], which consists in the computation of a sequence V l I l∈I N where the potential V l+1 I at step l + 1 is deduced from the potential V l I at step l by solving the following nonlinear equation:

-∂ 2 x V l+1 I (x) = q 2 ε n l I (x) exp V l I (x) -V l+1 I (x) /V ref -n D (x) , x ∈ (0, L) V l+1 I (0) = V l+1 I (L) = 0 (3.1)
for a fixed reference potential V ref . In (3.1), n l I is the density (2.10) where the wave functions are solution to (2.5) with V I = V l I . The repartition function g being exponentially decreasing at infinity, the integral in (2.10) can be restricted, in computations, to a domain [-κ, κ] where κ is chosen to be big enough. We consider a discretization x 0 = 0, x 1 , ..., x j , ..., x J = L of the interval [0, L] with uniform grid spacing ∆x and a discretization k 0 = -κ, k 1 , ..., k p , ..., k P = κ of the interval [-κ, κ] with non-necessarily uniform grid spacing ∆k p = k p+1 -k p . In order to verify, the condition (A.4) for all k ∈ [-κ, κ], the grid spacing for the space variable is chosen such that: ∆x < 1

κ 2 + 2m 2 B I . (3.2)
Then, for a given initial potential V 0 I , the algorithm writes:

Set

∆ = V l+1 I -V l I L 2 (0,L) V l+1 I L 2 (0,L) . (3.3)
Set l = l + 1. End Do.

Equation (3.1) is nonlinear and is solved with a Newton method where the Laplacian is discretized using finite difference at the points (x j ) 0≤j≤J . The approximation n l I,j of the density n l I (x j ) is computed by using the trapezoidal rule:

n l I,j = P -1 p=0 g(k p ) Φ l p,j 2 + g(k p+1 ) Φ l p+1,j 2 ∆k p 2 , j = 0, ..., J (3.4)
where Φ l p,j denotes the approximation of the wave function Φ l p (x j ). In our numerical experiments, the norm . L 2 (0,L) appearing in equation (3.3) is replaced by the discrete l 2 -norm

|v| 2 =   J j=0 |v j | 2   1 2
for v ∈ C J+1 , and V l I is replaced by the vector V l I,j 0≤j≤J where V l I,j is the approximation of the self-consistent potential V l I (x j ). Because of the peaked form of the transmission near resonances, the method used to compute the density n l I from the potential V l I (steps S1 and S2) is of major importance, and the Gummel algorithm may fail to converge if a non accurate method is used. In the case of the finite difference discretization, the strategy of [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF] can be described as follows. The step S1 is performed without particular treatment of resonances: the function Φ l p is computed on [0, L] by the resolution of the problem (2.5)-(2.8) with the finite difference scheme (A.3) and the discrete transparent boundary conditions (A.7), (A.10) for k ≥ 0 and (A.11), (A.12) for k < 0 where the potential Q I is replaced by U I + V l I . In that case, the convergence is provided by the choice of the frequency mesh {k p , p = 0, ..., P } for the trapezoidal rule (3.4). The reference method is obtained by taking the grid spacing ∆k p equal to a constant ∆k > 0 small enough. This method will be called Direct Resolution. In order to reduce the high numerical cost, which is due to the big number of Schrödinger equations to be solved, an Adaptive Method is proposed in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF]. For the Adaptive Method, the grid spacing ∆k p is taken equal to the small constant ∆k only when E kp is close to a resonant energy, otherwise ∆k p = ν∆k where ν ≥ 2. The frequency points where the mesh is refined are detected using the logarithmic derivative, with respect to the frequency k, of the transmission coefficient related to the wave function Φ k (see [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF][START_REF] Faraj | Méthodes asymptotiques et numériques pour le transport quantique résonnant[END_REF]). As noticed in the introduction, the Adaptive Method can not be generalized to the time-dependent case. In this context, the One Mode Approximation algorithm presented in section 3.3 is very useful. As proposed in [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF], this algorithm consists in decomposing the wave function in a non resonant part and a resonant one proportional to the first resonant mode. In [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF], using a WKB interpolation to compute each part of the wave function, an adapted treatment of the step S1 is realized.

In the present work, we do not use the WKB interpolation, however the accuracy required at step S1 is reached by an improvement of the computation of the first resonance. In particular, a precise computation of its imaginary part is essential. A simple reconstitution of the wave functions not allowing to modify the frequency mesh, the gain in computational time is reached by adapting the step S2. In the present work, two strategies are proposed for the step S2. The first method is inspired from [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] and consists in using different frequency scales for the different parts of the wave function. The frequency mesh is taken thin only for the coefficient of proportionality to the resonant mode and the number of Schrödinger equations to be solved is reduced. This method will be denoted One Mode Approximation and will be generalized to the time-dependent case. A variant of this method is obtained by an explicit integration of the coefficient of proportionality to the resonant mode instead of the trapezoidal rule. In that case all the resonant peak is integrated and the frequency mesh can be chosen coarse enough. This method will be denoted One Mode Approximation Integ. The algorithms Direct Resolution, Adaptive Method, One Mode Approximation and One Mode Approximation Integ are compared in section 5.2.

Accurate computation of resonances

In the founder work of [START_REF] Aguilar | A Class of Analytic Perturbations for One-body Schrödinger Hamiltonians[END_REF], see also [START_REF] Hislop | Introduction to Spectral Theory with Applications to Schrödinger Operators[END_REF], resonances of a self-adjoint operator are defined using analytic transformations, and it is a common fact that it corresponds to an eigenvalue in a modified L 2 space [START_REF] Fujiié | Matrice de scattering et résonances associées à une orbite hétérocline[END_REF]. Using the second approach, the authors in [START_REF] Patel | Développement de modèles macroscopiques pour des systèmes quantiques nonlinéaires hors-équilibre[END_REF] and [START_REF] Faraj | Méthodes asymptotiques et numériques pour le transport quantique résonnant[END_REF] define resonances for Hamiltonians of the form

- 2 2m ∂ 2 x + Q(x) , (3.5)
where the potential Q(x) verifies

Q(x) = 0 , x ≤ 0 and Q(x) = Q L , x ≥ L where Q L ≤ 0 . (3.6)
Starting from this definition, the resonances of the Hamiltonian (3.5) can be characterized as the numbers

z in C \ {iIR ∪ iIR + Q L } such that the problem [- 2 2m ∂ 2 x + Q(x)]u(x) = zu(x) , x ∈ IR ||u|| L 2 (0,L) = 1 (3.
M (z)u = 0 u H u = 1 , (3.8)
where M (z) is the nonlinear matrix valued function of z given by (B.10) and u H denotes the complex conjugate transpose of the vector u. A Newton-like method is introduced in [START_REF] Peters | Inverse Iteration, Ill-Conditioned Equations and Newton's Method[END_REF] to solve (3.8) when M (z) is a polynomial. Following the same line, we derive a method for the general case where M (z) is holomorphic. Note first that the term u H u appearing in (3.8) is not differentiable with respect to u, therefore a Newton method must be modified to be applied here. For a given iterate (u n , z n ) verifying (u n ) H u n = 1, we are looking for a direction (δu n , δz n ) such that (u n + δu n , z n + δz n ) is solution to the problem (3.8). Using (u n ) H u n = 1 the system:

M (z n + δz n )(u n + δu n ) = 0 (u n + δu n ) H (u n + δu n ) = 1
gives at order 2:

M (z n )δu n + δz n M ′ (z n )u n = -M (z n )u n (u n ) H δu n + (δu n ) H u n = 0 . (3.9)
As noticed earlier, the second equation in (3.9) does not have the suitable form. To deal with this problem, we remark that it is enough to impose (u n ) H δu n = 0 in order to verify this equation. This leads to the following linear system:

M (z n ) M ′ (z n )u n (u n ) H 0 δu n δz n = -r n 0 , where r n = M (z n )u n . (3.10)
Its resolution corresponds to an iteration of our Newton-like method to compute the resonant mode and the resonance. The assumption (u n ) H u n = 1, that we made to obtain the system (3.10), is verified at order 2 as long as (u 0 ) H u 0 = 1.

The One Mode Approximation

We start with the description of the step S1 of the computation of a wave function Φ k (x) solution to the stationary Schrödinger equation (2.5) for a given frequency k. To simplify the notations, the exponent l appearing in Algorithm 3.1 will be omitted. Following the works [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] and [START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF], the One Mode Approximation consists in the decomposition:

Φ k (x) = Φ nr k (x) + Φ r k (x) , (3.11)
where the non resonant part Φ nr k (x) solves the stationary Schrödinger equation:

[- 2 2m ∂ 2 x + U I,f ill (x) + V I (x)]Φ nr k (x) = E k Φ nr k (x) , x ∈ IR , (3.12)
with filled potential

U I,f ill (x) = U I (x) + v 0 1 [a3,b3] (x)
. Equation (3.12) comes with the scattering conditions (2.6), (2.7) and the relation (2.8). It is solved on [0, L] using the finite difference scheme (A.3) and the discrete transparent boundary conditions (A.7), (A.10) for k ≥ 0 and (A.11), (A.12) for k < 0 where the potential Q I is replaced by U I,f ill + V I . The statistic g(k) being rapidly decreasing at infinity, the resonant part Φ r k (x) is sought, on [0, L], proportional to the resonant mode u I (x) of minimal resonant energy Re (z

I ) solution to [- 2 2m ∂ 2 x + U I (x) + V I (x)]u I (x) = z I u I (x) and L 0 |u I (x)| 2 dx = 1 . (3.13)
In other words, we look for Φ r k (x) of the form:

Φ r k (x) = θ k u I (x), x ∈ [0, L] . (3.14)
Inserting (3.11) and (3.14) in the stationary Schrödinger equation (2.5), and using (3.12) and (3.13), we get the following explicit value of the proportionality coefficient:

θ k = 1 z I -E k v 0 b3 a3 Φ nr k u I dx . (3.15)
The resonant mode and the resonance are computed using the method presented in section 3.2 with the potential Q = U I + V I . This method is initialized at the fundamental mode and fundamental energy of the Hamiltonian:

[- 2 2m ∂ 2 x + U I (x) + V I (x)] (3.16)
equipped with homogeneous Dirichlet boundary conditions at a 2 and b 2 . It is shown in [START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells I[END_REF][START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells II[END_REF] that the real part of the resonances are close to the eigenvalues of the Dirichlet Hamiltonian (3.16).

The imaginary part of resonances being small, such an initialization insures the convergence of the algorithm to the resonance with smaller energy. This achieves the step S1.

For the One Mode Approximation, the step S2 consists in using the wave function decomposition (3.11) and (3.14) in order to compute accurately the density with a low numerical cost. First, using an argument of localization of support, see [START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF], the cross term in the development of

|(Φ nr k + θ k u I ) (x)|
2 can be neglected, which gives the approximation:

n I (x) = I R g(k) |Φ nr k (x)| 2 dk + I R g(k) |θ k | 2 dk |u I (x)| 2 .
(3.17)

The non resonant part of the wave function being regular with respect to k, the first integral above can be computed using a trapezoidal rule with coarse frequency mesh, reducing this way the number of Schrödinger equations to be solved. As remarked above, the resonance z I has a small imaginary part. Therefore, it appears from formula (3.15), that the coefficient θ k has a sharp peak when the wave function energy E k is close to the resonant energy Re (z I ). The key-point in the step S2 is then the approximation of the integral I R g(k) |θ k | 2 dk. For the sake of simplicity, the non resonant part of the wave function will be called non resonant wave function as well. The same remark holds for the time-dependent One Mode Approximation presented in section 4.2.

First method: One Mode Approximation

Consider a refined discretization k 0 = -κ, k 1 , ..., k p , ..., k P = κ of the interval [-κ, κ] with uniform grid spacing ∆k and suppose that there exist two integers P ′ and ν such that P P ′ = ν. Consider the large scale mesh kp ′ = k νp ′ for p ′ = 0, ..., P ′ with uniform grid spacing ∆ k = ν∆k. Then, if we note Φ nr p ′ ,j the approximation of the non resonant wave function Φ nr k νp ′ (x j ), u I,j the approximation of the resonant mode u I (x j ) and θ p the approximation of the coefficient θ kp , the formula for the density is

n I,j = P ′ -1 p ′ =0 g( kp ′ ) Φ nr p ′ ,j 2 + g( kp ′ +1 ) Φ nr p ′ +1,j 2 ∆ k 2 + P -1 p=0 g(k p ) |θ p | 2 + g(k p+1 ) |θ p+1 | 2 ∆k 2 |u I,j | 2 .
The crucial point for the reduction of the simulation time is that only P ′ Schrödinger equations have to be solved, instead of P equations for the Direct Resolution. However, this reduction implies that we only have access to the functions Φ nr p ′ ,j for p ′ = 0, ..., P ′ and the computation of the coefficient θ p , for p = νp ′ , requires an interpolation of the non resonant wave function in the r.h.s. of equation (3.15). It is provided by the piecewise constant interpolation formula below: for p = 0, ..., P

θ p = 1 z I -E p v 0 j∈w Φ nr p ′ ,j u I,j ∆x , (3.18)
where

w = {j | a 3 ≤ x j < b 3 }, E p = E kp and 0 ≤ p ′ ≤ P ′ is the integer such that p = νp ′ + r for some 0 ≤ r ≤ ν -1.
To simplify the presentation, we use in (3.18), and in what follows, the same notation for the resonance z I and its approximation. Since the evaluation of the coefficient θ p in formula (3.18) is numerically cheap, the number P can be chosen big enough to catch the resonance peak without a significant increase of the computational cost.

Second method: One Mode Approximation Integ

We only considered here the large scale frequency mesh k0 , ..., kP ′ introduced in section 3.3.1. As noticed above, the integral of the non resonant part can be computed with a simple trapezoidal rule. Then, if we set J p := kp+1 kp g(k) |θ k | 2 dk, we get the following formula for the density:

n I,j = P ′ -1 p=0 g( kp ) Φ nr p,j 2 + g( kp+1 ) Φ nr p+1,j 2 ∆ k 2 + P ′ -1 p=0 J p |u I,j | 2 . (3.19)
In that case, the required precision for the density is reached thanks to an adapted treatment of the integral in J p . It follows from (3.15) that:

J p = kp+1 kp S k k |E k -z I | 2 dk ,
where

S k = g(k) k v 2 0 b3 a3 Φ nr k u I dx 2 .
Using numerical observations, it is shown in [START_REF] Faraj | Méthodes asymptotiques et numériques pour le transport quantique résonnant[END_REF] that Φ nr k L ∞ (0,L) = O(k) when k tends to 0. This limit is also verified theoretically in the same work, under the condition that k > 0 and for particular potentials. Therefore, S k is a smooth function of the variable k in a small neighborhood of 0. Since the wave function Φ nr k has no resonance peak, S k is also smooth outside such a neighborhood. As a consequence, S k can be approximated by a constant on [ kp , kp+1 ] which gives:

J p = S p kp+1 kp k |E k -z I | 2 dk ,
where

S p := g kp kp v 2 0 j∈w Φ nr p,j u I,j ∆x 2
and w is defined in section 3.3.1. The resonance z I can be written

z I = E I -iΓ I , (3.20)
where

E I ∈ IR and Γ I > 0. Under the assumption 0 / ∈ [ kp , kp+1 ], it holds E k = 2 k 2 2m -B p on [ kp , kp+1 ],
where B p is a constant equal to 0 or B I . Then, it follows:

J p = S p kp+1 kp k 2 k 2 2m -B p -E I 2 + Γ 2 I dk .
Using the change of variable E = 2 k 2 2m -B p -E I in the above integral, we obtain the following approximation of J p :

J p = mS p 2 Γ I arctan Êp+1 -E I Γ I -arctan Êp -E I Γ I for p = 0, ..., P ′ , (3.21)
where Êp := E kp . We will use (3.21) regardless of the verification of the assumption 0 / ∈ [ kp , kp+1 ].

It provides an accurate treatment of the resonance peak even if ∆ k is large. The approximation of the density follows by replacing (3.21) in (3.19).

Remark 3.2. The computational time of the One Mode Approximation algorithm can be reduced when the method to compute the resonance is initialized at the fundamental mode and fundamental energy (u D , E D ), of the Dirichlet Hamiltonian (3.16), only for the first iteration l = 0 of the Gummel algorithm 3.1. Indeed, for l ≥ 1, the computation of (u D , E D ) is avoided by initializing the method at the resonant mode and resonance obtained at the previous iteration l -1.

Remark 3.3. The size of the grid spacing ∆k for the refined frequency mesh is chosen numerically.

It can be estimated analytically by imposing that, around the resonance energy E I , the difference E p+1 -E p is smaller than half the resonance peak width Γ I .

The time-dependent algorithm

This section deals with the time evolution of the system under a time-dependent applied bias B(t) as in Figure 1. We will consider only the case of the step like bias of the form:

B(t) = B I , t = 0 B ∞ , t > 0 , ( 4

The algorithm

As in the stationary regime, we first present the algorithm proposed in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF] to solve the timedependent nonlinear Schrödinger-Poisson system (2.1)(2.2)(2.3). Then, we give the details of our time-dependent One Mode Approximation algorithm. The space and frequency meshes are defined as in section 3.1. In order to verify the condition (C. [START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells I[END_REF]) for all k ∈ [-κ, κ], the time step ∆t is chosen such that:

∆t < π 2 2 κ 2 2m + sup t≥0 B(t) . (4.2)
Then, for a given exterior potential U (t, x), the algorithm corresponds to the computation of the sequence V l of approximations of the self-consistent potential at time t l = l∆t. It provides a sequence of densities n l at time t l . The initial potential V 0 , resp. wave function Ψ 0 p , is given by V I , and resp. Φ kp , solution to the stationary Schrödinger-Poisson system (2.5)(2.9)(2.10) with exterior potential equal to U I . The Schrödinger equation (2.1) is solved using a Crank-Nicolson scheme and, in our nonlinear framework, a choice of the intermediary potential V l+ 1 2 which preserves the unconditional stability of the semi-implicit scheme is required. In particular, it will be shown in section 5.3.1 that the extrapolation:

V l+ 1 2 = 3 2 V l - 1 2 V l-1 (4.3)
allows bigger time steps than the following extrapolation proposed in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF]:

V l+ 1 2 = 2V l -V l-1 2 . (4.4)
The extrapolation (4.3) is proposed in [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF] and corresponds to a modified Adams-Bashforth-Crank-Nicolson method. We can now give the time-dependent algorithm: Do For l = 0 to l tot -1:

Computation of the intermediary potential

Computation of the potential V l+ 1 2 from the potentials V l and V l-1 using (4.3) with the convention

V -1 = V 0 .
Computation of the density S1. For p = 0, ..., P : computation of the wave function Ψ l+1 p , solution to (2.1) at time t l+1 and frequency k p , using the potential V l+ 1 2 and the wave function Ψ l p . S2. Numerical integration: computation of the density n l+1 from the wave functions Ψ l+1 p 0≤p≤P by using equation (2.3).

Coupling, Poisson equation

Computation of the potential V l+1 from the density n l+1 by solving the Poisson equation (2.2):

-∂ 2 x V l+1 (x) = q 2 ε n l+1 (x) -n D (x) , x ∈ (0, L) V l+1 (0) = V l+1 (L) = 0 . (4.5)
End Do.

The approximation n l+1 j of the density n l+1 (x j ) is computed by using the trapezoidal rule:

n l+1 j = P -1 p=0 g(k p ) Ψ l+1 p,j 2 + g(k p+1 ) Ψ l+1 p+1,j 2 ∆k p 2 , j = 0, ..., J (4.6)
where Ψ l+1 p,j denotes the approximation of the wave function Ψ l+1 p (x j ). The linear equation (4.5) is discretized with finite difference and solved by a simple matrix inversion. As in the stationary case, the steps S1 and S2 to compute the density will be crucial for the treatment of the resonance peaks. In the algorithm proposed in [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF], that we will call Direct Resolution, the step S1 is performed as follows: the function Ψ l+1 p is computed on [0, L] using Ψ l p and V l+ 1 2 by the resolution of (2.1) with the Crank-Nicolson scheme (C.2), where the potential Q is equal to U + V . Since the support of the initial data Φ kp is not included in (0, L), the discrete transparent boundary conditions for (C.2) are the non-homogeneous ones (C.9)(C.12) or their simplified version (C.22)(C.23). In the Direct Resolution, no particular treatment of the resonances is performed in the step S1 and the accuracy of the method is provided in the step S2 by imposing of a small uniform grid spacing ∆k p = ∆k in the trapezoidal rule (4.6). Indeed, if the grid spacing is small only near the initial resonant frequencies, the refined mesh will loose its relevance because of the time evolution of the resonance. Therefore, the stationary Adaptive Method does not extend to the transient regime and, if no additional treatment is performed (a strategy would be to choose a small frequency grid spacing only in the regions where will live the resonance peaks), the number of Schrödinger equations to be solved have to be important in the time-dependent Direct Resolution. Moreover, in addition to a matrix inversion, which is equally required in the stationary case, each timedependent Schrödinger equation requires discrete convolutions in the boundary conditions. For the exact boundary conditions (C.9)(C.12), the number of operations required by these convolutions is proportional to the time iteration number l. For the simplified boundary conditions (C.22)(C.23), the size of the convolutions can be reduced but it can not be too small in order to preserve the accuracy of the method. This increases the numerical cost of each time step and makes the timedependent simulations very long: actually, due to the condition (4.2), the time step must be chosen much smaller than the scaling time of the system and the number l can become very big. In this context, it is important to look for an adapted treatment of the resonance peaks to reduce the number of frequency points. In the following section, we propose a One Mode Approximation method which extends the method proposed in section 3.3 to the time-dependent case.

The time-dependent One Mode Approximation

Let us start with the description of the decomposition of a wave function Ψ k (t, x) solution to the transient Schrödinger equation (2.1) for a given frequency k. We look for Ψ k (t, x) of the form:

Ψ k (t, x) = Ψ nr k (t, x) + Ψ r k (t, x) , (4.7)
where the non resonant part Ψ nr k (t, x) solves the transient Schrödinger equation:

i ∂ t Ψ nr k (t, x) = [- 2 2m ∂ 2 x + U f ill (t, x) + V (t, x)]Ψ nr k (t, x) , x ∈ IR Ψ nr k (0, x) = Φ nr k (x) , (4.8)
with filled potential

U f ill (t, x) = U (t, x) + v 0 1 [a3,b3] (x) and Φ nr k (x
) is solution to (3.12). Using the same argument as in the stationary regime, we suppose that Ψ r k (t, x) is proportional to the resonant mode u(t, x) corresponding to the first resonance z(t) of the Hamiltonian H(t) at time t:

in other words Ψ r k (t, x) = λ k (t)u(t, x) , x ∈ [0, L] , (4.9)
where [-

2 2m ∂ 2 x + U (t, x) + V (t, x)]u(t, x) = z(t)u(t,
Φ k (x) = Φ nr k (x) + θ k u I (x) , x ∈ [0, L] . ( 4
[i λ ′ k (t) -z(t)λ k (t)]u(t, x) + i λ k (t)∂ t u(t, x) = -v 0 1 [a3,b3] (x)Ψ nr k (t, x) .
Multiplying the previous equation by u(t, x) and integrating on (0, L), it follows:

λ ′ k (t) + [ i z(t) + L 0 ∂ t u(t, x)u(t, x)dx]λ k (t) = i v 0 b3 a3 Ψ nr k (t, x)u(t, x)dx λ k (0) = θ k . (4.12)
Remark 4.2. In the linear case, obtained by taking V (t, x) = 0 in (2.1), it holds for t > 0:

u(t, x) = u ∞ (x) and z(t) = E ∞ -iΓ ∞
, where E ∞ ∈ IR and Γ ∞ > 0. Then, under the assumption

Ψ nr k (t, x) = Ψnr k (x)e -i ε ∞ k t , where ε ∞ k = 2 k 2 2m , k ≥ 0 2 k 2 2m -B ∞ , k < 0 , (4.13)
equation (4.12) is an ODE which homogeneous solution oscillates at the angular frequency E∞ and which source term oscillates at the angular frequency

ε ∞ k . Therefore, λ k (t) has a peak at the frequencies k such that ε ∞ k = E ∞ .
Due to the initial condition in (4.12), one expects a coexistence of this peak with a second peak at the frequencies k such that E k = E I with decay rate Γ∞ , where E k and E I are given by (2.8) and (3.20). This is verified numerically in section 5.3.4.

Step S1

The aim here is to compute the elements of the decomposition (4.7)(4.9) at the time t = t l+1 . At this point of the algorithm, we do not have the value of the potential V l+1 and the corresponding resonant mode u l+1 can not be computed. Therefore, we will make the following approximation for the wave function at time t l+1 and frequency k: The resonant mode u l+ 1 2 and the resonance z l+ 1 2 are computed using the method presented in section 3.2 where the potential Q is equal to U l+ 1 2 + V l+ 1 2 and U l+ 1 2 (x) = U (t l+ 1 2 , x). The initial guess for (u l+ 1 2 , z l+ 1 2 ) are the first eigenfunction and eigenenergy of the Hamiltonian [-

Ψ l+1 k (x) = Ψ nr,l+1 k (x) + λ l+1 k u l+ 1 2 (x) . ( 4
2 2m ∂ 2 x + U l+ 1 2 (x) + V l+ 1 2 (x)] (4.15)
with homogeneous Dirichlet boundary conditions at a 2 and b 2 .

To achieve the step S1, we have to compute the coefficient λ l+1 k at time t l+1 and frequency k. It is obtained from λ l k , Ψ nr,l k , Ψ nr,l+1 k , z l+ 1 2 and u l+ 1 2 by the resolution of equation (4.12) with a trapezoidal rule. This leads to the following iteration:

(λ l+1 k -λ l k )/∆t + [ i z l+ 1 2 + L 0 ∂ t u l+ 1 2 u l+ 1 2 dx](λ l k + λ l+1 k )/2 = iv 0 2 b3 a3 Ψ nr,l k + Ψ nr,l+1 k u l+ 1 2 dx , (4.16)
where ∂ t u l+ 1 2 stands for ∂ t u(t l+ 1 2 ). By adequately fixing the resonant mode phase, the quantity

µ l+1/2 := L 0 ∂ t u l+1/2 u l+1/2 dx (4.17)
appearing in (4. [START_REF] Bonnaillie-Noël | Far from equilibrium steady states of 1D-Schrödinger-Poisson systems with quantum wells I[END_REF]) is fitted to zero. Given ũl+1/2 resonant mode of lower energy of the Hamiltonian (4.15) and verifying ũl+1/2 L 2 (0,L) = 1, we look for u l+1/2 of the form u l+1/2 = ũl+1/2 e iϕ l+1/2 where ϕ l+1/2 ∈ IR is such that µ l+1/2 can be neglected. We note first that we have the approximation

µ l+1/2 = 1 2∆t L 0 (u l+1/2 -u l-1/2 )(u l+1/2 + u l-1/2 )dx , which becomes µ l+1/2 = i ∆t Im[ L 0 u l+1/2 u l-1/2 dx] (4.18)
under the condition u l-1/2 L 2 (0,L) = u l+1/2 L 2 (0,L) = 1. Then, defining

ω l+1/2 = L 0 ũl+1/2 u l-1/2 dx , (4.19)
we choose

e iϕ l+1/2 = ω l+1/2 |ω l+1/2 |
and it follows:

L 0 u l+1/2 u l-1/2 dx = e iϕ l+1/2 L 0 ũl+1/2 u l-1/2 dx = ω l+1/2 e iϕ l+1/2 = |ω l+1/2 | ∈ IR .
As a consequence, Im[ L 0 u l+1/2 u l-1/2 dx] = 0 and equation (4.18) shows that u l+1/2 is such that µ l+1/2 is almost equal to 0. In the applications considered in section 5, the space integral in (4. [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF]) is computed using a rectangle method.

Remark 4.3. As it is done in the stationary case, see Remark 3.2, the algorithm can be accelerated by choosing (u l-1 2 , z l-1 2 ), for l ≥ 1, as initial guess of the method to compute (u l+ 1 2 , z l+ 1 2 ). When l = 0, the initial guess is the resonant mode and the resonance given by the stationary solution.

Remark 4.4. Although it is harder to implement than an explicit Euler method, the semi-implicit trapezoidal rule (4.16), used to solve the stiff-like equation (4.12), is required for stability purposes. Indeed, the imaginary part of the resonance z(t) being strictly negative, the ODE (4.12) is similar to the test equation: y ′ = αy where Re (α) < 0, for which the trapezoidal rule solution is bounded independently of the time step whereas the Euler explicit solution is conditionally bounded. In particular, as explained above, the quantity µ l+1/2 defined in (4.17) can been fitted to zero in (4.16), and a straightforward calculation shows that the resulting algorithm is unconditionally stable, in the sens that for all N ≥ 0

λ l k ≤ λ 0 k + T max 0≤n≤N -1 iv 0 2 b3 a3 Ψ nr,n k + Ψ nr,n+1 k u n+ 1 2 dx , 0 ≤ l ≤ N ,
where T = N ∆t.

Step S2

The step S2 is performed using the decomposition (4.14) to find an approximation of the density n l+1 j in (4.6) which provides an adapted treatment of the resonant peaks. Like in section 3.3, we make the approximation:

Ψ l+1 k (x) 2 = Ψ nr,l+1 k (x) 2 + λ l+1 k 2 u l+ 1 2 (x) 2 . (4.20)
The non resonant wave function Ψ nr k is regular with respect to the frequency k and the integral of the first term at the r.h.s. of equation (4.20) can be computed using a trapezoidal rule with coarse frequency mesh. As it is noticed in Remark 4.2, the coefficient λ k has sharp peaks at the resonant frequencies and the frequency mesh is taken thin for the integral of the second term at the r.h.s of equation (4.20). More precisely, we consider the frequency meshes { kp ′ , p ′ = 0, ..., P ′ } and {k p , p = 0, ..., P } introduced in section 3.3.1 and which have different scales adjusted using the ratio ν = P P ′ ∈ N. Then, if we note Ψ nr,l p ′ ,j the approximation of the non resonant wave function Ψ nr,l k νp ′ (x j ), λ l p the approximation of the coefficient λ l kp and u

l+ 1 2 j
the approximation of the resonant mode u l+ 1 2 (x j ), the formula for the density is:

n l+1 j = P ′ -1 p ′ =0 g( kp ′ ) Ψ nr,l+1 p ′ ,j 2 + g( kp ′ +1 ) Ψ nr,l+1 p ′ +1,j 2 ∆ k 2 + P -1 p=0 g(k p ) λ l+1 p 2 + g(k p+1 ) λ l+1 p+1 2 ∆k 2 u l+ 1 2 j 2 . (4.21)
The number of Schrödinger equations to solve is reduced, P ′ equations instead of P equations for the Direct Resolution, reducing this way the numerical cost. However, this reduction implies that we have only access to the functions Ψ nr,l p ′ ,j for p ′ = 0, ..., P ′ and the computation of the coefficient λ l+1 p , for p = νp ′ , requires an interpolation of the non resonant wave function to evaluate the source term in (4.16). If we set:

ε ∞ p = ε ∞ kp
where ε ∞ k is given by (4.13), it follows from Remark 4.2 that the coefficient λ l+1 p has a peak when the energy ε ∞ p is equal to E l+1 (the lower resonant energy corresponding to the potential V l+1 ) only if the approximation of the source term in (4.16) includes the time oscillations of the non resonant wave function. In particular, a polynomial interpolation of the non resonant wave function is not adapted. Solving equation (4.16) and performing an oscillatory interpolation of the non resonant wave function, we obtain the following suitable algorithm which extends the formula (3.18) to the time-dependent case: for p = 0, ..., P

λ l+1 p = 1 1 + i ∆tz l+ 1 2 2 1 -i ∆tz l+ 1 2 2 λ l p + iv 0 ∆t 2 j∈w Ψnr,l p ′ ,j e -iω ∞ p t l + Ψnr,l+1 p ′ ,j e -iω ∞ p t l+1 u l+ 1 2 j ∆x   , (4. 22
)
where w is defined in section 3.3.1 and 0 ≤ p ′ ≤ P ′ is the integer such that p = νp ′ + r for some 0 ≤ r ≤ ν -1. The amplitude-like function Ψnr,l p ′ ,j appearing in equation (4.22) is defined by

Ψnr,l p ′ ,j = Ψ nr,l p ′ ,j e iω ∞ νp ′ t l and ω ∞ p = ε ∞ p . (4.23)
To simplify the presentation, we used in (4.22) the same notation for the resonance z l+ 1 2 and its approximation.

Remark 4.5. The condition (4.2) insures that the time step is small enough with respect to one period of the oscillating term e -iω ∞ p t in (4.22). The time-dependent One Mode Approximation extends the method in section 3.3.1 and it is numerically verified that it is stable in time and allows long time simulations, see section 5.3.3. On the other hand, due to the function e -iω ∞ p t , a time-dependent extension of the method presented in section 3.3.2 is not obvious and gives a non-stable algorithm if no particular treatment is performed. Remark 4.6. In section C.2.2, it is noticed that the accuracy of the algorithm to solve the timedependent Schrödinger equation is improved by replacing the angular frequencies (C.8)(C.11) by the discrete ones (C.13)(C.14) in the boundary conditions (C.9)(C.12). Similarly, the accuracy of the time-dependent One Mode Approximation is improved when the angular frequency (4.23) is replaced by the discrete one below:

ω ∞ p = 2 ∆t arctan ∆tε ∞ p 2 . (4.24)

Results

The physical parameters used for the numerical computations are gathered in the following array: For all the tests, the two barriers have the same size which is equal to the size of the well. The external potential is determined by Figure 1 and the data below: In all the simulations, we fixed κ = 2m 2 (E F + 7k B T ) and the number of space points is equal to J = 300 which is such that the condition (3.2) is verified.

L (nm) a 1 (nm) a 2 (nm) a 3 (nm) b 3 (nm) b 2 (nm) b 1 (nm) v 0 (

Computation of resonances

In this section, we give the numerical values of the resonance of lower energy obtained using the algorithm presented in section 3.2 for different biases B I . The potential Q is equal to U I + V I where U I is the external potential corresponding to B I and V I the corresponding solution to the Schrödinger-Poisson system (2.5)(2.9)(2.10) given by the Direct Resolution, with P = 3750, as described in section 5.2. We obtain the following results: where the number of iterations before convergence N cv denotes the iteration n of the algorithm such that |M (z n )u n | 2 < 10 -15 , E D denotes the fundamental energy of the Dirichlet Hamiltonian (3.16) and the resonance is equal to z I = E I -iΓ I . The modulus of the normalized resonant mode for B I = 0.1 eV is represented in Figure 2.

B I (eV ) N cv E D (meV ) E I (
The present section provides a numerical verification of the following theoretical results, which were used to write the One Mode Approximation: E D is close to E I , Γ I is small with respect to E I and the restriction of the resonant mode to the interval [0, L] is small outside the island [a 2 , b 2 ].

The stationary regime

We give here a numerical comparison of the methods presented in section 3: the Direct Resolution, Adaptive Method, One Mode Approximation and One Mode Approximation Integ methods. The Direct Resolution plays the role of the reference method. For all the tests, the reference potential in (3.1) is given in electron-volts by V ref = T kB q . Remark 5.1. For the bias B I = 0, the Gummel algorithm 3.1 can be initialized at the potential V 0 I = 0. For B I > 0, such an initialization does not converge. In that case, the convergence is obtained by taking the initial potential equal to the solution corresponding to B I = 0.

As it was underlined in section 3.1, the choice of the frequency mesh determines the convergence and the numerical cost of the method. In particular, for the Direct Resolution, it is verified numerically that, P ≥ 3700 is required to insure the convergence of the method for the biases B I = 0 eV and B I = 0.1 eV . This fact is illustrated in Figure 3 for B I = 0 eV . Indeed, for P = 1500, after 40 iterations the method fails to converge for a CPU time equal to 1.8241 s. For P = 3750, the method converges: the difference ∆ becomes less than tol = 10 -14 after 34 iterations for a CPU time equal to 3.8562 s. For the Direct Resolution, it appears that a big number of frequency points, and therefore a relatively high numerical cost, is required to have convergence. This remark enhances the importance of the other methods. The results given by the four methods, for two different values of the bias B I , are given in Figure 4 and Table 1. For the Adaptive Method, the small frequency step size is equal to ∆k = 2κ 4500 and the large step size is equal to 12∆k. The frequency step size is equal to the small one when the logarithmic derivative of the transmission coefficient, see [START_REF] Pinaud | Transient simulations of a resonant tunneling diode[END_REF][START_REF] Faraj | Méthodes asymptotiques et numériques pour le transport quantique résonnant[END_REF] for its explicit formula, is greater than 1.3L. The number of frequency points changes from an iteration to the other, however it stays around a fixed number which is 587 for B I = 0 eV and 535 for B I = 0.1 eV . For the other methods, the frequency mesh is given by the integers P = 3750 and P ′ = 150.

In Table 1, the integer N cv is the total number of iterations of Algorithm 3.1 for the tolerance tol = 10 -14 . The column CPU denotes the CPU time of the method and the column Dist to reference, the L 2 -norm relative distance to the reference method potential. The last is given by: 100

V Ncv -V di L 2 (0,L) V di L 2 (0,L) ,
where V di is the potential of the Direct Resolution. We notice that the alternative methods converge as fast as the reference method in term of number of iterations, and give a solution very close to the reference solution, for a smaller CPU time. For the alternative methods, by adequately fixing the parameters, and possibly allowing a bigger value in the column Dist to reference in Table 1, the frequency mesh can be chosen such that the CPU time is reduced. Anyway, for the frequency meshes chosen here, the alternative methods are an improvement to the reference method since their CPU times are lower than the CPU time of the Direct Resolution for the lowest value of P such that the convergence of the method is ensured for the biases considered. Although the One Mode Approximation and One Mode Approximation Integ methods require less frequency points than the Adaptive Method, the computation of the second term on the r.h.s. of equation (3.17) increases the numerical cost of the two first methods, such that the three alternative methods have an equivalent CPU time. On the other hand, compared to the other alternative methods, the Adaptive Method provides a self-consistent potential which is closer to the Direct Resolution potential (this is due to the fact that the step S1 is performed in the same way for the two last methods). However, in the present work, the One Mode Approximation is more interesting than the other alternative methods since, as remarked in section 4, the last can not be generalized to the time-dependent case.

The transient regime

For all the tests, the external potential U (t, x) is the time-dependent external potential corresponding to the bias B(t) defined in (4.1) where B I = 0 eV and B ∞ = 0.1 eV .

Direct Resolution: Toward optimal parameters

The Direct Resolution presented in section 4.1 will play the role of the reference method and its CPU time is essential to measure the effective improvement provided by our alternative method. In this context, we show in the present section which choice of the parameters of the Direct Resolution provides the smallest CPU time while conserving a reasonable accuracy. The initial data of the method is given by the stationary Direct Resolution for the initial bias B I = 0 eV and the parameters given in section 5.2 where P = 3750. The stationary potential corresponding to the final bias B ∞ is computed in the same way and will be denoted V ∞ . In Figure 5, we depicted, for different parameters of the method, the time evolution of the L 2 -norm relative distance to the potential V ∞ :

100 V l -V ∞ L 2 (0,L) V ∞ L 2 (0,L) , (5.1)
where V l is the Direct Resolution self-consistent potential corresponding to the iteration number l of Algorithm 4.1. For each graphic of Figure 5, are also given in the corresponding array: the total number of iterations l tot , the CPU time of the method and the Final dist which corresponds to the value of the quantity (5.1) for l = l tot .

In the first graphic of Figure 5, the number of frequency points is P = 3750 and we use the boundary conditions (C.9)(C.12) for the time-dependent Schrödinger equation where the angular frequencies are given by (C.13)(C.14). When the extrapolation of the intermediary potential V l+ 1 2 is given by (4.4), curve titled Pinaud, the Direct Resolution verifies the convergence criterion to the stationary potential V ∞ only for the time step ∆t = 10 -15 s. Although ∆t = 2 × 10 -15 s verifies the time step limit (4.2) imposed by the linear Schrödinger equation, the method is clearly unstable for this time step and the loss of stability comes from (4.4). When the extrapolation is given by (4.3), curve titled ABCN, the Direct Resolution potential convergences to V ∞ for both time steps ∆t = 10 -15 s and ∆t = 2 × 10 -15 s. Therefore, the extrapolation (4.3) provides a reduction of the CPU time. We remark however that the increase of the time step makes the amplitude of the transient oscillations of the potential slightly bigger.

In the second graphic of Figure 5, the number of frequency points is P = 3750 and we use the boundary conditions (C.9)(C.12) and the extrapolation (4.

3) for the intermediary potential V l+ 1 2 . When the angular frequencies are the discrete ones (C.13)(C.14), curve titled Disc AF, the Direct Resolution potential convergences to V ∞ and Final dist is small enough for both time steps ∆t = 10 -15 s and ∆t = 2 × 10 -15 s. When the angular frequencies are the continuous ones (C.8)(C.11), curve titled Cont AF, no instability is noticed, however, the value Final dist is much bigger than in the case Disc AF and increases significantly with the time step. In the case Disc AF, the computational cost is largely reduced while conserving a very satisfactory accuracy thanks to the simplified boundary conditions (C.22)(C.23) with M = 500, curve titled Disc AF-M = 500. Indeed, the CPU time is divided by more than two for ∆t = 10 -15 s and, for the two time steps, the value Final dist is less than 10 -2 . In the third graphic of Figure 5, we use the extrapolation (4.3) for the intermediary potential V l+ 1 2 and the boundary conditions (C.9)(C.12) where the angular frequencies are given by (C.13)(C.14). The different curves correspond to different values of the parameter P , curve titled P = 3750, P = 2100 and P = 450. The case P = 2100 illustrates the fact that a reasonable accuracy can be obtained even if the number of frequency points is under the limit P = 3700 which insures the convergence of the stationary method. This fact is relatively stable with respect to the time step and allows a reduction of the CPU time of the method. However, the number of frequency points can not be chosen too small: in the case P = 450, the value Final dist is bigger than 10 -1 and the amplitude of the transient oscillations of the potential are large compared to the case P = 3750. As in the stationary case, a big number of frequency points, and therefore a relatively high numerical cost, is required for the Direct Resolution, which enhances the importance of our alternative method.

The One Mode Approximation

Before comparing the two methods presented in section 4, we show how to adjust the parameters of the One Mode Approximation method, as we did in section 5.3.1 for the Direct Resolution. The initial data of the method is given by the stationary One Mode Approximation with the initial bias B I = 0 eV the parameters given in section 5.2 where P = 3900 and P ′ = 300. The stationary potential V ∞ corresponding to B ∞ is computed in the same way. For different parameters of the method, we depicted in Figure 6 the time evolution of the L 2 -norm relative distance (5.1) to the potential V ∞ where V l is the One Mode Approximation self-consistent potential corresponding to the iteration number l of Algorithm 4.1. For each graphic of Figure 6, we give in the corresponding array: the total number of iterations l tot , the CPU time of the method and the number Final dist defined in section 5.3.1. For all the curves represented, the time-dependent Schrödinger equation boundary conditions are given by (C.9)(C.12), the angular frequencies by (C.13)(C.14) and the extrapolation of the intermediary potential by (4.3). The number of frequency points corresponding to the thin mesh of our double scale algorithm is equal to P = 3900. In the first graphic of Figure 6, the number of frequency points corresponding to the coarse mesh is equal to P ′ = 300. When the angular frequency ω ∞ p in (4.22) is the discrete one (4.24), curve titled Oscill Interpolation, the One Mode Approximation verifies the convergence criterion to the stationary potential V ∞ for both time steps ∆t = 10 -15 s and ∆t = 2 × 10 -15 s. This is also true when the angular frequency ω ∞ p is the continuous one (4.23), curve titled Oscill Interpolation-Cont AF. However, it is clear from the column Final dist that the discrete angular frequency provides a better convergence to the stationary potential than the continuous one. The constant interpolation of the non resonant wave function, curve titled Const Interpolation, is obtained by taking the angular frequency ω ∞ p equal to 0 in (4.22). In that case, the One Mode Approximation is not accurate: the value Final dist is bigger than 10 -1 and the potential oscillates when the time is bigger than 2 × 10 -12 s. In the second graphic of Figure 6, the angular frequency ω ∞ p is given by (4.24). The different curves correspond to different values of the parameter P ′ , curve titled P ′ = 300, P ′ = 150 and P ′ = 39. The case P ′ = 150 illustrates the fact that, if we allow a lower accuracy, the CPU time can be reduced by reducing the number P ′ of Schrödinger equations to be solved. This fact is relatively stable with respect to the time step. However, it is clear from the case P ′ = 39 that P ′ can not be too small. As in the Direct Resolution, the extrapolation equation ( 4 frequency points, of the thin mesh, below the limit P = 3700.

Comparison for a suitable set of parameters

We consider now Figure 7 and Table 2 where is given a numerical comparison of the methods presented in section 4: the Direct Resolution and the One Mode Approximation. The Direct Resolution plays the role of the reference method.

For the Direct Resolution, resp. One Mode Approximation, the initial data of the method and the stationary potential V ∞ and density n ∞ corresponding to the final bias B ∞ are obtained as in section 5.3.1, resp. 5.3.2, where P = 3900 and P ′ = 300.

For the two methods, the boundary conditions for the time-dependent Schrödinger equation are given by (C.22)(C.23) with M = 500, the corresponding angular frequencies by (C.13)(C.14) and the extrapolation of the intermediary potential by (4.3). The frequency mesh is given by the integers P = 3900 and P ′ = 300, the total number of iterations is equal to l tot = 5 × 10 3 and the time step is equal to ∆t = 2 × 10 -15 s. The last verifies the condition (4.2). For the One Mode Approximation the angular frequency ω ∞ p is given by (4.24). In the graphic of Figure 7 entitled Convergence to stationary density, we depicted the time evolution of the L 2 -norm relative distance to the density n ∞ :

100 n l -n ∞ L 2 (0,L) n ∞ -n D L 2 (0,L) , (5.2)
where n l is the electron density corresponding to the iteration number l of Algorithm 4.1. Since the relevant variation of the density occurs inside the device, we considered in (5.2) the relative distance from (n l -n D ) to (n ∞ -n D ). In the graphic entitled Evolution of the resonant energy, we depicted the time evolution of the real part of the first resonance z l computed using the method presented in section 3.2 where the potential Q is equal to U l + V l . In the last two graphics of Figure 7, we represented the self-consistent potential and the probability density at final time step provided by the Direct Resolution, curve titled Direct, and the One Mode Approximation, curve titled One Mode. They are compared with the stationary Direct Resolution potential V ∞ , curve titled Direct-V ∞ , and density n ∞ , curve titled Direct-n ∞ . The same was done for the One Mode Approximation: curves titled One Mode-V ∞ and One Mode-n ∞ . In Table 2, the column Final dist is the value of the quantity (5.1) for l = l tot . The two methods behave the same qualitatively: the potential (resp. density) converges, up to a small error, to V ∞ (resp. n ∞ ). At time t = 0 + , the current is close to the current corresponding to B I . This is not true for the energy due to the discontinuity of the bias at t = 0. After fast transient oscillations, the current stabilizes in a short time to a value relatively close to the current corresponding to B ∞ whereas the energy increases slowly to such a value after slightly longer transient oscillations. The comparison is also favorable quantitatively. Indeed, at each time step, the currents provided by the two methods are relatively close of each other which, from formula (5.3), signifies that the One Mode Approximation wave functions are computed accurately in comparison with the Direct Resolution. Similarly, at each time step, the small relative difference between the energies provided by the two methods validates the One Mode Approximation selfconsistent potential in comparison with the Direct Resolution. We conclude that for the parameters chosen in this section, the One Mode Approximation provides a solution accurate enough at each time step with a reduction of the CPU time by a factor close to 6 compared to the Direct Resolution, see Table 2. The above comparison was not performed with the parameters of the reference method which give the best accuracy/CPU time ratio. However, when the parameters of the Direct Resolution are adjusted, as explained in section 5.3.1, to reduce the CPU time within a reasonable range of accuracy, the One Mode Approximation provides a solution in the same range of accuracy with a lower CPU time. 3) for M = 500, P = 3900, P ′ = 300, ∆t = 2 × 10 -15 s and l tot = 5 × 10 3 .

Remark 5.2. As it is proposed in [START_REF] Bonnaillie-Noël | Simulation of resonant tunneling heterostructures: numerical comparison of a complete Schrödinger-Poisson system and a reduced nonlinear model[END_REF], the average current density at time t is defined by:

J (t) = q m I R g(k)Im 1 L L 0 ∂ x Ψ k (t, x)Ψ k (t, x)dx dk . (5.3)
In (5.3), the space derivative is computed using the finite difference method. For the Direct Resolution, the integral with respect to the frequency k is computed using the trapezoidal rule where the frequency mesh is the mesh used to compute the density. For the One Mode Approximation, the frequency integral is given by a double scale trapezoidal rule of the form (4.21). As we did for the density, it is obtained by injecting in (5.3) the decomposition (4.14) of the wave function and neglecting the cross term between the resonant and the non resonant part.

Time evolution of the resonance peaks

For biases of the form (4.1), since the change of potential at time t = 0 is abrupt, the adiabaticity hypothesis of [START_REF] Presilla | Transport properties in resonant tunneling heterostructures[END_REF] is not satisfied and one expects two peaks corresponding to both the resonant energy at time t = 0 -and the resonant energy at time t. This is what we verify numerically in the present section by using the Direct Resolution and the One Mode Approximation.

The parameters and the initialization of the two methods are the same than in section 5. In Figure 8, for the curve with title Direct, the number of electrons in the interval [a 2 , b 2 ] carried by one wave function is given by:

N l di,p = j∈ w |Ψ l p,j | 2 ∆x , (5.4)
where Ψ l p,j is the Direct Resolution wave function defined in section 4.1 and w = {j | a 2 ≤ x j < b 2 }. For the curve with title One Mode, the cross term between the resonant and the non resonant part of the wave function is neglected and equation (5.4) becomes:

N l om,p = j∈ w Ψ nr,l p ′ ,j 2 ∆x + λ l p 2 j∈ w u l-1 2 j 2 ∆x , (5.5)
where Ψ nr,l p ′ ,j , λ l p and u l-1 2 j are the components, defined in section 4.2.2, of the One Mode Approximation wave function. Like in section 4.2.2, 0 ≤ p ′ ≤ P ′ is the integer such that p = νp ′ + r for some 0 ≤ r ≤ ν -1, in other terms, we made in (5.5) a constant interpolation of the non resonant wave function for p = νp ′ . Consider the numbers E l ∈ IR and Γ l > 0 such that

z l = E l -iΓ l ,
where z l is the Direct Resolution resonance, at time t l , defined in section 5.3.3. Then, the resonant frequencies

k l -= - 2m 2 (E l + B ∞ ) , k l + = 2m 2 E l (5.6)
are the real numbers related to the resonant energy E l according to the relation (4.13).

We can now discuss the results in Figure 8. For the Direct Resolution, we observe that for small times t l , the localization of the resonance peak is not given by the resonant energy E l but by the resonant energy at time t = 0 -: E I . Indeed, for such times, the number of electrons N l di,p defined in (5.4) has only one peak which is located at the frequencies such that E k = E I where E k and E I are given by (2.8) and (3.20) respectively (the numerical value of E I is computed as in section 5.1 where B I = 0 eV and P = 3900). When time goes, the first peak persists but decays slowly while a second peak appears and then grows at the frequencies k l -and k l + defined in (5.6). Using ] for one wave function with respect to the frequency k for the One Mode Approximation, and Resonant frequencies given by (5.6) (a vertical line is represented at these frequencies). In the two last graphics, is depicted a zoom in on k l -of the preceding curves at time t = 2 × 10 -12 s and t = 10 -11 s.

the definition of k l -and k l + , the last statement means that the second peak lives at the frequencies such that ε ∞ k = E l . At time t l = 10 -11 s, the initial peak has almost vanished and the second peak almost reached its maximal amplitude. For the One Mode Approximation, the first peak is provided by the initial condition and, as stated in section 4.2.2, the accuracy of the second peak is reached only if the interpolation of the non resonant wave function in (4.16) is well chosen. In particular, the constant interpolation computed as explained in section 5.3.2, curve titled One Mode-CI, provides a piecewise constant function close to the resonant frequencies k l -and k l + and, therefore, does not describe correctly the second peak (this is clear in the two last graphics of Figure 8). Whereas in the case of the oscillatory interpolation, which corresponds to the iteration (4.22), the coefficient λ l p is such that, for large times, the number of electrons N l om,p defined in (5.5) has a sharp peak at the frequencies k l -and k l + with amplitude high enough and the One Mode Approximation is accurate compared to the Direct Resolution.

Acknowledgements.The second author acknowledges Prof. Norbert Mauser for his hospitality in WPI Vienna where a part of the present work was achieved.

A Resolution of the stationary Schrödinger equation using finite difference

Discrete transparent boundary conditions are proposed in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF] to solve equation (2.5) on the domain [0, L] using finite difference. In order to be able to justify the choice of the transparent boundary conditions in section B, we recall how the method in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF] is deduced under the assumption k > 0 and we give an extension of this method to the case k < 0. For a given potential Q I (x) such that:

Q I (x) = 0 , x ≤ 0 and Q I (x) = Q I,L , x ≥ L (A.1)
and for a given wave vector k ∈ IR, we consider the stationary Schrödinger equation

- 2 2m ∂ 2 x Φ(x) + Q I (x)Φ(x) = E k Φ(x) , x ∈ IR (A.2)
with the k-dependent scattering conditions (2.6), (2.7).

Using the space discretization defined in section 3.1 and noting Φ j the approximation of the wave function Φ(x j ) and Q I,j the approximation of the potential Q I (x j ), the finite difference method for equation (A.2) writes:

D 2 x Φ j = 2m 2 (Q I,j -E k )Φ j , j = 1, ..., J -1 , (A.3) where D 2 x Φ j = Φ j+1 -2Φ j + Φ j-1 ∆x 2 .
For k > 0, the energy is given by E k = 2 k 2 2m . Then, thanks to assumption (A.1), the solutions to equation (A.3) can be computed explicitly for j ≤ 1. Indeed, under the condition ∆x < 1

k 2 + 2m 2 |Q I,L | , (A.4)
they correspond to the linear combinations of α j 1 , α j 2 where:

α 1,2 = 1 -∆x 2 k 2 2 ± i ∆x 2 k 2 -∆x 4 k 4 4 (A.5)
In order to work with quantities of order 1, the condition u H u = 1 is imposed instead of the norm condition in (3.7) and the resonant mode have to be divided by ∆x

M (z)u = 0 u H u = 1 , (B.9)
where

M (z) =        α 2 (z) -1 -1 2 + 2m∆x 2 2 (Q 1 -z) -1 0 . . . . . . . . . 0 -1 2 + 2m∆x 2 2 (Q J-1 -z) -1 -1 β 2 (z)        . (B.10)
We note that the holomorphy of M (z) required in section 3.2 is verified in S.

We give now an interpretation of the solutions to the problem (B.9). In the strip S, resonances can be defined at the discrete level as the complex numbers z such that the whole-space scheme (B.1) considered on j ∈ Z has purely outgoing solutions, i.e. solutions verifying (B.4) and (B.7). Such a solution is a called a resonant mode. Then the following properties hold (proof left to the reader):

• When z is a resonance, the corresponding resonant mode is uniquely determined up to a phase factor.

• A complex number z is a resonance if and only if z is such that there exists a vector u verifying (B.9).

• If z is a resonance, a corresponding u solution to (B.9) is equal, up to a phase factor, to the restriction of the resonant mode to 0 ≤ j ≤ J.

Proposition B.1. If ∆x verifies the condition (B.2), then for all z in S the coefficients defined in (B.3) and (B.6) verify:

Im (α 1 (z)) > 0 , Im (α 2 (z)) < 0 and Im (β 1 (z)) > 0 , Im (β 2 (z)) < 0 . (B.11)
If in addition Im (z) < 0, then it holds:

|α 1 (z)| > 1 , |α 2 (z)| < 1 and |β 1 (z)| > 1 , |β 2 (z)| < 1 .
In the next proposition, we show that resonances have a negative imaginary part. This, combined with Proposition B.1 and equations (B.4)(B.7), implies that the corresponding resonant modes tend to the infinity when j → ±∞. It is the discrete transcription of the fact that the space L 2 (IR) has to be deformed in order to consider resonances as eigenvalues.

Proposition B.2. If ∆x verifies the condition (B.2) and if z ∈ S is such that the problem (B.9) has a solution u, then it holds Im (z) < 0.

Proof. We multiply equation (B.1) by u j and sum up for j = 1, ..., J -1. Then, using the summation by part rule:

J-1 j=1 g j D - x f j = - J-1 j=0 f j D + x g j + 1 ∆x (f J-1 g J -f 0 g 0 ) , (B.12) where D + x f j = f j+1 -f j ∆x , D - x f j = f j -f j-1 ∆x , D - x D + x f j = D 2 x f j , we get J-1 j=0 D + x u j 2 = - J-1 j=1 2m 2 (Q j -z)|u j | 2 + 1 ∆x D + x u J-1 u J -D + x u 0 u 0 . (B.13)
Inserting the boundary conditions (B.5)(B.8) in (B.13) and taking the imaginary part, it follows:

0 = J-1 j=1 2m 2 Im (z)|u j | 2 - Im (β 2 (z)) ∆x 2 |u J | 2 - Im (α 2 (z)) ∆x 2 |u 0 | 2 .
Using u H u = 1 and (B.11), the previous equation implies Im (z) < 0.

C Resolution of the time-dependent Schrödinger equation

Although in the applications considered here the potential depends nonlinearly on the wave function, this section is written in the case of the linear Schrödinger equation. As in [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF] and [START_REF] Zisowsky | Discrete artificial boundary conditions for nonlinear Schrödinger equations[END_REF], the extension to the nonlinear case is provided by an adapted choice of the potential at half time step.

C.1 The homogeneous case

In this section, we recall the scheme proposed in [START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF] to solve on the bounded domain [0, L] the time-dependent Schrödinger equation

i ∂ t Ψ(t, x) = - 2 2m ∂ 2 x Ψ(t, x) + Q(t, x)Ψ(t, x) , t > 0 , x ∈ IR Ψ(0, x) = Φ(x) , x ∈ IR , (C.1)
with the following hypothesis:

H1. The initial condition Φ(x) is supported in 0 < x < L.

H2. The potential Q(t, x) verifies: for t > 0

Q(t, x) = 0 , x ≤ 0 and Q(t, x) = Q L , x ≥ L .
Considering the time and space discretization defined in section 4.1, we note Q l+1/2 j the approximation of the potential Q(t l+1/2 , x j ) and Ψ l j the approximation of the solution Ψ(t l , x j ). Then, equation (C.1), is solved with the Crank-Nicolson method:

i D + t Ψ l j = - 2 2m D 2 x Ψ l+1/2 j + Q l+1/2 j
Ψ l 1 -s 0 0 Ψ l 0 = l-1 k=1 s l-k 0 Ψ k 0 -Ψ l-1 1 , l ≥ 1 (C.4) Ψ l J-1 -s 0 J Ψ l J = l-1 k=1 s l-k J Ψ k J -Ψ l-1 J-1 , l ≥ 1 , (C.5)
where we have for j = 0 and j = J:

s l j = 1 -i R 2 + σ j 2 δ 0 l + 1 + i R 2 + σ j 2 δ 1
l + α j exp (-ilϕ j ) P l (µ j ) -P l-2 (µ j ) 2l -1 (C.6) and

ϕ j = arctan 2R(σ j + 2) R 2 -4σ j -σ 2 j , µ j = R 2 + 4σ j + σ 2 j (R 2 + σ 2 j )[R 2 + (σ j + 4) 2 ]
,

σ j = 2m∆x 2 2 Q j , α j = i 2 ((R 2 + σ 2 j )[R 2 + (σ j + 4) 2 ]) 1/4 exp i ϕ j 2 , R = 4m∆x 2 ∆t .
Here P l denotes the Legendre polynomials with the convention P -1 = P -2 = 0, δ j l the Kronecker symbol related to the integers j, l and Q 0 = 0, Q J = Q L . The coefficients s l j are computed thanks to the recursion formula derived in [START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF] from the standard recursion formula for the Legendre polynomials. We remark that the scheme (C.2) with the boundary conditions (C.4)(C.5) has a unique solution. By construction, see [START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF], the unique solution Ψ l j in L 2 (Z) to the whole-space scheme (C.2) considered on j ∈ Z verifies the problem (C. x + Q(t, x)]ϕ(t, x) , x ≤ 0 and verifies ϕ(0, x) = 0. Thus, we can write the homogeneous boundary condition (C.4) for ϕ and obtain the following boundary condition at x = 0 for Ψ:

Ψ l 1 -s 0 0 Ψ l 0 = l-1 k=1 s l-k 0 Ψ k 0 -Φ 0 e -iω0t k -Ψ l-1 1 -Φ 1 e -iω0t l-1
+ Φ 1 e -iω0t l -s 0 0 Φ 0 e -iω0t l , l ≥ 1 , (C.9) where the coefficients s l j are given by (C.6) and Φ j is a solution to the stationary discrete problem (A.3). We proceed similarly at x = L by setting ϕ(t, x) = Ψ(t, x) -Φ(x)e -iωLt , (C.10)

where

ω L = 1 (E k + Q L -Q I,L ) . (C.11)
This leads to the following boundary condition at x = L:

Ψ l J-1 -s 0 J Ψ l J = l-1 k=1 s l-k J Ψ k J -Φ J e -iωLt k -Ψ l-1 J-1 -Φ J-1 e -iωLt l-1
+ Φ J-1 e -iωLt l -s 0 J Φ J e -iωLt l , l ≥ 1 . (C.12)

C.2.2 Discrete angular frequencies, stability

Since the homogeneous transparent boundary conditions (C.4)(C.5) are derived from the fully discrete scheme (C.2) where j ∈ Z, the accuracy of the non-homogeneous boundary conditions (C.9)(C.12) is improved under the condition that the functions in (C.7) and (C.10) are solutions to the discrete equation (C.2) for j ≤ 0 and, respectively, j ≥ J. As shown in [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF], this leads to the non-homogeneous boundary conditions (C.9)(C.12), where the angular frequencies given in (C. We remark that, for k > 0, the solution to the initial finite difference scheme (A.3) considered on j ∈ Z with the scattering conditions (A.6)(A.9) is in L ∞ (Z) but not in L 2 (Z). The same property is true for k < 0. Therefore, for such an initial condition, the suitable state space for the solutions to the scheme (C.2) considered on j ∈ Z is L 2 (Z) + L ∞ (Z). In particular, it is easy to show that this scheme has a unique solution Ψ l j such that Ψ l j j∈Z ∈ L 2 (Z) + L ∞ (Z) for l ≥ 0. By construction, Ψ l j verifies the problem (C.2) with the boundary conditions (C.9)(C.12) where the angular frequencies are given by (C.13)(C.14). It follows that the unique solution to (C.2)(C.9)(C.12) corresponds exactly with the restriction of Ψ l j to 0 ≤ j ≤ J. Moreover, we have the following stability result: Theorem C.2. Let Φ j be the solution to the scheme (A.3) with the boundary conditions (A.7)(A.10) for k > 0 and (A.11)(A.12) for k < 0. Let Ψ l j be the solution to the discrete Schrödinger equation (C.2) with the non-homogeneous boundary conditions (C.9)(C.12) and the initial condition Ψ 0 j = Φ j . We suppose in addition that the angular frequencies ω 0 and ω L in (C.9)(C.12) are given by (C.13)(C.14). Then, there exists a constant C > 0 which depends only on , m, L, Q I,L , sup x∈(0,L) |Q I (x)|, Q L , k such that for all ∆x verifying (A.4) it holds for all N ≥ 0:

Ψ l 2 ≤ C 1 + max 0≤n≤N Q n+1/2 ∞ T + 1 , 0 ≤ l ≤ N , (C.15)
where T = N ∆t and Q n+1/2 ∞ = max 1≤j≤J-1 |Q n+1/2 j |.

Even if no size limit on the time step is deduced from Theorem C.2, the time oscillating terms e -iω0t and e -iωLt in (C.9) and, respectively, (C.12) require the following condition for a reasonable resolution:

|ω j |∆t < π 2 , j = 0, L . (C.16)
We can now give the proof of Theorem C.2.

Proof. Without any practical restriction, we can suppose that J ≥ 3. Let us consider a smooth function χ(x) such that 0 ≤ χ(x) ≤ 1 , x ∈ IR , χ(x) = 1 , x ≤ L 3 , χ(x) = 0 , x ≥ 2L 3 .

Using the vector χ j = χ(x j ), we define for l ≥ 0 and 0 ≤ j ≤ J: θ l j = χ j Φ j e -iω0t l + (1 -χ j )Φ j e -iωLt l and ϕ l j = Ψ l j -θ l j . Using the scheme (C.2) for Ψ l j , we obtain the discrete Schrödinger equation below for the function ϕ l j : i D + t ϕ l j = - , where we used the summation by part rule (B.12). Taking the imaginary part, we get:

D + t ϕ l 2 2 = m -Im D - x ϕ l+1/2 J ϕ l+1/2 J + Im D + x ϕ l+1/2 0 ϕ l+1/2 0 + 2∆x J-1 j=1
Im F l j ϕ l+1/2 j .

The initial condition ϕ 0 = 0 obviously holds, therefore, we have by summing up with respect to l in the previous equation: for all N ≥ 0. Since the non-homogeneous boundary conditions (C.9)(C.12) hold for Ψ l j , it follows that the homogeneous ones (C.4)(C.5) hold for ϕ l j . Therefore, the boundary terms in (C.18) behave like in Theorem C.1. In particular, it is shown in [START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF] Then, it follows from (C.18) that:

ϕ N +1 2 2 ≤ 2∆t N l=0 F l 2 ϕ l+1/2 2 .
(C. [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF]

Algorithm 3 . 1 (

 31 Gummel algorithm). Fix l = 0. Set ∆ = tol + 1. Do While ∆ ≥ tol: Computation of the density S1. For p = 0, ..., P : computation of the wave function Φ l p at the frequency k p from the potential V l I by solving equation (2.5) with V I = V l I . S2. Numerical integration: computation of the density n l I from the wave functions Φ l p 0≤p≤P by using equation (2.10). Coupling to the Poisson equation : Gummel iteration Computation of the potential V l+1 I from the potential V l I and the density n l I by solving equation (3.1).

7 )

 7 has solutions which are purely outgoing outside the interval [0, L]. The function u(x) in (3.7) is called resonant mode associated to the resonance z. As described in Appendix B, the eigenproblem (3.7) is restricted (for numerical computations) to the domain [0, L] by the aid of transparent boundary conditions selecting only outgoing states. In particular, this problem is discretized using the finite difference scheme (B.1) inside the domain and the boundary conditions (B.5), (B.8) at x = 0 and x = L. Since the boundary conditions depend nonlinearly on the spectral parameter z, this process leads to the following nonlinear eigenvalue problem: find (u, z) ∈ C J+1 × C such that

  .1) where B I , B ∞ ≥ 0 and B I = B ∞ . The case of the general bias requires a different version of the discrete transparent boundary conditions for the resolution of the transient Schrödinger equation.

Algorithm 4 . 1 (

 41 Transient algorithm).

  x) and L 0 |u(t, x)| 2 dx = 1 . (4.10) If u I (x) and θ k are respectively the resonant mode and the proportionality coefficient given by (3.13) and (3.15), then it holds from (3.11)(3.14):

. 14 )

 14 Then the step S1 writes as follows: suppose that the quantities Ψ nr,l k and λ l k are known (at l = 0, it is given by the initial decomposition (4.11)). As it is done in the Direct Resolution for Ψ l+1 k , the non resonant wave function Ψ nr,l+1 k at time t l+1 and frequency k is computed on [0, L] using Ψ nr,l k and V l+ 1 2 by the resolution of (4.8) with the Crank-Nicolson scheme (C.2) where the potential Q is equal to U f ill + V and the initial data Φ is equal to Φ nr k . Equation (C.2) comes with the transparent boundary conditions (C.9)(C.12) or (C.22)(C.23).

Rel× 10 - 21 Jdensity n 1 D 10 24 m - 3 Donor density n 2 D 5 ×

 1021110325 Donor 10 21 m -3

Figure 2 :

 2 Figure 2: Representation of the potential U I + V I (dashed line) and the corresponding normalized resonant mode |uI (x)| |uI |2 (full line) for B I = 0.1 eV .

Figure 3 :

 3 Figure 3: Potential difference ∆ defined by (3.3) with respect to the number of iterations of Algorithm 3.1 for the Direct Resolution for B I = 0 eV and for the values P = 1500 and P = 3750.

  [START_REF] Ben Abdallah | Multiscale simulation of transport in an open quantum system: Resonances and WKB interpolation[END_REF] for B I = 0 eV and B I = 0.1 eV . 4e-08 6e-08 8e-08 1e-07 1.2e-07 1.4e-

Figure 4 :

 4 Figure 4: For the biases B I = 0 eV and B I = 0.1 eV and, for the methods considered, are represented from left to right: the external potential U I , the corresponding potential difference ∆ defined by (3.3) with respect to the number of iterations of Algorithm 3.1, the self-consistent potential V I and the electron density n I .

Figure 5 :

 5 Figure 5: Evolution of the potential relative distance (5.1) to the stationary solution for the Direct Resolution. From top to bottom: Impact of the potential extrapolation, Impact of the boundary conditions and Impact of the frequency mesh.

  of the non resonant wave function interpolation Oscill Interpolation, ∆t=10 -15 s Oscill Interpolation, ∆t=2x10 -15 s Oscill Interpolation-Cont AF, ∆t=10 -15 s Oscill Interpolation-Cont AF, ∆t=2x10 -15 s Const Interpolation, ∆t=10 -15 s ltot CPU(s) Final dist (%) Oscill Interpolation, ∆t = 10 -15 s coarse frequency mesh P'=300, ∆t=10 -15 s P'=300, ∆t=2x10 -15 s P'=150, ∆t=10 -15 s P'=150, ∆t=2x10 -15 s P'=39, ∆t=10 -15 s ltot CPU(s) Final dist (%) P ′ = 300, ∆t = 10 -15 s 10 4 496.03 5.5306 × 10 -3 P ′ = 300, ∆t = 2 × 10 -15 s 5 × 10 3 204.29 6.4993 × 10 -3 P ′ = 150, ∆t = 10 -15 s 10 4 344.06 0.12238 P ′ = 150, ∆t = 2 × 10 -15 s 5 × 10 3 152.50 0.12444 P ′ = 39, ∆t = 10 -15 s

Figure 6 :

 6 Figure 6: Evolution of the potential relative distance (5.1) to the stationary solution for the One Mode Approximation. From top to bottom: Impact of the interpolation of the non resonant wave function and Impact of the coarse frequency mesh.

Figure 7 :

 7 Figure7: For M = 500, P = 3900, P ′ = 300, ∆t = 2 × 10 -15 s, l tot = 5 × 10 3 and, for the methods considered, are depicted from left to right: Evolution of the resonant energy, Evolution of the potential relative distance (5.1) to the stationary solution, Evolution of the average current density defined by (5.3), Evolution of the density relative distance (5.2) to the stationary solution, Self-consistent potential at final iteration l tot , and Electron density at final iteration l tot .

3 . 3 .

 33 Since we are not interested in the reduction of the CPU time, we took the exact boundary conditions (C.9)(C.12) for the time-dependent Schrödinger equation, instead of the simplified ones (C.22)(C.23), and the time step equal to ∆t = 10 -15 s, instead of ∆t = 2 × 10 -15 s.

Figure 8 :

 8 Figure 8: For P = 3900, P ′ = 300, ∆t = 10 -15 s and at different values of the time t: Number (5.4) of electrons in the interval [a 2 , b 2 ] for one wave function with respect to the frequency k for the Direct Resolution, Number (5.5) of electrons in the interval [a 2 , b 2] for one wave function with respect to the frequency k for the One Mode Approximation, and Resonant frequencies given by (5.6) (a vertical line is represented at these frequencies). In the two last graphics, is depicted a zoom in on k l -of the preceding curves at time t = 2 × 10 -12 s and t = 10 -11 s.

1 2

 2 to verify J j=0 |u j | 2 ∆x = 1. Then, it follows from equations (B.1), (B.5) and (B.8) that the sought discrete problem is

Ψ l+1/ 2 j 2 )

 22 , j = 1, ..., J -1 , l ≥ 0 , comes with the discrete transparent boundary conditions:

  2)(C.4)(C.5). It follows that the solution to (C.2)(C.4)(C.5) corresponds exactly with the restriction of Ψ l j to 0 ≤ j ≤ J. Moreover, we have the following stability result: Theorem C.1 ([20]). The solution to the discrete Schrödinger equation (C.2) with the boundary conditions (C.4)(C.5) is uniformly bounded:Ψ l 2 ≤ Ψ 0 2 , l ≥ 1and therefore, the scheme is unconditionally stable.C.2 The non-homogeneous caseC.2.1 The non-homogeneous discrete transparent boundary conditionsWe consider the problem (C.1) where the initial condition Φ(x) is solution to (A.2). We suppose that the potentials Q I (x) and Q(t, x) are such that (A.1) and H2 hold and Q(0, x) = Q I (x). The initial condition does not verify assumption H1, therefore (C.4)(C.5) have to be replaced by suitable boundary conditions in the Crank-Nicolson scheme (C.2). Using equations (C.1) and (A.2), the function ϕ(t, x) = Ψ(t, x) -Φ(x)e -iω0t ,

  8) and (C.11) are replaced by the discrete ones below: + Q L -Q I,L ) . (C.14)

  j , j = 1, ..., J -1 , l ≥ 0 are defined as in (C.3). Multiplying equation (C.17) by ϕ l+1/2 j and summing up for j = 1, ..., J -1, it follows:

  that the boundary conditions (C.4)(C.5) imply:

Table 1 :

 1 Comparison of the Direct Resolution, Adaptive Method, One Mode Approximation and

			34	3.8642	/
		Adaptive One Mode	34 34	0.6040 0.6440	1.3859 × 10 -2 0.2561
		One Mode Integ	34	0.6040	0.3700
	B I = 0.1 eV	Direct Adaptive One Mode	31 31 31	3.3402 0.4800 0.5560	/ 4.9157 × 10 -3 0.7689
		One Mode Integ	31	0.4880	0.5126

One Mode Approximation Integ methods for the resolution of the stationary Schrödinger-Poisson system (2.5)(2.9)(2.

Table 2 :

 2 Comparison of the Direct Resolution and One Mode Approximation for the resolution of the time-dependent Schrödinger-Poisson system (2.1)(2.2)(2.

* Research supported by the project QUATRAIN (BLAN07-2 212988) funded by the French Agence Nationale de la Recherche and from the Marie Curie Project DEASE: MEST-CT-2005-021122 funded by the European Union.

are complex numbers verifying |α j | = 1 and α 1 = e i k∆x for some k > 0. Then α j 1 is an incoming plane wave, α j 2 = α -j 1 is the reflected wave and the discrete analogous to the scattering condition (2.6) is to select the linear combinations of α j 1 , α j 2 of the form:

Applying equation (A.6) to compute Φ 0 and Φ 1 , we get the following discrete transparent boundary condition at x = 0:

We note that inequality (A.4) implies k∆x < π 2 and the discretization provides enough space points in one wave length of the discrete plane wave. Similarly, for j ≥ J it holds Q I,j = Q I,L and, if k 2 -2m 2 Q I,L > 0, the solutions to equation (A.3) are the linear combinations of β j 1 , β j 2 where:

Inequality (A.4) implies |β j | = 1 and β 1 = e i k∆x for some k > 0. Then β j 1 is an outgoing wave and the scattering condition (2.6) gives:

Applying equation (A.9) to compute Φ J-1 and Φ J , we get the following discrete transparent boundary condition at x = L:

Again, condition (A.4) implies k∆x < π 2 . If k 2 -2m 2 Q I,L ≤ 0, the bounded solutions to (A.3) for j ≥ J -1 are given by (A.9), where Φ j is constant in the case of the equality and vanishing in the case of the strict inequality, and (A.10) still apply. For k < 0, the energy is given by E k = 2 k 2 2m + Q I,L . Then, using the assumptions (A.1), (A.4) and repeating the above calculations where the scattering condition (2.6) is replaced by (2.7), we obtain the following transparent boundary conditions for equation (A.3) at x = 0 and x = L:

where α 2 is given by equation (A.5) and

For k > 0, the scheme (A.3) with the boundary conditions (A.7)(A.10) has a unique solution. By construction, the unique solution Φ j to the whole-space scheme (A.3) considered on j ∈ Z with the scattering conditions (A.6)(A.9) verifies the problem (A.3)(A.7)(A.10). It follows that the solution to (A.3)(A.7)(A.10) corresponds exactly with the restriction of Φ j to 0 ≤ j ≤ J. A similar property holds for k < 0. Moreover we have the following estimate on the solution:

Lemma A.1. There exists a constant C > 0 which depends only on , m, L, Q I,L , sup x∈(0,L) |Q I (x)|, k such that for all ∆x verifying (A.4), the solution to the scheme (A.3) with the boundary conditions (A.7)(A.10) for k > 0 and (A.11)(A.12) for k < 0 verifies:

where

B Computation of resonances using finite difference

For a given potential Q(x) verifying (3.6), we give in this section a discrete version of the following problem: find (u, z) solution to (3.7) such that u is purely outgoing outside the interval [0, L].

Using the decay at infinity of the function g(k) in (2.4), it follows that only the resonances with real part smaller than

2m are important for the computation of the density, where κ is defined in section 3.1. Then, if we suppose that the resonances have a positive real part, which is not a restriction in the applications, the problem of the computation of resonances can be restricted to the strip

In order to work with analytic functions on S, the value of the square root function √ z will be the principal value. Proceeding as in section A, we write first the finite difference method for equation (3.7):

where u j denotes the approximation of the resonant mode u(x j ) and Q j the approximation of the potential Q(x j ). Thanks to assumption (3.6), the solutions to equation (B.1) can be computed explicitly for j ≤ 1. Indeed, under the condition ∆x < 1

they correspond to the linear combinations of α 1 (z) j , α 2 (z) j where:

Using the forthcoming equation (B.11), it holds α 2 (z) = ρe ik∆x with ρ > 0 and k < 0. Since we are looking for solutions which are outgoing outside [0, L], we impose:

Applying equation (B.4) to compute u 0 and u 1 , we get the following discrete transparent boundary condition at x = 0:

Similarly, for j ≥ J it holds Q j = Q L and the solutions to equation (B.1) are the linear combinations of β 1 (z) j , β 2 (z) j where:

Using equation (B.11), the outgoing wave is identified with the plus sign in (B.6), therefore we impose:

Applying equation (B.7) to compute u J-1 and u J , we get the following discrete transparent boundary condition at x = L:

Since β 1 (z)β 2 (z) = 1, the previous equation writes:

Using the equations (C.13)(C.14) and the equation (A.3), with the boundary conditions (A.7)(A.10) for k > 0 and (A.11)(A.12) for k < 0, and using the result of Lemma A.1, the following estimate holds for the source term:

where C denotes a constant which depends only on , m, L, Q I,L , sup x∈(0,L) |Q I (x)|, Q L , k. Then, after an easy computation, the inequalities (C. [START_REF] Chang | Difference Schemes for Solving the Generalized Nonlinear Schrödinger Equation[END_REF]) and (C.20) lead to:

where T = N ∆t. The inequality (C.15) is obtained by applying (C.21) and Lemma A.1 again, together with the following inequality

C.2.3 Simplified discrete transparent boundary conditions

As proposed in [START_REF] Ehrhardt | Discrete Transparent Boundary Conditions for the Schrödinger Equation[END_REF], the computational time required by our transparent boundary conditions can be reduced by using the decay of the coefficients s n j appearing in (C.9)(C.12). This improvement is obtained by replacing the coefficients s n j by 0 after an index M ≥ 1 chosen big enough to preserve the accuracy of the boundary conditions. In that case, the numerical cost of the convolutions in (C.9)(C.12) is reduced and the computations are performed as follows: for 1 ≤ l ≤ M + 1 the boundary conditions are given by (C.9)(C.12) and for l ≥ M + 2, they are given by:

-Φ 1 e -iω0t l-1 + Φ 1 e -iω0t l -s 0 0 Φ 0 e -iω0t l (C.22) and

J-1 -Φ J-1 e -iωLt l-1 +Φ J-1 e -iωLt l -s 0 J Φ J e -iωLt l .

(C.23)