On the classical limit of a time-dependent self-consistent field system: analysis and computation

We consider a coupled system of Schr\"odinger equations, arising in quantum mechanics via the so-called time-dependent self-consistent field method. Using Wigner transformation techniques we study the corresponding classical limit dynamics in two cases. In the first case, the classical limit is only taken in one of the two equations, leading to a mixed quantum-classical model which is closely connected to the well-known Ehrenfest method in molecular dynamics. In the second case, the classical limit of the full system is rigorously established, resulting in a system of coupled Vlasov-type equations. In the second part of our work, we provide a numerical study of the coupled semi-classically scaled Schr\"odinger equations and of the mixed quantum-classical model obtained via Ehrenfest's method. A second order (in time) method is introduced for each case. We show that the proposed methods allow time steps independent of the semi-classical parameter(s) while still capturing the correct behavior of physical observables. It also becomes clear that the order of accuracy of our methods can be improved in a straightforward way.


264
SHI JIN, CHRISTOF SPARBER AND ZHENNAN ZHOU 1. Introduction. The numerical simulation of many chemical, physical, and biochemical phenomena requires the direct simulation of dynamical processes within large systems involving quantum mechanical effects. However, if the entire system is treated quantum mechanically, the numerical simulations are often restricted to relatively small model problems on short time scales due to the formidable computational cost. In order to overcome this difficulty, a basic idea is to separate the involved degrees of freedom into two different categories: one, which involves variables that behave effectively classically (i.e., evolving on slow time-and large spatial scales) and one which encapsulates the (fast) quantum mechanical dynamics within a certain portion of the full system. For example, for a system consisting of many molecules, one might designate the electrons as the fast degrees of freedom and the atomic nuclei as the slow degrees of freedom.
Whereas separation of the whole system into a classical part and a quantum mechanical part is certainly not an easy task, it is, by now, widely studied in the physics literature and often leads to what is called time-dependent self-consistent field equations (TDSCF), see, e.g., [4,5,9,13,18,19,21,26] and the references therein. In the TDSCF method, one typically assumes that the total wave function of the system Ψ(X, t), with X = (x, y), can be approximated by Ψ(X, t) ≈ ψ(x, t)ϕ(y, t), (1.1) where x and y denote the degrees of freedom within a certain subsystem, only. The precise nature of this approximation thereby strongly depends on the concrete problem at hand (in particular, on the initial data and on the precise nature of the coupling between the two subsystems). Disregarding this issue for the moment, one might then, in a second step, hope to derive a self-consistently coupled system for ψ and ϕ and approximate it, at least partially, by the associated classical dynamics. In this article we will study a simple model problem for such a TDSCF system, motivated by [4,9,18,19,21], but one expects that our findings extend to other self-consistent models as well. We will be interested in deriving various (semi-) classical approximations to the considered TDSCF system, resulting in either a mixed quantum-classical model, or a fully classical model. As we shall see, this also gives a rigorous justification of what is known as the Ehrenfest method in the physics literature, cf. [5,9]. To this end, we shall be heavily relying on Wigner transformation techniques, developed in [11,20], which have been proved to be superior in many aspects to the more classical WKB approximations, see, e.g. [24] for a broader discussion. One should note that the use of Wigner methods to study the classical limit of nonlinear (self-consistent) quantum mechanical models is not straightforward and usually requires additional assumptions on the quantum state, cf. [22,20]. It turns out that in our case we can get by without them.
In the second part of this article we shall then be interested in designing an efficient and accurate numerical method which allows us to pass to the classical limit in the TDSCF system within our numerical algorithm. We will be particularly interested in the meshing strategy required to accurately approximate the wave functions, or to capture the correct physical observables (which are quadratic quantities of the wave function). To this end, we propose a second order (in time) method based on an operator splitting and a spectral approximation of the TDSCF equations as well as the obtained Ehrenfest model. These types of methods have been proven to be very effective in earlier numerical studies, see, e.g., [1,2,7,16,17] for previous results and [15] for a review of the current state-of-the art of numerical methods for semi-classical Schrödinger type models.The readers may also refer to [3,27] for some recent results on the numerical analysis of Born-Oppenheimer molecular dynamics with connections to the Ehrenfest model. In comparison to the case of a single (semi-classical) nonlinear Schrödinger equation with power law nonlinearities, where one has to use time steps which are comparable to the size of the small semi-classical parameter (see [1]), it turns out that in our case, despite of the nonlinearity, we can rigorously justify that one can take time steps independent of the semi-classical parameter and still capture the correct classical limit of physical observables.
The rest of this paper is now organized as follows: In Section 2, we present the considered TDSCF system and discuss some of its basic mathematical properties, which will be used later on. In Section 3, a brief introduction to the Wigner transforms and Wigner measures is given. In Section 4 we study the semi-classical limit, resulting in a mixed quantum-classical limit system. In Section 5 the completely classical approximation of the TDSFC system is studied by means of two different limiting processes, both of which result in the same classical model. The numerical methods used for the TDSCF equations and the Ehrenfest equations are then introduced in Section 6. Finally, we study several numerical tests cases in Section 7 in order to verify the properties of our methods.

2.1.
Basic set-up and properties. In the following, we take x ∈ R d , y ∈ R n , with d, n ∈ N, and denote by ·, · L 2 x and ·, · L 2 y the usual inner product in L 2 (R d x ) and L 2 (R n y ), respectively, i.e.
The total Hamiltonian of the system acting on L 2 (R d+n ) is assumed to be of the form where V (x, y) ∈ R is some (time-independent) real-valued potential. Typically, one has V (x, y) = V 1 (x) + V 2 (y) + W (x, y), (2.2) where V 1,2 are external potentials acting only on the respective subsystem and W represents an internal coupling potential in between the two subsystems. From now on, we shall assume that V satisfies where here and in the following, we denote by C 0 the set of continuous functions vanishing at infinity.
Remark 2.1. For potential bounded below, the requirement V 0 is not really an assumption, but merely corresponds to fixing the point 0 on the energy axis.
In (2.1), the Hamiltonian is already written in dimensionless form, such that only two (small) parameters ε, δ > 0 remain. In the following, they play the role of dimensionless Planck's constants. Dependence with respect to these parameters will be denoted by superscripts. The TDSCF system at hand is then (formally given by [9]) the following system of self-consistently coupled Schrödinger equations where we denote by the Hamiltonian of the subsystem represented by the x-variables (considered as the purely quantum mechanical variables) and in which y only enters as a parameter. It is obtained by substituting the ansatz (1.1) into the original Schrödinger equation and integrating over y and x respectively, see [9]. As a matter of fact, the TDSCF systems may take various forms, which are also equivalent to one another by certain gauge transformations, see [4,18,19,21] for broad discussions. Without loss of generality, we choose to study the specific TDSCF system (2.3).
For simplicity, we assume that at t = 0 the data ψ δ in only depends on δ, and that ϕ ε in only depends on ε, which means that the simultaneous dependence on both parameters is only induced by the time-evolution. Remark 2.2. A typical example of initial data which satisfies this assumption (and all upcoming requirements of our analysis) is Ψ(X, 0) ≈ ψ δ in (x)ϕ ε in (y) = a 1 (x)e iS1(x)/δ a 2 (y)e iS2(y)/ε , where S 1 , S 2 are some smooth, real-valued phases and a 1 , a 2 some (in general, complex-valued) amplitudes. In other words, Ψ(X, 0) is assumed to be approximated by a (two-scale) WKB type initial data in product form.
Finally, the coupling terms are explicitly given by and after formally integrating by parts Throughout this work we will always interpret the term ψ ε,δ , h δ ψ ε,δ L 2 x as above, i.e., in the weak sense. Both Υ ε,δ and Λ ε,δ are time-dependent, real-valued potentials, computed self-consistently via the dynamics of ϕ ε,δ and ψ ε,δ , respectively. Note that Here, the purely time-dependent part ϑ ε,δ (t), could in principle be absorbed into the definition of ϕ ε,δ via a gauge transformation, i.e., For the sake of simplicity, we shall refrain from doing so, but this nevertheless shows that the two coupling terms are in essence of the same form. Also note that this gauge transform leaves any H s (R d )-norm of ϕ ε,δ invariant (but clearly depends on the solution of the second equation within the TDSCF system). An important physical quantity is the total mass of the system, where m ε,δ 1 , m ε,δ 2 denote the masses of the respective subsystem. One can then prove that these are conserved by the time-evolution of (2.3).
Proof. Assuming for the moment, that both ψ ε,δ and ϕ ε,δ are sufficiently smooth and decaying, we multiply the first equation in (2.3) with ψ ε,δ and formally integrate with respect to x ∈ R d x . Taking the real part of the resulting expression and having in mind that Υ ε,δ (y, t) ∈ R, yields which, after another integration in time, is the desired result for m ε,δ 1 (t). By the same argument one can show the result for m ε,δ 2 (t). Integration in time in combination with a density argument then allows to extend the result to more general solutions in H 1 , respectively.
We shall, from now on assume that the initial data is normalized such that m ε,δ 1 (0) = m ε,δ 2 (0) = 1. Using this normalization, the total energy of the system can be written as Note that in view of our assumption (A1) on V this is well-defined and that E ε,δ (t) is, in fact, a sum of three non-negative terms.
Lemma 2.5. Let V satisfy (A1) and assume that ψ ε,δ ∈ C(R t ; H 1 (R d x )) and ϕ ε,δ ∈ C(R t ; H 1 (R n y )) solve (2.3). Then In other words, we have conservation of the total energy, which in itself implies a bound on the interaction energy (since V 0) and on the kinetic energies of the respective subsystems. Note however, that the energies of the respective subsystems are in general not conserved, unless W ≡ 0, i.e., V (x, y) = V 1 (x) + V 2 (y).

SHI JIN, CHRISTOF SPARBER AND ZHENNAN ZHOU
Proof. Assuming, as before that ψ ε,δ and ϕ ε,δ are sufficiently regular (and decaying), the proof is a lengthy but straightforward calculation. More precisely, using the shorthand notation where we denote We will now show that (I) + (III) = 0. By using (2.3), one gets where [A, B] := AB − BA denotes the commutator bracket. Similarly, one finds that Analogously, one can show that (II) + (IV) = 0 and hence, an integration in time yields E ε,δ (t) = E ε,δ (0). Using a density arguments allows to extend this result to more general solution in H 1 .

Existence of solutions.
In this subsection, we shall establish global in-time existence of solutions to the TDSCF system (2.3). Since the dependence on ε and δ does not play a role here, we shall suppress their appearance for the sake of notation.
Proposition 2.6. Let V satisfy (A1) and ψ in ∈ H 1 (R d x ), ϕ in ∈ H 1 (R n y ). Then there exists a global strong solution (ψ, ϕ) ∈ C(R t ; H 1 (R d+n )) of (2.3), satisfying the conservation laws for mass and energy, as stated above.
Proof. We shall first prove local (in-time) well-posedness of the initial value problem (2.3): To this end, we consider Ψ(·, t) = (ψ(·, t), ϕ(·, t)) : R d+n → C 2 and define the associated y , and consequently set H 1 (R d+n ) := {Ψ ∈ L 2 (R d+n ) : |∇Ψ| ∈ L 2 (R d+n )}. Using this notation, the TDSCF system (2.3) can be written as Clearly, H is the generator of a strongly continuous unitary Schrödinger group U (t) := e −itH , which can be used to rewrite the system using Duhamel's formula as Following classical semi-group arguments, cf. [8], it suffices to show that f (Ψ) is locally Lipschitz in H 1 (R d+n ) in order to infer the existence of a unique local in-time solution Ψ ∈ C([0, T ), H 1 (R d+n )). This is not hard to show, since: by assumption. Now we can estimate An analogous argument can be done for the second part of f (Ψ), since Combining all these estimates, we conclude The same reasoning can then be applied to ∇f (Ψ) L 2 , by noticing ∇ ϕ, V ϕ L 2 y = ϕ, ∇ x V ϕ L 2 y and ∇ ψ, hψ L 2 in view of (2.5). Since V satisfies (A1), these expressions are all well-defined. In summary, one gets that there exists a C = C( V L ∞ , Using this, [8,Theorem 3.3.9] implies the existence of a T = T ( Ψ H 1 ) > 0 and a unique solution Ψ ∈ C([0, T ), H 1 (R d+n )) of (2.9). It is then also clear, that this solution satisfies the conservation of mass and energy for all t ∈ [0, T ). Moreover, the quoted theorem also implies that if T < +∞, then However, having in mind the conservation laws for mass and energy stated in Lemmas 2.4 and 2.5 together with the fact that we assume w.l.o.g. V (x, y) 0, we immediately infer that Ψ(·, t) H 1 C for all t ∈ R and hence, the blow-up alternative (2.10) implies global in-time existence of the obtained solution.
Remark 2.7. Note that this existence result rests on the fact that the term Λ ε,δ (y, t) := ψ ε,δ , h δ ψ ε,δ L 2 x is interpreted in a weak sense, see (2.5). In order to interpret it in a strong sense, one would need to require higher regularity, in particular ψ ε,δ ∈ H 2 (R d x ).
3. Review of Wigner transforms and Wigner measures. The use of Wigner transformation and Wigner measures in the analysis of (semi)-classical asymptotic is, by now, very well established. We shall in the following, briefly recall the main results developed in [20,11] (see also [10,22,24] for further applications and discussions of Wigner measures): Denote by {f ε } 0<ε 1 a family of functions f ε ∈ L 2 (R d x ), depending continuously on a small parameter ε > 0, and by the corresponding Fourier transform. The associated ε-scaled Wigner transform is then given by [28]: Clearly, one has and thus Plancherel's theorem together with a simple change of variables yields . The real-valued function w ε (x, ξ) acts as a quantum mechanical analogue for classical phase-space distributions. However, w ε (x, ξ) 0 in general. A straightforward computation shows that the position density associated to f ε can be computed via Moreover, by taking higher order moments in ξ one (formally) finds In order to make these computations rigorous, the integrals on the r.h.s. have to be understood in an appropriate sense, since w ε ∈ L 1 (R m x × R m ξ ) in general, cf. [11,20] for more details.
It has been proved in [20,Proposition III.1], that if f ε is uniformly bounded in where C 0 is an ε-independent constant, then the set of Wigner functions {w ε } 0<ε 1 is uniformly bounded in A . The latter is the dual of the following Banach space denotes the space of continuous functions vanishing at infinity. More precisely, one finds that for any test function Denoting, we therefore obtain const., uniformly in ε. Thus, up to extraction of sub-sequences {ε n } n∈N , with ε n → 0 + as n → ∞, there exists a limiting object It turns out that the limit is in fact a non-negative, bounded Borel measure on Remark 3.1. One easily checks that the Schwartz space S is in fact dense in A. Thus, it would also be possible to state all the convergence results above in terms of convergence in S (R d x × R n y ). This is the framework used in [11]. If, in addition it also holds that f ε is ε-oscillatory, i.e., then one also gets (up to extraction of sub-sequences) i.e., for any test function σ ∈ C 0 (R d x ): Indeed, the Wigner measure µ is known to encode the classical limit of all physical observables. More precisely, for the expectation value of any Weyl-quantized operator Op ε (a), corresponding to a sufficiently "nice" classical symbol, say a( and hence where the right hand side resembles the usual formula from classical statistical mechanics.
In order to describe the dynamics of Wigner measures, we first recall that if is a pseudo-differential operator describing the influence of a potential U (x) ∈ R. Explicitly, Θ ε [U ] is given by [20]: Here, the symbol δU ε is found to be and, thus, under sufficient regularity assumptions on V , one consequently obtains It consequently follows that the measure µ(x, ξ, t) solves Liouville's equation on phase space, i.e.
in the sense of distributions. Here, µ in is the weak * limit of w ε in in A , along subsequences of (ε n ) n∈N (which, in principle, could all yield different limits).

Remark 3.2.
In fact, a more general formula for the asymptotic, or semi-classical expansion (in powers of ε) of any Wigner transformed Schrödinger-type equation is available in [11,Proposition 1.8]. This formula will be used in the numerical algorithm described below.
Finally, we note that for sufficiently regular potential V , one can improve the convergence statements and show that, indeed, where Φ t : R 2d → R 2d is the Hamiltonian flow associated to (3.5): This allows to prove uniqueness of the weak solution of (3.5), provided the initial measure µ in is the same for all sub-sequences 4. The mixed quantum-classical limit. In this section we will investigate the semi-classical limit of the TDSCF system (2.3), which corresponds to the case ε → 0 + and δ = O(1) fixed. In other words, we want to pass to the classical limit in the equation for ϕ ε,δ only, while retaining the full quantum mechanical dynamics for ψ ε,δ .
The standing assumption from now on, until the end of this work will be that the initial data ϕ ε in ∈ H 1 (R n y ) and In other words, the initial data are assumed to be such that the initial mass and the initial energy are uniformly bounded with respect to both ε and δ. In view of Lemma 2.4 and Lemma 2.5 this property consequently holds true for all times t 0, and hence, neither the mass, nor the energy can become infinite in the classical limit. In addition, we will assume, for simplicity, that the individual masses of the two sub-systems are initially normalized such that m ε,δ 1 (0) = m ε,δ 2 (0) = 1. This normalization is henceforth preserved by the time-evolution of (2.3).
Next, we introduce the ε-scaled Wigner transformation of ϕ ε,δ in the form (In this subsection, we thus could, in principle, suppress the dependence on δ completely, since it is fixed, but given that we will consider the subsequent δ → 0 + limit in Section 5, we shall keep its appearance within the superscript.) The assumption (A2), together with the a-priori estimates established in Lemma 2.4 and Lemma 2.5 then implies the uniform bound, for any t ∈ R, where C 0 is a constant independent of ε and δ. In other words, ϕ ε,δ (·, t) is εoscillatory for all times and we consequently infer the existence of a limiting Wigner measure µ 0,δ (y, η, t) ≡ µ δ (y, η, t) such that (up to extraction of sub-sequences) for all t ∈ [0, T ] it holds The measure µ δ encodes the classical limit of the subsystem described by the yvariables only.
In order to proceed, we will need to strengthen our convergence results with respect to the t-variable. To this end, we recall that since ϕ ε,δ solves the second equation in the TDSCF system (2.3), w ε,δ ≡ w ε [ϕ ε,δ ] solves the corresponding Wigner transformed equation where Θ[Λ ε,δ ] is explicitly given by The associated symbol δΛ ε,δ reads in view of the definition of Λ ε,δ , given in (2.5). Introducing the short hand notation . In particular, this shows that the purely time-dependent term ϑ ε,δ (t) appearing in (2.5) does not contribute to the symbol of the pseudodifferential operator Θ. can also be seen by using the time-dependent gauge transformation (2.6) from the beginning.
We can now prove the following lemma.
Lemma 4.2. Let Assumptions (A1) and (A2) hold. Then w ε,δ is equi-continuous in time and hence, up to extraction of sub-sequences, we have Proof. The proof follows along the lines of [20,11]. In order to infer the assertion of the Lemma it is sufficient to show that ∂ t w ε,δ ∈ L ∞ ((0, T ); A (R n y × R n η )). The latter implies time-equicontinuity of w ε,δ and hence, the Arzela-Ascoli Theorem guarantees that there exists a subsequence {ε n } n∈N , with ε n → 0 + as n → ∞, such that w εn,δ converges uniformly on compact subsets of R t .
In order to prove the uniform bound on ∂ t w ε,δ we consider the weak formulation of (4.2), i.e.
for any test function χ ∈ A(R n y × R n η ). We shall only show how to bound the term Θ[Λ ε,δ ]w ε,δ , χ , since the other term on the right hand of (4.4) can be treated similarly.
To this end, let χ ∈ A(R n y × R n η ) be a smooth test function with the property that its Fourier transform with respect to η, i.e.
has compact support with respect to both y and z. This kind of test functions are dense in A and hence it suffices to show the assertion for these type of χ only. A straightforward calculation (cf. the proof of [20, Theorem IV.1]) shows that having in mind that V ∈ L ∞ , by assumption, and that ψ ε,δ (·, t) L 2 x = 1, ∀ t ∈ R, due to mass conservation. Since V satisfies (A1), the same argument also applies to ∇ y V ε,δ , which by dominated convergence, is simply given by x . Having in mind the computation (3.2), we know that We thus conclude that uniformly in ε and t, since due to the compact support of χ.
A similar argument can be done to obtain a uniform bound on | η · ∇ y w ε,δ , χ |, and thus (4.4) implies that uniformly in ε and t, and we are done.
Next, we look at the nonlinear coupling term appearing in the first equation of our TDSCF system (2.3): We first note that, as before, ) and an analogous bound for ∂ α x Υ ε,δ , |α| 2. In addition, we have that Υ ε,δ is continuous in time, since by triangle inequality where we have again used the fact that ϕ ε,δ (·, t) L 2 y = 1, for all t ∈ R. In view of our existence result stated in Proposition 2.6, the right hand side is continuous in time, and thus Υ ε,δ (x, ·) is too.
Proof. Having in mind that test functions of the form V (x, y) = γ(x)σ(y) are dense in C 2 0 (R d x × R n y ), the weak measure convergence (4.1) implies that for all (x, t) it holds point-wise for all (x, t) ∈ R d+1 . In addition, the foregoing Lemma shows that the point-wise convergence, in fact, also holds uniformly on compact time-intervals. Using that µ δ 0, we find since ϕ ε,δ (·, t) L 2 = 1, for all t ∈ R. (Here we have used [20, Theorem III.1] in the second inequality.) An analogous bound also holds for ∂ α x Υ δ and |α| 2, since V satisfies (A1). By applying the push-forward formula (4.8) with χ(x, y) = V (x, y), it is easy to see that Υ δ (t, ·) is continuous in time, yielding Υ δ ∈ C b (R t ; C 2 0 (R d x )). The following Proposition then shows that the solution of the first equation within the TDSCF system (2.3) stays close to the one where the potential Υ ε,δ is replaced by its classical limit Υ δ .
Proposition 4.4. Let V satisfy (A1) and ψ ε,δ , ψ δ ∈ C(R t ; H 1 (R d x )) solve, respectively Proof. Denote the Hamiltonian operators corresponding to the above equations by In view of our assumptions on the potential V and the existence result given in Proposition 2.6, we infer that H 1 and H 2 are essentially self-adjoint on L 2 (R d x ) and hence they generate unitary propagators U ε,δ 1 (t, s) and U δ 2 (t, s), such that . Computing further, one gets By Minkowski's inequality, one therefore has which, firstly, implies continuity of the difference in L 2 norm w.r.t. to t ∈ R and, secondly, we also have In order to identify the limiting measure µ δ we shall derive the corresponding evolutionary system, by passing to the limit ε → 0 + in (4.2). The main difference to the case of a given potential V (as studied in, e.g., [11]) is that here V ε,δ itself depends on ε and is computed self-consistently from the solution of ψ ε,δ . We nevertheless shall prove in the following proposition that the limit of Θ[V ε,δ ] as ε → 0 + is indeed what one would formally expect it to be.
where in the second inequality we have used the Cauchy-Schwarz inequality together with the fact that ||a| 2 − |b| 2 | |a − b|(|a| + |b|) for any a, b ∈ C. The strong L 2convergence of ψ ε,δ stated in Proposition 4.4 therefore implies x (y, t), pointwise in y and uniformly on compact time-intervals. Analogously we infer x (y, t). Next, we recall from (4.5) that V ε,δ and ∇ y V ε,δ are uniformly bounded in ε. Moreover, by using the Mean-Value Theorem, we can estimate This shows that F ε,δ := −∇ y V ε,δ is equicontinuous in y, and hence the Arzela-Ascoli Theorem guarantees that there exists a subsequence, such that F ε,δ converges, as ε → 0 + , uniformly on compact sets in y, t. Recalling from before that for any For χ having compact support the uniform convergence of F ε,δ then allows us to conclude Ξ ε,δ ε→0+ −→ i∇ y V δ (y, t) · F −1 z→η (z χ(y, z))(y, η) = F δ (y, t) · ∇ η χ(y, η), and since these χ are dense in A the result follows. Remark 4.6. One should note that, even though Λ ε,δ is a self-consistent potential, depending nonlinearly upon the solution ψ ε,δ , the convergence proof given above is very similar to the case [20, Theorem IV.2], which treats the classical limit of nonlinear Hartree-type models with smooth convolution kernels. In particular, we do not require to pass to the mixed state formulation which is needed to establish the classical limit in other self-consistent quantum dynamical models such as [22].
In summary, this leads to the first main result of our work, which shows that the solution to (2.3), as ε → 0 + (and with δ = O(1) fixed) converges to a mixed quantum-classical system, consisting of a Schrödinger equation for the x-variables and a classical Liouville equation for the y-variables.

(4.7)
Here µ in is the initial Wigner measure obtained as the weak * limit of w ε [ϕ ε in ] and Proof. The result follows from Proposition 4.4 and Proposition 4.5.

SHI JIN, CHRISTOF SPARBER AND ZHENNAN ZHOU
and the mixed quantum-classical system becomes with y 0 , η 0 ∈ R n . This is a well-known model in the physics and quantum chemistry literature, usually referred to as Ehrenfest method. It has been studied in, e.g, [6,5] in the context of quantum molecular dynamics.
Remark 4.8. A closely related scaling-limit is obtained in the case where the time-derivatives in both equations of (2.3) are scaled by the same factor ε. At least formally, this leads to an Ehrenfest-type model similar to (4.9), but with a stationary instead of a time-dependent Schrödinger equation, cf. [5,9]. In this case, connections to the Born-Oppenheimer approximation of quantum molecular dynamics become apparent, see, e.g., [25]. From the mathematical point of view this scaling regime combines the classical limit for the subsystem described by the y-variables with a time-adiabatic limit for the subsystem described in x. However, due to the nonlinear coupling within the TDSCF system (2.3) this scaling limit is highly nontrivial and will be the main focus of a future work. 5. The fully classical limit. In order to get a better understanding (in particular for the expected numerical treatment of our model), we will now turn to the question of how to obtain a completely classical approximation for the system (2.3). There are at least two possible ways to do so. One is to consider the limit δ → 0 + in the obtained mixed quantum-classical system (4.7), which in itself corresponds to the iterated limit ε → 0 + and then δ → 0 + of (2.3), cf. Section 5.1. Another possibility is to take ε = δ → 0 + in (2.3), which corresponds to a kind of "diagonal limit" in the ε, δ parameter space, cf. Section 5.2.
5.1. The classical limit of the mixed quantum-classical system. In this section we shall perform the limit δ → 0 + of the obtained mixed quantum-classical system (4.7). To this end, we first introduce the δ-scaled Wigner transform of ψ δ : Assumption (A2) and the results of Lemma 2.4 and Lemma 2.5 imply that ψ δ is a family of δ-oscillatory functions, i.e, sup 0<δ 1 and thus there exists a limiting measure ν ∈ M + (R d x × R d ξ ), such that for all t ∈ [0, T ]: By Wigner transforming the first equation in the mixed quantum-classical system (4.7), we find that W δ [ψ δ ] ≡ W δ satisfies Having in mind that Υ δ ∈ C(R t ; C 2 0 (R d x )), the same arguments as in Lemma 4.2 then allow us to obtain a uniform bound on ∂ t W , and hence time-equicontinuity of W δ , which yields Furthermore, our assumptions one V together with the weak measure convergence (5.2) imply pointwise. By the same arguments as in the proof of Proposition 4.5, we find that this convergence holds uniformly on compact sets in y, t. With this in hand, we can prove the following result.
To obtain the convergence of the term Θ[Υ δ ]W δ , we note that with the convergence of the Wigner measure µ δ , which is obtained in Proposition 5.1, one gets point-wise, for all x, t. Similar to previous cases, one concludes that, up to extraction of possibly another sub-sequence, Υ δ converges, as δ → 0 + , uniformly on compact sets in x, t.
With the same techniques as in the proof of Proposition 4.5, one can then derive the equation for the associated Wigner measure ν. The classical limit of the mixed quantum-classical system can thus be summarized as follows.

(5.3)
Here ν in is the initial Wigner measure obtained as the weak * limit of W δ [ψ δ in ], and Remark 5.3. Note that system (5.3) admits a special solution of the form where x(t), y(t), ξ(t), η(t) solve the following Hamiltonian system: This describes the case of two classical point particles interacting with each other via V (x, y). Obviously, if V (x, y) = V 1 (x) + V 2 (y), the system completely decouples and one obtains the dynamics of two independent point particles under the influence of their respective external forces.
5.2. The classical limit of the TDSCF system. In this section we shall set ε = δ and consider the now fully semi-classically scaled TDSCF system where only 0 < ε 1 appears as a small dimensionless parameter: where, as in (2.4), we denote We shall introduce the associated ε-scaled Wigner transformations w ε [ϕ ε ](y, η, t) and W ε [ψ ε ](x, ξ, t) defined by (3.1). From the a-priori estimates established in Lemmas 2.4 and 2.5, we infer that both ψ ε and ϕ ε are ε-oscillatory and thus we immediately infer the existence of the associated limiting Wigner measures µ, ν ∈ M + , such that

The associated Wigner transformed system is
(5.5) By following the same arguments as before, we conclude that, up to extraction of sub-sequences, on compact sets in (x, t) and (y, t) respectively. Consequently, one can show the convergences of the terms Θ[Υ ε ]W ε and Θ[V ε ]w ε by the same techniques as in the proof of Proposition 4.5. In summary, we obtain the following result: Theorem 5.4. Let Assumptions (A1) and (A2) hold. Then, for any T > 0, we have that W ε and w ε converge as ε → 0 + , respectively, to µ ∈ L ∞ (R t ; M + (R n y × R n η )) and ν ∈ L ∞ (R t ; M + (R d x × R d ξ )), which solve the classical system (5.3) in the sense of distributions.
In other words, we obtain the same classical limiting system for ε = δ → 0 + , as in the iterated limit ε → 0 + and δ → 0 + . In summary, we have established the diagram of semi-classical limits as is shown in Figure 1. (It is very likely that the missing "upper right corner" within Fig. 1 can also be completed by using arguments similar to those given above.) Figure 1. The diagram of semi-classical limits: the iterated limit and the classical limit.
6. Numerical methods based on time-splitting spectral approximations.
In this section, we will develop efficient and accurate numerical methods for the semi-classically scaled TDSCF equations (2.3) and the Ehrenfest equations (4.9). The highly oscillatory nature of these models strongly suggests the use of spectral algorithms, which are the preferred method of choice when dealing with semi-classical models, cf. [15]. In the following, we will design and investigate time-splitting spectral algorithms, for both the TDSCF system and the Ehrenfest model, which will be shown to be second order in time. The latter is not trivial due to the self-consistent coupling within our equations and it will become clear that higher order methods can, in principle, be derived in a similar fashion. Furthermore, we will explore the optimal meshing strategy if only physical observables and not the wave function itself are being sought. In particular, we will show that one can take time steps independent of semi-classical parameters if one only aims to capture correct physical observables.
6.1. The SSP2 method for the TDSCF equations. In our numerical context, we will consider the semi-classically scaled TDSCF equations (2.3) in one spatial dimension and subject to periodic boundary conditions, i.e.
As before, we denote Υ ε,δ = ϕ ε,δ , V ϕ ε,δ L 2 y and Λ ε,δ = ψ ε,δ , h δ ψ ε,δ L 2 x . Clearly, a, b ∈ R have to be chosen such that the numerical domain [a, b] is sufficiently large in order to avoid the possible influence of boundary effects on our numerical solution. The numerical method developed below will work for all ε and δ, even if ε = o(1) or δ = o(1). In addition, we will see that it can be naturally extended to the multi-dimensional case. 6.1.1. The construction of the numerical method. We assume, on the computational domain [a, b], a uniform spatial grid in x and y respectively, x j1 = a + j 1 ∆x, y j2 = a + j 2 ∆y, where j m = 0, · · · N m − 1, N m = 2 nm , n m are some positive integers for m = 1, 2, and ∆x = b−a N1 , ∆y = b−a N2 . We also assume discrete time t k = k∆t, k = 0, · · · , K with a uniform time step ∆t.
The construction of our numerical method for (6.1) is based on the following operator splitting technique. For every time step t ∈ [t n , t n+1 ], we solve the kinetic step and the potential step possibly for some fractional time steps in a specific order. For example, if Strang's splitting is applied, then the operator splitting error is clearly second order in time (for any fixed value of ε). However, in the semi-classical regime ε → 0 + , a careful calculation shows that the operator splitting error is actually O(∆t 2 /ε), cf. [2,16]. Next, let U n j1 be the numerical approximation of the wave functions ψ ε,δ at x = x j1 and t = t n . Then, the kinetic step for ψ ε,δ can be solved exactly in Fourier space via: whereÛ n l1 are the Fourier coefficients of U n j1 , defined bŷ Similarly, the kinetic step for ϕ ε,δ can also be solved exactly in the Fourier space. On the other hand, for the potential step (6.3) with t 1 < t < t 2 , we formally find Λ ε,δ (y, s) ds ϕ ε,δ (y, t 1 ), (6.5) where 0 < t 2 − t 1 ≤ ∆t. The main problem here is, of course, that the mean field potentials Υ ε,δ and Λ ε,δ depend on the solution ψ ε,δ , ϕ ε,δ themselves. The key observation is, however, that within each potential step, the mean field potential Υ ε,δ is in fact time-independent (at least if we impose the assumption that the external potential V = V (x, y) does not explicitly depend on time). Indeed, a simple calculation shows

SHI JIN, CHRISTOF SPARBER AND ZHENNAN ZHOU
In other words, (6.4) simplifies to which is an exact solution formula for ψ ε,δ at t = t 2 .
The same argument does not work for the other self-consistent potential Λ ε,δ , since formally x . However, since Λ ε,δ (y, t) = ψ ε,δ , h δ ψ ε,δ L 2 x , the formula (6.6) for ψ ε,δ allows to evaluate Λ ε,δ (y, t) for any t 1 < t < t 2 . Moreover, the above expression for ∂ t Λ ε,δ , together with the Cauchy-Schwarz inequality and the energy estimate in Lemma 2.5, directly implies that ∂ t Λ ε,δ is O(1). Thus, one can use standard numerical integration methods to approximate the time-integral within (6.5). For example, one can use the trapezoidal rule to obtain Obviously, this approximation introduces a phase error of order O(∆t 2 /ε), which is comparable to the operator splitting error. This is the reason why we call the outlined numerical method SSP2, i.e., a second order Strang-spliting spectral method.
Remark 6.1. In order to obtain a higher order splitting method to the equations, one just needs to use a higher order quadrature rule to approximate the time-integral within (6.5).
6.1.2. Meshing strategy. In this subsection, we will analyze the dependence on the semi-classical parameters of the numerical error by applying the Wigner transformation onto the scheme we proposed above. In particular, this yields an estimate on the approximation error for (the expectation values of) physical observables due to (3.3). Our analysis thereby follows along the same lines as in Refs. [2,16]. For the sake of simplicity, we shall only consider the differences between their cases and ours. We denote the Wigner transforms W ε,δ ≡ W δ [ψ ε,δ ] and w ε,δ = w ε [ϕ ε,δ ], which satisfy the system   To this end, we are interested in analyzing two special cases: δ = O(1), and ε 1, or δ = ε 1. These correspond to the semi-classical limits we showed in Theorem 4.7 and Theorem 5.4.
We first consider the case δ = ε 1, where Wigner transformed TDSCF system reduces to (5.5). For convenience, we suppress the parameter δ, and write the potential step for ϕ ε as, and the associated Wigner transform w ε in the potential step satisfies, In view of (6.7), if we denote the approximation on the right hand side byφ ε , theñ ϕ ε is the exact solution to the following equation where G ε (y) = 1 2 (Λ ε (y, t 1 ) + Λ ε (y, t 2 )).
In summary, we conclude that for the SSP2 method, the approximation within the potential step results in an one-step error which is bounded by O(ε∆t + ∆t 3 ). Thus, for fixed ∆t, and as ε → 0 + , this one-step error in computing the physical observables is dominated by O(∆t 3 ) and we consequently can take ε-independent time steps for accurately computing semi-classical behavior of physical observables. By standard numerical analysis arguments, cf. [2,16], one consequently finds, that the SSP2 method introduces an O(∆t 2 ) error in computing the physical observables for ε 1 within an O(1) time interval. Similarly, one can obtain the same results when δ is fixed while ε 1. We remark that, if a higher order operator splitting is applied to the TDSCF equations, and if a higher order quadrature rule is applied to approximate formula (6.5), one obviously can expect higher order convergence in the physical observables.
6.2. The SVSP2 method for the Ehrenfest equations. In this section, we consider the one-dimensional Ehrenfest model obtained in Section 4.1. More precisely, we consider a (semi-classical) Schrödinger equation coupled with Hamilton's equations for a classical point particle, i.e with initial conditions ψ ε |t=0 = ψ ε in (x), y |t=0 = y 0 , η |t=0 = η 0 , and subject to periodic boundary conditions. Inspired by the SSP2 method, we shall present a numerical method to solve (6.10), which is second order in time and works for all 0 < δ 1.
As before, we assume a uniform spatial grid x j = a + j∆x, where N = 2 n0 , n 0 is an positive integer and ∆x = b−a N . We also assume uniform time steps t k = k∆t, k = 0, · · · , K for both the Schrödinger equation and Hamilton's ODEs. For every time step t ∈ [t n , t n+1 ], we split the system (6.10) into a kinetic step y = η,η = 0; (6.11) and a potential step We remark that, the operator splitting method for the Hamilton's equations may be one of the symplectic integrators. The readers may refer to [12] for a systematic discussion.
As before, the kinetic step (6.11) can be solved analytically. On the other hand, within the potential step (6.12), we see that ∂ t V (x, y(t)) = ∇ y V ·ẏ(t) = 0, (6.13) i.e., V (x, y(t)) is indeed time-independent. Moreover x . Now, we can use the first equation in (6.12) and the fact that V (x, y(t)) ∈ R to infer that the first two terms on the right hand side of this time-derivate cancel each other. We thus have in view of (6.13). In other words, also the semi-classical force is time-independent within each potential step. In summary, we find that for t ∈ [t 1 , t 2 ], the potential step admits the following exact solutions as well as This implies, that for this type of splitting method, there is no numerical error in time within the kinetic or the potential steps and thus, we only pick up an error of order O(∆t 2 /δ) in the wave function and an error of order O(∆t 2 ) in the classical coordinates induced by the operator splitting. Standard arguments, cf. [2,16], then imply that one can use δ-independent time steps to correctly capture the expectation values of physical observables. We call this proposed method SVSP2, i.e., a second order Strang-Verlet splitting spectral method. It is second order in time but can easily be improved by using higher order operator splitting methods for the Schrödinger equation and for Hamilton's equations.
7. Numerical tests. In this section, we test the SSP2 method for the TDSCF equations and the SVSP2 method for the Ehrenfest system. In particular, we want to test the methods after the formation of caustics, which generically appear in the WKB approximation of the Schrödinger wave functions, cf [24]. We also test the convergence properties in time and with respect to the spatial grids for the wave functions and the following physical observable densities i.e., the particle density and current densities associated to ψ ε (and analogously for ϕ ε ).
7.1. SSP2 method for the TDSCF equations. We first study the behavior of the proposed SSP2 method. In Example 1, we fix δ and test the SSP2 method for various ε. In Example 2 and Example 3, we take δ = ε and assume the same spatial grids in x and y.
Example 1. In this example, we fix δ = 1, and test the SSP2 method for various ε = o(1). We want to test the convergence in spatial grids and time, and whether ε-independent time steps can be taken to calculate accurate physical observables. Assume x, y ∈ [−π, π] and let V (x, y) = 1 2 (x + y) 2 . The initial conditions are of the WKB form, ψ δ in (x) = e −2(x+0.1) 2 e i sin x δ , ϕ ε in (y) = e −5(y−0.1) 2 e i cos y ε .
In the following all our numerical tests are computed till the stopping time T = 0.4.    These results show that, ε-independent time steps can be taken to obtain accurate physical observables, but not accurate wave functions. line represents the numerical solution with ε-independent ∆t, and the solid line represents the numerical solution with ε-dependent ∆t. From these figures, we observe the numerical convergence (in the weak sense) to the limit solutions after caustics formation, and that the numerical scheme can capture the physical observables with ε-independent ∆t.
Next, we come back to take the harmonic coupling potential V (x, y) = 1 2 (x + y) 2 , which ensures in a nontrivial coupling between the two sub-systems. We again want to test whether ε-independent ∆t can be taken to correctly capture the behavior of physical observables. We solve the TDSCF equations with resolved spatial grids, which means ∆x = O(ε). The numerical solutions with ∆t = O(ε) are used as the reference solutions. For ε = 1 256 , 1 512 , 1 1024 , 1 2048 , 1 4096 , we fix ∆t = 0.005, and compute till T = 0.54. The l 2 norm of the error for the wave functions and the error for the position densities is calculated. We see in Figure 7 that the former increases as ε → 0 + , but the error in the physical observables does not change noticeably.
Example 3. In this example, we want to test the convergence in the spatial grid ∆x and in the time step ∆t as ε = δ → 0 + . According to the previous analysis, the spatial oscillations of wavelength O(ε) need to be resolved. On the other hand, if the time oscillation with frequency O(1/ε) is resolved, one gets accurate approximation even for the wave functions itself (not only quadratic quantities of it). Unresolved time steps of order O(1) can still give correct physical observable densities. More specifically, one expects second order convergence with respect to time in both wave functions (and in the physical observables), and spectral convergence in the respective spatial variable.  , ∆x = ε 8 , respectively. The reference solution is computed with the same ∆x, but ∆t = 0.54ε 4 . These results show that, ε-independent time steps can be taken to obtain accurate physical observables, but not accurate wave functions. To test the spatial convergence, we take ε = 1 256 , and the reference solution is computed by well resolved mesh ∆x = 2πε 64 , ∆t = 0.4ε 16 until T = 0.4. Then, we compute with the same time step, but with different spatial grids. The results are illustrated in Figure 8. We observe that, when ∆x = O(ε), the error decays quickly to be negligibly small as ∆x decreases. However, when the spatial grids do not well resolve the ε-scale, the method would actually give solutions with O(1) error.
At last, to test the convergence in time, we take ε = 1 1024 , and the reference solution is computed through a well resolved mesh with ∆x = 2πε 16 , ∆t = 0.4 8192 till T = 0.4. Then, we compute with the same spatial grids, but with different time steps. The results are illustrated in Figure 9. We observe that the method is stable even if ∆t ε. Moreover, we get second order convergence in the wave functions as well as in the physical observable densities. 7.2. SVSP2 method for the Ehrenfest equations. Now we solve the Ehrenfest equations (6.10) by the SVSP2 method. Assume x ∈ [−π, π], and assume periodic boundary conditions for the electronic wave equation.
Example 4. In this example, we want to test if δ-independent time steps can be taken to capture correct physical observables and the convergence in the time step which is expected to be of the second order. The potential is again V (x, y) = These results show that, when δ = ε 1, the SSP2 method is unconditionally stable and is second order accurate in time.
First, we test whether δ-independent ∆t can be taken to capture the correct physical observables. We solve the equations with resolved spatial grids, which means ∆x = O(δ). The numerical solutions with ∆t = O(δ) are used as the reference solutions. For δ = 1/256, 1/512, 1/1024, 1/2048, 1/4096, we fix ∆t = 0.4 64 , and compute until T = 0.4. The l 2 norm of the error in wave functions, the error in position densities, and the error in the coordinates of the nucleus are calculated. We see in Figure 10 that the error in the wave functions increases as δ → 0 + , but the errors in physical observables and in the classical coordinates do not change notably.
Next, we test the convergences with respect to the time step in the wave function, the physical observables and the classical coordinates. We take δ = 1 1024 , and the reference solution is computed by well resolved mesh ∆x = 2πε 16 , ∆t = 0.4 8192 till T = 0.4. Then, we compute with the same spatial grids, but with difference time steps. The results are illustrated in Figure 11. We observe that, the method is stable even if ∆t ε, and clearly, we get second order convergence in the wave functions, the physical observable densities and the classical coordinates.   These results show that, the SVSP2 method is unconditionally stable and is second order accurate in time.