<figure><img src="image_1.png" /></figure>

HAL
open science

# On the approximation of the von Neumann equation in the semiclassical limit. Part II: numerical analysis

François Golse, Francis Filbet

► To cite this version:

François Golse, Francis Filbet. On the approximation of the von Neumann equation in the semiclassical limit. Part II: numerical analysis. 2025. hal-05045845v2

HAL Id: hal-05045845

https://hal.science/hal-05045845v2

Preprint submitted on 11 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# ON THE APPROXIMATION OF THE VON NEUMANN EQUATION IN THE SEMICLASSICAL LIMIT. PART II : NUMERICAL ANALYSIS

FRANCIS FILBET AND FRANÇOIS GOLSE

ABSTRACT. This paper is devoted to the numerical analysis of the Hermite spectral method proposed in [14], which provides, in the semiclassical limit, an asymptotic preserving approximation of the von Neumann equation. More precisely, it relies on the use of so-called Weyl's variables to effectively address the stiffness associated to the equation. Then by employing a truncated Hermite expansion of the density operator, we successfully manage this stiffness and provide error estimates by leveraging the propagation of regularity in the exact solution.

Keywords : Quantum mechanics, von Neumann equation, Semiclassical limit, Hermite spectral method.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Quantum dynamics in Weyl's variables</td><td>3</td></tr><tr><td>1.2. Semiclassical limit</td><td>5</td></tr><tr><td>2. Hermite Spectral method</td><td>6</td></tr><tr><td>2.1. Hermite approximation of the von Neumann equation</td><td>7</td></tr><tr><td>2.2. Hermite approximation of the semiclassical limit equation</td><td>9</td></tr><tr><td>2.3. Preliminary results</td><td>10</td></tr><tr><td>3. Analysis of the semiclassical model</td><td>11</td></tr><tr><td>3.1. Propagation of regularity for the semiclassical limit equation</td><td>11</td></tr><tr><td>3.2. Proof of Theorem 1.2</td><td>14</td></tr><tr><td>3.3. Proof of Theorem 2.2</td><td>14</td></tr><tr><td>4. Analysis of the von Neumann equation</td><td>17</td></tr><tr><td>4.1. Propagation of regularity on the von Neumann equation</td><td>18</td></tr><tr><td>4.2. Proof of Theorem 2.1</td><td>18</td></tr><tr><td>5. Numerical simulations</td><td>21</td></tr><tr><td>6. Conclusion</td><td>22</td></tr><tr><td>Appendix A. Proof of Proposition 4.1</td><td>23</td></tr><tr><td>References</td><td>29</td></tr></table>

## 1. INTRODUCTION

Quantum dynamics exhibits high-frequency waves and requires the numerical resolution of small wavelengths, which presents substantial computational difficulties [2, 13, 19, 23]. Hence, for such problems, it is essential to consider an appropriate mesh strategy, which includes the correlation between time steps, mesh size, and the physical wavelength $\lambda = 2\pi\hbar/p$, where $\hbar$ is the Planck constant and $p$ the particle momentum [24, section 17, Chapter III]. Similarly, numerical schemes

<footnote>Date: June 11, 2025.</footnote>

<footnote n="2020">2020 Mathematics Subject Classification. Primary 35Q40, Secondary 65N08, 65N35.</footnote>
used to solve the Schrödinger equation typically demand that both the time step $\Delta t$ and mesh size, in the semiclassical regime when $\hbar \ll 1$, be of order $O(\hbar)$, or may even $o(\hbar)$ [1, 30]. Conversely, time splitting methods have the potential to enhance large time steps when only the physical observables are concerned, as has been shown in [3]. A fundamental concept for understanding these mesh strategies is the Wigner transform, which serves as a useful tool for investigating the semiclassical limit of the Schrödinger equation.

Furthermore, in [16], uniform error estimates of time splitting methods have been established for the von Neumann equation in the semiclassical regime. This equation describes the evolution of mixed quantum states, and reduces to the Schrödinger equation in the case of pure quantum states. More precisely, an operator $\hat{\rho}^{\hbar}$, called the density-matrix operator, satisfies the Liouville-von Neumann or simply von Neumann equation, in the operator formulation

$$ \begin{cases} i \hbar \partial_t \hat{\rho}^{\hbar} = [H, \hat{\rho}^{\hbar}], \\ \hat{\rho}^{\hbar}(0) = \hat{\rho}_{\text{in}}^{\hbar}, \end{cases} $$

where $\hbar$ is the reduced Planck constant. The operator $\hat{\rho}^{\hbar}$ is acting on the Hilbert space $L^2(\mathbb{R}^d) := L^2(\mathbb{R}^d, \mathbb{C})$ equipped with the inner product

$$ (1.1) \qquad \langle f, g \rangle := \int_{\mathbb{R}^d} f(x) \overline{g(x)} dx, $$

with associated norm given by $\| \cdot \| := \sqrt{\langle \cdot, \cdot \rangle}$. Moreover, the quantum mechanical Hamiltonian $H$ is given, throughout this paper by

$$ H := -\frac{\hbar^2}{2m} \Delta + V(X), $$

where $V$ a real-valued function such that $H$ is a self-adjoint operator on $L^2(\mathbb{R}^d)$. On the other hand, for all $\phi \in L^2(\mathbb{R}^d)$, we set $V\phi(X) := V(X)\phi(X)$ while

$$ [H, \hat{\rho}^{\hbar}] := H \hat{\rho}^{\hbar} - \hat{\rho}^{\hbar} H $$

is the commutator. Here, we suppose that $\hat{\rho}^{\hbar}(t)$ is a density matrix or operator on $L^2(\mathbb{R}^d)$, assuming that

$$ \hat{\rho}^{\hbar}(t)^* = \hat{\rho}^{\hbar}(t) \ge 0, \quad \text{Tr}(\hat{\rho}^{\hbar}(t)) = 1, \quad t \ge 0. $$

In particular, $\hat{\rho}^{\hbar}(t)$ is trace-class and therefore Hilbert-Schmidt, hence by [12, Theorem 6.12], there exists a unique function $\rho^{\hbar}(t|X, Y) \in L^2_{X,Y} := L^2(\mathbb{R}^d \times \mathbb{R}^d, \mathbb{C})$, such that for each $\phi \in L^2(\mathbb{R}^d)$, we define $\hat{\rho}^{\hbar}(t)\phi$ as the function

$$ (1.2) \qquad X \mapsto \hat{\rho}^{\hbar}(t) \phi(X) := \int_{\mathbb{R}^d} \rho^{\hbar}(t|X, Y) \phi(Y) dY. $$

Thus, using this latter representation, the von Neumann equation can be written as follows for the function $(X, Y) \mapsto \rho^{\hbar}(t|X, Y)$:

$$ (1.3) \qquad \begin{cases} i \hbar \partial_t \rho^{\hbar}(t|X, Y) = -\frac{\hbar^2}{2m} (\Delta_X - \Delta_Y) \rho^{\hbar}(t|X, Y) + (V(X) - V(Y)) \rho^{\hbar}(t|X, Y), \\ \rho^{\hbar}(t=0|\cdot, \cdot) = \rho_{\text{in}}^{\hbar}. \end{cases} $$

First, it is important to emphasize that when $\hbar$ is significantly smaller than 1, the dynamics governed by the von Neumann equation exhibits stiffness, resulting in a multiscale problem. Consequently, one of the primary challenges in quantum dynamics is the development of efficient numerical techniques capable of describing a wide range of wavelengths. For example, we refer the
reader to [22, 33, 34, 35] and the associated literature for discussions concerning the numerical ap-
proximation of the Schrödinger equation within the semiclassical limit. Alternative methods based
on stochastic algorithms [41] may be also applied in this context.

Furthermore, the solution of the von Neumann equation is a density operator, i.e. a Hermit-
ian, positive semi-definite operator with trace equal to 1. Preserving these intrinsic properties at
the discrete level is obviously crucial to ensure the generation of realistic physical outcomes. In
this regard, an initial approach was introduced in [18] which involves the successive application
of short-time propagators evaluated using fast Fourier transforms. This method exhibits consider-
able efficiency for short time intervals, particularly when the density matrix is spatially localized.
However, it is worth noting that it requires a small time step of the order of ħ or potentially even
smaller. Later in [5], the authors proposed a numerical scheme for a related model, again based on
the Fourier pseudo-spectral method. It allows a description both in configuration as well as in mo-
mentum space. More recently, various structure preserving schemes have been studied to guarantee
trace conservation and positivity of the discrete density matrix allowing long time simulations.

Recently, we proposed a new approach based on the so-called Weyl's variables introduced in part I [14] of the present paper to address the stiffness associated with the equation, and applied a truncated Hermite expansion of the density operator : see [14] for more details. Actually, this approach is inspired by the link between the von Neumann equation and the Wigner equation, which share the same kind of numerical difficulties arising from the nonlocal and highly-oscillating pseudo-differential operator associated to the potential energy. In [39, 42], it has been shown that the Wigner function can be discretized in a robust way using adaptive pseudo-spectral methods, where the oscillatory components introduced by the Wigner kernel are solved exactly. Finally, the oscillatory quantum effects can also be mitigated by decomposing the potential into classical and quantum parts [37, 38] or by reformulating the Wigner equation using spectral components of the force field [40]. Furthermore, as we will see, the use of Hermite expansion allows for the partial preservation of the operator structure $\hat{\rho}^\hbar$, which will be helpful in ensuring the stability of this approach. It is important to emphasize that in the Wigner equation, the term involving the potential is nonlocal. However, when the equation is formulated in Weyl's variables, this term remains local since it still appears as a multiplication. Additionally, the dependence on $\hbar$ is significantly reduced. After Hermite discretization, our approach can be viewed as a form of discretization of the Wigner equation, but this is specific to the Hermite-Galerkin method since Hermite functions are eigenvalues of the Fourier transform. Using another discretization method would lead to a different analysis of the Wigner equation and of the von Neumann equation written in Weyl's variable (see (1.6) below).

Our aim here is to analyze the Hermite spectral method proposed in [14], for equation (1.3). We begin by rewriting the von Neumann equation in terms of Weyl's variables, which eliminates the stiffness associated with $\hbar \ll 1$. The resulting equation is well-suited for the semiclassical limit as $\hbar$ tends to zero (Theorem 1.2). Next, we expand the density matrix using Hermite polynomials, leading to a spectrally accurate approximation. We provide new estimates on the solution of (1.3), demonstrating the propagation of regularity in a weighted Sobolev space. Taking advantage of the regularity of the solution, we first derive error estimates in terms of the number of Hermite modes, which are uniform with respect to the parameter $\hbar$ (Theorem 2.1). Additionally, we establish error estimates for the Hermite spectral method applied to the semiclassical limit model (Theorem 2.2). Sections 3 and 4 are dedicated to the proofs of Theorems 2.2 and 2.1. Finally, we perform numerical simulations to illustrate the spectral accuracy of the Hermite method in the last section and refer the reader to [14] for more extensive numerical simulations.

1.1. **Quantum dynamics in Weyl's variables.** To remove the stiffness due to the parameter $\hbar$
in the von Neumann equation, the key idea consists in rewriting (1.3) in the new variables defined
below, called “Weyl’s variables”:

$$ (1.4) \qquad \left\{ \begin{array}{l} x := \frac{X + Y}{2} \\ y := \frac{X - Y}{\hbar}, \end{array} \right. $$

which can be explicitly inverted as

$$ X = x + \frac{\hbar}{2} y, \quad Y = x - \frac{\hbar}{2} y. $$

Although this terminology, introduced in [14], is not common, it is obviously suggested by Weyl’s quantization (see formula (1.1.9) of [27]). Hence, we set

$$ (1.5) \qquad R^h(t|x, y) := \rho^h(t|X, Y), $$

and find that $\rho^h$ is solution to (1.3) if and only if $R^h$ satisfies

$$ (1.6) \qquad \left\{ \begin{array}{l} \partial_t R^h(t|x, y) = i \sum_{j=1}^{d} \partial_{x_j} \partial_{y_j} R^h(t|x, y) - i \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} R^h(t|x, y), \\ R^h(t = 0|..) = R_{\text{in}}^h. \end{array} \right. $$

Let us first review some elementary properties of the operator $\hat{\rho}^h$ and their interpretation in terms of the function $R^h$. We first recall that $\hat{\rho}^h$ is a density operator on $L^2(\mathbb{R}^{2d})$, that is $\hat{\rho}^h(t)^* = \hat{\rho}^h(t)$, for all $t \in \mathbb{R}$, which means that $R^h$ satisfies

$$ (1.7) \qquad R^h(t|x, -y) = \overline{R^h(t|x, y)}, $$

whereas the trace property becomes

$$ (1.8) \qquad \operatorname{Tr} (\hat{\rho}^h(t)) = \int_{\mathbb{R}^d} \rho^h(t|X, X) \, dX = \int_{\mathbb{R}^d} R^h(t|x, 0) \, dx = 1. $$

Translating the positivity condition on $\hat{\rho}^h(t) \ge 0$ in terms of the function $R^h$ is more difficult : for all $\phi \in L^2(\mathbb{R}^d)$, we have

$$ \begin{align*} 0 \le \langle \phi | \hat{\rho}^h(t) | \phi \rangle &:= \iint_{\mathbb{R}^{2d}} \overline{\phi(X)} \rho^h(t|X, Y) \phi(Y) \, dX \, dY \\ &= \hbar^d \iint_{\mathbb{R}^{2d}} \overline{\phi(x + \frac{\hbar}{2} y)} R^h(t|x, y) \phi(x - \frac{\hbar}{2} y) \, dx \, dy. \end{align*} $$

We refer to [17] for a characterization of rank-1 density operators in Weyl’s variables. Finally, for all $t \ge 0$, we define $S_h(t)$ to be the one parameter group associated to the von Neumann equation (1.6), so that, for all $R_{\text{in}}^h \in L^2(\mathbb{R}^{2d})$,

$$ (1.9) \qquad R^h(t) = S_h(t) R_{\text{in}}^h. $$

A simple change of variables implies the conservation of the $L^2$ norm for equation (1.6), summarized in the following proposition.

**Proposition 1.1.** Consider $R^h$ the solution to the von Neumann equation (1.6). Then, we have for each $t \ge 0$,

$$ (1.10) \qquad \iint_{\mathbb{R}^{2d}} |R^h(t|x, y)|^2 dx dy = \|S_h(t) R_{\text{in}}^h\|_{L^2}^2 = \|R_{\text{in}}^h\|_{L^2}^2 < \infty. $$
Proof. On the one hand, since for all $t$,

$$ \hat{\rho}^{\hbar}(t) = \exp(-itH/\hbar) \hat{\rho}^{\hbar}(0) \exp(itH/\hbar), $$

one has

$$ \hat{\rho}^{\hbar}(t)^* = \hat{\rho}^{\hbar}(t) \quad \text{since } \hat{\rho}^{\hbar}(0)^* = \hat{\rho}^{\hbar}(0), \quad \text{and} \quad \text{Tr}(\hat{\rho}^{\hbar}(t)^2) = \text{Tr}(\hat{\rho}^{\hbar}(0)^2) < \infty. $$

Since the Hilbert-Schmidt norm of an operator is the $L^2$ norm of its integral kernel,

$$ \|\rho^{\hbar}(t|\cdot, \cdot)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)}^2 = \text{Tr}(\hat{\rho}^{\hbar}(t)^2) = \text{Tr}(\hat{\rho}^{\hbar}(0)^2) = \|\rho^{\hbar}(0|\cdot, \cdot)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)}^2. $$

The simple change of variables (1.4), leads to

$$ \|\rho^{\hbar}(t|\cdot, \cdot)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)}^2 = \hbar^{\hbar} \|R^{\hbar}(t|\cdot, \cdot)\|_{L^2(\mathbb{R}^d \times \mathbb{R}^d)}^2, $$

which implies in turn the conservation of the $L^2$ norm for equation (1.6) : for each $t \in [0, +\infty)$,

$$ \iint_{\mathbb{R}^{2d}} |R^{\hbar}(t|x, y)|^2 dx dy = \|R_{\text{in}}^{\hbar}\|_{L^2}^2 < \infty. $$

□

This latter property will be crucial in the stability analysis of our numerical method but, it is not enough in itself to prove convergence. We shall also need to establish some regularity estimates obtained by propagating moments and derivatives in $L^2$ for the solution of the von Neumann equation (1.6).

1.2. **Semiclassical limit.** The formulation of the von Neumann equation using Weyl's variables facilitates the semiclassical limit. Passing to the limit as $\hbar \to 0$ in the latter equation does not lead to any difficulty related to the stiffness in $\hbar$. More precisely, since $V \in \mathcal{C}^1(\mathbb{R}^d)$ we have :

$$ \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} \to \nabla V(x) \cdot y \quad \text{as } \hbar \to 0. $$

Therefore, setting

$$ (1.11) \qquad \mathcal{E}^{\hbar}(x, y) := \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} - \nabla V(x) \cdot y, $$

the function $R^{\hbar}$ now solves

$$ (1.12) \qquad \left\{ \begin{array}{l} \partial_t R^{\hbar}(t|x, y) = i \sum_{j=1}^{d} \partial_{x_j} \partial_{y_j} R^{\hbar}(t|x, y) - i \left( \nabla V(x) \cdot y + \mathcal{E}^{\hbar}(x, y) \right) R^{\hbar}(t|x, y), \\ R^{\hbar}(t=0|\cdot, \cdot) = R_{\text{in}}^{\hbar}, \end{array} \right. $$

where $\mathcal{E}^{\hbar}$ goes to zero as $\hbar \to 0$. The error made in approximating the solution $R^{\hbar}$ of (1.6), or equivalently (1.12), by the solution $R$ obtained after neglecting the source term $\mathcal{E}^{\hbar}$, corresponds to the semiclassical limit, given by

$$ (1.13) \qquad \left\{ \begin{array}{l} \partial_t R(t|x, y) = i \sum_{j=1}^{d} \partial_{x_j} \partial_{y_j} R(t|x, y) - i \nabla V(x) \cdot y R(t|x, y), \\ R(t=0|\cdot, \cdot) = R_{\text{in}}. \end{array} \right. $$

This limit was made rigorous by P.-L. Lions and T. Paul [28], and P. Markowich and C. Ringhofer [29, 31] for each smooth potential $V \in \mathcal{C}^{\infty}$. We also refer to [32, Section 1.4], where an asymptotic expansion of $\mathcal{E}^{\hbar}$ in powers of $\hbar$ is performed for smooth potentials $V \in \mathcal{C}^{\infty}$. In the sequel, we define $S_0(t)$ to be the one-parameter group corresponding to the semiclassical limit equation (1.13), so that

$$ (1.14) \qquad R(t) = S_0(t) R_{\text{in}}. $$
We now restrict ourselves to the one-dimensional case $d = 1$ and present an error estimate between the solutions $R^h(t|\cdot, \cdot)$ of the von Neumann equation in Weyl's variables (1.6) and $R(t|\cdot, \cdot)$ of the semiclassical equation (1.13) in terms of $\hbar$. The key point here is to establish some regularity estimates on the solution to (1.13). To this effect, we introduce the following quantity for any function $F$ and each integer $m \ge 0$:

$$ (1.15) \qquad \mathcal{N}_m(t) := \sum_{\substack{a+b+\alpha+\beta \le m \\ a,b,\alpha,\beta \in \mathbb{N}}} \|x^a \partial_x^b y^\alpha \partial_y^\beta F(t)\|_{L^2}. $$

Let us fix an integer $m \ge 0$ and suppose that the potential $V$ is such that

$$ (1.16) \qquad V \in \mathcal{C}^{m+1}(\mathbb{R}) \text{ with } m \ge 1, \text{ and } V(x) \to +\infty \text{ as } |x| \to \infty, $$

and for all $n = 2, \dots, m+1$, there exists a constant $C > 0$, independent of $\hbar$, such that

$$ (1.17) \qquad \|V^{(n)}\|_{L^\infty} \le C, \quad \text{with } V^{(n)} := \frac{d^n V}{dx^n}. $$

We prove the following result.

**Theorem 1.2.** Assume that $V$ satisfies (1.16)-(1.17) with $m = 3$ while for all integers $a, b, \alpha$ and $\beta$ such that $a+b+\alpha+\beta \le 3$, we have

$$ (1.18) \qquad \left\{ \begin{array}{l} \|x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}\|_{L^2} < \infty, \\ \|R_{\text{in}}^\hbar\|_{L^2} < \infty. \end{array} \right. $$

Then the solution $R^\hbar$ to the von Neumann equation (1.6) and the solution $R$ to (1.13) with the initial data $R_{\text{in}}^\hbar$ and $R_{\text{in}}$ satisfy

$$ \|R^\hbar(t) - R(t)\|_{L^2} \le \|R_{\text{in}}^\hbar - R_{\text{in}}\|_{L^2} + C(t) \hbar^2, $$

where $C(t)$ is given by the following expression

$$ C(t) := \frac{1}{24} \|V^{(3)}\|_{L^\infty} \frac{e^{C_3[V]^t} - 1}{C_3[V]} \mathcal{N}_3(0), $$

with $C_3[V]$ defined later in (3.1) and $\mathcal{N}_3(0)$ given in (1.15) applied to $R_{\text{in}}$.

The proof of this theorem is presented in Section 3. Actually, it is relatively elementary since it is a consequence of a simple Taylor expansion of the potential $V$, but the propagation of regularity of the solution to the semiclassical limit problem (1.13) is required. This will be achieved by studying the functional $\mathcal{N}_m(t)$, defined in (1.15), applied to the solution $R(t)$ to (1.13). This result can be interpreted as a modulated energy estimate where only regularity on the solution of the limit equation is required (see for instance [4] or [11]).

## 2. HERMITE SPECTRAL METHOD

The purpose of this section is to present a formulation of the von Neumann equation (1.6) written in Weyl's variables based on Hermite polynomials. We first use Hermite polynomials in the $y$ variable and write the von Neumann equation (1.6) as an infinite hyperbolic system for the Hermite coefficients depending only on time and space. The idea is to apply a Galerkin method keeping only a small finite set of orthogonal polynomials rather than discretizing the operator $R^\hbar$ in $y$ on a grid. The merit of using an orthogonal basis like the so-called scaled Hermite basis has been shown in [21, 36] and more recently in [7, 6, 10, 9, 8, 15] for the Vlasov-Poisson system. We define the weight function $\omega$ as

$$ \omega(y) := \pi^{-1/4} e^{-y^2/2}, \quad y \in \mathbb{R}, $$
and the sequence of Hermite polynomials (called the “physicist’s Hermite polynomials”) by

$$H_k(y) := (-1)^k \omega^{-2}(y) \partial_y^k \omega^2(y), \quad k \ge 0.$$

In this context the family of Hermite functions $(\Phi_k)_{k \in \mathbb{N}}$ defined as

$$\Phi_k(y) := \frac{1}{\sqrt{2^k k!}} \omega(y) H_k(y), \quad k \ge 0$$

is a complete, orthonormal system for the classical $L^2$ inner product, that is,

$$\int_{\mathbb{R}} \Phi_k(y) \Phi_l(y) \, dy = \delta_{k,l}, \quad \text{and} \quad \overline{\text{Span}\{\Phi_k, \ k \ge 0\}} = L^2(\mathbb{R}).$$

The sequence $(H_k)_{k \in \mathbb{N}}$ defined in (2.1) satisfies the following recursion relation :

$$(2.1) \qquad \begin{cases} H_{-1} = 0, & H_0 = 1, \\ H_{k+1}(y) = 2y H_k(y) - 2k H_{k-1}(y), & k \ge 0. \end{cases}$$

We now expand the solutions to the von Neumann and semiclassical limit equations in terms of Hermite functions.

**2.1. Hermite approximation of the von Neumann equation.** We consider the decomposition of $R^h$, solution to (1.6), into its components $R^h = (R_k^h)_{k \in \mathbb{N}}$ in the Hermite basis :

$$R^h(t|x, y) = \sum_{k \in \mathbb{N}} R_k^h(t|x) \Phi_k(y),$$

where

$$R_k^h(t|x) := \int_{\mathbb{R}} R^h(t|x, y) \Phi_k(y) \, dy, \quad k \in \mathbb{N}.$$

Since the Hermite functions verify the following relations [20] :

$$(2.2) \qquad \Phi'_k(y) = -\sqrt{\frac{k+1}{2}} \Phi_{k+1}(y) + \sqrt{\frac{k}{2}} \Phi_{k-1}(y), \quad \forall k \ge 0,$$

and

$$(2.3) \qquad y \Phi_k(y) = \sqrt{\frac{k+1}{2}} \Phi_{k+1}(y) + \sqrt{\frac{k}{2}} \Phi_{k-1}(y), \quad \forall k \ge 0,$$

we have

$$\int_{\mathbb{R}} \partial_y R^h(t|x, y) \Phi_k(y) \, dy = -\sqrt{\frac{k}{2}} R_{k-1}^h(t|x) + \sqrt{\frac{k+1}{2}} R_{k+1}^h(t|x),$$

whereas

$$\int_{\mathbb{R}} R^h(t|x, y) y \Phi_k(y) \, dy = \sqrt{\frac{k}{2}} R_{k-1}^h(t|x) + \sqrt{\frac{k+1}{2}} R_{k+1}^h(t|x),$$

with $R_{-1}^h = 0$. Furthermore, the contribution of the potential $V$ in (1.11), now becomes

$$\int_{\mathbb{R}} \left\{ \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} - V'(x)y \right\} R^h(t|x, y) \Phi_k(y) \, dy = \sum_{l \in \mathbb{N}} \mathcal{E}_{k,l}^h(x) R_l^h(t|x),$$

with

$$(2.4) \qquad \mathcal{E}_{k,l}^h(x) := \int_{\mathbb{R}} \mathcal{E}^h(x, y) \Phi_k(y) \Phi_l(y) \, dy, \quad k, l \in \mathbb{N},$$
where the function $\mathcal{E}^h(\cdot, \cdot)$ is defined in (1.11). Before going further, let us first review some properties of $\mathcal{E}_{k,l}^h$. Obviously $\mathcal{E}_{k,l}^h(x) \in \mathbb{R}$ and we also have

$$ (2.5) \qquad \mathcal{E}_{k,l}^h(x) = \mathcal{E}_{l,k}^h(x), \quad x \in \mathbb{R}, k, l \ge 0. $$

Furthermore, using that $H_k$ is proportional to the even function $\omega^{-2}$ times the $k$-th derivative of the even function $\omega^2$, it shows that

$$ (2.6) \qquad H_k(y) = (-1)^k H_k(-y) $$

and from (1.11), we also get that

$$ \mathcal{E}^h(x, y) = \mathcal{E}^h(x, -y). $$

Therefore, applying the change of variable $y \mapsto -y$ in (2.4), we have

$$ (2.7) \qquad \mathcal{E}_{k,l}^h(x) = -(-1)^{k+l} \mathcal{E}_{k,l}^h(x), \quad x \in \mathbb{R}, \quad k, l \ge 0. $$

This latter property implies that $\mathcal{E}_{k,l}^h$ is identically 0 when ever $k+l$ is even. We thus end up with the formulation of the von Neumann equation written in the Hermite basis for all integers $k \ge 0$,

$$ \left\{ \begin{array}{l} \partial_t R_k^h + i \left( \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^h + \sqrt{\frac{k+1}{2}} \mathcal{D}^* R_{k+1}^h \right) = -i \sum_{l \ge 0} \mathcal{E}_{k,l}^h R_l^h, \\ R_k^h(t=0|\cdot) = R_{\text{in},k}^h, \end{array} \right. $$

where $\mathcal{D}$ and $\mathcal{D}^*$ are given by

$$ \left\{ \begin{array}{l} \mathcal{D} u := +\partial_x u + V'(x) u, \\ \mathcal{D}^* u := -\partial_x u + V'(x) u. \end{array} \right. $$

Let us emphasize an important property satisfied by $\mathcal{D}$, which we will need to recover later on, in the discrete setting. The operator $\mathcal{D}^*$ is its adjoint operator in $L^2(\mathbb{R})$, meaning that for all $u$, $v \in \text{Dom}(\mathcal{D})$, it holds

$$ (2.8) \qquad \langle \mathcal{D}^* u, v \rangle = \langle u, \mathcal{D} v \rangle, $$

where $\langle \cdot, \cdot \rangle$ denotes the classical inner product (1.1) in $L^2(\mathbb{R})$. Then we define the (closed) linear space $\mathcal{P}_N(\mathbb{R}) \subset L^2(\mathbb{R})$ of dimension $N+1$:

$$ (2.9) \qquad \mathcal{P}_N(\mathbb{R}) := \text{Span}\{\Phi_k, \quad k = 0, \dots, N\}, $$

and construct an approximation $R^{h,N} = (R_k^{h,N})_{0 \le k \le N} \in \mathcal{P}_N(\mathbb{R})$, solution to the following system obtained after neglecting Hermite modes of order larger than $N$, that is, for $0 \le k \le N$

$$ (2.10) \qquad \left\{ \begin{array}{l} \partial_t R_k^{h,N} + i \left( \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^{h,N} + \sqrt{\frac{k+1}{2}} \mathcal{D}^* R_{k+1}^{h,N} \right) = -i \sum_{l=0}^{N} \mathcal{E}_{k,l}^h R_l^{h,N}, \\ R_k^{h,N}(t=0|\cdot) = R_{\text{in},k}^h, \end{array} \right. $$

where $R_{-1}^{h,N} = R_{N+1}^{h,N} = 0$. The right-hand side of (2.10) will later be designated as "the nonlocal term in (2.10)". (this may be slightly improper since the left-hand side also involves $R_{k \pm 1}^{h,N}$ which can be viewed as nonlocal in terms of $(R_k^{h,N})_{0 \le k \le N}$)

We are now ready to present our second main result on error estimates of between the Hermite-Galerkin approximation and the solution to the von Neumann equation (1.6). The use of Weyl variables allows us to obtain uniform accuracy with respect to the parameter $\hbar$, ensuring naturally the asymptotic preserving property of the Hermite approximation for the von Neumann equation. Again the key point here is to establish some regularity estimates on the solution to the von Neumann equation written in Weyl's variables.
**Theorem 2.1.** Let $p \ge 1$ and assume that $V$ satisfies (1.16)-(1.17) with $m = 2(p+1)$, while for all integers $a, b, \alpha$ and $\beta$ such that $a+b+\alpha+\beta \le m$,

$$ (2.11) \qquad \|x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}^\hbar\|_{L^2} < \infty. $$

Then the solutions $R^\hbar$ to (1.6) and $R^{\hbar,N}$ to (2.10) satisfy

$$ \begin{aligned} \|R^{\hbar,N}(t) - R^{\hbar}(t)\|_{L^2} & \le \frac{K_p}{(2N+3)^p} e^{C_{2p}[V]t} \mathcal{N}_{2p}(0) \\ & \quad + \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N+3)^{p-1/2} C_{2p+1}[V]} \mathcal{N}_{2p+1}(0) \\ & \quad + \frac{K_{p+1} \hbar^2}{4(2N+3)^{p-1/2}} \|V^{(3)}\|_{L^\infty} \frac{e^{C_{2p+2}[V]t} - 1}{C_{2p+2}[V]} \mathcal{N}_{2p+2}(0), \end{aligned} $$

where $\mathcal{N}_m(t)$ is defined in (1.15) and applied to $R^\hbar$, whereas $C_p, K_p$ and $K_{p+1}$ are constants depending only on $p$, while $C_l[V]$, for $l = 2p, \dots, 2p+2$, is the constant defined later in Proposition 4.1.

This result is particularly remarkable because, on the one hand, it provides spectral accuracy, that is, the convergence rate in $N$ depends only on the regularity of the solution to the continuous problem (1.6). On the other hand, these estimates are uniform with respect to the parameter $\hbar \ll 1$, which guarantees a uniform accuracy of the numerical approximation for different physical regimes. Thus we can guarantee the asymptotic preservation property in the limit where $\hbar$ tends to 0.

The proof is split in two steps. First, we establish stability estimates for the $L^2$ norm of $R^{\hbar,N}$. In other words we obtain a uniform in $\hbar$ error estimate for the approximation of the solution of the von Neumann equation (in Weyl's variables) by the Hermite Galerkin system (2.10). Then we prove the propagation of regularity of the solution to the von Neumann equation (1.6) written in Weyl's variables.

**2.2. Hermite approximation of the semiclassical limit equation.** As we have seen before, the formulation (2.10) is well adapted to the semiclassical limit as $\hbar \to 0$. Indeed, the Hermite formulation of the semiclassical limit equation (1.13) is simply obtained by neglecting the right-hand side term in (2.10). Thus, we construct an approximation $R^N = (R_k^N)_{0 \le k \le N}$, which satisfies, for all $k \in \{0, \dots, N\}$,

$$ (2.12) \qquad \left\{ \begin{array}{l} \partial_t R_k^N + i \left( \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^N + \sqrt{\frac{k+1}{2}} \mathcal{D}^* R_{k+1}^N \right) = 0, \\ R_k^N(t=0|.) = R_{\text{in},k}. \end{array} \right. $$

Following the same strategy as the one applied for the von Neumann equation, we provide error estimates of the Hermite-Galerkin method for smooth solutions of the semiclassical limit equation (1.13). Hence, using the propagation of regularity, we prove the following result.

**Theorem 2.2.** Let $p \ge 1$ and assume that $V$ satisfies (1.16)-(1.17) with $m = 2p+1$ while for all integers $a, b, \alpha$ and $\beta$ such that $a+b+\alpha+\beta \le m$, the initial data $R_{\text{in}}$ satisfies

$$ (2.13) \qquad \|x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}\|_{L^2} < \infty. $$
Then the solutions R to (1.13) and R^N to (2.12) satisfy

$$
\begin{aligned}
\|R(t) - R^N(t)\|_{L^2} &\le \frac{K_p e^{C_{2p}[V]t}}{(2N+3)^p} \mathcal{N}_{2p}(0) \\
&\quad + \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N+3)^{p-1/2}} \frac{C_{2p+1}[V]}{C_{2p+1}[V]} \mathcal{N}_{2p+1}(0),
\end{aligned}
$$

where $\mathcal{N}_m(t)$ is defined in (1.15) and applied to $R$, whereas $C_p$ and $K_p$ are constants depending only on $p$, while the constants $C_l[V]$, for $l = 2p$ and $2p+1$ are given later in (3.1).

First let us emphasize that this result requires slightly less regularity on the initial data than the one in Theorem 2.1. This is due to the absence of non-local in $k$ terms in the right-hand side of (2.12), at variance with (2.10). Furthermore, in order to offer a smooth presentation, we will first present the proof of this latter result in Section 3 and then focus on the estimate of the non-local term in the right-hand side of (2.10) to prove Theorem 2.1 in Section 4.

It is worth mentioning that the proof of these three results follows the lines of the famous Lax theorem in numerical analysis [26, 25], since we make use of the regularity of the solution to get an error estimate. A similar approach, either called the “relative entropy method” or the “modulated energy method”, has also been widely used in kinetic theory to compare the solutions of a kinetic and fluid equations [4] or [11]. In this work, we treat simultaneously the error estimates with respect to the physical parameter $\bar{h}$ and the discrete parameter $N$.

**2.3. Preliminary results.** Before starting our analysis on the discretization of the von Neumann equation using Hermite functions, let us briefly recall some basic properties on spectral methods.

Consider $G \in L^2(\mathbb{R})$ written in the Hermite basis

$$G(y) = \sum_{k \in \mathbb{N}} G_k \Phi_k(y),$$

where

$$G_k := \int_{\mathbb{R}} G(y) \Phi_k(y) \, dy, \quad k \ge 0.$$

The Hermite-Galerkin method consists in truncating the series, keeping only the $N$ first coefficients $G_k$ in the expansion of the function $G$ on the basis $(\Phi_k)_{k \in \mathbb{N}}$ of Hermite functions. Thus, let $\mathcal{P}_N$ be the orthogonal projection from $L^2(\mathbb{R})$ to $\mathcal{P}_N(\mathbb{R})$, where $\mathcal{P}_N(\mathbb{R})$ is defined in (2.9). Let us first evaluate the "error" of the projection operator $\mathcal{P}_N : L^2(\mathbb{R}) \to \mathcal{P}_N(\mathbb{R})$ when applied to smooth functions.

**Lemma 2.3.** Let $G \in L^2(\mathbb{R})$ satisfy, for some integer $p \ge 0$,

$$\|y^{2p} G\|_{L^2} + \|G\|_{H^{2p}} < \infty.$$

Then, we have

$$\|G - \mathcal{P}_N G\|_{L^2} \le \frac{1}{(2N+3)^p} \left( \int_{\mathbb{R}} \left| \left(y^2 - \frac{\mathrm{d}^2}{\mathrm{d}y^2}\right)^p G(y) \right|^2 \mathrm{d}y \right)^{1/2}.$$

*Proof.* On the one hand, since the Hermite functions is a complete family of orthonormal functions in $L^2(\mathbb{R})$, we deduce from the Parseval equality that

$$\|G - \mathcal{P}_N G\|_{L^2}^2 = \sum_{k>N} |G_k|^2.$$
On the other hand, for each integer $p \ge 0$, we apply again the Parseval equality to find that

$$ \int_{\mathbb{R}} \left| \left(y^2 - \frac{d^2}{dy^2}\right)^p G \right|^2 dy = \sum_{k \ge 0} \left| \left( \left(y^2 - \frac{d^2}{dy^2}\right)^p G \right)_k \right|^2 $$

and observe that

$$ \frac{1}{2} \left( y^2 - 1 - \frac{d^2}{dy^2} \right) \Phi_k(y) = k \Phi_k(y), \quad k \ge 0. $$

Hence,

$$ \left(y^2 - \frac{d^2}{dy^2}\right)^p \Phi_k(y) = (2k+1)^p \Phi_k(y), \quad k \ge 0, $$

and for all integer $p \ge 0$, we find that

$$ \left( \left(y^2 - \frac{d^2}{dy^2}\right)^p G \right)_k = \int_{\mathbb{R}} G(y) \left(y^2 - \frac{d^2}{dy^2}\right)^p \Phi_k(y) dy = (2k+1)^p G_k. $$

Therefore,

$$ \int_{\mathbb{R}} \left| \left(y^2 - \frac{d^2}{dy^2}\right)^p G(y) \right|^2 dy = \sum_{k \ge 0} (2k+1)^{2p} |G_k|^2. $$

This leads to

$$ \begin{aligned} \|G - \mathcal{P}_N G\|_{L^2}^2 &= \sum_{k>N} |G_k|^2 \\ &\le \frac{1}{(2N+3)^{2p}} \sum_{k>N} (2k+1)^{2p} |G_k|^2 \\ &\le \frac{1}{(2N+3)^{2p}} \int_{\mathbb{R}} \left| \left(y^2 - \frac{d^2}{dy^2}\right)^p G(y) \right|^2 dy, \end{aligned} $$

which is the sought inequality. $\square$

For sake of simplicity we first present our analysis on the semiclassical model, where the non-local term is neglected. Then in the second part, we treat the von Neumann equation and focus on the contribution of the non-local term $\mathcal{E}^\hbar$.

### 3. ANALYSIS OF THE SEMICLASSICAL MODEL

This section is mainly devoted to the propagation of regularity on the solution to the semiclassical limit equation (1.13). We first study the semiclassical limit equation in order to present the basic ideas. Then, we provide an error estimate with respect to $\hbar$ between $R^\hbar$ and $R$ when the solution to the semiclassical limit equation $R$ is regular enough. Finally, we give an error estimate between the solution $R$ to (1.13) and its approximation $R^N$ using the Hermite-Galerkin method (2.12).

**3.1. Propagation of regularity for the semiclassical limit equation.** In this section we consider $R$, the solution to the semiclassical limit equation (1.13). We first recall the conservation of $L^2$ norm, and then prove the propagation of regularity leading to the following result. Observe that our assumptions on the smooth potential $V$ include the harmonic oscillator for which $V(x) = x^2$.

**Proposition 3.1.** Let $m \ge 0$ and let $V$ be a potential satisfying (1.16), (1.17), and assume that, for each integer $a, b, \alpha$ and $\beta$ such that $a+b+\alpha+\beta \le m$, the condition (2.13) on the initial data $R_{\text{in}}$ is satisfied. Then the solution $R$ to (1.13) with $d=1$ satisfies, for all $t \ge 0$,

$$ \mathcal{N}_m(t) \le \mathcal{N}_m(0) e^{C_m[V]t}, $$
where $\mathcal{N}_m$ is defined in (1.15) and $C_m[V] > 0$ is given by

$$ (3.1) \quad C_m[V] := (m^2 + 2m + m \|V''\|_{L^\infty} + m |V'(0)|) + 2^m (1+m) m \max_{2 \le n \le m+1} \|V^{(n)}\|_{L^\infty}. $$

*Proof.* First, observing that for all $t \ge 0$

$$ (3.2) \quad \|R(t)\|_{L^2}^2 = \int_{\mathbb{R}^2} |R(t|x, y)|^2 dx dy = \|S_0(t) R_{\text{in}}\|_{L^2}^2 = \|R_{\text{in}}\|_{L^2}^2, $$

where $S_0(t)$ represents the one-parameter group (1.14), corresponding to the semiclassical limit equation (1.13), for which the $L^2$ norm is preserved by the evolution.

Now, for each integer $a, b, \alpha$ and $\beta \ge 0$, applying $x^a \partial_x^b y^\alpha \partial_y^\beta$ to both sides of (1.13), we have

$$ (3.3) \quad \begin{aligned} \partial_t (x^a \partial_x^b y^\alpha \partial_y^\beta R) = & (i \partial_x \partial_y - i V'(x)y) (x^a \partial_x^b y^\alpha \partial_\xi^\beta R) \\ & + i [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y] R. \end{aligned} $$

On the one hand, computing the last term of the right-hand side in the former equality and using that $[z^\gamma, \partial_z] = -\gamma z^{\gamma-1}$ for $z = x$ (resp. $y$) and $\gamma = \alpha$ (resp. $\beta$), we find

$$ \begin{aligned} [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y] &= [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y] - [x^a \partial_x^b y^\alpha \partial_y^\beta, V'(x)y] \\ &= [x^a, \partial_x] \partial_x^b y^\alpha \partial_y^{\beta+1} + \partial_x x^a \partial_x^b [y^\alpha, \partial_y] \partial_y^\beta \\ &\quad - x^a [\partial_x^b, V'(x)] y^\alpha \partial_y^\beta y - V'(x) x^a \partial_x^b y^\alpha [\partial_y^\beta, y]. \end{aligned} $$

Then, we have

$$ \begin{aligned} [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y] &= -a x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} - \alpha \partial_x x^a \partial_x^b y^{\alpha-1} \partial_y^\beta \\ &\quad - \beta V'(x) x^a \partial_x^b y^\alpha \partial_y^{\beta-1} - x^a [\partial_x^b, V'(x)] y^\alpha \partial_y^\beta y \\ &= -a x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} - \alpha x^a \partial_x^{b+1} y^{\alpha-1} \partial_y^\beta \\ &\quad - \alpha a x^{a-1} \partial_x^b y^{\alpha-1} \partial_y^\beta - \beta V'(x) x^a \partial_x^b y^\alpha \partial_y^{\beta-1} \\ &\quad - x^a [\partial_x^b, V'(x)] y^{\alpha+1} \partial_y^\beta - \beta x^a [\partial_x^b, V'(x)] y^\alpha \partial_y^{\beta-1}. \end{aligned} $$

Finally, applying the Leibniz formula

$$ [\partial_x^b, V'(x)] = \sum_{j=1}^{b} \binom{b}{j} V^{(j+1)}(x) \partial_x^{b-j}, $$

we get

$$ (3.4) \quad \begin{aligned} -[x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y] = & a x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} + \alpha x^a \partial_x^{b+1} y^{\alpha-1} \partial_y^\beta \\ & + \alpha a x^{a-1} \partial_x^b y^{\alpha-1} \partial_y^\beta + \beta V'(x) x^a \partial_x^b y^\alpha \partial_y^{\beta-1} \\ & + \sum_{j=1}^b \binom{b}{j} x^a V^{(j+1)}(x) \partial_x^{b-j} (y^{\alpha+1} \partial_y^\beta + \beta y^\alpha \partial_y^{\beta-1}). \end{aligned} $$

On the other hand, we deduce from (3.3) and the Duhamel formula that

$$ \begin{aligned} x^a \partial_x^b y^\alpha \partial_y^\beta R(t|x, y) = & S_0(t) (x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}) \\ & + i \int_0^t S_0(t-s) [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y] R(s|x, y) ds, \end{aligned} $$
where $S_0(t)$ again denotes the one parameter group (1.14). Hence, using the conservation of the $L^2$ norm (3.2) by $S_0(t)$, we obtain that

$$ \begin{aligned} \|x^a \partial_x^b y^\alpha \partial_y^\beta R(t)\|_{L^2} &\le \|x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}\|_{L^2} \\ &\quad + \int_0^t \left\| \left[x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x) y\right] R(s) \right\|_{L^2} \mathrm{d}s. \end{aligned} $$

Thus, the former inequality and the computation of $[x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y - V'(x)y]$ in (3.4) imply that

$$ \begin{aligned} \|x^a \partial_x^b y^\alpha \partial_y^\beta R(t)\|_{L^2} &\le \|x^a \partial_x^b y^\alpha \partial_y^\beta R_{\text{in}}\|_{L^2} \\ &\quad + \int_0^t \left\{ a \|x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} R\|_{L^2} + \alpha \|x^a \partial_x^{b+1} y^{\alpha-1} \partial_y^\beta R\|_{L^2} \right\} (s) \mathrm{d}s \\ &\quad + \int_0^t \left\{ \alpha a \|x^{a-1} \partial_x^b y^{\alpha-1} \partial_y^\beta R\|_{L^2} + \beta |V'(0)| \|x^a \partial_x^b y^\alpha \partial_y^{\beta-1} R\|_{L^2} \right\} (s) \mathrm{d}s \\ &\quad + \beta \|V''\|_{L^\infty} \int_0^t \left\{ \|x^{a+1} \partial_x^b y^\alpha \partial_y^{\beta-1} R\|_{L^2} \right\} (s) \mathrm{d}s \\ &\quad + \max_{2 \le n \le b+1} \|V^{(n)}\|_{L^\infty} \sum_{j=1}^b \binom{b}{j} \int_0^t \|x^a \partial_x^{b-j} y^{\alpha+1} \partial_y^\beta R\|_{L^2}(s) \mathrm{d}s \\ &\quad + \beta \max_{2 \le n \le b+1} \|V^{(n)}\|_{L^\infty} \sum_{j=1}^b \binom{b}{j} \int_0^t \|x^a \partial_x^{b-j} y^\alpha \partial_y^{\beta-1} R\|_{L^2}(s) \mathrm{d}s, \end{aligned} $$

where we used the mean-value inequality to prove that, for all $x \in \mathbb{R}$,

$$ (3.5) \qquad |V'(x)| \le |V(0)| + |x||V''||_{L^\infty}. $$

Then we set

$$ N_j(t) := \sum_{a,b \ge 0} \sum_{\alpha, \beta \ge 0} \|x^a \partial_x^b y^\alpha \partial_y^\beta R(t)\|_{L^2} \delta_{a+b+\alpha+\beta,j}, $$

where $\delta_{k,j}$ denotes the Kronecker symbol, which yields for $j \ge 1$,

$$ \begin{aligned} N_j(t) \le N_j(0) &+ (j^2 + 2j + j \|V''\|_{L^\infty}) \int_0^t N_j(s) \mathrm{d}s + j |V'(0)| \int_0^t N_{j-1}(s) \mathrm{d}s \\ &+ 2^j (1+j) \max_{2 \le n \le j+1} \|V^{(n)}\|_{L^\infty} \sum_{k=1}^j \int_0^t N_{j+1-k}(s) \mathrm{d}s. \end{aligned} $$

Thus, observing that

$$ N_m(t) = \sum_{j=0}^{m} N_j(t), $$

we deduce the following inequality on $N_m(t)$:

$$ \begin{aligned} N_m(t) \le N_m(0) &+ (m^2 + 2m + m \|V''\|_{L^\infty} + m |V'(0)|) \int_0^t N_m(s) \mathrm{d}s \\ &+ 2^m (1+m) m \max_{2 \le n \le m+1} \|V^{(n)}\|_{L^\infty} \int_0^t N_m(s) \mathrm{d}s. \end{aligned} $$

Then we obtain the expected upper bound by applying Gronwall's inequality. □
This latter result is particularly important since it allows us to establish several error estimates. On the one hand, we prove an error estimate with respect to $\hbar$ between the solution $R^\hbar$ of the von Neumann equation and its semiclassical limit $R$ (Theorem 1.2). On the other hand, we provide an estimate with respect to the number of Hermite modes $N$ for the semiclassical limit equation (Theorem 2.2).

### 3.2. Proof of Theorem 1.2.
Consider the solution $R^\hbar$ to the von Neumann equation (1.6) and the solution $R$ to (1.13). Then, the difference $R^\hbar - R$ satisfies

$$ \partial_t(R^\hbar - R) = i \left( \partial_x \partial_y - \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} \right) (R^\hbar - R) - i \mathcal{E}^\hbar R, $$

where $\mathcal{E}^\hbar$ is defined in (1.11).

By Taylor's formula at order 3,

$$ V(x+z) = V(x) + V'(x) z + \frac{1}{2} V''(x) z^2 + \frac{1}{2} \int_0^1 (1-\theta)^2 V^{(3)}(x+\theta z) z^3 d\theta, $$

so that

$$ \begin{aligned} \mathcal{E}^\hbar(x,y) &= \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} - V'(x) y \\ &= \frac{\hbar^2}{16} y^3 \int_0^1 (1 - \theta)^2 \left( V^{(3)}(x + \frac{\hbar}{2} \theta y) + V^{(3)}(x - \frac{\hbar}{2} \theta y) \right) d\theta. \end{aligned} $$

Hence,

$$ |\mathcal{E}^\hbar(x,y)| = \left| \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} - V'(x) y \right| \le \frac{\hbar^2}{24} |y|^3 \|V^{(3)}\|_{L^\infty}. $$

Using this upper bound and applying Proposition 3.1 on $R$, it yields that

$$ \begin{aligned} \|R^\hbar(t) - R(t)\|_{L^2} &\le \|R_{\text{in}}^\hbar - R_{\text{in}}\|_{L^2} + \frac{\hbar^2}{24} \|V^{(3)}\|_{L^\infty} \int_0^t \|y^3 R(s)\|_{L^2} ds \\ &\le \|R_{\text{in}}^\hbar - R_{\text{in}}\|_{L^2} + \frac{\hbar^2}{24} \|V^{(3)}\|_{L^\infty} \frac{e^{C_3[V]t} - 1}{C_3[V]} \mathcal{N}_3(0). \end{aligned} $$

### 3.3. Proof of Theorem 2.2.
We now consider $R$ (resp. $R^N$), the solution to the semiclassical equation (1.13) (resp. the Galerkin semiclassical equation (2.12)). We will first establish some properties on $R^N$ that are also satisfied at the continuous level and next show that the $L^2$ estimate given in (1.10) is preserved at the discrete level when the number of Hermite modes $N$ is finite.

**Proposition 3.2.** For a given $m \ge 0$, assume that $V$ satisfies (1.16)-(1.17) while $R_{\text{in}}$ satisfies (2.13), and consider the solution $R^N$ to (2.12). Then

(i) for all $k, l \in \{0, \dots, N\}$ and $t \ge 0$, we have

$$ (3.6) \qquad R_k^N(t|x) = (-1)^k \overline{R_k^N(t|x)}, $$

and, assuming that $N$ is even, the trace property in (1.8) becomes

$$ \int_{\mathbb{R}} R^N(t|x, 0) dx = \sum_{k=0}^{N} \frac{H_k(0)}{\sqrt{2^k k!}} \int_{\mathbb{R}} R_k^N(t|x) dx = \int_{\mathbb{R}} \mathcal{P}_N R_{\text{in}}(x, 0) dx; $$

(ii) for all $t \ge 0$, we have

$$ (3.7) \qquad \|R^N(t)\|_{L^2_{x,y}} = \|\mathcal{P}_N R_{\text{in}}\|_{L^2_{x,y}}^2. $$
Proof. On the one hand, using (2.6), we have

$$ \Phi_k(y) = (-1)^k \Phi_k(-y), \quad y \in \mathbb{R}, \quad k \ge 0, $$

hence,

$$ \begin{aligned} R_{\text{in},k}(x) &= \int_{\mathbb{R}} R_{\text{in}}(x,y) \Phi_k(y) \, dy = \int_{\mathbb{R}} \overline{R_{\text{in}}(x,-y)} \Phi_k(y) \, dy \\ &= \int_{\mathbb{R}} \overline{R_{\text{in}}(x,z)} \Phi_k(-z) \, dz = (-1)^k \overline{R_{\text{in},k}(x)}. \end{aligned} $$

On the other hand, we easily show that if for all $k \in \{0, \dots, N\}$, $R_k^N$ is solution to (2.12), then $(-1)^k \overline{R_k^N}$ also solves (2.12). By uniqueness of the solution of (2.12), we get the expected result

$$ R_k^N(t) = (-1)^k \overline{R_k^N}(t), \quad \forall t \ge 0. $$

Now let us turn to the second equality, for which we suppose that $N$ is even and compute the time derivative of the trace of $R^N$ to find that

$$ \frac{d}{dt} \int_{\mathbb{R}} R^N(t|x,0) \, dx = \sum_{k=0}^{N} \frac{H_k(0)}{\sqrt{2^k k!}} \int_{\mathbb{R}} \partial_t R_k^N(t|x) \, dx = \mathcal{J}_1(t) + \mathcal{J}_2(t), $$

with

$$ \left\{ \begin{aligned} \mathcal{J}_1(t) &= -i \sum_{k=0}^{N} \frac{H_k(0)}{\sqrt{2^k k!}} \int_{\mathbb{R}} \left( \sqrt{\frac{k}{2}} \partial_x R_{k-1}^N - \sqrt{\frac{k+1}{2}} \partial_x R_{k+1}^N \right) dx, \\ \mathcal{J}_2(t) &= -i \sum_{k=0}^{N} \frac{H_k(0)}{\sqrt{2^k k!}} \int_{\mathbb{R}} V'(x) \left( \sqrt{\frac{k}{2}} R_{k-1}^N + \sqrt{\frac{k+1}{2}} R_{k+1}^N \right) dx. \end{aligned} \right. $$

We first observe that the term $\mathcal{J}_1(t)$ vanishes since $R_k^N(t|\cdot)$ is sought in $L^2(\mathbb{R}_x)$, hence $R_k^N(t|x) \to 0$ along sequences $x_j \to \pm\infty$. Furthermore, using that $R_{-1}^N = R_{N+1}^N = 0$, $H_1(0) = 0$ while $H_{N-1}(0) = H_{N+1}(0) = 0$ since $N$ is even shows that $\mathcal{J}_2(t)$ may be written as,

$$ \begin{aligned} \mathcal{J}_2(t) &= -i \sum_{k=1}^{N-1} \int_{\mathbb{R}} V'(x) R_k^N \left( \sqrt{\frac{k+1}{2}} \frac{H_{k+1}(0)}{\sqrt{2^{k+1}(k+1)!}} + \sqrt{\frac{k}{2}} \frac{H_{k-1}(0)}{\sqrt{2^{k-1}(k-1)!}} \right) dx \\ &= -i \sum_{k=1}^{N-1} \int_{\mathbb{R}} V'(x) \frac{R_k^N}{\sqrt{2^{k+2}k!}} (H_{k+1}(0) + 2k H_{k-1}(0)) dx = 0, \end{aligned} $$

since the recursive relation on Hermite polynomials (2.1) gives $H_{k+1}(0) + 2kH_{k-1}(0) = 0$. Hence, the result follows when $N$ is even.

We next prove (ii) for any $N \in \mathbb{N}$ using that

$$ \|R^N(t)\|_{L_{x,y}^2}^2 = \sum_{k=0}^{N} \|R_k^N(t)\|_{L_x^2}^2, \quad t \ge 0. $$

Thus, we have

$$ \frac{1}{2} \frac{d}{dt} \sum_{k=0}^{N} \|R_k^N(t)\|_{L^2}^2 = \operatorname{Re} \left( i \sum_{k=0}^{N} \left\langle \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^N + \sqrt{\frac{k+1}{2}} \mathcal{D}^* R_{k+1}^N, R_k^N \right\rangle \right), $$

where $\langle \cdot, \cdot \rangle$ denotes the standard inner product (1.1) on $L^2(\mathbb{R}; \mathbb{C})$. Using the duality property (2.8) on the operators $\mathcal{D}$ and $\mathcal{D}^*$, we prove that the term on the right-hand side vanishes. Indeed, we
have

$$ \sum_{k=0}^{N} \left\langle \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^{N} + \sqrt{\frac{k+1}{2}} \mathcal{D}^{*} R_{k+1}^{N}, R_{k}^{N} \right\rangle = \sum_{k=0}^{N} \left( \sqrt{\frac{k}{2}} \langle R_{k-1}^{N}, \mathcal{D}^{*} R_{k}^{N} \rangle + \sqrt{\frac{k+1}{2}} \langle \mathcal{D}^{*} R_{k+1}^{N}, R_{k}^{N} \rangle \right). $$

Reordering and recalling that $R_{-1}^N = R_{N+1}^N = 0$, we find that

$$ \sum_{k=0}^{N} \left( \sqrt{\frac{k}{2}} \langle R_{k-1}^{N}, \mathcal{D}^{*} R_{k}^{N} \rangle + \sqrt{\frac{k+1}{2}} \langle \mathcal{D}^{*} R_{k+1}^{N}, R_{k}^{N} \rangle \right) = \Sigma + \bar{\Sigma} \in \mathbb{R}, $$

with

$$ \Sigma = \sum_{k=1}^{N} \sqrt{\frac{k}{2}} \langle R_{k-1}^{N}, \mathcal{D}^{*} R_{k}^{N} \rangle. $$

Hence

$$ \frac{d}{dt} \sum_{k=0}^{n} \| R_k^N(t) \|^2 = 2 \operatorname{Re} (i (\Sigma + \bar{\Sigma})) = 0, $$

so that, the $L^2$ norm for (2.12) is conserved. $\square$

<figure><img src="image_2.png" /></figure>

We are now ready to prove Theorem 2.2. We consider $R$ the solution to the semiclassical equation (1.13) with $d=1$ and $R^N$ the solution to the semiclassical Galerkin system (2.12) and recall that $\mathcal{P}_N R$ is the projection of the semiclassical solution $R$ onto the first $N$ Hermite modes. From the triangle inequality, the error estimate is obtained by bounding one term corresponding to the projection error and a second one coming from the error between $R^N$ and $\mathcal{P}_N R$, that is,

$$ (3.8) \qquad \|R^N(t) - R(t)\|_{L^2} \le \|R^N(t) - \mathcal{P}_N R(t)\|_{L^2} + \|\mathcal{P}_N R(t) - R(t)\|_{L^2}. $$

We first evaluate the propagation error between the solution $R^N$ to (2.12) and the projection $\mathcal{P}_N R$, where $R$ is the solution to semiclassical limit equation (1.13). For all $k \in \{0, \dots, N\}$,

$$ \partial_t (R_k^N - R)_k = -i \left( \sqrt{\frac{k}{2}} \mathcal{D} (R^N - R)_{k-1} + \sqrt{\frac{k+1}{2}} \mathcal{D}^* (R^N - R)_{k+1} \right), $$

and we deduce from the Duhamel formula that

$$ R^N(t) - \mathcal{P}_N R(t) = S_0(t) (R^N - \mathcal{P}_N R)(0) + i \sqrt{\frac{N+1}{2}} \int_0^t S_0(t-s) \mathcal{D}^* R_{N+1}(s) ds $$

where $S_0(t)$ is the one-parameter group (1.14). Hence, using the conservation of the $L^2$ norm by $S_0(t)$, we obtain that

$$ (3.9) \qquad \|R^N(t) - \mathcal{P}_N R(t)\|_{L^2} = \left( \sum_{k=0}^{N} \int_{\mathbb{R}} |R_k^N(t|x) - R_k(t|x)|^2 dx \right)^{1/2} \\ \le \sqrt{\frac{N+1}{2}} \int_0^t \left( \int_{\mathbb{R}} |\mathcal{D}^* R_{N+1}(\tau|x)|^2 dx \right)^{1/2} d\tau. $$

Let us estimate the right-hand side using the Parseval identity: since $\mathcal{D}^*$ acts on the variable $x$, while $\mathcal{P}_N$ acts on the variable $y$, these two operators commute, and applying Lemma 2.3 shows that

$$ \begin{align*} \int_{\mathbb{R}} |\mathcal{D}^* R_{N+1}(\tau|x)|^2 dx &\le \iint_{\mathbb{R}^2} |\mathcal{D}^* (R - \mathcal{P}_N R)(\tau|x, y)|^2 dx dy \\ &\le \frac{1}{(2N+3)^{2p}} \iint_{\mathbb{R}^2} |(y^2 - \partial_y^2)^p \mathcal{D}^* R(\tau|x, y)|^2 dx dy. \end{align*} $$
Then, using the mean-value inequality (3.5)

$$ \begin{aligned} |(y^2 - \partial_y^2)^p \mathcal{D}^* R(\tau |x, y)|^2 &\le 2 |(y^2 - \partial_y^2)^p \partial_x R(\tau |x, y)|^2 + 2 |V'(x)|^2 |(y^2 - \partial_y^2)^p R(\tau |x, y)|^2 \\ &\le 2 |(y^2 - \partial_y^2)^p \partial_x R(\tau |x, y)|^2 \\ &\quad + 4 (|V'(0)|^2 + x^2 \|V''\|_{L^\infty}^2) |(y^2 - \partial_y^2)^p R(\tau |x, y)|^2. \end{aligned} $$

Hence, there exists a constant $C_p > 0$, depending only on $p \ge 1$ such that, for all $\tau \ge 0$,

$$ \iint_{\mathbb{R}^2} |(y^2 - \partial_y^2)^p \mathcal{D}^* R(\tau |x, y)|^2 dx dy \le C_p \max(1, |V'(0)|, \|V''\|_{L^\infty}) \mathcal{N}_{2p+1}(\tau)^2, $$

where $\mathcal{N}_{2p+1}(\tau)$ is given in (1.15). Therefore, we conclude by applying Proposition 3.1 that

$$ \int_0^t \left( \int_{\mathbb{R}} |\mathcal{D}^* R_{N+1}(\tau |x)|^2 dx \right)^{1/2} d\tau $$

$$ (3.10) \qquad \le \frac{1}{(2N + 3)^p} \int_0^t \left( \iint_{\mathbb{R}^2} |(y^2 - \partial_y^2)^p \mathcal{D}^* R(\tau |x, y)|^2 dx dy \right)^{1/2} d\tau $$

$$ (3.11) \qquad \le \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N + 3)^p} \mathcal{N}_{2p+1}(0), $$

where $C_{2p+1}[V]$ is given in (3.1). Finally, gathering together (3.9)-(3.10), it yields the following inequality

$$ (3.12) \qquad \|R^N(t) - \mathcal{P}_N R(t)\|_{L^2} \le \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N + 3)^{p-1/2}} \frac{C_{2p+1}[V]}{C_{2p+1}[V]} \mathcal{N}_{2p+1}(0). $$

Now, it remains to evaluate the last term on the right-hand side in (3.8), which corresponds to the projection error. Thanks to Lemma 2.3, we have

$$ \|R(t) - \mathcal{P}_N R(t)\|_{L^2(\mathbb{R}^2)} \le \frac{1}{(2N + 3)^p} \|(y^2 - \partial_y^2)^p R(t)\|_{L^2(\mathbb{R}^2)}, $$

hence using the propagation of regularity in Proposition 3.1, we finally get

$$ (3.13) \qquad \|R(t) - \mathcal{P}_N R(t)\|_{L^2(\mathbb{R}^2)} \le \frac{K_p}{(2N + 3)^p} e^{C_{2p}[V]t} \mathcal{N}_{2p}(0), $$

where $K_p > 0$ is a constant depending only on $p$. Gathering together the former results (3.12) and (3.13), we get the desired inequality

$$ \begin{aligned} \|R^N(t) - R(t)\|_{L^2} &\le \frac{K_p}{(2N + 3)^p} e^{C_{2p}[V]t} \mathcal{N}_{2p}(0) \\ &\quad + \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N + 3)^{p-1/2}} \mathcal{N}_{2p+1}(0). \end{aligned} $$

## 4. ANALYSIS OF THE VON NEUMANN EQUATION

We now treat the von Neumann equation written in Weyl's variables (1.6) and focus on the contribution of the non-local term $\mathcal{E}^h$ in (1.11). We apply the same strategy as before, but the proof of the propagation of regularity now becomes more intricate.
**4.1. Propagation of regularity on the von Neumann equation.** It is worth mentioning that the latter result only requires that the solution to the von Neumann equation (1.6) is $R^h(t) \in L^2$ since we use the regularity of the solution $R$ of (1.13) proved in Proposition 3.1. However, when the initial datum $R_{\text{in}}^h$ is smooth, we prove that for any $m \ge 0$, the quantity $N_m(t)$ is also propagated for the solution $R^h$ to the von Neumann equation (1.6).

**Proposition 4.1.** Let $m \ge 0$ and $V$ be the potential such that (1.16)-(1.17) hold, and suppose that for all integer $a, b, \alpha$ and $\beta$ such that $a+b+\alpha+\beta \le m$, the assumption (2.11) on the initial data $R_{\text{in}}^h$ is satisfied. Then the solution $R^h$ to (1.6) is such that for all $t \ge 0$,

$$N_m(t) \le N_m(0) e^{C_m[V]t},$$

where $N_m(t)$ is defined in (1.15) and applied to $R^h$ while $C_m[V] > 0$ is a constant depending only on $m$ and on the potential $V$.

*Proof.* This proof follows the same strategy as the one presented in Proposition 3.1 but is rather more technical. See Appendix A for details. □

**4.2. Proof of Theorem 2.1.** Let us first show that the $L^2$ estimate given in (1.10) is preserved when the number of modes $N$ is finite in the approximation of the von Neumann equation using Hermite polynomials.

**Proposition 4.2.** For a given $m \ge 0$, assume that $V$ satisfies (1.16)-(1.17) while $R_{\text{in}}^h$ satisfies (2.11) and consider $R^{h,N}$ the solution to (2.10). Then

(i) for all $k, l \in \{0, \dots, N\}$ and $t \ge 0$, we have

$$(4.1) \qquad R_k^{h,N}(t) = (-1)^k \overline{R_k^{h,N}}(t);$$

(ii) for all $t \ge 0$, we have

$$(4.2) \qquad \|R^{h,N}(t)\|_{L^2} = \|\mathcal{P}_N R_{\text{in}}^h\|_{L^2}^2.$$

*Proof.* To prove (4.1), use (2.7) and proceed as in the proof of Proposition 3.2 to show that, if for all $k \in \{0, \dots, N\}$, $R_k^{h,N}$ is solution to (2.10), then $(-1)^k \overline{R_k^{h,N}}$ also solves (2.10), which gives the expected result.

We next prove the propagation of the $L^2$ norm for $R^{h,N}$ again as in Proposition 3.2. Thus, we have

$$\frac{1}{2} \frac{d}{dt} \|R^{h,N}(t)\|_{L^2}^2 = \mathcal{I}_1(t) + \mathcal{I}_2(t),$$

with

$$\left\{
\begin{aligned}
\mathcal{I}_1(t) &:= \operatorname{Re} \left( i \sum_{k=0}^{N} \left\langle \sqrt{\frac{k}{2}} \mathcal{D} R_{k-1}^{h,N} + \sqrt{\frac{k+1}{2}} \mathcal{D}^* R_{k+1}^{h,N}, R_k^{h,N} \right\rangle \right), \\
\mathcal{I}_2(t) &:= \operatorname{Re} \left( i \sum_{k,l=0}^{N} \langle \mathcal{E}_{k,l}^h R_l^{h,N}, R_k^{h,N} \rangle \right),
\end{aligned}
\right.$$

where $\langle \cdot, \cdot \rangle$ denotes the standard inner product (1.1) on $L^2(\mathbb{R}; \mathbb{C})$. On the one hand, observing that $\mathcal{I}_1(t)$ is the same term as the one of the right-hand side in the proof of Proposition 3.2, we get that $\mathcal{I}_1(t) = 0$. On the other hand, using that $\mathcal{E}_{k,l}^h \in \mathbb{R}$ is symmetric and satisfies (2.7) while $R^{h,N}$
verifies (4.1), we get that for all $k, l \in \{0, \dots, N\}$

$$ \begin{aligned} \langle \mathcal{E}_{k,l}^h R_l^{h,N}, R_k^{h,N} \rangle &= \int_{\mathbb{R}} \mathcal{E}_{k,l}^h R_l^{h,N} \overline{R_k^{h,N}} dx \\ &= \int_{\mathbb{R}} \left( (-1)^{l+k} \mathcal{E}_{k,l}^h \right) \left( (-1)^l \overline{R_l^{h,N}} \right) \left( (-1)^k R_k^{h,N} \right) dx \\ &= - \int_{\mathbb{R}} \mathcal{E}_{l,k}^h R_k^{h,N} \overline{R_l^{h,N}} dx = - \langle \mathcal{E}_{l,k}^h R_k^{h,N}, R_l^{h,N} \rangle. \end{aligned} $$

Therefore, $\langle \mathcal{E}_{k,l}^h R_l^{h,N}, R_k^{h,N} \rangle$ is antisymmetric in $(k, l)$ so that

$$ \sum_{k,l=0}^{N} \langle \mathcal{E}_{k,l}^{h} R_{l}^{h,N}, R_{k}^{h,N} \rangle = 0, $$

hence $\mathcal{I}_2(t)$ also vanishes, leading to the conservation of the $L^2$ norm for $R^{h,N}$. $\square$

We are now ready to complete the proof of Theorem 2.1 and follow the lines of the arguments provided in Section 3.2 for Theorem 2.2. Here we will only focus on the additional term $\mathcal{E}^h$ when evaluating the propagation error between the solution $R^{h,N}$ to (2.10) and the projection $\mathcal{P}_N R^h$, where $R^h$ is the solution to the von Neumann equation(1.6). We need to estimate the nonlocal term involving $\mathcal{E}_{k,l}^h$, namely,

$$ \mathcal{I}_N(t) = \sum_{k=0}^{N} \sum_{l>N} \int_{\mathbb{R}} \mathcal{E}_{k,l}^h(x) R_l^h(t|x) \overline{(R_k^h - R_k^{h,N})}(t|x) dx. $$

Let us re-write $\mathcal{I}_N(t)$ as

$$ \begin{aligned} \mathcal{I}_N(t) &= \iint \mathcal{E}^h(x, y) \overline{(\mathcal{P}_N R^h - R^{h,N})}(t|x, y) (I - \mathcal{P}_N) R^h(t|x, y) dy dx \\ &\le \|\mathcal{E}^h (\mathcal{P}_N R^h(t) - R^{h,N}(t))\|_{L^2} \|R^h(t) - \mathcal{P}_N R^h(t)\|_{L^2_{x,y}}. \end{aligned} $$

Using that

$$ |\mathcal{E}^h(x, y)| = \left| \frac{V(x + \frac{\hbar}{2} y) - V(x - \frac{\hbar}{2} y)}{\hbar} - V'(x) y \right| \le \frac{\hbar^2}{24} |y|^3 \|V^{(3)}\|_{L^\infty}, $$

we have

$$ \begin{aligned} \|\mathcal{E}^h (R^{h,N}(t) - \mathcal{P}_N R^h(t))\|_{L^2_{x,y}}^2 &= \iint |\mathcal{E}^h(x,y)|^2 |(R^{h,N} - \mathcal{P}_N R^h)(t|x,y)|^2 dy dx \\ &\le \left(\frac{\hbar^2}{24} \|V^{(3)}\|_{L^\infty}\right)^2 \sum_{j,k=0}^{N} \int (R_k^{h,N} - R_k^h) \overline{(R_j^{h,N} - R_j^h)}(t|x) dx \int |y|^6 \Phi_k(y) \Phi_j(y) dy. \end{aligned} $$

Now using repeatedly that

$$ y \Phi_k(y) = \sqrt{\frac{k+1}{2}} \Phi_{k+1}(y) + \sqrt{\frac{k}{2}} \Phi_{k-1}(y), \quad k \ge 0, $$

we have

$$ y^3 \Phi_k(y) = \alpha_k \Phi_{k+3}(y) + \beta_k \Phi_{k+1}(y) + \gamma_k \Phi_{k-1}(y) + \delta_k \Phi_{k-3}(y), $$
with

$$
\left\{
\begin{aligned}
\alpha_k &= \sqrt{\frac{(k+1)(k+2)(k+3)}{8}}, \\
\beta_k &= \frac{3(k+1)}{2} \sqrt{\frac{k+1}{2}}, \\
\gamma_k &= \frac{3k}{2} \sqrt{\frac{k}{2}}, \\
\tau_k &= \sqrt{\frac{k(k-1)(k-2)}{8}},
\end{aligned}
\right.
$$

so that

$$
\alpha_k, \beta_k, \gamma_k, \tau_k \le \frac{3}{2\sqrt{2}} (k+3)^{3/2}.
$$

Hence, we obtain

$$
\begin{align*}
\int_{\mathbb{R}} |y|^6 \Phi_j(y) \Phi_k(y) \, dy &= \int_{\mathbb{R}} |y|^3 \Phi_j(y) |y|^3 \Phi_k(y) \, dy \\
&= (\alpha_j^2 + \beta_j^2 + \gamma_j^2 + \tau_j^2) \delta_{j,k} \\
&\quad + (\alpha_k \beta_j + \beta_k \gamma_j + \gamma_k \tau_j) \delta_{j,k+2} + (\alpha_k \gamma_j + \beta_k \tau_j) \delta_{j,k+4} + \alpha_k \tau_j \delta_{j,k+6} \\
&\quad + (\alpha_j \beta_k + \beta_j \gamma_k + \gamma_j \tau_k) \delta_{j+2,k} + (\alpha_j \gamma_k + \beta_j \tau_k) \delta_{j+4,k} + \alpha_j \tau_k \delta_{j+6,k},
\end{align*}
$$

and for all $0 \le j, k \le N$, denoting $\max(j, k)$ by $j \lor k$,

$$
a_j b_k \le \frac{9}{8} (j \lor k + 9)^3, \quad \forall a, b \in \{\alpha, \beta, \gamma, \tau\}.
$$

Finally, we get

$$
\begin{align*}
\| \mathcal{E}^h (R^{h,N}(t) - \mathcal{P}_N R^h(t)) \|_{L^2}^2 &\le \frac{1}{8} \left( \frac{\hbar^2}{8} \|V^{(3)}\|_{L^\infty} \right)^2 \times \\
&\quad \left( \sum_{j=0}^{N} (j+9)^3 \int |R_j^{h,N} - R_j^h|^2 (t|x) \,dx \right. \\
&\qquad + 2 \sum_{j=0}^{N-2} (j+11)^3 \int |R_j^{h,N} - R_j^h| |R_{j+2}^{h,N} - R_{j+2}^h|(t|x) \,dx \\
&\qquad + 2 \sum_{j=0}^{N-4} (j+13)^3 \int |R_j^{h,N} - R_j^h| |R_{j+4}^{h,N} - R_{j+4}^h|(t|x) \,dx \\
&\qquad \left. + 2 \sum_{j=0}^{N-6} (j+15)^3 \int |R_j^{h,N} - R_j^h| |R_{j+6}^{h,N} - R_{j+6}^h|(t|x) \,dx \right) \\
&\le \frac{7(N+9)^3}{2} \left( \frac{\hbar^2}{8} \|V^{(3)}\|_{L^\infty} \right)^2 \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2}^2.
\end{align*}
$$
Now we are ready to bound the term $\mathcal{I}_N(t)$ by applying the latter estimate and Lemma 2.3, which gives for $N \ge 6$,

$$
\begin{align*}
\mathcal{I}_N(t) &\le \frac{\hbar^2}{4} \|V^{(3)}\|_{L^\infty} (N+9)^{3/2} \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2} \|R^h(t) - \mathcal{P}_N R^h(t)\|_{L^2} \\
&\le \frac{\hbar^2}{4} \|V^{(3)}\|_{L^\infty} \frac{(N+9)^{3/2}}{(2N+3)^{p+1}} \|(y^2 - \partial_y^2)^{p+1} R^h(t)\|_{L^2} \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2} \\
&\le \frac{\hbar^2}{4} \|V^{(3)}\|_{L^\infty} \frac{1}{(2N+3)^{p-1/2}} \|(y^2 - \partial_y^2)^{p+1} R^h(t)\|_{L^2} \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2} \\
&\le \frac{\hbar^2}{4} \|V^{(3)}\|_{L^\infty} \frac{K_p \mathcal{N}_{2p+2}(t)}{(2N+3)^{p-1/2}} \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2},
\end{align*}
$$

where $K_p$ is a constant depending on $p$ only and $\mathcal{N}_{2p+2}(t)$ is evaluated in Proposition 4.1.

Finally, gathering together this latter result and those obtained in the proof of Theorem 2.2 on the contribution of $\mathcal{D}^*R_{N+1}^h$, we are lead to the following inequality

$$
\begin{aligned}
\frac{d}{dt} \|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2} &\le \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})}}{(2N+3)^{p-1/2}} \mathcal{N}_{2p+1}(t) \\
&\quad + \frac{\hbar^2}{4} \|V^{(3)}\|_{L^\infty} \frac{2 K_p \mathcal{N}_{2p+2}(t)}{(2N+3)^{p-1/2}}.
\end{aligned}
$$

Thus, we have

$$
\begin{aligned}
\|R^{h,N}(t) - \mathcal{P}_N R^h(t)\|_{L^2} &\le \frac{\sqrt{C_p \max(1, |V'(0)|, \|V''\|_{L^\infty})} e^{C_{2p+1}[V]t} - 1}{(2N+3)^{p-1/2}} \frac{C_{2p+1}[V]}{\mathcal{N}_{2p+1}(0)} \\
&\quad + \frac{2 K_p \hbar^2}{4 (2N+3)^{p-1/2}} \|V^{(3)}\|_{L^\infty} \frac{e^{C_{2p+2}[V]t} - 1}{C_{2p+2}[V]} \mathcal{N}_{2p+2}(0).
\end{aligned}
$$

Finally from this latter inequality and the estimate on the norm of the operator $I - \mathcal{P}_N$ in Lemma 2.3, on the operator $I - \mathcal{P}_N$, we get the sought inequality.

5. NUMERICAL SIMULATIONS

We now perform a few numerical simulations to illustrate the spectral accuracy of the proposed Hermite method for the approximation of the solution to the von Neumann equation (1.6) in the case of a harmonic potential $V$. We refer to [14] for more extensive numerical simulations for various potential.

Thus, we consider the isotropic harmonic potential

$$V(x) := \frac{x^2}{2}$$

in order to benchmark the convergence rate of the proposed method. In this situation, for all $(t, x) \in [0, +\infty) \times \mathbb{R}^d$, we define the Wigner function $W(t, x, \cdot)$ by taking the Fourier transform of the function $R^h(t|x, \cdot)$

$$W(t, x, \xi) := \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} R^h(t|x, y) e^{-i y \cdot \xi} dy.$$

Then the function $W$ is formally a solution of the so-called Wigner equation

$$\partial_t W + \xi \partial_x W - \partial_x V \partial_\xi W = 0.$$
The exact solution is $W(t, x, \xi) = W_0(x(t), \xi(t))$, where

$$ egin{cases} x(t) = \cos(t) x(0) - \sin(t) \xi(0), \\ \xi(t) = \sin(t) x(0) + \cos(t) \xi(0). \end{cases} $$

Consider a Gaussian wavepacket

$$ W_{\text{in}}(x, \xi) := \frac{1}{2\pi} \exp \left( -\frac{(x-5)^2 + \xi^2}{2} \right) $$

as initial condition, which returns to the initial state at the final time $T = 2\pi$. It is worth mentioning that our Hermite spectral method is well suited to compute the Wigner transform. Indeed, $W$ is the Fourier transform of $R^h$ in the variable $y$, while Hermite functions $\Phi_k$ are eigenfunctions of the Fourier transform :

$$ \mathcal{F}(\Phi_k)(\xi) = (-i)^k \Phi_k(\xi), $$

where $\mathcal{F}$ represents the Fourier transform.

Here, we choose the $x$ domain as $[-8, 8]$ and mainly focus on the convergence with respect to $N$, hence the time step $\Delta t$ and the mesh size $\Delta x$ are chosen sufficiently small such that the error with respect to $N$ dominates. Numerical errors are presented in Table 5.1, where the numerical error is given as

$$ \mathcal{E}(N) := \max_{t \in [0, 2\pi]} \left( \sum_{k=0}^{N} \| R_k^{h,N}(t) - R_k^h(t^n) \|_{L_x^2}^2 \right)^{1/2} . $$

<figure><img src="image_3.png" /><figcaption>FIGURE 5.1. $L^2$ error in log scale.</figcaption></figure>

<table><caption>TABLE 5.1. $L^2$ error with respect to the number of Hermite modes $N$.</caption><thead><tr><th>Number of modes <i>N</i></th><th>Discrete <i>L</i><sup>2</sup> error <i>E</i>(<i>N</i>)</th><th>Order of accuracy</th></tr></thead><tbody><tr><td>10</td><td>6.0072 &times; 10<sup>-01</sup></td><td>X</td></tr><tr><td>20</td><td>4.0877 &times; 10<sup>-01</sup></td><td>0.5</td></tr><tr><td>30</td><td>8.3406 &times; 10<sup>-02</sup></td><td>3.9</td></tr><tr><td>40</td><td>6.0200 &times; 10<sup>-05</sup></td><td>25</td></tr><tr><td>50</td><td>1.8176 &times; 10<sup>-07</sup></td><td>36</td></tr><tr><td>60</td><td>7.0635 &times; 10<sup>-12</sup></td><td>43</td></tr></tbody></table>

## 6. CONCLUSION

In conclusion, we have introduced a new approach to discretize the von Neumann equation in the semiclassical limit. By using Weyl's variables and a truncated Hermite expansion of the density operator, we are able to address the stiffness due to the semiclassical scaling of the equation. Our method allows for error estimates by leveraging the propagation of regularity on the exact solution. Therefore, we believe that our asymptotic preserving numerical approximation, coupled with Hermite polynomials, represents an interesting approach for accurately solving the von Neumann equation in the semiclassical regime.
APPENDIX A. PROOF OF PROPOSITION 4.1

The case $m = 0$ corresponds to the conservation of the $L^2$ norm

$$ \|R^h(t)\|_{L^2} = \|R_{\text{in}}^h\|_{L^2}, \quad t \ge 0. $$

Let us begin with the propagation of quantities corresponding to phase-space derivatives in the Wigner equation, that is to say $\partial_x R^h$ and $y R^h$. Later, it will be convenient to use the notation

$$ L := \max(1, \|V''\|_{L^\infty}). $$

On the one hand, the equation for the $x$-derivative of $R^h$ is given by

$$ \begin{aligned} \partial_t \partial_x R^h(t|x, y) &= i \partial_x \partial_y \partial_x R^h(t|x, y) - i \partial_x \left( \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} R^h(t|x, y) \right) \\ &= i \partial_x \partial_y \partial_x R^h(t|x, y) - i \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} \partial_x R^h(t|x, y) \\ &\quad - i \frac{V'(x + \frac{1}{2}\hbar y) - V'(x - \frac{1}{2}\hbar y)}{\hbar} R^h(t|x, y), \end{aligned} $$

so that

$$ \partial_x R^h(t|x, y) = S_h(t) \partial_x R^h(0|x, y) - i \int_0^t S_h(t-s) \frac{V'(x + \frac{1}{2}\hbar y) - V'(x - \frac{1}{2}\hbar y)}{\hbar} R^h(s|x, y) \, ds, $$

where $S_h(t)$ represents the one-parameter group (1.9) associated to (1.6). Since $S_h(t)$ is an $L^2_{x,y}$ isometry,

$$ \begin{aligned} \|\partial_x R^h(t)\|_{L^2} &\le \|\partial_x R^h(0)\|_{L^2} + \int_0^t \left\| \frac{V'(x + \frac{1}{2}\hbar y) - V'(x - \frac{1}{2}\hbar y)}{\hbar y} y R^h(s|x, y) \right\|_{L^2} \, ds \\ &\le \|\partial_x R^h(0)\|_{L^2} + \|V''\|_{L^\infty} \int_0^t \|y R^h(s)\|_{L^2} \, ds. \end{aligned} $$

Similarly,

$$ \begin{aligned} \partial_t (y R^h(t|x, y)) &= i \partial_x \partial_y (y R^h(t|x, y)) - i \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} (y R^h(t|x, y)) \\ &\quad - i \partial_x R^h(t|x, y), \end{aligned} $$

so that, arguing as above,

$$ \|y R^h(t)\|_{L^2} \le \|y R^h(0)\|_{L^2} + \int_0^t \|\partial_x R^h(s)\|_{L^2} \, ds. $$

Summarizing, we have proved that

$$ \begin{aligned} \|\partial_x R^h(t)\|_{L^2} &+ \|y R^h(t)\|_{L^2} - \|\partial_x R^h(0)\|_{L^2} - \|y R^h(0)\|_{L^2} \\ &\le \int_0^t \left( \|V''\|_{L^\infty} \|y R^h(s)\|_{L^2} + \|\partial_x R^h(s)\|_{L^2} \right) \, ds \\ &\le L \int_0^t \left( \|\partial_x R^h(s)\|_{L^2} + \|y R^h(s)\|_{L^2} \right) \, ds. \end{aligned} $$

Applying Gronwall's lemma, implies that

$$ (A.1) \qquad \|\partial_x R^h(t)\|_{L^2} + \|y R^h(t)\|_{L^2} \le \left( \|\partial_x R^h(0)\|_{L^2} + \|y R^h(0)\|_{L^2} \right) e^{Lt}. $$
Here we observe that this estimate is uniform in (in fact independent of) $\hbar$. Next we discuss the propagation of quantities corresponding to moments in the Wigner equation, that is, $xR^\hbar$ and $\partial_y R^\hbar$. On the one hand, the equation for $xR^\hbar(t|x, y)$ is

$$ \partial_t(xR^\hbar(t|x, y)) = i \partial_x \partial_y(xR^\hbar(t|x, y)) - i \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} (xR^\hbar(t|x, y)) \\ - i \partial_y R^\hbar(t|x, y), $$

so that, arguing as above

$$ \|x R^\hbar(t)\|_{L^2} \le \|x R^\hbar(0)\|_{L^2} + \int_0^t \|\partial_y R^\hbar(s)\|_{L^2} ds. $$

On the other hand, the equation for $\partial_y R^\hbar$ is

$$ \begin{aligned} \partial_t \partial_y R^\hbar(t|x, y) &= i \partial_x \partial_y \partial_y R^\hbar(t|x, y) - i \partial_y \left( \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} R^\hbar(t|x, y) \right) \\ &= i \partial_x \partial_y \partial_y R^\hbar(t|x, y) - i \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} \partial_y R^\hbar(t|x, y) \\ &\quad - \frac{i}{2} \left( V' \left( x + \frac{1}{2}\hbar y \right) + V' \left( x - \frac{1}{2}\hbar y \right) \right) R^\hbar(t|x, y). \end{aligned} $$

Here it is worth mentioning that there are several options for the treatment of the last term on the last right-hand side. One could of course assume that $V'$ is bounded, but this is not very satisfying as it leaves out the case of a harmonic oscillator, where $V$ is quadratic in $x$. In the latter case, one can assume instead that $V''$ is bounded. This is the simpler case, allowing us to include a confining quadratic potential if needed. Besides, since $V(x) \to +\infty$ as $|x| \to +\infty$, it has a minimum point at $x_0 \in \mathbb{R}$. Without loss of generality, assume henceforth that $x_0 = 0$, so that $V'(0) = 0$. Therefore, by the mean value inequality

$$ |V'(z) - V'(0)| \le \|V''\|_{L^\infty} |z|, $$

and hence, having recast the last equality as

$$ \begin{aligned} \partial_t \partial_y R^\hbar(t|x, y) &= i \partial_x \partial_y \partial_y R^\hbar(t|x, y) - i \frac{V(x + \frac{1}{2}\hbar y) - V(x - \frac{1}{2}\hbar y)}{\hbar} \partial_y R^\hbar(t|x, y) \\ &\quad - i V'(x) R^\hbar(t|x, y) - i \frac{V'(x + \frac{1}{2}\hbar y) + V'(x - \frac{1}{2}\hbar y) - 2V'(x)}{2} R^\hbar(t|x, y), \end{aligned} $$

we conclude as above that

$$ \|\partial_y R^\hbar(t)\|_{L^2} \le \|\partial_y R^\hbar(0)\|_{L^2} + \|V''\|_{L^\infty} \int_0^t \left( \|x R^\hbar(s)\|_{L^2} + \frac{\hbar}{2} \|y R^\hbar(s)\|_{L^2} \right) ds. $$

Combining this inequality with the one for $x R^\hbar(t|x, y)$ leads to

$$ \begin{aligned} & \|x R^\hbar(t)\|_{L^2} + \|\partial_y R^\hbar(t)\|_{L^2} - \|x R^\hbar(0)\|_{L^2} - \|\partial_y R^\hbar(0)\|_{L^2} \\ & \le \int_0^t \left( \|V''\|_{L^\infty} \|x R^\hbar(s)\|_{L^2} + \|\partial_y R^\hbar(s)\|_{L^2} \right) ds + \frac{\hbar}{2} \|V''\|_{L^\infty} \int_0^t \|y R^\hbar(s)\|_{L^2} ds. \end{aligned} $$

The last term on the right-hand side is bounded as in (A.1), by

$$ \left( \|\partial_x R^\hbar(0)\|_{L^2} + \|y R^\hbar(0)\|_{L^2} \right) e^{t \max(1, \|V''\|_{L^\infty})}. $$
Thus

$$
\begin{align*}
& \|x R^h(t)\|_{L^2} + \|\partial_y R^h(t)\|_{L^2} - \|x R^h(0)\|_{L^2} - \|\partial_y R^h(0)\|_{L^2} \\
& \le L \int_0^t (\|x R^h(s)\|_{L^2} + \|\partial_y R^h(s)\|_{L^2}) \, ds + \frac{\hbar}{2} (\|\partial_x R^h(0)\|_{L^2} + \|y R^h(0)\|_{L^2}) (e^{Lt} - 1)
\end{align*}
$$

and by Gronwall's lemma, we conclude that

$$
\begin{equation}
\begin{aligned}
\|x R^h(t)\|_{L^2} + \|\partial_y R^h(t)\|_{L^2} \le {}& \left( \|x R^h(0)\|_{L^2} + \|\partial_y R^h(0)\|_{L^2} \right) e^{Lt} \\
& + \frac{\hbar}{2} \left( \|\partial_x R^h(0)\|_{L^2} + \|y R^h(0)\|_{L^2} \right) (e^{Lt} - 1) e^{Lt}.
\end{aligned}
\tag{A.2}
\end{equation}
$$

Unlike (A.1), this bound is not independent of $\hbar$, but its $\hbar$-dependence is not singular as $\hbar \to 0$.
Together with (A.1), the inequality (A.2) is easily transformed into a uniform in $\hbar$ bound for
$\hbar \in (0, 2]$, say:

$$
\begin{equation}
\begin{aligned}
\mathcal{N}_1(t) &\le \mathcal{N}_1(0) e^{Lt} + \frac{\hbar}{2} \left( \|\partial_x R^h(0)\|_{L^2} + \|y R^h(0)\|_{L^2} \right) (e^{Lt} - 1) e^{Lt} \\
(A.3) \qquad &\le \mathcal{N}_1(0) e^{Lt} + \frac{\hbar}{2} \mathcal{N}_1(0) (e^{Lt} - 1) e^{Lt} \\
&\le \mathcal{N}_1(0) e^{2Lt},
\end{aligned}
\end{equation}
$$

which yields the desired result for $n = 1$.

We now proceed in the same way, for each integer $a$, $b$, $\alpha$ and $\beta$ such that $a+b+\alpha+\beta=n$ and
seek an equation for

$$
x^a \partial_x^b y^\alpha \partial_y^\beta R^h(t | x, y).
$$

Starting from

$$
\partial_t R^h(t|x,y) = i \partial_x \partial_y R^h(t|x,y) - i \frac{V(x+\frac{1}{2}\hbar y) - V(x-\frac{1}{2}\hbar y)}{\hbar} R^h(t|x,y),
$$

we see that

$$
\begin{align*}
\partial_t (x^a \partial_x^b y^\alpha \partial_y^\beta R^h(t|x,y))
  &= i \partial_x \partial_y (x^a \partial_x^b y^\alpha \partial_y^\beta R^h(t|x,y)) \\
  &\quad - i \frac{V(x+\frac{1}{2}\hbar y) - V(x-\frac{1}{2}\hbar y)}{\hbar} (x^a \partial_x^b y^\alpha \partial_y^\beta R^h(t|x,y)) + i \mathcal{T}_1 - i \mathcal{T}_2,
\end{align*}
$$

where

$$
\mathcal{T}_1 := \left[ x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y \right] R^h(t|x,y), \\
\mathcal{T}_2 := \left[ x^a \partial_x^b y^\alpha \partial_y^\beta, \frac{V(x+\frac{1}{2}\hbar y) - V(x-\frac{1}{2}\hbar y)}{\hbar} \right] R^h(t|x,y).
$$

Hence, the idea is to treat $\mathcal{T}_j$ as source terms for $j = 1, 2$, to find

$$
\|x^a \partial_x^b y^\alpha \partial_y^\beta R^h(t|x,y)\|_{L^2} \le \|x^a \partial_x^b y^\alpha \partial_y^\beta R^h(0|x,y)\|_{L^2} \\
+ \int_0^t (\|\mathcal{T}_1(s|\cdot, \cdot)\|_{L^2} + \|\mathcal{T}_2(s|\cdot, \cdot)\|_{L^2}) ds.
$$
First, we consider $\mathcal{T}_1$ and use that

$$
\begin{align*}
[x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x \partial_y] &= [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_x] \partial_y + \partial_x [x^a \partial_x^b y^\alpha \partial_y^\beta, \partial_y], \\
&= [x^a, \partial_x] \partial_x^b y^\alpha \partial_y^\beta \partial_y + \partial_x x^a \partial_x^b [y^\alpha, \partial_y] \partial_y^\beta \\
&= -ax^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} - \alpha \partial_x x^a \partial_x^b y^{\alpha-1} \partial_y^\beta \\
&= -a x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} - \alpha x^a \partial_x^{b+1} y^{\alpha-1} \partial_y^\beta - \alpha [\partial_x, x^a] \partial_x^b y^{\alpha-1} \partial_y^\beta \\
&= -a x^{a-1} \partial_x^b y^\alpha \partial_y^{\beta+1} - \alpha x^a \partial_x^{b+1} y^{\alpha-1} \partial_y^\beta - \alpha a x^{a-1} \partial_x^b y^{\alpha-1} \partial_y^\beta.
\end{align*}
$$

In the first two terms on the right-hand side, the sum of exponents is

$$a - 1 + b + \alpha + \beta + 1 = a + b + 1 + \alpha - 1 + \beta = n,$$

equal to that in $x^a \partial_x^b y^\alpha \partial_y^\beta$, whereas the last term has a sum of exponents

$$a - 1 + b + \alpha - 1 + \beta = n - 2,$$

less than that in $x^a \partial_x^b y^\alpha \partial_y^\beta$. Hence we obtain

$$\|\mathcal{T}_1(s|\cdot, \cdot)\|_{L^2} \le 2n \mathcal{N}_n(s) + n^2 \mathcal{N}_{n-2}(s) \le (2n + n^2) \mathcal{N}_n(s).$$

Now for the second term $\mathcal{T}_2$, we have

$$\left[ x^a \partial_x^b y^\alpha \partial_y^\beta, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right] = x^a \left[ \partial_x^b, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right] y^\alpha \partial_y^\beta \\ + x^a \partial_x^b y^\alpha \left[ \partial_y^\beta, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right].$$

By Leibniz's formula

$$\left[ \partial_x^b, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right] = \sum_{j=1}^{b} \binom{b}{j} \frac{V^{(j)}(x + \frac{\hbar}{2}y) - V^{(j)}(x - \frac{\hbar}{2}y)}{\hbar} \partial_x^{b-j},$$

$$\left[ \partial_y^\beta, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right] = \sum_{k=1}^{\beta} \binom{\beta}{k} \frac{V^{(k)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \hbar^k \partial_y^{\beta-k},$$
so that

$$
\begin{align*}
& \left[ x^a \partial_x^b y^\alpha \partial_y^\beta, \frac{V(x + \frac{\hbar}{2}y) - V(x - \frac{\hbar}{2}y)}{\hbar} \right] \\
&= \sum_{j=1}^b \binom{b}{j} \frac{V^{(j)}(x + \frac{\hbar}{2}y) - V^{(j)}(x - \frac{\hbar}{2}y)}{\hbar} x^a \partial_x^{b-j} y^\alpha \partial_y^\beta \\
&\quad + \sum_{k=1}^\beta \binom{\beta}{k} \frac{V^{(k)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \hbar^k x^a \partial_x^b y^\alpha \partial_y^{\beta-k} \\
&\quad + \sum_{k=1}^\beta \binom{\beta}{k} x^a \left[ \partial_x^b, \frac{V^{(k)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \right] \hbar^k y^\alpha \partial_y^{\beta-k} \\
&= b \frac{V'(x + \frac{\hbar}{2}y) - V'(x - \frac{\hbar}{2}y)}{\hbar y} x^a \partial_x^{b-1} y^{\alpha+1} \partial_y^\beta \\
&\quad + \sum_{j=2}^b \binom{b}{j} \frac{V^{(j)}(x + \frac{\hbar}{2}y) - V^{(j)}(x - \frac{\hbar}{2}y)}{\hbar y} x^a \partial_x^{b-j} y^{\alpha+1} \partial_y^\beta \\
&\quad + \beta \frac{V'(x + \frac{\hbar}{2}y) + V'(x - \frac{\hbar}{2}y) - 2V'(x)}{2\hbar y} \hbar x^a \partial_x^b y^{\alpha+1} \partial_y^{\beta-1} \\
&\quad + \beta \frac{V'(x) - V'(0)}{x} x^{a+1} \partial_x^b y^\alpha \partial_y^{\beta-1} + \beta V'(0) x^a \partial_x^b y^\alpha \partial_y^{\beta-1} \\
&\quad + \sum_{k=2}^\beta \binom{\beta}{k} \frac{V^{(k)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \hbar^k x^a \partial_x^b y^\alpha \partial_y^{\beta-k} \\
&\quad + \sum_{k=1}^\beta \sum_{j=1}^b \binom{\beta}{k} \binom{b}{j} \frac{V^{(k+j)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k+j)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \hbar^k x^a \partial_x^{b-j} y^\alpha \partial_y^{\beta-k}.
\end{align*}
$$

The first term in the right-hand side involves the ratio

$$
\left| \frac{V'(x + \frac{\hbar}{2}y) - V'(x - \frac{\hbar}{2}y)}{\hbar y} \right| \le \|V''\|_{L^{\infty}} ,
$$

and the monomial $x^a \partial_x^{b-1} y^{\alpha+1} \partial_y^\beta$, whose exponents add up to

$a+b-1+\alpha+1+\beta=n,$

which is equal to the sum of exponents in the monomial $x^a \partial_x^b y^\alpha \partial_y^\beta$. The second term in the right-
hand side involves monomials whose exponents add up to

$a+b-j+\alpha+1+\beta<n$

since $j \ge 2$, multiplied by the ratio

$$
\left| \frac{V^{(j)}(x + \frac{\hbar}{2}y) - V^{(j)}(x - \frac{\hbar}{2}y)}{\hbar y} \right| \le \|V^{(j+1)}\|_{L^{\infty}} < \infty.
$$

The last three terms on the right-hand side involve monomials whose exponents add up to

$a+b-j+\alpha+\beta-k<n$
with $j+k \ge 1$, multiplied by $V'(0)$ or a ratio of the form

$$ \left| \frac{V^{(k+j)}(x + \frac{\hbar}{2}y) - (-1)^k V^{(k+j)}(x - \frac{\hbar}{2}y)}{2^k \hbar} \right| \hbar^k \le \|V^{(k+j)}\|_{L^\infty} \hbar^{k-1} $$

with $k \ge 1$ and $k+j \ge 2$. There remains the third and fourth terms, involving monomials whose exponents add up to $a+b+\alpha+\beta=n$, multiplied by ratios of the form

$$ \left| \frac{V'(x + \frac{\hbar}{2}y) + V'(x - \frac{\hbar}{2}y) - 2V'(x)}{2\hbar y} \right| \le \frac{1}{2} \|V''\|_{L^\infty} < \infty, $$

or

$$ \left| \frac{V'(x) - V'(0)}{x} \right| \le \|V''\|_{L^\infty} < \infty. $$

Summarizing, assuming that $\hbar \le 1$ while $a+b+\alpha+\beta=n$ with $a,b,\alpha,\beta \ge 0$,

$$ \begin{aligned} \|\mathcal{T}_2(s|\cdot, \cdot)\|_{L^2} &\le n \|V''\|_{L^\infty} \mathcal{N}_n(s) + \underbrace{\sum_{j=2}^b \binom{b}{j} \|V^{(j+1)}\|_{L^\infty} \mathcal{N}_{n-j+1}(s) + \frac{1}{2} \|V''\|_{L^\infty} n \hbar \mathcal{N}_n(s)}_{\le 2^n \max_{2 \le i \le n} \|V^{(i+1)}\|_{L^\infty} \mathcal{N}_{n-1}(s)} \\ &\quad + \|V''\|_{L^\infty} n \mathcal{N}_n(s) + n |V'(0)| \mathcal{N}_{n-1}(s) + \underbrace{\sum_{k=2}^\beta \binom{\beta}{k} \|V^{(k)}\|_{L^\infty} \hbar^{k-1} \mathcal{N}_{n-k}(s)}_{\le 2^n \max_{2 \le i \le n} \|V^{(i)}\|_{L^\infty} \mathcal{N}_{n-2}(s)} \\ &\quad + \underbrace{\sum_{k=1}^\beta \sum_{j=1}^b \binom{\beta}{k} \binom{b}{j} \|V^{(k+j)}\|_{L^\infty} \hbar^{k-1} \mathcal{N}_{n-k-j}(s)}_{\le 2^n \cdot 2^n \max_{2 \le i \le n} \|V^{(i)}\|_{L^\infty} \mathcal{N}_{n-2}(s)} \end{aligned} $$

Hence, it yields

$$ \begin{aligned} \|\mathcal{T}_2(s|\cdot, \cdot)\|_{L^2} &\le \max \left( |V'(0)|, \max_{2 \le k \le n+1} \|V^{(k)}\|_{L^\infty} \right) \left( n + 2^n + \frac{1}{2}n\hbar + n + n + 2^n\hbar + 2^n \cdot 2^n \right) \mathcal{N}_n(s) \\ &\le \max \left( |V'(0)|, \max_{2 \le k \le n+1} \|V^{(k)}\|_{L^\infty} \right) \left( \frac{7}{2} + 2^n \right) (2^n + 2) \mathcal{N}_n(s) \end{aligned} $$

since $0 < \hbar \le 1$.

Thus, we have proved that

$$ \begin{aligned} \|x^a \partial_x^b y^\alpha \partial_y^\beta R^\hbar(t)\|_{L^2} &\le \|x^a \partial_x^b y^\alpha \partial_y^\beta R^\hbar(0)\|_{L^2} \\ &\quad + \left( 1 + \max \left( |V'(0)|, \max_{2 \le k \le n+1} \|V^{(k)}\|_{L^\infty} \right) \right) \int_0^t \left( \frac{7}{2} + 2^n \right) (2^n + 2) \mathcal{N}_n(s) \,ds, \end{aligned} $$

so that

$$ \begin{aligned} \mathcal{N}_n(t) &\le \mathcal{N}_n(0) \\ &\quad + \left( 1 + \max \left( |V'(0)|, \max_{2 \le k \le n+1} \|V^{(k)}\|_{L^\infty} \right) \right) \left( \frac{7}{2} + 2^n \right) (2^n + 2)(n+1)^4 \int_0^t \mathcal{N}_n(s) \,ds \end{aligned} $$

and one concludes by Gronwall's inequality.
REFERENCES

[1] A. G. Athanassoulis. Smoothed Wigner transforms in the numerical simulation of semiclassical (high-frequency) wave propagation. Preprint, arXiv:0704.3404, 2007.

[2] P. Bader, A. Iserles, K. Kropielnicka, and P. Singh. Effective approximation for the semiclassical Schrödinger equation. *Foundations of Computational Mathematics*, 14:689–720, 2014.

[3] W. Bao, S. Jin, and P. A. Markowich. On time-splitting spectral approximations for the Schrödinger equation in the semiclassical regime. *Journal of Computational Physics*, 175(2):487–524, 2002.

[4] C. Bardos, F. Golse, and D. Levermore. Fluid dynamic limits of kinetic equations. I. Formal derivations. *Journal of Statistical Physics*, 63(1-2):323–344, 1991.

[5] M. Berman, R. Kosloff, and H. Tal-Ezer. Solution of the time-dependent Liouville-von Neumann equation: dissipative evolution. *Journal of Physics A: Mathematical and General*, 25(5):1283, 1992.

[6] M. Bessemoulin-Chatard and F. Filbet. On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. *Journal of Computational Physics*, 451:Paper No. 110881, 28, 2022.

[7] M. Bessemoulin-Chatard and F. Filbet. On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. *SIAM Journal on Numerical Analysis*, 61(4):1664–1688, 2023.

[8] A. Blaustein. Structure preserving solver for Multi-dimensional Vlasov-Poisson type equations. Preprint https://hal.science/hal-04440391v1, 2024.

[9] A. Blaustein and F. Filbet. A structure and asymptotic preserving scheme for the Vlasov-Poisson-Fokker-Planck model. *Journal of Computational Physics*, 498:112693, 2024.

[10] A. Blaustein and F. Filbet. On a discrete framework of hypocoercivity for kinetic equations. *Mathematics of Computation*, 93(345):163–202, 2024.

[11] Y. Brenier. Convergence of the Vlasov-Poisson system to the incompressible Euler equations. *Communications in Partial Differential Equations*, 25(3-4):737–754, 2000.

[12] H. Brézis. *Functional Analysis, Sobolev Spaces and Partial Differential Equations*. Springer Science+Business Media, LLC, 2011.

[13] D. Fang, S. Jin, and C. Sparber. An efficient time-splitting method for the Ehrenfest dynamics. *Multiscale Modeling and Simulation*, 16(2):900–921, 2018.

[14] F. Filbet and F. Golse. On the approximation of the von-Neumann equation in the semi-classical limit. Part I: Numerical algorithm. *Journal of Computational Physics*, 527:113810, 2025.

[15] F. Filbet and T. Xiong. Conservative discontinuous Galerkin/Hermite spectral method for the Vlasov-Poisson system. *Communications on Applied Mathematics and Computation*, 4(1):34–59, 2022.

[16] F. Golse, S. Jin, and T. Paul. On the convergence of time splitting methods for quantum dynamics in the semiclassical regime. *Foundations of Computational Mathematics*, 21(3):613–647, 2021.

[17] F. Golse and J. Möller. Velocity Averaging for the Mixed State Wigner Transform. preprint, 2024.

[18] B. Helling and H. Metiu. An efficient method for solving the quantum Liouville equation: Applications to electronic absorption spectroscopy. *Chemical Physics Letters*, 127(1):45–49, 1986.

[19] M. Hochbruck and C. Lubich. Exponential integrators for quantum-classical molecular dynamics. *BIT Numerical Mathematics*, 39:620–645, 1999.

[20] H. Hochstadt. *The functions of mathematical physics*. Courier Corporation, 2012.

[21] J. P. Holloway. Spectral velocity discretizations for the Vlasov-Maxwell equations. *Transport Theory and Statistical Physics*, 25(1):1–32, 1996.

[22] S. Jin. Asymptotic-preserving schemes for multiscale physical problems. *Acta Numerica*, 31:415–489, 2022.

[23] S. Jin, P. Markowich, and C. Sparber. Mathematical and computational methods for semiclassical Schrödinger equations. *Acta Numerica*, 20:121–209, 2011.

[24] L. D. Landau and E. M. Lifshitz. *Quantum mechanics: non-relativistic theory*, volume 3. Elsevier, 2013.

[25] P. D. Lax. The scope of the energy method. *Bulletin of the American Mathematical Society*, 66:32–35, 1960.

[26] P. D. Lax and R. D. Richtmyer. Survey of the stability of linear finite difference equations. *Communications on Pure and Applied Mathematics*, 9:267–293, 1956.

[27] N. Lerner. *Metrics on the phase space and non-self adjoint pseudo-differential operators*, volume 3. Springer Science & Business Media, 2011.

[28] P.-L. Lions and T. Paul. Sur les mesures de Wigner. *Revista Matemática Iberoamericana*, 9(3):553–618, 1993.

[29] P. A. Markowich and H. Neunzert. On the equivalence of the Schrödinger and the quantum Liouville equations. *Mathematical Methods in the Applied Sciences*, 11(4):459–469, 1989.

[30] P. A. Markowich, P. Pietra, and C. Pohl. Numerical approximation of quadratic observables of Schrödinger-type equations in the semi-classical limit. *Numerische Mathematik*, 81:595–630, 1999.

[31] P. A. Markowich and C. Ringhofer. An analysis of the quantum Liouville equation. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik*, 69(3):121–127, 1989.
[32] P. A. Markowich, C. A. Ringhofer, and C. Schmeiser. *Semiconductor equations*. Springer Science & Business Media, 2012.

[33] B. Miao, G. Russo, and Z. Zhou. A novel spectral method for the semiclassical Schrödinger equation based on the Gaussian wave-packet transform. *IMA Journal of Numerical Analysis*, 43(2):1221–1261, 2023.

[34] G. Russo and P. Smereka. The Gaussian wave packet transform: Efficient computation of the semi-classical limit of the Schrödinger equation. Part 1—Formulation and the one dimensional case. *Journal of Computational Physics*, 233:192–209, 2013.

[35] G. Russo and P. Smereka. The Gaussian wave packet transform: efficient computation of the semi-classical limit of the Schrödinger equation. Part 2. Multidimensional case. *Journal of Computational Physics*, 257:1022–1038, 2014.

[36] J. W. Schumer and J. P. Holloway. Vlasov simulations using velocity-scaled Hermite representations. *Journal of Computational Physics*, 144(2):626–661, 1998.

[37] J. Sellier, M. Nedjalkov, I. Dimov, and S. Selberherr. A comparison of approaches for the solution of the Wigner equation. *Mathematics and Computers in Simulation*, 107:108–119, 2015.

[38] D. Sels, F. Brosens, and W. Magnus. Wigner distribution functions for complex dynamical systems: A path integral approach. *Physica A: Statistical Mechanics and its Applications*, 392(2):326–335, 2013.

[39] S. Shao, T. Lu, and W. Cai. Adaptive conservative cell average spectral element methods for transient Wigner equation in quantum transport. *Communications in Computational Physics*, 9(3):711–739, 2011.

[40] M. L. Van de Put, B. Sorée, and W. Magnus. Efficient solution of the Wigner-Liouville equation using a spectral decomposition of the force field. *Journal of Computational Physics*, 350:314–325, 2017.

[41] Y. Xie and Z. Zhou. Frozen Gaussian sampling: a mesh-free Monte Carlo method for approximating semiclassical Schrödinger equations. *Communications in Mathematical Sciences*, 22(4):1133–1166, 2024.

[42] Y. Xiong, Z. Chen, and S. Shao. An advective-spectral-mixed method for time-dependent many-body Wigner simulations. *SIAM Journal on Scientific Computing*, 38(4):B491–B520, 2016.

(F.F.) IMT, UNIVERSITÉ DE TOULOUSE, 31062 TOULOUSE CEDEX, FRANCE
*Email address: francis.filbet@math.univ-toulouse.fr*

(F.G.) CMLS, ÉCOLE POLYTECHNIQUE, 91128 PALAISEAU CEDEX, FRANCE
*Email address: francois.golse@polytechnique.edu*