# Decoherence ensures classicality beyond the Ehrenfest time as $\hbar\to 0$

Felipe Hernández, Daniel Ranard, C. Jess Riedel

###### Abstract

In closed quantum systems, wavepackets can spread exponentially in time due to
chaos, forming long-range superpositions in just seconds for ordinary
macroscopic systems. A weakly coupled environment is conjectured to decohere
the system and restore the quantum-classical correspondence while necessarily
introducing diffusive noise—but for what coupling strength, and under what
conditions? For Markovian open systems with Hamiltonians of the form
$\hat{H}=\hat{p}^{2}/2m+V(\hat{x})$ and linear Lindblad operators, we prove
the quantum and classical evolutions are close whenever the strength of the
environment-induced diffusion exceeds a threshold $\hbar^{4/3}A_{c}$, were
$A_{c}$ is a characteristic scale of the classical dynamics. (A companion
paper treats more general Hamiltonians and Lindblad operators.) The bound
applies for all observables and for times exponentially longer than the
Ehrenfest timescale, which is when the correspondence can break down in closed
systems. The strength of the diffusive noise can vanish in the classical limit
to give the appearance of reversible dynamics. The $4/3$ exponent may be
optimal, as Toscano et al. have found evidence that the quantum-classical
correspondence breaks down in some systems when the diffusion is any weaker.

We study the macroscopic emergence of classical mechanics from quantum
mechanics in the limit when $\hbar$ is small compared to the characteristic
features of the system. This is a well-studied topic in the case of _closed_
quantum systems – that is to say, in Hamiltonian systems isolated from any
influence from an external environment – where classical and quantum
observables are known Egorov (1969); Zworski (2022); Combescure and Robert
(1997); Hagedorn and Joye (2000); Silvestrov and Beenakker (2002); Robert and
Combescure (2021) by Egorov’s theorem Egorov (1969) to match closely for times
up to the Ehrenfest time
$\tau_{\mathrm{E}}\sim\lambda_{\mathrm{L}}^{-1}\log(S_{\mathrm{c}}/\hbar)$.
This timescale is governed by the dominant Lyapunov exponent
$\lambda_{\mathrm{L}}$ and the characteristic action scale $S_{\mathrm{c}}$ of
the classical dynamics, and it quantifies the time for a minimal uncertainty
wavepacket to spread significantly due to chaos. The above correspondence
arises essentially because of Ehrenfest’s theorem Ehrenfest (1927), which
implies that a localized wavepacket will approximately follow a classical
equation of motion as long as it remains well-localized.

Figure 1:  (a) An initial pure quantum Gaussian state $\rho(t\mathord{=}0)$
evolves in phase space. (b) At short times the dynamics admit a local harmonic
approximation, squeezing the distribution via classical flow and broadening it
via diffusion $D$. The state, now mixed, can be approximated by a mixture
$\tilde{\rho}(t)$ of pure squeezed Gaussian states (ellipses). (c) Over
intermediate times, the true distribution $\rho(t)$ bends with radius of
curvature $L$ due to anharmonicity, while maintaining minimum thickness
$\sqrt{D/\lambda_{\mathrm{L}}}$, for classical local Lyapunov exponent
$\lambda_{\mathrm{L}}$. Then $\rho(t)$ can still be approximated by
$\tilde{\rho}(t)$ if the distribution $\rho(t)$ can be snugly covered by
ellipses of width $w$ and length $v=\hbar/w$, which geometrically requires
$v\lesssim\sqrt{wL}$ and $w\lesssim\sqrt{D/\lambda_{\mathrm{L}}}$, or
$D\gtrsim\hbar^{4/3}\lambda_{\mathrm{L}}L^{-2/3}$. (d) Over long times the
squeezing of the Gaussian decomposition varies locally but remains bounded as
$\rho(t)$ spreads over the accessible phase space.

For real macroscopic systems, the Ehrenfest time can be quite short—even
seconds or minutes—because the dependence on $\hbar^{-1}$ is only logarithmic
Berman and Zaslavsky (1978); Berry  _et al._ (1979); Chirikov  _et al._
(1988). In closed systems beyond the Ehrenfest time, even as $\hbar\to 0$, the
correspondence between the classical and quantum evolution breaks down. In
particular,

  1. (1)

superpositions over macroscopic distances are generated, detectable through
delicate interference experiments Haake  _et al._ (1987); Zurek and Paz
(1994); Zurek (1998); Habib  _et al._ (1998); Zurek and Paz (1999); Berry
(2001); Silvestrov and Beenakker (2002) (cf. Casati and Chirikov (1995); Zurek
and Paz (1995); Wiebe and Ballentine (2005); Schlosshauer (2008)); relatedly,
the Wigner function develops negativity Berry and Balazs (1979); Berry  _et
al._ (1979); Korsch and Berry (1981); Zurek (2001); Toscano  _et al._ (2006)
(cf. Tomsovic and Heller (1991, 1993)).

  2. (2)

large differences possibly arise between quantum expectation values and the
corresponding classical predictions even for smooth observables like
$\hat{x}^{2}$, according to some numerical studies Karkuszewski  _et al._
(2002); Carvalho  _et al._ (2004) (cf. Emerson (2002)).

In other words, the $\hbar\to 0$ and $t\to\infty$ limits do not commute in
closed systems: if one fixes a time duration $t>0$ and takes $\hbar\to 0$, the
quantum state trajectory approaches the classical state trajectory, but if one
fixes arbitrarily small $\hbar>0$ and takes $t\to\infty$, then the quantum
trajectory may develop superpositions over macroscopic distances. (The limit
is _singular_ in the sense of Berry Berry (1983, 1995, 2001, 2002).)

Despite this theoretical breakdown, macroscopic systems appear to obey the
laws of classical mechanics for much longer times. To theoretically justify a
quantum-classical correspondence beyond the Ehrenfest time, one longstanding
suggestion is to consider decoherence effects from the environment,
acknowledging that macroscopic systems are rarely well-isolated Zurek and Paz
(1994); Shiokawa and Hu (1995); Zurek (1998); Zurek and Paz (1999); Perthame
and Ryzhik (2004) (cf. Ballentine  _et al._ (1994); Fox and Elston (1994a, b);
Schlautmann and Graham (1995); Casati and Chirikov (1995); Zurek and Paz
(1995); Brun  _et al._ (1996); Emerson (2002); Wiebe and Ballentine (2005);
Kofler and Brukner (2007, 2008)). Numerical simulations and analytical
arguments suggest that decoherence successfully restores the quantitative
agreement between the quantum and classical evolution of many specific systems
Spiller and Ralph (1994); Habib  _et al._ (1998); Bhattacharya  _et al._
(2000, 2003); Toscano  _et al._ (2005). In particular, it was conjectured that
there is a regime where the system-bath coupling is (i) large enough for
decoherence to ensure classicality by inhibiting long-range coherence in phase
space, but still (ii) small enough that the noise introduced by the bath does
not significantly alter the classical dynamics on fixed time scales 111More
precisely, the diffusion term $D$ in Eq. (2) may be taken to vanish as
$\hbar\to 0$, and thus for any fixed timescale its effect on the trajectory
vanishes as $\hbar\to 0$. However, for any fixed nonzero $D$, at sufficiently
large time, the effect of $D$ on the classical solution to $\eqref{eq:fp-
simple}$ becomes large. For chaotic evolution, this effect already occurs
after time $t\sim\frac{1}{\lambda_{\mathrm{L}}}\log(D^{-1})$, or equivalently
the Ehrenfest time if $D\sim\hbar$, because the diffusive noise is
exponentially amplified., so the system _appears_ well-isolated Bhattacharya
_et al._ (2000, 2003); Toscano  _et al._ (2005).

In this work we prove several aspects of the quantum-classical correspondence
suggested above using decoherence modeled by the Lindblad equation. The
results presented here apply to Hamiltonians of the common form
$\hat{H}=\hat{p}^{2}/2m+V(\hat{x})$ with linear Lindblad operators and no
friction. Such a system may represent multiple particles in multiple spatial
dimensions, with general position-dependent interactions 222The kinetic term
could be replaced with a general quadratic function of momenta, allowing
particles with distinct masses, but we omit this variation for simplicity..
(With $N$ particles in $n$ spatial dimensions, we have $d=Nn$ degrees of
freedom.)

The case of general Hamiltonians $\hat{H}(\hat{x},\hat{p})$ and Lindblad
operators $\hat{L}_{k}(\hat{x},\hat{p})$ will be presented in a companion
paper. The linear Lindblad operators used here may be seen as the first-order
approximation of a general environmental interaction, and their effect is to
introduce both decoherence and noise. The noise manifests as diffusion in
phase space, or “environment-induced diffusion.”

We prove that given sufficient diffusion, for times $t\ll\hbar^{-1/2},$ the
quantum evolution under the Lindblad equation is well-approximated by the
above classical evolution for all possible observables, in a sense we will
make precise. This timescale is exponentially longer than the Ehrenfest time.
Indeed, for a macroscopic system with characteristic action scale
$S_{\mathrm{c}}\sim 1\,\mathrm{kg}\cdot\mathrm{m}^{2}/\mathrm{s}$ and Lyapunov
exponent $\lambda_{\mathrm{L}}\sim 1\,\mathrm{s}^{-1}$, the Ehrenfest
timescale is
$\tau_{\mathrm{E}}\sim\lambda_{\mathrm{L}}^{-1}\log(S_{\mathrm{c}}/\hbar)\sim
1$ minute, while the bound presented below extends to the timescale
$\tau\sim\lambda_{\mathrm{L}}^{-1}\sqrt{S_{\mathrm{c}}/\hbar}\sim
10^{17}\,\mathrm{s}\sim 1$ billion years.

For a quantum system of finite $d$ degrees of freedom, assuming initial state
$\rho(t\mathord{=}0)$ that can be well-approximated by a mixture of coherent
Gaussian states, we consider the quantum evolution $\rho(t)$ under the
Lindblad equation $\partial_{t}\rho(t)=\hat{\mathcal{L}}[\rho(t)]$ with
Lindbladian

$\displaystyle\hat{\mathcal{L}}[\rho]$
$\displaystyle=\frac{-i}{\hbar}[\hat{H},\rho]+\frac{1}{\hbar}\sum_{k}(\hat{L}_{k}\rho\hat{L}_{k}^{\dagger}-\frac{1}{2}\\{\hat{L}_{k}^{\dagger}\hat{L}_{k},\rho\\}).$
(1)

We treat the case of linear Lindblad operators
$\hat{L}_{k}=\ell_{\mathrm{x}}\hat{x}^{k}$ and
$\hat{L}_{k+d}=\ell_{\mathrm{p}}\hat{p}^{k}$ for $k=1,\ldots,d$. We compare
the quantum evolution $\rho(t)$ to the evolution of a probability distribution
$f(t)$ that obeys a corresponding classical dynamics, making use of the Wigner
transform $\mathcal{W}[\cdot]$ (whose needed properties are recalled in
Appendix B). The initial classical distribution is taken to be the one
naturally associated to the initial quantum state:
$f(t\mathord{=}0)=\mathcal{W}[\rho(t\mathord{=}0)]$. The corresponding
classical dynamics are given by a frictionless Fokker-Planck equation
$\partial_{t}f=\mathcal{L}[f]$ using the Liouvillian Risken (1984); Gardiner
(2009)

$\mathcal{L}[f]=-\omega^{ab}(\partial_{a}f)(\partial_{b}H)+\frac{1}{2}D^{ab}\partial_{a}\partial_{b}f,$
(2)

where the first term above produces classical Hamiltonian flow using
$H=p^{2}/2m+V(x)$, while the second term above produces diffusion with
diffusion matrix

$D=\left(\begin{matrix}|\ell_{\mathrm{p}}|^{2}\mathds{1}_{d}&0\\\
0&|\ell_{\mathrm{x}}|^{2}\mathds{1}_{d}\end{matrix}\right).$ (3)

Above $\omega=\big{(}\begin{smallmatrix}0&\mathds{1}_{d}\\\
-\mathds{1}_{d}&0\end{smallmatrix}\big{)}$ is the (antisymmetric) symplectic
form and $\partial_{a}$ denotes the partial derivative in phase space, where
the indices
$a,b\in\\{\mathrm{x}_{1},\ldots,\mathrm{x}_{d},\mathrm{p}_{1},\ldots,\mathrm{p}_{d}\\}$
range over the $2d$ directions in phase space and repeated indices are summed.

Just as the diffusion equation arises from a random walk, the Fokker-Planck
equation in (2) could also be interpreted as arising from a Langevin equation
Kerr and Graham (2000), describing noisy trajectories on phase space. Thus (2)
describes an ensemble of noisy classical trajectories, each roughly following
the classical equation of motion when $D$ is small.

In order to state the main result, we now identify some key scales in any
Hamiltonian $\hat{H}=\hat{p}^{2}/2m+V(\hat{x})$. At each position $x$ there is
a local harmonic approximation to the dynamics depending on the Hessian
$\nabla^{2}V(x)$ of the potential. The harmonic time
$\tau_{H}:=\sqrt{m/\|\nabla^{2}V\|}$ is the shortest timescale associated with
such local dynamics, where

$\displaystyle\|\nabla^{j}V\|:=\sup_{x}\sup_{\|w_{i}\|=1}\left|\left[\prod_{i=1}^{j}(w_{i}\cdot\nabla)\right]V(x)\right|$
(4)

is the maximum over phase space of $j$ directional derivatives of the
potential. Similarly, the aspect parameter
$\eta_{H}:=\sqrt{m\|\nabla^{2}V\|}=m/\tau_{H}=\|\nabla^{2}V\|\tau_{H}$ (with
units of [momentum/length] $\sim$ kg/s) picks out a preferred “aspect ratio”
in phase space with respect to which a Gaussian distribution can be considered
circular or unsqueezed. The anharmonic action
$s_{H}:=\eta_{H}^{3}/(\tau_{H}\|\nabla^{3}V\|)^{2}=m^{1/2}\|\nabla^{2}V\|^{5/2}\|\nabla^{3}V\|^{-2}$
is the action scale above which anharmonicities of the potential (quantified
by the max _third_ derivative $\|\nabla^{3}V\|$) are important. (It is the
unique action scale constructible from $m$, $\|\nabla^{2}V\|$, and
$\|\nabla^{3}V\|$.) The anharmonic length
$x_{H}:=\sqrt{s_{H}/\eta_{H}}=\|\nabla^{2}V\|/\|\nabla^{3}V\|$ and anharmonic
momentum
$p_{H}:=\sqrt{s_{H}\vphantom{/}\eta_{H}}=\sqrt{m\|\nabla^{2}V\|}/\|\nabla^{3}V\|$
are then preferred scales satisfying $x_{H}p_{H}=s_{H}$. Finally, the pure
Gaussian quantum states with the preferred covariance matrix

$\displaystyle\sigma_{*}:=\frac{\hbar}{2}\left(\begin{matrix}\eta^{-1}_{H}\mathds{1}_{d}&0\\\
0&\eta_{H}\mathds{1}_{d}\end{matrix}\right)\propto\left(\begin{matrix}x_{H}\mathds{1}_{d}&0\\\
0&p_{H}\mathds{1}_{d}\end{matrix}\right)$ (5)

are the (unsqueezed) coherent states.

We can now state our main result.

######  Theorem 1.

Let $\rho(t)$ solve the Lindblad equation (1) in $d$ variables with
Hamiltonian $\hat{p}^{2}/2m+V(\hat{x})$ and $2d$ Lindblad operators
$\hat{L}_{k}=\ell_{\mathrm{x}}\hat{x}^{k}$ and
$\hat{L}_{k+d}=\ell_{\mathrm{p}}\hat{p}^{k}$ for $k=1,\ldots,d$, corresponding
to isotropic diffusion in position and momentum with respective diffusion
constants $D_{(\mathrm{x})}=\hbar|\ell_{\mathrm{p}}|^{2}$ and
$D_{(\mathrm{p})}=\hbar|\ell_{\mathrm{x}}|^{2}$. Assume the initial state
$\rho(t\mathord{=}0)$ is given by a mixture of Gaussian states that are
squeezed relative to the coherent states by no more than a factor
$z:=\max\\{\frac{\hbar/s_{H}}{D_{0}},1\\}\geq 1$ (i.e., with covariance matrix
$\sigma$ s.t. $z^{-1}\sigma_{*}\leq\sigma\leq z\sigma_{*}$) where

$\displaystyle
D_{0}:=\min\left\\{\frac{D_{(\mathrm{x})}}{x_{H}^{2}/\tau_{H}},\frac{D_{(\mathrm{p})}}{p_{H}^{2}/\tau_{H}}\right\\}$
(6)

is a dimensionless measure of the diffusion strength. Finally, let $f(t)$
solve the Fokker-Planck equation (2) with initial condition
$f(0)=\mathcal{W}[\rho(0)]$, and assume the potential $V(x)$ has bounded
second and third derivatives. Then there exists a quantum trajectory
$\tilde{\rho}(t)$ which is a mixture of Gaussians satisfying

$\displaystyle\begin{split}\left\lVert{\widetilde{\rho}}(t)-\rho(t)\right\rVert_{\mathrm{Tr}}&\leq\epsilon\\\
\left\lVert\mathcal{W}[{\widetilde{\rho}}(t)]-f(t)\right\rVert_{L^{1}}&\leq\epsilon\end{split}$
(7)

with error

$\displaystyle\epsilon$
$\displaystyle=d^{\frac{3}{2}}\frac{t}{\tau_{H}}\sqrt{\frac{\hbar}{s_{H}}}\max\left\\{\frac{\hbar/s_{H}}{D_{0}},1\right\\}^{\frac{3}{2}}.$
(8)

Recall $\tau_{H}$, $s_{H}$, $x_{H}$, $p_{H}$ above are just characteristic
scales set by the classical Hamiltonian. In (7),
$\|\hat{A}\|_{\mathrm{Tr}}:=\operatorname{Tr}[(\hat{A}^{\dagger}\hat{A})^{1/2}]$
is the trace norm on quantum operators and
$\|f\|_{L^{1}}:=\int\\!\mathrm{d}\alpha|f(\alpha)|$ is the analogous $L^{1}$
norm on classical phase-space functions; as recalled in Appendix A, they
constrain the probability of a discrepancy being revealed by any observation
or measurement.

We turn to the proof after some interpretation. How much environment-induced
diffusion is necessary for Theorem 1 to ensure a close quantum-classical
correspondence? We find Eq. (7) holds with error $\epsilon$ for any

$\epsilon\geq d^{\frac{3}{2}}\frac{t}{\tau_{H}}\sqrt{\frac{\hbar}{s_{H}}},$
(9)

if the diffusion satisfies

$D_{0}\geq\left(\frac{d^{\frac{3}{2}}t}{\epsilon\,t_{H}}\right)^{\frac{2}{3}}\left(\frac{\hbar}{s_{H}}\right)^{\frac{4}{3}}.$
(10)

Alternatively, if $D_{0}\gtrsim(\hbar/S_{H})^{\frac{4}{3}-p}$ for some power
$p>0$, we find that the error $\epsilon$ in the correspondence (7) is small
for times $t\lesssim t_{H}(s_{H}/\hbar)^{q}$, for power
$q=\min\\{\frac{1}{2},\frac{3p}{2}\\}$. Note this time range is exponentially
longer than the Ehrenfest time.

See Fig. 1 for a heuristic argument suggesting the adequacy of condition (10).
We illustrate how the bound might be applied to physical examples in Appendix
F.

###### Proof of Theorem 1.

We will build a trajectory $\tilde{\rho}(t)$ that approximates the true
evolution $\rho(t)$ by using a mixture of Gaussian states
$\hat{\tau}_{\alpha,\sigma}$, each centered at a point
$\alpha=(\alpha^{\mathrm{x}},\alpha^{\mathrm{p}})\in\mathbb{R}^{2d}$ in phase
space with covariance matrix $\sigma\in\mathbb{R}^{2d}\times\mathbb{R}^{2d}$,
that evolve under a local harmonic approximation to $\hat{\mathcal{L}}$. Thus
our trajectory will be defined through a time-dependent distribution
$p_{\alpha,\sigma}(t)\geq 0$,

$\displaystyle\tilde{\rho}(t)=\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(t)\hat{\tau}_{\alpha,\sigma},$
(11)

but with the crucial proviso that $p_{\alpha,\sigma}(t)$ only supports
covariance matrices $\sigma$ that are not too squeezed (NTS) in the sense
$\sigma\leq z\sigma_{*}$. (The distribution $p_{\alpha,\sigma}$ is a
generalization of the Glauber-Sudarshan P function Glauber (1963); Cahill and
Glauber (1969); Sudarshan (1963); Agarwal and Wolf (1970), which uses a single
fixed $\sigma$.) Importantly, we have assumed the true initial state
$\rho(t\mathord{=}0)=\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(t\mathord{=}0)\hat{\tau}_{\alpha,\sigma},$
is a (possibly trivial) mixture of such states, so that the true trajectory
and our approximation initially coincide:
$\tilde{\rho}(t\mathord{=}0)=\rho(t\mathord{=}0)$.

Completing the proof requires three steps.

_Step 1. Constructing the Gaussian mixture $\tilde{\rho}$:_ If
$\tilde{\rho}(t)$ is to approximate $\rho(t)$, which satisfies the Lindblad
equation (1), then we would like to find $p_{\alpha,\sigma}(t)$ such that
$\partial_{t}\tilde{\rho}(t)\approx\hat{\mathcal{L}}[\tilde{\rho}(t)]=\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(t)\hat{\mathcal{L}}[\hat{\tau}_{\alpha,\sigma}]$.
To accomplish this, we will use a _harmonic_ approximation
$\hat{\mathcal{L}}^{(\alpha)}$ to $\hat{\mathcal{L}}$ near the point $\alpha$
in phase space, where harmonic dynamics are characterized by a quadratic
Hamiltonians and linear Lindblad operators. The benefit of using this
approximation is that harmonic dynamics exactly preserve the set of Gaussian
states Isar  _et al._ (1994); Brodier and Ozorio de Almeida (2004); Graefe
_et al._ (2018) (see also Strunz and Percival (1998); Halliwell and Zoupas
(1995)). We will find a function $p_{\alpha,\sigma}(t)$ that exactly solves

$\partial_{t}\tilde{\rho}(t)=\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(t)\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}].$
(12)

In the present case of linear Lindblad operators and Hamiltonian
$\hat{H}=\hat{p}^{2}/2m+V(\hat{x})$, the harmonic approximation
$\hat{\mathcal{L}}^{(\alpha)}$ to $\hat{\mathcal{L}}$ at $\alpha$ is obtained
simply by replacing the potential with its local quadratic approximation about
$x=\alpha^{\mathrm{x}}$, i.e., by replacing $V(\hat{x})$ with
$V^{[\alpha^{\mathrm{x}},2]}(\hat{x})$ where
$V^{[\alpha^{\mathrm{x}},2]}(x)=V(\alpha^{\mathrm{x}})+(x-\alpha^{\mathrm{x}})\cdot\nabla
V(\alpha^{\mathrm{x}})+[(x-\alpha^{\mathrm{x}})\cdot\nabla]^{2}V(\alpha^{\mathrm{x}})$.
Before we proceed to analyze the error between $\tilde{\rho}(t)$ and
$\rho(t)$, we first show that a function $p_{\alpha,\sigma}(t)$ solving (12)
can be found.

Under the harmonic approximation about $\alpha$, a Gaussian evolves such that
Isar  _et al._ (1994); Brodier and Ozorio de Almeida (2004); Graefe  _et al._
(2018) its centroid $\alpha$ follows the classical (diffusionless) flow on
phase space with flow vector

$\displaystyle U(\alpha):=$
$\displaystyle\left(\begin{matrix}\alpha^{\mathrm{p}}/m\\\ -\nabla
V(\alpha^{\mathrm{x}})\end{matrix}\right)$ (13)

in the sense that $\mathrm{d}\alpha(t)/\mathrm{d}t=U(\alpha)$, while the
covariance matrix $\sigma$ evolves by
$\mathrm{d}\sigma(t)/\mathrm{d}t=S(\alpha,\sigma)$ where

$\displaystyle S(\alpha,\sigma)$ $\displaystyle:=F(\alpha)\sigma+\sigma
F^{\top}(\alpha)+D,$ (14) $\displaystyle F(\alpha)$
$\displaystyle:=\left(\begin{matrix}0&-\mathds{1}_{d}/m\\\
\nabla^{2}V(\alpha^{\mathrm{x}})&0\end{matrix}\right),$ (15)

and where
$D=\operatorname{diag}(D_{(\mathrm{x})}\mathds{1}_{d},D_{(\mathrm{p})}\mathds{1}_{d})$
is the diffusion matrix. This describes the skewing (by $F$) and broadening
(by $D$) of the Gaussian; only the latter increases mixedness. Equivalently,
the Gaussian obeys

$\displaystyle\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]$
$\displaystyle=\left[U(\alpha)\cdot\partial_{\alpha}+S(\alpha,\sigma)\cdot\partial_{\sigma}\right]\hat{\tau}_{\alpha,\sigma}.$
(16)

With (11) and (16), our desired condition (12) becomes

$\displaystyle\begin{split}&\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,\hat{\tau}_{\alpha,\sigma}\frac{\mathrm{d}}{\mathrm{d}t}p_{\alpha,\sigma}(t)\\\
&=\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(t)\left[U(\alpha)\cdot\partial_{\alpha}+S(\alpha,\sigma)\cdot\partial_{\sigma}\right]\hat{\tau}_{\alpha,\sigma}.\end{split}$
(17)

We could integrate the right-hand side of (17) by parts in $\sigma$ and
$\alpha$ to obtain a transport equation for $p_{\alpha,\sigma}$, but we would
quickly lose control of the covariance matrix $\sigma$, which could be
stretched arbitrarily long by the evolution. Instead, we observe that any
component of the flow in the “positive” $\sigma$ direction (which increases
mixedness of the state) can also be re-interpreted as diffusion in the
$\alpha$ direction Weidlich  _et al._ (1967a, b); Haken  _et al._ (1967); Lax
and Louisell (1967); Isar  _et al._ (1991); Diósi and Kiefer (2000).

For any choice of decomposition $S=S_{\mathrm{D}}+S_{0}$ we have
$S\cdot\partial_{\sigma}\hat{\tau}_{\alpha,\sigma}=[S_{0}\cdot\partial_{\sigma}+\frac{1}{2}S_{\mathrm{D}}\cdot\partial_{\alpha}\partial_{\alpha}]\hat{\tau}_{\alpha,\sigma}$
by the Gaussian derivative identity
$\partial_{\sigma}\hat{\tau}_{\alpha,\sigma}=\frac{1}{2}\partial_{\alpha}\partial_{\alpha}\hat{\tau}_{\alpha,\sigma}$
(reviewed in Appendix E.1). Plugging this into (17) and integrating by parts
333We will show below that the support of $p_{\alpha,\sigma}$ is contained
within the set of NTS states, so that there is no boundary term., we see that
(12) is satisfied (with $p_{\alpha,\sigma}$ guaranteed to be non-negative when
$S_{\mathrm{D}}\geq 0$) so long as $p_{\alpha,\sigma}(t)$ solves
$\partial_{t}p_{\alpha,\sigma}(t)=\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.5pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.5pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\mathcal{L}\hss$\crcr}}}\limits}^{(\alpha)}[p_{\alpha,\sigma}(t)]$
with

$\displaystyle\begin{split}\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.5pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.5pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\mathcal{L}\hss$\crcr}}}\limits}^{(\alpha)}[p_{\alpha,\sigma}]:=&\,\Big{[}-\partial_{\alpha}\cdot
U(\alpha)-\partial_{\sigma}\cdot S_{0}(\alpha,\sigma)\\\
&\qquad+\frac{1}{2}\partial_{\alpha}\partial_{\alpha}\cdot
S_{\mathrm{D}}(\alpha,\sigma)\Big{]}p_{\alpha,\sigma}\end{split}$ (18)

where the partial derivatives in (18) are understood to act also on
$p_{\alpha,\sigma}$.

We have some limited freedom in choosing $S_{0}$ and $S_{\mathrm{D}}$; that
is, the decomposition of $S$ into skewing and broadening parts is not unique.
We require (i) $S_{\mathrm{D}}\geq 0$ so that $p_{\alpha,\sigma}$ undergoes
non-negative diffusion, and we must choose $S_{0}$ so that both (ii)
$p_{\alpha,\sigma}$ remains supported on covariance matrices of pure Gaussian
states and (iii) $p_{\alpha,\sigma}$ remains supported on NTS covariance
matrices. Below, we choose a decomposition fulfilling these three
requirements.

It is useful to work with the rescaled matrices
${\bar{X}}:=\sigma_{*}^{-1/2}X\sigma_{*}^{-1/2}$ for
$X=\sigma,S,S_{0},S_{\mathrm{D}},D$. We also use the asymmetric
${\bar{F}}:=\sigma_{*}^{-1/2}F\sigma_{*}^{1/2}$ (note $\pm\,1/2$ exponents) so
that
${\bar{S}}={\bar{F}}{\bar{\sigma}}+{\bar{\sigma}}{\bar{F}}^{\top}+{\bar{D}}$.
Then we make the choice

$\displaystyle{\bar{S}}_{0}(\alpha,{\sigma})$
$\displaystyle:=[{\bar{F}}(\alpha)-g({\bar{\sigma}})]{\bar{\sigma}}+{\bar{\sigma}}[{\bar{F}}^{\top}(\alpha)-g({\bar{\sigma}})],$
(19) $\displaystyle{\bar{S}}_{\mathrm{D}}(\alpha,{\sigma})$
$\displaystyle:={\bar{D}}+[g({\bar{\sigma}}){\bar{\sigma}}+{\bar{\sigma}}g({\bar{\sigma}})],$
(20)

satisfying ${\bar{S}}_{0}+{\bar{S}}_{\mathrm{D}}={\bar{S}}$, where

$\begin{split}g({\bar{\sigma}})&=\left(\frac{D_{0}s_{H}}{\hbar\tau_{H}}\right)\frac{{\bar{\sigma}}-{\bar{\sigma}}^{-1}}{1-z^{-2}}.\end{split}$
(21)

First, a bit of algebra shows that ${\bar{S}}_{0}={\bar{S}}_{0}{}^{\top}$ and
$({\bar{\sigma}}^{-1}{\bar{S}}_{0})^{\top}=-\omega^{\top}{\bar{\sigma}}^{-1}{\bar{S}}_{0}\omega$
(because $\omega^{\top}{\bar{\sigma}}\omega={\bar{\sigma}}^{-1}$). As recalled
in Appendix C, this ensures that a covariance matrix evolving by $S_{0}$
remains a legal covariance matrix for a pure Gaussian state under the dynamics
(18).

Next, we need to show that with these choices the distribution
$p_{\alpha,\sigma}(t)$ never develops support on covariance matrices violating
the NTS condition $\sigma\leq z\sigma_{*}$, i.e., that
${\bar{\sigma}}=\sigma_{*}^{-1/2}\sigma\sigma_{*}^{-1/2}\leq z\mathds{1}_{d}$
for any $\sigma$ such $p_{\alpha,\sigma}\neq 0$. Note that an equivalent NTS
condition is $\sigma\geq z^{-1}\sigma_{*}$ because $\sigma/(\hbar/2)$ and
$\sigma_{*}/(\hbar/2)$ are symplectic matrices. (See Appendix C, Eq. (56) for
an elementary demonstration.) This equivalent lower-bound condition will be
preserved if $\langle v|{\bar{S}}_{0}(\alpha,\sigma)|v\rangle\geq 0$ whenever
$v$ is an eigenvector of ${\bar{\sigma}}$ with eigenvalue $\lambda\leq z^{-1}$
(${\bar{\sigma}}v=\lambda v$). From (20) we compute

$\begin{split}\langle
v|{\bar{S}}_{0}(\alpha,\sigma)|v\rangle&\geq-2\|{\bar{F}}(\alpha)\|_{\mathrm{op}}z^{-1}+2D_{0}s_{H}\tau_{H}^{-1}\geq
0\end{split}$ (22)

because $z:=\max\\{\hbar/s_{H}D_{0},1\\}\geq\hbar/s_{H}D_{0}$ and
$\|{\bar{F}}(\alpha)\|_{\mathrm{op}}\leq\sqrt{\|\nabla^{2}V\|/m}=\tau_{H}^{-1}$
by (15). (Here, $\|\,\cdot\,\|_{\mathrm{op}}$ denotes the operator norm, i.e.,
the largest singular value of a matrix.)

This ensures that $p_{\alpha,\sigma}(t)$ is only supported on NTS covariance
matrices ($z^{-1}\sigma_{*}\leq\sigma\leq z\sigma_{*}$) for all time. Then
using
$\|{\bar{D}}\|_{\mathrm{op}}=\|\sigma_{*}^{-1/2}D\sigma_{*}^{-1/2}\|_{\mathrm{op}}=\min\\{D_{(\mathrm{x})}\eta_{H},D_{(\mathrm{p})}/\eta_{H}\\}/(\hbar/2)=2D_{0}s_{H}/\hbar\tau_{H}$
we have by (20) that
$S_{\mathrm{D}}(\alpha,\sigma)=\sigma_{*}^{1/2}{\bar{S}}_{\mathrm{D}}(\alpha,\sigma)\sigma_{*}^{1/2}\geq
0$ for all allowed $\sigma$. Therefore the diffusion term in (18) ensures that
$p_{\alpha,\sigma}(t)\geq 0$ for all $t\geq 0$, i.e., $\tilde{\rho}(t)$ is
always a true mixture of squeezed Gaussians.

To summarize, we have constructed the trajectory $\tilde{\rho}(t)$ defined
through (12) with the probability distribution $p_{\alpha,\sigma}(t)$ over
covariance matrices satisfying the NTS condition $\sigma\leq z\sigma_{*}$
defined through the dynamics
$\partial_{t}p_{\alpha,\sigma}(t)=\mathchoice{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\displaystyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{5.00003pt}{3.8889pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\textstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{3.5pt}{2.72223pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptstyle\mathcal{L}\hss$\crcr}}}\limits}{\mathop{\vbox{\halign{#\cr\kern
0.80002pt$\hss\leavevmode\resizebox{2.5pt}{1.94444pt}{\rotatebox[origin={c}]{90.0}{(}}\hss$\crcr\nointerlineskip\cr$\hss\scriptscriptstyle\mathcal{L}\hss$\crcr}}}\limits}^{(\alpha)}[p_{\alpha,\sigma}(t)]$
of (18) using choices (19) and (20).

_Step 2. Duhamel bound:_ Having defined our trajectory $\tilde{\rho}$, we can
compare it to the exact evolution $\rho(t)$ using the Duhamel formula:

$\displaystyle\tilde{\rho}(t)-\rho(t)$
$\displaystyle=\int_{0}^{t}\\!\mathrm{d}s\,e^{(t-s)\hat{\mathcal{L}}}\left(\partial_{s}-\hat{\mathcal{L}}\right)[\tilde{\rho}(s)]$
(23)

Then

$\displaystyle\left\lVert\tilde{\rho}(t)-\rho(t)\right\rVert_{\mathrm{Tr}}$
(24)
$\displaystyle=\left\lVert\int_{0}^{t}\\!\mathrm{d}se^{(t-s)\hat{\mathcal{L}}}\left(\partial_{s}-\hat{\mathcal{L}}\right)[\tilde{\rho}(s)]\right\rVert_{\mathrm{Tr}}$
(25) $\displaystyle\leq\int_{0}^{t}\\!\mathrm{d}s\left\lVert
e^{(t-s)\hat{\mathcal{L}}}\left(\partial_{s}-\hat{\mathcal{L}}\right)[\tilde{\rho}(s)]\right\rVert_{\mathrm{Tr}}$
(26)
$\displaystyle\leq\int_{0}^{t}\\!\mathrm{d}s\left\lVert\left(\partial_{s}-\hat{\mathcal{L}}\right)[\tilde{\rho}(s)]\right\rVert_{\mathrm{Tr}}$
(27)
$\displaystyle=\int_{0}^{t}\\!\mathrm{d}s\left\lVert\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(s)\left(\hat{\mathcal{L}}^{(\alpha)}-\hat{\mathcal{L}}\right)[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}$
(28)
$\displaystyle\leq\int_{0}^{t}\\!\mathrm{d}s\iint\\!\mathrm{d}\alpha\mathrm{d}\sigma\,p_{\alpha,\sigma}(s)\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}$
(29) $\displaystyle\begin{split}&\leq\sup_{\sigma\leq
z\sigma_{*}}\sup_{\alpha}\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}\int_{0}^{t}\\!\mathrm{d}s\end{split}$
(30) $\displaystyle=t\sup_{\sigma\leq
z\sigma_{*}}\sup_{\alpha}\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}},$
(31)

where (27) follows from the fact that $e^{(t-s)\hat{\mathcal{L}}}$ is a
completely positive map and so cannot increase the trace norm, (28) follows
from (12) and the normalization of $p_{\alpha}(t)$, and in (29) we have
defined
$\delta\hat{\mathcal{L}}^{(\alpha)}:=\hat{\mathcal{L}}-\hat{\mathcal{L}}^{(\alpha)}$.

If we make the replacements $\tilde{\rho}\to\mathcal{W}[\tilde{\rho}]$,
$\rho\to f$, $\hat{\tau}_{\alpha,\sigma}\to\tau_{\alpha,\sigma}$,
$\hat{\mathcal{L}}\to\mathcal{L}$, and
$\left\lVert\,\cdot\,\right\rVert_{\mathrm{Tr}}\to\left\lVert\,\cdot\,\right\rVert_{L^{1}}$
in (23–31), then identical manipulations give

$\displaystyle\left\lVert\mathcal{W}[\tilde{\rho}(t)]-f(t)\right\rVert_{L^{1}}\leq
t\sup_{\sigma\leq
z\sigma_{*}}\sup_{\alpha}\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}.$
(32)

Here $\delta\mathcal{L}^{(\alpha)}:=\mathcal{L}-\mathcal{L}^{(\alpha)}$ where,
analogously to $\hat{\mathcal{L}}^{(\alpha)}$, $\mathcal{L}^{(\alpha)}$ is the
harmonic approximation to the classical Liovillian obtained by replacing $V\to
V^{[\alpha^{\mathrm{x}},2]}$ in (2).

_Step 3. Harmonic approximation error:_ The error from the harmonic
approximation will be small when the spatial size of the Gaussian state is
small compared to the anharmonicity of the potential $V$. Indeed, the only
modification necessary to obtain the harmonic approximation to the dynamics
involves replacing the potential with its second-order Taylor approximation,
$V\to V^{[\alpha^{\mathrm{x}},2]}$, so the error is proportional to the max
leading-order (i.e., third-order) correction $\Delta x^{3}\|\nabla^{3}V\|$ at
the characteristic spatial width $\Delta
x\sim\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{\frac{1}{2}}$ of the
wavepacket with covariance matrix $\sigma$. (The $\mathrm{x}\mathrm{x}$
superscript is used to denote the upper left block of $\sigma$.) More
precisely, in Appendix D we prove the bounds

$\displaystyle\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}$
$\displaystyle\leq\mu,\quad\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}\leq\mu$
(33)

where

$\displaystyle\mu=\sqrt{3}\,d^{\frac{3}{2}}\hbar^{-1}\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{\frac{3}{2}}\|\nabla^{3}V\|.$
(34)

In the quantum case, the $\hbar^{-1}$ factor above can be seen to arise from
the $\hbar^{-1}$ in the Schrödinger equation; in the classical case, it is
related to the fact that the coherent state has area $\hbar$ in phase space.

Because $\tilde{\rho}(t)$ is a mixture of Gaussian states that are NTS, the
harmonic errors appearing in the quantum error (31) and classical error (32)
are taken only over covariance matrices satisfying $\sigma\leq z\sigma_{*}$,
so $\|\sigma^{\mathrm{x}\mathrm{x}}\|\leq
z\|\sigma_{*}^{\mathrm{x}\mathrm{x}}\|=z\hbar/2\eta_{H}$. Thus applying (33)
we get that quantum and classical errors are both upper bounded by

$\displaystyle\epsilon$
$\displaystyle=t\mu\leq\sqrt{3}\,t\,d^{\frac{3}{2}}\hbar^{-1}(z\hbar/2\eta_{H})^{\frac{3}{2}}\|\nabla^{3}V\|$
(35)

giving (8). ∎

###  .1 Acknowledgements

We thank Yonah Borns-Weil for valuable discussion and collaboration on related
efforts, Eugene Tang for detailed discussion of the manuscript, and Wojciech
Zurek for foundational insights that inspired this work. We thank Joseph
Emerson, Fabricio Toscano, and Diego Wisniacki for helpful discussion of
numerical evidence for the quantum-classical correspondence in closed and open
systems.

## References

  * Egorov (1969) Y. V. Egorov, Uspekhi Matematicheskikh Nauk 24, 235 (1969).
  * Zworski (2022) M. Zworski,  _Semiclassical analysis_ , Vol. 138 (American Mathematical Society, 2022).
  * Combescure and Robert (1997) M. Combescure and D. Robert, Asymptotic Analysis 14, 377 (1997).
  * Hagedorn and Joye (2000) G. Hagedorn and A. Joye, Annales Henri Poincaré 1, 837 (2000).
  * Silvestrov and Beenakker (2002) P. G. Silvestrov and C. W. J. Beenakker, Physical Review E 65, 035208 (2002).
  * Robert and Combescure (2021) D. Robert and M. Combescure,  _Coherent States and Applications in Mathematical Physics_ , 2nd ed. (Springer, Cham, Switzerland, 2021).
  * Ehrenfest (1927) P. Ehrenfest, Zeitschrift für Physik 45, 455 (1927).
  * Berman and Zaslavsky (1978) G. P. Berman and G. M. Zaslavsky, Physica A: Statistical Mechanics and its Applications 91, 450 (1978).
  * Berry  _et al._ (1979) M. V. Berry, N. L. Balazs, M. Tabor, and A. Voros, Annals of Physics 122, 26 (1979).
  * Chirikov  _et al._ (1988) B. Chirikov, F. Izrailev, and D. Shepelyansky, Physica D: Nonlinear Phenomena 33, 77 (1988).
  * Haake  _et al._ (1987) F. Haake, M. Kuś, and R. Scharf, Zeitschrift für Physik B Condensed Matter 65, 381 (1987).
  * Zurek and Paz (1994) W. H. Zurek and J. P. Paz, Physical Review Letters 72, 2508 (1994).
  * Zurek (1998) W. H. Zurek, Physica Scripta T76, 186 (1998).
  * Habib  _et al._ (1998) S. Habib, K. Shizume, and W. H. Zurek, Physical Review Letters 80, 4361 (1998).
  * Zurek and Paz (1999) W. H. Zurek and J. P. Paz, in  _Epistemological and Experimental Perspectives on Quantum Physics_, Vienna Circle Institute Yearbook [1999] No. 7, edited by D. Greenberger, W. L. Reiter, and A. Zeilinger (Springer Netherlands, 1999) pp. 167–177.
  * Berry (2001) M. V. Berry, in  _Vatican Observatory CTNS Publications_ (2001).
  * Casati and Chirikov (1995) G. Casati and B. V. Chirikov, Physical Review Letters 75, 350 (1995).
  * Zurek and Paz (1995) W. H. Zurek and J. P. Paz, Physical Review Letters 75, 351 (1995).
  * Wiebe and Ballentine (2005) N. Wiebe and L. E. Ballentine, Physical Review A 72, 022109 (2005).
  * Schlosshauer (2008) M. Schlosshauer, Foundations of Physics 38, 796 (2008).
  * Berry and Balazs (1979) M. V. Berry and N. L. Balazs, Journal of Physics A: Mathematical and General 12, 625 (1979).
  * Korsch and Berry (1981) H. J. Korsch and M. V. Berry, Physica D: Nonlinear Phenomena 3, 627 (1981).
  * Zurek (2001) W. H. Zurek, Nature 412, 712 (2001).
  * Toscano  _et al._ (2006) F. Toscano, D. A. R. Dalvit, L. Davidovich, and W. H. Zurek, Physical Review A 73, 023803 (2006).
  * Tomsovic and Heller (1991) S. Tomsovic and E. J. Heller, Physical Review Letters 67, 664 (1991).
  * Tomsovic and Heller (1993) S. Tomsovic and E. J. Heller, Physical Review E 47, 282 (1993).
  * Karkuszewski  _et al._ (2002) Z. P. Karkuszewski, J. Zakrzewski, and W. H. Zurek, Physical Review A 65, 042113 (2002).
  * Carvalho  _et al._ (2004) A. Carvalho, R. de Matos Filho, and L. Davidovich, Physical Review E 70, 026211 (2004).
  * Emerson (2002) J. Emerson, arXiv:quant-ph/0211035 (2002), arxiv:quant-ph/0211035 .
  * Berry (1983) M. V. Berry, Les Houches lecture series 36, 171 (1983).
  * Berry (1995) M. Berry, in  _Studies in Logic and the Foundations of Mathematics_ , Logic, Methodology and Philosophy of Science IXProceedings of the Ninth International Congress of Logic, Methodology and Philosophy of Science, Vol. 134, edited by B. S. a. D. W. Dag Prawitz (Elsevier, 1995) pp. 597–607.
  * Berry (2002) M. Berry, Physics Today 55, 10 (2002).
  * Shiokawa and Hu (1995) K. Shiokawa and B. L. Hu, Physical Review E 52, 2497 (1995), chao-dyn/9503009 .
  * Perthame and Ryzhik (2004) B. Perthame and L. Ryzhik, Methods and Applications of Analysis 11, 447 (2004).
  * Ballentine  _et al._ (1994) L. E. Ballentine, Y. Yang, and J. P. Zibin, Physical Review A 50, 2854 (1994).
  * Fox and Elston (1994a) R. F. Fox and T. C. Elston, Physical Review E 49, 3683 (1994a).
  * Fox and Elston (1994b) R. F. Fox and T. C. Elston, Physical Review E 50, 2553 (1994b).
  * Schlautmann and Graham (1995) M. Schlautmann and R. Graham, Physical Review E 52, 340 (1995).
  * Brun  _et al._ (1996) T. A. Brun, I. C. Percival, and R. Schack, Journal of Physics A: Mathematical and General 29, 2077 (1996), arxiv:quant-ph/9509015 .
  * Kofler and Brukner (2007) J. Kofler and Č. Brukner, Physical Review Letters 99, 180403 (2007).
  * Kofler and Brukner (2008) J. Kofler and Č. Brukner, Physical Review Letters 101, 090403 (2008).
  * Spiller and Ralph (1994) T. P. Spiller and J. F. Ralph, Physics Letters A 194, 235 (1994).
  * Bhattacharya  _et al._ (2000) T. Bhattacharya, S. Habib, and K. Jacobs, Physical Review Letters 85, 4852 (2000).
  * Bhattacharya  _et al._ (2003) T. Bhattacharya, S. Habib, and K. Jacobs, Physical Review A 67, 042103 (2003).
  * Toscano  _et al._ (2005) F. Toscano, R. de Matos Filho, and L. Davidovich, Physical Review A 71, 010101 (2005).
  * Note (1) More precisely, the diffusion term $D$ in Eq. (2) may be taken to vanish as $\hbar\to 0$, and thus for any fixed timescale its effect on the trajectory vanishes as $\hbar\to 0$. However, for any fixed nonzero $D$, at sufficiently large time, the effect of $D$ on the classical solution to $\eqref{eq:fp-simple}$ becomes large. For chaotic evolution, this effect already occurs after time $t\sim\frac{1}{\lambda_{\mathrm{L}}}\log(D^{-1})$, or equivalently the Ehrenfest time if $D\sim\hbar$, because the diffusive noise is exponentially amplified.
  * Note (2) The kinetic term could be replaced with a general quadratic function of momenta, allowing particles with distinct masses, but we omit this variation for simplicity.
  * Risken (1984) H. Risken,  _The Fokker-Planck Equation_, Springer Series in Synergetics, Vol. 18 (Springer, Berlin, Heidelberg, 1984).
  * Gardiner (2009) C. W. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , 4th ed. (Springer, Berlin, 2009).
  * Kerr and Graham (2000) W. C. Kerr and A. Graham, The European Physical Journal B-Condensed Matter and Complex Systems 15, 305 (2000).
  * Glauber (1963) R. J. Glauber, Physical Review 131, 2766 (1963).
  * Cahill and Glauber (1969) K. E. Cahill and R. J. Glauber, Physical Review 177, 1882 (1969).
  * Sudarshan (1963) E. C. G. Sudarshan, Physical Review Letters 10, 277 (1963).
  * Agarwal and Wolf (1970) G. S. Agarwal and E. Wolf, Physical Review D 2, 2161 (1970).
  * Isar  _et al._ (1994) A. Isar, A. Sandulescu, H. Scutaru, E. Stefanescu, and W. Scheid, International Journal of Modern Physics E 3, 635 (1994).
  * Brodier and Ozorio de Almeida (2004) O. Brodier and A. M. Ozorio de Almeida, Physical Review E 69, 016204 (2004).
  * Graefe  _et al._ (2018) E. Graefe, B. Longstaff, T. Plastow, and R. Schubert, Journal of Physics A: Mathematical and Theoretical 51, 365203 (2018).
  * Strunz and Percival (1998) W. T. Strunz and I. C. Percival, Journal of Physics A: Mathematical and General 31, 1801 (1998).
  * Halliwell and Zoupas (1995) J. Halliwell and A. Zoupas, Physical Review D 52, 7294 (1995).
  * Weidlich  _et al._ (1967a) W. Weidlich, H. Risken, and H. Haken, Zeitschrift für Physik 201, 396 (1967a).
  * Weidlich  _et al._ (1967b) W. Weidlich, H. Risken, and H. Haken, Zeitschrift für Physik A Hadrons and nuclei 204, 223 (1967b).
  * Haken  _et al._ (1967) H. Haken, H. Risken, and W. Weidlich, Zeitschrift für Physik 206, 355 (1967).
  * Lax and Louisell (1967) M. Lax and W. Louisell, IEEE Journal of Quantum Electronics 3, 47 (1967).
  * Isar  _et al._ (1991) A. Isar, W. Scheid, and A. Sandulescu, Journal of Mathematical Physics 32, 2128 (1991).
  * Diósi and Kiefer (2000) L. Diósi and C. Kiefer, Physical Review Letters 85, 3552 (2000).
  * Note (3) We will show below that the support of $p_{\alpha,\sigma}$ is contained within the set of NTS states, so that there is no boundary term.
  * Curtright  _et al._ (2014) T. L. Curtright, D. B. Fairlie, and C. K. Zachos,  _A Concise Treatise on Quantum Mechanics in Phase Space_ (WORLD SCIENTIFIC, 2014).
  * Simon  _et al._ (1994) R. Simon, N. Mukunda, and B. Dutta, Physical Review A 49, 1567 (1994).
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Reviews of Modern Physics 84, 621 (2012).
  * Petersen and Pedersen (2012) K. B. Petersen and M. S. Pedersen, “The Matrix Cookbook,” Internet Archive (2012), uRL: https://archive.org/details/K_B_Petersen_and_M_S_Peders__The_Matrix_Cookbook (version: 2013-10-5).
  * Note (4) Some sources will express this for a not-necessarily invertible matrix $Z$ using the matrix adjugate $\mathrm{adj}(Z)$. A property of the adjugate is that $Z\,\mathrm{adj}(Z)=(\det Z)\mathds{1}$ so that, when $Z$ is invertible, $\mathrm{adj}(Z)=(\det Z)Z^{-1}$.
  * Joos  _et al._ (2013) E. Joos, H. D. Zeh, C. Kiefer, D. J. Giulini, J. Kupsch, and I.-O. Stamatescu,  _Decoherence and the appearance of a classical world in quantum theory_ (Springer Science & Business Media, 2013).

##  Appendix A Norms

For a function $f(\alpha)$ over phase space variable
$\alpha=(x,p)\in\mathbb{R}^{2d}$, the (Lebesgue) $L^{q}$ norm is

$\displaystyle\left\lVert
f\right\rVert_{L^{q}}:=\left(\int\\!\mathrm{d}\alpha\left|f(\alpha)\right|^{q}\right)^{1/q}.$
(36)

Our classical error bound on the difference in two distributions $f$ and $g$
is stated with $L^{1}$ norm: $\left\lVert
f-g\right\rVert_{L^{1}}:=\int\\!\mathrm{d}\alpha|f(\alpha)-g(\alpha)|$. The
analogous norm on the quantum side is the trace norm
$\|\hat{A}\|_{\mathrm{Tr}}=\operatorname{Tr}[(\hat{A}^{\dagger}\hat{A})^{1/2}]=\sum_{i}\sigma_{i}(\hat{A})$
of an operator $\hat{A}$, i.e., the sum of the singular values
$\sigma_{i}(\hat{A})$.

The importance of these two norms follows from their standard variational
characterization. For quantum states $\rho$ and $\tilde{\rho}$, the respective
probabilities $q_{\mathrm{q}}=\operatorname{Tr}[\hat{Q}\rho(t)]$ and
$\tilde{q}_{\mathrm{q}}=\operatorname{Tr}[\hat{Q}\tilde{\rho}(t)]$ for the
measurement outcome associated with a projector $\hat{Q}\leq\hat{I}$ are
constrained by

$\displaystyle|q_{\mathrm{q}}-\tilde{q}_{\mathrm{q}}|=|\operatorname{Tr}[Q(\rho(t)-\tilde{\rho}(t))]|\leq\|\hat{Q}\|_{\mathrm{op}}\left\lVert\rho(t)-\tilde{\rho}(t)\right\rVert_{\mathrm{Tr}}\leq\left\lVert\rho(t)-\tilde{\rho}(t)\right\rVert_{\mathrm{Tr}}$
(37)

Likewise, for classical states $f$ and $\mathcal{W}[\tilde{\rho}(t)]$ the
respecive probabilities
$q_{\mathrm{c}}=\int\\!\mathrm{d}\alpha\,Q(\alpha)f(t)$ and
$\tilde{q}_{\mathrm{c}}=\int\\!\mathrm{d}\alpha\,Q(\alpha)\mathcal{W}[\tilde{\rho}(t)]$
for the measurement of any classical indicator variable $Q(\alpha)\leq 1$ are
constrained by

$\displaystyle|q_{\mathrm{c}}-\tilde{q}_{\mathrm{c}}|=\left|\int\\!\mathrm{d}\alpha\,Q(\alpha)(f(t)-\mathcal{W}[\tilde{\rho}(t)])\right|\leq\left[\sup_{\alpha}|Q(\alpha)|\right]\left\lVert\rho(t)-\tilde{\rho}(t)\right\rVert_{L^{1}}\leq\left\lVert\rho(t)-\tilde{\rho}(t)\right\rVert_{L^{1}}.$
(38)

Thus, two classical states cannot be readily distinguished when they are close
in $L^{1}$ norm, and two quantum states cannot be readily distinguished when
they are close in trace norm, no matter what measurement is performed.

##  Appendix B Wigner functions

In order to construct our classical approximation to the quantum trajectory,
we use the Wigner function. The Wigner function $\mathcal{W}[\rho]$ of a
density matrix $\rho$ is a real-valued function on phase space. In this paper,
we will only discuss the Wigner function of mixtures of Gaussian states. Hence
we are using $\mathcal{W}[\cdot]$ essentially as a compact notation for
discussing Gaussian states and Gaussian distributions, without the Wigner
functions playing a key role in the argument. In particular, we only need to
know three facts about the Wigner function:

  1. 1.

When $\rho$ is a Gaussian state, $\mathcal{W}[\rho]$ is a Gaussian probability
distribution on phase space with the same respective mean and variance:
$\mathcal{W}[\hat{\tau}_{\alpha,\sigma}]=\tau_{\alpha,\sigma}$.

  2. 2.

When $\rho$ is a mixture of Gaussian states, $\mathcal{W}[\rho]$ is a mixture
of the respective Gaussian distributions, because $\mathcal{W}[\cdot]$ is a
linear mapping.

  3. 3.

The Weyl trace formula in the special case of the expectation value of a
function $V$ of position for a Gaussian state

$\displaystyle\operatorname{Tr}[\hat{\tau}_{\alpha,\sigma}V(\hat{x})]=\int\\!\mathrm{d}\alpha\,\tau_{\alpha,\sigma}(\alpha)V(\alpha^{\mathrm{x}}).$
(39)

For completeness, we note that for a more general density matrix $\rho$, the
Wigner function is defined as

$\displaystyle\mathcal{W}[\rho](\alpha):=$
$\displaystyle\,(2\pi\hbar)^{d}\int_{\mathbb{R}^{2d}}\\!\mathrm{d}\chi\,\operatorname{Tr}\left[e^{i\chi\cdot(\hat{r}-\alpha)}\rho\right],$
(40) $\displaystyle=$
$\displaystyle\,\int_{\mathbb{R}^{d}}\\!\mathrm{d}y\,e^{iy\cdot
p/\hbar}\langle x+y/2|\rho|x-y/2\rangle.$ (41)

where $\hat{r}=(\hat{x},\hat{p})$ and $|y\rangle$ denotes a position
eigenstate. In general it takes both positive and negative values, often
called a “quasiprobability distribution.” For more about the Wigner function,
see e.g. the review in Curtright  _et al._ (2014).

##  Appendix C Symplectic and Hamiltonian matrices

In this section we recall basic facts about symplectic and Hamiltonian
matrices. Our starting point is the symplectic form

$\displaystyle\omega=\begin{pmatrix}0&\mathds{1}_{d}\\\
-\mathds{1}_{d}&0\end{pmatrix}$ (42)

associated with the phase space $\mathbb{R}^{2d}$ for $d$ classical degrees of
freedom. A symplectic matrix $A$ satisfies

$A^{\top}\omega A=\omega.$ (43)

Note that the symplectic matrices $\mathrm{Sp}(2d,\mathbb{R})$ form a group
which is a subgroup of the special linear group $\mathrm{SL}(2d,\mathbb{R})$.

For an one-parameter family $A(T)$ of symplectic matrices,

$\begin{split}0&=\frac{\mathrm{d}}{\mathrm{d}t}(A^{\top}\omega A)\\\
&=\dot{A}^{\top}\omega A+A^{\top}\omega\dot{A}.\end{split}$ (44)

Taking a family starting at the origin, $A(0)=\mathds{1}_{2d}$, we see that
Lie algebra $\mathrm{sp}(2d,\mathbb{R})$ for the symplectic group
$\mathrm{Sp}(2d,\mathbb{R})$ is the set of Hamiltonian matrices, which satisfy

$F^{\top}\omega+\omega F=0.$ (45)

(Hamiltonian matrices should not to be confused with the classical Hamiltonian
variable [energy] or the quantum Hamiltonian operator.) Using
$\omega^{2}=\mathds{1}_{2d}$ and $\omega^{\top}=-\omega$, we can rearrange
this as

$F^{\top}=\omega F\omega=-\omega^{\top}F\omega.$ (46)

Using the identities $\omega=A^{\top}\omega A=A^{-\top}\omega A^{-1}$ we can
see that if $F$ is Hamiltonian and $A$ is symplectic, then $A^{-1}FA$ is also
symplectic.

When $A(0)$ is not the identity we can rearrange (44) by left-multiplying by
$A^{-\top}$ and right-multiplying by $A^{-1}$:

$A^{-\top}\dot{A}^{\top}\omega+\omega\dot{A}A^{-1}=0.$ (47)

Thus $\dot{A}A^{-1}$ is Hamiltonian when $A(t)$ is symplectic. Conjugating by
a $A$ we can see that $A^{-1}\dot{A}$ is also Hamiltonian in this case.

When $A$ is additionally symmetric, it follows that $\dot{A}$ is also
symmetric so we can write

$\dot{A}=FA+AF^{\top}$ (48)

for some Hamiltonian matrix $F$.

In order that $\sigma$ is the covariance matrices of a pure Gaussian state, it
most be positive definite ($\sigma>0$), and $\sigma/(\hbar/2)$ must be
symplectic ($\sigma^{\top}\omega\sigma=(\hbar/2)^{2}\omega$) Simon  _et al._
(1994); Weedbrook  _et al._ (2012). As can be checked, the same is true for

$\displaystyle\sigma_{*}:=\frac{\hbar}{2}\left(\begin{matrix}\eta_{H}^{-1}\mathds{1}_{d}&0\\\
0&\eta_{H}\mathds{1}_{d}\end{matrix}\right).$ (49)

In the proof of Theorem 1 we work with such covariance matrices $\sigma$ that
additionally satisfy the matrix inequality

$\displaystyle\sigma\leq z\sigma_{*}.$ (50)

for some $z\geq 1$. Inverting this (and using the positivity of $\sigma$ and
$\sigma_{*}$) we obtain the inequality

$\sigma^{-1}\geq z^{-1}\sigma_{*}^{-1}.$ (51)

Now rewriting the symplectic conditions on $\sigma/(\hbar/2)$ and
$\sigma_{*}/(\hbar/2)$ as

$\displaystyle\sigma^{-1}$
$\displaystyle\,=(\hbar/2)^{-2}\omega^{\top}\sigma\omega,$ (52)
$\displaystyle\sigma_{*}^{-1}$
$\displaystyle\,=(\hbar/2)^{-2}\omega^{\top}\sigma_{*}\omega,$ (53)

we derive from (51) the inequality

$\omega^{\top}\sigma\omega\geq z^{-1}\omega^{\top}\sigma_{*}\omega.$ (54)

De-conjugating by the symplectic form $\omega$ gives our desired equivalent
NTS condition:

$\displaystyle\sigma\geq z^{-1}\sigma_{*}$ (55)

That is, we conclude that for pure Gaussian states being “not too long” is
equivalent to being “not too thin”:

$\sigma\geq z^{-1}\sigma_{*}\iff\sigma\leq z\sigma_{*}.$ (56)

More generally, for a pure Gaussian state, the eigenvalues of the covariance
matrix $\sigma$ come in pairs $\lambda$ and $(\hbar/2)^{2}\lambda^{-1}$ Simon
_et al._ (1994); Weedbrook  _et al._ (2012).

##  Appendix D Harmonic approximation error bound

######  Lemma 1 (Error in harmonic approximation).

For a Lindblad equation (1) with Hamiltonian
$\hat{H}=\hat{p}^{2}/2m+V(\hat{x})$ and linear lindblad operators
$\hat{L}_{k}$, the error
$\delta\hat{\mathcal{L}}^{(\alpha)}:=\hat{\mathcal{L}}-\hat{\mathcal{L}}^{(\alpha)}$
for the local harmonic approximation $\hat{\mathcal{L}}^{(\alpha)}$ to the
quantum dynamics at $\alpha$ acting on the pure Gaussian quantum state
$\hat{\tau}_{\alpha,\sigma}$ satisfies

$\displaystyle\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}$
$\displaystyle\leq\sqrt{\frac{5d^{3}}{3}}\|\nabla^{3}V\|\frac{\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{\frac{3}{2}}}{\hbar}$
(57)

Likewise for a Fokker-Planck equation (2), the error
$\delta\mathcal{L}^{(\alpha)}:=\mathcal{L}-\mathcal{L}^{(\alpha)}$ for the
local harmonic approximation $\mathcal{L}^{(\alpha)}$ to the classical
dynamics at $\alpha$ acting on the Gaussian classical state
$\tau_{\alpha,\sigma}$ satisfies

$\displaystyle\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}$
$\displaystyle\leq\sqrt{3d^{3}}\|\nabla^{3}V\|\frac{\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{\frac{2}{3}}}{\hbar}$
(58)

###### Proof.

We first bound the quantum error

$\displaystyle\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}$
$\displaystyle=\left\lVert(\hat{\mathcal{L}}-\hat{\mathcal{L}}^{(\alpha)})[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}=\left\lVert-\frac{i}{\hbar}\left[\delta\hat{H}^{[\alpha,2]},\hat{\tau}_{\alpha,\sigma}\right]\right\rVert_{\mathrm{Tr}}\leq\frac{2}{\hbar}\left\lVert\delta\hat{H}^{[\alpha,2]}\hat{\tau}_{\alpha,\sigma}\right\rVert_{\mathrm{Tr}}$
(59)

for the Gaussian quantum state
$\hat{\tau}_{\alpha,\sigma}=|\alpha,\sigma\rangle\langle\alpha,\sigma|$ with
covariance matrix $\sigma$ and mean $\alpha$. Here, $\hat{I}$ denotes the
identity operator and
$\delta\hat{H}^{[\alpha,2]}=\hat{H}-\hat{H}^{[\alpha,2]}=\delta
V^{[\alpha^{\mathrm{x}},2]}(\hat{x})$ is the operator error from the harmonic
approximation. By Taylor’s theorem we have the classical remainder from the
quadratic approximation

$\displaystyle\begin{split}\delta
V^{[\alpha^{\mathrm{x}},2]}(\alpha^{\mathrm{x}}+\Delta
x)&=V(\alpha^{\mathrm{x}}+\Delta
x)-V^{[\alpha^{\mathrm{x}},2]}(\alpha^{\mathrm{x}}+\Delta
x)=\frac{1}{3!}\left[(\Delta
x\cdot\nabla)^{3}V\right](\alpha^{\mathrm{x}}+\xi\Delta x)\end{split}$ (60)

for some choice of $\xi\in[0,1]$ (depending on $\alpha^{\mathrm{x}}$ and
$\Delta x$).

Recalling that
$\left\lVert\mathinner{|{\psi}\rangle}\mathinner{\langle{\phi}|}\right\rVert_{\mathrm{Tr}}^{2}=\|\psi\|^{2}\|\phi\|^{2}$,
we have

$\displaystyle\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}^{2}\leq$
$\displaystyle\
\frac{4}{\hbar^{2}}\left\lVert\delta\hat{H}^{[\alpha,2]}|\alpha,\sigma\rangle\right\rVert^{2}$
(61) $\displaystyle=$ $\displaystyle\
\frac{4}{\hbar^{2}}\operatorname{Tr}\left[\hat{\tau}_{\alpha,\sigma}(\delta
V^{[\alpha^{\mathrm{x}},2]}(\hat{x}))^{2}\right]$ (62) $\displaystyle=$
$\displaystyle\
\frac{4}{\hbar^{2}}\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)\left[\delta
V^{[\alpha^{\mathrm{x}},2]}(\alpha^{\mathrm{x}}+\beta^{\mathrm{x}})\right]^{2}$
(63)

where
$\tau_{\alpha,\sigma}(\alpha+\beta)=\mathcal{W}[\hat{\tau}_{\alpha,\sigma}](\alpha+\beta)=\exp(-\beta^{a}\sigma^{-1}_{ab}\beta^{b}/2)/[(2\pi)^{d}\sqrt{\det\sigma}]$
(a positive-valued function on phase space) is the Wigner function of the pure
Gaussian state
$\hat{\tau}_{\alpha,\sigma}=|\alpha,\sigma\rangle\langle\alpha,\sigma|$. In
Eq. 63 we have made use of the Weyl trace formula in the particularly simple
case of the expectation value of an operator that is a function of position:
$\operatorname{Tr}[\rho
V(\hat{x})]=\int\\!\mathrm{d}\alpha\,\mathcal{W}[\rho](\alpha)V(\alpha^{\mathrm{x}})$.
Next we apply the approximation (60) from Taylor’s theorem:

$\displaystyle\left\lVert\delta\hat{\mathcal{L}}^{(\alpha)}[\hat{\tau}_{\alpha,\sigma}]\right\rVert_{\mathrm{Tr}}^{2}\leq$
$\displaystyle\
\frac{4}{\hbar^{2}}\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)\left[\frac{1}{3!}[(\beta^{\mathrm{x}}\cdot\nabla)^{3}V](\alpha^{\mathrm{x}}+\xi(\beta^{\mathrm{x}})\beta^{\mathrm{x}})\right]^{2}$
(64) $\displaystyle\leq$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{9\hbar^{2}}\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)|\beta^{\mathrm{x}}|^{6}$
(65) $\displaystyle=$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{9\hbar^{2}}\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)(\beta^{\top}P_{\mathrm{x}}\beta)^{3}$
(66) $\displaystyle=$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{9\hbar^{2}}\left[(\operatorname{Tr}\sigma^{\mathrm{x}\mathrm{x}})^{3}+6(\operatorname{Tr}\sigma^{\mathrm{x}\mathrm{x}})\operatorname{Tr}((\sigma^{\mathrm{x}\mathrm{x}})^{2})+8\operatorname{Tr}((\sigma^{\mathrm{x}\mathrm{x}})^{3})\right]$
(67) $\displaystyle\leq$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{9\hbar^{2}}\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{3}\left(d^{3}+6d^{2}+8d\right)$
(68) $\displaystyle\leq$
$\displaystyle\frac{5\|\nabla^{3}V\|^{2}}{3\hbar^{2}}\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{3}d^{3}$
(69)

where $P_{\mathrm{x}}:=\left(\begin{smallmatrix}\mathds{1}_{d}&0\\\
0&0\end{smallmatrix}\right)$ projects onto the $\mathrm{x}$ block. In (65) we
have used the Cauchy-Schwartz inequality, in (67) we have performed the
Gaussian integral (as recalled in Appendix E.2) and in (68) we have used
$\operatorname{Tr}[A^{n}]\leq\|A\|_{\mathrm{op}}^{n}d^{n}$ for $d\times d$
positive semidefinite matrix $A$. Eq. (68) implies Eq. (57) because $d\geq 1$.

Now we turn to the classical harmonic error:

$\displaystyle\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}$
$\displaystyle=\left\lVert(\mathcal{L}-\mathcal{L}^{(\alpha)})[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}=\left\lVert-\omega^{ab}\partial_{a}(\tau_{\alpha,\sigma}\partial_{b}\delta
H^{[\alpha,2]})\right\rVert_{L^{1}}=\left\lVert(\nabla_{\mathrm{p}}\tau_{\alpha,\sigma})\cdot(\nabla_{\mathrm{x}}\delta
V^{[\alpha^{\mathrm{x}},2]})\right\rVert_{L^{1}}$ (70)
$\displaystyle=\left\lVert-\tau_{\alpha,\sigma}(\sigma^{-1}\beta)^{\mathrm{p}}\cdot(\nabla_{\mathrm{x}}\delta
V^{[\alpha^{\mathrm{x}},2]})\right\rVert_{L^{1}}$ (71)
$\displaystyle=\|\tau_{\alpha,\sigma}^{1/2}\|_{L^{2}}^{2}\,\|\tau_{\alpha,\sigma}^{1/2}(\sigma^{-1}\beta)^{\mathrm{p}}\cdot(\nabla_{\mathrm{x}}\delta
V^{[\alpha^{\mathrm{x}},2]})\|_{L^{2}}^{2}$ (72)

where the Gaussian derivative
$\partial_{a}\tau_{\alpha,\sigma}=-\sigma^{-1}_{ab}\beta^{b}\tau_{\alpha,\sigma}$
is recalled in Appendix E.1 and where in the last line we have used the
Cauchy-Schwartz inequality. (Here, $\beta$ in the norm is understood to
represent the function $f(\alpha+\beta)=\beta$.) Then because
$\|\tau_{\alpha,\sigma}^{1/2}\|_{L^{2}}^{2}=\int\\!\mathrm{d}\beta\,|\tau_{\alpha,\sigma}(\alpha+\beta)|=1$
we have

$\displaystyle\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}^{2}=\,$
$\displaystyle\int\\!\mathrm{d}\beta\,\tau_{\alpha,\sigma}(\alpha+\beta)\left|(\sigma^{-1}\beta)^{\mathrm{p}}\cdot[(\nabla_{\mathrm{x}}\delta
V^{[\alpha^{\mathrm{x}},2]})(\alpha^{\mathrm{x}}+\beta^{\mathrm{x}})]\right|^{2}$
(73) $\displaystyle=\,$
$\displaystyle\int\\!\mathrm{d}\beta\,\tau_{\alpha,\sigma}(\alpha+\beta)\left|(\sigma^{-1}\beta)^{\mathrm{p}}\cdot\frac{1}{2!}[(\beta^{\mathrm{x}}\cdot\nabla_{\mathrm{x}})^{2}\nabla_{\mathrm{x}}V](\alpha^{\mathrm{x}}+\xi(\beta^{\mathrm{x}})\beta^{\mathrm{x}})\right|^{2}$
(74) $\displaystyle=\,$
$\displaystyle\int\\!\mathrm{d}\beta\,\tau_{\alpha,\sigma}(\alpha+\beta)\left|\frac{1}{2!}[(\beta^{\mathrm{x}}\cdot\nabla_{\mathrm{x}})^{2}((\beta\sigma^{-1})^{\mathrm{p}}\cdot\nabla_{\mathrm{x}})V](\alpha^{\mathrm{x}}+\xi(\beta^{\mathrm{x}})\beta^{\mathrm{x}})\right|^{2}$
(75) $\displaystyle\leq\,$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{4}\int\\!\mathrm{d}\beta\,\tau_{\alpha,\sigma}(\alpha+\beta)\left[|\beta^{\mathrm{x}}|^{2}|(\beta\sigma^{-1})^{\mathrm{p}}|\right]^{2}$
(76) $\displaystyle\leq\,$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{4}\int\\!\mathrm{d}\beta\,\tau_{\alpha,\sigma}(\alpha+\beta)(\beta^{\top}P_{\mathrm{x}}\beta)^{2}(\beta^{\top}\sigma^{-1}P_{\mathrm{p}}\sigma^{-1}\beta)$
(77) $\displaystyle\leq\,$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{4}\left[(\operatorname{Tr}\sigma^{\mathrm{x}\mathrm{x}})^{2}\operatorname{Tr}(\sigma^{-1}_{\mathrm{p}\mathrm{p}})+2\operatorname{Tr}(\sigma^{-1}_{\mathrm{p}\mathrm{p}})\operatorname{Tr}[(\sigma^{\mathrm{x}\mathrm{x}})^{2}]\right]$
(78) $\displaystyle=\,$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{4}d^{2}(d+2)\|\sigma^{-1}_{\mathrm{p}\mathrm{p}}\|_{\mathrm{op}}\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{2}$
(79)

where we have again used the Cauchy-Schwartz inequality and the Gaussian
integrals reviewed in Appendix E.2. (In particular,
$\operatorname{Tr}[\sigma\sigma^{-1}P_{\mathrm{p}}\sigma^{-1}]=\operatorname{Tr}[\sigma^{-1}_{\mathrm{p}\mathrm{p}}]$
and $\operatorname{Tr}[\sigma
P_{\mathrm{x}}\sigma\sigma^{-1}P_{\mathrm{p}}\sigma^{-1}]=\operatorname{Tr}[P_{\mathrm{x}}P_{\mathrm{p}}]=0$.)
When $\sigma$ is the covariance of a pure Gaussian quantum state, its the
eigenvalues come in pairs $\lambda$, $\hbar^{2}/(4\lambda)$ that are
associated with symplectically conjugate directions (Appendix C). This means
$\|\sigma^{-1}_{\mathrm{p}\mathrm{p}}\|_{\mathrm{op}}=(4/\hbar^{2})\|\sigma_{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}$
so that

$\displaystyle\left\lVert\delta\mathcal{L}^{(\alpha)}[\tau_{\alpha,\sigma}]\right\rVert_{L^{1}}^{2}\leq\,$
$\displaystyle\frac{\|\nabla^{3}V\|^{2}}{\hbar^{2}}d^{2}(d+2)\|\sigma^{\mathrm{x}\mathrm{x}}\|_{\mathrm{op}}^{3}$
(80)

which implies (58) because $d\geq 1$. ∎

For $N$ particles in $n$ spatial dimensions with $k$-wise interactions, the
dimension is $d=Nn$ but the $d^{\frac{3}{2}}$ dimensional factor in Lemma 1
can be replaced with a factor of only $k^{\frac{3}{2}}n^{\frac{3}{2}}$, a more
favorable scaling. One therefore expects extensive error for many-body
systems. This growth in error is analogous to that of the orthogonality
catastrophe, and therefore a more local notion of error seems to be needed to
study the quantum-classical correspondence in many-body systems.

##  Appendix E Gaussian derivatives and integrals

###  E.1 Gaussian derivatives

The Gaussian probability distribution with mean $\alpha$ and covariance matrix
$\sigma$ is

$\displaystyle\tau_{\alpha,\sigma}(\alpha+\beta)=\frac{e^{-\beta^{\top}\sigma^{-1}\beta/2}}{(2\pi)^{d}\sqrt{\det\sigma}}=\frac{1}{(2\pi)^{d}\sqrt{\det\sigma}}\exp\left(-\frac{1}{2}\beta^{a}\sigma^{-1}_{ab}\beta^{b}\right)$
(81)

Let us consider this a real-valued function of any vector $\beta$ and any
invertible matrix $\sigma$, including non-symmetric ones, so that
$\sigma_{ab}$ and $\sigma_{ba}$ are independent variables for the purposes of
partial derivatives. However, at the end we will evaluate these derivatives on
the subspace where $\sigma$ is symmetric. Recalling our notation
$\partial_{c}=\partial/\partial\beta^{c}$ so
$\partial_{c}\beta^{a}=\delta_{c}^{\phantom{a}a}$, we have

$\displaystyle\partial_{d}\left(\beta^{a}\sigma^{-1}_{ab}\beta^{b}\right)\,$
$\displaystyle=\sigma^{-1}_{db}\beta^{b}+\beta^{a}\sigma^{-1}_{ad}$ (82)
$\displaystyle\partial_{c}\partial_{d}\left(\beta^{a}\sigma^{-1}_{ab}\beta^{b}\right)\,$
$\displaystyle=\sigma^{-1}_{dc}+\sigma^{-1}_{cd}.$ (83)

We also deploy the standard Petersen and Pedersen (2012) matrix derivative
identities444Some sources will express this for a not-necessarily invertible
matrix $Z$ using the matrix adjugate $\mathrm{adj}(Z)$. A property of the
adjugate is that $Z\,\mathrm{adj}(Z)=(\det Z)\mathds{1}$ so that, when $Z$ is
invertible, $\mathrm{adj}(Z)=(\det Z)Z^{-1}$.

$\displaystyle\frac{\partial\det Z}{\partial y}$ $\displaystyle=(\det
Z)\operatorname{Tr}\left[Z^{-1}\frac{\partial Z}{\partial y}\right],$ (84)
$\displaystyle\frac{\partial Z^{-1}}{\partial y}$
$\displaystyle=-Z^{-1}\frac{\partial Z}{\partial y}Z^{-1}$ (85)

for an invertible matrix $Z$, so in particular

$\displaystyle\frac{\partial\det Z}{\partial Z^{ab}}$ $\displaystyle=(\det
Z)Z^{-1}_{ba},$ (86) $\displaystyle\frac{\partial Z^{-1}_{cd}}{\partial
Z^{ab}}$ $\displaystyle=-Z^{-1}_{ca}Z^{-1}_{bd}.$ (87)

Combining these we get

$\displaystyle\partial_{a}\partial_{b}\tau_{\alpha,\sigma}(\alpha+\beta)=(\sigma^{-1}_{ac}\beta^{c}\sigma^{-1}_{bd}\beta^{d}-\sigma^{-1}_{ab})\tau_{\alpha,\sigma}(\alpha+\beta)=2\frac{\partial}{\partial\sigma^{ab}}\tau_{\alpha,\sigma}(\alpha+\beta),$
(88)

when evaluated for symmetric $\sigma$. (As expected, this is singular when
$\sigma$ is non-invertible.) Weyl quantizing both sides with
$\operatorname{Op}_{\hbar}=\mathcal{W}^{-1}$ gives the corresponding quantum
expression
$\partial_{a}\partial_{b}\hat{\tau}_{\alpha,\sigma}=2\frac{\partial}{\partial\sigma_{ab}}\hat{\tau}_{\alpha,\sigma}$.

###  E.2 Gaussian integrals

Here we recall the evaluation of some Gaussian integrals, as can be done with
Wick’s theorem. We define the shorthand:

$\displaystyle\begin{split}\langle(\beta^{\top}A\beta)\rangle_{\sigma}:=&\int\mathrm{d}\beta\tau_{0,\sigma}(\beta)(\beta^{\top}A\beta)\\\
=&\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)(\beta^{\top}A\beta)\\\
=&A_{ab}\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)\beta^{a}\beta^{b}\\\
=&A_{ab}\sigma^{ab}\\\ =&\operatorname{Tr}[\sigma A].\end{split}$ (89)

for any positive semidefinite matrix $A$. (The covariance matrix $\sigma$ is
also positive semidefinite, of course.) Likewise, for $B$ and $C$ also
positive semidefinite, we have

$\displaystyle\begin{split}\langle(\beta^{\top}A\beta)(\beta^{\top}B\beta)\rangle_{\sigma}:=&\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)(\beta^{\top}A\beta)(\beta^{\top}B\beta)\\\
=&A_{ab}B_{cd}\left[\sigma^{ab}\sigma^{cd}+2\sigma^{ad}\sigma^{bc}\right]\\\
=&\operatorname{Tr}[\sigma A]\operatorname{Tr}[\sigma
B]+2\operatorname{Tr}[\sigma A\sigma B]\end{split}$ (90)

and

$\displaystyle\begin{split}\langle(\beta^{\top}A\beta)(\beta^{\top}B\beta)(\beta^{\top}C\beta)\rangle_{\sigma}:=&\int\mathrm{d}\beta\tau_{\alpha,\sigma}(\alpha+\beta)(\beta^{\top}A\beta)(\beta^{\top}B\beta)(\beta^{\top}C\beta)\\\
=&A_{ab}B_{cd}C_{ef}\Big{[}\sigma^{ab}\sigma^{cd}\sigma^{ef}+2\left(\sigma^{ab}\sigma^{cf}\sigma^{de}+\sigma^{af}\sigma^{cd}\sigma^{be}+\sigma^{ad}\sigma^{bc}\sigma^{ef}\right)\\\
&\qquad\qquad\qquad+4\left(\sigma^{ad}\sigma^{be}\sigma^{cf}+\sigma^{af}\sigma^{bc}\sigma^{de}\right)\Big{]}\\\
=&\operatorname{Tr}[\sigma A]\operatorname{Tr}[\sigma
B]\operatorname{Tr}[\sigma C]+2\operatorname{Tr}[\sigma
A]\operatorname{Tr}[\sigma B\sigma C]+2\operatorname{Tr}[\sigma
B]\operatorname{Tr}[\sigma C\sigma A]\\\ &\qquad+2\operatorname{Tr}[\sigma
C]\operatorname{Tr}[\sigma B\sigma A]+8\operatorname{Tr}[\sigma A\sigma
B\sigma C].\end{split}$ (91)

##  Appendix F Applications to physical examples

Previous studies have estimated the strength of the environment-induced
diffusion for real physical systems: a grain of dust being decohered by the
cosmic microwave background, a large molecule being decohered by sunlight, and
so on.

We encounter a difficulty when plugging these physical estimates into our
bound. In many naturally occurring interactions between system and bath, only
the position variables are coupled to the bath, with no direct coupling to
momentum. Thus $D_{(\mathrm{x})}=\hbar|\ell|^{2}_{\mathrm{p}}=0$ while
$D_{(\mathrm{p})}=\hbar|\ell_{\mathrm{x}}|^{2}\neq 0$, i.e., there is
diffusion in momentum but not position. On the other hand, if we unpack the
constants in Eq. (8),

$\displaystyle\begin{split}\epsilon&=t\hbar^{2}d^{\frac{3}{2}}\frac{\|\nabla^{3}V\|}{\left(\|\nabla^{2}V\|m\right)^{\frac{3}{4}}}\max\left\\{\frac{1/m}{D_{(\mathrm{x})}},\frac{\|\nabla^{2}V\|}{D_{(\mathrm{p})}},\frac{1}{\hbar}\right\\}^{\frac{3}{2}},\end{split}$
(92)

one can see our bound deteriorates if either $D_{(\mathrm{x})}$ or
$D_{(\mathrm{p})}$ become too small. At first glance, this problem may seem
irreparable. The diffusion in position variables, which is missing in most
natural models of system-bath interactions, corresponds to decohering
superpositions in the system’s momentum variables. Yet if superpositions of
momentum are allowed to persist, why should we expect classical behavior?

Physically, the remedy is that even if a superposition
$|\psi_{1}\rangle+|\psi_{2}\rangle$ of two coherent states with different
momenta but overlapping positions is not _immediately_ decohered by purely
positional decoherence, they will ultimately evolve into a superposition of
two states with different positions as well, at which point the positional
decoherence takes effect. Thus heuristically, we can model an effective
positional diffusion $D_{(\mathrm{x}),\mathrm{eff}}$ as
$D_{(\mathrm{x}),\mathrm{eff}}\propto D_{(\mathrm{p})}$ even when
$D_{(\mathrm{x})}$ is actually zero. In particular, simple dimensional
analysis suggests we can model $D_{(\mathrm{x}),\mathrm{eff}}$ as

$D_{(\mathrm{x}),\mathrm{eff}}\sim\frac{D_{(\mathrm{p})}}{\|\nabla^{2}V\|m}.$
(93)

In that case, we ignore the first argument in the $\min\\{\cdot\\}$ of (92),
or equivalently take $D_{0}=D_{(\mathrm{p})}\tau_{H}/p_{H}^{2}$ in (6). We
emphasize this substitution is heuristic, though it is justified by work in
progress.

With this substitution, we can apply our bound to examples using Table 3.1 of
Joos  _et al._ (2013), which provides estimates for the “localization rate”
$\Lambda_{(\mathrm{x})}\sim\hbar^{-2}D_{(\mathrm{p})}$, quantifying positional
decoherence. As an example, we consider a dust particle with diameter
$10\mathrm{\mu m}$ and mass $10^{-11}$ kg decohering in the sunlight, for
which $\Lambda_{(\mathrm{x})}\sim 10^{25}\mathrm{m}^{-2}\mathrm{s}^{-1}$.

Consider a particle with characteristic energy $E\sim mv^{2}$ for
characteristic velocity $v$, and characteristic length scale $s$ for the
variation of the potential. That is, we take $\|\nabla^{2}V\|\sim Es^{-2}$ and
$\|\nabla^{3}V\|\sim Es^{-3}$. Then we have

$\left\lVert{\widetilde{\rho}}(t)-\rho(t)\right\rVert_{\mathrm{Tr}}\lesssim
t\hbar^{-1}v^{\frac{7}{2}}m\Lambda^{-\frac{3}{2}}_{(\mathrm{x})}s^{-\frac{9}{2}},$
(94)

i.e. the bound is useful for time

$t\lesssim\hbar v^{-\frac{7}{2}}m^{-1}\Lambda^{\frac{3}{2}}s^{\frac{9}{2}}.$
(95)

For the $10$ micron dust particle in sunlight, with characteristic velocity
$v=1\mathrm{m/s}$ in a potential varying on the length scale of
$s=1\,\mathrm{m}$, the error
$\left\lVert\rho(t)-\tilde{\rho}(t)\right\rVert\ll 1$ is small for three
million years.