An Efficient Time-splitting Method for the Ehrenfest Dynamics

TL;DR

This paper introduces a stable and efficient time-splitting numerical method for Ehrenfest dynamics, enabling accurate quantum-classical simulations with larger time steps independent of the semiclassical parameter.

Contribution

A novel time-splitting scheme for Ehrenfest dynamics that is stable, preserves the semiclassical limit, and allows larger time steps for efficient simulations.

Findings

Scheme is stable uniformly with respect to the semiclassical parameter.

Method preserves the discrete semiclassical limit.

Numerical examples confirm the effectiveness of the meshing strategy.

Abstract

The Ehrenfest dynamics, representing a quantum-classical mean-field type coupling, is a widely used approximation in quantum molecular dynamics. In this paper, we propose a time-splitting method for an Ehrenfest dynamics, in the form of a nonlinearly coupled Schr\"odinger-Liouville system. We prove that our splitting scheme is stable uniformly with respect to the semiclassical parameter, and, moreover, that it preserves a discrete semiclassical limit. Thus one can accurately compute physical observables using time steps induced only by the classical Liouville equation, i.e., independent of the small semiclassical parameter - in addition to classical mesh sizes for the Liouville equation. Numerical examples illustrate the validity of our meshing strategy.

Tables

Table 1. Table 1: The relative ℓ 2 − limit-from superscript ℓ 2 \ell^{2}- difference (defined as ‖ μ 1 − μ 2 ‖ ℓ 2 ‖ μ 2 ‖ ℓ 2 subscript norm subscript 𝜇 1 subscript 𝜇 2 superscript ℓ 2 subscript norm subscript 𝜇 2 superscript ℓ 2 \frac{\left\|\mu_{1}-\mu_{2}\right\|_{\ell^{2}}}{\left\|\mu_{2}\right\|_{\ell^{2}}} ) for various h ℎ h .

$h$	$1/256$	$1/1024$	$1/4096$
$\frac{{\Vert {\mu }_1-{\mu }_2\Vert }_{{\ell }^2}}{{\Vert {\mu }_2\Vert }_{{\ell }^2}}$	1.65e-03	1.69e-03	1.70e-03

Equations

$$\{\begin{array}{rl}ih{\partial }_t\psi & =-\frac{h^2}{2}{\Delta }_x\psi +V\left(x,y(t)\right)\psi ,\psi (0,x)={\psi }_{\mathrm{in}}(x),\\ \dot{y}(t)& =\eta (t),y(0)=y_0,\\ \dot{\eta }(t)& =-{\nabla }_y{V}_{\mathrm{E}}(y(t)),\eta (0)={\eta }_0.\end{array}$$

$$V_{\mathrm{E}}(y,t)={\int }_{{\mathbb{R}}^d}V\left(x,y\right){\left|\psi (x,t)\right|}^{2}dx.$$

$$\ddot{y}(t)=-{\nabla }_y{V}_{\mathrm{E}}(y(t)),$$

$$\{\begin{array}{rl}& ih{\partial }_t{\psi }^h=-\frac{h^2}{2}{\Delta }_x{\psi }^h+{{\rm Y}}^h\left(x,t\right){\psi }^h,{\psi }^h(0,x)={\psi }_{\mathrm{in}}^h(x)\\ & {\partial }_t{\mu }^h+\eta \cdot {\nabla }_y{\mu }^h+F^h\left(y,t\right)\cdot {\nabla }_{\eta }{\mu }^h=0,{\mu }^h(0,x,\eta )={\mu }_{\mathrm{in}}(y,\eta ).\end{array}$$

$${{\rm Y}}^h\left(x,t\right)$$

$$\mu (t,y,\eta )=\delta (y-y(t),\eta -\eta (t)).$$

$$V\in C_0^2\left({\mathbb{R}}_x^d\times {\mathbb{R}}_y^n\right)\text{and }V\left(x,y\right)\ge 0,\forall \left(x,y\right)\in {\mathbb{R}}_x^d\times {\mathbb{R}}_y^n,$$

$$J=\frac{d-c}{\Delta y},K=\frac{\beta -\alpha }{\Delta \eta },M=\frac{b-a}{\Delta x},y_j=c+j\Delta y,{\eta }_k=\alpha +k\Delta \eta ,x_j=a+j\Delta x.$$

$$\{\begin{array}{rl}& ih{\partial }_t{\psi }^h=-\frac{h^2}{2}{\Delta }_x{\psi }^h,\\ & {\partial }_t{\mu }^h=-\eta \cdot {\nabla }_y{\mu }^h-F^h\left(y,t\right)\cdot {\nabla }_{\eta }{\mu }^h,\end{array}$$

$$\{\begin{array}{rl}& ih{\partial }_t{\psi }^h={{\rm Y}}^h\left(x,t\right){\psi }^h,\\ & {\partial }_t{\mu }^h=0,\end{array}$$

$${\psi }_j^{h,n+1}=\exp \left(-\frac{i}{h}{{\rm Y}}^h\left(x_j,t_{\ast }\right)\Delta t\right){\psi }_j^{h,\ast }.$$

$$\{\begin{array}{rl}& ih{\partial }_t{\psi }^h=-\frac{h^2}{2}{\partial }_{xx}{\psi }^h,\\ & \frac{d}{dt}{\mu }_{jk}^h=-{\eta }_k{\left(D_y{\mu }^h\right)}_{jk}-F_j^h{\left(D_{\eta }{\mu }^h\right)}_{jk},\end{array}$$

$$\{\begin{array}{rl}& {\psi }_j^{h,\ast }=\frac{1}{M}\sum _{\ell =-M/2}^{M/2-1}{e}^{-ih{\omega }_{\ell }^2/2}{\hat{\psi }}_{\ell }^{h,n}{e}^{i{\omega }_{\ell }\left(x_j-a\right)},j=0,\dots ,M-1,\\ & \frac{{\mu }_{jk}^{h,\ast }-{\mu }_{jk}^{h,n}}{\Delta t}=-{\eta }_k{\left(D_y{\mu }^{h,n}\right)}_{jk}-F_j^{h,n}{\left(D_{\eta }{\mu }^{h,n}\right)}_{jk},\end{array}$$

$$\{\begin{array}{rl}& {\eta }_k{\left(D_y{\mu }^{h,n}\right)}_{jk}=\frac{1}{2}\left({\eta }_k+\left|{\eta }_k\right|\right)\frac{{\mu }_{jk}^{h,n}-{\mu }_{j-1,k}^{h,n}}{\Delta y}+\frac{1}{2}({\eta }_k-\left|{\eta }_k\right|)\frac{{\mu }_{j+1,k}^{h,n}-{\mu }_{j,k}^{h,n}}{\Delta y},\\ & F_j^{h,n}{\left(D_{\eta }{\mu }^{h,n}\right)}_{jk}=\frac{1}{2}(F_j^{h,n}+\left|F_j^{h,n}\right|)\frac{{\mu }_{jk}^{h,n}-{\mu }_{j,k-1}^h}{\Delta \eta }+\frac{1}{2}(F_j^{h,n}-\left|F_j^{h,n}\right|)\frac{{\mu }_{j,k+1}^{h,n}-{\mu }_{jk}^{h,n}}{\Delta \eta }.\end{array}$$

$$\{\begin{array}{rl}& ih{\partial }_t{\psi }^h={{\rm Y}}_d^h\left(x,t\right){\psi }^h,\\ & \frac{d}{dt}{\mu }_{jk}^h=0,\end{array}$$

$$\begin{array}{r}{\psi }_j^{h,n+1}=\exp \left(-i{{\rm Y}}_d^{h,\ast }\left(x_j\right)\Delta t/h\right){\psi }_j^{h,\ast },{\mu }_{jk}^{h,n+1}={\mu }_{jk}^{h,\ast }\end{array}$$

$${{\rm Y}}_d^{h,\ast }\left(x\right)=\sum _{j=0}^{J-1}\sum _{k=0}^{K-1}V\left(x,y_j\right){\mu }_{jk}^{h,\ast }\Delta y\Delta \eta =\sum _{j=0}^{J-1}\sum _{k=0}^{K-1}V\left(x,y_j\right){\mu }_{jk}^{h,n+1}\Delta y\Delta \eta ,$$

$${\Vert \varphi \Vert }_{L^2}={\left({\int }_a^b{\left|\varphi \left(x\right)\right|}^{2}dx\right)}^{1/2},{\Vert \psi \Vert }_{{\ell }^2}={\left(\frac{b-a}{M}\sum _{j=0}^{M-1}{\left|{\psi }_j\right|}^2\right)}^{1/2}.$$

$${\Vert f_{\mathrm{I}}\Vert }_{L^2}^2={\Vert f\Vert }_{{\ell }^2}^2=\frac{b-a}{M}\sum _{j=0}^{M-1}{\left|f\left(x_j\right)\right|}^2$$

$$f_{\mathrm{I}}\left(x\right)=\frac{1}{M}\sum _{\hat{ȷ}=-\frac{M}{2}}^{\frac{M}{2}-1}{\hat{f}}_{\hat{ȷ}}{e}^{i{\omega }_{\hat{ȷ}}\left(x-a\right)},{\omega }_{\hat{ȷ}}=\frac{2\pi \hat{ȷ}}{b-a}.$$

$${\Vert {\psi }^{h,n}\Vert }_{{\ell }^2}={\Vert {\psi }^{h,0}\Vert }_{{\ell }^2},{\Vert {\psi }_{\mathrm{I}}^{h,n}\Vert }_{L^2}={\Vert {\psi }_{\mathrm{I}}^{h,0}\Vert }_{L^2},\text{ for n=1,2,…,}$$

$$\sum _{j=0}^{J-1}\sum _{k=0}^{K-1}{\mu }_{jk}^{h,n}=\sum _{j=0}^{J-1}\sum _{k=0}^{K-1}{\mu }_{jk}^{h,0}.$$

$$\frac{1}{b-a}{\Vert {\psi }^{h,n+1}\Vert }_{{\ell }^2}^2$$

$$\frac{1}{M}\sum\limits_{j=0}^{M-1}e^{i\left(\omega_{q}-\omega_{p}\right)\left(x_{j}-a\right)}=\frac{1}{M}\sum\limits_{j=0}^{M-1}e^{i2\pi\left(q-p\right)j/M}=\left\{\begin{array}[]{cc}0,&q-p\neq mM\\ 1,&q-p=mM\end{array}\right.,\ m\in\mathbb{Z}.$$

$$\frac{1}{b-a}{\Vert {\psi }^{h,n+1}\Vert }_{{\ell }^2}^2$$

$$\frac{1}{M}\sum\limits_{\hat{\jmath}=-\frac{M}{2}}^{\frac{M}{2}-1}e^{i\omega_{\hat{\jmath}}\left(x_{q}-x_{p}\right)}=\frac{1}{M}\sum\limits_{\hat{\jmath}=-\frac{M}{2}}^{\frac{M}{2}-1}e^{i2\pi\left(q-p\right)\hat{\jmath}/M}=\left\{\begin{array}[]{cc}0,&q-p\neq mM\\ 1,&q-p=mM\end{array}\right.,\ m\in\mathbb{Z}.$$

$$\underset{k}{\max }|{\eta }_k|\frac{\Delta t}{\Delta y}+{\Vert {\partial }_{y}V\Vert }_{L^{\infty }}\frac{\Delta t}{\Delta \eta }\le 1.$$

$$w^h[f^h]\left(x,\xi \right)=\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{f}^h\left(x-\frac{h}{2}y\right)\overline{f^h}\left(x+\frac{h}{2}y\right)e^{i\xi \cdot y}dy.$$

$${\Vert w^h[f^h]\Vert }_{L^2\left({\mathbb{R}}^{2d}\right)}=\frac{1}{{\left(2\pi \right)}^{\frac{d}{2}}{h}^d}{\Vert f^h\Vert }_{L^2\left({\mathbb{R}}^d\right)}^2.$$

Full text

1 Introduction

Ab initio methods have played a fundamental role in the numerical simulation of large quantum systems, in particular in quantum molecular dynamics. Different from classical approaches based on pre-defined potentials, the underlying idea of ab initio molecular dynamics is to compute the forces acting on the nuclei as a feedback of the electronic structures. This procedure is also known as the “on-the-fly” calculation in the chemistry literature (for detailed reviews, see, e.g., [4, 23, 22, 28]). One of the most widely used of these methods is the so-called Ehrenfest dynamics, a mean-field treatment named in honor of Paul Ehrenfest who was among the first to address the problem of how to derive classical dynamics from the underlying quantum mechanical equations [10]. His idea is to separate the whole system into two parts: a fast varying, quantum mechanical part (for, say, electrons) and a slowly varying part (for the much heavier nuclei) in which one can pass to the (semi-)classical limit. In quantum chemistry, this is usually possible by taking advantage of the large mass difference between electrons and nuclei.

Typically, the Ehrenfest molecular dynamics refers to a Schrödinger equation, coupled with a classical Newtonian flow, cf. [6, 9, 24, 26, 3, 28]. The simplest such model reads

[TABLE]

Here, we denote by $0<h\ll 1$ a dimensionless rescaled Planck’s constant, and by $\psi=\psi\left(x,t\right)$ with $x\in\mathbb{R}^{d}$, $t\in\mathbb{R}_{+}$, the wave function of the fast, quantum mechanical degrees of freedom, which is assumed to be normalized such that $\|{\psi}(\cdot,t)\|_{L^{2}}=1$ for all $t\geq 0$. In addition, the slow degrees of freedom are described, for any time $t\in\mathbb{R}_{+}$, by their classical position $y(t)\in\mathbb{R}^{n}$ and momentum $\eta(t)\in\mathbb{R}^{n}$. We thereby allow for $n\in\mathbb{N}$ and $d\in\mathbb{N}$ to be not necessarily equal, depending on the physical application. Finally, for a given coupling potential $V=V(x,y)\in\mathbb{R}$, the force describing the back-reaction of the quantum part onto the slow degrees of freedom is given by the gradient in $y\in\mathbb{R}^{n}$ of the so-called Ehrenfest potential

[TABLE]

Clearly, one obtains a version of Newton’s second law for $y(t)$ by eliminating the momentum variable $\eta(t)$ and writing

[TABLE]

instead of the first order Hamiltonian system above.

Regarding the derivation of Ehrenfest dynamics, the majority of literature available today invokes WKB asymptotics for the slow degrees of freedom, leading to a Hamilton-Jacobi equation which suffers from the appearance of caustics, see, e.g., [6, 24]. To circumvent this problem and derive a semiclassical limit which is valid globally in time, a, by now classical, tool is the Wigner transform [29]. The latter gives rise to a Liouville equation for the associated semi-classical phase-space measure (or, Wigner measure) which “unfolds the caustics”, see [11, 19, 21, 25]. In the context of Ehrenfest dynamics, such an analysis was carried out in [15]. In there, the authors start from a system of time-dependent, self-consistent field equations, motivated by [5, 16, 17, 20], and derive (among other things) the following mixed quantum-classical system:

[TABLE]

Here, $\mu^{h}(\cdot,\cdot,t)\in\mathcal{M}^{+}(\mathbb{R}^{n}_{y}\times\mathbb{R}^{n}_{\eta})$ denotes the phase-space probability density for the slowly varying degrees of freedom at time $t$, $F^{h}=-\nabla_{y}V_{\rm E}$, i.e., the force obtained from the Ehrenfest potential, and

[TABLE]

We call this system the Schrödinger-Liouville-Ehrenfest (SLE) System and from now on represent the dependence on the small semi-classical parameter $h>0$ by superscripts. Note that the dependence of $\mu^{h}$ on $h$ stems purely from the forcing through the Ehrenfest potential appearing in the Liouville equation. The latter is an Eulerian description of the classical Hamiltonian flow. In particular, one formally obtains (1.1), from (1.2), in the case where $\mu$ corresponds to a single particle distribution concentrated on the classical trajectories $(y(t),\eta(t))$, i.e.,

[TABLE]

Such kind of Wigner measures can be obtained as the classical limit of a particular type of wave functions, called semi-classical wave packets, or coherent states, see [19].

Given the dispersive nature of Schrödinger’s equation, the main numerical difficulty for $h\ll 1$ is that one needs to resolve oscillations of frequency of order $\mathcal{O}(1/h)$ in both time and space, as they are present in the solution $\psi^{h}$, see [14] for a broad review of this problem. Naively, this requires one to use time-steps of order $\Delta t=o(h)$ as well as a spatial grid with $\Delta x=o(h)$. However, it was proved in [1], using a Wigner measure analysis, that for a single linear Schrödinger equation, a time-splitting spectral method can still correctly capture physical observables, i.e., real-valued quadratic quantities in $\psi^{h}$, even for time-steps much larger than $h$. Thus one only needs to resolve the high frequency oscillations spatially, which is a huge numerical advantage. For nonlinear Schrödinger equations, in general, this is no longer true, as was numerically demonstrated in [2]. The SLE system (1.2) is a nonlinearly coupled system, and one therefore expects the same type of problem at first glance. Nevertheless, we shall in the following develop an efficient numerical method for the SLE system which allows large (compared with $h$) computational mesh-sizes in both $y$ and $\eta$ and a large time step for both the Schrödinger and the Liouville equations, while still correctly capturing physical observables. While large meshes in $y$ and $\eta$ do not seem so surprising, the possibility of large time steps for solving the Schrödinger equation is far from obvious, due to the nonlinear nature of the SLE system.

Our numerical algorithm is inspired by, but different from the time-splitting method used in [15]. Based on a spectral method for the Schrödinger equation and an upwind scheme for the Liouville part of (1.2), we shall first prove stability for our algorithm, uniformly in $h$. Furthermore, by utilizing the Wigner analysis developed in [15] and adopting it to our particular setting, we shall also prove that physical observables (which can be characterized by the moments of the Wigner distribution), are captured correctly even if $\Delta y$, $\Delta\eta$ and $\Delta t$, i.e., the time step for the entire SLE system, are $\mathcal{O}(1)$ and thus independent of $h$. To this end, we follow the strategy of [1], and prove that the semi-discretized SLE-system, with $\Delta y$, $\Delta\eta$, $\Delta t$ fixed, converges to the correct semiclassical limiting system, as $h\to 0$. In this analysis we shall, for simplicity, consider $x$ to be continuous, since, as already stated above, $\Delta x\to 0$, as $h\to 0$, even for a single linear Schrödinger equation. In summary, our scheme can be seen to be asymptotic-preserving in $t$, $y$, and $\eta$, which is a well-established numerical concept for multi-scale kinetic equations, cf. [12, 13]. To our knowledge, this is the first work that proves the existence of a global in-time $h$-independent meshing strategy for physical observables associated to a nonlinear Schrödinger-type system.

In this context, we note that the authors of [7, 8] study a time-splitting scheme for nonlinear Schrödinger equations with cubic nonlinearity. Using a WKB type representation of the solution, they are able to prove a similar asymptotic preserving property. However, the main drawback of their method is, that it is only valid before the formation of caustics in the Hamilton-Jacobi equation for the WKB phase function. The system studied in the present paper has a weaker nonlinear structure which allows the use of Wigner transformation techniques which are valid for all time $t$.

The rest of this paper is now organized as follows: In Section 2 we present the time-splitting method for the SLE system and briefly discuss some of the inherent numerical difficulties. The stability, uniformly in $h$, is then proved in Section 2.2 for the fully discretized system. In Section 3, we shall give a brief review of Wigner transformation methods and the classical limit of the SLE system. The spatial meshing strategy announced above is then studied in Section 4 by deriving the classical limit of a semi-discrete SLE system. In Section 5 we focus on the time-discretization and prove that our scheme allows for time-steps independent of $h$. Finally, Section 6 presents some numerical examples illustrating our analytical results.

2 The time-splitting scheme and its basic properties

We shall, from now on, consider the Schrödinger-Liouville-Ehrenfest (SLE) System (1.2) with the following assumption on the coupling potential $V$:

[TABLE]

where $C^{2}_{0}$ denotes the set of twice continuously differentiable functions which vanish at infinity together with all their derivatives.

Remark 2.1**.**

This is the same assumption as in [15], where it is used to furnish a rigorous Wigner analysis of the self-consistent field equation. Note that, in particular, it implies $V\in W^{2,\infty}\left(\mathbb{R}_{x}^{d}\times\mathbb{R}_{y}^{n}\right)$. It is conceivable that the regularity requirement and the decay at infinity can be lowered at the expense of more technicalities. The assumption $V(x,y)\geq 0$ is in fact not very restrictive for the potentials bounded from below. It corresponds to a proper choice of the zero point of the potential axis.

We aim for an algorithm which fully utilizes the quantum-classical coupling. Thus, while it makes sense to use a finer (i.e., smaller than $h$) spatial discretization in $x$ to solve the Schrödinger equation, we want to use much larger (than $h$) meshes in $y$ and $\eta$ when solving the Liouville equation. That this is indeed possible is not obvious, since the potential $\Upsilon^{h}\left(x,t\right)$ appearing in the Schrödinger equation is time-dependent, and moreover nonlinearly coupled to the Liouville equation (hence it inherits the computational error obtained from discretizing in $y$ and $\eta$).

Remark 2.2**.**

In our discussion, we will only consider compactly supported initial data $\psi^{h}_{\rm in}$, $\mu_{\rm in}$, in order to simulate the SLE system problem based on an infinitely large spatial domain within a sufficiently large, but finite box with periodic boundary conditions.

2.1 A new time-splitting scheme for the SLE system

In order to describe our scheme, we henceforth assume that we are given a sufficiently small $\Delta x\sim\mathcal{O}(h)$, used to solve the quantum mechanical part of (1.2), while the larger grid meshes $\Delta y,\Delta\eta\sim\mathcal{O}\left(1\right)$ are applied for the classical part. With this in mind, let

[TABLE]

The time-splitting spectral scheme can then described as follows: From time $t=t_{n}=n\Delta t$ to $t=t_{n+1}=\left(n+1\right)\Delta t$, with $\Delta t$ given, the SLE system is solved in two steps. First, solve

[TABLE]

from $t=t_{n}$ to an intermediate time $t_{\ast}$. Then, solve

[TABLE]

with initial data obtained from Step 1, to obtain the solution at time $t=t_{n+1}$.

In (2.1), the Schrödinger equation will be discretized in space by a spectral method and integrated in time exactly using a Fast Fourier Transform. The Liouville equation can be solved either by a spectral method, or by a finite difference (e.g., upwind) scheme in space, and then marching the corresponding ODE system forward in time. An advantage of our splitting method is that in the second step, $\Upsilon^{h}\left(x,t\right)$ defined in (1.3) is indeed independent of time, since obviously $\mu^{h}$ is. In view of this, the time integration in (2.2) can also be solved exactly, which yields

[TABLE]

For the convenience of our later discussions, we shall now state our numerical scheme using an upwind spatial discretization of $\mu$ in more detail: The problem is solved in one spatial dimension $d=n=1$ from time $t=t_{n}$ to time $t=t_{n+1}$ using the following two steps:

In the first step, we solve

[TABLE]

where both $D_{y}\mu^{h}$ and $D_{\eta}\mu^{h}$ represent the numerical derivatives in our algorithm, which are treated using a standard conservative (for example, the upwind type) discretization. To solve the Liouville equation we shall we apply a forward-in-time Euler scheme for the time discretization. Explicitly, we thus have

[TABLE]

where ${w_{\ell}=\frac{2\pi\ell}{b-1}}$ and, for the upwind spatial discretization,

[TABLE]

The second step is then given by

[TABLE]

where $\Upsilon_{d}^{h}\left(x,t\right)$ is the quadrature approximation of $\Upsilon^{h}\left(x,t\right)$. Thus, we explicitly have

[TABLE]

where

[TABLE]

which can be viewed as a trapezoidal rule for $\mu$ with compact support, cf. Remark 2.2.

Remark 2.3**.**

It is straightforward to obtain an algorithm second order in time using the Strang splitting, which is omitted here.

2.2 Conservation property and stability of the scheme

We shall now prove the stability of the scheme given by (2.4) and (2.6). To this end, let $\mathbf{\psi=}\left(\psi_{0},\ldots,\psi_{M-1}\right)^{T}$. Let $\left\|\cdot\right\|_{L^{2}}$ and $\left\|\cdot\right\|_{\ell^{2}}$ be the usual $L^{2}$ and $\ell^{2}$ norm on the interval $\left(a,b\right)$ respectively, i.e.

[TABLE]

Notice that, for any periodic function $f$, the equality

[TABLE]

holds, where $f_{\rm I}$ denotes the trigonometric interpolant of $f$ on $\left\{x_{0},x_{1},\ldots,x_{M}\right\}$, i.e.

[TABLE]

Using this we can prove the following theorem.

Theorem 2.4**.**

The time-splitting spectral scheme conserves the mass. More precisely, it holds

[TABLE]

where, as before, $\psi_{\rm I}^{h,n}$ denotes the trigonometric interpolant of $\psi^{h,n}$. In addition,

[TABLE]

Proof.

First note that the last identity for $\mu_{jk}^{h,n}$ is a straightforward consequence of the fact that the discretized derivatives $D_{y}\mu$ and $D_{\eta}\mu$ are conservative.

It suffices to prove the first identity stated above due to (2.8). Noting our numerical algorithm (2.4), (2.6) and the definition of the norms (2.7), one computes

[TABLE]

Changing the order of summation, this is equal to

[TABLE]

where the second equality of the above follows from the fact that

[TABLE]

Similarly, by changing the order of summation again, we arrive at

[TABLE]

where the following identity has been used

[TABLE]

∎

Remark 2.5**.**

Theorem 2.4 implies that the scheme is stable uniformly in $h$, provided the positivity of $\mu$ under the following CFL condition, cf. [18]:

[TABLE]

3 Classical limit of the SLE system

As a preparatory step to the discussion of Section 4, we will now briefly review the results of [15] concerning the classical limit (via Wigner transforms) of the SLE system as $h\to 0$.

3.1 Wigner transform and Wigner measure

Let us first recall that $h$-scaled Wigner transform associated to any continuously parametrized family $f^{h}\equiv\{f^{h}\}_{0\leq h\leq 1}\in L^{2}\left(\mathbb{R}^{d}\right)$ is given by, cf. [11, 19, 21, 25]:

[TABLE]

By Plancherel’s Theorem and a change of variables one easily finds

[TABLE]

The real-valued function $w^{h}(x,\xi)$ acts as a quantum mechanical analogue for classical phase-space distributions. However, $w^{h}(x,\xi)\not\geq 0$ in general.

It has been proved in [19], that if the family of functions $f^{h}=\{f^{h}\}_{0\leq h\leq 1}$ is uniformly bounded in $L^{2}({\mathbb{R}}^{d})$ as $h\to 0_{+}$, i.e., if

[TABLE]

then the set of Wigner functions $\{w^{h}\}_{0<h\leq 1}$ is uniformly bounded in $\mathcal{A}^{\prime}$. The latter is the dual of the following Banach space

[TABLE]

where $C_{0}({\mathbb{R}}^{d})$ denotes the space of continuous functions vanishing at infinity and $\mathcal{F}_{\xi}$ denotes the Fourier transform with respect to the velocity $\xi$, only. More precisely, one finds that for any test function $\chi\in\mathcal{A}({\mathbb{R}}^{d}_{x}\times{\mathbb{R}}^{d}_{\xi})$,

[TABLE]

uniformly in $h$. Thus, up to extraction of sub-sequences $\{h_{n}\}_{n\in{\mathbb{N}}}$, with $h_{n}\to 0_{+}$ as $n\to\infty$, there exists a limiting object $w^{0}\equiv w\in\mathcal{A}^{\prime}({\mathbb{R}}^{d}_{x}\times{\mathbb{R}}^{d}_{\xi})$ such that

[TABLE]

It turns out that the limit is in fact a non-negative, bounded Borel measure on phase-space $w\in\mathcal{M}^{+}({\mathbb{R}}^{d}_{x}\times{\mathbb{R}}^{d}_{p})$, called the Wigner measure of $f^{h}$.

3.2 The Classical limit of the SLE system

Let $\psi^{h}$ and $\mu^{h}$ be the solution of the SLE system (1.2) and denote the Wigner function of $\psi^{h}\left(x,t\right)$ by

[TABLE]

A straightforward computation shows that the position density associated to $\psi^{h}\in L^{2}({\mathbb{R}}^{d})$ can be computed via

[TABLE]

where we recall, that due to our normalization,

[TABLE]

Moreover, by taking higher order moments in $\xi$ one (formally) finds the current density

[TABLE]

and the kinetic energy density

[TABLE]

Remark 3.1**.**

In order to make these computations rigorous, the integrals on the r.h.s. have to be understood in an appropriate sense, since $w^{h}\not\in L^{1}({\mathbb{R}}^{m}_{x}\times{\mathbb{R}}^{m}_{\xi})$ in general, see [19] for more details.

After Wigner transforming the Schrödinger equation, one finds that $w^{h}\left(x,\xi,t\right)$ satisfies the following nonlocal kinetic equation (see, e.g., [19]):

[TABLE]

where $w_{\rm in}^{h}\equiv w^{h}[\psi^{h}_{\rm in}]$ and

[TABLE]

with $\widehat{w}^{h}$ denoting the Fourier transformation of $w^{h}$ w.r.t. the second variable only.

Now, in order to utilize the weak-$\ast$ compactness properties of the Wigner function, we shall impose from now on that the initial mass and the initial kinetic energy are uniformly bounded with respect to $h$, i.e.,

[TABLE]

Remark 3.2**.**

In other words, we assume that

[TABLE]

This assumption is easily satisfied by initial data of WKB type, or by semi-classical wave packets.

It is proved in [15] that these uniform bounds on the initial mass and kinetic energy are propagated by the SLE system (1.2), which in turn implies that for all times $t\in{\mathbb{R}}_{+}$, the wave function $\psi^{h}(\cdot,t)$ is:

uniformly bounded in $L^{2}\left(\mathbb{R}^{d}\right)$ as $h\rightarrow 0^{+}$, i.e.

[TABLE] 2. 2.

$h$-oscillatory, i.e.

[TABLE]

where $C_{1}$ and $C_{2}$ are some constants independent of $0<h<1$.

In particular this implies the existence of a limiting Wigner measure $\nu(\cdot,\cdot,t)\in\mathcal{M}^{+}(\mathbb{R}_{x}^{d}\times\mathbb{R}_{\xi}^{d})$, such that for all $T>0$

[TABLE]

up to the extraction of subsequences. Moreover, on the same time-interval, one has

[TABLE]

Under our assumption (A1) on $V$, this can be used to prove that (see [15] for more details):

[TABLE]

uniformly on compact intervals in $y$ and $t$.

Similarly, one can pass to the limit $h\to 0_{+}$ in the equation for $\mu^{h}$ to find that there exists a limiting measure $\mu^{0}\equiv\mu\in\mathcal{M}^{+}({\mathbb{R}}^{n}_{y}\times{\mathbb{R}}^{n}_{\eta})$ which consequently solves (in the sense of distributions):

[TABLE]

Moreover, one can prove that

[TABLE]

In view of the definition (3.1), one also finds that

[TABLE]

and thus, the Wigner measure associated to $\psi^{h}$ satisfies the following Liouville equation (in the sense of distributions):

[TABLE]

In summary, one finds a system of two coupled Liouville equations (3.2)–(3.3) in the classical limit (we refer to [15] for a rigorous proof and further details).

4 The spatial meshing strategy

4.1 The semi-discretized SLE system and its energy

The analysis in this section will focus on the spatial meshing strategy. In order to show that it is possible to use a grid with $\Delta y,\Delta\eta\sim O(1)$, and thus, independent of $h$, we will consider a semi-discretized version of the SLE system (1.2) in one spatial dimension $d=n=1$ where the Liouville is discretized using an upwind scheme:

[TABLE]

Here, $\Upsilon_{d}^{h}\left(x,t\right)$ stands for the trapezoidal quadrature approximation of $\Upsilon^{h}\left(x,t\right)$, as before, whereas $F^{h}\left(y,t\right)$ includes exact derivative of the known function $V\left(x,y\right)$. We shall refer to (4.1) as the semi-discretized SLE system (s-SLE) and show that it yields the “correct” classical limit, i.e., the semi-discretized version of (3.2)–(3.3).

Before doing so, we will need to prove an a-priori estimate and the energy associated to (4.1). To this end, we define the semi-discrete energy as

[TABLE]

Here, and in the following, because of the periodicity of $\mu$, we shall use a cyclic index for $\mu_{jk}$, such that $\mu_{jk}=\mu_{j+J,k}=\mu_{j-J,k}=\mu_{j,k+K}=\mu_{j,k-K}$.

Theorem 4.1**.**

Under the assumptions (A1) and (A2), the energy $E_{d}(t)$ is bounded by a constant independent of $h$ for all $t\geq 0$.

Proof.

We start by showing that the initial energy is bounded. This is easily seen from

[TABLE]

where the first two integrals are clearly bounded by assumptions (A1) and (A2) and the last term is just a quadrature approximation of $\iint\mu_{\rm in}\,dy\,d\eta=1$, and hence bounded.

Next, we compute the time-derivative of $E_{d}$ as

[TABLE]

where

[TABLE]

First, a straightforward calculation shows $($I$)+($II$)=0$, since

[TABLE]

For simplicity we will, from now on, denote

[TABLE]

as well as

[TABLE]

A key observation is that $G^{h}$ is in fact Lipschitz with a Lipschitz constant $L>0$ independent of $h$, since

[TABLE]

for some $\xi\in\left(y_{j},y_{j+1}\right)$. Thus

[TABLE]

since $\left\|\psi^{h}\right\|_{L^{2}}=\left\|\psi_{\rm in}^{h}\right\|_{L^{2}}=1$, in view of mass conservation established in Theorem 2.4. In addition, we have that $F^{h}_{j}$ is uniformly bounded, i.e.

[TABLE]

Coming back to $($III$)$, we first note

[TABLE]

where the last equality follows from the fact that we have a telescoping series in $k$ with zero boundary conditions. Recalling (4.1), $($III$)$ can be written as

[TABLE]

where summation by parts is used in the last equality. In view of the Lipschitz property above, it is then straightforward to estimate $($III$)$ via

[TABLE]

Similarly, one proves that

[TABLE]

where summation by parts is used as before. Combining the coefficients of $F_{j}^{h}$ and $\left|F_{j}^{h}\right|$ respectively, this is equal to

[TABLE]

In summary, we thus find

[TABLE]

Using the fact that

[TABLE]

is a conserved quantity with respect to time, we consequently find the following estimate

[TABLE]

By Gronwall’s inequality, this yields

[TABLE]

which gives the desired bound independent of $h$. ∎

Remark 4.2**.**

It is easy to find a sharper bound of the energy by considering times $t\leq e$ and $t>e$, respectively, but the estimate above is sufficient for our purposes.

4.2 The classical limit of the s-SLE system

In this section, we shall perform the limit $h\rightarrow 0^{+}$ of the s-SLE system (4.1). By proving that it converges, as $h\to 0_{+}$, to the semi-discretized version of the coupled Liouville-system (3.2)–(3.3), we infer that it is possible to choose a spatial meshing strategy such that $\Delta y,\Delta\eta\sim\mathcal{O}(1)$.

To this end, we first note that the a-priori bounds on the mass and energy obtained in Theorems 2.4 and 4.1, together with our assumptions on $V\geq 0$ imply that the solution $\psi^{h}$ of (4.1) is uniformly bounded in $L^{2}\left(\mathbb{R}\right)$, and $h-$oscillatory. Thus, there exists an associated Wigner measure $\nu(,\cdot,\cdot,t)\in\mathcal{M}^{+}(\mathbb{R}_{x}\times\mathbb{R}_{\xi})$ and we directly infer that

[TABLE]

by the same arguments as in [15] (recall that $x,t$ are taken to be continuous in (4.1)). In the following we shall use the short-hand notation $F^{h}_{j}(t)$ and $F_{j}^{0}(t)$, respectively, as given in (4.1).

4.3 Convergence of $\mu^{h}$

Next, we turn to the solution $\mu^{h}$ within (4.1), which we recall to be discretized via an upwind scheme. We shall prove the following result about its limiting behavior as $h\to 0_{+}$. In the following, $C_{\rm b}(\mathbb{R})$ denotes the space of continuous and bounded functions on $\mathbb{R}$.

Proposition 4.3**.**

Let$\ \mu_{jk}^{h}\left(t\right)\in C_{\rm b}\left(\mathbb{R}_{t}\right)$ be a solution of

[TABLE]

and $\mu_{jk}\left(t\right)\in C_{\rm b}\left(\mathbb{R}_{t}\right)$ be a solution of

[TABLE]

where $j=0,\cdots,J-1$ and $k=0,\cdots,K-1$, such that initially $\mu_{jk}^{h}\left(0\right)=\mu_{jk}\left(0\right)$. Then for any given $T>0$,

[TABLE]

up to the extraction of subsequences.

Proof.

Denote the difference between $\mu_{jk}^{h}$ and its limit by

[TABLE]

which solves the following system of equations:

[TABLE]

subject to initial data $e_{jk}^{h}\left(0\right)=0$, since $\mu_{jk}^{h}(0)=\mu_{jk}(0)$. For simplicity we shall write the system above in vector form, i.e.,

[TABLE]

Here $E^{h}\left(t\right)$ and $b^{h}\left(t\right)$ are both $JK$-dimensional vectors, $A^{h}\left(t\right)$ is a continuous $JK\times JK$ matrix-valued function of $t\in\left(0,T\right]$, and

[TABLE]

for $j=1,\ldots,J$, and $k=1,\ldots,K$. Clearly, $\left\|A^{h}\left(t\right)\right\|_{\infty}\leq C$, where $C$ is a constant independent of $h$, due to the fact that $\left|F_{j}^{h}\left(t\right)\right|$ and $\left|F_{j}^{0}\left(t\right)\right|$ are both bounded by some constants independent of $h$, see the proof of Theorem 4.1. Classical ODE theory then implies (see, e.g., [27]), that there is a matrix-valued function $S^{h}\left(t,s\right)$ such that the solution of (4.3) is given by Duhamel’s principle:

[TABLE]

Moreover, there exists a bound on the propagator $S$ of the

[TABLE]

where $C_{3}$ is a constant independent of $h$.

Next, we recall that $F_{j}^{h}\left(t\right)=F^{h}\left(y_{j},t\right)$ is uniformly bounded, by equation (4.2), and

[TABLE]

pointwise (up to the extraction of subsequences). In addition, $F_{j}^{h}\left(t\right)$ is easily seen to be equi-continuous in time, by the same type of argument as in [15]. Namely, by using Schrödinger’s equation, one finds

[TABLE]

Integrating by parts, it reads

[TABLE]

where the last inequality follows from the $h-$oscillatory nature of $\psi^{h}$.

This consequently implies that

[TABLE]

locally uniformly in $t$, up to extraction of some subsequence, which in turn yields convergence of $E^{h}\left(t\right)$ itself, as can be seen by considering its $m$-th component, for $m=1,\dots,JK$:

[TABLE]

where we have used (4.4). ∎

4.4 Equation for $\nu$ and the main result

With the convergence theorem of $\mu^{h}$ in hand, we can now state the following result, which represents the final step in our analysis.

Proposition 4.4**.**

Assume (A1) and (A2) and let $\Theta[\Upsilon_{d}^{h}]$ be the pseudo-differential operator defined in (3.1) applied to the trapezoidal quadrature approximation of $\Upsilon^{h}\left(x,t\right)$. Then it holds

[TABLE]

where

[TABLE]

The proof of this proposition follows from the same arguments as given in the proof of Lemma 4.5 in [15], and we therefore omit it here.

We are now in the position to state the main result of this section.

Theorem 4.5**.**

Let Assumption (A1) and (A2) hold. Then, for any $T>0$, the solution of semi-discretized SLE system (4.1) satisfies, up to extraction of sub-sequences,

[TABLE]

where $j=0,\cdots,J-1$ and $k=0,\cdots,K-1$. In addition, $\nu$ and $\mu_{jk}$ solve the semi-discretized Liouville-system

[TABLE]

Remark 4.6**.**

Numerical experiments show that the same type of behavior is true not only for mixed spectral-finite difference schemes, but also purely spectral schemes, see [15]. Our proof, however, only works for the former case due to the required positivity of the energy.

5 Time-discretization

We finally turn to the time-discretization of our splitting scheme (in one dimension $d=n=1$) as given by (2.3), (2.5). In this section we want to show that it is asymptotic preserving in the sense that in the limit $h\to 0_{+}$, it yields the corresponding time-splitting scheme of (3.2)–(3.3), i.e.,

[TABLE]

and

[TABLE]

In turn, this shows that $\Delta t\sim\mathcal{O}(1)$ can be chosen independent of the small parameter $h$. To this end, it suffices to show that in our time-spitting method $\psi^{h}$ is $h-$oscillatory, i.e.

[TABLE]

where the constant $C$ depends only on the final time $T$. Then, following the arguments given in the previous section, one has convergence of the forcing term $F^{h}_{j}$ as $h\to 0$. In turn, this yields convergence of our numerical scheme towards the corresponding scheme of the limiting equation, as stated in (5.1) and (5.2).

We consequently consider the splitting scheme (2.3), (2.5) and recall that in both splitting steps, the first equation, i.e., the quantum part is solved exactly in time. A straightforward calculation then shows that $\left\|h\partial_{x}\psi^{h}\right\|^{2}_{L_{x}^{2}}$ is conserved in the first splitting step (2.3), i.e.

[TABLE]

Next, we shall show that $\left\|h\partial_{x}\psi^{h}\right\|^{2}_{L_{x}^{2}}$ remains bounded during the second splitting step (2.5): Recall that $\Upsilon_{d}^{h}\left(x,t\right)$ is in fact independent of $t$, due to the fact that $\frac{d}{dt}\mu_{jk}^{h}=0$ in this step. Since

[TABLE]

we find

[TABLE]

where

[TABLE]

Since $\left\|\psi^{h}\right\|_{L_{x}^{2}}=1$ is conserved by our scheme, we thus have

[TABLE]

where $C_{0}=\left\|\partial_{x}V\left(x,y\right)\right\|_{L^{\infty}}$ is some constant independent of $h$. Hence, in the second splitting step (2.5) one has

[TABLE]

where we used the fact that $\left\|h\partial_{x}\psi^{h}\right\|_{L_{x}^{2}}$ is conserved during (2.3). In summary, this yields

[TABLE]

where the right hand side is some constant independent of $h$ thanks to the assumption on initial data (A2). This shows that $\psi^{h,n}$ is $h-$oscillatory for any $n\in{\mathbb{N}}$, with $0\leq t_{n}\leq T$ and the result follows from the arguments in the previous section.

Remark 5.1**.**

Note that all estimates above remain valid in the context of a Strang splitting scheme.

6 Numerical examples

In this final section, we shall report on a few numerical examples, which illustrate the validity of our algorithm and meshing strategy. To this end, we choose an interaction potential of the form

[TABLE]

and solve the one-dimensional SLE system on the interval $x\in\left[-\pi,\pi\right]$ and $y,\eta\in\left[-2\pi,2\pi\right]$ with periodic boundary conditions.

Example 6.1** ($\Delta t$ independent of $h$).**

We choose initial conditions for the SLE system (1.2) as follows:

[TABLE]

and

[TABLE]

Here, $C_{\rm N}>0$ is the normalization factor such that $\iint_{\mathbb{R}^{2}}\mu_{in}dy\,d\eta$=1. Here we use the time-splitting method with spectral-upwind scheme (i.e., with an upwind scheme for the Liouville’s equation). For $h=$ $\frac{1}{256},\frac{1}{1024},\frac{1}{4096}$, we fix the stopping time $T=0.5$ and choose $\Delta x=\frac{2\pi h}{16}$, $\Delta y=\Delta\eta=\frac{4\pi}{128}$. For each choice of $h$, we shall solve the SLE system first with $\Delta t\text{ independent of }h$ and, second, with $\Delta t=o\left(h\right)$. To be more specific, we compare the two cases where $\Delta t=0.01$ and $\Delta t=\frac{h}{10}$. It can be observed from Figure 1f, that the macroscopic position and current densities associated to the solution of Schrödinger’s equation agree well with each other.

In addition, we compare the numerical values of $\mu$ computed by $\Delta t=0.01$ and $\Delta t=\frac{h}{10}$ (denoted as $\mu_{1}$ and $\mu_{2}$, respectively). As shown in Table 1, the error is insensitive in $h$, showing a uniform in $h$ convergence in $\Delta t$.

Example 6.2** (Numerical error as $h$ decreases).**

In this example, we choose the same initial data for $\mu_{\rm in}$ as before and

[TABLE]

Now, we fix $\Delta t=$ $0.01$, a stopping time $T=0.4$, and $\Delta y=\Delta\eta=\frac{4\pi}{128}$. We choose $\Delta x=\frac{2\pi h}{16}$, for $h=\frac{1}{64},\frac{1}{128},\frac{1}{256},\frac{1}{512},\frac{1}{1024},\frac{1}{2048}$, respectively. The reference solution is computed with $\Delta t=\frac{h}{10}$. From the $\ell^{2}$-error plotted in Figure 2, one can see that although the error in the wave function increases as $h$ decreases, the error for the position density $|\psi^{h}|^{2}$ as well as for the macroscopic quantity $\mu$ does not change noticeably. This shows that $h-$independent time steps can be taken to accurately obtain physical observables, but not the wave function $\psi^{h}$ itself.

Example 6.3** (Convergence in time).**

Finally, to examine the convergence in time of our scheme, let the initial data be as in the example before. Fix $h=$ $\frac{1}{8192}$, a stopping time $T=0.4$ and a spatial discretization with $\Delta x=\frac{2\pi h}{16}$, $\Delta y=\Delta\eta=\frac{4\pi}{128}$. Choose $\Delta t=\frac{0.4}{32},\frac{0.4}{64},\frac{0.4}{128},\frac{0.4}{256},\frac{0.4}{512},\frac{0.4}{1024}$. The reference solution is computed with $\Delta t=\frac{0.4}{81920}$. The $\ell^{2}$-error is plotted in Figure 3, which shows first order accuracy in time of our scheme. Again, we see that the wave function exhibits errors several orders of magnitude larger than the physical observable densities.