P Godignon

J Montserrat
email: bjosep.montserrat@cnm.es

J Rebollo
email: cjose.rebollo@cnm.es

D Planson

Edge terminations for 4H-SiC power devices: a critical issue

Keywords: edge termination, JTE, breakdown voltage, OBIC

Edge termination is a critical part of a power devices. Numerous edge termination types have been developed for silicon devices. Implementation of these termination architectures are not straightforward in SiC due to physical and processing specificities: lower junction depths, higher electric field, trench depth and shaping limitations, etc. Two main families of terminations are currently used in commercial devices, pure Field Guard Rings, and JTE + Rings combination. The increasing number of trench commercial devices requires new approaches based on etched rings filled with dielectrics or polysilicon. For epitaxied bipolar devices, MESA with bevel angle termination combined with JTE based architecture are also suitable. In any case, and especially regarding avalanche capability requirements, not only the termination architecture is relevant, but also the passivation type, the channel stopper design, the 3D design. As modelling using conventional tools is not fully reliable, specific characterization methods are needed. For instance, micro-OBIC can be very effective to determine the electric field distribution in the periphery of the power devices.

Introduction

Most of planar and trench power devices are based of one or several n-type/p-type junctions. When the active part of the junction ends on the lateral side of the device, a very high electric field appear under reverse voltage biasing, due to junction curvature or trench corners. Then, a periphery protection, called edge termination is necessary. Reliable and robust high voltage devices need effective edge termination structures to protect the device periphery, so that blocking values close to the ideal 1D avalanche voltage value can be achieved. Many types of terminations have been developed in Silicon technology, the main ones being MESA termination for rectifiers and thyristors, field plates for low voltage range, field guard rings (FGR) and junction termination extension (JTE) for planar rectifiers and transistors. Silicon Carbide is not as friendly as Silicon in terms of processing technology for power device fabrication. Several physical limitations make SiC edge terminations more complex to design in a same way to Silicon. The main limitation is the low diffusion coefficients of dopant atoms inside the semiconductor. Deep junction (> 1μm) are difficult to implement and lateral dopants diffusion is very limited. On the other hand, one of the theoretical advantages of SiC versus Si, its high critical electric field, is also a drawback when looking for compatible dielectric layers capable of operating under such high fields in a reliable way. It clearly limits the use of field plates for instance. Finally, the higher interface traps level observed in the SiC/SiO2 interface compared to Si is also affecting the efficiency of the standard Si terminations architectures in the SiC power devices. Charges accumulate at the interface, generating breakdown voltage instabilities, especially with temperature. Finally, the termination design is also relevant to reach high avalanche ruggedness in order to comply with applications prone to unclamped inductive switching. In this case, the edge termination must reach breakdown at higher voltage than the main active junction [1].

Edge termination structures

Termination optimization is relevant in power device design since it may affect the reliability of the component, but also because it may consume a large area of the die, especially in the case of very high voltage devices (Fig. 1). This is particularly relevant for SiC semiconductor, considering the cost of the starting material and the defects density, which usually limit the die area of the components. In addition, in some cases, edge termination integration increases the fabrication process complexity (increasing the number of photolithography mask levels), like in trench JFETs, pure Schottky diodes, BJTs or thyristors. These aspects have a impact on cost and reliability. Usual terminations architectures. The more common periphery protection used in SiC commercial devices is the field guard ring (FGR) termination (Fig. 2a). It is formed with P + -type implantation wells, having a width of few microns. The key design parameters are the rings distances, which usually increase with the ring position, and the first ring distance with the device main junction. This termination usually consumes a large area and require very precise photolithography and mask etching resolution. The optimal first ring distance depends on the epilayer doping and must be lower than 1μm for devices voltage ranges below 1.7kV. Rings termination also depends on the surface charges at the SiC/dielectric interface. In this sense, buried rings configuration have been proposed in order to shift the high electric field peaks from the surface to the SiC bulk [START_REF] Mihaila | th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD06)[END_REF], limiting interface parasitic effects. This is efficient but more complex to implement, since epilayer regrowth or very high energy implantation are required.

The second common type of edge protection is the Junction Termination Extension (JTE), which enable smaller area, but which typically requires an extra Al ion implantation step, with an accurate control of the implanted dose. JTE is also sensible to surface charges [START_REF] Matsushima | th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD17)[END_REF]. Then, single implantation JTE is not used anymore, as multiple JTE combination, having decreasing implantation dose, have shown to be more efficient in terms of process control. SiC etching process can be also use to mitigate the JTE dose (see Fig. 2b) and create a double or triple JTE structure without the need for 3 Al implantation steps [START_REF] Dai | 33 th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD21)[END_REF]. Even if the number of photolithographic masks is not reduced, the global technological process results simpler and cheaper. However, a combination of JTE + FGR, or JTE + FGR + JTE outer rings (see Fig. 3) is the most stable, reproducible and efficient combination for most of the SiC devices structures, especially for blocking voltages higher than 1.7kV. For instance, this architecture has been used in [START_REF] Kaji | [END_REF] for 27 kV diodes demonstration. Other combination of JTE and rings have been proposed in order to reduce termination length [START_REF] Mihaila | th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD16)[END_REF]. In a previous paper [START_REF] Soler | [END_REF], we analyzed different field guard rings, JTE and JTE + rings combinations for 1.7kV and 4.5kV planar MOSFETs fabrication. Main results are summarized in Table 1. We intended to use the MOSFET P-well doping to form deep FGR termination (Fig. 1). The efficiency (calculated using simulated ideal 1D breakdown) of these FGR terminations (designs D3 to D8) is lower than the JTE + Rings one (designs D9 to D11). We believe this difference is due to the surface charges type present in our devices. We also tried to use the P-well doping profile as the 1 st JTE in a double JTEs configuration (design D2), without clear success. The lateral extension of the P-well in the periphery has basically no impact on the termination efficiency. Other 2D and 3D parameters must be defined when designing a JTE (+ rings) termination. JTE length or corner radius must be optimized to minimize the consumed area. In table 2 are reported breakdown voltages of 1.7kV JBS diodes with the termination architecture shown in Fig. 3a. When comparing design D5 to the reference, we can see that a reduction of the corner radius decreases the breakdown capability. The corner radius will also have a strong impact on the device avalanche ruggedness. On the other hand, the recent availability of SiC trench devices also required adaptation of edge termination. In this sense, 3D edge termination like trench guard rings terminations have been recently proposed for SiC devices [8]. This termination configuration has shown to be efficient in Silicon technologies. It allows to reduce the consumed area, and is compatible with the trench fabrication process. Additionally, the bottom of the trench can be implanted with self-aligned buried P + ring (see Fig 4a), allowing a reduction of the electric field at the SiC surface. In Fig. 4b, we can see another alternative using a wide trench filled with Benzocyclobutene (BCB) in an as JTE configuration [9][START_REF] Liu | 31 st Int. Symp. on Power Semiconductor Devices and ICs (ISPSD19)[END_REF]. The critical part of this termination is the trench side wall charges displacement and the possible parasitic channel that can be created with the electric field peaks. The field plate geometry on top of the termination structure is also a key parameter. For high voltage components, especially for bipolar SiC device, positive and negative beveled MESA terminations can also be used. These terminations, formerly used in thyristor and early power diodes components presents the advantage of lower dimensions than planar terminations [START_REF] Sung | [END_REF]. The main issue, especially for the negative bevel, is the SiC etching process necessary to reach low aspect ratio angles. However, solution playing with the etching mask topology (photoresist) or using multiple etching steps are possible, even if not fully compatible with mass production requirements. Other termination elements. One element of the termination never mentioned in the literature is the channel stopper (see Fig. 5a). Nevertheless, it is part of the termination design, as the distance between the end of the JTE or the rings and this channel stopper may affect the breakdown capability, the dynamic behavior and directly contribute to the total area required for the periphery protection of the power devices. To determine the optimal position of the channel stopper by modelling is tricky as it depends on a wide range of parameters such as the interface traps and fixed charge of the different passivation layers, the type of packaging technique used, the position of the top metal electrode, etc. We studied the optimal position of the channel stopper in 1.7 kV JBS diodes passivated with 5µm of polyamide. Diodes with a distance of 100µm between the end of the JTE and the channel stopper were able to withstand more than 2 kV while the diodes with a JTE-CS distance of 60µm exhibit a premature breakdown around 1.6 kV after arcing between the anode and channel stopper metals. A picture of a failing diode at 1.6 kV is shown in Fig. 5b, where the damaged area is clearly seen. Last but not least, the passivation of the surface in the termination area is a critical aspect for the wide band gap semiconductor high voltage devices. When no passivation is used, premature breakdown due to arcing in air can occur. However, the addition of a passivation layer (for instance polyamide) may change the surface charges state, and impact the breakdown voltage. This is particularly visible in JTE terminations, as the surface charges directly impact the optimal value of JTE implantation dose. Again, a combination of JTE and rings is the most suitable solution to mitigate these charges issues. For high voltage devices (> 2.5 kV), additional effects due to ambient factors like humidity are enhanced and may seriously affect the reliability of power modules if not properly passivated [START_REF] Hoffman | th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD18)[END_REF]. In this case, solutions pass through a proper packaging technology, for instance using organic passivant like GlobTop. Stability and reliability evaluation of the termination can be done using High Temperature Reverse Bias test (HTRB) [START_REF] Hoffman | th Int. Symp. on Power Semiconductor Devices and ICs (ISPSD18)[END_REF], including Humidity (H3TRB) [START_REF] Jormanainen | [END_REF][14], if a proper evaluation of the passivation is to be done. It is also relevant to test the termination under Unclamped Inductive Switching (UIS) to check the avalanche capability of the device [1].

Reference

Termination optimization through advanced characterization:

A good understanding of termination behavior also requires advanced characterization tools such as lock-in thermography, photoluminescence, Electron Beam Induced Current (EBIC). Optical Beam Induced Current technique (OBIC) applied to WBG semiconductor devices have shown to be useful to study the efficiency of different edge protection techniques [15]. In order to illustrate some OBIC experiments on different terminations, PiN diodes with several JTE protections (1.7 kV-class) has been measured with increasing reverse voltage up to 400V. Since OBIC setup is under air, as described in [16], it limits the maximum voltage bias for unpassivated devices. Scans were performed from the center of the diode towards the outer periphery. For each results of Fig 6, one can find the drawing of the edge termination at the bottom. Diodes (#1-#3) exhibit increasing length of JTE. Diode #4 has a JTE length of 35 µm + 3 JTE rings. As it can be seen, the OBIC signal appears at least on the JTE length even with 0V. As the reverse voltage increases, the OBIC signal increases and extends laterally towards the outside. For diode #4, none of the JTE rings are depleted at 0V. However, these rings are quickly depleted with a reverse voltage as low as 50V.

Conclusion

Current edge termination in SiC can reach 100% efficiency, which is relevant for the avalanche ruggedness capability. Termination design is still improving to adapt to the new Trench and Super-Junction technologies or to reduce the consumed area, especially in high and very high voltage devices. However, physical limits as the minimum distance between main top electrode and channel stopper will fix the maximum size reduction achievable. Passivation and packaging are also two of the main concerns for the next generation of high voltage devices.

Fig. 1 :

 1 Fig. 1: 6.5kV commercial device picture and evaluation of its respective active and termination area

Fig. 2 :

 2 Fig. 2: Schematic of typical termination types including (a) guard rings termination using VDMOS pwell (b) triple JTE created by SiC dry etching

Fig. 3 :

 3 Fig. 3: Schematic of JBS diode termination including (a) JTE with JTE outer rings and P-well FGR (b) comparison of simulated breakdown voltage of a 6.5kV wafer for different JTE + rings designs

Fig. 4 :

 4 Fig. 4: Schematic of a) a trench rings termination [8] (b) trench termination filled withBenzocyclobutene BCB[START_REF] Liu | 31 st Int. Symp. on Power Semiconductor Devices and ICs (ISPSD19)[END_REF]

Fig. 5 :

 5 Fig. 5: (a) schematic of the termination with channel stopper (b) A picture of a failed 1.7 kV diode after premature breakdown at 1.6 kV due to short distance between anode and channel stopper.

6 :

 6 PiN Diode #1 protected with a 15 µm JTE PiN Diode #2 protected with a 40 µm JTE PiN Diode #3 protected with a 50 µm JTE PiN Diode #4 with a 35 µm JTE and 3 JTE rings Fig. OBIC measured current versus distance in experiments on 1.7 kV Class PiN diodes with several JTE protections geometries. Basic drawing of each termination is shown as inset.

Table 1 :

 1 Measured breakdown voltage and termination efficiency of 1.7kV & 4.5kV MOSFET[START_REF] Soler | [END_REF]

			Wafer#1		Wafer#2
	Design	Schematic cross section	Average VBR	Efficiency (%)	Averag e VBR	Efficiency (%)
			(V)		(V)
	D1 (130 μm)		1930	85	4500	77
	D2 (100 + 30μm)		1870	82	4000	68
	D3 (12 rings)		2000	88	4870	83
	D4		1360	60	4550	78
	D5		960	42	2880	49
	D6 (20 rings)		2080	91	4550	78
	D7		1300	57	4230	72
	D8		890	39	2700	46
	D9		2260	99	5050	86
	D10		2260	99	5130	88
	D11		2270	99	5230	89

Table 2 :

 2 Impact of the JTE length and corner radius in 1.7kV JBS diode avalanche voltage

		design	D2	D3	D4	D5
	JTE length	150μm	100μm	50μm	35μm	150μm
	Corner radius	150μm	150μm	150μm	150μm	100μm
	Avalanche voltage	> 2000V	> 2000V	1450V	1320V	1910V

Acknowledgments: Autors want to acknowledge Dr. V. Soler for his contribution to this work. This work was partially supported by the DiamMOS project funded by the Spanish Government (TEC2017-86347-C2-2-R), cofounded by the EU-ERDF (FEDER), and by the Agència de Gestió d'Ajuts Universitaris i Recerca Funds under Grant 2017-SGR-1384.


Arnaud Leclerc
email: arnaud.leclerc@univ-lorraine.fr

Georges Jolicard
email: georges.jolicard@utinam.cnrs.fr

Global integration of the Schrödinger equation: a short iterative scheme within the wave operator formalism using discrete Fourier transforms

Keywords: numbers: 31.15.p, 02.70.-c, 02.30.Tb, 33.80.-b

come

Global integration of the Schrödinger equation: a short iterative scheme within the wave operator formalism using discrete Fourier transforms

Georges Jolicard, Arnaud Leclerc

Introduction

The numerical solution of the Schrödinger equation i ∂Ψ/∂t = HΨ plays a key role in the understanding of molecular dynamics processes. Molecular inelastic [1] and reactive collisions [2] have long been treated using quantal wave packet methods. In procedures such as the semiclassical treatment of multiphoton processes [3] explicitly time-dependent Hamiltonians have to be used to describe the molecule-field interaction. A popular and efficient wave-packet propagation algorithm is the Multi-Configuration Time-Dependent Hartree method, which is able to handle a large variety of high-dimensional problems [4]. Here, our aim is not to deal with many dimensions but to develop an efficient way of treating problems with rapidly oscillating time dependencies in the Hamiltonian operator or in the wavefunctions. Several algorithms which are used for time-resolved experiments with time-independent Hamiltonians can be adapted to the time-dependent case and are described in [START_REF] Leforestier | [END_REF] and references therein. However for such methods the length of the integration steps must be reduced to handle the high frequency terms. The propagation scheme is then based on the decomposition of the evolution operator using small time increments.

U (t, 0) = N -1 n=0 U ((n + 1)∆t, n∆t) (1)
where ∆t = t/N and

U (t + ∆t, t) ≃ exp[-(i/ )H(t + ∆t/2) ∆t] (2)
These step by step integration schemes will be referred to in later sections as "continuous methods". Second order differencing schemes (SOD) [6], split operator methods [7] and the short iterative Lanczos propagation [8] all give cumulative propagation errors which are proportional to (∆t) 3 in (1) and also errors due to the approximate calculation of the action of exp[-(i/ )H∆t] on the wave function. Other high-accuracy integrators which can be used for multi-dimensional systems include the generalization of the sympletic partitioned Runge-Kutta method developed by Sanz-Serna and Portillo [9] and the method proposed by Kormann [10] in which the Hamiltonian H is replaced by a suitable truncation of the Magnus series.

Unfortunately, these continuous integrators suffer from limitations which could prevent their use in some cases. SOD cannot handle non-Hermitian Hamiltonians, while the split operator scheme requires that the kinetic operator does not mix coordinates and their associated momenta. The calculation of the Magnus development is only tractable if the couplings have separate time and coordinate dependencies [10].

A recently proposed global treatment, the Constrained Adiabatic Trajectory Method (CATM), [11,12,13,14] introduces the Floquet Hamiltonian H F = Hi ∂/∂t defined in an extended Hilbert space including the time coordinate. The dynamical Schrödinger equation is transformed into an equivalent partial eigenvalue problem in the extended space and the wave function is forced to conform to consistent boundary conditions over a short artificial time extension by using a time-dependent absorbing potential. In [14], the CATM was formulated as a global integrator of the Schrödinger equation which does not belong to the category of methods described by (1) and (2) and it was concluded that the CATM is well suited to the description of systems driven by Hamiltonians with explicit and complicated time variations. The method does not suffer from cumulative errors (because of its global character) and the only error sources are the non-completeness of the finite molecular and temporal basis sets used and the imperfection of the time-dependent absorbing potential. The dynamics of a quantum system can also be obtained in a global way (as in the CATM) by using the time-dependent wave operator concept (TDWO) [15]. In this formalism the dynamics is formally separated into a simple evolution within a given active subspace driven by an effective low-dimensional Hamiltonian plus a secondary evolution from the subspace to the complete space. The TDWO concept has also been generalized [16] to treat almost adiabatic quantum dynamics. This generalization is based on a time-dependent adiabatic deformation of the active space and is useful for dynamics which do not escape too far from an adiabatic subspace. However the effectiveness of both TDWO and CATM hinges on the capacity to integrate some nonlinear differential equations and this is the central problem treated in this work. A solution of this problem is essential to preserve the principal advantage of these formalisms, i.e. their capacity to give a global solution over the whole interaction time. To reach this goal an iterative solution seems appropriate. Such a solution was proposed in ref. [17], using a time-dependent version of the Recursive Distorted Wave Approximation (RDWA) [15] together with non-linear transformations such as Padé approximants [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF]. Nevertheless, this iterative solution is not fully satisfactory and often fails because of a small radius of convergence and a high sensitivity to the choice of the initial guessed solution.

The present paper proposes a new iterative solution of the time-dependent wave operator equations. This is equivalent to solving the time-dependent Schrödinger equation without using the approximation of ( 1) and ( 2) at any stage of the calculation. In section 2 the basic wave operator equations are presented together with the iterative integration procedure. In practice we limit ourselves to the non-degenerate case (one-dimensional active subspace) which describes the dynamics issuing from a given quantum state. This is equivalent to solving the Schrödinger equation for the wavefunction. The iteration procedure involves calculating numerous integrals over time of matrix and vector elements. Section 3 explores the possibility of using Fast Fourier Transforms (FFT) to compute these integrals. In section 4 we illustrate the algorithm by studying the transition of a vibrational wave packet between the two wells of a model potential energy surface under the influence of two laser fields. Various numerical features are analysed in details. Section 5 gives a discussion and conclusion.

Iterative calculation of the time dependent wave operator

Wave operator equations

Let H be the Hilbert space associated with a molecular system and let S o be a model subspace of rank m. It is assumed that the initial molecular state is included within S o . The orthogonal projector corresponding to the model space is called P o , with

P 2 o = P o , P † o = P o , tr(P o ) = m.
During the evolution, the model space is continuously transformed into successive active spaces S(t) whose projectors P (t) are solutions of the Schrödinger-von Neumann equation:

i ∂P (t) ∂t = [H(t), P (t)]. (3)
The wave operator associated with the two subspaces S o and S(t) is defined as [15] Ω(t) = P (t

)(P o P (t)P o ) -1 = U (t, 0; H)(P o U (t, 0; H)P o ) -1 (4)
where U represents the quantum evolution operator associated with the Hamiltonian H(t), i.e. i ∂U(t,0;H) ∂t = H(t)U (t, 0; H). In (4), (P o P P o ) -1 is the inverse of P within S o . It exists only if the Fubini-Study distance between P o and P is small: dist F S (P o , P ) ≤ π 2 [16], i.e. if the dynamics does not escape too far from the initial subspace. In this framework the time evolution can be written as

U (t, 0; H)P o = Ω(t)U (t, 0; H ef f ) H ef f (t) = P o H(t)Ω(t) (5)
In equation 5 the evolution issuing from the initial state is separated into two terms. The first one, U (t, 0; H ef f ), describes the dynamics within the model space S o . The second one can be written as Ω(t) = P o + X(t) and it induces transitions from the model space to the full space, the off diagonal part (X = Q o ΩP o ) inducing transitions to the complementary space exclusively (Q o being the projector associated with the complementary space). It is evident that the factorisation of U (H ef f ), which includes all the fast variations within S o , is an advantage for the dynamical integration process if the model space is well chosen. It can be proven that the reduced wave operator

X(t) = Q o ΩP o satisfies [15] i ∂X(t) ∂t = Q o (1 -X(t))H(t)(1 + X(t))P o . (6)
Solving this equation is the central point of the wave operator formalism. The solution gives Ω(t) as well as H ef f (t) (cf equation 5), from which a partial evolution operator can be obtained.

A continuous integration of this non-linear differential equation has been proposed in section VI of the review [15]. Unfortunately, many approximations were necessary to make the calculation tractable and the final result was not satisfactory in the strong coupling regime. To derive a new global solution, we can use the fact that

Q o (1 -X) ∂ ∂t (1 + X)P o = Q o ∂ ∂t XP o to rewrite (6) as Q o (1 -X(t))H F (t)(1 + X(t))P o = 0 ( 7
)
where H F (t) is the Floquet Hamiltonian,

H F (t) = H(t) -i ∂ ∂t . (8)
We now derive an iterative treatment, assuming that eq (7) is not perfectly satisfied at a finite iteration order n,

Q o (1 -X (n) (t))H F (t)(1 + X (n) (t))P o = ∆ (n) (t) = 0. (9)
We are looking for the increment δX (n) which exactly solves the problem:

Q o (1 -X (n) (t) -δX (n) (t))H F (t)(1 + X (n) (t) + δX (n) (t)) = 0. ( 10
)
Expanding (10) gives four terms of different significance,

∆ (n) -Q o δX (n) (t)H F (t)δX (n) (t)P o -Q o δX (n) (t)H F (t)(1 + X (n) (t))P o +Q o (1 -X (n) (t))H F (t)δX (n) (t)P o = 0 (11)
In (11), the first term is the error at the previous iteration defined in (9). The second term is quadratic in the increment (δX (n) (t)) and can be neglected; this is an approximation commonly used in such iterative treatments. The third term can be rewritten as

Q o δX (n) H F (1 + X (n) )P o = Q o δX (n) P o H(P o + Q o X (n) P o ) -P o i ∂ ∂t (P o + Q o X (n) P o ) = Q o δX (n) P o H(P o + Q o X (n) P o ) = Q o δX (n) (t)H (n) ef f (t)P o , (12)
with

H (n) ef f (t) = P o H(t)(P o + X (n) (t)). ( 13
)
The last term in (11) needs more detailed attention:

Q o (1 -X (n) (t))H F δX (n) (t)P o = Q o H δX (n) P o -Q o X (n) H δX (n) i ∂ ∂t Q o X (n) P o (14)
We cannot handle differential equations in which (HδX (n) ) or (X (n) H) couple the unknown components of δX (n) to all the complementary space spanned by Q o . We shall thus assume that the non-diagonal elements of H in the complementary space, given by Q o (H -H diag )Q o , are small and make two last approximations,

Q o H(t)Q o ≈ Q o [H(t)] diag Q o Q o X (n) (t)H(t)Q o ≈ Q o X (n) (t)H(t) diag Q o ( 15
)
This simplification is important but does not have important consequences if the procedure still converges. The neglected terms are gradually taken into account during the iterative process and finally equation ( 7) is satisfied. In summary, the increment defined in (10) can be computed by solving 16) with H ef f defined in (13) and

i ∂ ∂t δX (n) (t) = ∆ (n) (t) -δX (n) (t)H (n) ef f (t) + H(n) diag (t) δX (n) (t) (
H(n) diag (t) = Q o H(t) -X (n) (t)H(t) diag Q o ( 17
)
In the following we limit ourselves to the case of a one-dimensional model subspace, with P o = |i i| and Q o = k =i |k k|, |i being the initial state of the system. X (n) and δX (n) becomes vectors and H ef f (t) is a pure scalar function. Equation ( 16) becomes an ordinary differential equation whose rigorous solution for the component

k|δX (n) |i = δX (n) k
can be expressed using exponential evolution operators,

δX (n) k (t) = exp 1 i t 0 H(n) diag,kk (t ′ ) -H (n) ef f (t ′ ) dt ′ (18)
× t 0 ∆ (n) k (t ′ ) exp - 1 i t ′ 0 H(n) diag,kk (t ′′ ) -H (n) ef f (t ′′ ) dt ′′ dt ′ .
This solution is consistent with the initial value of the wave operator

Ω(t = 0) = P o , (19)
hence X (n) (0) = 0. Equations ( 13), (17) and [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF] are the central equations of this paper which are to be solved numerically.

Periodic basis set and constraints on the boundaries

Similar equations to those given above were derived in [START_REF] Jolicard | [END_REF] by expanding the evolution operator using successive interaction pictures. However they were not applied numerically because of the difficulty of calculating the many time integrals which are similar to those present in [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF]. Only adiabatic limits could be calculated, with very poor results. Here we show that it is possible to design an efficient implementation by observing that the global approach represented by equations ( 7) and ( 18) is consistent with the use of a discrete Fourier basis set to span the whole propagation time interval, together with the associated periodic quadrature. These tools give effective ways to express the matrices H(t)X (n) (t) and ∂X (n) (t)/∂t present in ∆ (n) (t) (equation ( 9)) and to calculate the integrals in ( 18), as will be explained in section 3. However the accurate use of Fast Fourier Transforms (FFT) is only possible if all the timedependent functions are made perodic and smooth.

The Hamiltonian H(t) of a molecule submitted to a laser pulse is the sum of a time-independent part H o which represents the unperturbed molecule and of a perturbation V (t) with V (0) = V (T o ) = 0, where T o is the laser pulse duration. In equation ( 9) this corresponds to X (n) (t = 0) = 0 and thus ∆ (n) (t = 0) = 0. However at the end of the pulse duration, X (n) (t = T o ) can differ considerably from its initial value. To impose the continuity and the initial conditions of X (n) (t) one should use the CATM procedure [11,12]. To this end, a time-dependent absorbing potential is introduced over an artificial time extension [T o , T ]. In the present case, this imaginary potential V abs (t) destroys all the components of the wave operator (and the wavefunction), except the one corresponding to the initial state i|Ω(t) ,

V abs (t) = -iV opt (t) k =i |k k| (20)
where V opt (t) is a real positive function localised on the interval [T o , T ]. This supplementary term produces the boundary value X (n) (T ) → 0 and hence the periodicity of δX (n) and ∆ (n) (equation ( 9)):

δX (n) (t = 0) = δX (n) (t = T ) = 0, ∆ (n) (t = 0) = ∆ (n) (t = T ) = 0. ( 21
)

The problem of nested integrals over time

The n th order iteration requires solving the equations ( 13), ( 17) and ( 18) and the passage to the (n + 1) th order, via the incrementation of X using X (n+1) = X (n) + δX (n) and the construction of ∆ (n+1) (cf equation ( 9)). We want to use the same discretization at regularly distributed grid points to describe all the time dependencies. Suppose that the molecular basis set is composed of N v orthonormal states with

P o = |i i| and Q o = k =i |k k| and that the discrete time-grid which spans the time-interval [0, T ] is composed of N discrete time values t j = jT /N, j = 0, . . . , N -1. The calculation of δX (n) (t) in (18) calls for the calculation of N v -1 integrals tj 0 H(n) diag,kk -H (n) ef f dt ′ ,
where k takes all the N v -1 index values in the complementary space, followed by

N v -1 integrals tj 0 ∆ (n) j (t ′ )e -1 i t ′ 0 H(n) diag,kk -H (n) ef f dt ′′
dt ′ . This is obviously a significant task, which must be undertaken at each interation order (n). The calculation of each integral must be fast and very accurate. The selected procedure is explained in the next section with a preliminary discussion about integrating an arbitrary complex function of time.

Discrete Fourier integration at the n th iteration order

In this section we explain how to compute efficiently a numerical definite integral of a time-dependent function described on regularly distributed sampling points. The function is not necessarily continuous at the boundaries of the grid. This preliminary discussion is important in facilitating the numerical implementation of equation [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF]. We require that the procedure should determine with the same accuracy, not only the integral at the final time value t = T (which is easily obtained using the standard Fourier quadrature rule) but also its value for each of the N intermediate discrete time values t j .

Continuous integrals and numerical FFT-integrators

In the present subsection the Fourier transform F (ν) is defined as

F (ν) = +∞ -∞ f (t) exp(-i2πνt)dt. ( 22
)
Let I(t) be the integral

I(t) = t -∞ f (t ′ )dt ′ = +∞ -∞ f (t ′ )h(t -t ′ )dt ′ (23)
where h(t) is the Heaviside step function,

h(t) = 0 t < 0 1 t ≥ 0.
The Fourier transform of h(t) is the distribution

H(ν) = 1 2 δ(ν) -P V i 2πν . ( 24
)
PV denotes the Cauchy principal value which becomes relevant when the above expression is integrated. The convolution theorem gives

I(t) = f * h(t) = +∞ -∞ F (ν)H(ν) exp(2iπν)dν. ( 25
)
Using H(ν) from equation ( 24) then gives

I(t) = 1 2 F (0) - i 2π P V +∞ -∞ F (ν)e 2iπνt ν dν. (26)
In molecule-field interaction problems, the field term is frequently an oscillating function of time. Let us consider for a moment the special case of an oscillating function of the form f (t) = g(t) exp(i2πνt). Its Fourier transform is the translated function

F (ν) = G(ν -ν),
where G(ν) is the transform of g(t). Using this F (ν) in ( 26) and using the variable change ν ′ = ν -ν, we obtain

I(t) = t -∞ g(t) exp(i2πνt)dt = 1 2 G(-ν) - i 2π e i2π νt P V +∞ -∞ G(-ν ′ )e -2iπν ′ t ν -ν ′ dν ′ . ( 27
)
When the frequency ν is very high, we can neglect the Fourier coefficient G(-ν) and approximate the denominator ν -ν ′ by ν, giving the approximation

t -∞ g(t ′ ) exp(2iπνt ′ )dt ′ ∼ ν→∞ exp[i(2πνt -π/2)] g(t) 2πν . ( 28
)
We can now go back to the general case of equation (26) to derive a numerical FFT-integrator. The detailed derivation is given in Appendix A. This appendix explains how equation ( 26) can be implemented for periodic or non-periodic functions known as N sample values distributed over the interval [0, T ]. The integration algorithm needs only two FFT, with a computational cost scaling as 2N log N to obtain the N values of I(t j ), with t j = jT /N, j = 0, 1, . . . , N -1. Here we recall only the main result,

I(t j ) = F F T -1 j (µ.F F T (f )) + j N × T √ N F F T k=0 (f ). ( 29
)
where µ.F F T (f ) is a component by component vector multiplication, with

µ ℓ = i 2π T ℓ ℓ = 1, . . . N 2 -1 , i 2π T ℓ-N ℓ = N 2 , . . . N -1. (30)
In addition, the ℓ = 0 component of (µ.F F T (f )) in ( 29) must be replaced by the following quantity:

a = - N -1 ℓ=1 µ ℓ .F F T ℓ (f ). ( 31
)
In the next subsection we perform some numerical tests using the algorithm derived in appendix Appendix A with rapidly oscillating functions.

Numerical test of the FFT-based integration algorithm

We study the arbitrary case of a function composed of six gaussians multiplied by complex exponential functions. The gaussians have various widths and are distributed on the interval [0, T = 180] to give the function

f (t) = 6 j=1 exp[-a j (t -t j ) 2 ] exp(i2πx j t/T ). ( 32
)
The various parameters appearing in (32) are given in table 1. The a j are chosen so that the modulus of the function f (t) is equal to or smaller than 10 -20 at the two boundaries (see figure 1). In this example, the FFT procedure uses the function f (t) defined by (A.6). This modified function is characterized by a spectrum with a central component ℓ = 0 equal to zero. This can be seen in figure 2, which shows the modulus of the spectrum of f (t). As the functions converge rapidly to zero at t = 0 and at t = T , the continuity of the function f (t) (assumed to be a periodic function by the FFT procedure) and also the continuity of the first two derivatives are ensured without the need to extend the time interval or to use an asymptotic polynomial (eq (A.10)). Figure 2 confirms that the functions are well represented by the Fourier expansion; the components omitted in the spectrum are smaller than 10 -13 . We first compare the speed of convergence of two algorithms, the numerical FFTintegrator (this work) and the standard Simpson rule. As the exact values of the integrals I j = I(t j ) are not known, we perform two series of calculations by doubling the number N of sampling points, comparing for each of the two integration procedures the result for 2N points with that for N points. The convergence factor CF N is calculated by using the N points (with even indices) that the 2N grid has in common with the N grid. At each of the N sampling points the absolute value of the difference between the two results for the integral is calculated and the convergence factor is defined as their maximum value, namely

CF N = N -1 max j=0 (|I (2N ) 2j -I (N ) j |). (33)
The corresponding numerical results are shown in figure 3. As the number of sampling points increases, the FFT convergence factor decreases rapidly by many orders of magnitude up to a critical N value, N cri F F T = 1024. Beyond this point the convergence factor reaches a constant plateau with a precision close to machine precision (10 -14 ). This is the usual behaviour expected from the Nyquist-Shannon theorem [26] when using Fourier basis expansions. The Simpson rule results show a strongly different behaviour. The convergence factor decreases regularly up to a minimum value for a critical value N cri Simp , which is about 100 times larger than for the FFT case. Moreover the minimum convergence factor reached is about two orders of magnitude greater than for the FFT results. After this minimum value, the convergence factor increases again. In this new regime the roundoff errors become predominant. Since no exact result is available, we next compare the results obtained with the best precision for each method, to be sure that the two algorithms are consistent and converge to the same value. This corresponds to N = 1024 for the FFT algorithm and to about N = 150000 for the Simpson algorithm. Figure 4 presents for each discrete time value the modulus of the difference between the two integrals calculated with the FFT algorithm and with the Simpson rule,

|I F F T (t j ) -I Simp (t j )|. ( 34
)
This figure reveals that the difference on the whole interval [0, T ] is never larger than 10 -13 . Similar results, not presented here, are obtained for examples in which it was necessary to impose a continuity of the two first derivatives by using a 6 th order polynomial (see Appendix A.2). These applications prove that the FFT integrator possesses the expected qualities. Figure 4 shows that the precision of the results is very high and virtually independent of the discrete instant t j taken as the upper bound. Moreover, figure 3 reveals that the number of sampling points necessary to give convergence for the FFT is 100 times less than that for the convergence of the Simpson algorithm. Finally the operations necessary for the FFT algorithm are reduced to the calculation of two FFT, the first being made in any case during the propagation scheme [START_REF] Kosloff | Quantum Molecular Dynamics on Grids chapter V Dynamics of molecules and chemical reactions[END_REF]. For the selected example the CPU time associated with the FFT procedure is about 20 times smaller than that associated with the Simpson rule. However, the degree of advantage of the FFT approach is closely related to the shape of the spectrum of f (t). In the present cases as well as for the Schrödinger equation integration procedure treated in the next section, the spectrums are narrow and are thus well represented by a Fourier basis set.

An asymmetric double-well transition experiment

The illustrative example is that of a model diatomic molecule submitted to two laser pulses. Throughout this section we use arbitrary units with the convention = 1.

The various numerical parameters have been adjusted so that the relative orders of magnitude of the vibrational level spacing, the strength of the dipole couplings and the pulse duration produce realistic dynamics with significant transition probabilities.

Two potential curves S = 1, 2 are defined, which can refer to the two first electronic states of a one-dimensional vibrational Hamiltonian in the framework of the Born-Oppenheimer approximation. These two curves are shown in figure 5. The lower surface is a double well similar to those used in references [START_REF] Bavli | [END_REF]29],

ǫ 1 (R) = -5R 2 + 0.5R 3 + R 4 . ( 35
)
The upper surface is a single quartic well,

ǫ 2 (R) = 0.2R 4 . ( 36
)
The R coordinate is chosen as a relative position centred on the barrier of ǫ 1 . All the eigenstates are pure bound states and were calculated using a Fourier grid method on the radial coordinate [30,31]. They are shown in figure 5. These model potentials are also similar to those used to describe, for instance, effective isomerization problems along a reaction coordinate [32]. The evolution of the vibrational evolution operator

(R) = -5R 2 + 0.5R 3 + R 4 ; ǫ 2 (R) = 0.2R 4 .
is driven by the following equation within the framework of the dipole approximation:

i ∂ ∂t U (t, 0) = T N + ǫ 1 -µ 1,2 . E(t) -µ 1,2 . E(t) ǫ 2 U (t, 0) (37)
T N is the relative kinetic energy of the two atoms, i.e.

T N = -1 2m ∂ 2
∂R 2 with an arbitrary mass equal to 10. ǫ S=1,2 are the two energy surfaces. There is no rotation in the model and the dipole transition moment µ is assumed to be constant along the radial axis R and equal to 1.

The laser field amplitude projected on the µ direction is defined as the sum of two pulses with gaussian envelopes

E(t) = 2 j=1 E j cos (ω j (t -T j )) exp - t -T j τ j 2 (38)
Vibrational states (v = 1, S = 1) and (v = 6, S = 1) have in common a strong overlap with state (v = 7, S = 2) of the excited potential curve. We have therefore chosen the carrier frequencies as ω 1 = E(v = 7, S = 2) -E(v = 1, S = 1) ≃ 9.98449 and ω 2 = E(v = 7, S = 2) -E(v = 6, S = 1) ≃ 4.77725. This choice induces resonant transitions between the first and the sixth eigenstate of the lower surface:

ω 1 -ω 2 = E(v = 6, S = 1) -E(v = 1, S = 1
). This will ensure strong transitions between the two asymmetric wells and a rather large depletion of the initial state. Calculations were made using τ 1 = 3.90 and τ 2 = 4.50 and two different electric field amplitudes and centres,

E 1 = 0.05, T 1 = 23.5 E 2 = 0.08, T 2 = 21.5 (39)
in the first example and

E 1 = 0.09, T 1 = 22.5 E 2 = 0.05, T 2 = 21.5 (40)
in the second example. The corresponding electric fields are represented in figure 6.

Those choices for the laser parameters place the system in a non optimal STIRAP configuration [33]. In both cases the initial state is the first eigenstate of the first surface (v = 1, S = 1). In terms of the wave operator, the first global iteration starts with a guess chosen as Ω (0 (i.e. a total molecular basis set of N v = 60 states) is sufficient to obtain converged results. This molecular basis set is small but here the principal aim is to test the time propagation algorithm. Here we need N t = 4096 Fourier states |p with t|p = exp(i2πpt/T ) (or equivalently N t = 4096 time grid points) to describe the time evolution correctly. In what follows the vibrational states are numbered from v = 1 to v = 30 for the first surface and from v = 31 to v = 60 for the second surface. The total basis set describing the extended Hilbert space is then composed of N = N v × N t = 245760 states. In this framework, the dynamics is calculated by using the iterative scheme explained in section 2 (equations ( 5), ( 13), ( 17) and [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF] with Ω (n+1) = P o + X (n+1) and X (n+1) = X (n) + δX (n) .) The integrals appearing in (18) are calculated using the numerical integrator presented in section 3 and Appendix A. This scheme is used to calculate all the components of the off-diagonal part of the wave operator for each of the N t discrete t k values:

) (t) = |v = 1, S = 1 v = 1, S = 1|, i.e.
X (n) v (t k ), v = 1, . . . , N v , k = 1, .
. . , N t , at each iteration order (n). Our main concern is the speed of convergence of the global iterative scheme. We define the following convergence factor:

F = Nv v=1 Nt k=1 |δX (n) v (t k )| 2 Nv v=1 Nt k=1 |X (n+1) v (t k )| 2 , (41)
with δX (n) defined in [START_REF] Killingbeck | Microcomputer Algorithms Action from Algebra[END_REF]. Figure 7 shows that the two calculations converge. A plateau of about 10 -17 for the first example and of about 10 -16 for the second one is reached at the 22 th iteration. Between the 10 th and the 22 th iteration, the accuracy is improved by one order of magnitude at each iteration. The convergence factor of equation ( 41) adds the residues δX (n) (t) over the whole propagation interval. This means that the plateau seen in figure 7 indicate a very high accuracy at each time grid point. It is essential to note that, by contrast with a continuous propagation scheme, in which errors are accumulated at each time step, the global scheme does not accumulate the errors between two iteration orders. Nevertheless our iterative procedure possesses, as does every iterative treatment, a finite radius of convergence. The small increase of the convergence factor observed in figure 7 during the first four iterations for the second example is the consequence of a large Fubini-Study distance between P o and P in this case (the survival probability at the end of the interaction | v = 1|Ψ(t f inal ) | 2 is 0.2129 in the first example but only 0.02864 in the second one). This increases the distance between P o and P (t f inal ) in this last case. This can be seen in figure 8 which illustrates the progressive separation between the active subspace S o and the dynamic subspace S(t). In the non-degenerate case, the Fubini-Study distance is simply

dist F S (P o , P (t)) = arccos ( v i |Ψ(t) ) = arccos 1 Ω(t) . ( 42
)
The right-hand side is obtained using the normalization properties of the wave operator. The final distance for the second example of ( 40) is clearly larger. The physical reason explaining this strong difference is that example (40) is closer to an efficient STIRAP situation than example (39) (but still non-optimal, the final transition probability to state |v = 6 being here only 0.22, because of many nonadiabatic transitions to other vibrational states).

The global character of the integration is evident in figures 9 and 10 which present the transition moduli to the different vibrational states at the final time | v|Ψ (n) (t f inal ) | at four different iteration orders. It is interesting to note that in figure 9, an approximate vibrational distribution close to the exact one is obtained as early as at the fourth iteration. The small convergence difficulties during the first iteration orders in the second example are evident in figure 10. The 4 th order shows some transition amplitude larger than 50, very far from the final results. In this case more iterations are needed to achieve a stabilized distribution. Probabilities can thus exceed one in the global scheme if the results are observed at an early intermediate order, before the completion of the iterative process, such as in Fig. 10. There is no fixed norm conservation in the structure of the algorithm but the norm conservation is recovered after convergence. Figures 11 and12 also illustrate the global character of the time-dependent propagation. In these figures we present the evolution of the transition amplitude to state |v = 37 . This particular state is strongly coupled with the initial state |v = 1 .

Figure 11 shows that each iteration order produces a projected amplitude which converges towards the exact solution over the whole interaction time interval. At the fourth iteration, the difference

| v = 37|Ψ (n=22) (t) -Ψ (n) (t)
| is already lower than 10 -1 . At the 9 th order, the same error becomes lower than 10 -3 . Figure 12 reveals the difficulties due to the large Fubini-Study distance between P o and P (t f inal ) (here the final survival probability is | v = 1|Ψ(t f inal | 2 = 0.02864). As a consequence, the amplitude of the projected wave function onto the state v = 37 at the 4 th iteration order is about 10 times larger than the exact one for t > 30. Surprisingly, the iterative procedure converges again between the 4 th and the 9 th order. At the 9 th order the error is smaller than 10 everywhere. We think that this efficiency is essentially due to the use of a time-dependent effective Hamiltonian in the basic equations of the integration scheme (( 13), ( 17) and ( 18)), especially when it is compared with previous iterative schemes. Figure 13 reveals large oscillation of the effective Hamiltonian defined in [START_REF] Leforestier | [END_REF]. The imaginary part imposes the unitarity of the evolution in the wave operator formulation. A contribution to the evolution operator of the form exp 1 i is present and gives the correct survival amplitude. The magnitude of these oscillations (≈ 4 for the real part) is even much larger than the distance between the first closest vibrational eigenvalues of each surface (for example

t f inal 0 ℑm(H ef f )dt ′
| v = 37|Ψ (n=22) (t) -Ψ (n) (t) | for (n = 4) and (n = 9).
E(v = 2) -E(v = 1) = 1.451).
Other iterative methods in which these time variations are neglected completely fail to converge on this same numerical example. For example, the adiabatic approximation would be justified if in equation ( 18) the amplitude | v|

H(n) diag (t)|v -v i |H (n)
ef f (t)|v i | is large at each instant, |v i being the initial vibrational state. In such a case the introduction of the approximation (28) into equation ( 18) leads to:

δX (n) (t) = i ∆ (n) (t) H(n) diag (t) -H (n) ef f (t) (43)
Equation ( 43) is very similar to the RDWA iterative method which has already been used to compute time-dependent wavepacket propagations [15,14]. However it cannot be used in the present case. It is evident that the strong time-dependency of H ef f (t) (see figure 13) and the very small survival probability at t f inal are inconsistent with the use of any adiabatic approximation. Actually the use of (43) in the integration scheme produces a very fast divergence of the iterative series caused by the crossings of H(n) diag and H

(n) ef f . Another approximate (but more sophisticated) integration scheme was proposed in refs [17] and [15]. By using a discrete Fourier basis set: |p with t|p = exp(i2πpt/T ) derived from the discrete time grid basis set by a FFT transformation, a generalised RDWA approximation leads to:

v ′ , p ′ |X (n+1) |v i , p = 0 = v ′ , p ′ | (H F -H(n) )X (n) + H |v i , p = 0 v i , p = 0|H (n) ef f |v i , p = 0 -v ′ , p ′ | H(n) |v ′ , p ′ . ( 44
)
Introducing (44) into the iterative scheme also produces a divergence, albeit slower than the one produced by (43). Undoubtedly this defect comes from the fact that equation (44) only uses an averaged H ef f , represented by the first Brillouin component p = 0|H ef f |p = 0 , rather than the exact and strongly time-dependent expression.

Finally it is crucial to verify that the global scheme gives the correct results by a comparison with some continuous wave-packet propagation algorithms. We first compare our results with those obtained for the same system by using a split operator [7] procedure by separating the role of the diagonal unperturbed Hamiltonian and the role of the laser-molecule coupling. This last contribution is not diagonal and its role is approximated by using a Second Order Differencing (SOD) scheme. We also calculate the same dynamics using the short iterative Lanczos (SIL) method [8]. To be significant, the same vibrational basis set of N v = 60 states is used in the continuous propagations (for the SIL propagations, the size of the Lanczos subspace was set up to 10). A convergence factor is defined as

F C (N t ) = Nv i=1 |Ω C,Nt (i, t f inal ) -Ω G,4096 (i, t f inal )| 2 . (45)
Ω G,4096 is the wave operator calculated at the final time using the global propagator with 4096 time grid points after 22 iterations and Ω C is the wave operator reconstructed by using equation ( 4) after a complete propagation with the continuous propapagators with N time steps. To show the convergence of the global method itself we have also computed a similar convergence factor between global results,

F G (N t ) = Nv i=1 |Ω G,Nt (i, t f inal ) -Ω G,8192 (i, t f inal )| 2 . ( 46
)
Figure 14 shows the evolution of F C and F G with respect to the number of discrete time steps. This figure reveals that the three procedures converge to the same final values (same amplitudes and same phases). Figure 14 also confirms that the number of time steps needed to converge to a given accuracy is much smaller with the global scheme than with the continuous schemes. The convergence is faster for the first pulse than for the second one. The number of time-steps of length dt for the continuous integrator is quite large especially if a high accuracy is needed. The CPU time needed to compute δX (n) (equation ( 18)) at a given iteration order is of course much larger than that required to construct one step in the continuous algorithms. Nevertheless the convergence plateau in figure 7 (39). The corresponding number of significant digits on the largest transition probabilities are also given.

if a low accuracy is sufficient. The SIL algorithm and the global method then give comparable CPU times for intermediate accuracies, while the split operator stays a little faster. The global method becomes slightly faster if a high accuracy is desired.

Conclusion

The theoretical analysis and the simple illustrative examples presented in this paper demonstrate that a short global iterative scheme can solve the time-dependent Schrödinger equation with an explicitly time-dependent Hamiltonian operator. This global method is competitive with a continuous step by step algorithm. The following points should be highlighted: i) A global integrator does not accumulate discretization errors along the time interval and consequently can generate high precision results. The numerical experiments have shown that each new iteration order increases the precision by one order of magnitude. The structure of the algorithm also allows for a rigourous estimation of the error associated with each iteration order by following the quantity H F Ω -ΩHΩ, the calculation being stopped when the expected accuracy is reached.

ii) The use of an FFT procedure to calculate the numerous time integrals appears as a very efficient choice which is compatible with producing a highly accurate solution to the Schrödinger equation. Nevertheless it is essential that the Fourier spectrum of the perturbation should be well reproduced by the Fourier basis set and it is certainly more difficult to solve problems with broader spectra.

iii) The global iterative scheme does not introduce limitations concerning the nature or the structure of the Hamiltonian which drives the dynamics. Future applications could usefully be aimed at checking how well it can reproduce dissipative dynamics.

iv) The inability of previous global iterative techniques to determine the timedependent wave-operator for similar problems is partly due to a wrong estimation of the effective Hamiltonian acting in the selected model space. In the selected examples H ef f (t) exhibits large oscillations which must correctly be taken into account.

v) In the present applications the starting input for the wave-operator is simply X (n=1) = 0. The iterative nature of the method makes possible a rapid generation of new solutions, after a slight modification of the Hamiltonian parameters, by simply starting with the wave-operator obtained from the previous solution. Many applications are possible, among which we can mention the study of laser control processes using exceptional points of parametric Hamiltonians [34,35,36,37] or of vibrational cooling using zero-width resonances [38].

The algorithm proposed in this paper appears to be useful for investigating new formulations of quantum control in molecular physics or for treating complicated nearadiabatic evolutions, for which the use of an effective Hamiltonian can efficiently reduce the dimension of the Hilbert space [39]. The global algorithm add the time as a new quantum coordinate, which can be a disadvantage to compute the dynamics of high dimensional problems over long times. On the other hand the method is not dependent on the structure of the Hamiltonian as it can be the case for the split operator algorithm or Magnus expansions. The wave operator approach is in principle transposable to dynamical processes like dissociation, ionization, dynamics involving degenerate eigenvalues or non-adiabatic interactions. For dissipative processes, a special difficulty can appear in situations for which the dynamics escapes too far from the selected model subspace. In such cases, the algorithm can diverge. The best option is to generalise the iterative algorithm of section 2 to handle a degenerate active space, by incorporating all the states which are strongly coupled into the model space. Study of this generalisation is in progress, as the adaptation of the algorithm to the calculation of Floquet eigenstates. where

F F T k (f ) = 1 √ N N -1 j=0
f j e -i2πjk/N (A.4) with f j = f (t j ). The integral in equation ( 26) is then apparently approximated by the sum

I(t j ) = 1 2 T √ N F F T k=0 (f ) - i 2π 1 √ N N -1 ℓ=0 F F T ℓ (f ) ν ℓ e 2iπℓj/N . (A.5)
However, the passage from the continuous expression (26) to the FFT equivalent (equation (A.5)) has not been carried out with sufficient rigour, since there is a division by 0 for ℓ = 0 in (A.5). We must adjust the term ℓ = 0 before doing the inverse transform. The use of a Cauchy principal value and the presence of a singularity in (26) at ν = 0 if F (ν = 0) = 0 are inconsistent with the use of a Fourier series expansion. The equations should thus be modified to take into account this discrepancy. However, we can preserve the feature that only two Fourier transforms are required to obtain the integral. The calculation of the N integrals I(t j ), j = 0, . . . , N -1 is made as follows. The Fourier transform F F T (f ) is first calculated and then a new function f (t) is defined:

f (t) = f (t) - 1 √ N F F T k=0 (f ). (A.6)
Both f and f have the same discrete spectrum, except for the ℓ = 0 component, which for f is equal to zero, so that the integral of f (t) over the interval [0, T ] is equal to zero. This is equivalent to using a translation of the ordinate axis. We work temporarily with f instead of f . We assume that the spectrum F (ν) (equal to zero for ν = 0) has a Taylor expansion near ν = 0: F (ν) = aν + bν 2 + . . .

We are looking for the finite coefficient a = lim ν→0 F (ν) ν , so as to apply equation (A.5) at ℓ = 0 in a consistent way. The divergence produced by the term ν ℓ ∝ 1 ℓ for ℓ = 0 in (A.5) has disappeared but the value of a is still unknown. One can guess it by continuity at the boundaries. Let f (0), f ′ (0), f ′′ (0) and f (T 0 ), f ′ (T 0 ), f ′′ (T 0 ) be the values of the function and of the first two derivatives at the two boundaries. It will be impossible to obtain accurate results if these values do not match. One can extend the time interval by adding a narrow extra interval [T 0 , T ] spanned by an extension of the grid t j = jT 0 /N 0 , j = N 0 , . . . , N -1 and then constructing on this new interval a 6 th order polynomial: Seven conditions are imposed on the polynomial. The first three are related to the continuity of f (t) at the point t j = T 0 and give the first three coefficients a j :

f (t) =
a o = f (T 0 ), a 1 = f ′ (T 0 ), a 2 = f ′′ (T 0 )/2.
(A.11)

The three values on the right hand sides of (A.11) are known exactly if the function f (t) has an analytical expression. Alternatively, the derivatives can be approximated by a FFT algorithm if f (t) is known at a sample of N 0 values. The last four coefficients are related to the continuity between t = T and t = 0 of f (t) and to the condition

T 0 f (t)dt = 0: f (T ) = f (0), f ′ (T ) = f ′ (0), f ′′ (T ) = f ′′ (0), T T0 f (t)dt = 0.
(A.12)

With ∆ = T -T 0 , the last four unknowns are found by solving the small linear system:

   ∆ 3 ∆ 4 ∆ 5 ∆ 6 3∆ 2 4∆ 3 5∆ 4 6∆ 5 6∆ 12∆ 2 20∆ 3 30∆ 4 ∆ 4 4 ∆ 5 5 ∆ 6 6 ∆ 7 7         a 3 a 4 a 5 a 6     =     f (0) -a 0 -a 1 ∆ -a 2 ∆ 2 f ′ (0) -a 1 -2a 2 ∆ f ′′ (0) -2a 2 -a 0 ∆ -a 1 ∆ 2 2 -a 2 ∆ 3 3     (A.13)
Note that the last (integral) condition (A.12) is not absolutely necessary and can be relaxed by first defining the above polynomial and then applying the zero-integral condition on the whole interval using equation (A.6).

In summary, the FFT-based integration algorithm is implemented using (A.9), completed by a preliminary treatment on the additional interval [T 0 , T ] using (A.10)-(A.13).

Figure 3 .

 3 Figure 3. Convergence factor (33) for the FFT algorithm ( ) and the Simpson algorithm ( ).

Figure 4 .

 4 Figure 4. Difference between the FFT and the Simpson result: |I F F T (t j ) -I Simp (t j )| (full line) and the modulus of the FFT result |I F F T (t j )| (dashed line).

Figure 5 .

 5 Figure5. The two potential energy curves and some of the first unperturbed vibrational eigenstates of the first potential surface (full lines) and of the second one (dashed lines). The potentials are quartic polynomials defined as: ǫ 1 (R) = -5R 2 + 0.5R 3 + R 4 ; ǫ 2 (R) = 0.2R 4 .

Figure 6 .

 6 Figure 6. The electric field E(t) in the first example of (39) (left frame) and in the second example of (40) (right frame).

Figure 7 .

 7 Figure 7. Convergence factor (41) with the parameters defined in (39) ( ) and those defined in 40 ( ) with respect to the iteration number.

Figure 8 .

 8 Figure 8. Evolution of the Fubini-Study distance dist F S (Po, P (t)) during the pulse defined in (39) (solid line) and (40) (dashed line). The distance has been computed using the wave operator at iteration 22.

Figure 9 .Figure 10 .

 910 Figure 9. Final transition amplitudes to the different vibrational states | v|Ψ (n) (t f inal ) | after the pulse defined in (39), obtained from the wave operator at successive iteration orders. + : (n = 2), × : (n = 4), * : (n = 9), • : (n = 22).

Figure 11 .

 11 Figure 11. The dashed line represents the modulus of the projected wave function | v = 37|Ψ (n) (t) | as a function of time at the end of the iterative process (n = 22) for the first example (parameters given in (39)). The full lines are the defects at intermediate iteration orders: | v = 37|Ψ (n=22) (t) -Ψ (n) (t) | for (n = 4) and (n = 9).

Figure 12 .Figure 13 .

 1213 Figure 12. The dashed line represents the modulus of the projected wave function | v = 37|Ψ (n) (t) | as a function of time at the end of the iterative process (n = 22) for the second example (parameters given in (40)). The thin continuous lines are the defects at intermediate iteration orders : | v = 37|Ψ (n=22) (t) -Ψ (n) (t) | for (n = 4) and (n = 9). The deep black line shows the survival amplitude | v = 1|Ψ (n=22) (t) |.

Figure 14 .

 14 Figure 14. Convergence factors F C and F G defined in equations (45) and (46) for the global method (circles), the short iterative Lanczos scheme (triangles) and the split-operator algorithm (squares) using the field parameters of equation (39) ( , △ and •) and of equation (40) ( , and •), as a function of the number of time steps.

2 + 1 , 0 f

 210 j = -ν N -j , j = N . . . , N -1.(A.2)If the time t and frequencies ν are selected as one of the points of the above discrete grid, a discrete version of the integral can be computed. The continuous Fourier transform is approximated by a discrete Fast-Fourier Transform (FFT):+∞ -∞ f (t)e -2iπν k t dt = T (t)e -2iπν k t dt ≈ T √ N F F T k (f ) (A.3)

6 j=0a

 6 j (t -T 0 ) j . (A.10)

Table 1 .

 1 Parameters used in equation(32).

	Gaussian	a j	t j	x j
	index j
	1	6.790 27.000	12.000
	2	3.819 36.000	135.600
	3	1.018 90.000	1.750
	4	1.591 108.000 154.700
	5	2.118 135.000 3.250
	6	3.310 144.000 18.150

time Figure 1. |ℜef (t)|, see (32) and table 1.

Table 2 .

 2 are reached after only 22 iterations. With the field parameters of equation(39), passing from iteration 22 to iteration 23 indicates that the transition amplitudes are converged with 7 stable and significant digits. More detailed informations about the CPU time are given in table 2 with varying required accuracies. This table shows the ratio between the CPU times required to reach a given threshold of convergence using the global or the continuous strategies, using the convergence factors F G and F C defined in (45) and (46). The results of the global method have been observed after 22 iterations. The continuous algorithms are faster Convergence factor < 10 -6 10 -8 10 -10 10 -12 10 -14 Approximative ratios of the CPU times required to reach a given convergence threshold between the global (G) and the continuous methods (SIL or split-operator), for example of equation

	Stable digits ≃ T CP U G /T CP U SIL T CP U G /T CP U Split-Op.	3 3.83 3.11 0.775 0.757 0.211 4 5 6 7 15.3 6.52 3.31 1.62 0.912

Acknowledgments

We thank John P. Killingbeck for fruitful discussions and for a careful reading of the English manuscript. Simulations have been executed on computers of the Utinam Institute of the Université de Franche-Comté, supported by the Région de Franche-Comté and Institut des Sciences de l'Univers (INSU).

Appendix A. Numerical FFT integration

Appendix A.1. Integration of periodic functions Consider equation (26),

where F (ν) is the Fourier transform of f (t). If the function f (t) is assumed to be nonzero only over a finite interval [0, T ], then a discrete finite time-grid can be introduced to span this time-interval and also the corresponding finite frequency representation.

For even N we have:

requiring that I(t 0 = 0) = 0 0 f (t)dt = 0. Equation (A.5) in this case leads to the value

where µ.F F T (f ) denotes a component by component vector multiplication, with

This value of a replaces the ℓ = 0 component of (µ.F F T (f )) in equation (A.9) below. So far we have ignored the zero frequency component of the function by working with f instead of f . The final result should then be corrected by a linear correction function, giving

The first term on the right hand side represents the integral of f at the discrete time values t j and the second term restores the correct value of tj 0 f (t)dt. Equations (A.9), (A.8) and (A.7) are the basic equations needed to calculate the N integrals I(t j ) using two FFT calls (with a computational cost scaling as 2N log N ).

The above equations have close similarities with those of the FFT-based integration algorithms proposed in ref. [20,[START_REF] Yaroslavsky | Fast discrete sinc-interpolation: a gold standard for image resampling in[END_REF] in the framework of image resampling. The derivation starting from (26) emphasizes the difficulties introduced by the Cauchy principal value term for small frequencies. Here the principal value integral is approximated by first cancelling F (ν = 0) and guessing its first derivative at ν = 0 in a consistent and unambiguous way. Our expressions also work well for functions with integrals over the interval [0, T ] which are not necessarily equal to zero.

Appendix A.2. Integration of non-periodic functions

In the algorithm developed in section 2 the smoothness of the wave operator and the wavefunction at the time boundary t = T is ensured by the presence of a timedependent absorbing potential (cf. section 2.2). However the intermediate matrix elements of H (n) ef f (t) and H(n) diag,kk (t) to be integrated are not necessarily T-periodic. FFT-based differentiation and integration algorithms are very accurate but suffer from boundary effects (Gibb's phenomenon). The present algorithm will give accurate results only with periodic continuous functions, ideally those with time derivatives which are themselves also continuous functions. There are different ways of adapting a Fourier-based algorithm to deal with non-periodic signals. For example, extension methods have been developed to solve partial differential equations on a complex domain [22,23], to represent surfaces of complex bodies [24] and to treat imaging problems [20].

Here we use the alternative procedure of Hermite interpolants to ensure the continuity of the function f (t) and of the first derivatives [25]. Let f (t) be a complex continuous L 2 function defined on a large time interval [0, T 0 ], and represented by a sample of N 0 equally distributed values f (t j ), t j = jT 0 /N 0 , j = 0, . . . , N 0 -1 on this interval. The procedure described in (A.6) has been applied to f (t) so that F F T k=0 (f ) = 0. As the function f (t) is assumed to be continuous with continuous first and second derivatives inside the interval [0, T 0 ], the principal problem is to ensure