Semiclassical Analysis of Discretizations of Schr 6 dinger-type Equations

We apply Wigner-transform techniques to the analysis of difference methods for Schr6dinger-type equations in the case of a small Planck constant. In this way we are able to obtain sharp conditions on the spatial-temporal grid which guarantee convergence for average values of observables as the Planck constant tends to zero. The theory developed in this paper is not based on local and global error estimates and does not depend on whether caustics develop or not. Numerical examples are presented to help interpret the theory.

INTRODUCTION Many problems of solid-state physics require the (numerical) solution of the Schr6dinger equation in the case of a small (scaled) Planck constant c" C2 euet i--f Au + iV(x)u O, x m t 0) (x) x (1.1b) Here V is a given electrostatic potential, 0 < e << and u= u(x, t) is the wave function.By classical quantum physics [10] the wave function is an auxil- iary quantity used to compute the primary physical quantities, which are quadratic function(al)s of u(t), e.g., the position density n (x, t) ]u (x, t) [2,   (1.2a) the current density J(x,t) Im(u(x,t)Vu(x,t)), ( (where "-" denotes complex conjugation).It is well known that the Eq.(1.1a) propagates oscilla- tions of wave length e, which inhibit u from con- verging strongly in, say, L (L 2x).Clearly, the weak *Tel.: +43 (0)732-2468 9190, e-mail: markowic@math.tu-berlin.de Tel." +39 0 382 529600, e-mail: pietra@dragon.ian.pv.cnr.itCorresponding author.Tel." +43 (0) 732-2468 9187, e-mail: carsten.pohl@jk.uni-linz.ac.at convergence of u is, for example, not sufficient for passing to the limit in the macroscopic densities (1.2), which implies that the analysis of the so- called semi-classical limit is a mathematically rather complex issue.Recently, much progress has been made in this area, particularly by the introduction of tools from microlocal analysis (defect measures [5], H-mea- sures [19] and Wigner measures [6,11,13,7]).These techniques, which go far beyond classical WKBmethods, have shown the right way to exploit prop- erties of the Schr6dinger equation which, despite the oscillatory nature, allow the passage to the limit e 0 in the macroscopic densities by revealing an underlying kinetic structure.
Exactly the same problem, i.e., the highly oscillatory nature of the solutions, has to be coped with when the Schr6dinger equation with small c is solved numerically.Even for stable discretiza- tion schemes (or under mesh size restrictions which guarantee stability) the oscillations may very well pollute the solution in such a way that the quadratic macroscopic quantities (1.2) and other physical observables come out completely wrong when the scaled Planck constant is small. [3]).
However, the classical strategies cannot be employed to analyze properties of discretization schemes for e-+ 0.
In this paper we adapt the microlocal techniques used to analyze the semi-classical limit for the IVP (1.1) to the analysis of finite difference discretiza- tions.We choose a sample of time discretizations: the Crank-Nicolson scheme and the Leap-Frog scheme (both belonging to the set of mostly used methods for the Schr6dinger equation).Corresponding spatial discretizations are general arbi- trary-order symmetric schemes.For these methods we identify the semiclassieal Wigner measure (on the scale ) for all (sensible) combinations of e and of the time and space mesh sizes.We clearly have convergence for the average values of all (regular) observables in those cases, for which the Wigner measure of the numerical scheme is identical to the Wigner measure of the Schr6dinger equation it- self.Thus, from this theory we obtain sharp (i.e., necessary and sufficient) conditions on the mesh sizes which guarantee good approximation quality of all (smooth) observables for e small.We remark that the approach is sharp because we do not use local and global error estimates.
The theory and the presented test calculations clearly show the big risk in doing Schr6dinger calculations for small Planck constant.Even highly stable schemes might produce completely wrong observables under seemingly reasonable meshing strategies (i.e., asymptotic resolution of the oscilla- tion is not always enough).Worse enough, in these cases there is no warning from the scheme (like blow-up) that something went wrong (since local error control cannot be used anyway).It seems that the only safety anchor here lies in the analysis and, for more difficult problems than the linear Schr6dinger equation, in physical insight.
The paper is organized as follows.In Section 2 we present the basic analytical tools (e.g., Wigner transforms) needed to carry out the program on the discretization schemes, Section 3 is concerned with formulating the finite difference schemes in the language of pseudo-differential calculus and Section 4 contains the identification of the Wigner measures of the discretizations.Numerical sample computations and interpretations of the theory are given in Section 5.

Iklg
where k (k km) N m denotes a multiindex, K is the order of the differential operator (2.2) and kl "-k +... + km the order of the multiindex k.
The DO (2.2) can now be written as e(x,)( We denoted D -iOn. The convenience in the Weyl-calculus lies in the fact that a Weyl-operator is formally selfadjoint iff it has a real valued symbol (@ [8]).
Since we are interested in generalizations of the Schr6dinger-equation we assume for the following (i) Q is real valued for 0 Ik K, (ii) Vk, a N with Ikl KCk, > o.
(A1) IOQ(x)l c,x m.This implies in particular that Q(x, sD) w is selfadjoint on L 2(Nm).By Stone's Theorem exp(-i(t/s)Q(.,sD)w) is a strongly continuous group of unitary operators on L2(Nm).Thus we conclude the Ll(Nm)-conser- vation in time of the so-called position-density ,(x, 0 I(x, t)l: , , n(x, t) dx n (x) dx Vt N, (2.5)   where we set n "-ul 2. In quantum mechanics the wave function u c= uc(x, t) (i.e., the solution of the Schr6dinger- equation) is usually considered an auxiliary quantity.It facilitates the calculation of physical observables of the system under consideration [10]  corresponding to actual measurements.An obser- vable A , which depends on the position variable x and on the momentum operator eD, is given by the Weyl-operator A a(., eD) w (2.6) with the real valued symbol a(x, s()..7)Here (.,.) stands for the L2(Nm)-scalar product and, of course, it is assumed that u(t) lies in the domain of a(., eD) w.
A good framework for manipulating quantities which are quadratic in the wave function (e.g.(2.7)), is given by the Wigner-transform [7, 22].For given functions f, g ,9' (m) and a given scale e (0, 0] we define the Wigner-transform (on the scale ) by we(f,g)(x, C) (2.8)For fixed e this defines a bilinear continuous mapping from $,([m) x $,(Nm)into $'(N x N).

(A2)
Then there exists w L([t; A/l + (x m x n)) (where All + stands for the cone of positive Borel- measures) such that after selection of a subsequence and w satisfies the transport equation after selection of a subsequence.Since the limit process (2.12) is actually locally uniform in [7], the convergence (2.14) takes place in Cloc(Nt).Easy calculations give for the position density f t) / w (x, t)d. (2.15) Additional assumptions on the initial datum u have to be imposed in order to be able to pass to the limit e --, 0 in (2.15) (note that (2.14) cannot be applied directly since n corresponds to a momen- tum independent observable, whose symbol is not in S(x m x )).A complete account of this is given in [7].Here we only remark that holds if u is e-oscillatory, i.e., lim f Iu-()12d--0 (2.17) e0 Jll-> 0 w0 w0 Oi +{Q' } --0, m subject to the initial condition (2.13a) for every continuous, compactly supported func- tion .The assumption (2.17) says that u oscil- lates with wavelength at least O(e).
In order to avoid taking subsequences we shall assume for the following w/ is the unique Wigner-measure of u. (A3) In (2.13a) {., .}denotes the Poisson bracket: {f g} Of Oxg Oxf Og.

FINITE DIFFERENCE SCHEMES
Note that the convergence (2.12) holds for the whole sequence w e if w/ is the unique Wigner- measure of u (i.e., independent of the subscale The unique solution of (2.13) then allows the calculation of the limit e --, 0 of the average value of an observable A determined by a symbol a a(x, ) $ in the state u (cf.(2.7) and (2.10)).
We obtain Let F {# llal +'" + lmam ljG 7/for _<j_< m} [m be the lattice generated by the linearly independent vectors al,..., am m.For a multi- index k N' we construct a discretization of the order N of the operator O k as follows: Ea --gO a f a(x, )w(dx, d, t) (2.14)Here h6(0, ho] is the mesh-size, FkGF is the finite set of discretization points and a,k [ are coefficients satisfying a#,klZ ) where 61,k if k and 0 otherwise.It is an easy exercise to show that the local discretization error of (3.1) is O(hN) for all smooth functions if (D1) holds.For a detailed discussion of the linear prob- lem (D1) (i.e., possible choices of the coefficients a,k) we refer to [14].
Given the discretization (3.1) for 0 < Ik < K we now define the corresponding finite difference discretization of Q(., cD) w by applying (3.1) (with 0 iD) directly to (2.4).
Denoting au,ke Qk(x), we obtain the finite difference discretization of (2.4)   in the form Qh,e(x, cD) Wg(x) oikl(--O 'k' ZaMQk(x +?)(x + h#).Since Qh,e(x,D) w is a bounded operator on L 2(m), it is selfadjoint if i11 au,keiU. is real valued for 0 <_ Ikl K. (D2) As temporal discretizations we consider the one- step schemes with time step At > 0: eun+lAt-un + iQh#(x, eD)W(aun+ + flun) (3.4a) with a>_0, fl>_0, a+fl= 1.Here (and in the sequel) we denote the vector of small parameters by cr (, h, At).u is the obtained approximation of u(&) where t nAt, n E No.In this paper we will restrict ourselves to the analysis of the Crank- Nicolson scheme (a fl (1/2)).The (non time- reversible and non mass-conserving) implicit/ explicit Euler schemes are not investigated here but the interested reader may find the analysis in [15].Also we shall analyze the Leap-Frog method 2At + iQh'e(x'eD)WU+l The choice of will be discussed later on.
Note that the selfadjointness of Qh,(x, eD) w implies that the operator Id+iwQ,(x, eD) w is boundedly invertible on L2() for all w .T here- fore also the scheme (3.4) for > 0 gives well- for n=l 2. if defined approximations u u L2().Moreover we remark that it is suffici- ent to evaluate (3.4) and, resp., (3.5) at x hP in order to obtain discrete equations for {u(h>)l > P}.Clearly, artificial 'far out' boundary condi- tions have to be imposed for practical computa- tions.Their impact will not be taken into account in the subsequent analysis.
We now collect properties of the finite difference schemes.We start with the spatial discretization:  (3.6) Proof see [1 5].
Choosing h such that p (e/h) corre- sponds to asymptotically resolving the oscillations of wave-length O(h) of the solution u(t) of (2.1).
Iklg (3.7)   n+2 Rn In the case ,ho/__+ 0 (which corresponds to not resolving the oscillations asymptotically) we have h,e0 #cF0 and, thus, Qh,(x, eD) w does not approximate Q(x, eD) w.Therefore, we cannot expect reason- able numerical results in this case (which will not be investigated further).
The Crank-Nicolson scheme is unconditionally stable.The Leap-Frog scheme is stable if for some d > 0 sufficiently small.

CONVERGENCE OF WIGNER MEASURES
The consistency-stability concept of classical nu- merical analysis provides a framework for the convergence analysis of finite difference discretiza- tions of linear partial differential equations.Thus, for > 0 fixed it is easy to prove that the schemes (3.4) and, resp., (3.5) (with an appropriate choice of fi) are convergent of order N in space and order 2 in time if the solution u is sufficiently smooth and if the stability constraints on the mesh- sizes At and h of Lemma 3.1 are taken into account.
Therefore, again for fixed > 0 we conclude convergence of the same order for averages of the observables defined in (2.7) assuming that a is smooth.Due to the oscillatory nature of the solutions of (2.1) the local discretization error of the finite difference schemes and, consequently, also the global discretization error, generally tend to infinity as e tends to 0. The situation is still complicated by the fact that again due to the occurrence of O(e)-wavelength oscillations the solution u of (2.1) and their discrete approxima- tions u, which solve (3.4) or, resp., (3.5), generally only converge weakly in L 2(m) as ---+ 0 and, resp., a 0. The limit processes -0, a 0 do not commute with the quadratically nonlinear opera- tion which has to be carried out to compute the average values of observables.In practice one is interested in finding (mild) conditions on the mesh sizes h and At, in dependence of and the used discretization such that the average values of the observables in the discrete state converge (for e small) to E given by (2.7).
(4. la) The function E(t),tE +, then is defined by piecewise linear interpolation of the values E(t).
We want to find conditions on h, At (in the form, and again defining w(t), [, by piecewise linear interpolation of the values (4.2), we conclude from (2.10) that (4.1 b) is equivalent to proving liw(t)-w(t) in S'(x m x [?), (4.3) crGB locally uniformly in t, where w w(u , u) is the Wigner-transform of the solution u of (2.1).Note that w(tn) denotes the Wigner-transform of the finite difference solution u on the scale e.
In this Section we shall compute the accumula- tion points of w as cr --+ 0. We shall see that the set of Wigner-measures of the difference schemes {wl subsequence (,)of A (4.4) W 0 lim w ' " depends decisively on the discretization method and on the relative sizes of e,h and At.In those cases, in which W w( lim0w) holds, (4.3)   follows for rl B, while (4.3) does not hold if the measures W and w are different.
(i) (At/e)--O.Then the assertion of (1) remains true if Q is replaced by Qe in (4.9).
We remark that the convergence of w to W o is locally uniform in in Theorem 4.1, Theorem 4.2, Cases (1) and (2) (i) while it is only in L(Rt, $') weak-, in the Case (2) (ii) of Theorem 4.2.In the latter case we cannot generally exclude the occurrence of more than one accumulation point of w.
The two Theorems give necessary and sufficient conditions for a-{w w[u]} i.e., for the property of the difference schemes that their Wigner-measures are equal to the Wigner- measure associated to the IVP (2.1).Summarizing, these conditions are: (1) Crank-Nicolson scheme: (At/e) --O, (h/e) -+ O.
In all other cases there are initial data u for which either instabilities occur or such that the Wigner-measure ofthe difference scheme under con- sideration is different from the Wigner-measure of the Schr6dinger-type IVP.It can be shown that the term 1-w Q (appearing in (4.10a)) is strictly positive if the stability condition (3.8) is satisfied.
The proofs of the Theorems can be found in [15].

EXAMPLES AND NUMERICAL RESULTS
We consider the linear Schr6dinger equation (already in scaled form) which describes the trans- port of a charged particle under the influence of an electrostatic potential V(x)E C o(Rm) with uni- formly bounded derivatives c 2 eu -i-Au + iV(x)u e --0, x Rm, .(5.1a) (5.1b) Here e is the scaled Planck constant.In this case the symbol Q is simply 2 O(x, ) -Il + V(x). (5. 2) The initial condition (5.1 b) is chosen in WKB form u(x) v/ni(x) exp ( x [ m, (5.3)   with tti(X) and Si(x) independent of c, real valued, regular and with nl(x) decaying to zero sufficiently fast as [xl---, oc.
As already pointed out, the main goal when solving the Schr6dinger equation is to compute macroscopic quantities associated to the wave function u.The most important macroscopic quantities are the position density n and the current density J: n (x, t) lu (x, t)12, Je(x, t)"-Im (ue(x, t)Vue(x, t)).
The test problem (5.1) is considered in the one dimensional case and it is discretized in space with the usual three-point symmetric scheme, and in time with the Crank-Nicolson scheme given by (3.4) or with the Leap-Frog scheme given by (3.5).
We denote by u, the approximation of u (x#,tn), where t nAt, x# =jh is a given mesh point and cr (e, h, At).We recall the definition of the discrete position density and current density J+(xj, t,) e Im . . .u]+l'n j,n j,n h n O, 1,.,xj c P. (5.8b) We want to investigate in which situations the computed quantities n and J converge (weakly) to n o and j0, as cr 0. Therefore, in the pictures pre- sented in the following n and J are plotted and compared with n o and j0.As an example, Fig- ures 4-6 show a sequence of n and Ja converging to n o and j0.Most of the other figures illustrate cases where the sequence does not converge to n o and j0.The answer is (partially) given by Theorems 4.1 and 4.2.They give conditions on the mesh sizes h and At to guarantee that the Wigner measure of the discrete scheme W (defined in (4.4)) coincides with the Wigner measure of the continuous problem w (defined in (5.4)).As a direct consequence, all the observables of the finite difference scheme under consideration converge to the exact observables.
Cases where W does not converge to w , but n converges to n o are possible.In the continuous case, n o and j0 are recovered as moments of the Wigner measure (see (5.5)).In the discrete case the situation is more delicate.For the Crank-Nicolson scheme, Corollary 4.1 states that n---v := (x,d(,t) if the initial datum u s e-oscillatory.
As already mentioned, the WKB initial datum (5.3) implies the e-oscillatory property, so if W coincides with w , then convergence of n to n o is guaranteed, otherwise n can converge to something else.The main ingredient to prove Corollary 4.1 is that the e-oscillatory property of the initial datum is preserved by the scheme (and for Crank-Nicolson this follows from the L2-conservation).For the Leap-Frog scheme the e-oscillatory property is not preserved in general, because of the lack of L2-conservation.However, in the case of constant coefficient operators the c-oscillatory property of the solution is preserved in time.Thus, an analogous result as in Corollary 4.1 holds true in this case for all the time discretization schemes considered here.Moreover, for constant coefficient operators a complete characterization of the convergence of the current density is possible.(see [12]) In all the numerical tests presented here, the com- putations are carried out in the interval [0, 1] and periodic boundary conditions are used.The poten- tial V-0 is chosen, unless otherwise specified.
All the pictures presented in the following show n and J at time 0.54 (after caustics have developed).The reference n o and jo are the ones of Figure 3.
Crank-Nieolson The first set of tests refers to the Crank-Nicolson time discretization (3.4) with a =/3 0.5.The Crank-Nicolson scheme is one of the most widely used finite difference schemes for the discretization of the Schr6dinger equation.
The main reason lies in its conservation properties: the position density and the total energy ((c2/2) f Vx u'(x, tn) lZdx +f V(x) lu'(x, tn)l 2 dx)are pre- served.These conservation properties are important, but by far not enough to guarantee good convergence for small c. scheme verifies (4.5) and consequently the transport equation (2.13) (see also (5.4)) of the continuous Wigner measure.The functions n and J' oscillate about n o and j0 (resp.)with wave length e.It is evident that when e is halved, also the wave length is halved.The amplitude of the oscillations does not grow as e becomes smaller, except the first and the last one which increase with c.Thus, the pictures indicate that the sequences {n} and {J} converge weakly to n o and j0.We can also say that n and J" in Figures 4-6 are very good approximations of n and J for the selected e's.When the step sizes violates (1) (i), a patholo- gical behavior is expected, since the Wigner measure of the scheme identified in Theorem 4.1 does not coincide with the continuous one.Here we observe the effect on macroscopic quantities.
We consider first the case when At is still small enough (At e5), but the mesh size h is equal to e (Case (2) (i)).that n and J converge to a function with smaller support than n and j0.Moreover, comparing these plots with those obtained with h e.2, we observe that, for the same fixed e, the results are qualitatively different.The support of n and J in   ) is smaller than the one of the functions in Figure 5, (6, resp.) and consequently, due to charge conservation, the oscillation ampli- tude is larger, n and J of Figures 7, 8 are bad approximations of n and JC for a fixed small e.The situation is similar when h e 1"2 and At 3e are selected (Case (1) (ii)).Again the solutions exhibit smaller support and higher amplitude density density FIGURE 16 Leap-Frog, e 10-3, At 0.4e, h e, V 0.
1.5 0 0.1 0 0.2 FIGURE 18 Leap-Frog, e 10-3, At (1/30)e, h e, V 10, 0.2 < x < 0.8, 0 < < 0.1.0.8 oscillations (see Fig. 9, for e-0.5.10-3).The effect is amplified when both h and At are too large.Figure 10 illustrates case (2) (ii) with h e, At 3e, for e 0.5.10-3.In order to illustrate case (1) (iii), condition (4.6a) must be fulfilled.A constant potential V > 0 is considered, so that Q is bounded away from zero.Figures 11-13 show the results when V 0.25, h e2, At e6, for a sequence of e's.When e is large, (the plot for e 8.10 -3 is presented in Fig. 11) it is evident that n and J are far away from the correct solution, but the pictures do not suggest what the limit of the sequence is.As e becomes smaller, n does not exhibit oscillations and its maximum decreases with e (see Fig. 12) which refers to a computation with e-10-3.Figure 13 displays (for e-0.5.10-4) n and J very close to the initial value n1 and J:=n(d/dx)Si, in agreement with the theoretical result W wi.
Leap-Frog The Leap-Frog scheme (3.5) has some interesting properties when being applied to Schr6dinger type equations.It is an explicit scheme, but with a not too restrictive stability condition.where gma maxx IV(x) I.Moreover, although it is not a conservative scheme it has similar properties: 'Charge conservation': Re (u(x, tn+l)bla(X, tn))dx const 'Energy conservation" 2 Re(VxU(X'tn)VxU(X'tn+l))dx+ V(x)Re (u(x, tn)u(x, tn+l))dx const.
All the tests presented here are carried out with -uu, so that z/ w/ (see Theorem 4.2).
Let us discuss case (1) ofTheorem 4.2 first.Only a condition on h is stated, since the stability condition (5.10) implies that (At/e) goes to zero in this case.
For instance, if V 0 is selected, then the stability condition forces (At/e)< 0.5(h2/e2).For h e 1"2 and At 0.49e TM the stability condition is satis- fied and the results for e 2.10-3, e 10 -3 e 0.5.10 -3 are very similar to those of Figures [4][5][6] (obtained with the Crank-Nicolson scheme) and are not reported.Note that in this case the stability condition allows a At very close to At e 5 used in the tests with the Crank-Nicolson scheme.Of course, if h is much smaller than e, (5.10) requires a much smaller At than the one requested by condition (i) of Theorem 4.1 for Crank-Nicolson.
All the pathological situations are possible for h O(e).As in the Crank-Nicolson case, when Qp appears in the transport Eq. (4.9) (case (2) (i)), n and J exhibit smaller support and higher oscilla- tions than in the case when the exact Q is present.
Figures 14 and 15 show the results for h e, /t e 1"25 with e 10 .3 and e 0.5.10-3.Comparing Figure 15 and Figure 8, which corresponds to case (2) (i) for Crank-Nicolson, no difference is observed.
In the case (2) (ii), we remark that the term cd2Q],e in (4.10a), with c At/e, is always posi- tive when the stability condition (5.10) is satisfied.Instead, when this coefficient is small, Z and W are significantly different and they evolve accord- ing to two different laws.The effect on the macroscopic quantities is weird.If, for example, e= 10-3,At=(1/30)e,h=eand V= 10thenthe maximum of n oscillates almost every time step.
Figure 18 shows the evolution of n'(x,t) for 0_< < 0.1.Although the density n is smooth in this time intervall (we recall that caustics develop at time 0.2) the wave function u (solution of the Schr6dinger equation) has oscillations in the phase.A choice of the discretization parameters that does not take care of this fact can give the wrong results.Actually, Theorems 4.1 and 4.2 apply regardless of the occurance of caustics.

correspond to e 10 - 3
and e 0.5.10-3.The pictures indicate

Figures 16 and 17
Figures16 and 17 show results with e 10 -3 and e 0.5.10 -3 for At 0.4e, h e and V 0 which are very similar to those of Figures14 and 15, since in this case 1-co 2g12 is almost h,e