Wigner Analysis of Particle Dynamics
Decoherence in Wide Nonharmonic Potentials

andAndreu Riera-Campeny1,2 , Marc Roda-Llordes1,2 , Piotr T. Grochowski1,2,3 , andOriol Romero-Isart1,2,4,5
1

Institute for Quantum Optics and Quantum Information of the Austrian Academy of Sciences, A-6020 Inns-

arXiv:2307.14106v6 [quant-ph] 2 Sep 2024

bruck, Austria
2Institute for Theoretical Physics, University of Innsbruck, A-6020 Innsbruck, Austria

3

Center for Theoretical Physics, Polish Academy of Sciences, Aleja Lotników 32/46, 02-668 Warsaw, Poland

4

ICFO - Institut de Ciències Fotòniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels(Barcelona), Spain

5

ICREA, Pg. Lluis Companys 23, 08010 Barcelona, Spain

We derive an analytical expression of a Wigner function that approximately
describes the time evolution of the one-dimensional motion of a particle in anonharmonic potential. Our method involves two exact frame transformations,accounting for both the classical dynamics of the centroid of the initial state andthe rotation and squeezing about that trajectory. Subsequently, we employ two
crucial approximations, namely the constant-angle and linearized-decoherence
approximations, upon which our results rely. These approximations are effective in the regime of wide potentials and small fluctuations, namely potentials
that enable spatial expansions orders of magnitude larger than the one of theinitial state but that remain smaller compared to the relevant dynamical lengthscale (e.g., the distance between turning points). Our analytical result elucidates the interplay between classical and quantum physics and the impact of
decoherence during nonlinear dynamics. This analytical result is instrumental
to designing, optimizing, and understanding proposals using nonlinear dynamics to generate macroscopic quantum states of massive particles.

1 Introduction
The one-dimensional position and momentum of a particle is perhaps the most fundamental degree of freedom in physics and of paramount importance in the development of
quantum mechanics [1–3]. The state of such a continuous-variable degree of freedom can be
described by the Wigner function [4–6], a phase-space quasi-probability distribution. The
Wigner function and its time evolution allow us to study the interplay between classicaland quantum mechanics as well as the impact of noise and decoherence [7–9]. Usually, thedynamics of the Wigner function is studied in scenarios of small phase-space areas, i.e.,
scenarios where the phase-space surface occupied by the initial state is of the same orderas the space available to be explored through the dynamics. From a theoretical viewpoint,
scenarios with small phase-space areas are advantageous as they circumvent the numericalinstabilities associated with the nonlocal and highly oscillatory character of the Wigner
equation [10, 11], which appear for large coherent expansions due to the presence of arbitrarily high-order derivatives. Experimentally, this is the natural regime for quantum
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

1

experiments performed with photons [12–14] and atomic systems [15, 16], whose state,even at near-zero temperatures, is delocalized over scales comparable to the length scale ofthe confining potential. Motivated by recent experimental progress in controlling the quantum center-of-mass motion of a single nanoparticle [17–24], which has approximately nineorders of magnitude more mass than a single atom, here we focus on the rather unexplored
regime of large-scale quantum dynamics. That is, the dynamics generated when a highly
localized phase-space probability distribution evolves in a wide, nonharmonic potentialthat enables the particle to explore an orders-of-magnitude larger phase-space area.
In this paper, we present a method to obtain an approximated analytical expression
for the time evolution of the Wigner function in a nonharmonic potential, in the presenceof decoherence, caused by weak coupling to a high-temperature bath [25–27], and whitenoise fluctuations in both the amplitude and position of the potential [28–30]. Our method
implements two exact frame transformations. The first one accounts for the classical trajectory of the centroid of the initial state. The second incorporates the rotation and squeezingaround the centroid classical trajectory. After those exact transformations, we perform two
approximations, namely the constant-angle and linearized-decoherence approximations.
On the one hand, the constant-angle approximation is applicable when the local phasespace rotation angle (defined below) around the centroid classical trajectory changes slowly
over time. This approximation is expected to hold in the regime of wide potentials, where
the relevant length scale of the potential is much larger than the initial spatial extent of theparticle’s position, thereby allowing large coherent expansions. Importantly, the approximation breaks down close to the turning points, where the quantum state compresses andthe local phase-space rotation angle changes rapidly. On the other hand, the linearizeddecoherence approximation is accurate in the limit of small fluctuations, where the de
Broglie wavelength of the quantum state remains significantly smaller than the relevant
length scale of the potential. Both, the breakdown of the approximation near the turning
points and the small de Broglie wavelength evoke properties reminiscent of the so-called
semiclassical methods [31–34], such as the Wentzel–Kramers–Brillouin approximation, thatstudy the asymptotic limit ℏ → 0 of the Schrödinger equation. Among semiclassical methods, we highlight the connection of our method with the Gaussian wave packet dynamics
developed by Heller and co-workers in a series of works [35–40]. This method relies on
the so-called thawed Gaussian approximation [41], in which the potential is replaced by its
quadratic expansion along the classical trajectory of the wave packet. Consequently, ourmethod reproduces this approximation in the limit in which all nonlinear and decoherence
effects are disregarded. However, the precise connection between the method discussed
here and previously studied semiclassical methods remains unclear, and we leave it for
future investigation.
Here, we show that the analytical expression obtained with this method provides an
excellent approximation to the generated nonlinear open dynamics in the regime of wide
potentials and small fluctuations. Our results are timely as they synergize with recently
developed numerical tools [11] that allow to design, optimize, and understand protocols
where a quantum ground-state-cooled nanoparticle in a wide nonharmonic potential rapidly
evolves into a macroscopic quantum superposition state [42].
The rest of the paper is organized as follows. In Sec. 2, we derive an analytical approach
using Wigner analysis to describe the dynamics of a particle in a nonharmonic potential.
We define the regime in which a key approximation, referred to as the constant-angle
approximation, is used to integrate the nonlinear dynamics, thereby providing an analyticalform of the time-evolved Wigner function. In Sec. 3, we apply our analytical method tothe example of a particle evolving in a wide double-well potential. We explicitly show how

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

2

our analytical approach reproduces the numerical results obtained using the split-method
operator [43] and the numerical method presented in [11]. We draw our conclusions in
Sec. 4 and provide further details of our analysis in the Appendix.

2 Wigner Analysis of Particle Nonlinear Dynamics
In this Section, we provide an analytical treatment of the one-dimensional dynamics ofa particle evolving in a nonharmonic potential in the presence of decoherence. Since theclassical equations of motion associated with a nonharmonic potential are nonlinear, we
say that the particle undergoes nonlinear dynamics. Initially (Sec. 2.1), we will delve
into the detailed description of the dynamical problem under study, presenting both the
Liouville-von Neumann equation and the equation of motion of the Wigner function. Subsequently, we will conduct two frame transformations that are exact: the classical centroid
frame transformation (Sec. 2.2), and the Gaussian frame transformation (Sec. 2.3). These
two transformations yield a valuable and exact reformulation of the dynamical problem.
Following this, we will perform two approximations (Sec. 2.4): the constant-angle approximation and the linearized-decoherence approximation. These will lead to an easily
evaluable expression of the time-evolved Wigner function. Finally (Sec. 2.5), we will discuss the regime where we expect the constant-angle approximation to provide an accurate
description of the nonlinear particle dynamics.
Throughout the section, we will use the following nomenclature for the different typesof dynamics. First, the term linear dynamics refers to the dynamics generated by a Hamiltonian whose associated classical equations of motion are linear. That is, whose solutions
are closed under linear combinations. This implies that the corresponding Hamiltonian is
at most quadratic in the position and momentum variables. The term coherent dynamics
refers to the dynamics that preserves the coherence and purity of the state over time. In
the wave function picture, those are the dynamics generated by the Schrödinger equation.
Finally, we use the term Gaussian dynamics to refer to the dynamics that preserves theGaussianity of the state. Namely, the combination of linear dynamics plus diffusion-likedissipation. Finally, given the definitions above, we use interchangeably the terms linear
dynamics and coherent Gaussian dynamics.

2.1 Description of the dynamical problem
We consider the one-dimensional dynamics of a particle with mass m, where the position
and momentum operators are given by X̂ and P̂ , respectively, and they fulfill the canonicalcommutation rule [X̂, P̂ ] = iℏ. The particle evolves under a potential denoted by V (X̂),
and in the presence of decoherence. The time evolution of the particle’s state, specified byits density operator ρ̂(t), is governed by the Liouville-von Neumann equation:
"

#

∂ ρ̂
1 P̂ 2(t) =+ V (X̂), ρ̂(t) + D[ρ̂(t)],
∂t
iℏ 2m

(1)

Here, D[·] represents the decoherence superoperator, which we specify below. We areinterested in studying the evolution of a particle that is initially cooled to a low-temperaturethermal state in a harmonic potential of frequency Ω (for example, an
p optical trap). At zerotemperature, the particle has zero-point fluctuations, with XΩ = ℏ/(2mΩ) in position,
and PΩ = ℏ/(2XΩ ) in momentum.
Motivated by experiments with levitated nanoparticles in ultra-high vacuum [17–24,
44], we will model the decoherence as follows. First, we consider position localization
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

3

decoherence [45–47], which is defined by the dissipator:
Dloc [·] = −Γloc[X̂, [X̂, ·]],
2XΩ2

(2)

where Γloc is the position localization decoherence rate. This type of decoherence emerges
in the weak-coupling high-temperature limit of Brownian motion [25–27, 48] as well as
after performing an ensemble average of a white-noise stochastic force [49]. This decoherence model mimics the recoil heating due to laser light [50, 51] and the emission ofthermal photons [46, 52–54], among other sources of decoherence. The dissipator in Dloc [·]corresponds to the long-wavelength limit of the decoherence superoperator [52, 55], andit overestimates the actual decoherence of the system. Hence, even in the case where theform in Eq. (2) is not accurate, it can be used to upper bound the effect of environmental
decoherence. Additionally to Dloc [·], we take into account the decoherence due to fluctuations in the potential’s center and its amplitude [28–30]. In particular, we consider thefluctuating potential:
Vfluc (X, t) = (1 + ξ2 (t))V (X − XΩ ξ1 (t)) ≈ V (X) − ξ1 (t)XΩ∂V(X) + ξ2 (t)V (X),∂X

(3)

where ξ1 (t) and ξ2 (t) are zero-average stochastic processes representing fluctuations of thepotential’s center and amplitude, respectively. In Eq. (3), we have assumed small potential
fluctuations, which justify the use of a Taylor series expansion in each ξi (t), truncated to
the first order. For simplicity, we assume ξi (t) to be independent Gaussian white noise
processes with the correlation function ⟨⟨ξi (t)ξj (t′ )⟩⟩ = 2πAi δij δ(t − t′ ), where ⟨⟨·⟩⟩ denotes
the average over trajectories and Ai is a measure of the noise amplitude. The average overthese stochastic processes can be computed using the cumulant expansion method to thesecond order [56], and leads to the decoherence superoperator:
Dfluc [·] = −∂V
πA2
πA1 XΩ2 ∂V[(X̂), [(X̂), ·]] − 2 [V (X̂), [V (X̂), ·]].
2
ℏ
∂X
∂X
ℏ

(4)

Putting everything together, the decoherence model used in this paper is given by D[·] =Dloc [·] + Dfluc [·]. This dissipator can be compactly written as:
D[·] = −

∞
X

∞
X
Γnm
n
m
Dnm [·].
[
X̂
,
[
X̂
,
·]]
≡
n+m
n,m=1 2XΩ
n,m=1

(5)

This form is obtained by performing a Taylor expansion of the potential, i.e., V (X) =nnnn=1 (∂ V /∂X )(0)X /n!. The expression of the decoherence rates Γnm can be obtainedby collecting the corresponding terms from Dloc [·] and Dfluc [·], and they are explicitly given
in App. A. It should be noted that the dissipator Dnm [·] for n or m strictly larger than one,which originates from the potential fluctuations Dfluc [·], generates non-Gaussian dissipative
dynamics.
We chose to use the Wigner function formalism [4, 5] to analytically handle the dynamics generated by Eq. (1). Moreover, we use the dimensionless position x = X/XΩ ,momentum p = P/PΩ , and time τ = ωt variables, where ω is an arbitrary frequency scaleassociated with the potential V (X). The corresponding dimensionless potential in theseunits is defined as U (x) ≡ V (xXΩ )/(mω 2 XΩ2 ). In accordance with Eq. (1), the equationof motion for the dimensionless Wigner function W (r, τ ) is given byP∞∂W(r, τ ) = (Lc + Lq + Ld ) W (r, τ ),∂τ
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(6)4

where r ≡ (x, p)T denotes a point in the dimensionless phase space. The first term generates classical (i.e., Liouville) dynamics and is expressed as:
ω ∂U
∂
Ω ∂
+
(x) .
Lc = − pω ∂x Ω ∂x∂p(7)

The second term generates genuine quantum dynamics and is written as:∞
ω X∂(−1)n ∂ 2n+1 U
Lq =
(x)
2n+1
Ω n=1 (2n + 1)! ∂x
∂p



2n+1

.

(8)

Note that Lc + Lq generates the same evolution as the Schrödinger equation and that
Lq = 0 for up to quadratic potentials. Lastly, the third term generates decoherence, whichPin accordance with Eq. (5), can be written as Ld = ∞n,m=1 Ld,nm , where:
Ld,nm =

n+m
X
k=2

cnmk xn+m−k



∂
∂p

k

.

(9)

The specific expression of the coefficients cnmk is provided in App. A. Note that while
we have focused on the one-dimensional motion of massive particles, our results, especiallywhen written in this dimensionless form, can be applied to describe the nonlinear dynamics
of other continuous-variable degrees of freedom (e.g., a single electromagnetic field mode).
The aim of this paper is to solve Eq. (6) for an initial Gaussian state, whose Wigner
function can be written as:
W (r, 0) = G[C(0)](r − µ(0)).

(10)

Here, G[C](r) denotes the two-dimensional Gaussian distribution:
1
r T C −1 r
G[C](r) ≡exp −2π det C
2

!

,

(11)

and µ(0) = (⟨x̂⟩(0), ⟨p̂⟩(0))T and C(0) refer to the initial mean values and the covariance matrix of ρ̂(0), respectively. The elements Cij (0), where i, j ∈ {x, p}, of the always
symmetric covariance matrix are given by:1
Cij (0) = ⟨r̂i r̂j + r̂j r̂i ⟩(0) − ⟨r̂i ⟩(0)⟨r̂j ⟩(0).
2(12)

For instance, a thermal state with a mean phonon occupation number n̄ serves as a particularly relevant initial state. This state corresponds to C(0) = (2n̄ + 1)12 , where 12
represents the two-by-two identity matrix.
In the following section, we will derive an expression for W (r, τ ), which approximately
solves Eq. (6). We obtain this expression by performing two consecutive frame transformations (i.e., the classical centroid frame and the Gaussian frame) and applying a key
approximation, namely the constant-angle approximation. We will further simplify thisexpression using the linearized-decoherence approximation.

2.2 Classical centroid frame
The first frame transformation aims to remove the dynamics of the mean position andmomentum values, effectively moving to a comoving frame that aligns with the approximate centroid of the Wigner function [35]. In our Wigner function context, such frame
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

5

transformations follow from the sometimes called disentangling theorem [57, 58], and are aconsequence of the properties of exponential dynamical maps. Essentially, they are equivalent to the transformation between different pictures in standard quantum mechanics.
In order to implement the frame transformation, we define W (C) (τ ) ≡ MC (τ )−1 W (r, τ ),where MC (τ ) is a time-dependent phase-space displacement map, expressed as:∂∂
MC (τ ) ≡ exp −rc (τ ) ∇ = exp −xc (τ ).
− pc (τ )
∂x
∂p
h

T

i





(13)

where ∇ = (∂/∂x, ∂/∂p)T is the phase space gradient. In Eq. (13), rc (τ ) ≡ (xc (τ ), pc (τ ))T
represents the classical centroid trajectory, which is the solution to dimensionless Hamilton’s equations∂xcΩ(τ ) = pc (τ ),
∂τ
ω
∂pc
ω ∂U
(τ ) = −(xc (τ )),
∂τ
Ω ∂x

(14)
(15)

with initial conditions given by rc (0) = (xc (0), pc (0))T = (⟨x̂⟩(0), ⟨p̂⟩(0))T . The actionof the map MC (τ ) on an arbitrary function f (r) translates its phase-space variables as
MC (τ )f (r) = f (r −rc (τ )). In this frame, we describe the dynamics of the Wigner function
centered at the classical trajectory that is followed by its initial position and momentum
expected value. This frame transformation does not imply any approximation. It will
be particularly useful for nonlinear dynamics that maintain a small distance between thevector rc (τ ) = (xc (τ ), pc (τ ))T and µ(τ ) = (⟨x̂⟩(τ ), ⟨p̂⟩(τ ))T in comparison to the available
phase-space dimensions. As we will demonstrate later, this will be the case for wide
potentials and small fluctuations.
In this classical centroid frame, one can show that the state W (C) (r, τ ) evolves according
to an effective time-dependent potential, represented byUeff (x, τ ) ≡

∞
X
1 ∂nU

n! ∂xnn=2(xc (τ ))xn .

(16)This effective potential, at least quadratic in position, exhibits time dependence determined by the local derivatives of the potential along the classical trajectory. Thus, in the
comoving frame with the approximate centroid of the Wigner function, the particle experiences a time-dependent potential while traversing the static potential. More specifically,
the evolution equation for W (C) (r, τ ) can be expressed as:hi∂W (C)(C)(C)(C)(r, τ ),(r, τ ) = L(C)c (τ ) + Lq (τ ) + Ld (τ ) W
∂τ
(C)

(17)

(C)

Here, Lc (τ ) and Lq (τ ) are obtained by replacing U (x) with Ueff (x, τ ) in Eq. (7) andEq. (8), respectively. The generator of decoherence also becomes time-dependent and is
P
(C)
(C)
given by Ld (τ ) ≡ MC (τ )−1 Ld MC (τ ) = ∞n,m=1 Ld,nm (τ ), with(C)

Ld,nm (τ ) =

n+m
X
k=2

cnmk [x + xc (τ )]n+m−k



∂
∂p

k

.

(18)

2.3 Gaussian frame(C)

In the centroid frame, the coherent dynamics are generated by the sum of Lc (τ ) +(C)
Lq (τ ). The second frame transformation aims to remove the coherent Gaussian dynamics
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

6

generated by the harmonic part of the effective potential Eq. (16), which we define as:
1
1 ∂2UUG (x, τ ) ≡ α(τ )x2 =(xc (τ ))x2 .
2
2 ∂x2

(19)

Here, α(τ ) represents an effective (dimensionless) spring constant, which can be either
positive or negative. To remove the coherent Gaussian dynamics, we define W (G) (r, τ ) ≡MG (τ )−1 W (C) (r, τ ), where MG (τ ) is given by:Z τ

MG (τ ) ≡ exp+

0



(C)

dτ ′ LG (τ ′ ) ,

(20)
(C)

where exp+ (·) represents the time-ordered exponential and LG (τ ) is the generator of
coherent Gaussian dynamics in the centroid frame, defined as:
Ω ∂
ω ∂UG
∂
(C)
LG (τ ) = − p+(x, τ ) .
ω ∂x Ω ∂x∂p(21)

The generator of the coherent non-Gaussian dynamics, attributed to the nonharmonic(C)(C)
terms of the potential, is given by the complementary generator LnG (τ ) ≡ Lc (τ ) +(C)(C)
Lq (τ ) − LG (τ ). While assessing if a given state of motion qualifies as a quantum nonGaussian state lies beyond the scope of this paper, we remark that it is possible to findsufficient criteria based on the observation of squeezing in nonlinear variables [59].
The map MG (τ ) generates coherent Gaussian physics, namely squeezing and phasespace rotations. Its action on an arbitrary function f (r) is MG (τ )f (r) = f (S(τ )−1 r),
where S(τ ) is a symplectic matrix. The matrix S(τ ) is obtained as the solution to thedifferential equation!

∂S(τ ) =∂τ0Ω/ω
S(τ ),−(ω/Ω)α(τ )0(22)with the initial condition S(0) = 12 . One can verify that det[S(τ )] = 1 for all τ . Therefore,
denoting Sij (τ ) with i, j ∈ {x, p} the components of S(τ ), its inverse is given by
!

−1

S(τ )

=

Spp (τ ) −Sxp (τ ).
−Spx (τ ) Sxx (τ )(23)Importantly, the action of the map MG (τ ) on the phase space variable x is given by
MG (τ )−1 xMG (τ ) ≡ η(τ )xφ(τ ) = η(τ )[cos(φ(τ ))x + sin(φ(τ ))p],(24)

where rφ(τ ) ≡ R[φ(τ )]r = (xφ(τ ) , pφ(τ ) )T are rotated phase-space variables according to
the rotation matrix!

R[φ(τ )] ≡cos[φ(τ )] sin[φ(τ )],− sin[φ(τ )] cos[φ(τ )](25)

and we have defined the two key variables
η(τ ) ≡q

Sxx (τ )2 + Sxp (τ )2 ,tan[φ(τ )] ≡Sxp (τ ).
Sxx (τ )Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(26)
(27)

7

From Eq. (24), the angle φ(τ ) can be regarded as the local rotation angle of the position
quadrature in phase space. Similar transformation rules can be derived for ∂/∂p, involving
the same functions η(τ ) and φ(τ ), and for p and ∂/∂x using similar functions that combine
Spp (τ ) and Spx (τ ) instead of Sxx (τ ) and Sxp (τ ). In summary, the map MG (τ ) generates
time-dependent phase-space rotations with an angle φ(τ ) and time-dependent squeezing
with a squeezing parameter η(τ ).
In the centroid and Gaussian frame, the evolution equation of W (G) (r, τ ) is given byhi∂W (G)(G)(G)(r, τ ) = LnG (τ ) + Ld (τ ) W (G) (r, τ ),
∂τ
(G)

(C)

(28)

(G)

(C)

!

∂
∂pφ(τ )

where LnG (τ ) ≡ MG (τ )−1 LnG (τ )MG (τ ) and Ld (τ ) ≡ MG (τ )−1 Ld (τ )MG (τ ) that also
P
(G)
(G)
can be expanded as Ld (τ ) = ∞n,m=1 Ld,nm . The explicit form of the generator of coherent
non-Gaussian dynamics is given by
∞
X
∂
(G)
LnG (τ ) =βn+1 (τ )xnφ(τ )∂pφ(τ )
n=2
∞
X

+

n

(−1)

n=1
m=0

2n + m β2n+m+1 (τ ) mxφ(τ )
2n + 1
m

!2n+1

,

(29)

where the first and second terms account for classical and quantum non-Gaussian dynamics,
respectively, and we have definedβn (τ ) ≡ω1∂nU
(xc (τ ))η(τ )n .
Ω (n − 1)! ∂xn(30)

Note that the magnitude of the parameters βn (τ ) will determine whether non-Gaussian
dynamics is generated. This shows that large squeezing, that is, a large value of η(τ ),enhances the effect of the nonharmonicities in the potential [60]. The specific form of thegenerators of decoherence is given byn+mhin+m−kX(G)
Ld,nm (τ ) =cnmk η(τ )xφ(τ ) + xc (τ )
η(τ )kk=2∂∂pφ(τ )

!k

.

(31)

Note that decoherence is also enhanced by squeezing.
We emphasize that up to this point, the analytical approach is exact; that is, noapproximations have been made. We have, however, singled out the challenging partof solving a nonlinear open quantum dynamical problem, which is to integrate Eq. (28)with the initial condition given by W (G) (r, 0) = MG (0)−1 MC (0)−1 W (r, 0). If one obtains
W (G) (r, τ ), the exact solution of Eq. (6) is given by
W (r, τ ) = MC (τ )MG (τ )W (G) (r, τ ).

(32)

Interestingly, replacing W (G) (r, τ ) by W (G) (r, 0) in Eq. (32), one obtains approximate
Gaussian dynamics that are equivalent to performing the so-called Gaussian thawed approximation [35, 41] but in the Wigner, instead of the wave function, representation. Notethat such a replacement is equivalent to ignoring all the nonlinear and decoherence effects. Moreover, as we show below, the exact reformulation of the problem presented here
allows us to identify a crucial approximation, yielding an approximated Wigner function(G)
Wa (r, τ ), that incorporates nonlinear effects while being significantly easier to calculate
than its exact counterpart.
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

8

2.4 Constant-angle and linearized-decoherence approximation
Ultimately, the reason a closed-form solution for W (G) (r, τ ) is not possible stems from the(G)fact that the angle φ(τ ) changes over time. This results in the generators LnG (τ ) and(G)
Ld (τ ) being noncommutative with themselves at different times, and the solution of the
Wigner equation being a time-ordered instead of a simple exponential map. As we shallsee, in the regime of wide potentials and small fluctuations, where xc (τ ) ≫ η(τ ) ≫ 1, it
becomes feasible to approximate the angle φ(τ ) as a piecewise constant function. Thisapproximation significantly simplifies the integration of the evolution equation for the
Wigner function. Let us focus on the simplest case, where φ(τ ) can be replaced by a single
constant angle ϕ over the integration regime of interest, while the more general case isdiscussed in App. B. In this scenario, the generator of Eq. (28) commutes with itself at
different times, and furthermore, all its summands also commute with each other. We will
first implement this approximation and then, in the subsequent subsection, discuss theregime where we expect it to provide an accurate description of the particle’s nonlinear
dynamics.
Within the constant-angle approximation, the approximated solution of W (G) (r, τ ) is
given by an exponential map (rather than a time-ordered exponential one, see Ch. 2 of [61]),
which is time-dependent and can be factorized. This can be expressed as follows:
W (G) (r, τ ) ≈ Wa(G) (r, τ ) = MnG (τ )D(G) (τ )W (G) (r, 0).

(33)Here, MnG (τ ) generates coherent non-Gaussian dynamics and is given by the exponential
map MnG (τ ) = exp[ΛnG (τ )], where∞X

ΛnG (τ ) ≡∞X∂
2n + m κ2n+m+1 (τ ) m
+
(−1)n
xϕ
m
∂pϕ n=1
2n + 1!

κn+1 (τ )xnϕ

n=2

∂
∂pϕ

!2n+1

.

(34)

m=0

The first and second terms account for classical and quantum non-Gaussian dynamics,
respectively. The time-dependent coefficients quantifying the strength of the generator ofthese coherent non-Gaussian dynamics are given byZ τ

ω
1κn (τ ) ≡dτ βn (τ ) =Ω (n − 1)!
0
′

′

Z τ
0

dτ ′

∂nU
(xc (τ ′ ))η(τ ′ )n .
∂xn

(35)

Similarly, D(G) (τ ) ≡ exp[Λd (τ )] accounts for the decoherence in the Gaussian frame,whereΛd (τ ) ≡

Z τ

dτ
0

′

∞ n+m
X
X
n,m=1 k=2

′ n+m

cnmk xc (τ )



n+m−k
η(τ ′ )
η(τ ′ ) k
1+xϕxc (τ ′ )xc (τ ′ )







∂
∂pϕ

!k

. (36)

Note that by using the constant-angle approximation, the solution to the dynamics
in Wigner space has been significantly simplified. This is because one has replaced theintegration of a complicated partial differential equation with the more efficient integrationof a set of time-dependent functions.
In addition to the constant-angle approximation, a second approximation can be performed on the decoherence map under the assumption of small fluctuations. This meansthat the following condition is met:
η(τ )
≪ 1.
xc (τ )
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(37)9This is the regime in which the position fluctuations of the particle are smaller than its
mean value. Even in the case where xc (τ ) = 0, the linearization of the decoherence
generator may still be possible, albeit for different underlying reasons (see App. C). Under
this regime, one can perform the linearized-decoherence approximation, which consists ofexpanding Eq. (36) in powers of η(τ )/xc (τ ) and keeping the lowest order terms up to thesecond order. The lowest-order terms can be written compactly as follows:σb (τ )2
Λd (τ ) ≈ Λd,a (τ ) ≡
2

∂
∂pϕ

!2

.

(38)

The time-dependent parameter σb (τ ) quantifies the strength of the accumulated decoherence and is given byσb (τ )2 ≡ 4

Z τ

dτ ′

0

Γeff (τ ′ )
η(τ ′ )2 .
ω

(39)

The effective time-dependent decoherence rate is given by Γeff (τ ) ≡ Γloc + Γfluc (τ ), with
Γloc defined below Eq. (2) and with the contribution arising from the fluctuations of thepotential


#2

"

Γfluc (τ )πω 4 A1  ∂ 2 U=(xc (τ ))
Ω
2Ω3
∂x2

2
A2 ∂U
(xc (τ ))  .
A1 ∂x





+

(40)
(G)

Note that this approximation corresponds to replacing the decoherence generator Ld (τ ) =P∞(G)n,m=1 Ld,nm (τ ) in Eq. (31), by its linear approximation defined by2Γeff (τ )(G)
Ld (τ ) ≈η(τ )2
ω

∂
∂pϕ

!2

.

(41)

Putting everything together, the constant-angle and linearized-decoherence approximations allow us to obtain an approximate solution to the open nonlinear quantum mechanicalproblem described by Eq. (6). This can be represented as follows:
Wa (r, τ ) = MC (τ )MG (τ )MnG (τ )Da(G) (τ )W (G) (r, 0)(42)= Da (τ )MC (τ )MG (τ )MnG (τ )W (G) (r, 0),

(43)

(G)

(G)

Here, we have defined Da (τ ) ≡ exp[Λd,a (τ )] and Da (τ ) ≡ MG (τ )Da (τ )MG (τ )−1 . Also,
recall that W (G) (r, 0) = MG (0)−1 MC (0)−1 W (r, 0). In the second equation, we have used
the fact that D(G) (τ ) commutes with MnG (τ ) and D(τ ) with MC (τ ). The expression in
the second equation is particularly convenient as the decoherence map is applied to thecoherently evolved state. In this case, the decoherence map is given by:"#∇T Cb (τ )∇Da (τ ) ≡ exp
,
2

(44)

T
Cb (τ ) = σb (τ )2 S(τ )R(ϕ)T (ep eTp )R(ϕ)S(τ )

(45)

where

is the blurring covariance matrix and ep = (0, 1)T is a unit vector along the momentum
component. The action of the decoherence map Da (τ ) corresponds to a convolution witha Gaussian distribution of covariance matrix Cb (τ ), namely:Z

Da (τ )f (r) =dr ′ G[Cb (τ )](r − r ′ )f (r ′ ).

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(46)

10

Here, the Gaussian distribution G[C](r) is defined as per Eq. (11). Physically, σb (τ ) determines the smallest length scale that remains unaffected by decoherence and is sometimes
referred to as the blurring distance [62–64].
We can now use the expression Wa (r, τ ) in Eq. (43) to calculate the first and second
moments. As shown in App. D.1 and D.2, the first moments µa (τ ) can be expressed asµa (τ ) = rc (τ ) − S(τ )R(ϕ)T ep

X(G)κn (t)⟨x̂n−1
(0).
ϕ ⟩

(47)

n≥3

where ⟨·⟩(G) (τ ) is the average in the Gaussian frame; that is, computed with W (G) (r, τ ).
Note that the first moments are not affected by decoherence. The second moments, namely
the covariance matrix evaluated using Wa (r, τ ), can be expressed asCa (τ ) = Cb (τ ) + S(τ )R(ϕ)⟨x̂nϕ ⟩(G) (0)C(0) −
κn (t)⟨x̂nϕ ⟩(G) (0) ⟨x̂ϕn−1 p̂ϕ ⟩(G) (0)
n≥3

T

X

0

!

!

+

X
n,m≥3

0
0
 R(ϕ)S(τ )T .
κn (t)κm (t)n+m−2 (G)n−1 (G)0 ⟨x̂ϕ⟩ (0) − ⟨x̂ϕ ⟩ (0)⟨x̂m−1⟩(G) (0)ϕ(48)

Note that decoherence appears in the covariance matrix through the addition of theblurring covariance matrix defined in Eq. (45). Nonlinear dynamics is manifested by thedependence of Ca (τ ) on the initial moments higher than the second [65]. In the next
section, we argue that the accuracy of µa (τ ) and Ca (τ ) in approximating µ(τ ) and C(τ )is a useful measure of the merit of our analytical approach.

2.5 Validity of the constant-angle approximation
Our analytical approach relies on the constant-angle and linearized-decoherence approximations, which lead to the expression Wa (r, τ ), given in Eq. (43). We note that the linearizeddecoherence approximation, performed after the constant-angle approximation, is convenient but not essential, and its validity is well characterized by the small-fluctuations
condition in Eq. (37). Regarding the constant-angle approximation, we need to compare
how closely Wa (r, τ ) approximates the solution W (r, τ ). One can show that an upperbound for the norm of the difference of the exact and approximated Wigner functions is
proportional to the functionχ(τ ; ϕ) ≡XZ τdτ ′ |βn (τ ′ )||φ(τ ′ ) − ϕ|.

(49)

n≥3 0

Indeed, this quantity is minimized if the angle φ(τ ) is constant and close to ϕ on the timescales where the nonharmonicities are large, namely when the βn (τ ) are large. In thiscontext, from the definition of the angle φ(τ ) in Eq. (27), it follows thatdφΩ/ω(τ ) =,dτη(τ )2(50)

indicating that a large delocalization η(τ ) ≫ 1, reduces the rate of change of the angle. Together with the small-fluctuations condition Eq. (37) required for the linearizeddecoherence approximation, our analytical approach requires the regime xc (τ ) ≫ η(τ ) ≫ 1
in the time scales where the generators of nonlinear dynamics cannot be neglected.
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

11

Ultimately, the best way to quantify the performance of our analytical approach is
to compare Wa (r, τ ) with the numerically evaluated W (r, τ ). In this case, we suggest
comparing observables rather than the overlap between the two states. The reason isthat large-scale dynamics generates states with very small phase-space features, even subPlanckian [8, 11] that can immediately provide a very low overlap if not captured accurately. These small features can nevertheless be irrelevant in capturing the observables ofinterest and may immediately disappear in the presence of unavoidable sources of noiseand decoherence. Therefore, we suggest quantifying the performance of the constant-angle
approximation by comparing global observable properties of the state such as their firstand second moments. This can be done by evaluating their relative error, namely,2∥µ(τ ) − µa (τ )∥,∥µ(τ )∥ + ∥µa (τ )∥(51)

2∥C(τ ) − Ca (τ )∥2,∥C(τ )∥2 + ∥Ca (τ )∥2(52)

ϵ1 (τ ) ≡with the vector norm ∥y∥ =py† y, andϵ2 (τ ) ≡q

where ∥A∥2 = tr(A† A) is the Hilbert-Schmidt norm. In summary, we recommend evaluating χ(τ ; ϕ), ϵ1 (τ ) and ϵ2 (τ ) to quantify the performance of the constant-angle approximation, something we will do in the example shown in the following section.
Finally, let us emphasize that around a turning point, say τm , the state is expected
to compress, namely to minimize η(τ ). Hence, according to Eq. (50), around a turning
point the angle φ(τ ) will not change slowly. By performing a Taylor expansion of φ(τ )around the time τm , recall Eq. (26), it follows that the time derivative of the angle arounda compression time has a Lorentzian shapedφγ(τm )(τ ) ≈,2dτ
γ(τm ) + (τ − τm )2(53)

with width γ(τm ) ≡ (ω/Ω)η(τm )2 . Therefore, by integrating this Lorentzian, we see thataround a turning point, the angle should change approximately by π. This implies that
if the constant-angle approximation is performed with an angle ϕ during the time scale
before a turning point, a second angle around ϕ + π will need to be employed after theturning point. Note that this phase change across a turning point is in agreement withthe change in the Maslov index [66] that often appears in semiclassical methods. We willshow and use this explicitly in the example provided in the following section.

3 Example: dynamics in a wide double-well potential
Let’s demonstrate how the formalism developed in the previous section can be applied
to describe the nonlinear dynamics of a particle using a relevant example. Inspired by arecent proposal to prepare a macroscopic quantum superposition of a levitated nanoparticle
via the dynamics in a wide nonharmonic potential [42], we will focus on the double-wellpotential given by:
1
X4
Vdw (X) = mω 2 −X 2 +
2
2D2

!

.

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(54)

12

0.5

4

0

2

-0.5

0
0

0.4

0.8

1.2

1

102
101
0

0

0

103

2

1

100

2

103

109

101

106

10-1

103

10-3

0

1

2

0

1

2

100
0

1

2

Figure 1: Relevant functions for our analytical approach computed for the parameters in Table 1. Panel(a) shows the classical trajectory rc (τ ) in units of d. Panels (b),(c),(d),(e) and (f) show, respectively,α(τ ) (cf. Eq. (19)), η(τ ) (Eq. (26)), φ(τ ) (Eq. (27)), βn (τ ) (Eq. (30)), and σb (τ )Ω/Γu both for theupper bounded (Eq. (62)) and the fluctuating decoherence (Eq. (39)) with Γloc = 0 and 25A1 = 2A2 d2such that both terms in Eq. (61) contribute equally to Γu . The polygons in panel (a) indicate relevant
instances of time, namely times when α(τ ) becomes zero (τ = τ3 ), times when η(τ ) is maximum(τ = τd ) and the time halfway through the classical trajectory (τ = τm ). These instances of time
appear as grid lines in the rest of the plots.

This potential is characterized by the frequency ω and the distance scale D. To adapt theresults of the previous section, we will employ the dimensionless potential:
1x4
Udw (x) =−x2 + 2
2
2d

!

,

(55)

where D = dXΩ . We will consider the wide-potential regime defined by d ≫ 1. We
consider the initial condition given by a thermal state of a harmonic potential of frequencyΩ centered at x = 0 with phonon mean number occupation n̄ and mean position andmomentum values denoted by µ(0) = (⟨x̂⟩(0), ⟨p̂⟩(0))T = (xs , 0)T , where d ≫ xs ≫ 1, suchthat the dynamics is constrained to the right side of the double-well potential. Whenevernecessary, we will utilize the parameters listed in Table 1, which correspond to the size L
dynamics defined in [42].
To transition to the classical centroid frame, we need to evaluate the classical trajectoryrc (τ ) = (xc (τ ), pc (τ ))T with rc (0) = µ(0). This classical trajectory can be solved analytically [67, 68] using the Jacobi elliptic functions [69]. In Fig. 1a, we display the classicaltrajectory, which is periodic and has an avocado shape. The analytical solution allows us
to obtain analytical expressions for the timescales governing the classical dynamics, which
are solely functions of xs /d. We will be particularly interested in the time scale√ !
4 2dτm = log+ O[(xs /d)2 ],(56)
xs
which corresponds to the turning point of the classical trajectory of the centroid. It alsoParameter:
Value:

ω/Ω
10−2

d
104

xs /d
10−1

n̄
0

Table 1: Parameters used in the double-well potential example.

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

13

corresponds to the half-period of the evolution, andpthe time at which the distance fromthe origin reaches its maximum, given by xc (τm ) = 2d2 − x2s .
In the classical centroid frame, the effective potential from Eq. (16) becomes timedependent. The harmonic part of this time-dependent effective potential, given by Eq. (19),has a dimensionless spring constant:xc (τ )α(τ ) = −1 + 3
d


2

,

(57)

which we illustrate in Fig. 1(b). Initially, this effective harmonic
potential is inverted√(α(τ ) < 0) until the time τ3 , which is defined by xc (τ3 ) = d/ 3. It then transforms into aharmonic potential (α(τ ) > 0) until time 2τm − τ3 , where it reverts to being inverted. Thissequence of inverted, harmonic, and inverted potentials echoes the dynamics leveraged in
the loop protocol proposed in [70].
Using the effective Harmonic potential, we can transition to the Gaussian frame. In
this frame, we can compute both the delocalization length η(τ ) and angle φ(τ ), which are
depicted in Fig.1(c) and Fig.1(d), respectively. The state initially expands (i.e., squeezes)
until it reaches a maximum at τ = τd , which for large Ω/ω is approximately given by1 Ω xsη⋆ = √
.
2ω d

(58)

The state then compresses (i.e., anti-squeezes) until the turning point τm . From here, thedynamics is mirrored until 2τm . Regarding the angle φ(τ ), we observe the anticipatedbehavior discussed in Sec. 2.5, i.e., it remains approximately constant (after the initialvariation) until the first turning point, at which point it undergoes an approximated πshift.
In the Gaussian frame, the strength of the nonlinear dynamics induced by the nonharmonicities in the potential is parameterized by the βn (τ ) (n > 2) variables, as given in
Eq. (30). For the double-well potential, we haveωβ3 (τ ) =3xc (τ )η(τ )3 ,
(59)
Ωd2
ω
β4 (τ ) =η(τ )4 .
(60)Ωd2These variables are illustrated in Fig.1(e). They demonstrate how the strength of these
nonlinearities is determined by the amount of squeezing η(τ ), which can be comparedwith Fig.1(c). Furthermore, it is evident that the cubic nonharmonicity is larger than thequartic nonharmonicity in the small-fluctuations regime, denoted by η(τ ) ≪ xc (τ ).
Regarding the impact of decoherence, the key variable within the linearized-decoherence
approximation, which is valid within the regime η(τ ) ≪ xc (τ ), is the blurring distanceσb (τ )2 , as shown in Eq. (39). The blurring distance depends on the effective timedependent decoherence rate Γeff (τ ) = Γloc + Γfluc (τ ), which includes the contribution fromthe potential’s √fluctuations as given in Eq. (40). For the case of the double-well potential
where xc (τ ) < 2d, Γfluc (τ ) is upper bounded by the time-independent decoherence rateΓfluc (τ )Γuπω 4 
2
<
≡
25A
+
2A
d
.
(61)
1
2
Ω
Ω
2Ω3
This allows us to conveniently bound the time-dependent decoherence rate by a timeindependent rate, namely Γeff (τ ) < Γ ≡ Γloc + Γu , leading to an upper-bound blurring
distance given by4Γσb (τ ) <ω2Z τdτ ′ η(τ ′ )2 .

(62)

0

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

14

10-1

100

10-3

10-2

10-5

10-4

10-7

0

0.25

0.5

0.75

1

10-6

0

0.25

0.5

0.75

1

Figure 2: Quantities related to the validity of the constant-angle approximation for dynamics in a
double-well with the parameters in Table 1. Panel (a) shows χ(τ ; ϕ), and panel (b) shows the relative
errors ϵ1 (τ ) (see Eq. (51)) and ϵ2 (τ ) (see Eq. (52)).

In Fig.1(f), we plot both the original blurring distance (dashed line) and the upper bound(solid line), showing that in the relevant time scale τ ≈ τm the upper bound is a factor of
10 larger than the original blurring distance.
To implement the constant-angle approximation, we need to set two angles: one, ϕ, forτ ∈ [0, τm ], and a second one, ϕ2 , for τ ∈ (τm , 2τm ]. As discussed in Sec. 2.5, and observedin Fig.1d, ϕ − ϕ2 = π − δ, where 0 < δ ≪ 1 necessitates fine tuning. We choose ϕ = φ(τd );that is, at the time when dφ(τ )/dτ is at its first minimum (recall Eq. (50) and that η(τ )is at its maximum at τ = τd ), and we select δ so that the state at τ = 2τm is closest to
the one numerically calculated. For the parameters given in Table 1, this corresponds to
ϕ/π = φ(τd )/π ≈ 0.499969 and δ/π = 7.63944 × 10−6 .
We are now in a position to analyze how well our analytical approach describes thenumerically exact nonlinear dynamics of this particular example.

3.1 Comparison of Analytical Approach with Numerically Exact Results
Following the discussion in Sec. 2.5, we first display in Fig. 2(a) the function χ(τ ; ϕ), as
defined in Eq. (49), evaluated for τ ∈ [0, τm ]. The fact that this quantity is small is
an indicator of the validity of the constant-angle approximation. More importantly, inFig. 2(b), we demonstrate the relative errors ϵn (τ ) for the first (n = 1) and second (n = 2)
moments, as defined in Eq. (51) and Eq. (52), for the case of dynamics without decoherence,
which is the most sensitive case. The numerical calculations, performed using the splitoperator method [43], reveal relative errors well below one percent, with the exceptionof the second moments around the turning point τ = τm . This is consistent with thefact that the constant-angle approximation is less accurate around a turning point. Thisis also evidenced in Fig. 2(a), where the increase of χ(τ ; ϕ) towards the turning point isobservable. These results suggest that our analytical approach should yield an excellent
approximation of nonlinear open dynamics.
To demonstrate this explicitly, we plot the Wigner function in the centroid frame at
the six specific instances of time indicated in Fig. 1a, as shown in Fig. 3. We present boththe numerically exact results W (r + rc (τ ), τ ) (which require several hours of computation)and the results using our analytical approach WnG (r + rc (τ ), τ ) (which only require a
few minutes). The agreement between the two is remarkably excellent. Furthermore, thecubic-phase states [64, 71–74], generated during the times τ ∈ [0, τm ], can be obtained
analytically. Indeed, one can show, as detailed in App. D.3, that the state in both thecentroid and Gaussian frames, in the absence of decoherence during this time scale, canAccepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

15

10

0

1
0.5

-10
0

10

-0.5
0

-1

-10
-10

0

10

-10

3

0

10

3

-10

3

0

10

3

-15 -10

-5

0

5

-10

3

0

10

3

-10

3

0

10

3

-10

0

10

Figure 3: Wigner function of the state of a particle evolving in a double-well potential with parameters
given in Table 1 at different instances of time. These instances of time are indicated by polygons
and correspond to the times indicated in Fig. 1(a). The first row shows the numerically exact Wigner
function W (r + rc (τ ), τ ) obtained using a numerically exact method whereas the second row shows
the approximated Wigner function WnG (r + rc (τ ), τ ) obtained using our analytical approach.
0.08
0.06
0.04
0.02
0

-14

-12

-10

-8

-6

-4

-2

0

2

Figure 4: Position probability distribution at time τm for a state evolving in a double-well potential forthe parameters in Table 1 and for different values of Γ. The lines are computed using the analyticalmethod described in this paper, whereas dots correspond to a numerically exact computation using
Q-Xpanse [11].

be expressed as follows:
Wa(G) (r, τ ) = √κ3 (τ )x2ϕ + pϕ
6κ3 (τ )pϕ + 1
1
1
exp
Ai
+
2
1/3
12|κ3 (τ )||κ3 (τ )|4|κ3 (τ )|4/3
2π|κ3 (τ )|1/3!

. (63)

Here, Ai(z) denotes the Airy function [75]. For simplification, we have omitted the contribution of κ4 (τ ) in Eq. (63), which is negligible for τ ≤ τm . The complete expression canbe found in App. D.3. Equation (63) indeed represents the Wigner function of a cubicphase state [64, 71–74]. Interestingly, the state at τ = 2τm exhibits a pronounced quartic
character. This might seem surprising given that, as we mentioned before, κ3 (τ ) ≫ κ4 (τ ).
However, upon inspecting the map MnG (τ ), one realizes that the two contributions of κ3 (τ )and κ4 (τ ), corresponding to the first half-orbit (0 < τ ≤ τm ) and the second half-orbit(τm < τ ≤ 2τm ), would exactly cancel out or add up for a change in the angle of precisely
π (i.e., δ = 0). When δ is small but nonzero, the contribution from κ4 (τ ) still accumulates,
while the difference of the cubic terms, to first order, contributes another quarticlike term,
thereby giving the state at 2τm a quartic character.
Finally, an interesting feature of the cubic-phase state generated at τ = τm is thatit
produces an interference pattern in the position probability distribution P(x, τ ) =R
dpW (r, τ ). This is depicted in Fig. 4 for different values of Γ/Ω, using both the numerically exact Q-Xpanse method [11] (dots) and our analytical approach (solid line).
This demonstrates not only excellent qualitative agreement but also quantitative, parAccepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

16

ticularly for the first and most relevant interference fringes. With our analytical approach, it canbe shown that the interference pattern in the presence of decoherence,R
Pa (x, τ ) = dpWa (r, τ ), can be expressed as a Gaussian convolution of the interference
pattern in the absence of decoherence [52], Pcoh,a (x, τ ), as follows:1

Pa (x, τ ) = q
2πCb,xx (τ )

Z

"

#

(x − x′ )2
dx exp −Pcoh,a (x′ , τ ).
2Cb,xx (τ )′(64)Here, Cb,xx (τ ) is the xx component of the blurring covariance matrix given in Eq. (45),andZ
Pcoh,a (x, τ ) ≡ dpMC (τ )MG (τ )MnG (τ )W (G) (r, 0),(65)
corresponds to the probability distribution in the absence of decoherence, which is explicitly
provided in App. D.4. In addition, the interference pattern in the absence of decoherence
can be derived from Eq. (63) and is given by:x − xc (τ )Pcoh,a (x, τ ) ∝ exp −2mxp (τ )κ3 (τ )!

!2x − xc (τ )Ai
+ z(τ )
mxp (τ )(4κ3 (τ ))1/3

.

(66)

In this equation, mxp (τ ) = Sxx (τ ) sin(ϕ)−Sxp (τ ) cos(ϕ), and z(τ ) ∈ C is a time-dependentfunction that is given in App. D.4. This analytical expression allows us to derive the scalingof the fringe separation, xf , which is the distance between the largest interference peak in
Pcoh,a (x, τm ) and its second one. As shown in App. D.4, the fringe separation scales as
xf ∝ |κ3 (τ )|1/3 , or in terms of the physical parameters as xf ∝ (Ω/ω)2/3 (1/d)1/3 , for a fixed
xs /d.

4 Conclusions
We have provided an analytical treatment of the Wigner function dynamics for a continuousvariable degree of freedom, such as the position and momentum of a massive particle, in anonharmonic potential and in the presence of decoherence caused by, among other things,
fluctuations of the nonharmonic potential. Our analytical treatment has demonstrated
that it can yield an expression of the time-evolved Wigner function that approximates theexact dynamics very accurately, particularly in the case of wide nonharmonic potentials
and small fluctuations.
Our analytical method specifically entails performing an exact, suitable reformulationof the open nonlinear dynamical problem via two frame transformations: the classical centroid and the Gaussian frame transformations. In this transformed frame, two approximations can be applied: the key constant-angle approximation and the linearized-decoherence
approximation. These approximations facilitate the integration of the dynamical problem
and are argued to provide excellent approximations for dynamics that allow the state toexpand over scales many orders of magnitude larger than the initial spatial expansion, provided these fluctuations are smaller than the available phase space that can be explored.
In other words, the fluctuations of the state should not occupy the entire available phasespace surface. These requirements define the concept of wide nonharmonic potentials andsmall fluctuations.
We have tested our analytical method using an example of a massive particle evolving
in a wide double-well potential. The method shows excellent agreement with numericalsimulations, not just qualitatively, but also quantitatively. In the wide regime, especially in

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

17

the presence of decoherence, the numerical calculations can be challenging. For these calculations, we utilized a numerical tool, Q-Xpanse, that we recently developed [11]. Thesenumerical calculations often require several hours of computation, while our analytical
approach reproduces the results with simple evaluations requiring only a few minutes ofcalculation. This example is particularly relevant as it has been recently proposed for
preparing macroscopic quantum states of a levitated nanoparticle via the generated nonlinear quantum dynamics [42]. The combination of the numerical tool and the analyticalmethod has enabled us to design, optimize, and understand this experimental proposal.
Beyond its utility in modeling experiments, our analytical method also paves the way
for interesting research questions that warrant further investigation. Specifically, it would
be intriguing to better characterize the dynamical problems in which our method proveseffective. It could potentially be utilized to define a class of semiclassical quantum dynamical problems (see e.g., [32, 35]). Furthermore, one could concentrate on the classicalregime, either by taking the limit as ℏ → 0 and studying the efficacy of this method in understanding classical nonlinear dynamics, or by augmenting the strength of the generatorsof decoherence. In the latter case, it would also be interesting to focus on different typesof decoherence, including non-Gaussian decoherence, and apply our method to understandand characterize their effect on the dynamics.
In conclusion, we hope that our analytical approach to understanding open quantum
nonlinear dynamics and its classical limit will significantly contribute to this captivating,albeit long-standing, topic.

Acknowledgments
ARC and MRL contributed equally to this work. We thank the Q-Xtreme synergy groupfor fruitful discussions. This research was supported by the European Union’s Horizon 2020
research and innovation programme under grant agreement No. [863132] (IQLev) and by
the European Research Council (ERC) under the grant Agreement No. [951234] (Q-Xtreme
ERC-2020-SyG). ARC acknowledges funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant agreement
No. [101103589] (DecoXtreme, HORIZON-MSCA-2022-PF-01-01). PTG was partially
supported by the Foundation for Polish Science (FNP).

A Explicit expression of the decoherence rates
In the main body of the text, we have reformulated the dissipator due to the fluctuation
potential (cf. Eq. (5)) as a double sum with coefficients denoted by Γnm . To achieve this,
we performed a Taylor expansion of the potential, V (X), as well as its derivative in Eq. (4).
This led us to the explicit form of the equation:"#

2π X m+n∂ m+1 V∂ n+1 V∂mV∂nV
Γnm = 2 Ω
A1 XΩ2(0)(0) + A2
(0)(0) +δn,1 δm,1 Γloc . (67)
m+1
n+1
m
ℏ m!n!
∂X
∂X
∂X
∂X n
These coefficients Γnm are particularly convenient for expressing the decoherence superoperator within the Hilbert space. However, their translation to the Wigner equation isconsiderably different, necessitating the introduction of a new set of coefficients denoted as
cnmk in Eq. (9). These can be derived by applying the Wigner transform to the decoherence

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

18

superoperator D[·], resulting in the following equation:hcnmk = (−i)k 1 + (−)k



iΓ

nm 

2ω

k
X

(−)q

q=0

n
k

!

!

!

m
m+n 
−
.
k−q
k

(68)

B Details on the constant-angle approximation
Here, we discuss in detail the constant-angle approximation. The dynamics in a wide
potential is such that the coherent expansion η(τ ) ≫ 1, and therefore the angle φ(τ ) is
approximately constant except close to compression times. Hence, one can approximateφ(τ ) as a piecewise time-dependent function. Let τm be the time at which the only turning
point between 0 and τf occurs. Then, we first use the exact factorization property of theevolution map Z τf

M+ (τf ) = exp+ g

0

 Z τ
f



L(τ1 )dτ1 = exp+ gτm Z τm

L(τ1 )dτ1 exp+ g

L(τ1 )dτ1 ,

(69)

0

for a generic gL(τ ). Now in each time interval, [0, τm ] and [τm , τf ] the angle is approximately
constant. Namely, our analytic calculation is based on the replacement of M+ (τf ) by
Z τ

Z τm



f

L(τ1 )dτ1 exp

M (τf ) = expτm

L(τ1 )dτ1 ,

(70)

0

together with the approximation(for τ ∈ [0, τm ),.
for τ > τm

ϕ1
ϕ2

φ(τ ) ≈

(71)

for two constant angles ϕ1 and ϕ2 . Note that we have reduced the problem to the scenario
in which there are no compression times in the integration region. In what follows, let us
assume that there are no compression times from between 0 and τ .
The constant angle approximation described in the main text has the effect of replacing
a time-ordered exponential, which is difficult to calculate, by an easier simple exponential
map. Mathematically, the reason why this is possible can be justified as follows. Consider
the expansion of the time-ordered exponential map Z τ

M+ (τ ) = exp+ g

dτ1 L(τ1 ) =

0

∞
X

gn

Z τ

dτ1 · · ·dτn L(τ1 ) · · · L(τn ),

(72)

0

0

n=0

Z τn−1

with the time-dependent generator gL(τ ). The corresponding “standard” exponential map
is instead Z τM (τ ) = exp gdτ1 L(τ1 ) =

0

Z
∞
X
gn τ
n=0

n!

0

dτ1 · · ·Z τdτn L(τ1 ) · · · L(τn ).

(73)

0

Note that if L(τ ) commutes with itself at different times, both expressions coincide after
redefining the integration region. In the general case, the difference between the two maps
is given by
M+ (τ ) − M (τ ) =

g2
2

Z τ

Z τ1

dτ1
0

dτ2 [L(τ1 ), L(τ2 )] + O(g 3 ),

(74)

0

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

19

which explicitly shows that the first correction depends on the out-of-time commutator of
the generator.
(G)(G)
In our case of interest L(τ ) = LnG (τ )+Ld (τ ), where L(τ ) includes only multiplication
with respect to xφ(τ ) and (multiple) derivatives with respect to pφ(τ ) . Hence, the term∂x(τ2 ) = sin(φ(τ1 ) − φ(τ2 )),∂pφ(τ1 )(75)

controls the strength of the commutator. The explicit calculation is lengthy and therefore
we do not include it here. However, let us show how to proceed by considering only the“classical” term defined as
L(G)c (τ ) ≡X

βn+1 (τ )xnφ(τ )

n≥2

∂
,
∂pφ(τ )

(76)

since the rest are analogous. We find that the first-order correction in Eq. (74) is given by
1
2

Z τ1

Z τ

dτ1
0

0

(G)
dτ2 [L(G)c (τ1 ), Lc (τ2 )] =

1
2

Z τ1

Z τ

dτ1

dτ2 A(τ1 , τ2 ) sin(φ(τ1 ) − φ(τ2 )), (77)00where we have defined the differential operator"

A(τ1 , τ2 ) =X

n−1 m−1xmxφ(τ1 )βn+1 (τ1 )βm+1 (τ2 )xφ(τ
1 ) φ(τ2 )

n,m≥2

∂
∂pφ(τ2 )+ nxφ(τ2 )∂∂pφ(τ1 )#. (78)In the regime of large expansions η(τ ) ≫ 1 allowed by wide potentials, the angle φ(τ ) is
approximately constant (see Eq. (50)). Hence, replacing the time-dependent angle φ(τ )by a constant ϕ corresponds to replacing the time-ordered exponential by a standard
exponential map, which in turn allows to continue with the analytical calculation.
The approximation still depends on the value of the angle ϕ that one chooses. We decide
to take ϕ as the angle φ(τ⋆ ) evaluated at the time of largest expansion τ⋆ = argmaxτ η(τ ),since this corresponds to the time at which the derivative of φ(τ ) is the smallest.

C An alternative linearized-decoherence approximation requirement
The small fluctuations condition xc (τ ) ≫ η(τ ) ensures that decoherence differential op(G)erator Ld (τ ) can be linearized (see Eq. (41)). However, a drawback of this conditionis that xc (τ ) is defined relative to an arbitrary origin. While the coefficients cnmk (andconsequently Γnm ) appropriately account for this arbitrariness, it can potentially obscure
the more general condition for the linearized-decoherence approximation to hold. Here,
we provide a more intricate derivation of the operator Λd (τ ) in Eq. (36), but that has theadvantage of eliminating the arbitrariness associated with the choice of origin.
Let f (X̂/XΩ ) be a smooth function of the position operator X̂, and let us compute
the Wigner transform∂∂f (X̂/XΩ )ρ̂f (X̂/XΩ ) 7→ f x + i
f x−i
W (r).
∂p
∂p








(79)

Using this relation, a double commutator term reads∂∂
f x + is′−[f (X̂), [f (X̂), ·]] 7→ Ξ· ≡(−ss )f x + is∂p∂ps,s′ =±
X

′







Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.



·.

(80)

20

In the dissipator due to fast fluctuating forces, we have two such terms: f (x) ∝ (∂U/∂x)(x)and f (x) ∝ U (x). The transformation to the classical centroid and Gaussian frames
amounts to transforming the expression Ξ intoX

Ξ(G) (τ ) =(−ss′ )f xc (τ ) + η(τ )Zs (φ(τ )) f xc (τ ) + η(τ )Zs′ (η(τ )) · .




(81)

s,s′ =±1

with the short-hand notation∂
xϕ + is∂pϕ

Zs (ϕ) =

!

.

(82)

One can now expand Eq. (81) around xc (τ ) to obtain
Ξ(G) (τ ) =

∞  n
X
∂ f

(xc (τ ))η(τ )n
∂xn

n,m=0

 m
∂ f

(xc (τ ))η(τ )m
∂xm

 X

(−ss′ )Zs (ϕ(t))n Zs′ (ϕ(t))m .

ss′ =±(83)
If we label the terms in the series as (n, m), the first nonvanishing term is (1, 1). The
linearized decoherence approximation corresponds to keeping the term (1, 1) and assumingthat the rest are subdominant. In essence, we require that the terms of the sequence {an }with∂nfan = η(τ )n n (xc (τ )),(84)∂xdecreases for increasing n, where f (x) is proportional to either the potential U (x) or its
derivative. This is expected in the limit of wide potentials. Indeed, let d be the length-scaleof the potential, then we haveη(τ ) n
η(τ )
(x
(t))
∼
f (xc (τ )).
(85)
c
∂xn
d
Therefore, the term with the least number of factors η(τ )/d is the term (1, 1), and we expect
it to dominate the series. Note that the above logic does not exploit that xc (τ ) ≫ η(τ ),but rather that η(τ ) is small as compared to the derivatives of the potential U (x); that is,that the potential is wide.
n∂

nf





D Derivations using the analytical expression of the Wigner function
This appendix is devoted to providing additional details regarding the derivations of analytical expressions for the approximated first and second moments, along with the approximated Wigner function, using the constant-angle approximation; that is, using our main
result in Eq. (43). In the derivation, we use the properties
ZZ

drf (r){G[C](r) ⋆ g(r)]} =Z

drf (r)g(T −1 r) =dr{G[C](r) ⋆ f (r)}g(r),

(86)

Z

drf (T r)g(r)JT (r),(87)

where ⋆ denotes the convolution,Zf (r) ⋆ g(r) =dr ′ f (r − r ′ )g(r ′ ),(88)

and the Jacobian JT (r) is the determinant of the matrix∂(T r)=∂(r)∂x (T r)x ∂p (T r)x∂x (T r)p ∂p (T r)p

!

(89)

associated to the transformation T .
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

21

D.1 First moments
In the Wigner representation, expectation values of symmetrically-ordered (also Weylordered) operators are computed as phase space integrals. Using our main result in
Eq. (43), we find thatµa (τ ) ≡ ⟨r̂⟩a (τ ) =Z

drrDa (τ )MC (τ )MG (τ )MnG (τ )W (G) (r, 0).

(90)

Using the property in Eq. (86) together with the convolution representation of the decoherence map in Eq. (46), one can show that the decoherence map does not contribute to
the evaluation of the first moments. Moreover, using the property in Eq. (87) for the maps
MC (τ ) and MG (τ ), and noting that both transformations have unit Jacobian, one arrivestoZµa (τ ) = rc (τ ) + S(τ )drMnG (τ )W (G) (r, 0).

(91)

Finally, we note that, in our approximation, the quantum terms included in MnG (τ ) do
not contribute to the calculation of the first or second moments. This is a consequenceof the fact that the quantum terms include derivatives of order higher than three, which
after integration by parts would give no contribution to the integral. Hence, the action of
the map MnG (τ ) amounts to a nonlinear transformation
T r = r − R(ϕ)T ep

Xn−1
κn (τ )(eTx R(ϕ)r)

(92)

n≥3

which also has a unit Jacobian. Hence, combining our previous results and taking the
average with the state Wg (r, 0) we arrive at the result quoted in Eq. (47).

D.2 Covariance matrix
The covariance matrix is computed as the average
Ca (τ ) ≡⟨(r̂ − µ(τ ))(r̂ − µ(τ ))T ⟩(τ )Z=dr(r − µ(τ ))(r − µ(τ ))T Da (τ )MC (τ )MG (τ )MnG (τ )W (G) (r, 0).

(93)

Taking advantage of the property in Eq. (86), we see that the action of the decoherence
map is to increase the value of the covariance matrix by Cb (τ ), as compared to the casewithout decoherence. Namely, one arrives atZ

Ca (τ ) = Cb (τ ) +dr(r − µ(τ ))(r − µ(τ ))T MC (τ )MG (τ )MnG (τ )W (G) (r, 0).

(94)

Moreover, similarly to the case of the first moments, taking advantage of Eq. (87) andEq. (92), one can manipulate the result into the form displayed in Eq. (48).

D.3 State before the turning point in the double-well caseHere, we give more details on how to obtain an analytical expression for the Wigner functionin the centroid and Gaussian frame under the constant-angle approximation. To this end,
we assume a thermal initial state Wg (r, 0), which we decompose into Fourier modes as
W (G) (r, 0) =hi
exp[−x2ϕ /(2(2n̄ + 1))] Z dkϕpexp −(2n̄ + 1)kϕ2 /2 + ipϕ kϕ ,2π
2π(2n̄ + 1)

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(95)

22

Applying the evolution map given by the exponential of Λ(τ ) in Eq. (34) to the state in
Eq. (95) leads to an integral expression of the form
1
2π

Z

"

k3
k2
dk exp i c3 + ic2 + c1 k
32!#

c32 + 6c3 c2 c1
1
exp=
12|c3 |2
|c3 |1/3!

c2 + 4c3 c1
Ai 2
4|c3 |4/3

!

,
(96)

where c3 ≡ κ3 (τ ) + 3κ4 (τ )xϕ , c2 ≡ 2n̄ + 1 + σb (τ )2 , and c1 ≡ pϕ + κ3 (τ )x2ϕ + κ4 (τ )x3ϕ andwe have introduced the Airy function
1
Ai(z) ≡
2π

"

!#

u3
du exp i
+ zu
3
−∞

Z ∞

,

(97)

for z ∈ C. Multiplying by the remaining Gaussian in xϕ , one can simplify the expression
into the final result
1
Wa(G) (r, τ ) = p2π(2n̄ + 1)|κ3 (τ ) + 3κ4 (τ )xϕ |1/3
(2n̄ + 1 + σb (τ )2 )[6(κ3 (τ ) + 3κ4 (τ )xϕ )(pϕ + κ3 (τ )x2ϕ + κ4 (τ )x3ϕ )]× exp
12|κ3 (τ ) + 3κ4 (τ )xϕ |2"x2ϕ(2n̄ + 1 + σb (τ )2 )3
× exp
exp
−
12|κ3 (τ ) + 3κ4 (τ )xϕ |2
2σ 2
"

× Ai

#

(κ3 (τ )x2ϕ + κ4 (τ )x3ϕ + pϕ )|κ3 (τ ) + 3κ4 (τ )xϕ |1/3

"

#

#

(2n̄ + 1 + σb (τ )2 )2+
4|κ3 (τ ) + 3κ4 (τ )xϕ |4/3

!

.

(98)

If one sets into the above expression κ4 (τ ) 7→ 0, σb (τ ) 7→ 0, and n̄ 7→ 0 it yields thesimplified result
Wa(G) (r, τ ) = √κ3 (τ )x2ϕ + pϕ
6κ3 (τ )pϕ + 1
1
1
Ai
+
exp
2
1/3
1/3
12|κ3 (τ )||κ3 (τ )|4|κ3 (τ )|4/3
2π|κ3 (τ )|!

, (99)

which we show in the main text. Finally, one can move back to the centroid frame by applying MG (τ ), and to the original frame by applying MC (τ ). It follows that the approximatedWigner function yields
Wa (r, τ ) = Wa(G) (S(τ )−1 (r − rc (τ )), τ ).

(100)

D.4 Interference pattern at the turning point in the double-well caseIn this Appendix, we compute the coherent interference pattern of the probability distribution Pcoh,a (x, τ ) within our approximation. To obtain the desired probability distribution,
we first rewrite xϕ and pϕ in terms of x and p using the rotation matrix R(ϕ). Then, weapply the map MG (τ ) to move to the classical trajectory frame, which mixes the coordinates x and p according to the propagator S(τ )−1 . If we denote f (rϕ ) the function of theright-hand-side of Eq. (99), the Wigner function in the classical centroid frame is f (m(τ )r),where m(τ ) ≡ R(ϕ)S(τ )−1 , with components mij (τ ) for i, j = x, p. For conciseness of thenotation, hereafter we denote mij (τ ) simply as mij .
In order to obtain the probability distribution, we have to perform the integral overthe variable p, which can be done using the following property of the Airy function (seeCh.3 of [75])ZdyAi(y 2 + y0 ) exp(iky) = 22/3 πAi[2−2/3 (y0 + k)]Ai[2−2/3 (y0 − k)].

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

(101)

23

In order to use the property in Eq. (101), we have to rewrite the argument of the Airyfunction in Eq. (63) without a linear term. To this end, we define1y=κ3 (τ )1/6

κ3 (τ )

1/2

2κ3 (τ )mxp mxx x + mpp
mxp p +
2mxp κ3 (τ )1/2

!

,

(102)

m2xp − m2pp − 4mxp κ3 (τ )x1y0 (x) =,4m2xp κ3 (τ )
κ3 (τ )1/3
mpp.
k = −i2mxp κ3 (τ )4/3
!

(103)
(104)

Noting that mxx mpp − mxp mpx = 1, and using y as an integration variable, we can use
formula Eq. (101) to compute(C)Pcoh,a (x, τ ) =m2xp − 3m2pp − 6mxp κ3 (τ )x
π/2
exp
12m2xp κ3 (τ )2|mxp |(κ3 (τ )/2)2/3
p

"

"

m2pp
1
x
+
1− 2
−
mxp 4κ3 (τ )
mxx

1
× Ai
(4κ3 (τ ))1/3

!

#

mpp
+i
2mxp κ3 (τ )

!# 2

.
(105)

Introducing back the expressions for mxx , mxp , mpx and mpp , and shifting x 7→ x − xc (τ )yields the probability Pcoh,a (x, τ ). The effect of decoherence can be included with Eq. (64).
Let xf be the distance between the first and second maximum in the fringe pattern
defined by Eq. (105). At τm , the dependence on the position variable x is scaled withxf ∼ |κ3 (τm )|1/3 mxp .

(106)

with mxp ≈ Sxx (τm ), which is only a function of xs /d. Moreover, one has thatω1
κ3 (τm ) = 3
Ωd

Z τm

′
′ 3 xc (τ )′dτ η(τ )
0

d

ω1 3
=3
η
Ωd ⋆

Z τm
0

dτ ′

η(τ ′ )3 xc (τ ′ )
.
η⋆3
d

(107)

For sufficiently large Ω/ω the integral in Eq. (107) becomes, in very good approximation,
also a function of only xs /d. In that regime, we find that the fringe separation scales with


xf ∼

Ω
ω

2/3  1/3
1df (xs /d),(108)

where f (xs /d) is an unspecified function of only xs /d.

References[1] E. Schrödinger, Phys. Rev. 28, 1049 (1926).
[2] M. Born, Z. Phys. 38, 803 (1926).
[3] W. Heisenberg, Z. Phys. 43, 172 (1927).
[4] E. Wigner, Phys. Rev. 40, 749 (1932).
[5] W. P. Schleich, Quantum optics in phase space (John Wiley & Sons, Ltd, 2001).
[6] W. B. Case, Am. J. Phys. 76, 937 (2008).
[7] S. Habib, K. Shizume, and W. H. Zurek, Phys. Rev. Lett. 80, 4361 (1998).
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

24

[8] W. H. Zurek, Nature 412, 712 (2001).
[9] W. H. Zurek, “Decoherence and the transition from quantum to classical — revisited”, in Quantum decoherence: poincaré seminar 2005 (Birkhäuser Basel, Basel,2007), pp. 1–31.
[10] R. Cabrera, D. I. Bondar, K. Jacobs, and H. A. Rabitz, Phys. Rev. A 92, 042122(2015).
[11] M. Roda-Llordes, D. Candoli, P. T. Grochowski, A. Riera-Campeny, T. Agrenius,
J. J. García-Ripoll, C. Gonzalez-Ballestero, and O. Romero-Isart, Phys. Rev. Res.
6, 013262 (2024).
[12] S. Deléglise, I. Dotsenko, C. Sayrin, J. Bernu, M. Brune, J.-M. Raimond, and S.
Haroche, Nature 455, 510 (2008).
[13] M. Hofheinz, H. Wang, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D.
O’Connell, D. Sank, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature 459, 546(2009).
[14] G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, Nature 495, 205 (2013).
[15] D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano, and D. J.
Wineland, Phys. Rev. Lett. 77, 4281 (1996).
[16] C. Flühmann and J. P. Home, Phys. Rev. Lett. 125, 043602 (2020).
[17] C. Gonzalez-Ballestero, M. Aspelmeyer, L. Novotny, R. Quidant, and O. RomeroIsart, Science 374, eabg3027 (2021).
[18] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Science 367, 892 (2020).
[19] L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong,
N. Kiesel, A. Kugi, and M. Aspelmeyer, Nature 595, 373 (2021).
[20] F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Nature
595, 378 (2021).
[21] M. Kamba, R. Shimizu, and K. Aikawa, Opt. Express 30, 26716 (2022).
[22] A. Ranfagni, K. Børkje, F. Marino, and F. Marin, Phys. Rev. Res. 4, 033051 (2022).
[23] J. Piotrowski, D. Windey, J. Vijayan, C. Gonzalez-Ballestero, A. de los Ríos Sommer,
N. Meyer, R. Quidant, O. Romero-Isart, R. Reimann, and L. Novotny, Nat. Phys.
19, 1009 (2023).
[24] M. Kamba, R. Shimizu, and K. Aikawa, Nat. Commun. 14, 7943 (2023).
[25] A. Caldeira and A. Leggett, Physica A 121, 587 (1983).
[26] W. G. Unruh and W. H. Zurek, Phys. Rev. D 40, 1071 (1989).
[27] H. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, 2002).
[28] M. E. Gehm, K. M. O’Hara, T. A. Savard, and J. E. Thomas, Phys. Rev. A 58, 3914
(1998).
[29] S. Schneider and G. J. Milburn, Phys. Rev. A 59, 3766 (1999).
[30] C. Henkel, S. Pötting, and M. Wilkens, App. Phys. B 69, 379 (1999).
[31] M. Berry and K. E. Mount, Rep. Prog. Phys. 35, 315 (1972).
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

25

[32] M. V. Berry, Philos. Trans. Royal Soc. A 287, 237 (1977).
[33] J. B. Keller, SIAM Review 27, 485 (1985).
[34] M. Zworski, Semiclassical analysis, Vol. 138 (American Mathematical Society, 2022).
[35] E. J. Heller, J. Chem. Phys. 62, 1544 (1975).
[36] E. J. Heller, J. Chem. Phys. 65, 1289 (1976).
[37] E. J. Heller, J. Chem. Phys. 67, 3339 (1977).
[38] M. J. Davis and E. J. Heller, J. Chem. Phys. 80, 5036 (1984).
[39] R. Littlejohn, Phys. Rep. 138, 193 (1986).
[40] D. Huber, E. J. Heller, and R. G. Littlejohn, J. Chem. Phys. 89, 2003 (1988).
[41] E. J. Heller, The semiclassical way to dynamics and spectroscopy (Princeton University Press, 2018).
[42] M. Roda-Llordes, A. Riera-Campeny, D. Candoli, P. T. Grochowski, and O. RomeroIsart, Phys. Rev. Lett. 132, 023601 (2024).
[43] C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff,
J. Comput. Phys. 94, 59 (1991).
[44] L. Dania, D. S. Bykov, F. Goschin, M. Teller, A. Kassid, and T. E. Northup, Phys.
Rev. Lett. 132, 133602 (2024).
[45] E. Joos and H. D. Zeh, Z. Phys. B 59, 223 (1985).
[46] M. Schlosshauer, Decoherence and the quantum-to-classical transition (SpringerVerlag, Berlin, Heidelberg, 2007).
[47] O. Romero-Isart, Phys. Rev. A 84, 052121 (2011).
[48] B. L. Hu, J. P. Paz, and Y. Zhang, Phys. Rev. D 45, 2843 (1992).
[49] N. Van Kampen, Phys. Rep. 24, 171 (1976).
[50] V. Jain, J. Gieseler, C. Moritz, C. Dellago, R. Quidant, and L. Novotny, Phys. Rev.
Lett. 116, 243601 (2016).
[51] P. Maurer, C. Gonzalez-Ballestero, and O. Romero-Isart, Phys. Rev. A 108, 033714(2023).
[52] O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer,
and J. I. Cirac, Phys. Rev. Lett. 107, 020405 (2011).
[53] J. Bateman, S. Nimmrichter, K. Hornberger, and H. Ulbricht, Nat. Commun. 5, 4788
(2014).
[54] T. Agrenius, C. Gonzalez-Ballestero, P. Maurer, and O. Romero-Isart, Phys. Rev.
Lett. 130, 093601 (2023).
[55] M. R. Gallis and G. N. Fleming, Phys. Rev. A 42, 38 (1990).
[56] N. Van Kampen, Physica 74, 239 (1974).
[57] R. P. Feynman, Phys. Rev. 84, 108 (1951).
[58] N. Van Kampen, Physica 74, 215 (1974).
[59] D. W. Moore and R. Filip, Comm. Phys. 5, 128 (2022).
[60] C. A. Rosiek, M. Rossi, A. Schliesser, and A. S. Sørensen, arXiv: 2312.12986 (2023).
Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

26

[61] S. Mukamel, Principles of nonlinear optical spectroscopy (Oxford University Press,
New York, 1999).
[62] O. Romero-Isart, New J. Phys. 19, 123029 (2017).
[63] H. Pino, J. Prat-Camps, K. Sinha, B. P. Venkatesh, and O. Romero-Isart, Quantum
Sci. Technol. 3, 025001 (2018).
[64] L. Neumeier, M. A. Ciampini, O. Romero-Isart, M. Aspelmeyer, and N. Kiesel, Proc.
Natl. Acad. Sci. U.S.A. 121, e2306953121 (2024).
[65] T. Weiss and O. Romero-Isart, Phys. Rev. Res. 1, 033157 (2019).
[66] V. P. Maslov and M. V. Fedoriuk, Semi-classical approximation in quantum mechanics, Vol. 7 (D. Reidel publishing company, London, England, 1981).
[67] C. S. Hsu, Q. Appl. Math. 17, 393 (1960).
[68] A. J. Brizard, Eur. J. Phys. 30, 729 (2009).
[69] E. T. Whittaker and G. N. Watson, A course of modern analysis, 4th ed., Cambridge
Mathematical Library (Cambridge University Press, 1996).
[70] T. Weiss, M. Roda-Llordes, E. Torrontegui, M. Aspelmeyer, and O. Romero-Isart,
Phys. Rev. Lett. 127, 023601 (2021).
[71] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro,and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
[72] M. Brunelli and O. Houhou, Phys. Rev. A 100, 013831 (2019).
[73] V. Kala, R. Filip, and P. Marek, Opt. Express 30, 31456 (2022).
[74] A. A. Rakhubovsky and R. Filip, npj Quantum Inf. 7, 120 (2021).
[75] O. Vallée and M. Soares, Airy functions and applications to physics, 2nd ed. (Imperial
College Press, 2010).

Accepted in Quantum 2024-06-24, click title to verify. Published under CC-BY 4.0.

27