\documentclass[JEP,XML,SOM,Unicode,NoEqCountersInSection,published]{cedram}
\datereceived{2022-04-08}
\dateaccepted{2023-03-20}
\dateepreuves{2023-04-03}

\makeatletter
\def\@settitle{%
\vspace*{-8mm}
\raggedleft\includegraphics[scale=.5]{titre-jep}
\vtop to 50 mm{%
 \parindent=0pt
 {\abstractfont\article@logo\par}
 \medskip
 \hrule
 \vfil
 \begin{center}
 \def\baselinestretch{1.2}\large\vfil
   {\didottitraille\MakeUppercase\@title\par}
 \vfil\vfil
 \begin{minipage}{.8\textwidth}\centering
   \ifx\@empty\smfbyname\else
   {\smf@byfont\smfbyname\ifsmf@byauthor\enspace\else\ \fi}%
   \fi {\smf@authorfont \edef\smfandname{{\noexpand\smf@andfont
         \smfandname}} \andify\authors\authors\par}
 \end{minipage}
 \vfil \vrule height .4pt width .3\textwidth \vfil
 \end{center}}%
 \par\enlargethispage{.5\baselineskip}%
}
\def\@setthanks{\def\thanks##1{\par##1\@addpunct{{\upshape.}}}\vspace*{-5pt}\thankses}
\makeatother

\usepackage[scr=rsfs,cal=euler]{mathalfa}
\multlinegap0pt
\newcommand{\spfrac}[2]{\sfrac{#1}{(#2)}}
\newcommand{\psfrac}[2]{\sfrac{(#1)}{#2}}
\newcommand{\pspfrac}[2]{\sfrac{(#1)}{(#2)}}
\newcommand{\Psfrac}[2]{(\sfrac{#1}{#2})}
\newcommand {\aTr}[1] {\tr(#1)}
\newcommand {\bTr}[1] {\tr\bigl( #1 \bigr)}
\newcommand\mto{\mathchoice{\longmapsto}{\mapsto}{\mapsto}{\mapsto}}
\newcommand {\acom}[1] {[{#1}]}
\newcommand {\aNrm}[2] {\lVert{#1}\rVert_{#2}}
\newcommand {\an}[1] {\lvert #1\rvert}
\let\epsilon\varepsilon

\usepackage{bbm}
\let\mathds\mathbbm
\usepackage{braket}

\newtheorem{lem}{Lemma}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{cor}{Corollary}[section]
\newtheorem{prop}{Proposition}[section]

\theoremstyle{definition}
\newtheorem{definition}{Definition}[section]
\newtheorem{hyp}{Hypothesis}
\newtheorem{remark}{Remark}[section]

\newcommand{\step}[1] {\subsubsection*{#1}}

\newcommand {\N} {\mathbb N}
\newcommand {\RR} {\mathbb R}
\let\R\RR
\newcommand {\Rd} {\RR^3}
\newcommand {\RRd} {\RR^6}
\let\Rdd \RRd

\renewcommand {\L} {\mathcal L}
\newcommand {\cW} {\mathcal W}
\newcommand {\cL} {\mathcal L}

\newcommand \sfm {\mathsf m}

\newcommand {\ecF} {\mathscr{F}}

\let\lt\left
\let\rt\right
\renewcommand {\(} {\lt(}
\renewcommand {\)} {\rt)}
\newcommand {\bangle}[1] {\lt\langle #1\rt\rangle}
\newcommand {\weight}[1] {\bangle{#1}}
\newcommand {\com}[1] {\lt[{#1}\rt]}
\newcommand {\bcom}[1] {\big[{#1}\big]}
\newcommand {\Bcom}[1] {\Big[{#1}\Big]}
\newcommand {\n}[1] {\lt\lvert #1 \rt\rvert}
\newcommand {\norm}[1] {\lVert #1 \rVert}
\newcommand {\nrm}[1] {\lt\lVert #1\rt\rVert}
\newcommand {\bnrm}[1] {\big\lVert #1\big\rVert}
\newcommand {\Bnrm}[1] {\Big\lVert #1\Big\rVert}
\newcommand {\Nrm}[2] {\nrm{#1}_{#2}}
\newcommand {\bNrm}[2] {\bnrm{#1}_{#2}}
\newcommand {\BNrm}[2] {\Bnrm{#1}_{#2}}
\newcommand {\llp}[1] {\norm{#1}}
\newcommand {\Lp}[2] {\Nrm{#1}{#2}}

\let\bd\partial
\let\eps\varepsilon
\let\grad\nabla
\let\lapl\Delta
\renewcommand {\d} {\mathop{}\!\mathrm{d}}
\newcommand {\dpt} {\partial_t}
\newcommand {\dt} {\frac{\d}{\d t}\!}
\newcommand {\Dx} {\nabla_x}
\newcommand {\Dv} {\nabla_\xi}
\newcommand {\conj}[1] {\overline{#1}}
\newcommand {\Id} {\mathrm{Id}}

\DeclareMathOperator{\tr} {Tr}
\DeclareMathOperator{\diag} {diag}

\newcommand {\F}[1] {\ecF\!\( #1 \)}
\newcommand {\Tr}[1] {\tr\!\( #1 \)}
\newcommand {\Diag}[1] {\diag\!\( #1 \)}

\newcommand {\intd} {\int_{\Rd}}
\newcommand {\intdd} {\int_{\RRd}}
\newcommand {\iintd} {\iint_{\RRd}}
\newcommand {\iintdd} {\iint_{\RRd\times\RRd}}

\newcommand {\ii} {\mathrm{i}}
\newcommand {\init} {\mathrm{in}}

\newcommand {\cC} {\mathcal{C}}

\newcommand {\Inprod}[2] {\Braket{#1 | #2}}
\renewcommand {\r} {\op}
\newcommand {\op} {\boldsymbol{\rho}}
\newcommand{\bp}{\boldsymbol{p}}
\let\opm\sfm
\newcommand {\opmu} {\boldsymbol{\mu}}
\newcommand {\opnu} {\boldsymbol{\nu}}
\newcommand {\opp} {\boldsymbol{p}}
\newcommand {\Dh} {\boldsymbol{\nabla}}
\newcommand {\Dhx}[1] {\Dh_{\!x} #1}
\newcommand {\Dhv}[1] {\Dh_{\!\xi} #1}

\datepublished{2023-04-07}
\begin{document}
\frontmatter
\title[On the $L^2$ rate of convergence]{On the $L^2$ rate of convergence in the limit from the Hartree to the Vlasov--Poisson equation}

\author[\initial{J.} \lastname{Chong}]{\firstname{Jacky} \middlename {J.} \lastname{Chong}}
\address{Department of Mathematics, The University of Texas at Austin\\ Austin, TX 78712, USA\\
\& Department of Mathematics, School of Mathematical Sciences, Peking University\\
Beijing, China}
\email{jwchong@math.utexas.edu}
\email{jwchong@math.pku.edu.cn}
\urladdr{https://web.ma.utexas.edu/users/jwchong/}

\author[\initial{L.} \lastname{Lafleche}]{\firstname{Laurent} \lastname{Lafleche}}
\address{Department of Mathematics, The University of Texas at Austin\\ Austin, TX 78712, USA\\
\& Institut Camille Jordan, CNRS, Université Claude Bernard Lyon 1\\
43 boulevard du 11 novembre 1918, F-69622 Villeurbanne Cedex, France}
\email{lafleche@math.univ-lyon1.fr}
\urladdr{https://laurent-lafleche.perso.math.cnrs.fr/}

\author[\initial{C.} \lastname{Saffirio}]{\firstname{Chiara} \lastname{Saffirio}}
\address{Department of Mathematics and Computer Science, University of Basel\\
4051 Basel, Switzerland}
\email{chiara.saffirio@unibas.ch}
\urladdr{}

\thanks{L.L. has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 865711). J.C. was supported by the NSF through the RTG grant DMS-RTG 184031. C.S. acknowledges the NCCR SwissMAP and the support of the SNSF through the Eccellenza project PCEFP2\_181153}

\subjclass{81Q20, 35Q55, 35Q83, 82C10, 82C05}

\keywords{Semiclassical limit, Hartree equation, Vlasov equation, Coulomb potential, gravitational potential.}

\begin{abstract}
Using a new stability estimate for the difference of the square roots of two solutions of the Vlasov--Poisson equation, we obtain the convergence in the $L^2$ norm of the Wigner transform of a solution of the Hartree equation with Coulomb potential to a solution of the Vlasov--Poisson equation, with a rate of convergence proportional to $\hbar$. This improves the $\hbar^{3/4-\varepsilon}$ rate of convergence in $L^2$ obtained in [L.\,Lafleche, C.\,Saffirio: \textit{Analysis \& PDE}, to appear]. Another reason of interest of this paper is the new method, reminiscent of the ones used to prove the mean-field limit from the many-body Schrödinger equation towards the Hartree--Fock equation for mixed states.
\end{abstract}

\altkeywords{Limite semi-classique, équation de Hartree, équation de Vlasov, potentiel de Coulomb, potentiel gravitationnel}

\alttitle{Sur le taux de convergence dans $L^2$ dans la limite de l'équation de Hartree à l'équation de Vlasov-Poisson}

\begin{altabstract}
Grâce à une nouvelle estimée de stabilité pour la différence entre les racines carrées de deux solutions de l'équation de Vlasov-Poisson, nous obtenons la convergence en norme $L^2$ de la transformée de Wigner d'une solution de l'équation de Hartree avec potentiel de Coulomb vers une solution de l'équation de Vlasov-Poisson, avec un taux de convergence proportionnel à la constante de Planck $\hbar$. Ceci améliore le taux de convergence $h^{3/4 - \varepsilon}$ dans $L^2$ obtenu dans [L.\,Lafleche, C.\,Saffirio: \textit{Analysis \& PDE}, à paraître]. Un autre intérêt de cet article est la nouvelle méthode, réminiscente de celles utilisées pour prouver la limite de champ moyen de l'équation de Schrödinger à $N$ corps vers l'équation de Hartree-Fock pour des états mixtes.
\end{altabstract}
\maketitle
\vspace*{-1.5\baselineskip}
\tableofcontents
\mainmatter

\vspace*{-4\baselineskip}\mbox{}
\section{Introduction and main result}

In this paper, we study the limit of solutions to the Hartree equation towards solutions to the Vlasov--Poisson equation in the semiclassical regime, that is when the Planck constant converges to $0$. We focus on the cases of the three dimensional Coulomb and gravitational interaction potentials
\begin{equation}\label{eq:K}
K(x) = \frac{\pm 1}{4\pi\n{x}},
\end{equation}
which are the most relevant singular potentials from a physical viewpoint and the most challenging from a mathematical viewpoint.

We consider the Hartree equation
\begin{equation}\label{eq:Hartree}
i\hbar\,\dpt \op = \com{H_{\op},\op} \quad \text{ with } \quad \op(0)= \op^{\init},
\end{equation}
where $\hbar=\sfrac{h}{2\pi}$ is the reduced Planck constant, $\op = \op(t)$ is a time-dependent positive self-adjoint trace class operator acting on $L^2(\Rd)$ satisfying the normalization
\begin{equation}\label{def:normalization}
h^3\Tr{\op}=1\quad \text{ and } \quad \Nrm{\op}{\infty}=\cC_\infty,
\end{equation}
where $\Nrm{\cdot}{\infty}$ is the operator norm and $\cC_\infty$ is a fixed independent-of-$\hbar$ constant.
Here $H_{\op} = -\sfrac{\hbar^2\Delta}{2} + V_{\op}$ denotes the Hartree Hamiltonian where the operator $V_{\op}$ is the operator of multiplication by the mean-field potential $V_{\op}(x) = (K*\rho)(x)$ associated to the two-body interaction potential $K$ defined in our case by Equation~\eqref{eq:K}, and $\rho(x)$ is the quantum spatial density defined as the scaled restriction to the diagonal of the integral kernel of $\op$,
\begin{equation*}
\rho(x)=\diag(\op)(x) := h^3\op(x,x),
\end{equation*}
where we adopt the same notation to denote both the operator $\op$ and its integral kernel $\op(x,y)$. The choice of the normalization \eqref{def:normalization} is motivated by semiclassical analysis and is natural in the analysis of many-body fermionic systems, although the results of this paper apply to both the bosonic and fermionic settings.

\noindent The classical analogue of Equation~\eqref{eq:Hartree} is the Vlasov--Poisson equation given by
\begin{equation}\label{eq:Vlasov}
\dpt f + \xi\cdot\Dx f + E_f \cdot\Dv f = 0 \quad \text{ with } \quad f(0, x, \xi)= f^{\init}(x, \xi)\ge 0,
\end{equation}
where $E_f = -\nabla V_f$ is the self-consistent force field associated to the mean-field poten\-tial $V_f(x)=(K*\rho_f)({t,}x)$ with $\rho_f$ the spatial density given by
\begin{equation*}
\rho_f({t,}x)=\intd f(t,x,\xi) \d\xi.
\end{equation*}

For the well-posedness of the Hartree and the Vlasov--Poisson equations we refer to \cite{castella_$l^2$_1997} and \cite{pfaffelmoser_global_1992, lions_propagation_1991} respectively and the references therein.

Equation~\eqref{eq:Vlasov} can be seen as the semiclassical approximation of a system of many interacting quantum particles, as pointed out in the pioneering works by Narnhofer and Sewell~\cite{narnhofer_vlasov_1981} and by Spohn~\cite{spohn_vlasov_1981} where the Vlasov equation was obtained directly from the many-body Schrödinger equation with smooth interaction in the combined mean-field and semiclassical regime. This has been reconsidered in~\cite{graffi_mean-field_2003} and more recently in~\cite{chong_many-body_2021, chen_combined_2021}, where the case of the Coulomb potential with a $N$ dependent cut-off has been addressed. Moreover, a combined mean-field and semiclassical limit for particles interacting via the Coulomb potential has been treated in~\cite{golse_mean-field_2021} for factorized initial data whose first marginal is given by a monokinetic Wigner measure (that can be seen as the Klimontovich solutions to the Vlasov equation), which leads to the pressureless Euler--Poisson system.

Most of the above mentioned works rely on compactness methods that do not allow for an explicit bound on the rate of convergence, which is essential for applications. For this reason the Hartree equation~\eqref{eq:Hartree} has been considered as an intermediate step to decouple the problem into two separate parts, namely to prove the convergence of the mean-field limit from the many-body Schrödinger equation towards the Hartree equation, and then the semiclassical limit from the Hartree equation to the Vlasov equation. In this paper, we are interested in the latter problem, that has been largely studied in different settings. It was first proved by Lions and Paul in~\cite{lions_sur_1993}, and later in~\cite{markowich_classical_1993,figalli_semiclassical_2012}, that the Wigner transforms of the solutions of the Hartree equation~\eqref{eq:Hartree} converge in some weak sense to solutions of the Vlasov--Poisson equation. Quantitative rates of convergence were then obtained, first in the case when the Coulomb potential is replaced by a smoother potential, in Lebesgue-type norms~\cite{athanassoulis_strong_2011, amour_semiclassical_2013, benedikter_hartree_2016} and in a quantum analogue of the Wasserstein distances \cite{golse_schrodinger_2017}. The case of singular interactions was then treated in~\cite{lafleche_propagation_2019, lafleche_global_2021} with the same quantum Wasserstein distances, and in~\cite{saffirio_semiclassical_2019, saffirio_hartree_2020, lafleche_strong_2021} in Lebesgue-type norms. In particular, for $K$ as in Equation~\eqref{eq:K}, the explicit rate has been established in~\cite{lafleche_propagation_2019} for the weak topology and in~\cite{saffirio_hartree_2020, lafleche_strong_2021} for the Schatten norms.

In a different setting, the semiclassical limit has also been studied for local perturbations of stationary states in the case of infinite gases in~\cite{lewin_hartree_2020}.

Our goals here are twofold.~On the one hand, we want to improve the rate of convergence in the $L^2$ norm in terms of $\hbar$ for the convergence of the Wigner transform of the solution of the Hartree equation to the Vlasov--Poisson equation. The rate $\hbar^{3/4-\eps}$, for $\eps>0$ small enough, was obtained in~\cite{lafleche_strong_2021}, the reason being that the method used in that previous paper was relying on a $L^1$ weak-strong uniqueness principle, thus leading to the expected correct rate of order $\hbar$ in trace norm, but lower order rates for higher Schatten norms. In this paper we recover a rate of order~$\hbar$ for the $L^2$ convergence, that is expected to be optimal.

On the other hand, the ideas used in this paper are closer to those used in~\cite{benedikter_mean-field_2016-1, chong_many-body_2021} to prove the semiclassical mean-field limit from the many-body Schrödinger equation towards the Hartree--Fock equation for mixed states, in the sense that it considers the square roots of the phase space densities in the $L^2$ setting. Inspired by these ideas, we introduce a new $L^2$ weak-strong stability estimate for the Vlasov--Poisson and Hartree equations, which can be compared with the $L^1$ method used in~\cite{lafleche_strong_2021}. It serves as a guiding principle to our main theorem.

\Subsection{Notations}

\subsubsection{Functional spaces} Let us fix some notations before stating our main result. A point in the one-particle phase space $\Rd_x\times \Rd_{\xi}$ is denoted by $z = (x, \xi)$. We denote the phase space Lebesgue norm and its mixed norm variant by
\begin{equation*}
\Nrm{f}{L^r} = \Big(\intdd \n{f(z)}^r\d z\Big)^{\!\sfrac{1}{r}}, \quad \Nrm{f}{L^p_x L^q_{\xi}} = \bNrm{\!\Nrm{f(x, \cdot)}{L^q_{\xi}(\Rd)}\!}{L^p_x(\Rd)}
\end{equation*}
for $1\le r<\infty$ and with the usual modification when $r=\infty$. The corresponding weighted Sobolev spaces $W^{k,p}_n$ are defined by the norm
\begin{equation*}
\Nrm{f}{W^{k,p}_{n}} = \Big(\sum_{\substack{ \n{\alpha} \le k}}\Nrm{\weight{\xi}^n\partial^\alpha f}{L^{p}}^2\Big)^{\!\sfrac{1}{2}}\quad \text{ with }\quad
\begin{cases}
\alpha \in {\N_{0}^6},\, \n{\alpha}= \sum^6_{\ii=1} \alpha_{\ii},
\\
\bd^\alpha = \bd_{x_1}^{\alpha_1}\cdots \bd_{\xi_3}^{\alpha_6},
\end{cases}
\end{equation*}
where $\weight{z}:=\sqrt{1+|z|^2}$ is the standard bracket notation. We also use the standard notation $H^{k}_{n} = W^{k,2}_{n}$.

Due to our choice of the potential, it is natural to introduce the Lorentz spaces $L^{p,q}_x$ with $(p,q)\in[1,\infty]^2$ that are intermediate spaces between Lebesgue spaces such that $\nabla K \in L^{3/2,\infty}$. They can indeed be defined by real interpolation of Lebesgue spaces (see \eg \cite[Th.\,5.3.1]{bergh_interpolation_1976}), and as such they verify for any $\eps\in(0,p-1]$, $L^{p-\eps} \cap L^{p+\eps} \subset L^{p,q}$. Moreover, their dual is given by $(L^{p,q})' = L^{p',q'}$ where $p' = \spfrac{p}{p-1}$ denotes the Hölder conjugate (see \eg \cite{hunt_lpq_1966}). This implies the analogue of Hölder's inequality for these spaces.

For any bounded linear operator $\op$ acting on $L^2(\Rd)$, we denote its operator norm by $\Nrm{\op}{\infty}$, and for $1\le p < \infty$, we define its Schatten-$p$ norm by
\begin{equation*}
\Nrm{\op}{p}= \aTr{\n{\op}^p}^{\sfrac1p},
\end{equation*}
where $\n{\op} = \sqrt{\op^\ast \op}$. To obtain meaningful quantities in the semiclassical limit \hbox{$\hbar\to 0$}, we also define the semiclassical analogue of the Schatten spaces via the rescaled Schatten norm
\begin{equation*}
\Nrm{\op}{\L^p} = h^{\sfrac3p} \Nrm{\op}{p}= h^{\sfrac3p}\Tr{\n{\op}^p}^{\sfrac1p}
\end{equation*}
and $\Nrm{\op}{\cL^\infty}=\Nrm{\op}{\infty}$. Notice that, unlike its unscaled counterparts, the semiclassical Schatten norms do not enjoy the same inclusion inequality in the limit as $\hbar\to 0$, that is, $\Nrm{\op}{\cL^q}$ is not bounded uniformly in $\hbar$ by $\Nrm{\op}{\cL^p}$ for any $1\le p< q\le \infty$. These semiclassical spaces play the roles of the Lebesgue spaces on the phase space. We will denote by $C$ a $\hbar$-independent constant that can change from line to line.

\subsubsection{Weyl quantization and Wigner transform} The Fourier transform is defined with the convention that
\begin{equation*}
\widehat g(\xi) = {\F{g}(\xi) =} \intd e^{-2\pi i\, x\cdot \xi}\,g(x)\d x.
\end{equation*}
Consistently with our normalization~\eqref{def:normalization}, we associate to every function \hbox{$f \!\in\! L^1(\Rd_x\!\times\!\Rd_{\xi})$} the operator $\op_f$ with integral kernel
\begin{equation*}
\op_f(x, y) = \intd e^{-2\pi i\, (y-x)\cdot \xi} f(\psfrac{x+y}{2}, h\xi)\d \xi.
\end{equation*}
The operator $\op_f$ is called the Weyl quantization associated to $f(z)$. By the Fourier inversion theorem, we could also define the inverse mapping called the Wigner transform. More precisely, for any operator $\op$ with sufficiently regular kernel $\op(x, y)$, we~define its Wigner transform by
\begin{equation*}
f_{\op}(x, \xi) = \intd e^{-2\pi i\, \xi\cdot y/h} \op(x+\sfrac{y}{2}, x-\sfrac{y}{2})\d y.
\end{equation*}

\subsection{Main result}
Our main result reads as follows.
\begin{theorem}\label{thm:main}
Let $\op\geq 0$ be a solution to the Hartree equation~\eqref{eq:Hartree} with initial condition $\op^\init\in\L^1\cap\L^\infty$ verifying~\eqref{def:normalization} and $f\geq 0$, $f\in L^1(\R^6)$ be a solution to the Vlasov--Poisson equation~\eqref{eq:Vlasov} with initial condition $f^\init$ verifying
\begin{equation*}
f^\init,\,\sqrt{f^\init}\in W^{4,\infty}_4 \cap H^4_4 \quad \text{ and } \quad \iintd f^\init \n{\xi}^{n_1}\d x\d\xi <\infty
\end{equation*}
for $n_1 > 6$. Then there exist time-dependent functions $\Lambda, C_1, C_2 \in C^0(\R_+,\R_+)$, inde\-pendent of $\hbar$ and depending on the initial conditions of equations~\eqref{eq:Hartree} and~\eqref{eq:Vlasov} such~that
\begin{equation}\label{eq:Wigner_main_estimate}
\Nrm{f_{\op}-f}{L^2} \leq \cC_{\infty}^{\sfrac12} \(\bNrm{f_{\sqrt{\op^\init}} - \sqrt{f^\init}}{L^2} + C_1(t) \,\hbar\) e^{\Lambda(t)} + C_2(t)\,\hbar.
\end{equation}
\end{theorem}

\begin{remark}
The behavior of the time-dependent functions $\Lambda$, $C_1$ and $C_2$ appearing in the theorem depends strongly on the propagation of weighted Sobolev norms for the Vlasov--Poisson equation, and might lead in the worse case to functions growing faster in time than $e^{e^t}$. To simplify the exposition, more details about these functions are given in the next section.
\end{remark}

\begin{remark}
Our method also allows us to treat the Hartree--Fock equation.
In~that case, the exchange term vanishes in the semiclassical limit and can be treated as an error term as in~\cite[Prop.\,5.1]{lafleche_strong_2021}.
\end{remark}

\subsubsection{Strategy and explicit constants} Before giving more precise upper bounds on the functions $\Lambda$, $C_1$ and $C_2$, let us explain our strategy. The rationale of our result is a stability estimate for the square root of the solutions of the Vlasov equation in~$L^2$, as explained in Section~\ref{sec:stability-est}. Rephrasing this estimate in the quantum context to estimate the difference between the solutions to the Hartree equation and the Vlasov equation, we have to deal with the fact that the Weyl quantization of a non-negative function is not always non-negative. Thus, we have to consider an alternative quantization of $f$, sometimes called the Wick quantization, the anti-Wick quantization, the Töplitz quantization, the Husimi quantization or the coherent state quantization. We~introduce the coherent state at the point $z$ and its corresponding density operator
\begin{equation*}
\psi_z(y):= \(\pi \hbar\)^{-3/4} e^{-\n{y-x}^2/(2\hbar)}\,e^{i y\cdot \xi/\hbar} \quad \text{ and } \quad \op_z := h^{-3}\ket{\psi_z}\bra{\psi_z}.
\end{equation*}
Then, to every function $f$ on the phase space, we associate the operator $\widetilde{\op}_f$ defined by taking averages of the coherent states $\op_z$ against $f(z)$, that is
\begin{equation}\label{eq:wick}
\widetilde{\op}_f = \intdd f(z)\,\op_z\d z = \frac{1}{h^3}\intdd f(z) \ket{\psi_z}\bra{\psi_z}\d z.
\end{equation}
It follows from this formula that $\widetilde{\op}_f$ is a positive operator whenever $f\ge 0$. We~summarize some of its other well known properties which we will need in this work in the following lemma. We refer the reader for example to \cite{lions_sur_1993, lerner_fefferman-phong_2007} and the references therein for additional properties of this quantization and the proof of the following lemma.

\begin{lem}\label{lem:properties_of_Wick}
Let $f$ be a function of the phase space. Then the quantization~\eqref{eq:wick} is obtained from the Weyl quantization by a convolution with a Gaussian
\begin{equation}\label{eq:Wick_convolution}
\widetilde{\op}_f =\op_{\tilde f}, \quad \text{ where } \quad \widetilde f = g_h\ast f \text{ with } g_h(z):= (\pi\hbar)^{-3}\,e^{-\n{z}^2/\hbar}.
\end{equation}
By \eqref{eq:Wick_convolution}, we have the estimate
\begin{equation}\label{eq:Wick_Weyl_error}
\Nrm{\op_f-\widetilde{\op}_f}{\L^2} = \Nrm{f-g_h\ast f}{L^2} \le \frac{3\, \hbar}{2}\,\aNrm{\nabla^2 f}{L^2}.
\end{equation}
For $1\le p\le \infty$, we have that
\begin{equation}\label{eq:Wick_Schatten}
\aNrm{\widetilde{\op}_f}{\L^p}\le \Nrm{f}{L^p}.
\end{equation}
\end{lem}

This quantization will allow us to consider an intermediate equation between the Hartree and the Vlasov--Poisson equations, depending on the mean-field force of the Vlasov equation, but whose solution is a positive operator. More precisely, we consider~$\widetilde{\op}$ a solution to the equation
\begin{equation*}
i\hbar\,\dpt\widetilde{\op} = \com{H_f,\widetilde{\op}}
\end{equation*}
with Hamiltonian $H_f=-\Psfrac{\hbar^2}{2}\,\Delta+V_{f}$ and initial data $\widetilde{\op}^\init := \widetilde{\op}_{\sqrt{f^\init}}^2$, and we define
\begin{equation*}
\tilde v := \sqrt{\widetilde{\op}},
\end{equation*}
which satisfies the same equation.

\subsubsection{Quantum Sobolev spaces} To define a semiclassical version of Sobolev spaces on the phase space, we introduce the following operators
\begin{equation*}
\Dhx \op := \com{\nabla,\op} \quad \text{ and } \quad
\Dhv \op := \Bcom{\frac{x}{i\hbar},\op},
\end{equation*}
which can be seen as an application of the correspondence principle of quantum mechanics. More precisely, one observes that these operators correspond to the gradients of the Wigner transform, that is
\begin{equation*}
f_{\Dhx \op} = \Dx f_{\op}\quad \text{ and } \quad f_{\Dhv \op} = \Dv f_{\op}.
\end{equation*}
We will refer to $\Dhx\op$ and $\Dhv\op$ as the quantum gradients. The semiclassical analogues of the kinetic homogeneous Sobolev norms are defined by
\begin{subequations}
\begin{align*}
\Nrm{\op}{\dot{\cW}^{k,p}}^p &:= \sum_{\n{\alpha}=k} \Nrm{\Dh_{\!z}^\alpha\op}{\L^p}^p\quad \text{ with } \quad \Dh_{\!z}^\alpha = \Dh_{\!x_1}^{\alpha_1}\cdots\Dh_{\!\xi_3}^{\alpha_6},
\\
\Nrm{\op}{\dot{\cW}^{k,\infty}} &:= \sup_{\n{\alpha}= k} \aNrm{\Dh_{\!z}^\alpha\op}{\L^\infty},
\end{align*}
\end{subequations}
and the inhomogeneous version by
\begin{align*}
\Nrm{\op}{\cW^{k,p}}^p &:= \Nrm{\op}{\L^p}^p + \sum_{1\le \ell \le k}{\Nrm{\op}{\dot{\cW}^{\ell,p}}^p},
\end{align*}
for any $k\in\N$ with the usual modification when $p=\infty$. Notice that for $p=2$, we have $\Nrm{\op}{\cW^{k,2}} = \Nrm{f_{\op}}{H^k}$ which is suggestive of the strong connection between $\cW^{k, p}$ and the classical kinetic Sobolev spaces. Moreover, we will also consider weighted versions of these norms with the particular case of the weights defined for $n\in \R$ by
\begin{equation}\label{eq:weight}
\opm := \weight{\opp}^n = \big(1 + \n{\opp}^2\big)^{\!\sfrac{n}{2}},
\end{equation}
where $\opp = -i\hbar\nabla$ so that $\n{\opp}^2 = -\hbar^2\Delta$.

\subsubsection{A refined version of the main theorem} With the above definitions, we now give another version of the main theorem with more explicit upper bounds on the constants appearing in Inequality~\eqref{eq:Wigner_main_estimate}.

\begin{theorem}\label{thm:main_2}
Under the assumptions of Theorem~\ref{thm:main}, the following estimate holds
\begin{equation}\label{eq:Weyl_main_estimate}
\Nrm{\op-\op_f}{\L^2} \leq \cC_{\infty}^{\sfrac12} \(\Nrm{\widetilde{\op}_{\sqrt{f^\init}} - \sqrt{\op^\init}}{\L^2} + C_1(t) \,\hbar\) e^{\Lambda(t)} + C_2(t)\,\hbar.
\end{equation}
An upper bound for $\Lambda$ is given by $\Lambda(t) \leq C \int_0^t \lambda(s)\d s$ where $\lambda$ is defined for some $\eps\in(0,1)$ by
\begin{equation}\label{eq:lambda}
\lambda(t) = \Nrm{\Dhv \tilde{v}}{\cW^{1,2}} \Nrm{\rho_f}{L^\infty_x}^{\sfrac12} + \cC_{\infty}^{\sfrac12} \Nrm{\Dhv \tilde{v}\, \opm}{\L^{3\pm\eps}}
\end{equation}
with $n > 2$, and where $\Nrm{\cdot}{\L^{3\pm \eps}}$ is a shorthand notation for $\Nrm{\cdot}{\L^{3+\eps}}+\Nrm{\cdot}{\L^{3-\eps}}$. The functions $C_1(t)$ and $C_2(t)$ are bounded from above by
\begin{align}\label{eq:bound-c}
C_1(t)^2 &\leq C \int_0^t e^{C_E(s)} \Nrm{\Dhv \tilde{v}(s)}{\cW^{1,2}}^2 \n{\int_0^s C^\init + C_f(\tau) \bNrm{\op_{f(\tau,\cdot)} \weight{\opp}^2}{\cW^{2,2}}\d\tau}^2\d s
\\
C_2(t) &\leq C \biggl(C^\init + \int_0^t C_{f}(s) \Nrm{\nabla_{\xi}^2 f(s,\cdot)}{L^2} \d s\biggr),\label{eq:bound-c2}
\end{align}
which remain bounded at any time $t\geq 0$, and where $C_f(t) = \Nrm{\rho_{f(t,\cdot)}}{W^{1,\infty}_x\cap L^1_x}$,
\begin{equation}\label{eq:C-t}
C_E(t) = 2\int_0^t\Nrm{E_f(s,\cdot)}{L^\infty_x}\d s
\end{equation}
and
\begin{equation}\label{eq:C_init}
C^\init = \bNrm{\sqrt{f^\init}}{W^{1,4}_1\cap H^3}^2 + \bNrm{f^\init}{H^1_1\cap H^2}.
\end{equation}
\end{theorem}

The fact that these quantities~\eqref{eq:lambda}, \eqref{eq:bound-c} and \eqref{eq:bound-c2} remain bounded uniformly in $\hbar$ at any time is proved in the last section in Proposition~\ref{prop:propagation_regularity}.

\begin{remark}\label{rem:initial_data}
Notice that the quantity
\begin{equation*}
\Nrm{\widetilde{\op}_{\sqrt{f^\init}} - \sqrt{\op^\init}}{\L^2}= \Nrm{\op_{g_h*\sqrt{f^\init}} - \sqrt{\op^\init}}{\L^2} = \Nrm{g_h*\sqrt{f^\init} - f_{\sqrt{\op^\init}}}{L^2}
\end{equation*}
that appears in the right-hand side of Equation~\eqref{eq:Weyl_main_estimate} is $0$ in the case when $\op^\init = \widetilde{\op}_{\sqrt{f^\init}}^2$, that is when $\op^\init$ is the square of a Wick quantization. Due to the strategy of the proof, that mimics the proof of the stability estimates presented in Section~\ref{sec:stability-est}, the initial datum for the auxiliary problem \eqref{eq:linear-Hartree} has to be chosen positive, close to $\op_{f^\init}$ in~$\L^2$, its square root is close to $\sqrt{\op^\init}$ in $\cL^2$ and regular in the sense that the quantity
$\bNrm{\Dhv{\sqrt{\smash[t]{\widetilde{\op}}^{\init}}}\,\opm}{\L^p} + \bNrm{\Dhx{\sqrt{\smash[t]{\widetilde{\op}}^{\init}}}\,\opm}{\L^p}$ is uniformly bounded with respect to~$\hbar$. The choice $\widetilde{\op}^\init = \widetilde{\op}_{\sqrt{f^\init}}^2$ guarantees these properties.
\end{remark}

\section{Stability estimates for the Vlasov--Poisson and Hartree equations}\label{sec:stability-est}

\subsection{Classical case}

We start by explaining our method through examining the classical case and obtaining a stability estimate for the Vlasov--Poisson equation.

\begin{theorem}
Let $f_1$ and $f_2$ be two positive solutions of the Vlasov--Poisson equation verifying $\Nrm{f_1}{L^\infty}, \Nrm{f_2}{L^\infty} \leq \cC_{\infty}$. Then we have the following bound
\begin{equation*}
\bNrm{\sqrt{f_1} - \sqrt{f_2}}{L^2} \leq \bNrm{\sqrt{f_1^{\smash{\init}}} - \sqrt{f_2^{\smash{\init}}}}{L^2} e^{\Lambda(t)},
\end{equation*}
where $\Lambda(t) = C\,\int_0^t \lambda(s)\d s$ with
\begin{equation*}
\lambda(t) = \Nrm{{\rho_{f_2}}}{L^{\infty}_x}^{\sfrac12} \bNrm{\Dv \sqrt{f_2}}{L^3_xL^2_\xi} + \cC_{\infty}^{\sfrac12} \bNrm{\Dv \sqrt{f_2}}{L^{3,1}_xL^1_\xi}.
\end{equation*}
In particular, the following estimate holds
\begin{equation*}
\Nrm{f_1-f_2}{L^2} \leq 2\, \cC_{\infty}^{\sfrac12} \Nrm{f_1^\init-f_2^\init}{L^1}^{\sfrac12} e^{\Lambda(t)}.
\end{equation*}
\end{theorem}

\begin{proof}
Let $v_1=\sqrt{f_1}$, $v_2=\sqrt{f_2}$ and $v := v_1-v_2$. Then $v$ satisfies the equation
\begin{equation*}
\big(\dpt+\xi\cdot\Dx+E_{v_1^2}\cdot\Dv\big)v = \big(E_{v_2^2}-E_{v_1^2}\big)\cdot\Dv v_2 = \(\nabla K * \rho_{\(v_2+v_1\)v}\)\cdot\Dv v_2.
\end{equation*}
Then, by direct computation, we have that
\begin{align*}
\frac{1}{2}\frac{\d}{\d t}\Nrm{v}{L^2}^2 &= \intdd v \(\nabla K * \rho_{\(v_2+v_1\)v}\)\cdot\Dv v_2\d x\d\xi
\\
&= \int_{\R^{12}} \big(\n{v(z')}^2 + 2\, v_2(z')\,v(z')\big) \nabla K(x-x') \cdot \(v(z)\, \Dv v_2(z)\) \d z\d z'\\
&=: I_1 + 2 \, I_2.
\end{align*}
The first term is bounded by writing first
\[
I_1\leq \BNrm{\nabla K * \intd v\, \Dv v_2 \d\xi}{L^\infty_x} \Nrm{v}{L^2}^2 \leq \BNrm{\n{\nabla K} * \intd \n{\Dv v_2} \d\xi}{L^\infty_x} \Nrm{v}{L^\infty} {\Nrm{v}{L^2}^2}
\]
and then by applying Hölder's inequality for the Lorentz spaces to get
\begin{equation*}
I_1 \leq \Nrm{\nabla K}{L^{3/2,\infty}_x} \Nrm{\Dv v_2}{L^{3,1}_x L^1_\xi} \Nrm{v}{L^\infty} {\Nrm{v}{L^2}^2}.
\end{equation*}
The second term is bounded by the Hardy--Littlewood--Sobolev inequality and then by the Cauchy--Schwarz inequality, leading to
\begin{align*}
I_2 &\leq C\, \BNrm{\intd v\,v_2\d \xi\,}{L^2_x}\,\BNrm{\intd v\,\Dv v_2 \d \xi\,}{L^{6/5}_x}
\\
&\leq C \Nrm{\Nrm{v}{L^2_\xi}\Nrm{v_2}{L^2_\xi}}{L^{2}_x} {\Nrm{\Nrm{v}{L^2_\xi}\Nrm{\Dv v_2}{L^2_\xi}}{L^{6/5}_x}}.
\end{align*}
Finally, by Hölder's inequality, we obtain
\begin{align*}
I_2 &\leq C \Nrm{v_2}{L^\infty_xL^2_\xi}\,\Nrm{\Dv v_2}{L^3_xL^2_\xi}\, \Nrm{v}{L^2}^2 = C \Nrm{{\rho_{f_2}}}{L^\infty_x}^{\sfrac12}\, \Nrm{\Dv v_2}{L^3_xL^2_\xi}\, {\Nrm{v}{L^2}^2}.
\end{align*}
Combining the bounds for $I_1$ and $I_2$ leads to
\begin{equation*}
\frac{\d}{\d t} \Nrm{v}{L^2}^2 \leq C\, \Bigl(\Nrm{{\rho_{f_2}}}{L^\infty_x}^{\sfrac{1}{2}} \Nrm{\Dv v_2}{L^3_xL^2_\xi} + \Nrm{v}{L^\infty} \Nrm{\Dv v_2}{L^{3,1}_x L^1_\xi}\Bigr) {\Nrm{v}{L^2}^2}.
\end{equation*}
We conclude by using the fact that $\Nrm{v}{L^\infty} \leq \Nrm{f_1}{L^\infty}^{\sfrac12} + \Nrm{f_2}{L^\infty}^{\sfrac12} = 2\, \cC_{\infty}^{\sfrac12}$ and Grönwall's lemma.
\end{proof}

\subsection{Quantum case}\label{sect:Qcase}

Using the quantum analogue of gradients and Lebesgue norms of the phase space, we deduce a similar estimate for the Hartree equation.

Let us begin by recalling the following useful semiclassical commutator estimate.
\begin{lem}[Proposition 4.3 in \cite{lafleche_strong_2021}]\label{lem:semiclassical_commutator_est}
Let $p \!\in\! [1, \sfrac{3}{2})$ and $q$ satisfying \hbox{$p^{-1}\!\!=\!q^{-1}\!+\!\sfrac23$}. Then for any $\varepsilon \in (0, q-1)$ and $n>2$, there exists a constant $C>0$ such that
\begin{align*}
\Nrm{[K(\cdot-x), \op]}{\cL^p} \le C\, h\, \Nrm{\Dhv \op\, \opm}{\cL^{q+\varepsilon}}^{\Psfrac{1}{2}+\tilde \epsilon}\aNrm{\Dhv \op\, \opm}{\cL^{q-\varepsilon}}^{\Psfrac{1}{2}+\tilde \varepsilon},
\end{align*}
where $\tilde \varepsilon = \varepsilon/q$ and $\opm=\weight{\opp}^n$.
\end{lem}
Notice that our scaling is different from the scaling used in \cite[Prop.\,4.3]{lafleche_strong_2021}, but one can easily pass from that scaling to our current scaling by multiplying both side of the inequality in \cite{lafleche_strong_2021} by $h^d$ with $d=3$.
\begin{theorem}\label{thm:Q-uniqueness}
Let $\op_1$ and $\op_2$ be two solutions of the Hartree equation~\eqref{eq:Hartree} such that $\Nrm{\op_1}{\L^\infty} \leq \cC_{\infty}$ and $\Nrm{\op_2}{\L^\infty} \leq \cC_{\infty}$ and let $v_1 = \sqrt{\op_1}$ and $v_2 = \sqrt{\op_2}$. Then there exists a universal constant $C>0$ and $T>0$ such that for any $t\in [0,T]$,
\begin{equation}\label{eq:L^2-L^2}
\nrm{v_1-v_2}_{\L^2} \leq {\Nrm{v_1^\init-v_2^\init}{\L^2} e^{\Lambda(t)}},
\end{equation}
where $\Lambda(t) = C \int_0^t \lambda(s)\d s$ with $\lambda$ given for some $\eps\in(0,1)$ by
\begin{equation}\label{eq:lambda_2}
\lambda(t) = \Nrm{\Dhv v_2}{\cW^{1,2}} \bNrm{{\rho_{v_2^2}}}{L^\infty_x}^{\sfrac12} + \cC_{\infty}^{\sfrac12} {\Nrm{\Dhv v_2\,\opm}{\L^{3\pm\eps}}},
\end{equation}
{where $n>2$}. This implies the following estimate
\begin{equation}\label{eq:L^2-L^1}
\Nrm{\op_1-\op_2}{\L^2} \leq 2\, \cC_{\infty}^{\sfrac12} \Nrm{\op_1^\init-\op_2^\init}{\L^1}^{\sfrac12} e^{\Lambda(t)}.
\end{equation}
\end{theorem}

\begin{remark}
It is actually not difficult to see from the interpolation argument in the proof of Lemma~\ref{lem:commutator_estimate_HS} that $\Nrm{\Dhv v_2}{\cW^{1,2}} = \Nrm{\Dv f_{v_2}}{H^1}$ can actually be replaced by $\Nrm{\Dv f_{v_2}}{H^{1/2}}$.
\end{remark}

\begin{remark}
Contrarily to our main theorem, the result here is only local-in-time since it is not yet known whether the global-in-time uniform in $\hbar$ propagation of regularity holds for the Hartree equation in the case of the Coulomb potential. For slightly less singular potentials, we have proved it in our recent paper~\cite{chong_global--time_2022}.
\end{remark}

\begin{proof}
Let $v = v_1-v_2$. Then $v$ satisfies the equation
\begin{equation*}
i\hbar\,\dpt v = \com{H_{\op_1}, v} + \bcom{V_{\op_1}-V_{\op_2}, v_2}.
\end{equation*}
Differentiating its Hilbert--Schmidt norm with respect to time, we obtain
\begin{align*}
\frac{1}{4\pi}\, \frac{\d}{\d t} \Nrm{v}{\L^2}^2 &= h^2 \Tr{\com{V_{\op_1}-V_{\op_2}, v_2} v}
= h^2 \bTr{\bcom{V_{v_1^2-v_2^2}, v_2} v}.
\end{align*}
Since $v_1^2 - v_2^2 = \(v+v_2\)^2 - v_2^2 = v^2 + v\,v_2 + v_2\,v$ and $\Diag{v\,v_2} = \conj{\Diag{v_2\, v}}$, we get
\begin{align*}
\frac{1}{4\pi}\,\frac{\d}{\d t} \Nrm{v}{\L^2}^2 = h^2 \Tr{\com{V_{v^2}, v_2} v} + 2\,h^2 \Tr{\com{V_{v_2\,v}, v_2} v} =: I_1 + 2\, I_2.
\end{align*}
To bound the first term, we introduce the notation $K_x := K(x-y)$ and write
\begin{multline*}
I_1 = h^2 \intd \Tr{\com{K_x,v_2} v} \rho_{v^2}(x) \d x \leq \sup_{x\in\Rd} \(h^{-1}\Nrm{\com{K_x,v_2}}{\L^1} \Nrm{v}{\L^\infty}\) \intd \n{\rho_{v^2}} \d x
\\
\leq C {\Nrm{\Dhv{v_2}\,\opm}{\L^{3\pm\eps}}} \Nrm{v}{\L^\infty} {\Nrm{v}{\L^2}^2},
\end{multline*}
where we have used Lemma \ref{lem:semiclassical_commutator_est} to bound the commutator in trace norm. To bound the second term, we use Hölder's inequality to get
\begin{equation*}
I_2 \leq h^{-1}\Nrm{\com{V_{v\,v_2}, v_2}}{\L^2} \Nrm{v}{\cL^2}
\end{equation*}
and then we apply the following inequality which we will prove in the subsequent section (see Lemma~\ref{lem:commutator_estimate_HS})
\begin{equation*}
\Nrm{\com{V_{v\,v_2},v_2}}{\L^2} \leq C \Nrm{\Dhv v_2}{\cW^{1,2}} {\Nrm{\rho_{v\,v_2}}{L^2_x}}.
\end{equation*}
Also, notice that
\begin{align*}
\Nrm{\rho_{v\,v_2}}{L^2_x} &= h^3 \(\intd \Big|\intd v(x,y)\,v_2(y,x)\d y\Big|^2\d x\)^{\!\sfrac{1}{2}} \leq h^{\sfrac{3}{2}} \Nrm{v_2}{L^\infty_xL^2_y} {\Nrm{v}{\L^2}}.
\end{align*}
Hence, since $h^3\Nrm{v_2}{L^2_y}^2 = \Diag{\op_2} = \rho_2$, we deduce that
\begin{equation*}
I_2 \leq C\Nrm{\Dhv v_2}{\cW^{1,2}} \Nrm{\rho_2}{L^\infty_x}^{\sfrac12} {\Nrm{v}{\L^2}^2},
\end{equation*}
and, as in previous proposition, applying Grönwall's lemma proves~\eqref{eq:L^2-L^2}. Inequality~\eqref{eq:L^2-L^1} follows from the identity $\op_1-\op_2=\frac{1}{2}(v_1-v_2)(v_1+v_2)+\frac{1}{2}(v_1+v_2)(v_1-v_2)$ and the Powers--St{\o}rmer inequality
\begin{align*}
\bNrm{\sqrt{A}-\sqrt{B}}{\cL^2}^2\le {\aNrm{A-B}{\cL^1}},
\end{align*}
which holds for any positive operators $A, B$
(see~\cite[Lem.\,4.1]{powers_free_1970}).
\end{proof}

\Subsubsection{A Semiclassical inequality for commutators}

\begin{lem}\label{lem:commutator_estimate_HS}
There exists a constant $C>0$ such that for any self-adjoint trace class operators $\op$ and $\opmu$, we have the following estimate
\begin{equation*}
\frac{1}{\hbar} \Nrm{\com{V_{\op},\opmu}}{\L^2} \leq C \Nrm{\rho}{L^2_x} \Nrm{\Dhv \opmu}{\cW^{1,2}}
\end{equation*}
where $\rho = \Diag{\op}$.
\end{lem}

\begin{proof}[Proof of Lemma~\ref{lem:commutator_estimate_HS}]
We use a first order Taylor expansion of $V_{\op}$ of the form
\begin{equation*}
V_{\op}(x) = V_{\op}(y) + \(x-y\)\cdot \int_0^1 \nabla V_{\op}(z_\theta)\d\theta,
\end{equation*}
where $z_\theta = \(1-\theta\)x+\theta\, y$. This implies that $\frac{1}{i\hbar}\com{V_{\op},\opmu}$ is the operator with kernel
\begin{equation*}
\frac{1}{i\hbar}\com{V_{\op},\opmu}(x,y) = \int_0^1 \nabla V_{\op}(z_\theta) \cdot \Dhv \opmu(x,y) \d\theta.
\end{equation*}
Its Hilbert--Schmidt norm being given by the $L^2$ norm of its kernel, and doing the change of variables $(x,y) \mto (x+\theta y, x-(1-\theta) y)$ with Jacobian equal to $1$, we obtain
\begin{align*}
\Nrm{\tfrac{1}{i\hbar}\com{V_{\op},\opmu}}{\L^2}^2 &\leq {h^3}\int_0^1\!\!\! \iintd \n{\nabla V_{\op}(x) \cdot \Dhv \opmu(x+\theta y,x-(1-\theta) y)}^2\d x\d y \d\theta
\\[-3pt]
&\leq {h^3} \Nrm{\nabla V_{\op}}{L_x^6}^2 \int_0^1\!\!\! \intd \Nrm{\Dhv \opmu(x+\theta y,x-(1-\theta) y)}{L^3_x}^2 \d y \d\theta,
\end{align*}
where the second inequality follows by Hölder's inequality. The first factor on the right-hand side is controlled using the Hardy--Littlewood--Sobolev inequality by $\Nrm{\rho}{L^2_x}$. To~control the second factor, we perform the change of variable $x\mto x-\theta y$ to get\vspace*{-3pt}
\begin{equation*}
\int_0^1\!\!\! \intd \Nrm{\Dhv \opmu(x+\theta y,x-(1-\theta) y)}{L^3_x}^2 \d y \d\theta = \intd \Nrm{\Dhv \opmu(x,x-y)}{L^3_x}^2 \d y.
\end{equation*}
The $L^3$ norm in this term is now bounded using the Gagliardo--Nirenberg--Sobolev inequality associated to the embedding $H^1\subset L^3$ which yields
\begin{equation*}
\Nrm{\Dhv \opmu(x,x-y)}{L^3_x}^2 \leq C \Bigl(\Nrm{\Dhv\opmu(x,x-y)}{L^2_x}^2 + \Nrm{\Dx(\Dhv \opmu(x,x - y))}{L^2_x}^2\Bigr).
\end{equation*}
However, notice now that the gradient with respect to $x$ appearing in the above formula is actually also a quantum gradient of the operator since
\begin{equation*}
\nabla_x(\Dhv\opmu(x,x-y)) = ((\nabla_1+\nabla_2)\Dhv\opmu)(x,x-y) = \com{\nabla,\Dhv{\opmu}}(x,x-y)
\end{equation*}
and this is exactly $\Dhx{\Dhv{\opmu}}(x,x-y)$. Using the fact that for both $\opnu = \Dhv\opmu$ and for $\opnu = \Dhx\Dhv\opmu$, by the Fubini theorem and the change of variables $y\mto x-y$ it holds\vspace*{-3pt}
\begin{align*}
{h^3} \intd \Nrm{\opnu(x,x-y)}{L^2_x}^2 \d y = {h^3} \iintd \n{\opnu(x,x-y)}^2 \d x\d y = {\Nrm{\opnu}{\L^2}^2},
\end{align*}
we therefore arrive at the inequality\vspace*{-3pt}
\begin{equation*}
\int_0^1\!\!\! \intd \Nrm{\Dhv \opmu(x+\theta y,x-(1-\theta) y)}{L^3_x}^2 \d y \d\theta \leq C {\Nrm{\Dhv \opmu}{\cW^{1,2}}},
\end{equation*}
which concludes the proof.
\end{proof}
\medskip

\section{Semiclassical limit in \texorpdfstring{$\L^2$}{L2} Schatten norm}

We consider $\op$ a solution of the Hartree equation~\eqref{eq:Hartree} and $f$ a solution of the Vlasov--Poisson equation~\eqref{eq:Vlasov}. Then, in the same spirit as~\cite{lafleche_strong_2021}, we will use the fact that the Weyl transform $\op_f$ of $f$ solves the equation\vspace*{-3pt}
\begin{equation}\label{eq:Weyl-Vlasov}
i\hbar\,\dpt \op_f = \com{H_{f},\op_f} - B_f(\op_f) \quad \text{ with } \quad \op_f(0) = \op_{f^\init},
\end{equation}
where $B_f(\op_f)$ is the operator with kernel\vspace*{-3pt}
\begin{equation}\label{def:B_operator}
B_f(\op_f)(x,y) = \(V_{f}(x) - V_{f}(y) - \(x-y\)\cdot\nabla V_{f}\Big(\frac{x+y}{2}\Big)\!\) \op_f(x,y).
\end{equation}
A second order Taylor expansion (see the proof of \cite[Prop.\,4.4]{lafleche_strong_2021}) leads to the following estimate\vspace*{-3pt}
\begin{equation}\label{eq:B_bound}
\frac{1}{\hbar}\Nrm{B_f(\op_f)}{\L^2} \leq C\,\hbar\Nrm{\nabla E_f}{L^\infty_x} {\Nrm{\nabla_\xi^2 f}{L^2}},
\end{equation}
where $\Nrm{\nabla E_f}{L^\infty_x}$ is controlled by $\Nrm{\rho_f}{W^{1,\infty}_x}$ and $\Nrm{\rho_f}{L^1_x}$.

Let $\widetilde{\op}$ be a solution to the following Cauchy problem for the linear Hartree equation
\begin{equation}\label{eq:linear-Hartree}
i\hbar\,\dpt\widetilde{\op} = \com{H_f,\widetilde{\op}} \quad \text{ with } \quad \widetilde{\op}(0) = \widetilde{\op}^\init := \widetilde{\op}_{\sqrt{f^\init}}^2,
\end{equation}
with Hamiltonian $H_f=-\sfrac{\hbar^2\Delta}{2}+V_{f}$.

\begin{remark}
The choice $\widetilde{\op}^\init = \widetilde{\op}_{\sqrt{f^\init}}^2$ guarantees
\begin{equation*}
\bNrm{\Dh_{\!\eta} \widetilde{\op}_{\sqrt{f^\init}} \,\opm}{\L^p} = \bNrm{\op_{g_h*\nabla_\eta\sqrt{f^\init}}\,\opm}{\L^p},
\end{equation*}
which is bounded by a Sobolev norm of $\sqrt{f^{\init}}$ by \cite[Prop.\,3.2]{lafleche_strong_2021}. See also Remark~\ref{rem:initial_data}.
\end{remark}

Our aim is to show that the bound on the difference between $\op$ and $\op_f$ in $\L^2$ is of order $\hbar$. To this end, we consider $\widetilde{\op}$ the solution to the auxiliary problem \eqref{eq:linear-Hartree} with initial datum $\widetilde{\op}^\init = \widetilde{\op}_{\sqrt{f^\init}}^2$. Then, applying Minkowski's inequality yields
\begin{equation}\label{eq:op-opf}
\Nrm{\op-\op_f}{\L^2} \leq \Nrm{\op-\widetilde{\op}}{\L^2} + {\Nrm{\widetilde{\op}-\op_f}{\L^2}}.
\end{equation}
The second term on the right-hand side is needed in our method because of the possible lack of positivity of the Weyl quantization of $f$, which prevents us to take its square root. Hence we will call it the term due to the lack of positivity. The equations verified by $\widetilde{\op}$ and $\op_f$ yields
\begin{equation*}
i\hbar\dpt\(\widetilde{\op}-\op_f\) = \com{H_f,\widetilde{\op}-\op_f} + B_f(\op_f).
\end{equation*}
Hence, Inequality~\eqref{eq:B_bound} leads to
\begin{equation}\label{eq:optilde-opf}
\Nrm{\widetilde{\op}-\op_f}{\L^2} \leq \bNrm{\widetilde{\op}^\init-\op_{f^\init}}{\L^2} + C\, \hbar \int_0^t\Nrm{\nabla E_f(s,\cdot)}{L^\infty_x}\Nrm{\nabla_{\xi}^2 f(s,\cdot)}{L^2} \d s.
\end{equation}
The first term on the right-hand side of Inequality~\eqref{eq:op-opf} is the true nonlinear error that corresponds to the stability estimate for the Hartree equation proved in the previous section. It can be readily estimated in a similar manner as in Section~\ref{sect:Qcase} by
\begin{equation}\label{eq:sqrt-op_optilde}
\Nrm{\op-\widetilde{\op}}{\L^2} \leq \bNrm{\sqrt{\op}-\sqrt{\widetilde{\op}}\,}{\L^2} \bigl(\Nrm{\sqrt{\op}}{\L^\infty} + \bNrm{\sqrt{\widetilde{\op}}\,}{\L^\infty}\bigr).
\end{equation}
We now estimate each term separately.

\subsection{An estimate on the term due to the lack of positivity}

For the second term on the right-hand side of Equation \eqref{eq:optilde-opf}, $\Nrm{\nabla E_f}{L^\infty_x}$ and $\aNrm{\nabla^2_\xi f}{L^2}$ are bounded by \cite[Prop.\,A.1]{lafleche_strong_2021}. As for the first term, the difference between the initial data $\widetilde{\op}^\init$ and $\op_{f^\init}$ can be written as follows
\begin{equation*}
\widetilde{\op}^\init - \op_{f^\init} = \widetilde{\op}_{\sqrt{f^\init}}^2 - \widetilde{\op}_{f^\init} + \widetilde{\op}_{f^\init} - \op_{f^\init}.
\end{equation*}
Since, by Inequality~\eqref{eq:Wick_Weyl_error}, we have that
\begin{equation}\label{eq:Wick_Weyl_error_init}
\bNrm{\widetilde{\op}_{f^\init} - \op_{f^\init}}{\L^2} \leq \frac{3\,\hbar}{2}\, {\Nrm{\nabla^2 f^\init}{L^2}},
\end{equation}
it remains to estimate $\widetilde{\op}_{\sqrt{f^\init}}^2 - \widetilde{\op}_{f^\init}$ in $\L^2$. This is done via the following lemma.

\begin{lem}\label{lem:commut_square_wick}
{Let $g\in W^{1,2p}$ with $p\in[1,\infty]$.} Then
\begin{equation*}
\bNrm{\widetilde{\op}_{g^2} - \widetilde{\op}_g^2}{\L^p} \leq 48\,\hbar\, {\Nrm{\nabla g}{L^{2p}}^2},
\end{equation*}
where $\widetilde{\op}_g^2$ means $(\widetilde{\op}_g)^2$.
\end{lem}

\begin{proof}
Notice that for any $(f,g)\in (L^1+ L^\infty)^2$,
\begin{equation*}
\widetilde{\op}_{f}\,\widetilde{\op}_{g} = h^{-6} \iintd f\,g' \,G(z,z') \, \ket{\psi_{z}}\bra{\psi_{z'}} \d z \d z',
\end{equation*}
where $f = f(z), g' = g(z')$ and
\begin{equation*}
G(z,z') = \Inprod{\psi_z}{\psi_{z'}} = \Bigl(\frac{2}{h}\Bigr)^{\sfrac32} \intd e^{-(\n{y-x}^2+\n{y-x'}^2)/(2\hbar)}\, e^{i\,y\cdot\(\xi'-\xi\)/\hbar} \d y.
\end{equation*}
Using the parallelogram identity and the formula for the Fourier transform of a Gaussian, one obtains the following expression for $G$,
\begin{equation*}
G(z,z')
= e^{-\an{\psfrac{z-z'}{2}}^2/\hbar}\, e^{i\,\(x+x'\)\cdot\(\xi'-\xi\)/(2\hbar)}.
\end{equation*}
Notice also that $\widetilde{\op}_1 = \mathds{1}$ is the identity operator. Hence, we write
\[
\widetilde{\op}_{fg} = \frac{\widetilde{\op}_{fg}\, \widetilde{\op}_1 + \widetilde{\op}_1\, \widetilde{\op}_{fg}}{2} = h^{-6} \iintd \frac{f'g'+fg}{2} \,G(z,z') \ket{\psi_{z}}\bra{\psi_{z'}} \d z \d z'.
\]
Therefore, since $f\,g + f'\,g' - f\,g' - f'\,g = \(f-f'\) \(g-g'\)$, we obtain
\begin{align*}
\Lambda(f,g) &:= \widetilde{\op}_{fg} - \frac{\widetilde{\op}_f\,\widetilde{\op}_g + \widetilde{\op}_g\,\widetilde{\op}_f}{2}
\\
&\hphantom{:}= h^{-6} \iintd \frac{\(f-f'\)\(g-g'\)}{2} \,G(z,z') \ket{\psi_{z}}\bra{\psi_{z'}} \d z \d z'.
\end{align*}
To obtain bounds on $\Lambda(g,g) = \widetilde{\op}_{g^2} - \widetilde{\op}_g^2$, we will now bound $\Lambda(f,g)$ by bilinear interpolation.

\step{Trace norm estimate} Using Minkowski's inequality and the fact that the trace norm of $\ket{\psi_{z}}\bra{\psi_{z'}}$ is 1, we obtain the following trace norm estimate
\begin{multline*}
\Nrm{\Lambda(f,g)}{\L^1} \leq \frac{h^{-3}}{2} \iint_{[0,1]^2}\iintdd \n{\nabla f(z_\theta)} \n{\nabla g(z_{\theta'})} G_2(z-z') \d z \d z' \d \theta\d \theta'
\\
\leq \frac{h^{-3}}{2} \(\iintdd \n{\nabla f(u)}^2 G_2(v) \d u \d v\)^{\sfrac{1}{2}} \(\iintdd \n{\nabla g(u)}^2 G_2(v) \d u \d v\)^{\sfrac{1}{2}}
\\
\leq 48\,\hbar \Nrm{\nabla f}{L^2} {\Nrm{\nabla g}{L^2}},
\end{multline*}
where $G_2(z) = \n{z}^2 e^{-\an{\sfrac{z}{2}}^2/\hbar}$ and $z_\theta = \(1-\theta\) z + \theta\, z'$.
\step{Operator norm estimate} Let $(\varphi,\phi)\in (L^2)^2$ and define $\psi_\varphi(z) := \Inprod{\psi_z}{\varphi}$. Then
\begin{align*}
\Inprod{\phi}{\Lambda(f,g)\,\varphi} &= h^{-6} \iintdd \frac{\(f-f'\)\(g-g'\)}{2} \,G(z,z') \, \psi_{\phi}(z)\,\psi_{\varphi}(z') \d z \d z'
\\
&\leq h^{-6} \Nrm{\nabla f}{L^{\infty}}\Nrm{\nabla g}{L^{\infty}}\iintdd G_2(z-z') \n{\psi_{\phi}(z)}\n{\psi_{\varphi}(z')} \d z \d z'
\\
&\leq h^{-6} \Nrm{\nabla f}{L^{\infty}}\Nrm{\nabla g}{L^{\infty}}\Nrm{G_2}{L^1} \Nrm{\psi_{\phi}}{L^2}{\Nrm{\psi_{\varphi}}{L^2}},
\end{align*}
where the last inequality follows from Young's inequality. Observing that $\psi_\varphi(z) = \F{\psi_0(x-\cdot)\,\varphi}(\xi/h)$, we deduce that $\Nrm{\psi_\varphi}{L^2} = h^{\sfrac32} \Nrm{\varphi}{L^2}$ and so
\begin{equation*}
\Nrm{\Lambda(f,g)}{\L^\infty} \leq 48\,\hbar \Nrm{\nabla f}{L^{\infty}}{\Nrm{\nabla g}{L^{\infty}}}.
\end{equation*}

\step{General Schatten norms} By the complex bilinear interpolation (see \eg \cite[\S 4.4]{bergh_interpolation_1976}), we deduce that for any $p\in [1,\infty]$,
\begin{equation*}
\Nrm{\Lambda(f,g)}{\L^p} \leq 48\,\hbar \Nrm{\nabla f}{L^{2p}}{\Nrm{\nabla g}{L^{2p}}},
\end{equation*}
which leads to the desired result by taking $f=g$.
\end{proof}

\subsection{Bound on \texorpdfstring{$\Nrm{\op-\protect\widetilde{\op}}{\L^2}$}{auxiliary and Hartree} }\label{sec:op-opt}

In Equation \eqref{eq:sqrt-op_optilde}, $\aNrm{\sqrt{\op}}{\L^\infty}$ and $\aNrm{\sqrt{\widetilde{\op}}}{\L^\infty}$ are bounded because $\op\in\L^\infty$. The main result of this section is Proposition~\ref{prop:sqrt-op-vs-opt} which gives an estimate for $\aNrm{\sqrt{\op}-\sqrt{\widetilde{\op}}}{\L^2}$ in terms of $\hbar$. To prove the proposition, we come back to the equations satisfied by $\op$ and $\tilde{\op}$ and proceed similarly as in the proof of Theorem~\ref{thm:Q-uniqueness} about the stability for the Hartree equation. However, the equation of $\tilde{\op}$ is different from the Hartree equation as its force field is given by the classical force field. This makes appear several additional error terms. To bound these error terms, we first establish some preliminary results.

As a first tool, it will be useful to exchange the role of $x$ and $\xi$ to manipulate weights of the form $\weight{x}$ instead of weights of the form $\weight{\opp}$. The quantum analogue of exchanging the $x$ and $\xi$ variables is obtained by conjugation with the semiclassical Fourier transform $\ecF_h\varphi(\xi):= h^{\sfrac32}\,\widehat{\varphi}(h\xi)$. More precisely, we define $\op^\star := \ecF^{-1}_h\op\, \ecF_h$. Then it holds
\begin{equation}\label{eq:property_op_star}
\op_f^\star = \op_{f^\star} \quad \text{ with } \quad f^\star(x, \xi) = f(\xi, x).
\end{equation}
This exchange operation is a linear automorphism which preserves the $\L^p$ norms. Moreover, from the definition of a Fourier multiplier, $\weight{\opp}\ecF_h\!=\!\ecF_h\weight{x}$ and so \hbox{$\weight{x}^\star \!=\! \weight{\opp}$} and $\weight{\opp}^\star = \weight{x}$.

\begin{proof}[Proof of~\eqref{eq:property_op_star}]
Computing the kernel of $\op_f^\star$ yields
\begin{equation*}
\op_f^\star(x,y) = \frac{1}{h^3}\int_{\RR^9} e^{2\pi i\, \(x\cdot x'-y'\cdot y\)/h}\,e^{-2\pi i\, (y'-x')\cdot\xi}\, f\!\bigl(\tfrac{x'+y'}{2}, h\,\xi\bigr) \d x'\d y'\d \xi
\end{equation*}
Therefore, by the change of variables $\zeta = \psfrac{x'+y'}{2\,h}$ and $\eta = x'-y'$
\begin{equation*}
\op_f^\star(x,y) = \intdd e^{-2\pi i\, \(y-x\)\cdot\zeta/h} \, e^{2i\pi\eta\cdot\psfrac{x+y}{2\,h}}\,\ecF\!\(f(h\,\zeta, h\,\cdot)\)\!(\eta) \d \eta\d \zeta
\end{equation*}
and we deduce Equation~\eqref{eq:property_op_star} by the Fourier inversion theorem.
\end{proof}

We can now use the above defined operation to prove the following lemma.
\begin{lem}\label{lem:Weyl_and_weight}
{ Let $f$ be a phase space function. Then there exist $r, r_2 \in L^2(\Rdd)$ such that\vspace*{-3pt}
%\begin{subequations}
\begin{align*}
\op_{\weight{\xi} f}-\op_{f}\weight{\opp} &= \op_{r},
\\
\weight{\opp}\op_f \weight{\opp} -\op_{\weight{\xi}^2f} &= \op_{r_2},
\end{align*}
%\end{subequations}
satisfying the following estimates\vspace*{-3pt}
\begin{subequations}
\begin{align}\label{eq:p-weight_remainder}
\Nrm{r}{L^2} = \bNrm{\op_{\weight{\xi} f}-\op_{f}\weight{\opp}}{\L^2} &\le \frac{\hbar}{2}\,\aNrm{\grad_{x} f}{L^2},
\\\label{eq:p-weight_remainder_2}
\Nrm{r_2}{L^2} = \bNrm{\weight{\opp}\op_f \weight{\opp} -\op_{\weight{\xi}^2f}}{\L^2} &\le \frac{3\hbar^2}{4}\,\aNrm{\Delta_x f}{L^2}.
\end{align}
Notice that a similar result holds when we replace $\op_{f}\weight{\opp}$ by $\weight{\opp}\op_{f}$, which follows by taking the adjoint.
Furthermore, we also have
\begin{equation*}
\Nrm{r}{H^2} \leq C\,\hbar\Nrm{f}{H^3}
\end{equation*}
for some $C>0$ independent of $\hbar$.
\end{subequations}
}
\end{lem}

\begin{proof}
Notice the kernel of the difference $\op_{\weight{x}f}-\op_{f}\weight{x}$ is given by\vspace*{-5pt}
\begin{multline*}
\intd e^{-2\pi i \(y-x\)\cdot \xi}\(\weight{\tfrac{x+y}{2}}-\weight{y}\)f\(\tfrac{x+y}{2}, h\xi\)\d \xi
\\
= \frac{i\hbar}{2}\frac{(\frac{x+y}{2})+ y}{\weight{\tfrac{x+y}{2}}+\weight{y}}\cdot \intd e^{-2\pi i \(y-x\)\cdot \xi}\, \grad_{\xi}f\(\tfrac{x+y}{2}, h\xi\)\d \xi = \op_r,
\end{multline*}
where $r \in L^2$ is defined by $r(z)=\tfrac{i\hbar}{2}\,(a_x(i\hbar\, \grad_\xi)\cdot\grad_{\xi}f)(z)$ where $a_x(i\hbar\, \grad_\xi)$ is a bounded, $x$-dependent Fourier multiplier of the $\xi$ variable associated to the function $a_x(\eta) = \pspfrac{2\,x- \eta/2}{\weight{x}+\weight{x-\eta/2}}$. More explicitly, we have
\begin{equation*}
r(z) = \frac{i\hbar}{2}\intdd e^{2\pi i\(\xi-\xi'\)\cdot \eta}\, a_x(h\eta)\cdot \grad_{\xi}f(x, \xi')\d\xi'\d\eta.
\end{equation*}
In particular, it follows that\vspace*{-3pt}
\begin{equation}\label{eq:x-weight_remainder}
\bNrm{\op_{\weight{x} f}-\op_{f}\weight{x}}{\L^2} = \Nrm{\op_{r}}{\L^2} = \Nrm{r}{L^2} \le \frac{\hbar}{2}\, \aNrm{\grad_{\xi}f}{L^2}.
\end{equation}
We now use Inequality~\eqref{eq:x-weight_remainder} together with the properties of the exchange operation~\eqref{eq:property_op_star} to get\vspace*{-3pt}
\begin{align*}
\bNrm{\op_{\weight{\xi} f}-\op_{f}\weight{\opp}}{\L^2} &= \bNrm{\op_{\weight{\xi} f}^\star-\op_{f}^\star\weight{x}}{\L^2} = \bNrm{\op_{\weight{x} f^\star}-\op_{f^\star}\weight{x}}{\L^2}
\\
&\le \hbar \Nrm{\grad_{\xi}f^\star}{L^2} = \hbar \Nrm{(\grad_{x}f)^\star}{L^2} = \hbar\,{\Nrm{\grad_{x}f}{L^2}},
\end{align*}
which completes the proof of Inequality~\eqref{eq:p-weight_remainder}. To get the second inequality, a similar computation gives that
\begin{equation*}
\weight{x}\op_f \weight{x} -\op_{\weight{x}^2f} = \op_{r_2}
\end{equation*}
with $r_2(z)=\tfrac{-\hbar^2}{4}\,(b_x(i\hbar\, \grad_\xi)\Delta_{\xi}f)(z)$, where $b_x$ is the bounded function given by
\begin{equation*}
b_x(2\eta) = \frac{{4\,(\frac{\eta}{\n{\eta}}\cdot x)^2-2\n{x}^2}-2-\n{\eta}^2}{\weight{x}^2+\weight{x+\eta}\weight{x-\eta}} = \frac{{4\,(\frac{\eta}{\n{\eta}}\cdot x)^2-3\n{x}^2}-2+\(x+\eta\)\(x-\eta\)}{\weight{x}^2+\weight{x+\eta}\weight{x-\eta}},
\end{equation*}
and the same reasoning leads to Inequality~\eqref{eq:p-weight_remainder_2}.
\end{proof}

Another simple observation is that the commutator between convolution with the Gaussian function $g_h(z) = \(\pi\hbar\)^{-3}e^{-\n{z}^2/\hbar}$ and multiplication by a weight is also of order $\hbar$.
\begin{lem}\label{lem:p-weight_remainder_gaussian}
Let $h\in(0,1)$. Then there exists $C>0$ such that for any $p\in[1,\infty]$ and any function of the phase space it holds
\begin{subequations}
\begin{align}\label{eq:p-weight_remainder_gaussian}
\Nrm{\com{\weight{\xi},\sim}f}{W^{2,p}} \leq C\,\hbar\, {\Nrm{f}{W^{3,p}}},
\\\label{eq:p-weight_remainder_gaussian_2}
\Nrm{\com{\langle\xi\rangle^2,\sim}f}{L^p} \leq C\,\hbar\, {\Nrm{f}{W^{1,p}_1}},
\end{align}
\end{subequations}
where $\com{m,\sim} f := m\(g_h*f\) - \(g_h*(mf)\)$.
\end{lem}

\begin{proof}
For any $z=(x,\xi)\in\Rdd$ and any general $m(\xi)$, it holds
\begin{equation*}
\com{m,\sim} f(z) = {\intdd} g_h(z') \(m(\xi)-m(\xi-\xi')\) f(z-z')\d z',
\end{equation*}
where $z' = (x',\xi')$. By a second order Taylor expansion, this can be written as
\begin{align*}
\com{m,\sim} f &= {\intdd} g_h(z')\, \xi'\cdot\(\nabla m(\xi)- \int_0^1 {(1-\theta)}\, \xi'\cdot \nabla^2m(\xi-\theta\xi')\d\theta \) f(z-z')\d z'
\\
&= I_1 - I_2.
\end{align*}
Next, using the fact that that $g_h$ is even, one can replace $f(z-z')$ by $f(z-z')-f(z)$ in $I_1$, leading to
\begin{equation*}
I_1 = \nabla m(\xi)\cdot \int_0^1\!\!\!{\intdd} \xi'\,g_h(z')\, \xi' \cdot\,\Dv f({x-x',\xi-\theta \xi'})\d z' \d\theta.
\end{equation*}
By Minkowski's inequality, we obtain
\begin{equation*}
\Nrm{I_1}{L^p} \le \Nrm{\nabla m}{L^\infty} \,\bNrm{\!\n{\xi}^2g_h}{L^1} \Nrm{\Dv f}{L^p} = \frac{3}{2}\,\hbar\,\Nrm{\nabla m}{L^\infty}{\Nrm{\Dv f}{L^p}}.
\end{equation*}
Similarly, $I_2$ is bounded by
\begin{equation*}
\Nrm{I_2}{L^p} \le \frac{1}{2}\Nrm{\nabla^2 m}{L^\infty} \,\bNrm{\!\n{\xi}^2g_h}{L^1}\Nrm{f}{L^p} = \frac{3}{4}\,\hbar\,\Nrm{\nabla^2 m}{L^\infty}{\Nrm{f}{L^p}}.
\end{equation*}
In particular, if $m = \weight{\xi}$, then it yields $\Nrm{\com{m,\sim} f}{L^p} \leq 3\,\hbar \Nrm{f}{W^{1,p}}$. Noticing that $\Dx^2\com{m,\sim} f = \com{m,\sim} \Dx^2 f$ also yields $\Nrm{\Dx^2\com{m,\sim} f}{L^p} \leq 3\,\hbar \Nrm{\Dx^2 f}{W^{1,p}}$. Finally, noticing that
\begin{equation*}
\Dv^2\com{m,\sim} f = \com{m,\sim} \Dv^2 f + 2\com{\nabla m,\sim} \Dv f + \com{\nabla^2 m,\sim} f
\end{equation*}
and using the above estimates with $m$ replaced by $\nabla m$ and $\nabla^2 m$ leads to
\begin{equation*}
\Nrm{\Dv^2\com{m,\sim} f}{L^p} \leq 36 \,\hbar \Nrm{f}{H^3}
\end{equation*}
and so to Inequality~\eqref{eq:p-weight_remainder_gaussian}. To get Inequality~\eqref{eq:p-weight_remainder_gaussian_2}, notice that $\acom{\weight{\xi}^2,\sim}f = \acom{\n{\xi}^2,\sim}f$ and write
\begin{equation*}
\bcom{\n{\xi}^2,\sim}f = \intdd g_h(z')\(\xi'\)^{\otimes 2} : \(f(z-z')\,\Id + 2\, \xi\otimes \int_0^1 \Dv f(z-\theta \,\xi')\d\theta\)\d z',
\end{equation*}
where $A:B=\tr(A^\top B) = \sum_{i, j}A_{ij}B_{ij}$ is the Frobenius inner product for square matrices $A$ and $B$.
Then it follows that
\begin{equation*}
\bNrm{\bcom{\n{\xi}^2,\sim}f}{L^p}\le \frac{3\,\hbar}{2} \(\Nrm{f}{L^p} + 2 \Nrm{\n{\xi}\Dv f}{L^p}\) + \frac{16\,\hbar^{\sfrac32}}{\sqrt{\pi}} {\Nrm{\Dv f}{L^p}},
\end{equation*}
which concludes the proof.
\end{proof}

\begin{lem}\label{lem:init_diff_diag}
Let $h\in(0,1)$. Then there exists $C>0$ such that for any function $g$ on the phase space
\begin{equation*}
\bNrm{\weight{\opp}\bigl(\widetilde{\op}_{g}^2 - \widetilde{\op}_{g^2}\bigr) \weight{\opp}}{\L^2} \leq C\,\hbar\, \bigl(\Nrm{g}{W^{1,4}_1\cap H^3}^2 + \Nrm{g^2}{H^1_1\cap H^2}\bigr).
\end{equation*}
\end{lem}

\begin{proof}
Let $\opm := \weight{\opp}$ and $m(\xi) = \weight{\xi}$. By the definition of the absolute value for operators, by the proof of Lemma~\ref{lem:Weyl_and_weight} and by Lemma~\ref{lem:p-weight_remainder_gaussian}, we can write
\begin{equation*}
\opm\,\widetilde{\op}_{g}^2\,\opm = \an{\widetilde{\op}_{g} \opm}^2 = \an{\op_{\tilde{g}} \opm}^2 = \an{\op_{m\,\tilde{g}} + \op_r}^2 = \an{\widetilde{\op}_{m\,g} + \op_{r_1}}^2,
\end{equation*}
where $r_1 = \com{m,\sim}g+r$, and
\begin{equation*}
\opm\,\widetilde{\op}_{g^2}\,\opm = \op_{m^2\widetilde{g^2}} + \op_{r_2} = \widetilde{\op}_{m^2\,g^2} + \op_{r_0},
\end{equation*}
where $r_0 = \com{m^2,\sim}g^2+r_2$. It follows that
\begin{multline}\label{est:difference_toeplitz_square}
\bNrm{\opm\big(\widetilde{\op}_{g}^2 - \widetilde{\op}_{g^2}\big)\opm}{\L^2} = \bNrm{\an{\widetilde{\op}_{m\,g} +\op_{r_1}}^2 - \widetilde{\op}_{m^2\,g^2} - \op_{r_0}}{\L^2}
\\
\leq \bNrm{\widetilde{\op}_{m\,g}^2 - \widetilde{\op}_{m^2\,g^2}}{\L^2} + \aNrm{\op_{r_1}}{\L^4}^2 + 2 \aNrm{\widetilde{\op}_{m\,g}\, \op_{r_1}}{\L^2} + \aNrm{\op_{r_0}}{\L^2}.
\end{multline}
The first term on the right-hand side of Inequality \eqref{est:difference_toeplitz_square} is bounded by Lemma~\ref{lem:commut_square_wick} with $p=2$, leading to
\begin{equation*}
\bNrm{\widetilde{\op}_{m\,g}^2 - \widetilde{\op}_{m^2\,g^2}}{\L^2} \leq 48\,\hbar\, {\Nrm{\nabla(m\,g)}{L^4}^2}.
\end{equation*}
The second term can be bounded by
\begin{equation*}
\aNrm{\op_{r_1}}{\L^4} \leq \aNrm{\op_{r_1} - \widetilde{\op}_{r_1}}{\L^4} +\aNrm{\widetilde{\op}_{r_1}}{\L^4}
\end{equation*}
and then noticing that the Schatten-$4$ norm is controlled by the Hilbert--Schmidt norm, one obtains $\aNrm{\op_{r_1} - \widetilde{\op}_{r_1}}{\L^4} \leq h^{-\sfrac34} \aNrm{\op_{r_1} - \widetilde{\op}_{r_1}}{\L^2} \leq h^{\sfrac14} \aNrm{\nabla^2 r_1}{L^2}$, where we used Inequality~\eqref{eq:Wick_Weyl_error} in the last inequality. Hence we deduce that
\begin{equation*}
\Nrm{\op_{r_1}}{\L^4} \leq \Nrm{r_1}{L^4} + h^{\sfrac14} \aNrm{\nabla^2 r_1}{L^2} \leq C\,\hbar\, \aNrm{g}{H^3}.
\end{equation*}
The third term can be bounded by Hölder's inequality and Inequality~\eqref{eq:Wick_Schatten} as follows
\begin{equation*}
\aNrm{\widetilde{\op}_{m\,g} \, \op_{r_1}}{\L^2} \leq \aNrm{\widetilde{\op}_{m\,g}}{\L^4} \aNrm{\op_{r_1}}{\L^4} \leq C\,\hbar \Nrm{m\,g}{L^4} \aNrm{g}{H^3}.
\end{equation*}
Finally, the last term can be bounded by
\begin{align*}
\aNrm{\op_{r_0}}{\L^2} \le C\, \hbar\,\bigl(\aNrm{g^2}{H^1_1} + \hbar\,\aNrm{\Delta_x g^2}{L^2}\bigr).
\end{align*}
Combining all the inequalities leads to the desired result.
\end{proof}

\begin{lem}\label{lemma:diag-opt-opf}
Let $\op_f$ be a solution to Equation~\eqref{eq:Weyl-Vlasov} and $\widetilde{\op}$ be a solution to Equation~\eqref{eq:linear-Hartree}. Then
\begin{equation*}
\Nrm{\diag(\widetilde{\op}-\op_f)}{L^2_x} \leq C\, e^{C_E}(t)\, \hbar \biggl(C^\init+ \int^t_0 \aNrm{\rho_{f(s,\cdot)}}{W^{1,\infty}_x\cap L^1_x} \bNrm{\op_{f(s,\cdot)}\, \weight{\opp}^2}{\cW^{2,2}}\d s\biggr),
\end{equation*}
where $C^\init$ is given by Equation~\eqref{eq:C_init} and $C_E(t)$ by Equation~\eqref{eq:C-t}.
\end{lem}

\begin{proof}
Let $\opm = \weight{\opp}$. We will use the fact that
\begin{equation*}
\Nrm{\diag(\widetilde{\op}-\op_f)}{L^2_x} \leq C\bNrm{\opm\,(\widetilde{\op}-\op_f)\,\opm}{\L^2}.
\end{equation*}
as follows for example from the proof of \cite[Prop.\,6.4]{chong_many-body_2021}. By the equation verified by $\widetilde{\op}-\op_f$ and the cyclicity of the trace, one gets
\begin{multline*}
i\hbar\,\dt \Nrm{\opm\,(\widetilde{\op}-\op_f)\,\opm}{\L^2}^2 = h^d\Tr{2\n{\opm\(\widetilde{\op}-\op_f\)\opm}^2\opm^{-1} \bcom{\n{\opp}^2,V_f}\opm^{-1}}
\\
+ h^d \bTr{\opm\(\widetilde{\op}-\op_f\)\opm^2\,B_f(\op_f)\,\opm}.
\end{multline*}
By the fact that $\tfrac{1}{i\hbar}\bcom{\n{\opp}^2,V_f} = E_f\cdot\opp + \opp\cdot E_f$ and Hölder's inequality for the trace, we deduce that
\begin{align*}
\dt\Nrm{\opm\,(\widetilde{\op}-\op_f)\,\opm}{\L^2} &\leq 2 \Nrm{\opm\,(\widetilde{\op}-\op_f)\opm}{\L^2}\! \Nrm{\opm^{-1} E_f\cdot\opp\,\opm^{-1}}{\L^\infty}\! + \frac{1}{\hbar}\Nrm{\opm\,B_f(\op_f)\,\opm}{\L^2}
\\
&\leq 2 \Nrm{E_f}{L^\infty} \Nrm{\opm\,(\widetilde{\op}-\op_f)\,\opm}{\L^2} + \frac{1}{\hbar}\,\bNrm{B_f(\op_f)\,\opm^2}{\L^2}.
\end{align*}
By Grönwall's lemma, this implies
\begin{equation*}
\Nrm{\opm\,(\widetilde{\op}-\op_f)\,\opm}{\L^2} \leq e^{C_E(t)}\biggl( \aNrm{\opm\,(\widetilde{\op}^\init-\op_{f^\init})\,\opm}{\L^2} +\frac{1}{\hbar}\int_0^t \Nrm{B_f(\op_f)\,\opm^2}{\L^2}\biggr).
\end{equation*}
with $C_E(t)$ given in \eqref{eq:C-t}. With our choice of the initial datum $\widetilde{\op}^\init = \widetilde{\op}^2_{\sqrt{f^\init}}$, the first term on the right-hand side is bounded in Lemma~\ref{lem:init_diff_diag}. For the second term, we compute explicitly the $L^2$ norm of the kernel of the operator $B_f(\op_f)\,\opm^2$. Recall the definition of $B_f(\op)$ given by Formula~\eqref{def:B_operator} and define
\begin{equation*}
\delta^2 V(x,y) := V_{f}(x) - V_{f}(y) - \(x-y\)\cdot\nabla V_{f}\Big(\frac{x+y}{2}\Big)
\end{equation*}
so that the integral kernel of $B$ can be written $B_f(\op)(x,y) = \delta^2 V(x,y)\,\op(x,y)$. Then by the chain rule, the integral kernel of the operator $B_f(\op_f)\n{\opp}^2$ is given by
\begin{multline*}
\big(B_f(\op_f)\n{\opp}^2\big)(x,y)
= -\hbar^2\, \delta^2V(x,y)\, \lapl_y\op_f(x,y) -\hbar^2 \lapl_y\!\(\delta^2V(x,y)\) \op_f(x,y)
\\
- 2\,\hbar^2\, \grad_y\delta^2V(x, y)\cdot \grad_y\op_f(x,y) =: J_1+J_2+J_3.
\end{multline*}
The $J_1$ term is simply $B_f(\op_f\n{\opp}^2)$ which can be bounded in $\L^2$ by Inequality~\eqref{eq:B_bound}, leading to
\begin{align*}
\bNrm{B_f(\op_f\n{\opp}^2)}{\L^2} \le C \, \hbar^2 \Nrm{\rho_f}{{W^{1,\infty}_x\cap L^1_x}} \bNrm{\Dhv^2\!\(\op_f\n{\opp}^2\)}{\L^2}.
\end{align*}
For the $J_2$ term, since $V$ is the solution of the Poisson equation, we have that
\begin{align*}
-\hbar^2\,\lapl_y \delta^2V(x, y) = \pm\hbar^2 \Big(
\rho_{f}\Big(\frac{x+y}{2}\Big)-\rho_{f}(y)\Big) - \frac{\pm \hbar^2}{4}\, \grad \rho_f\Big(\frac{x+y}{2}\Big)\cdot(x-y),
\end{align*}
which means
\begin{align*}
J_2(x,y) = \pm\hbar^2 \Big(\rho_{f}\Big(\frac{x+y}{2}\Big)-\rho_{f}(y)\Big)\op_{f}(x, y)
\pm\frac{\hbar^3}{4i} \grad \rho_f\Big(\frac{x+y}{2}\Big)\cdot \Dhv \op_{f}(x,y).
\end{align*}
Then we deduce that
\begin{equation*}
\Nrm{J_2}{\L^2} \le 2\,\hbar^2 \Nrm{\rho_f}{L^\infty_x}\Nrm{\op_f}{\L^2}+\frac{\hbar^3}{4} \Nrm{\grad \rho_f}{L^\infty_x}\aNrm{\Dhv\op_f}{\L^2}.
\end{equation*}
Lastly, for the $J_3$ term, observe that
\begin{align*}
-2\,\hbar^2\, \nabla_y \delta^2 V(x, y)= 2\,\hbar^2\Bigl(E_f\Bigl(\frac{x+y}{2} \Bigr)-E_f(y)\Bigr)+\hbar^2 \(x-y\)\cdot \grad E_f\Bigl(\frac{x+y}{2}\Bigr),
\end{align*}
which means
\begin{multline*}
J_3(x,y) = i\hbar^3 \int^1_0 \nabla E_f(y+t\,\psfrac{x-y}{2}) : \nabla_y\Dhv\op_f(x,y) \d t
\\
+ i\hbar^3 \grad E_f\(\psfrac{x+y}{2}\) : \nabla_y\Dhv\op_f(x,y).
\end{multline*}
Hence it follows that
\begin{align*}
\Nrm{J_3}{\L^2} \le C\,\hbar^2 \Nrm{\grad E_f}{L^\infty_x}\aNrm{\Dhv\op_f\, \opp}{\L^2}.
\end{align*}
This completes our proof of the lemma.
\end{proof}

\begin{prop}\label{prop:sqrt-op-vs-opt}
Let $\op$ be a solution to the Hartree equation \eqref{eq:Hartree} with initial datum~$\op^\init$ and $\widetilde{\op}$ be a solution to \eqref{eq:linear-Hartree} with initial datum $\widetilde{\op}^\init$. We denote by $v_1 := \sqrt{\op}$, $\tilde{v} := \sqrt{\widetilde{\op}}$ and $v_1^\init$ and $\tilde{v}^\init$ the corresponding initial data. Then for any $t\in\R_+$,
\begin{equation*}
\Nrm{v_1-\tilde{v}}{\L^2}^2\leq \Nrm{v_1^\init - \tilde{v}^\init}{\L^2}^2 e^{2\,\Lambda(t)} + \hbar^2 \int_0^t c(s)^2 \,e^{2\(\Lambda(t)-\Lambda(s)\)}\d s,
\end{equation*}
where $\Lambda(t) = C \int_0^t \lambda(s)\d s$ with $\lambda$ given by~\eqref{eq:lambda} and $c(t)$ is given by
\begin{equation*}
c(t) = C\,e^{C_E(t)} \Nrm{\Dhv \tilde{v}}{\cW^{1,2}} \int_0^t C^\init + \aNrm{\rho_{f(s,\cdot)}}{W^{1,\infty}_x\cap L^1_x} \bNrm{\op_{f(s,\cdot)} \weight{\opp}^2}{\cW^{2,2}} \d s
\end{equation*}
with $C^\init$ given in Equation~\eqref{eq:C_init} and $C_E(t)$ in Equation~\eqref{eq:C-t}.
\end{prop}

\begin{proof}
Define $v:=v_1-\tilde{v}$. Observe that $v_1$ solves
\begin{equation*}
i\hbar\,\dpt v_1=\acom{H_1,v_1},\quad\quad H_1:=H_{v_1^2}
\end{equation*}
and $\tilde{v}$ solves
\begin{equation*}
i\hbar\,\dpt \tilde{v}=\acom{H_f,\tilde{v}}.
\end{equation*}
Now we proceed by mimicking the proof of Theorem~\ref{thm:Q-uniqueness}. By direct computation, we have that
\begin{equation*}
i\hbar\,\dpt v = \com{H_1,v}+\bcom{V_{\op}-V_{\op_f},\tilde{v}},
\end{equation*}
which implies
\begin{equation}\label{eq:2norm-v}
\frac{1}{4\pi}\,\dt\Nrm{v}{\L^2}^2 = h^2 \bTr{\bcom{V_{v_1^2}-V_{\op_f},\tilde{v}}v}.
\end{equation}
Now adding and subtracting $h^2 \Tr{\com{V_{\tilde{\op}},\tilde{v}}v}$ to the right-hand side of Equation~\eqref{eq:2norm-v} yields
\begin{equation}\label{eq:L2-estimate}
\begin{split}
\frac{1}{4\pi}\,\frac{\d}{\d t}\Nrm{v}{\L^2}^2 &= h^2 \Tr{\com{V_{\op}-V_{\tilde{\op}},\tilde{v}}v} + h^2 \bTr{\bcom{V_{\tilde{\op}}-V_{\op_f},\tilde{v}}v}
\\
&= h^2 \bTr{\bcom{V_{v_1^2-\tilde{v}^2},\tilde{v}}v} + h^2 \bTr{\bcom{V_{\tilde{\op}-\op_f},\tilde{v}}v}
=: J_1+J_2.
\end{split}
\end{equation}
We bound $J_1$ in the same manner as in Theorem~\ref{thm:Q-uniqueness} and get that
\begin{equation*}
J_1 \leq h^2 \Tr{\com{V_{v^2},\tilde{v}}v} + 2\, h^2 \Tr{\com{V_{\tilde{v}v},\tilde{v}}v} =: J_{1,1}+2\,J_{1,2},
\end{equation*}
where the two terms are readily bounded as follows
\begin{align*}
J_{1,1} &\leq C \Nrm{\Dhv\tilde{v}\,\opm}{\L^{3\pm\eps}}\Nrm{v}{\L^\infty}\aNrm{v}{\L^2}^2,
\\
J_{1,2} &\leq C\, { h^{-1}\Nrm{\com{V_{v\tilde{v}},\tilde{v}}}{\L^2}\aNrm{v}{\L^2}}.
\end{align*}
The right-hand side of $J_{1,1}$ is controlled by $\Nrm{v}{\L^2}^2$ since $\Nrm{\Dhv \tilde{v}\,\opm}{\L^{3\pm\eps}}$ is bounded thanks to Proposition~\ref{prop:propagation_regularity} and $\Nrm{v}{\L^\infty} \leq \Nrm{v_1}{\L^\infty} + \Nrm{\tilde{v}}{\L^\infty}$, which are both bounded. By Lemma~\ref{lem:commutator_estimate_HS}, similarly as in the proof of Theorem~\ref{thm:Q-uniqueness}, we get
\begin{equation*}
J_{1,2} \leq C \Nrm{\Dhv \tilde{v}}{\cW^{1,2}} {\Nrm{\rho_{\tilde v}}{L^\infty_x}^{\sfrac12}}\, \aNrm{v}{\L^2}^2,
\end{equation*}
where $\rho_{\tilde v}\in L^\infty(\R^3)$ and $\Nrm{\Dhv \tilde{v}}{\cW^{1,2}}$ is bounded in time thanks to Proposition~\ref{prop:propagation_regularity}. We~are now left to bound the term $J_2$. By Hölder's inequality for Schatten norms and Young's inequality for products we obtain
\begin{equation*}
J_2 = h^2 \Tr{\bcom{V_{\tilde{\op}-\op_f},\tilde{v}}v} \leq \bNrm{\tfrac{1}{h}\bcom{V_{\tilde{\op}-\op_f},\tilde{v}}}{\L^2}^2+ \aNrm{v}{\L^2}^2.
\end{equation*}
We aim at showing that $\tfrac{1}{h}\bNrm{\bcom{V_{\tilde{\op}-\op_f},\tilde{v}}}{\L^2}$ is small. To this end, apply Lemma~\ref{lem:commutator_estimate_HS} to get
\begin{equation}\label{eq:V-opt-opf}
\frac{1}{h}\,\bNrm{\bcom{V_{\tilde{\op}-\op_f},\tilde{v}}}{\L^2} \leq C \Nrm{\Dhv \tilde{v}}{\cW^{1,2}} \bNrm{\diag(\widetilde{\op}-\op_f)}{L^2_x}.
\end{equation}
The term $\aNrm{\diag(\widetilde{\op}-\op_f)}{L^2_x}$ in Equation~\eqref{eq:V-opt-opf} is controlled in Lemma~\ref{lemma:diag-opt-opf}, leading to the statement of the proposition.
\end{proof}

\subsection{Proof of the main theorems}

Now that we have a bound on $\sqrt{\op}-\sqrt{\widetilde{\op}}$ in $\L^2$ by the above proposition, we are ready to finish the proof of the semiclassical limit.
\begin{proof}[Proof of Theorem~\ref{thm:main} and Theorem~\ref{thm:main_2}]
As explained in the beginning of the section, by Lemma~\ref{lem:commut_square_wick} with $p=2$, we get that
\begin{equation}\label{eq:initial-state-f}
\bNrm{\widetilde{\op}_{f^\init} - \widetilde{\op}^\init}{\L^2} \leq 48\,\hbar\,\bNrm{\nabla \sqrt{f^\init}}{L^4}^2,
\end{equation}
and so, together with Inequality~\eqref{eq:optilde-opf}, Inequality~\eqref{eq:Wick_Weyl_error_init} and the fact that $\Nrm{\nabla E_f}{L^\infty_x}$ can be controlled by $\Nrm{\rho_f}{{W^{1,\infty}_x\cap L^1_x}}$, the error due to the lack of positivity can be estimated above by\vspace*{-3pt}
\begin{equation}\label{eq:bound-rho-tilde-rhof}
\Nrm{\widetilde{\op}-\op_f}{\L^2} \leq C\,\hbar\,\biggl({\Nrm{\nabla^2 f^\init}{L^2} +} \bNrm{\nabla \sqrt{f^\init}}{L^4}^2 + \int_0^t \Nrm{\rho_f}{W^{1,\infty}_x} \Nrm{\nabla_{\xi}^2 f}{L^2}\biggr).
\end{equation}
On the other side, the nonlinear stability error~\eqref{eq:sqrt-op_optilde} is controlled, using Proposition~\ref{prop:sqrt-op-vs-opt} and recalling that $\tilde{\op}^\init = \widetilde{\op}_{\sqrt{f^\init}}^2$, by
\begin{equation}\label{eq:bound-eq:bound-rho-rho-tilde}
\Nrm{\op-\widetilde{\op}}{\L^2} \leq \cC_{\infty}^{\sfrac{1}{2}}\(\bNrm{\sqrt{\op^\init} - \widetilde{\op}_{\sqrt{f^\init}}}{\L^2} + \hbar \Nrm{c}{L^2([0,t])}\) e^{\Lambda(t)}.
\end{equation}
Summing up \eqref{eq:bound-rho-tilde-rhof} and \eqref{eq:bound-eq:bound-rho-rho-tilde} and using the estimate on the initial state \eqref{eq:initial-state-f} prove Theorem~\ref{thm:main_2}.

Moreover, we can control the difference of the square roots using Inequality~\eqref{eq:Wick_Weyl_error} which tells us that the Wick and Weyl quantization are close since $\sqrt{f^\init}$ is sufficiently regular, and then the fact that the Wigner transform is an isometry from $\L^2$ to $L^2$ to get\vspace*{-3pt}
\begin{equation*}
\bNrm{\sqrt{\op^\init} - \widetilde{\op}_{\sqrt{f^\init}}}{\L^2} \leq \bNrm{\sqrt{f^\init} - f_{\sqrt{\op^\init}}}{L^2} + \frac{3\,\hbar}{2} \bNrm{\nabla^2 \sqrt{f^\init}}{L^2},
\end{equation*}
which leads to our main inequality~\eqref{eq:Wigner_main_estimate}.
\end{proof}

\subsection{Propagation of regularity}

In this section, we prove the following proposition about the semiclassical propagation of regularity for $\tilde{v}$ assuming sufficient regularity on the solution to the Vlasov--Poisson equation.
\begin{prop}\label{prop:propagation_regularity}
Let $h\in(0,1)$, $q\in[1,\infty)$, $n\in\N$, $f$ be a sufficiently smooth solution of the Vlasov--Poisson equation~\eqref{eq:Vlasov} and $\tilde{v}$ be the solution to the linear equation\vspace*{-3pt}
\begin{equation*}
i \hbar\,\dpt \tilde{v} = \acom{H_f, \tilde{v}},
\end{equation*}
where $H_f = -\tfrac{\hbar^2}{2}\lapl + V_f$. Let $\opm = \weight{\opp}^{2n}$. Then for any $\eps\in(0,1)$, there exists $C_{q,n}$ depending only on $q$, $n$ and $\eps$, such that for every $t>0$ it holds\vspace*{-3pt}
\begin{equation*}
\Nrm{\tilde{v}\, \opm}{\cW^{k, q}} \le \aNrm{\tilde{v}^\init\, \opm}{\cW^{k, q}} \exp \biggl(C_{q,n} \int^t_0 \Nrm{\rho_f(s,\cdot)}{W_x^{2n,3\pm\eps}} \d s\biggr).
\end{equation*}
\end{prop}

\begin{remark}
To obtain the result of Theorem~\ref{thm:main}, we need bounds on the norm of $\tilde{v}$ in $\cW^{2,2}\cap \cW^{1,3\pm\eps}(\weight{\opp}^{2n})$ with $2n>2$. In particular, by the inequalities in~\cite{lafleche_strong_2021}, these norms are bounded for the initial data $\tilde{v}^\init = \tilde{\op}_{\sqrt{f^\init}}$ if
\begin{equation*}
\bNrm{\sqrt{f^\init}}{H^2\cap W^{1,3\pm\eps}_{2n}} +
\hbar^2\, \bNrm{\sqrt{f^\init}}{W^{4,3\pm\eps}}
\end{equation*}
is bounded uniformly in ${\hbar\in(0,1)}$.
\end{remark}

\begin{proof}
We define $\Dhx^k$ and $\Dhv^k$ to be the composition of quantum gradients. From this definition, we see that\vspace*{-3pt}
\[
\Dhx^k \tilde v = \underbrace{[\grad, [\grad, \cdots [\grad, \tilde v]]]}_{k\textup{-times}} \quad \text{ and } \quad \Dhv^k \tilde v = \frac{1}{(i\hbar)^k}\underbrace{[x, [x, \cdots [x, \tilde v]]]}_{k\textup{-times}}
\]
are tensor-valued operators acting on $L^2(\R^3)$.
Here the multiplication of vectors is defined to be their tensor products, that is, for any two vectors $x, y\in \R^3$ we have that $xy= x\otimes y$.

By the Jacobi identity, we arrive at the equation\vspace*{-3pt}
\begin{equation}\label{eq:qgrad_x-equation}
\bd_t \Dhx^k \tilde{v} = \frac{1}{i\hbar} \bcom{H_f, \Dhx^{k} \tilde{v}} -\sum^k_{j=1}\binom{k}{j}\frac{1}{i\hbar}\bcom{\Dhx^{j-1} E_f, \Dhx^{k-j} \tilde{v}}.
\end{equation}
Likewise, using the fact that $\Dhv^j H_f = 0$ for $j\ge 3$, we have that\vspace*{-3pt}
\begin{equation*}
\bd_t \Dhv^k \tilde{v} = \frac{1}{i\hbar} \bcom{H_f, \Dhv^k \tilde{v}} - k\,\Dhv^{k-1}\Dhx \tilde{v}
\end{equation*}
and\vspace*{-3pt}
\begin{multline*}
\bd_t \Dhx^\ell \Dhv^{k-\ell} \tilde{v} = \frac{1}{i\hbar} \bcom{H_f, \Dhx^\ell \Dhv^{k-\ell} \tilde{v}} - \(k-\ell\)\Dhv^{k-\ell-1}\Dhx^{\ell+1} \tilde{v}
\\[-5pt]
-\sum^\ell_{j=1}\binom{\ell}{j}\frac{1}{i\hbar} \bcom{\Dhx^{j-1} E_f, \Dhx^{\ell-j} \Dhv^{k-\ell} \tilde{v}},
\end{multline*}
where $\ell <k$. Let us focus on propagating moments for $\Dhx^k \tilde{v}$ since Equation~\eqref{eq:qgrad_x-equation} only depends on $\Dhx^\ell \tilde{v}$ for $\ell \le k$. Define a weight equivalent to $\opm$ by $\widetilde \opm= 1+\sum^3_{\ii=1}\bp^{2n}_\ii$ as in~\cite[Lem.\,6.3]{chong_many-body_2021}. Then, by~\cite[Lem.\,6.2]{chong_many-body_2021}, we arrive at the estimate\vspace*{-3pt}
\begin{equation*}
\dt\,\bnrm{\Dhx^k \tilde{v}\, \widetilde \opm}_{\cL^q} \le \frac{1}{\hbar}\,\bnrm{\com{V_f, \widetilde \opm} \Dhx^k \tilde{v}}_{\cL^q}+\frac{1}{\hbar}\sum^k_{j=1}\binom{k}{j}\bnrm{\bcom{\Dhx^{j-1}E_f, \Dhx^{k-j} \tilde{v}} \widetilde\opm}_{\cL^q}
\end{equation*}
for $q\ge 2$. Let us estimate the first term. To simplify the notation, we write $V=V_f$, $E=E_f$, and $\opmu = \Dhx^{k} \tilde{v}$. Recall the identity $\com{V, \bp} = i\hbar\, E$. Then it follows\vspace*{-3pt}
\begin{equation*}
\frac{1}{i\hbar} \com{V, \bp^{2n}_{\ii}} = \sum^{2n-1}_{k=0} \bp^k_\ii\, E_\ii\, \bp_\ii^{2n-1-k} =\sum^{2n-1}_{k=0} \sum^k_{\ell = 0} \binom{k}{\ell} g_\ell\, \bp_\ii^{2n-1-\ell},
\end{equation*}
where $g_\ell = \(-i\hbar\)^\ell \bd_\ii^\ell E_\ii$. Using the above identity yields the estimate\vspace*{-3pt}
\[
\frac{1}{\hbar} \Nrm{\com{V, \bp_\ii^{2n}} \opmu}{\cL^q}
\le \sum^{2n-1}_{\ell =0} \binom{2n}{\ell+1} {\Nrm{g_\ell}{L^\infty_x}} \aNrm{\opmu\,\opm}{\cL^q}.
\]
Notice that\vspace*{-3pt}
\begin{equation*}
\Nrm{g_\ell}{L^\infty_x} = \hbar^\ell \Nrm{\grad K * \bd_\ii^{\ell} \rho_f}{L^\infty_x} \leq C \,\hbar^\ell \Nrm{\bd_\ii^{\ell} \rho_f}{L^{3\pm\eps}_x}
\end{equation*}
for $\eps\in(0,1)$, then it follows that\vspace*{-3pt}
\begin{equation*}
\frac{1}{\hbar} \Nrm{\com{V, \bp_\ii^{2n}} \opmu}{\cL^q} \leq C_{f}(t) \Bigl(\frac{(1+\hbar)^{2n}-1}{\hbar}\Bigr) \aNrm{\opmu\,\opm}{\cL^q}.
\end{equation*}
Let us now estimate the second term. Write $\opmu_{k-j} = \Dhx^{k-j} \tilde{v}$. Using the fact that
\[
\com{\Dhx^{j-1}E, \opmu_{k-j}} \bp_\ii^{2n} = \com{\Dhx^{j-1}E, \opmu_{k-j}\, \bp_\ii^{2n}} - \opmu_{k-j}\, \bcom{\Dhx^{j-1}E, \bp_\ii^{2n}},
\]
then it follows\vspace*{-3pt}
\begin{multline*}
\frac{1}{\hbar}\nrm{\com{\Dhx^{j-1}E, \opmu_{k-j}} \bp_\ii^{2n}}_{\cL^q}
\\
\le \frac{1}{\hbar}\nrm{\com{\Dhx^{j-1}E, \opmu_{k-j}\, \bp_\ii^{2n}}}_{\cL^q}+\frac{1}{\hbar}\nrm{\opmu_{k-j}\com{\Dhx^{j-1}E, \bp_\ii^{2n}}}_{\cL^q}=:I_1+I_2.
\end{multline*}
By \cite[Prop.\,6.5]{chong_many-body_2021}, we have that
\begin{equation*}
I_1 \le C\Nrm{\grad^{j-1}\rho_f}{ W^{1,r}_x}\Nrm{\Dhv\opmu_{k-j}\, \bp_\ii^{2n}}{\cL^q}\quad \text{ with } \quad \frac{1}{r}=\frac{1}{2}-\frac{1}{q} < \frac{1}{3}\text{ and $k>j$.}
\end{equation*}
Term $I_2$ is handled in the same manner as in the case of $\tfrac{1}{i\hbar}\com{V, \bp_\ii^{2n}} \opmu$, which yields
\begin{equation*}
I_2 \le C_q \sup_{1\le \ell\le 2n} \nrm{\bd_\ii^{\ell} \rho_f}_{L^{3\pm\eps}_x} \Bigl(\frac{(1+\hbar)^{2n}-1}{\hbar}\Bigr) \aNrm{\opmu\,\opm}{\cL^q}.
\end{equation*}
Hence it follows
\begin{subequations}\label{est:v2_propagation_of_regularity}
\begin{equation*}
\dt\, \bNrm{\Dhx^k \tilde{v}\, \widetilde \opm}{\cL^q} \le C_{q, f, n}(t) \biggl(\bNrm{\Dhx^k \tilde{v}\, \widetilde \opm}{\cL^q} + \sum^k_{j=1}\binom{k}{j} \bNrm{\Dhv\Dhx^{k-j} \tilde{v}\, \widetilde \opm}{\cL^q} \biggr).
\end{equation*}
In fact, recycling the above argument yields the estimates
\begin{equation*}
\dt\, \bNrm{\Dhv^k \tilde{v}\, \widetilde \opm}{\cL^q} \leq\, C_{q, f, n}(t)\( \bnrm{\Dhv^k \tilde{v}\, \widetilde \opm}_{\cL^q}+\bnrm{\Dhv^{k-1} \Dhx \tilde{v}\, \widetilde \opm}_{\cL^q} \)
\end{equation*}
and
\begin{multline*}
\dt\,\bNrm{\Dhx^\ell \Dhv^{k-\ell} \tilde{v}\, \widetilde \opm}{\cL^q}
\le C_{q, f, n}(t) \biggl( \bNrm{\Dhx^\ell \Dhv^{k-\ell} \tilde{v}\, \widetilde \opm}{\cL^q}
\\
+ (k-\ell) \bNrm{\Dhv^{k-1-\ell} \Dhx^{\ell+1} \tilde{v}\, \widetilde \opm}{\cL^q} + \sum^\ell_{j=1} \binom{\ell}{j} \bNrm{\Dhx^{\ell-j} \Dhv^{k-\ell} \tilde{v}\, \widetilde \opm}{\cL^q} \biggr).
\end{multline*}
\end{subequations}
Combining the above estimates~\eqref{est:v2_propagation_of_regularity} gives us a bound of the form
\begin{equation*}
\dt\, \Nrm{\tilde{v}\, \opm}{\cW^{k, q}} \le C_{q, f, n}'(t) \aNrm{\tilde{v}\, \opm}{\cW^{k, q}},
\end{equation*}
which then by Grönwall's lemma leads to the result.
\end{proof}

\backmatter
\bibliographystyle{jepplain+eid}
\bibliography{chong-et-al}
\end{document}