On the Convergence of Time Splitting Methods for Quantum Dynamics in the Semiclassical Regime

TL;DR

This paper proves that time splitting methods for quantum dynamics converge uniformly with respect to the Planck constant, providing explicit error estimates for first and second order algorithms.

Contribution

It introduces a pseudo-metric to analyze convergence and establishes uniform error bounds for splitting methods in the semiclassical regime.

Findings

Uniform convergence in Planck constant $$ for splitting algorithms

Explicit error estimates for Lie-Trotter and Strang methods

Applicable to quantum and classical density comparisons

Abstract

By using the pseudo-metric introduced in [F. Golse, T. Paul: Archive for Rational Mech. Anal. 223 (2017) 57-94], which is an analogue of the Wasserstein distance of exponent $2$ between a quantum density operator and a classical (phase-space) density, we prove that the convergence of time splitting algorithms for the von Neumann equation of quantum dynamics is uniform in the Planck constant $\hbar$. We obtain explicit uniform in $\hbar$ error estimates for the first order Lie-Trotter, and the second order Strang splitting methods.

Equations

$$R=R^{\ast }\ge 0,{\mathrm{trace}}_{\mathfrak{H}}(R)=1.$$

$${\mathrm{trace}}_{\mathfrak{H}}(R^{1/2}(-{\hslash }^2{\Delta }_y+|y{|}^2)R^{1/2})<\infty .$$

$${\mathrm{trace}}_{\mathfrak{H}}((-{\hslash }^2{\Delta }_y+|y{|}^2{)}^{1/2}R(-{\hslash }^2{\Delta }_y+|y{|}^2{)}^{1/2})={\mathrm{trace}}_{\mathfrak{H}}(R^{1/2}(-{\hslash }^2{\Delta }_y+|y{|}^2)R^{1/2})<\infty$$

$${\mathrm{trace}}_{\mathfrak{H}}(T^{1/2}A{T}^{1/2})={\mathrm{trace}}_{\mathfrak{H}}(A^{1/2}T{A}^{1/2})\in [0,+\infty ].$$

$${\iint }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}(|x{|}^2+|\xi {|}^2)f(x,\xi )dxd\xi <\infty .$$

$$Q(x,\xi )=Q(x,\xi {)}^{\ast }\ge 0,{\mathrm{trace}}_{\mathfrak{H}}(Q(x,\xi ))=f(x,\xi ),{\iint }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}Q(x,\xi )dxd\xi =R.$$

$$E_{\hslash }(f,R):=\underset{Q\in \mathcal{C}(f,R)}{\inf }{\left({\iint }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}{\mathrm{trace}}_{\mathfrak{H}}(Q(x,\xi {)}^{1/2}c(x,\xi ,y,\hslash D_y)Q(x,\xi {)}^{1/2})dxd\xi \right)}^{1/2},$$

$$c(x,\xi ,y,\hslash D_y):=|x-y{|}^2+|\xi -\hslash D_y{|}^2,D_y:=-i{\nabla }_y.$$

$$W_{\hslash }(R)(x,\xi ):=\frac{1}{(2\pi {)}^d}{\int }_{{\mathbf{R}}^d}r(x+\frac{1}{2}\hslash y,x-\frac{1}{2}\hslash y)e^{-i\xi \cdot y}dy,$$

$$W_{\hslash }(R)(0,0)=-\frac{1}{(2\pi {)}^d}{\int }_{{\mathbf{R}}^d}|\psi (\frac{1}{2}\hslash y){|}^{2}dy<0.$$

$${\tilde{W}}_{\hslash }(R):=e^{\hslash {\Delta }_{x,\xi }/4}{W}_{\hslash }(R).$$

$$|q,p⟩(x):=(\pi \hslash {)}^{-d/4}{e}^{-|x-q{|}^2/2\hslash }{e}^{ip\cdot (x-q/2)/\hslash },x,q,p\in {\mathbf{R}}^d.$$

$${\mathrm{OP}}_{\hslash }^T(\mu ):=\frac{1}{(2\pi \hslash {)}^d}{\int }_{{\mathbf{R}}^d\times {\mathbf{R}}^d}|q,p⟩⟨q,p|\mu (dqdp).$$

$${\mathrm{OP}}_{\hslash }^T(1)=I.$$

$$\mathrm{dis}{\mathrm{t}}_{\mathrm{MK},2}(\rho ,{\rho }^{\prime }):=\underset{\pi \in \Pi (\rho ,{\rho }^{\prime })}{\inf }{\left({\int }_{{\mathbf{R}}^{2d}}(|q-q^{\prime }{|}^2+|p-p^{\prime }{|}^2)\pi (dqdpd{q}^{\prime }d{p}^{\prime })\right)}^{1/2},$$

$${\int }_{{\mathbf{R}}^{2d}\times {\mathbf{R}}^{2d}}$$

$$={\int }_{{\mathbf{R}}^{2d}}\varphi (q,p)\rho (dqdp)+{\int }_{{\mathbf{R}}^{2d}}{\varphi }^{\prime }(q^{\prime },p^{\prime }){\rho }^{\prime }(d{q}^{\prime }d{p}^{\prime })$$

$$E_{\hslash }(f,R{)}^2\ge d\hslash .$$

$$E_{\hslash }(f,R{)}^2\ge \mathrm{dis}{\mathrm{t}}_{\mathrm{MK},2}(f,{\tilde{W}}_{\hslash }(R){)}^2-d\hslash$$

$$E_{\hslash }(f,{\mathrm{OP}}_{\hslash }^T((2\pi \hslash {)}^d\mu ){)}^2\le \mathrm{dis}{\mathrm{t}}_{\mathrm{MK},2}(f,\mu {)}^2+d\hslash .$$

$$V^{-}\in L^{d/2}({\mathbf{R}}^d),\text{ and }V\in C^{1,1}({\mathbf{R}}^d).$$

$$H_{\lambda }(x,\xi ):=\frac{1}{2}\lambda |\xi {|}^2+V(x).$$

$$(X(t;x,\xi ),\Xi (t;x,\xi ))$$

$$\dot{X}=\lambda \Xi ,\dot{\Xi }=-\nabla V(X),(X(0;x,\xi ),\Xi (0;x,\xi ))=(x,\xi ).$$

$${\partial}_{t}f+\{H_{\lambda},f\}=0\,,\qquad f{\big{|}}_{t=0}=f^{in}\,.$$

$$\{x_j,x_k\}=\{{\xi }_j,{\xi }_k\}=0,\{{\xi }_j,x_k\}={\delta }_{jk}.$$

$${\mathcal{H}}_{\lambda }:=-\frac{1}{2}{\hslash }^2\lambda {\Delta }_y+V(y).$$

$$i\hbar{\partial}_{t}R=[\mathcal{H}_{\lambda},R]\,,\qquad R{\big{|}}_{t=0}=R^{in}\,.$$

$$E_{\hslash }(f^{in}\circ {\Phi }_{-t},U(t{)}^{\ast }{R}^{in}U(t))\le E_{\hslash }(f^{in},R^{in})\exp \left(\frac{1}{2}t(\lambda +\max (1,\mathrm{Lip}(\nabla V{)}^2))\right).$$

$$\left\{\begin{aligned} {}&R^{0}_{\hbar}=R^{in}_{\hbar}\,,\\ &i\hbar{\partial}_{t}A_{\hbar}=[-\tfrac{1}{2}{\hbar}^{2}{\Delta}_{x},A_{\hbar}]\,,&&\qquad A_{\hbar}{\big{|}}_{t=0}=R^{n}_{\hbar}\,,\qquad n\in\mathbf{N}\,,\\ &i\hbar{\partial}_{t}B_{\hbar}=[V(x),B_{\hbar}]\,,&&\qquad B_{\hbar}{\big{|}}_{t=0}=A_{\hbar}({\Delta}t)\,,\\ &R^{n+1}_{\hbar}=B_{\hbar}({\Delta}t)\,.\end{aligned}\right.$$

$$\mathrm{dis}{\mathrm{t}}_{\mathrm{MK},2}({\tilde{W}}_{\hslash }(R^n),{\tilde{W}}_{\hslash }(R(n\Delta t)))$$

Full text

1. Introduction

One of the main challenges in quantum dynamics and high frequency waves is that one needs to numerically resolve the small wave length which is computationally prohibitive [1, 6, 13, 11]. When a numerical method is developed one would like to know its mesh strategy, namely, the dependence of the time step and mesh size on the wave length $\hbar$ (for a misuse of notation in this article we will not distinguish the difference between the reduced Planck constant $\hbar$ and the wave length).

Finite difference schemes for the Schrödinger equation typically require both time step and mesh size in the semiclassical regime (i.e. for $\hbar\ll 1$) to be of order $O(\hbar)$ (see [16]), or even $o(\hbar)$. On the other hand, the time splitting spectral method can improve the time step to be of order $O(1)$ if only the physical observables are of interest [2]. An important mathematical object to understand these mesh strategies is the Wigner transform [18], which is a convenient tool to study the semiclassical limit of the Schrödinger equation [7, 15]. In fact, the mesh strategy of $\Delta t=O(1)$, for the time step $\Delta t$, of the time-splitting spectral method can only be understood in the Wigner framework, and not in terms of the wave function [2].

Since the solution to the Schrödinger equation is oscillatory with wave length of order $\hbar$, if one uses a standard metric, such as the $L^{2}$ or Sobolev norm, one would end up with an numerical error of order $O((\Delta t/\hbar)^{m})$ for some integer $m$ which depends on the order of the method. This will not allow one to see an $\hbar$ independent mesh strategy. The argument of an $\hbar$ independent time-step strategy in [2] for the time splitting discretization to the linear Schrödinger equation, which was also useful in establishing a similar mesh strategy for the nonlinear Erhenfest dynamics [5], was made at a formal level without quantifying the numerical error.One would be interested in finding a suitable metric which allows one to establish such a mesh strategy at the rigorous level. In the present paper, we use the pseudo-metric introduced in [9] to establish a uniform (in $\hbar$) error estimate of the time splitting methods for the von Neumann equation (which describes the evolution of mixed quantum states, and reduces to the Schrödinger equation in the case of pure quantum states [3]) in the semiclassical regime.

2. A Pseudo-Metric for the Classical Limit

Definition 2.1**.**

A density operator on $\mathfrak{H}:=L^{2}(\mathbf{R}^{d})$ is an operator $R\in\mathcal{L}(\mathfrak{H})$ such that

[TABLE]

The set of all density operators on $\mathfrak{H}$ will be denoted by $\mathcal{D}(\mathfrak{H})$.

In the definition above, the notation $\mathcal{L}(\mathfrak{H})$ designates the algebra of bounded linear operators defined on $\mathfrak{H}$. Henceforth, we also denote by $\mathcal{L}^{p}(\mathfrak{H})$ for all $p\geq 1$ the two-sided ideal of $\mathcal{L}(\mathfrak{H})$ whose elements are the operators $T\in\mathcal{L}(\mathfrak{H})$ such that $|T|^{p}=(T^{*}T)^{p/2}$ is a trace-class operator on $\mathfrak{H}$. For instance $\mathcal{L}^{1}(\mathfrak{H})$ and $\mathcal{L}^{2}(\mathfrak{H})$ are respectively the sets of trace-class and of Hilbert-Schmidt operators on $\mathfrak{H}$. The notation $\operatorname{trace}_{\mathfrak{H}}(T)$ designates the trace of $T\in\mathcal{L}^{1}(\mathfrak{H})$.

We denote by $\mathcal{D}^{2}(\mathfrak{H})$ the set of density operators on $\mathfrak{H}$ such that

[TABLE]

If $R\in\mathcal{D}^{2}(\mathfrak{H})$, one has

[TABLE]

as can be seen from the lemma below (applied to $A={\lambda}^{2}|y|^{2}-\hbar^{2}{\Delta}_{y}$ and $T=R$).

Lemma 2.2**.**

Let $T\in\mathcal{L}(\mathfrak{H})$ satisfy $T=T^{*}\geq 0$, and let $A$ be an unbounded operator on $\mathfrak{H}$ such that $A=A^{*}\geq 0$. Then

[TABLE]

Proof.

The definition of $T^{1/2}$ and $A^{1/2}$ can be found in Theorem 3.35 in chapter V, §3 of [14], together with the fact that $A^{1/2}$ and $T^{1/2}$ are self-adjoint.

If $\operatorname{trace}_{\mathfrak{H}}(T^{1/2}AT^{1/2})<\infty$, then $A^{1/2}T^{1/2}\in\mathcal{L}^{2}(\mathfrak{H})$ and the equality holds by formula (1.26) in chapter X, §1 of [14].

If $\operatorname{trace}_{\mathfrak{H}}(T^{1/2}AT^{1/2})=\infty$, then $\operatorname{trace}_{\mathfrak{H}}(A^{1/2}TA^{1/2})=+\infty$, for otherwise $T^{1/2}A^{1/2}$ and its adjoint $A^{1/2}T^{1/2}$ would belong to $\mathcal{L}^{2}(\mathfrak{H})$, so that $T^{1/2}AT^{1/2}\in\mathcal{L}^{1}(\mathfrak{H})$, which would be in contradiction with the assumption that $\operatorname{trace}_{\mathfrak{H}}(T^{1/2}AT^{1/2})=\infty$. ∎

Let $f\equiv f(x,\xi)$ be a probability density on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ such that

[TABLE]

Definition 2.3**.**

Let $f\equiv f(x,\xi)$ be a probability density on $\mathbf{R}^{2d}$ and let $R\in\mathcal{D}(\mathfrak{H})$. A coupling of $f$ and $R$ is a measurable operator-valued function $(x,\xi)\mapsto Q(x,\xi)$ such that, for a.e. $(x,\xi)\in\mathbf{R}^{d}\times\mathbf{R}^{d}$,

[TABLE]

The second condition above implies that $Q(x,\xi)\in\mathcal{L}^{1}(\mathfrak{H})$ for a.e. $(x,\xi)\in\mathbf{R}^{d}\times\mathbf{R}^{d}$. Since $\mathcal{L}^{1}(\mathfrak{H})$ is separable, the notion of strong and weak measurability are equivalent for $Q$. The set of couplings of $f$ and $R$ is denoted by $\mathcal{C}(f,R)$. Notice that the operator-valued function $(x,\xi)\mapsto f(x,\xi)R$ belongs to $\mathcal{C}(f,R)$.

In [9], one considers the following “pseudometric”: for each probability density $f$ on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ and each $R\in\mathcal{D}^{2}(\mathfrak{H})$,

[TABLE]

where the quantum transportation cost is the quadratic differential operator in $y$, parametrized by $(x,\xi)\in\mathbf{R}^{d}\times\mathbf{R}^{d}$:

[TABLE]

Let $R\in\mathcal{D}(\mathfrak{H})$. The Wigner transform of $R$ is

[TABLE]

where $r\equiv r(x,y)$ is the integral kernel of $R$. It is a well known fact that $W_{\hbar}(R)$ is real-valued (since $R=R^{*}$). It is also well known that $W_{\hbar}(R)$ is not necessarily nonnegative a.e. on $\mathbf{R}^{d}\times\mathbf{R}^{d}$. For instance, if $r(X,Y)=\psi(X)\overline{\psi(Y)}$ with $\psi$ odd, one has

[TABLE]

The Husimi transform of $R$ henceforth denoted $\tilde{W}_{\hbar}(R)$ is defined in terms of its Wigner transform by the formula

[TABLE]

Finally, we recall the definition of a Töplitz operator. The family of Schrödinger coherent states is

[TABLE]

Let $\mu$ be a positive Borel measure on $\mathbf{R}^{d}\times\mathbf{R}^{d}$; the Töplitz operator with symbol $\mu$ is

[TABLE]

One easily checks that, if $\mu$ is the Lebesgue measure (denoted by $1$), then

[TABLE]

Moreover, one easily checks that, if $\mu$ is a Borel probability measure on $\mathbf{R}^{d}\times\mathbf{R}^{d}$, then $\operatorname{OP}^{T}_{\hbar}((2\pi\hbar)^{d}\mu)\in\mathcal{D}(\mathfrak{H})$. In addition, if $\mu$ has finite second order moment as in (3), then one has $\operatorname{OP}^{T}_{\hbar}((2\pi\hbar)^{d}\mu)\in\mathcal{D}^{2}(\mathfrak{H})$.

The pseudo-metric $E_{\hbar}$ satisfies the following fundamental properties. Henceforth, we denote by $\operatorname{dist_{MK,2}}$ the Monge-Kantorovich or Wasserstein distance with exponent $2$ defined on the set of Borel probability measures satisfying the finite second order moment condition (3) (see chapter 7 in [17]), whose definition is recalled below.

Definition 2.4**.**

For all $\rho$ and $\rho^{\prime}$, Borel probability measures on $\mathbf{R}^{2d}$, we set

[TABLE]

where $\Pi(\rho,\rho^{\prime})$ designates the set of couplings of $\rho$ and $\rho^{\prime}$. More precisely, $\Pi(\rho,\rho^{\prime})$ is the set of Borel probability measures on $\mathbf{R}^{2d}\times\mathbf{R}^{2d}$ with first and second marginals $\rho$ and $\rho^{\prime}$ resp., i.e. such that

[TABLE]

for all $\phi,\phi^{\prime}\in C_{b}(\mathbf{R}^{2d})$.

Specifically, the proposition below explains how the pseudo-metric $E_{\hbar}$ compares to the Wasserstein distance $\operatorname{dist_{MK,2}}$.

Proposition 2.5**.**

Let $R\in\mathcal{D}^{2}(\mathfrak{H})$ and let $f$ be a probability distribution on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ with finite second order moment (3).

(a) One has

[TABLE]

(b) One has

[TABLE]

(c) For each Borel probability measure $\mu$ on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ with finite second order moment as in (3), one has

[TABLE]

The pseudo-metric $E_{\hbar}$ can be used to obtain a quantitative formulation of the classical limit of quantum mechanics, as explained in [9]. Henceforth, we denote by $V$ a real-valued function satisfying

[TABLE]

(Here, the notation $V^{-}$ designates the function $x\mapsto V^{-}(x):=\max(-V(x),0)$.)

Let ${\lambda}\geq 0$, and set

[TABLE]

(From the physical point of view, the parameter ${\lambda}\geq 0$ which appears in the definition of the Hamiltonian $H_{\lambda}$ can be thought of as the reciprocal mass of the particle whose dynamics is defined in terms of the Hamiltonian flow associated to $H_{\lambda}$, whose definition is recalled below. In the present paper, the parameter ${\lambda}$ is used only as a convenient notation for defining the various time-splitting algorithms considered.)

Since $V$ satisfies the second condition in (4), we deduce from the Cauchy-Lipschitz theorem that the Hamiltonian $H_{\lambda}$ generates a globally defined flow denoted

[TABLE]

on $\mathbf{R}^{d}\times\mathbf{R}^{d}$. In other words, $t\mapsto(X(t;x,\xi),\Xi(t;x,\xi))$ is the solution to the Cauchy problem

[TABLE]

Equivalently, for each probability distribution $f^{in}$ on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ satisfying (3), the function $f^{in}\circ\Phi_{-t}$, where $\Phi_{t}$ is the map $(x,\xi)\mapsto\Phi_{t}(x,\xi):=(X(t;x,\xi),\Xi(t;x,\xi))$, is the solution to the Cauchy problem for the Liouville equation

[TABLE]

Here, the notation $\{\cdot,\cdot\}$ designates the Poisson bracket defined on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ by the relations

[TABLE]

Likewise consider the quantum Hamiltonian

[TABLE]

The parameter ${\lambda}$ that appears in the definition of the operator $\mathcal{H}_{\lambda}$ has the same meaning, and is used similarly as in the classical setting.

The first condition in (4) implies that $\mathcal{H}_{\lambda}$ has a self-adjoint extension (still denoted by $\mathcal{H}_{\lambda}$) on $\mathfrak{H}$ (see Lemma 4.8b in chapter VI, §4 of [14]). By the Stone theorem, $U(t):=e^{it\mathcal{H}_{\lambda}/\hbar}$ is a unitary group on $\mathfrak{H}$, and, for each $R^{in}\in\mathcal{D}(\mathfrak{H})$, the density operator $R(t):=U^{*}(t)R^{in}U(t)$ is the generalized solution to the Cauchy problem for the von Neumann equation

[TABLE]

Theorem 2.6**.**

Let $R^{in}\in\mathcal{D}^{2}(\mathfrak{H})$ and let $f^{in}$ be a probability density on $\mathbf{R}^{d}\times\mathbf{R}^{d}$ satisfying (3). Then

[TABLE]

This is a straightforward variant of Theorem 2.7 in [9] with an external potential and without interaction potential (i.e. in the special case $N=n=1$). The parameter ${\lambda}\geq 0$ appearing in the statement above is the other (unessential) difference with the situation discussed in [9].

3. Main Result

The simple time-splitting method for the von Neumann equation is

[TABLE]

Theorem 3.1**.**

Let $V$ satisfy (4), and assume that $R^{in}\in\mathcal{D}^{2}(\mathfrak{H})$ is a Töplitz operator on $\mathfrak{H}$. Let $t\mapsto R_{\hbar}(t)$ be the solution of the Cauchy problem (6), and let $R_{\hbar}^{n}$ be the sequence of density operators constructed by the simple splitting method (7). Let $T>0$, and pick a time step ${\Delta}t\in(0,\tfrac{1}{2})$. Then, for each $n=0,\ldots,[T/{\Delta}t]$, the simple splitting method satisfies the following error estimate, stated in terms of the quadratic Monge-Kantorovich or Wasserstein distance between the Husimi functions of the approximate and the exact quantum density operators:

[TABLE]

where the constant $C_{T}$ depends only on $T$, ${\nabla}V(0)$ and $\operatorname{Lip}({\nabla}V)$, and is defined in formula (12) below.

Instead of the simple splitting method, one can instead consider the Strang splitting method

[TABLE]

Strang splitting is a second order (in ${\Delta}t$) method, so that the convergence rate obtained in the previous theorem can be improved as indicated below.

Theorem 3.2**.**

Let $V$ satisfy (4) and

[TABLE]

Let $R^{in}\in\mathcal{D}^{2}(\mathfrak{H})$ be a Töplitz operator on $\mathfrak{H}$, and let $t\mapsto R_{\hbar}(t)$ be the solution of the Cauchy problem (6). On the other hand, let $R_{\hbar}^{n}$ be the sequence of density operators constructed by the Strang splitting method (8). Let $T>0$, and pick a time step ${\Delta}t\in(0,\tfrac{1}{2})$. Then, for each $n=0,\ldots,[T/{\Delta}t]$, the Strang splitting method satisfies the following error estimate, stated in terms of the quadratic Monge-Kantorovich or Wasserstein distance between the Husimi functions of the approximate and the exact quantum density operators:

[TABLE]

where the constant $D_{T}$ depends only on $T$ and $\|{\nabla}^{m}V\|_{L^{\infty}}$ for $m=1,2,3$, and is defined in formula (16) below.

Our strategy is the following. First, Theorem 2.6 gives the error between the solution of the von Neumann solution (6) and that of the classical Liouville equation (5). Then we obtain an analogous error between the time split von Neumann and the time split Liouville. Finally we estimate the time splitting error of the classical Liouville equation, measured in distance $\text{dist}_{MK2}$. Then a triangle type inequality leads to the results in Theorem 3.1 and 3.2. This strategy is best illustrated by Figure 1.

The error estimates Theorems 3.1 and 3.2 do not provide a uniform in $\hbar$ error estimate, since they contain an $O(\hbar^{1/2})$ term in their right hand side. In particular, these error estimates are useful only in the vanishing $\hbar$ limit. Yet these two theorems contain all the new information on the time splitting methods for quantum dynamics in the semiclassical regime that can be obtained with our approach. Besides, these two theorems are of independent interest, and lead to better convergence rates that the uniform error estimate given below in the vanishing $\hbar$ regime. In contrast, a classical $L^{2}$ norm estimate gives an error of order $O(\Delta t/\hbar)^{m}$ [1] (for some positive integer $m$ that depends on the order of the splitting), which blows up as $\hbar\to 0$.

In order to obtain uniform in $\hbar$ error estimates for the simple and the Strang splitting methods, we need to optimize these estimates with the error estimates for the time splitting method in the case of the Schrödinger equation with fixed $\hbar$ (or equivalently for $\hbar=1$). Such error estimates have been studied in detail and can be found for instance in [4]. The idea of combining and optimizing the error estimates in the asymptotic (macroscopic) regime and the microscopic regime is often used in numerical methods for kinetic and hyperbolic equations involving multiple scales, a computational methodology known as Asymptotic-Preserving Schemes [8, 12].

Our final uniform error estimate will be formulated in terms of an optimal transport distance denoted $\operatorname{dist}_{1}$, already used in [10] (see formula (13) in [10]). All the convergence statements in Theorems 3.1 and 3.2 are ultimately formulated in terms of the Monge-Kantorovich-Wasserstein distance $\operatorname{dist_{MK,2}}$, to which the “pseudo-metric” $E_{\hbar}$ can be conveniently compared (see Proposition 2.5), the uniform in $\hbar$ error estimates stated below as Corollaries 3.4 and 3.5 are all based on some optimization procedure comparing the $L^{1}$ and the $\operatorname{dist_{MK,2}}$ distances between the Husimi transforms of the exact and of the approximate solutions of the quantum dynamical problem. This optimization procedure is precisely the reason for using the distance $\operatorname{dist}_{1}$, a weaker variant of the Monge-Kantorovich-Wasserstein distance of exponent $1$, with a transportation cost that is truncated at infinity. The definition of $\operatorname{dist}_{1}$ is recalled below for the reader’s convenience.

Definition 3.3**.**

For all $\rho$ and $\rho^{\prime}$, Borel probability measures on $\mathbf{R}^{2d}$, we set

[TABLE]

Here, the notation $\Pi(\rho,\rho^{\prime})$ designates the set of couplings of $\rho,\rho^{\prime}$ already used to define the Monge-Kantorovich-Wasserstein distance $\operatorname{dist_{MK,2}}$ (Definition 2.4).

Our uniform estimates for the simple splitting method is given in the following statement, which is a consequence of Theorem 3.1 and of the error estimate in Theorem 2 of [4].

Corollary 3.4**.**

Let $V\in C^{2}(\mathbf{R}^{d})$ satisfy (4), and let $R^{in}=\operatorname{OP}^{T}_{\hbar}((2\pi\hbar)^{d}\mu^{in})$, where $\mu^{in}$ is a Borel probability measure on $\mathbf{R}^{2d}$ such that

[TABLE]

Let $t\mapsto R_{\hbar}(t)$ be the solution of the Cauchy problem (6), and let $R_{\hbar}^{n}$ be the sequence of density operators constructed by the simple splitting method (8). Let $T>0$, and pick a time step ${\Delta}t\in(0,\tfrac{1}{2})$. Then, for each $n=0,\ldots,[T/{\Delta}t]$, the simple splitting method satisfies the following uniform in $\hbar$ error estimate:

[TABLE]

where $C[T,V,\mu^{in}]$ is defined in (20). In particular, the constant $C[T,V,\mu^{in}]$ is independent of $\hbar$.

Likewise, Theorem 3.2 and the error estimate in Theorem 3 of [4] lead to the following statement.

Corollary 3.5**.**

Let $V\in W^{4,\infty}(\mathbf{R}^{d})$ satisfy (4), and let $R^{in}=\operatorname{OP}^{T}_{\hbar}((2\pi\hbar)^{d}\mu^{in})$, where $\mu^{in}$ is a Borel probability measure on $\mathbf{R}^{2d}$ such that

[TABLE]

Let $t\mapsto R_{\hbar}(t)$ be the solution of the Cauchy problem (6), and let $R_{\hbar}^{n}$ be the sequence of density operators constructed by the Strang splitting method (7). Let $T>0$, and pick a time step ${\Delta}t\in(0,\tfrac{1}{2})$. Then, for each $n=0,\ldots,[T/{\Delta}t]$, the simple splitting method satisfies the following uniform in $\hbar$ error estimate:

[TABLE]

where $D[T,V,\mu^{in}]$ is defined in (21). In particular, the constant $D[T,V,\mu^{in}]$ is independent of $\hbar$.

As will be clear from the proofs, the uniform in $\hbar$ estimates obtained in these two corollaries involve the $O(\hbar^{1/2})$ term in the convergence rates in Theorems 3.1 and 3.2, and the nonuniform bounds in Theorems 2 and 3 resp. of [4].

Specifically, the uniform $O({\Delta}t^{1/3})$ error bound in Corollary 3.4 is obtained as the minimum of the $O(\hbar^{1/2})$ term in Theorem 3.1 and of the (nonuniform) $O({\Delta}t/\hbar)$ error bound in Theorem 2 of [4]. This $O({\Delta}t^{1/3})$ uniform error estimate corresponds to the “worst” possible distinguished asymptotic regime $\hbar\sim{\Delta}t^{2/3}$. Although the $O({\Delta}t)$ term in Theorem 3.1 is smaller than the $O({\Delta}t^{1/3})$ uniform error estimate in Corollary 3.4, the $O({\Delta}t+\hbar^{1/2})$ error estimate in Theorem 3.1 is still of independent interest in all cases where $\hbar$ is small and satisfies $\hbar=O({\Delta}t^{2})$.

Likewise, the uniform $O({\Delta}t^{2/3})$ error bound in Corollary 3.5 comes as the minimum of the $O(\hbar^{1/2})$ term in the convergence rates in Theorem 3.2 and of the (nonuniform) $O({\Delta}t^{2}/\hbar^{2})$ error bound in Theorem 3 of [4]. In this case, the “worst” possible distinguished asymptotic regime is $\hbar\sim{\Delta}t^{4/5}$. Here again, the $O({\Delta}t^{2})$ term in Theorem 3.2 is smaller than the $O({\Delta}t^{2/3})$ uniform error estimate in Corollary 3.5. Nevertheless, the $O({\Delta}t^{2}+\hbar^{1/2})$ error estimate in Theorem 3.2 is of interest independently of the uniform $O({\Delta}t^{2/3})$ bound in Corollary 3.5 whenever $\hbar$ is small and satisfies $\hbar=O({\Delta}t^{4})$. Observe that the Strang splitting method is of second order in time in that regime, for the quantum dynamics as well as for the classical dynamics.

4. The Simple Splitting Algorithm

4.1. The Simple Splitting Algorithm for the von Neumann Equation in the Semiclassical Regime

In this subsection we estimate the error between the time split von Neumann and the time split Liouville equations. By analogy with the simple splitting method (7) for the von Neumann equation, consider the simple time-splitting method for the Liouville equation:

[TABLE]

Applying Theorem 2.6 to one time step of the free dynamics, i.e. with $V\equiv 0$ and ${\lambda}=1$ shows that

[TABLE]

Next we apply the same Theorem 2.6 to the Hamiltonian dynamics defined by the potential $V$, with ${\lambda}=0$: thus

[TABLE]

Putting both estimates together shows that

[TABLE]

Let $T>0$; then for each $n=0,\ldots,[T/{\Delta}t]+1$, one has

[TABLE]

Observe that the amplification rate $\exp(\tfrac{1}{2}T(1+\max(1,\operatorname{Lip}({\nabla}V)^{2})))$ in this estimate is uniform in (i.e. independent of) $\hbar$. This is the key point in our analysis.

4.2. The Simple Splitting Algorithm for the Liouville Equation

In this subsection, we estimate the distance between the classical Liouville equation and its time split approximation. The error analysis for the simple splitting method is well known in general. However, for our purpose in this paper, we shall formulate this error analysis for this splitting method applied to the Liouville equation in terms of the quadratic Monge-Kantorovich or Wasserstein distance.

One expresses the solutions $a$ and $b$ of the kinetic and the potential part of the Liouville evolution as follows, by using the method of characteristics. For the kinetic phase

[TABLE]

As for the potential phase

[TABLE]

Hence, one step of simple splitting corresponds to setting

[TABLE]

with

[TABLE]

Since the transformation $P_{{\Delta}t}\circ K_{{\Delta}t}$ has Jacobian one, the formula (10) means that $f^{n+1}(y,\eta)dyd\eta$ is the image of the measure $f^{n}(y,\eta)dyd\eta$ by $P_{{\Delta}t}\circ K_{{\Delta}t}$.

Next we seek to estimate the splitting error

[TABLE]

where $q^{n}$ is any coupling of $f(n{\Delta}t,\cdot,\cdot)$ and $f^{n}$.

For this, we seek to bound

[TABLE]

where

[TABLE]

is the numerical particle trajectory. First we derive the dynamic equations for $(Y,H)$. Inverting these relations, and denoting $Y_{t}:=Y(t;y,\eta)$ and $H_{t}:=H(t;y,\eta)$ for simplicity, we see that

[TABLE]

Equivalently

[TABLE]

Differentiating in time, we find that

[TABLE]

or equivalently

[TABLE]

Thus, we seek to compare the trajectories of the two following differential systems:

[TABLE]

Therefore, we set

[TABLE]

with

[TABLE]

Set $E:=|{\nabla}V(0)|$. Then, by the mean value inequality,

[TABLE]

On the other hand

[TABLE]

so that

[TABLE]

By Gronwall’s inequality, setting ${\Lambda}:=\max(1,E,\|{\nabla}^{2}V\|_{L^{\infty}}))$, one has

[TABLE]

Assume that $0\leq t\leq\tfrac{1}{2}$ for simplicity; then

[TABLE]

Choosing at this point an optimal coupling $q^{n}$ of $f(n{\Delta}t,\cdot,\cdot)$ and $f^{n}$ (see Theorem 1.3 in [17] for the existence of an optimal coupling), one has

[TABLE]

and it remains to control the last term in the right hand side.

Since $f^{n}(y,\eta)dyd\eta$ is the image of the measure $f^{n-1}(y,\eta)dyd\eta$ by the transformation $P_{{\Delta}t}\circ K_{{\Delta}t}$, one has

[TABLE]

(by substitution in the left hand side), so that

[TABLE]

Denoting

[TABLE]

we easily check that

[TABLE]

Therefore

[TABLE]

Thus, we arrive at the inequality

[TABLE]

Iterating in $n$, we conclude that, for $n=0,1,\ldots,[T/{\Delta}t]$

[TABLE]

where

[TABLE]

4.3. Error Estimate for the Simple Splitting Method

According to Theorem 2.6, for each $n=0,1,\ldots$, one has

[TABLE]

and in particular

[TABLE]

for $n=0,\ldots,[T/{\Delta}t]$. Putting together (9), (11) and (13) shows that

[TABLE]

According to Proposition 2.5 (b) and using the triangle inequality for $\operatorname{dist_{MK,2}}$, we conclude that

[TABLE]

In particular, if $R^{in}$ is the Töplitz operator with symbol $(2\pi\hbar)^{d}f^{in}$, we conclude from Proposition 2.5 (c) that

[TABLE]

5. The Strang Splitting Algorithm

In this subsection we estimate the error between the time split von Neumann and the time split Liouville equations. The Strang time-splitting method for the Liouville equation is

[TABLE]

Applying Theorem 2.6 to one time step of the free dynamics, i.e. with $V\equiv 0$ and ${\lambda}=1$ shows that

[TABLE]

Next we apply the same Theorem 2.6 to the Hamiltonian dynamics defined by the potential $V$, with ${\lambda}=0$: thus

[TABLE]

Finally, we apply again Theorem 2.6 to the last time step of the free dynamics, so that

[TABLE]

Hence the uniform in ${\hbar}$ estimate (9) also holds for the Strang splitting method.

Next we analyze the Strang splitting method for the Liouville equation in terms of the Monge-Kantorovich or Wasserstein distance. With the same notation as in the previous section, we seek to bound

[TABLE]

In order to do so, we seek to bound

[TABLE]

where the numerical particle trajectory or bi-characteristic flow of the Liouville equation is

[TABLE]

Writing $Z_{t}:=Z(t,z,{\omega})$ and ${\Omega}_{t}:={\Omega}(t,z,{\omega})$ for simplicity, we first observe that

[TABLE]

Hence

[TABLE]

so that

[TABLE]

Likewise

[TABLE]

Summarizing, the numerical bi-characteristic field for the Strang splitting method is

[TABLE]

where

[TABLE]

Here again, we seek to compare the solution $(Z_{t},{\Omega}_{t})$ to the Strang splitting differential equation with the solution $(X_{t},\Xi_{t})$ of the Newton system of motion equations, i.e.

[TABLE]

Arguing as in the case of the simple splitting method, we observe that

[TABLE]

so that

[TABLE]

One easily checks that

[TABLE]

Setting

[TABLE]

we see that

[TABLE]

Choosing an optimal coupling $q^{n}$ of $f(n{\Delta}t,\cdot,\cdot)$ and $f^{n}$, one has

[TABLE]

Arguing as in the case of the simple splitting algorithm, one has

[TABLE]

by substitution in the integral on the left hand side, since $f^{n}(y,\eta)dyd\eta$ is the image of the measure $f^{n-1}(y,\eta)dyd\eta$ by the transformation $K_{\frac{1}{2}{\Delta}t}\circ P_{{\Delta}t}\circ K_{\frac{1}{2}{\Delta}t}$. Since

[TABLE]

and

[TABLE]

one has

[TABLE]

so that

[TABLE]

Hence

[TABLE]

so that, iterating in $n$,

[TABLE]

for $n=0,1,\ldots,[T/{\Delta}t]$ with $0<{\Delta}t\leq\tfrac{1}{2}$. In other words

[TABLE]

for $n=0,1,\ldots,[T/{\Delta}t]$ with $0<{\Delta}t\leq\tfrac{1}{2}$, with

[TABLE]

Putting together (9), (15) and (13) shows that

[TABLE]

By Proposition 2.5 (b) and the triangle inequality for $\operatorname{dist_{MK,2}}$,

[TABLE]

for $n=0,1,\ldots,[T/{\Delta}t]$ with $0<{\Delta}t\leq\tfrac{1}{2}$.

In particular, if $R^{in}$ is the Töplitz operator with symbol $(2\pi\hbar)^{d}f^{in}$, we conclude from Proposition 2.5 (c) and the inequality above that

[TABLE]

6. Uniform in $\hbar$ Error Estimates

Proof of Corollary 3.4.

Throughout this section, we denote

[TABLE]

and

[TABLE]

For the first order time splitting, one has

[TABLE]

Hence

[TABLE]

At this point, we apply Theorem 2 from [4] for the error of the simple splitting scheme:

[TABLE]

where

[TABLE]

One has

[TABLE]

so that

[TABLE]

Next we apply Lemmas 8.2 and 8.1 in [10]:

[TABLE]

Using (18), (14) to bound the right hand side of (19) shows that

[TABLE]

so that

[TABLE]

with

[TABLE]

∎

Proof of Corollary 3.5.

Arguing as in the proof of Corollary 3.4 for the Strang splitting, we write

[TABLE]

and

[TABLE]

so that

[TABLE]

By Theorem 3 from [4]:

[TABLE]

where the constant $M^{\prime}$ depends on the final time $T$, on $\|V\|_{W^{4,\infty}(\mathbf{R}^{d})}$, and on

[TABLE]

since

[TABLE]

Thus

[TABLE]

where

[TABLE]

Optimizing in $\hbar$ leads to

[TABLE]

corresponding to the choice $\hbar={\Delta}t^{4/3}$. ∎

Acknowledgements. The work of François Golse and Thierry Paul was partly supported by LIA LYSM (co-funded by AMU, CNRS, ECM and INdAM). The work of Shi Jin was supported by NSFC grants Nos. 31571071 and 11871297.