Stochastic Representation of the Quantum Quartic Oscillator

Gennaro Tucci1 a, Stefano De Nicola2 b, Sascha Wald3,4 c, and Andrea Gambassi5
d

1 Max Planck Institute for Dynamics and Self-Organization,

37077 Göttingen, Germany

2 ISTA, Am Campus 1, 3400 Klosterneuburg, Austria

3 Statistical Physics Group, Centre for Fluid and Complex Systems,

Coventry University, Coventry, England

4 $\mathbb{L}^{4}$ Collaboration & Doctoral College for the Statistical
Physics of Complex Systems, Leipzig-Lorraine-Lviv-Coventry, Europe

5 SISSA - International School for Advanced Studies and INFN, via Bonomea 265,

I – 34136 Trieste, Italia

a<EMAIL_ADDRESS>
b<EMAIL_ADDRESS>
c<EMAIL_ADDRESS>
d<EMAIL_ADDRESS>
## Abstract

Recent experimental advances have inspired the development of theoretical
tools to describe the non-equilibrium dynamics of quantum systems. Among them
an exact representation of quantum spin systems in terms of classical
stochastic processes has been proposed. Here we provide first steps towards
the extension of this stochastic approach to bosonic systems by considering
the one-dimensional quantum quartic oscillator. We show how to exactly
parameterize the time evolution of this prototypical model via the dynamics of
a set of classical variables. We interpret these variables as stochastic
processes, which allows us to propose a novel way to numerically simulate the
time evolution of the system. We benchmark our findings by considering
analytically solvable limits and providing alternative derivations of known
results.

###### Contents

  1. 1 Introduction
  2. 2 Stochastic Representation of the Quantum Quartic Oscillator
  3. 3 Physical Observables
    1. 3.1 Dynamics of a Gaussian Wave Packet
    2. 3.2 Exactly Solvable Cases
      1. 3.2.1 Harmonic case
      2. 3.2.2 Commuting limit
  4. 4 Stochastic interpretation and numerical benchmark
  5. 5 Perturbative Expansion
  6. 6 Semiclassical Approximation
    1. 6.1 Propagator
    2. 6.2 Partition Function
  7. 7 Summary and outlook
  8. A Stochastic integral and time-dependent quartic Hamiltonian
    1. A.1 Gaussian integral
    2. A.2 Time-dependent quartic Hamiltonian
  9. B Action of $\exp(w\,\hat{S}^{\pm,z})$ on a Gaussian wave packet
  10. C Time evolution of a Gaussian wave packet
  11. D Derivations of the Dyson Series
  12. E Semiclassical limit

##  1 Introduction

Recent advances in cold atom experiments have motivated great theoretical
interest in the non-equilibrium dynamics of isolated many-body quantum systems
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. With the notable exception of integrable
models [11], analytical insights into the time evolution of these systems are
scarce. This motivates the search for novel numerical and analytical tools to
analyse many-body quantum dynamics and to complement other widely used
numerical approaches [12, 13, 14]. Some recent works [15, 16, 17, 18]
considered an exact representation of the dynamics of many-body quantum spin
systems in terms of classical stochastic processes. This approach is based on
a series of exact transformations, through which an interacting many-body
quantum system can be exactly represented as non-interacting system under the
action of a set of classical stochastic fields. In this stochastic approach,
sometimes referred to as disentanglement formalism111This terminology
originates from the fact that the approach described here “disentangles” a
time-ordered exponential, featured in the time-evolution operator, into a
product of ordinary exponentials, see, e.g., Refs. [17, 19]. Hence, this term
does not actually point to a relation to quantum entanglement. Here we mostly
stick with the terminology “stochastic approach”., physical observables can be
expressed as averages of classical functions over realizations of suitably
constructed stochastic processes. This not only breaks down the complexity of
many-body quantum interactions to that of classical stochastic differential
equations but also allows one to bridge the gap between the quantum realm and
classical stochastic processes, for which powerful numerical tools have been
developed [20]. Beside this numerical application, the formalism can also be
used as a field-theoretical tool to develop analytic expansions for
observables [21]. The stochastic formalism has been mainly applied to many-
body quantum spin systems [15, 16, 17, 18, 22, 21, 23], where interactions are
mapped to classical stochastic fields. Here, we explore a generalization of
the stochastic approach to systems of interacting bosons. In particular we
shall study the time-evolution of a zero-dimensional bosonic system with a
non-linear interaction potential. These systems have attracted significant
attention for their rich dynamical behavior and relevant experimental
applications [24, 25, 26, 27, 28, 29]. In particular, we shall consider the
quantum quartic oscillator as a paradigmatic, non-linear bosonic quantum
system [30, 31, 32, 33, 34, 35, 36, 37]. While any system with harmonic
interactions can be mapped to non-interacting bosons, this is no longer the
case for anharmonic potentials, whose simplest representative is arguably the
quartic one.

In spite of its apparent simplicity, the quantum quartic oscillator is
notoriously hard to study and involves subtle conceptual mathematical issues,
even in its equilibrium formulation: notably, the divergence of the
perturbation theory for the ground state energy in powers of the quartic
coupling $\lambda$ and the fact that the energy levels have an infinite number
of branch points in the complex $\lambda$-plane around the origin $\lambda=0$,
i.e., in the harmonic limit. From a physical perspective, these singularities
are caused by level crossings in the eigenenergies [30, 31], leading to a non-
trivial perturbation theory of the energy spectrum [36, 38]. The quantum
quartic oscillator is also of practical relevance, as it can be used to
approximate the low-energy behavior of more general, real-world quantum
systems [39, 40, 41].

The quantum quartic oscillator has been investigated by approximation
techniques, such as the semiclassical evaluation of its propagator [32, 33,
34, 37, 42] and, more recently, by looking at the evolution of relevant time-
dependent observables [43, 29]. Interestingly, the exact form of the wave
functions of the quartic oscillator has been reported only recently [44] but
there is still no analytic expression for the corresponding quantized energy
levels. Here, we show how the stochastic approach can be used to represent the
dynamics of a quantum quartic oscillator in terms of ensembles of
stochastically evolving harmonic oscillators. The representation we introduce
is formally exact and can be used to develop analytical approximations or to
evaluate numerically expectation values, providing a different viewpoint as
well as a practical alternative to existing techniques.

This manuscript is organized as follows. In section 2 we derive a field-
theoretical representation of the quantum quartic oscillator following earlier
applications of the stochastic approach. This is done in two steps. First, the
quartic part of the quantum quartic oscillator is decoupled by means of a
Hubbard-Stratonovich transformation [45, 46], which casts the problem into
that of the evolution of a quantum harmonic oscillator with time-dependent
frequency. Second, a Lie-algebraic transformation is used to express the time-
evolution operator in terms of ordinary exponentials featuring time-dependent
classical coefficients [47, 48, 19], which we refer to as stochastic
variables. The transformation yields a set of differential equations for the
stochastic variables, which effectively describe the exact quantum dynamics of
the quantum quartic oscillator. We also show that the formalism can be applied
to quartic Hamiltonians with time-dependent coefficients. In section 3 we
explain how expectation values of operators are calculated within this
formalism, providing a general recipe and exact equations for a range of
observables and initial states. The classical formulas we obtain are then
benchmarked by considering two exactly solvable limits: the quantum harmonic
oscillator, and the inherently classical limit in which all operators in the
Hamiltonian commute with each other, corresponding to the disappearance of the
kinetic energy in the Hamiltonian. In section 4, we elaborate on the
stochastic interpretation of the present approach and discuss the conditions
under which the field theory we obtain can be numerically simulated by means
of stochastic processes. Building on this discussion, we show that for the
quantum quartic oscillator with generic parameters (i.e., away from the two
exactly solvable limits) our method can be benchmarked by applying a numerical
stochastic scheme to evaluate time-evolved observables. In section 5 we adopt
a field-theoretical viewpoint and develop a functional expansion of the time-
evolution operator about the harmonic limit; we show that this is equivalent
to the standard perturbative Dyson series. In section 6, we discuss how the
semiclassical propagator and the partition function of the quantum quartic
oscillator can be analytically recovered from the proposed field-theoretical
picture. We present our conclusions in section 7, summarizing our results and
outlining directions for further research. Several appendices cover
mathematical details.

##  2 Stochastic Representation of the Quantum Quartic Oscillator

In this section, we introduce the stochastic approach briefly described in the
introduction alongside with its novel application to bosonic systems. In
particular, we focus on the quantum quartic oscillator, which is described by
the Hamiltonian

$\hat{H}\equiv\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}\hat{x}^{2}+\frac{\lambda}{4}\hat{x}^{4}\
.$ (1)

Here, the position operator $\hat{x}$ and the momentum operator $\hat{p}$
satisfy the canonical commutation relation $[\hat{x},\hat{p}]=i\hbar$. The
mass of the oscillator is denoted by $m$, the harmonic frequency by $\omega$
and the quartic potential is parametrized by the coupling constant
$\lambda\geq 0$. Since an exact solution is not available for $\lambda\neq 0$,
different approaches have been used in order to obtain insights on the quantum
quartic oscillator, e.g., perturbation theory [31, 30, 36] and semiclassical
approximations [42, 33, 32].

Here, we develop an alternative exact theoretical formulation of the problem,
following the disentanglement approach recently applied to an ensemble of
interacting quantum spin systems [15, 16, 17, 18]. In particular, we
investigate the unitary dynamics of the quantum quartic oscillator by exactly
mapping it to stochastically driven operators, which have spin-like
properties. This is done by employing a functional representation of the time-
evolution operator

$\displaystyle\hat{U}(t)\equiv\exp\left(-\frac{i}{\hbar}\hat{H}t\right).$ (2)

First, we decouple the quartic term in the exponent of Eq. (2) via a Hubbard-
Stratonovich transformation [46, 45] and trotterize [49, 50] the time-
evolution operator on the time interval $\tau_{n}\equiv t/n$, in the limit
$n\to\infty$, i.e.,

$\hat{U}(t)=\lim_{n\to\infty}\left(\exp\left[-\frac{i\tau_{n}}{\hbar}\frac{\hat{p}^{2}}{2m}\right]\exp\left[-\frac{i\tau_{n}}{\hbar}\left(\frac{m}{2}\omega^{2}\hat{x}^{2}+\frac{\lambda}{4}\hat{x}^{4}\right)\right]\right)^{n}.$
(3)

To each of the $n$ Suzuki-Trotter factors appearing in Eq.(3), we apply a
Hubbard-Stratonovich transformation which allows us to replace the quartic
interaction with a quadratic one by introducing a real-valued auxiliary field
$\phi$, i.e.,

$\exp\left(-i\frac{\tau_{n}}{\hbar}\frac{\lambda}{4}\hat{x}^{4}\right)=\sqrt{\frac{\tau_{n}}{i\lambda\hbar\pi}}\int_{-\infty}^{\infty}\mathrm{d}\phi\,\exp\left[\frac{i\tau_{n}}{\hbar}\left({\frac{\phi^{2}}{\lambda}-\hat{x}^{2}\,\phi}\right)\right],$
(4)

where we fix $\sqrt{i}=e^{i\pi/4}$ due to the positivity of $\lambda$, in
order to ensure convergence. Equation (4) is derived in Appendix A.1. We can
substitute the expression in Eq. (4) to each slice of the Suzuki-Trotter
decomposition in Eq. (3), labeling the corresponding integration variable
$\phi_{k}$ by $k\in\\{1,\dots,n\\}$. In the continuum limit
$n\rightarrow\infty$, we can express $\hat{U}(t)$ as a functional integral
with respect to the Hubbard-Stratonovich field $\phi(t)$. The corresponding
measure is given by
$\mathcal{D}\phi(t)\equiv\prod_{k}^{n}\sqrt{\tau_{n}/(i\lambda\hbar\pi)}\,\mathrm{d}\phi_{k}$
in the limit $n\rightarrow\infty$, and we find

$\hat{U}(t)=\int\mathcal{D}\phi\,e^{iS_{0}[\phi]}\,\mathbb{T}\exp\left[-\frac{i}{\hbar}\int_{0}^{t}\mathrm{d}\tau\left(\frac{\hat{p}^{2}}{2m}+\frac{m}{2}\Omega^{2}(\tau)\,\hat{x}^{2}\right)\right],$
(5)

where $\mathbb{T}\exp\left(\cdot\right)$ denotes the time-ordered exponential.
The quadratic coupling is absorbed into an effective time-dependent real-
valued frequency $\Omega^{2}(t)$, defined by

$\Omega^{2}(t)\equiv\omega^{2}+\frac{2\phi(t)}{m},$ (6)

and $S_{0}[\phi]$ denotes the Gaussian scalar action

$S_{0}[\phi]=\frac{1}{\hbar\lambda}\int_{0}^{t}\mathrm{d}\tau\,\phi^{2}(\tau).$
(7)

Equation (5) casts the time-evolution operator of the quartic oscillator as an
expectation value with respect to the Gaussian action $S_{0}[\phi]$ of the
propagator of a harmonic oscillator with time-dependent frequency $\Omega(t)$.
The associated effective Hamiltonian reads

$\displaystyle\hat{H}_{S}(t)\equiv\frac{\hat{p}^{2}}{2m}+\frac{m}{2}\Omega^{2}(t)\hat{x}^{2}.$
(8)

Despite the apparent simplification of dealing with the time-evolution
operator of a quadratic Hamiltonian, time-ordering prevents the direct
evaluation of the action of the operator in Eq. (5). This difficulty can be
circumvented by means of a Lie-algebraic disentanglement transformation [15,
16, 17]. We define a new set of operators, whose linear combination with
suitable coefficients reproduces the Hamiltonian in Eq. (1), i.e.,

$\displaystyle\hat{S}^{+}\equiv\frac{\hat{x}^{2}}{2\hbar},\qquad\hat{S}^{z}\equiv
i\frac{\\{\hat{x},\hat{p}\\}}{4\hbar},\quad\text{and}\quad\hat{S}^{-}\equiv\frac{\hat{p}^{2}}{2\hbar}.$
(9)

These operators satisfy the commutation relations of the SU(2) algebra, viz.,

$[\hat{S}^{+},\hat{S}^{-}]=2\hat{S}^{z}\;\;\;\text{and}\;\;\;[\hat{S}^{z},\hat{S}^{\pm}]=\pm\hat{S}^{\pm}.$
(10)

However, these operators differ from the conventional spin operators, since
$\hat{S}^{+}=(\hat{S}^{+})^{\dagger}$ and
$\hat{S}^{-}=(\hat{S}^{-})^{\dagger}$ are Hermitian. The operators defined in
Eq. (9) allow one to write the effective Hamiltonian as

$\hat{H}_{S}(t)=\hbar\left[\frac{\hat{S}^{-}}{m}+m\Omega^{2}(t)\,\hat{S}^{+}\right].$
(11)

Following Refs. [47, 48, 19], we can express the time-ordered exponential in
Eq. (5) as

$\hat{U}_{S}\equiv\mathbb{T}\exp\left[-\frac{i}{\hbar}\int_{0}^{t}\mathrm{d}\tau\left(\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\Omega^{2}(\tau)\,\hat{x}^{2}\right)\right]=e^{\xi^{+}(t)\,\hat{S}^{+}}e^{\xi^{z}(t)\,\hat{S}^{z}}e^{\xi^{-}(t)\,\hat{S}^{-}},$
(12)

with suitable “stochastic” variables $\xi^{+},\xi^{z}$ and $\xi^{-}$. In
section 4 we shall discuss the interpretation of these variables as stochastic
processes, which motivates their naming.

The time evolution of $\xi^{+},\xi^{z}$ and $\xi^{-}$ is obtained by imposing
that the factorized expression on the right-hand side of Eq. (12) satisfies
the same Heisenberg equation as $\hat{U}_{S}$ [19], leading to

$i\frac{\mathrm{d}\xi^{+}}{\mathrm{d}t}+\frac{1}{m}(\xi^{+})^{2}=m\Omega^{2}(t),\qquad
i\frac{\mathrm{d}\xi^{z}}{\mathrm{d}t}+\frac{2}{m}\xi^{+}=0,\qquad
i\frac{\mathrm{d}\xi^{-}}{\mathrm{d}t}-\frac{e^{\xi^{z}}}{m}=0.$ (13)

The initial condition $\hat{U}(0)=\mathds{1}$ imposes
$\xi^{+}(0)=\xi^{z}(0)=\xi^{-}(0)=0$. Note that, for a real-valued field
$\phi$, $\xi^{+}$ and $\xi^{-}$ are purely imaginary while $\xi^{z}$ is real,
implying that the exponential operators in Eq. (12) are unitary. We shall see
that the stochastic variables are generally complex and thus the product of
exponential operators in Eq. (12) is not unitary. The unitarity of the time-
evolution operator is however recovered upon averaging, as described in Eq.
(5). Substituting Eq. (12) into Eq. (5), allows one to express the time-
evolution operator as

$\hat{U}(t)=\left<e^{\xi^{+}(t)\hat{S}^{+}}e^{\xi^{z}(t)\hat{S}^{z}}e^{\xi^{-}(t)\hat{S}^{-}}\right>_{\phi},$
(14)

where $\left<\dots\right>_{\phi}$ denotes the average with respect to the
Gaussian field $\phi$, as defined in Eq. (5). The representation in Eq. (14)
is exact and allows us to map the quantum dynamics on the time evolution of
the stochastic parameters, see Eq. (13). Since this mapping is exact, the
ensemble of trajectories $\xi^{+}(t)$, $\xi^{z}(t)$, $\xi^{-}(t)$ determined
by the fields $\phi(t)$ encodes all the information about the underlying
quantum problem. Moreover, Eq. (14) suggests that the time evolution is given
by a weighted statistical average of successive actions of the exponential
operators.

It is worth noting that the operators in Eq. (12) are the matrix elements of
the covariance matrix

$\begin{pmatrix}\hat{x}^{2}&\frac{1}{2}\\{\hat{x},\hat{p}\\}\\\
\frac{1}{2}\\{\hat{x},\hat{p}\\}&\hat{p}^{2}\\\ \end{pmatrix},$ (15)

customarily used in the study of the dynamics of Gaussian wave packets under a
quantum oscillator Hamiltonian [51, 52]. Indeed, as we shall further discuss
below, the operators in Eq. (14) preserve the Gaussianity of a wave-packet and
justifies why it is convenient to work within this setting.

We can physically understand the action of the individual exponential
operators in Eq. (12), and hence of $\hat{U}_{S}$, by studying their effect on
a Gaussian wave packet $\ket{\psi}$, generally given by

$\ket{\psi}=\int_{-\infty}^{\infty}\frac{\mathrm{d}x}{(\pi\sigma^{2})^{1/4}}\exp\left(-\frac{(x-a)^{2}}{2\sigma^{2}}+i(x-a)k\right)\ket{x}.$
(16)

The wave packet in Eq. (16) is parametrized by its average position $\langle
x\rangle\equiv\bra{\psi}\hat{x}\ket{\psi}=a$, variance $\langle
x^{2}\rangle_{c}\equiv\bra{\psi}(\hat{x}-a)^{2}\ket{\psi}=\sigma^{2}/2$ (where
$\langle\dots\rangle_{c}$ denotes the connected component of the expectation
value) and average momentum $\langle
p\rangle\equiv\bra{\psi}\hat{p}\ket{\psi}=k$. The variance of the momentum
operator is given by $\langle
p^{2}\rangle_{c}\equiv\bra{\psi}(\hat{p}-k)^{2}\ket{\psi}=(2\sigma^{2})^{-1}$.
Here $\ket{x}$ denotes an eigenstate of the position operator and we set
$\hbar=1$. In order to visualize the action of the operators introduced in Eq.
(12), it is convenient to see how the average and variance of the momentum and
position operators transform. Namely, an operator of the type

  * •

$\exp(\xi^{-}\hat{S}^{-})$ acts on a Gaussian wave packet $\ket{\psi}$ by
leaving the momentum cumulants $\langle p\rangle$ and $\langle
p^{2}\rangle_{c}$ unaltered, while shifting, respectively, the average and the
variance of the position operator

$\langle x\rangle=a-\operatorname{Im}(\xi^{-})k,\qquad\langle
x^{2}\rangle_{c}=\sigma^{2}/2+(\operatorname{Im}(\xi^{-})/\sigma)^{2}/2,$

by terms that depend on the imaginary part of $\xi^{-}$;

  * •

$\exp(\xi^{z}\hat{S}^{z})$ rescales the $n-$th cumulant $\langle
x^{n}\rangle_{c}$ of the position by a homogeneous constant $e^{-n\xi^{z}/2}$,
while the $n-$th momentum cumulant $\langle\hat{p}^{n}\rangle_{c}$ by
$e^{n\xi^{z}/2}$. For the first cumulants we have $\langle
x\rangle=a\,e^{-\xi^{z}/2}$, $\langle
x^{2}\rangle_{c}=\sigma^{2}e^{-\xi^{z}}/2$, $\langle p\rangle=ke^{\xi^{z}/2}$
and $\langle p^{2}\rangle_{c}=e^{\xi^{z}}/(2\sigma^{2})$. Thanks to the
homogeneity of the transformation, the products $\langle x\rangle\langle
p\rangle$ and $\langle x^{2}\rangle_{c}\langle p^{2}\rangle_{c}$ are
preserved;

  * •

$\exp(\xi^{+}\hat{S}^{+})$ leaves unaltered the position cumulants $\langle
x\rangle$ and $\langle x^{2}\rangle_{c}$, while shifting, respectively the
average and the variance of the momentum operator

$\langle p\rangle=k+\operatorname{Im}(\xi^{+})a,\qquad\langle
p^{2}\rangle_{c}=(2\sigma^{2})^{-1}+(\operatorname{Im}(\xi^{+})\sigma)^{2}/2,$

by terms that depend on the imaginary part of $\xi^{+}$.

In order to further understand the physical significance of the action of the
exponential operators in Eq. (12) on a Gaussian wave packet $\ket{\psi}$, it
is useful to consider the exact phase space representation given by the Wigner
function [53, 54, 55], defined as

$W(x,p)\equiv\frac{1}{\pi}\int\mathrm{d}y\,e^{-i2yp}\langle
x+y|\psi\rangle\langle\psi|x-y\rangle.$ (17)

The knowledge of $W(x,p)$ allows one to compute the expectation value of any
operator $O(\hat{x},\hat{p})$, expressed as a function of $\hat{x}$ and
$\hat{p}$, as
$\bra{\psi}O(\hat{x},\hat{p})\ket{\psi}=\int\mathrm{d}x\int\mathrm{d}p\,W(x,p)O(x,p)$.
This generalizes the results obtained above for the first cumulants of the
wave packet, which are retrieved by identifying the operator
$O(\hat{x},\hat{p})$ with $\hat{x}$, $\hat{p}$, $\hat{x}^{2}$, $\hat{p}^{2}$,
and their connected expectation value. Figure 1 shows how the Wigner function

$W(x,p)=\frac{1}{\pi}\exp\left[-\frac{(x-a)^{2}}{\sigma^{2}}-\sigma^{2}(p-k)^{2}\right],$
(18)

of a Gaussian wave packet is transformed upon the action of the three
exponential operators in Eq. (12). In general, the exponential operators in
Eq. (14) transform the parameters of the Gaussian wave packet by altering its
original Heisenberg uncertainty relation
$\langle\hat{x}^{2}\rangle_{c}\langle\hat{p}^{2}\rangle_{c}=1/2$ [56], while
they preserve the Gaussian structure of the relative Wigner function. For
further details, see Appendix B.

Figure 1:  Contour plots of different Wigner functions. Panel (a) shows the
Wigner function $W(x,p)$ for a Gaussian wave packet, see Eq. (18). In the $x$
and $p$ direction, $W(x,p)$ is a Gaussian centered around the phase space
point $(a,k)$ with variance $\sigma^{2}/2$ and $(2\sigma^{2})^{-1}$,
respectively. We chose $a=0,\,k=1$, and $\sigma=1$. Panel (b) shows the Wigner
function $W_{+}(x,p)$ for the wave function
$\ket{\psi_{+}}\equiv\exp(\xi^{+}\hat{S}^{+})\ket{\psi}$, see Eq. (82). This
transformation shifts the $p$ variable by an $x$-dependent linear term, i.e.,
$p\rightarrow p-x\,\operatorname{Im}(\xi^{+})$; here we choose $\xi^{+}=1.5i$.
Panel (c) shows the Wigner function $W_{z}(x,p)$ for the wave function
$\ket{\psi}\equiv\exp(\xi^{z}\hat{S}^{z})\ket{\psi}$, see Eq. (77). This
transformation uniformly rescales the variables $(x,p)$ to
$(x\,e^{\xi^{z}/2},p\,e^{-\xi^{z}/2})$; here we chose $\xi^{z}=1.5$. Panel (d)
shows the Wigner function $W_{-}(x,p)$ for the wave function
$\ket{\psi}\equiv\exp(\xi^{-}\hat{S}^{-})\ket{\psi}$, see Eq. (71). This
transformation shifts the $x$ variable by a $p$-dependent linear term
according to $x\rightarrow x+p\,\operatorname{Im}(\xi^{-})$; we chose
$\xi^{-}=1.5i$.

##  3 Physical Observables

In this section, we illustrate how expectation values of observables can be
expressed in the stochastic formalism by considering the position and momentum
operator. In general, for a system prepared in a state $\ket{\psi}$, the
expectation value of an observable $\hat{O}$ is given by
$\bra{\psi}\hat{O}\ket{\psi}$. In particular, we refer to
$\bra{\psi}\hat{O}^{n}\ket{\psi}$ as the $n$-th moment of the operator
$\hat{O}$ with respect to the state $\ket{\psi}$. We can then express quantum
expectation values as functional averages by replacing each time-evolution
operator by its exact representation given by Eq. (14). This requires
introducing independent Hubbard-Stratonovich fields $\phi$ and $\bar{\phi}$
for the two time-evolution operators. We refer to these fields as forward
($\phi$) and backward ($\bar{\phi}$) fields in analogy to the nomenclature of
the Schwinger-Keldysh formalism [57, 58, 59].

###  3.1 Dynamics of a Gaussian Wave Packet

The dynamics of a particle in the presence of a quartic potential can be
studied via the time-dependent moments of its position and momentum. We model
the particle as a Gaussian wave packet, see Eq. (16), whose time evolution is
governed by the time-evolution operator in Eq. (14). The evolved state is thus
simply obtained by the subsequent action of the exponential operators in Eq.
(14) and we find

$\ket{\psi(t)}=\left<\int\mathrm{d}y\frac{\sqrt{\beta(t)\sigma}}{\pi^{1/4}}\exp\left(-\frac{\sigma^{2}k^{2}}{2}-\frac{y^{2}}{2}\gamma(t)+i\sigma^{2}\beta(t)\mu\,y+\frac{\sigma^{4}\mu^{2}}{2}\alpha(t)\right)\ket{y}\right>_{\phi}.$
(19)

Here, we have introduced the generalized initial momentum $\mu$ of the wave
packet as $\mu\equiv k-ia\sigma^{-2}$ and the variables $\alpha(t),\beta(t)$
and $\gamma(t)$, which depend on the stochastic variables as

$\begin{split}\alpha(t)\equiv\frac{1}{\sigma^{2}-\xi_{\vphantom{0}}^{-}(t)},\qquad\beta(t)\equiv\frac{e^{\xi^{z}(t)/2}}{\sigma^{2}-\xi_{\vphantom{0}}^{-}(t)},\qquad\gamma(t)\equiv\frac{e^{\xi^{z}(t)}}{\sigma^{2}-\xi_{\vphantom{0}}^{-}(t)}-\xi^{+}(t).\end{split}$
(20)

The time evolution of $\alpha,\beta$, and $\gamma$ is readily obtained from
Eq. (13), i.e.,

$\begin{split}i\frac{{\rm d}\gamma}{{\rm
d}t}=\frac{\gamma^{2}}{m}-m\Omega^{2},\qquad i\frac{{\rm d}\beta}{{\rm
d}t}=\frac{\beta\gamma}{m},\qquad i\frac{{\rm d}\alpha}{{\rm
d}t}=\frac{\beta^{2}}{m},\end{split}$ (21)

with initial conditions $\alpha(0)=\beta(0)=\gamma(0)=\sigma^{-2}$. These
parameters are generally complex. Notably, the presence of $\sigma^{2}>0$ in
Eq. (20) prevents divergences otherwise occurring for $\xi^{-}(t)=0$, e.g.,
for $t=0$.

The moments of the position operator on the Gaussian wave packet in Eq. (16)
can be explicitly obtained by computing the expectation values
$\bra{\psi(t)}\hat{x}^{n}\ket{\psi(t)}$. To this end, we replace the time-
evolved state $\ket{\psi(t)}$ by its representation involving the functional
average $\langle\dots\rangle_{\phi}$, derived in Eq. (14). A similar procedure
is applied to $\hat{U}^{\dagger}(t)$, associated with the field $\bar{\phi}$
with action $-S_{0}[\bar{\phi}]$. We denote by $\alpha$, $\beta$ and $\gamma$
the solutions of Eqs. (21) that depend on the forward field $\phi$, and
similarly we write $\bar{\alpha}$, $\bar{\beta}$, and $\bar{\gamma}$ for the
solutions of the complex conjugates of Eqs. (21) associated with the backward
field $\bar{\phi}$. Note that the fields $\phi$ and $\bar{\phi}$ are
independent. The $n$-th moment of the position operator is finally found as

$\bra{\psi(t)}\hat{x}^{n}\ket{\psi(t)}=\Biggl{\langle}\frac{2\sigma
i^{n}\sqrt{\beta\bar{\beta}}}{(2\Gamma)^{(n+1)/2}}\,\mathrm{H}_{n}\left(\frac{\sigma^{2}\Delta}{\sqrt{2\Gamma}}\right)\exp\left[-\sigma^{2}k^{2}+\frac{\sigma^{4}}{2}\left(A-\frac{\Delta^{2}}{\Gamma}\right)\right]\Biggl{\rangle}_{\phi,\,\bar{\phi}},$
(22)

where we defined the auxiliary variables $\Gamma\equiv\gamma+\bar{\gamma}$,
$\Delta\equiv\mu\beta-\mu^{*}\bar{\beta}$,
$A\equiv\mu^{2}\alpha+(\mu^{*})^{2}\bar{\alpha}$, and $\mathrm{H}_{n}$ denotes
the $n$-th degree Hermite polynomial and $\mu^{*}$ the complex conjugate of
$\mu$. Similarly, the expression of the moments of the momentum operator are
found to be

$\displaystyle\bra{\psi(t)}\hat{p}^{n}\ket{\psi(t)}$
$\displaystyle=\Biggl{\langle}2\sigma\left[2\left(\frac{1}{\gamma}+\frac{1}{\bar{\gamma}}\right)\right]^{-(n+1)/2}i^{n}\sqrt{\frac{\beta\bar{\beta}}{\gamma\bar{\gamma}}}\,\mathrm{H}_{n}\left(-i\frac{\sigma^{2}\left(\mu\frac{\beta}{\gamma}+\mu^{*}\frac{\bar{\beta}}{\bar{\gamma}}\right)}{\sqrt{2\left(\frac{1}{\gamma}+\frac{1}{\bar{\gamma}}\right)}}\right)$
(23)
$\displaystyle\times\exp\left\\{-\sigma^{2}k^{2}+\frac{\sigma^{4}}{2}\left[A-\mu^{2}\frac{\beta^{2}}{\gamma}-(\mu^{*})^{2}\frac{\bar{\beta}^{2}}{\bar{\gamma}}+\frac{\left(\mu\frac{\beta}{\gamma}+\mu^{*}\frac{\bar{\beta}}{\bar{\gamma}}\right)^{2}}{\left(\frac{1}{\gamma}+\frac{1}{\bar{\gamma}}\right)}\right]\right\\}\Biggl{\rangle}_{\phi,\,\bar{\phi}}.$

Note that these moments are formally retrieved from Eq. (22) by substituting

$\displaystyle\gamma\rightarrow\gamma^{-1},\quad\beta\rightarrow-i\beta/\gamma,\quad\text{and}\quad\alpha\rightarrow\alpha-\beta^{2}/\gamma.$

As in the case of the time-evolution operator and the wave packet evolution,
it is possible to express the dynamics of an observable as the expectation
value of functions of the auxiliary parameters $\alpha$, $\beta$ and $\gamma$,
which depend on the field $\phi$. Below, we will see how we may use these
expressions for numerical calculations. As a final remark, we note that the
convergence of Eq. (19) requires that $\operatorname{Re}(\gamma)>0$, which
follows naturally from the unitarity of the exponential operators in Eq. (12),
see Appendix C for further details.

###  3.2 Exactly Solvable Cases

To the best of our knowledge, the non-linear evolution of the system of Eqs.
(13) cannot in general be solved exactly. However, exact solutions can be
found in two cases: the harmonic limit, in which Eqs. (13) become purely
deterministic, and the commuting limit, in which the solutions of Eqs. (13)
can be expressed in terms of the time integral of the field $\phi$. We use
these exactly solvable instances as benchmarks for the stochastic formalism as
well as for developing a physical intuition of its significance.

####  3.2.1 Harmonic case

In the case $\lambda=0$ of the harmonic oscillator, the absence of the quartic
term implies that Eqs. (13) reduce to the system of ordinary differential
equations

$\displaystyle
i\frac{\mathrm{d}\xi^{+}}{\mathrm{d}t}+\frac{1}{m}(\xi^{+})^{2}=m\omega^{2},\qquad
i\frac{\mathrm{d}\xi^{z}}{\mathrm{d}t}+\frac{2}{m}\xi^{+}=0,\qquad
i\frac{\mathrm{d}\xi^{-}}{\mathrm{d}t}-\frac{e^{\xi^{z}}}{m}=0,$ (24)

which can be solved explicitly, i.e.,

$\xi^{+}=-im\omega\tan(\omega t),\qquad\xi^{-}=-\frac{i}{m\omega}\tan(\omega
t),\qquad\xi^{z}=-\log\cos^{2}(\omega t).$ (25)

This amounts to a known, exact parameterization of the quantum harmonic
oscillator in terms of classical variables [19]. The time evolution of the
stochastic variables $\xi^{+,-,z}$ for the harmonic oscillator shows periodic
divergences at $t=t_{n}=\pi(1+2n)/(2\omega)$ with integer $n$, which, however,
cancel out in the analytic computations of observables. This issue can be
circumvented by equivalently considering the time evolution of the variables
$\alpha$, $\beta$, $\gamma$, introduced in section 3.1. The solutions (25) can
be plugged in the expressions of the observables, obtained in section 3, in
order to compute exactly the corresponding dynamics. For instance, by
inserting Eqs. (25) into Eq. (22), we retrieve the expressions of the moments
of the position and momentum operator for the Gaussian wave packet:

$\displaystyle\left<\hat{x}^{n}\right>$
$\displaystyle=\left(\frac{i}{2\sigma}\right)^{n}\left[\sigma^{4}\cos^{2}(\omega
t)+x_{0}^{4}\sin^{2}(\omega
t)\right]^{\frac{n}{2}}\,\mathrm{H}_{n}\left(-\frac{i\sigma[a\cos(\omega
t)+kx_{0}^{2}\sin(\omega t)]}{\sqrt{\sigma^{4}\cos^{2}(\omega
t)+x_{0}^{4}\sin^{2}(\omega t)}}\right),$ (26)
$\displaystyle\left<\hat{p}^{n}\right>$
$\displaystyle=\left(\frac{i}{2x_{0}^{2}\sigma}\right)^{n}\left[x_{0}^{4}\cos^{2}(\omega
t)+\sigma^{4}\sin^{2}(\omega
t)\right]^{\frac{n}{2}}\,\mathrm{H}_{n}\left(-\frac{i\sigma[kx_{0}^{2}\cos(\omega
t)-a\sin(\omega t)]}{\sqrt{x_{0}^{4}\cos^{2}(\omega
t)+\sigma^{4}\sin^{2}(\omega t)}}\right),$

where we introduced the typical harmonic oscillator length
$x_{0}\equiv(m\omega)^{-1/2}$.

####  3.2.2 Commuting limit

We study a particular case of the Hamiltonian in Eq. (1) in which the time
evolution of the stochastic variables in Eq. (13) can be exactly solved. We
consider the limit $m\rightarrow\infty$ with $m\omega^{2}$ held constant in
order to ensure a constant finite energy. In this limit, the kinetic part of
the Hamiltonian is suppressed and it coincides with the quartic potential.
Accordingly, the Hamiltonian is only a function of the position operator and
does not contain the momentum operator. In this sense we refer to this
scenario as the commuting limit. Correspondingly, the non-linear terms in the
differential equations (13) vanish, allowing us to express the explicit
solution as

$\xi^{+}=-i\omega^{2}mt-2i\int_{0}^{t}\mathrm{d}\tau\,\phi(\tau),\qquad\xi^{-}=\xi^{z}=0,$
(27)

where $\xi^{+}$ depends on the variable
$\mathcal{W}\equiv\int_{0}^{t}\mathrm{d}s\,\phi(s)$, and allows us to express
the average $\langle\cdots\rangle_{\phi}$ in Eq. (14) as an average with
respect to the Gaussian weight of $\mathcal{W}$, i.e.,

$\hat{U}(t)=\int_{-\infty}^{+\infty}\mathrm{d\mathcal{W}}\,\frac{e^{i\mathcal{W}^{2}/(\lambda
t)}}{\sqrt{i\pi\lambda
t}}\exp\left[-i\left(2\mathcal{W}+m\omega^{2}t\right)\hat{S}^{+}\right].$ (28)

From Eq. (28), it is apparent that, as expected in this case, the time
evolution operator commutes with operators that depend only on the position
operator $\hat{x}$. This implies the absence of dynamics for the particle
position, compatibly with the vanishing kinetic energy. On the other hand, the
moments of the momentum $\hat{p}$ grow in time, as a consequence of the
Heisenberg uncertainty principle. We compute these moments by setting the
auxiliary variables in Eq. (20) to $\alpha=\beta=\sigma^{-2}$ and
${\gamma=\sigma^{-2}-i(2\mathcal{W}+\omega^{2}mt)}$ in Eq. (23) and by
substituting the average $\langle\cdots\rangle$ with the integral
$\int_{-\infty}^{+\infty}\mathrm{d\mathcal{W}}\,e^{i\mathcal{W}^{2}/(\lambda
t)}/\sqrt{i\pi\lambda t}$ and similarly for the dual variables $\bar{\alpha}$,
$\bar{\beta}$ and $\bar{\gamma}$. These expectation values can then be
evaluated numerically, as we discuss in the next section.

Figure 2: Panels (a) and (b) display a numerical estimate of the average value
$\langle t_{\gamma}\rangle$ and the standard deviation
$\sigma_{\gamma}\equiv\sqrt{\langle t_{\gamma}^{2}\rangle-\langle
t_{\gamma}\rangle^{2}}$ of the random variable $t_{\gamma}$, respectively, for
the quartic oscillator, where $t_{\gamma}$ is the time at which
$\operatorname{Re}(\gamma(t_{\gamma}))=0$ for the first time. We have used the
Euler scheme with time step $\Delta t=2/\pi\times 10^{-3}$ to simulate the
time evolution of $\gamma$ with physical parameters are $\omega=1$, $m=10$,
and extracted the value of $t_{\gamma}$ for $N=10^{4}$ trajectories. Panels
(c) and (d) show the same quantities and parameters except for the presence of
a time-dependent potential $\omega^{2}(t)=\sin^{2}(t)$ and
$\lambda(t)=\lambda_{0}\sin^{2}(t)$. Note that, in general, $\langle
t_{\gamma}\rangle\sim\sigma_{\gamma}$ at large values of $\lambda$: these
strong fluctuations increase the probability of
$\operatorname{Re}(\gamma(t_{\gamma}))=0$ at smaller values of $t_{\gamma}$ as
the coupling $\lambda$ increases.

##  4 Stochastic interpretation and numerical benchmark

In this section we show how the formalism presented above can be interpreted
in terms of stochastic processes, which also allows us to benchmark our
approach in cases where the model is not exactly solvable. In particular, by
rotating the integration contour of the variable $\phi$ in the Hubbard-
Stratonovich transformation in Eq. (4), after a generalization to the case of
time-dependent couplings as in Eq. (68) of, cf., Appendix A.2, it is possible
to show that

$\exp\left(-i\frac{\tau_{n}}{\hbar}\frac{\lambda_{n}}{4}\hat{x}^{4}\right)=\sqrt{\frac{\tau_{n}}{\hbar\pi}}\int_{-\infty}^{\infty}\mathrm{d}\phi\,\exp\left[-\frac{\tau_{n}}{\hbar}\left(\phi^{2}+i\sqrt{i\lambda_{n}}\,\hat{x}^{2}\,\phi\right)\right],$
(29)

which follows from the change of variable $\phi/\sqrt{i}\rightarrow\phi$ in
Eq. (68). Equation (29) allows one to represent $\hat{U}$ as

$\hat{U}(t)=\int\mathcal{D}\phi\,e^{-S_{0}[\phi]}\,\hat{U}_{S}[\phi],$ (30)

where $S_{0}$ is given by Eq. (7) and the time-evolution operator
$\hat{U}_{S}$ in Eq. (12) displays the corresponding effective frequency
$\Omega^{2}(t)\equiv\omega^{2}(t)+2\sqrt{i\lambda(t)}\phi(t)/m$. The
exponential of the action $S_{0}$ in Eq. (7) can be identified as a Gaussian
probability measure for the field $\phi(t)$, whose time integral can be
interpreted as a Wiener process [20]. Accordingly, Eqs. (13) can be understood
as complex stochastic differential equations with a Gaussian white noise
$\phi$. In this reformulation, however, the stochastic variables are generally
complex, and do not preserve the unitarity of $\hat{U}_{S}$, which is only
retrieved upon averaging. As a consequence, the expression inside the average
$\langle\cdots\rangle_{\phi}$ of the observable expressions in Eq. (19) does
not have to be convergent at all times and for all values of the quartic
coupling. More precisely, since the effective frequency $\Omega^{2}(t)$ is
complex, $\operatorname{Re}(\gamma)$ may attain a negative value after a
certain time $t_{\gamma}>0$, even though the initial value
$\gamma(0)=\sigma^{-2}>0$ is positive and real. One can numerically check that
the average value $\langle t_{\gamma}\rangle$ of $t_{\gamma}$, interpreted as
a random variable, decreases upon increasing the strength of the quartic
coupling $\lambda$, see Fig. 2. Indeed, $t_{\gamma}$ is the first-passage time
to the origin for the random variable $\operatorname{Re}(\gamma)$. This
constitutes a limitation to the numerical application of the stochastic
approach in the large-$\lambda$ regime, where $\langle t_{\gamma}\rangle$ is
comparable with the standard deviation $\sigma_{\gamma}\equiv\sqrt{\langle
t_{\gamma}^{2}\rangle-\langle t_{\gamma}\rangle^{2}}$. In this stochastic
description, as reported in Appendix C, the divergences originate from the
non-commutativity of the average $\langle\cdots\rangle_{\phi}$ over
trajectories and the action of the operator $\hat{U}_{S}$ on a prescribed
initial state $\ket{\psi}$ since, due to the non-unitarity of $\hat{U}_{S}$,
the relation
$\langle\hat{U}_{S}\rangle_{\phi}\ket{\psi}=\langle\hat{U}_{S}\ket{\psi}\rangle_{\phi}$
may not be satisfied. In the commuting limit, the above considerations
translate in the simple change of variable
$\mathcal{W}/\sqrt{i}\equiv\mathcal{W}^{\prime}$ in Eq. (28), where this new
$\mathcal{W}^{\prime}$ can be interpreted a Gaussian random number with
variance $\lambda t/2$.

Figure 3: Time evolution of the first and second moment of the
$\hat{p}$-operator for a Gaussian wavepacket with $\sigma=0.5,\,k=1.0,\,a=0.5$
and evolving according to Eq. (1) with $\lambda=0.2$ in the commuting limit.
The dashed line is computed by integrating numerically the Schrödinger
equation with the Crank-Nicholson method [60] with lattice spacing $\Delta
x=6\times 10^{-4}$ and time step $\Delta t=\pi\times 10^{-4}$. Blue dots are
computed with the stochastic method with $\Delta t=10^{-4}$, and by sampling
$N=1.2\times 10^{6}$ Gaussian random numbers. Error bars, corresponding to the
statistical standard deviation over the sampling average, are not visible on
the scale of the plot.

In spite of this limitation, the possibility to evaluate observables
numerically by simulating classical stochastic dynamics allows us to further
benchmark our approach. As a first check, we determine numerically the
dynamics of the average momentum $\langle p\rangle$ for a Gaussian wave packet
in the commuting limit. As we have shown, the stochastic differential
equations are exactly solvable in this limit. It is thus possible to obtain
directly the expressions for observables at a given time $t$ without having to
integrate the time evolution numerically. These expressions are known
functions of the time integral
$\xi^{+}=-2\sqrt{i}\mathcal{W}^{\prime}-i\omega^{2}mt$, which can be
numerically simulated by drawing Gaussian random numbers with zero mean and
variance given by $\lambda t/2$. Figure 3 shows the comparison of the
numerical prediction of the dynamics of the first two moments of the position
of a wave packet between the Crank-Nicholson method (dashed line) [60] used to
integrate the Schrödinger equation numerically, and the prediction based on
the stochastic interpretation of Eq. (30) (dots), finding good agreement. As a
further validation of the presented stochastic description, we evaluate the
dynamics of the expectation values in Eq. (22) within a range of parameters
where no exact solutions are available. We determine our numerical results up
to a time $t<\langle t_{\gamma}\rangle$, where no divergences are actually
detected. For this purpose, we use an Euler discretization scheme [20] with
time step $\Delta\tau=10^{-5}$ to solve the complex-valued stochastic
differential equations (13) for a given realization of the Wiener process
$\phi(t)$. Once a sufficiently large number of realizations for the stochastic
variables $\xi^{+},\,\xi^{z}$, and $\xi^{-}$ or $\alpha,\,\beta$, and $\gamma$
are known, by averaging with respect to them, it is possible to compute the
expectation value of a desired observable, see e.g., Eq. (22) or (23). Figure
4 shows the time evolution of the first moments of the position operator for a
Gaussian wave packet for various choices of the parameters. In particular, we
compare the numerical prediction of the stochastic method (dots) with standard
integration of the Schrödinger equation with the Cranck-Nicholson method
(dashed line). Numerically, the proposed stochastic method has the advantage
that the time-evolution of the many trajectories of the stochastic parameters,
e.g., $\xi^{+,-,z}$, can be straightforwardly parallelized. On the other hand,
an increasing large number of realizations is needed in order to have accurate
predictions for observables at longer times or larger quartic coupling
strength $\lambda$, since fluctuations due to the noise grow correspondingly.
This is similar to the behavior found for quantum spins systems [22].

Figure 4: Time evolution of the expected (a) position $\langle\,x\,\rangle$
and (b) squared position $\langle\,x^{2}\,\rangle$ of a Gaussian wave packet
with $\sigma=0.5,\,k=1.0,$, and $a=0.5$, evolving according to Eq. (1) with
$m=10,\,\omega=1,\,\lambda=0.2$. The dashed line corresponds to the numerics
performed with the Crank-Nicholson method with space bin $\Delta x=10^{-5}$
and time step $\Delta t=10^{-3}$, while the blue dots are computed with the
stochastic method with $\Delta t=10^{-4}$, and sampling $N=1.2\times 10^{5}$
trajectories. The time $t$ has been chosen such that $t<\langle
t_{\gamma}\rangle$ with $\langle t_{\gamma}\rangle\simeq 27$, and
$\sigma_{\gamma}\simeq 19$, see Fig. 2a-b . Similarly, panels (c) and (d)
display the evolution of $\langle x\rangle$ and $\langle x^{2}\rangle$,
respectively, in the case of $\omega^{2}(t)=\omega_{0}^{2}\sin^{2}(t)$ and
$\lambda(t)=\lambda_{0}\sin^{2}(t)$ with $\omega_{0}=1$, $\lambda_{0}=0.2$,
and $\langle t_{\gamma}\rangle\simeq 31$, with $\sigma_{\gamma}\simeq 18$, see
Fig. 2c-d. As above, the dashed line corresponds to the numerics performed
with the Crank-Nicholson method with space bin $\Delta x=10^{-6}$ and time
step $\Delta t=10^{-3}$, while blue dots are computed with the stochastic
method with $\Delta t=10^{-4}$, and sampling $N=1.2\times 10^{5}$
trajectories. Error bars are given by statistical standard deviation over the
trajectories average and are only visible in panel (d).

In summary, we have demonstrated that our formalism allows one to compute
quantum observables for an interacting bosonic system from averaging classical
stochastic processes, but we found that there are limitations to the
applicability of this numerical technique. The mapping to stochastic processes
also made it possible to further benchmark our approach in non-solvable cases.
Since the stochastic description discussed in this section is formally
identical to the field theory introduced in section 2, we refer to the present
method as the “stochastic approach", although only its numerical application
discussed here makes explicit use of stochastic processes.

##  5 Perturbative Expansion

In the stochastic approach, the time-evolution operator $\hat{U}$ in Eq. (14)
is represented as a functional average over classical fields $\phi$. In this
section, we show how this can be used to derive a perturbative expansion of
$\hat{U}$ for the quartic oscillator in terms of the time-evolution operator
$\hat{U}_{0}(t)$ of the harmonic case obtained for $\lambda=0.$ We start by
rescaling the Hubbard-Stratonovich field as $\phi=\sqrt{\lambda}\,\varphi$ in
the functional integral representation of $\hat{U}(t)$ in Eq. (5), yielding

$\hat{U}(t)=\int\mathcal{D}\varphi\,e^{i\int_{0}^{t}\mathrm{d}s\,\varphi^{2}(s)}\,\hat{U}_{S}\left[\sqrt{\lambda}\varphi\right].$
(31)

Here, the functional $\hat{U}_{S}[\phi]$ (which is also a function of time) is
identified with the time-evolution operator of a harmonic oscillator with
time-dependent frequency, given by Eq. (12). By Taylor-expanding the
functional $\hat{U}_{S}$ around $\phi=0$, corresponding to the harmonic time-
evolution operator $\hat{U}_{0}$, and calculating the resulting Gaussian
integrals, we get an asymptotic series for the propagator

$\displaystyle\hat{U}(t)=\sum_{n=0}^{\infty}\left(i\frac{\lambda}{4}\right)^{n}\left(\prod_{m=1}^{n}\int_{0}^{t}\mathrm{d}t_{m}\right)\frac{\delta^{2n}\hat{U}_{S}[\phi]}{\delta\phi(s_{1})\cdots\delta\phi(s_{2n})}\Bigg{|}_{\begin{subarray}{c}\phi=0\\\
s_{2n-1}=s_{2n}=t_{n}\end{subarray}},$ (32)

where, on the right-hand side, only even orders of functional derivatives
appear as a consequence of Wick’s theorem, leaving the functional derivative
of $\hat{U}_{S}$ evaluated at $\phi=0$ as the only unknown. The series in Eq.
(32) can be shown to be equivalent term by term to the Dyson series, see
Appendix D. The equivalence of the functional expansion about $\phi=0$ with
the Dyson series allows us to use this functional formulation to calculate
perturbative approximations of observables by field-theoretical means: we
express the time evolution operators in the stochastic formalism, such that
all operators are replaced by classical functionals, and then functionally
expand about the non-interacting case. As it is usually the case in
perturbative calculations, the asymptotic series in Eq. (32) is expected to
fail whenever we consider states for which the quartic term
$\lambda\,\hat{x}^{4}/4$ is not negligible relative to the harmonic
Hamiltonian. Indeed, it is a well-known fact that the Dyson series of the
quartic oscillator has a vanishing radius of convergence [30, 31].

##  6 Semiclassical Approximation

In this section, we show how the semiclassical approximation for the
propagator associated with $\hat{U}(t)$ in the representation of Eq. (14) and
the partition function of the system in Eq. (1) can be expressed within the
present formalism. Other than giving us an additional benchmark for the
theory, this shows that it is possible to find an alternative description of
the stationary trajectories contributing to the semiclassical approximation
for the quartic oscillator.

###  6.1 Propagator

The propagator $G(x_{f},t|x_{i},0)$, which gives the probability amplitude for
a particle located at $x_{i}$ at the initial time $t_{i}=0$ to reach the
position $x_{f}$ at time $t_{f}=t$, is defined by

$G(x_{f},t|x_{i},0)\equiv\bra{x_{f}}\hat{U}(t)\ket{x_{i}}.$ (33)

The semiclassical approximation of the propagator for the quantum quartic
oscillator has been extensively studied in the literature, see, e.g., Ref.
[32] for an overview. Here, we show how to express $G(x_{f},t|x_{i},0)$ in
terms of the stochastic variables. This expression can be derived by inserting
the representation of $\hat{U}$ in Eq. (14), by acting on an eigenstate of the
position $\ket{x_{i}}$ according to Eq. (86) in Appendix C, and finally by
projecting on $\bra{x_{f}}$. This leads to

$G(x_{f},t|x_{i},0)=\Biggl{\langle}\frac{1}{\sqrt{-2\pi\hbar\xi^{-}}}\exp\left[\frac{\xi^{z}}{4}+\frac{\xi^{+}x_{f}^{2}}{2\hbar}+\frac{(x_{f}e^{\xi^{z}/2}-x_{i})^{2}}{2\hbar\xi^{-}}\right]\Biggl{\rangle}_{\phi},$
(34)

where $\xi^{+,-,z}$ are evaluated at time $t.$ Note that the stochastic
variables in the above expression are functions of $\phi$ evaluated at the
final time $t=t_{f}-t_{i}$, while the initial and final position $x_{i}$ and
$x_{f}$ are fixed parameters.

In the harmonic case $\lambda=0$, by explicit substitution of Eq. (25) into
Eq. (34), we retrieve the known expression of the harmonic propagator $G_{\rm
HO}$, i.e.,

$G_{\rm HO}(x_{f},t|x_{i},0)=\sqrt{\frac{m\omega}{2i\pi\hbar\sin(\omega
t)}}\exp\left[\frac{im\omega}{2\hbar\sin(\omega
t)}\left((x_{f}^{2}+x_{i}^{2})\cos(\omega t)-2x_{i}\,x_{f}\right)\right].$
(35)

In the quartic case, inside the average in Eq. (34) we recognize a different
way of representing the propagator of an harmonic oscillator with time-
dependent frequency $\Omega^{2}(t)$. It is well-known that in the case of
quadratic interactions, even in the time-dependent case, the propagator can be
expressed in a closed form through the contributions arising from classical
paths, see, e.g., Ref. [32]. Accordingly, the propagator of a harmonic
oscillator with generic time-dependent frequency can be reformulated as

$G_{\rm HO}(x_{f},t|x_{i},0)=G_{\rm
HO}(0,t|0,0)\exp\left[\frac{i}{\hbar}S_{\rm HO}(x_{f},t|x_{i},0)\right].$ (36)

Here, the classical action $S_{\rm HO}$ of the harmonic oscillator is given by

$S_{\rm
HO}(x_{f},t|x_{i},0)\equiv\frac{m}{2}\int_{0}^{t}\mathrm{d}\tau\left[\dot{x}^{2}(\tau)-\Omega^{2}(\tau)x^{2}(\tau)\right],$
(37)

and is computed along the classical path $x(\tau)$ which satisfies the Euler-
Lagrange equation

$\ddot{x}(\tau)+\Omega^{2}(\tau)\,x(\tau)=0,$ (38)

with boundary conditions $x(0)=x_{i}$ and $x(t)=x_{f}$. The prefactor $G_{\rm
HO}(0,t|0,0)$ is given by

$G_{\rm HO}(0,t|0,0)=\sqrt{\frac{m}{2\pi i\hbar f(t)}},$ (39)

where the density of paths $f(t)$ is obtained, according to Gelfand-Yaglom
formula [61], as a solution of the differential equation in Eq. (38) with
$x\rightarrow f$, with initial conditions $f(0)=0$ and $\dot{f}(0)=1$. Note
that the propagator $G_{\rm HO}(0,t|0,0)$ can be represented via the Feynman
path integral associated to the quadratic action in Eq. (37) with boundary
conditions $x_{f}=x_{i}=0$, which entails that the function $f(t)$ is
proportional to the determinant of the linear operator
$d_{t}^{2}+\Omega^{2}(t)$ expressed through Eq. (39) [32].

The solutions of Eq. (38) depend on the realization of the field $\phi$ which
enters $\Omega$ according to Eq. (6). Alternatively, they can be expressed in
terms of the stochastic variables according to

$\displaystyle x(\tau|t)$
$\displaystyle=\frac{e^{-\xi^{z}(\tau)/2}}{\xi^{-}(t)}\left\\{x_{f}\xi^{-}(\tau)e^{\xi^{z}(t)/2}+x_{i}\left[\xi^{-}(t)-\xi^{-}(\tau)\right]\right\\},$
(40) $\displaystyle f(\tau)$
$\displaystyle=im\,\xi^{-}(\tau)\,e^{-\xi^{z}(\tau)/2},$

where we emphasize that $t,\,\,x_{i}$, and $x_{f}$ are fixed parameters, and
$\tau\in[0,t]$ is a variable. In turn, $\xi^{+,-,z}$ depend on $\phi$ via Eq.
(13). In order to simplify the notation, the dependence of $x$ and $f$ on
$x_{i}$ and $x_{f}$ and the functional dependence on $\phi$ are omitted.

Finally, by inserting Eq. (36) into (34), the propagator $G(x_{f},t|x_{i},0)$
for the quartic oscillator reads

$\displaystyle G(x_{f},t|x_{i},0)$
$\displaystyle=\int\mathcal{D}\phi\,\,\sqrt{\frac{m}{2\pi i\hbar
f(t)}}\exp\left[\frac{i}{\hbar\lambda}\int_{0}^{t}\phi^{2}(\tau)d\tau+\frac{i}{\hbar}S_{\rm
HO}(x_{f},t|x_{i},0)\right]$ (41)
$\displaystyle=\int\mathcal{D}\phi\,\,\exp\left[\frac{i}{\hbar\lambda}\int_{0}^{t}\phi^{2}(\tau)d\tau\right]G_{\rm
HO}[\\{\xi(t)\\}].$

As anticipated, Eq. (41) illustrates the fact that the propagator for the
quartic oscillator is given by an infinite collection of classical path
contributions of harmonic oscillators with different time dependent
frequencies. Moreover, Eq. (41) provides the starting point to perform the
semiclassical approximation, corresponding to the limit
$\hbar\lambda\rightarrow 0$. We start by rescaling spatial coordinates as
$y\equiv x\sqrt{\lambda}$, and we set $y_{f}\equiv\sqrt{\lambda}\,x_{f}$,
$y_{i}\equiv\sqrt{\lambda}\,x_{i}$. Due to the homogeneity of $S_{\rm HO}$
with respect to $x_{i}$ and $x_{f}$, this rescaling allows us to cast the
propagator $G$ as

$\displaystyle G(x_{f},t|x_{i},0)$
$\displaystyle=\int\mathcal{D}\phi\,\,\sqrt{\frac{m}{2\pi i\hbar
f(t)}}\exp\left[\frac{i}{\hbar\lambda}\left(\int_{0}^{t}\phi^{2}(\tau)d\tau+\,S_{\rm
HO}(y_{f},t|y_{i},0)\right)\right].$ (42)

In the limit $\hbar\lambda\rightarrow 0$, we can approximate the functional
integral by applying the stationary phase method [62]. We obtain

$\bar{\phi}(\tau)=-\frac{1}{2}\frac{\delta S_{\rm
HO}(y_{f},t;y_{i},0)}{\delta\phi(\tau)}\Big{|}_{\bar{\phi}}=\frac{\bar{y}^{2}(\tau)}{2}=\lambda\frac{\bar{x}^{2}(\tau)}{2}.$
(43)

By inserting Eq. (43) in Eq. (38) we retrieve the equation for the classical
trajectories of the quartic oscillator

$\ddot{\bar{x}}(\tau)+\omega^{2}\bar{x}(\tau)+\frac{\bar{x}^{3}(\tau)}{m}=0,$
(44)

with boundary conditions $\bar{x}(0)=x_{i}$ and $\bar{x}(t)=x_{f}$. The
solution of Eq. (44) can be expressed in terms of Jacobi elliptic functions
[42]. Note that, according to Eqs. (43) and (44), the stationary field
$\bar{\phi}$ is continuous and twice differentiable, meaning that among all
possible realizations of $\phi$ only a subset with sufficient regularity
contributes in the semiclassical limit. As reported in Ref. [33], the
associated stationary trajectories $\bar{x}$ can be classified in terms of the
sign of the momentum $\hat{p}$ of the particles $\bar{p}=m\dot{\bar{x}}$ at
the boundary points $x_{i}$ and $x_{f}$.

The semiclassical approximation is obtained by considering terms of the
expansion in $\sqrt{\hbar\lambda}$ up to the second order around the
stationary phase solution. In this spirit, we introduce the change of variable
$\phi=\bar{\phi}+\sqrt{\hbar\lambda}\,\varphi$ and truncate the expansion
around $\bar{\phi}$ at second order in $\varphi$, leading to

$\displaystyle G(x_{f},t|x_{i},0)$ $\displaystyle=\sum_{k}\sqrt{\frac{m}{2\pi
i\hbar
f_{k}(t)}}\,\,e^{\frac{i}{\hbar}S_{k}(x_{f},t|x_{i},0)}\int\mathcal{D}\varphi\exp\left[i\iint_{0}^{t}\mathrm{d}t_{1}\mathrm{d}t_{2}\,\varphi(t_{1})H_{k}(t_{1},t_{2})\varphi(t_{2})\right]$
(45)

with $\mathcal{D}\varphi\equiv\prod_{n}\mathrm{d}\phi_{n}/\sqrt{i\pi}$, where
the second functional derivative of the action computed at $\bar{\phi}$
corresponds to the operator

$H_{k}(t_{1},t_{2})\equiv\delta(t_{1}-t_{2})+\frac{1}{2}\frac{\delta^{2}\,S_{\rm
HO}(y_{f},t|y_{i},0)}{\delta\phi(t_{1})\,\delta\phi(t_{2})}\Big{|}_{\bar{\phi}_{k}};$
(46)

the zeroth order of the expansion renders the classical action
$S_{k}(x_{f},t|x_{i},0)$ of the quartic oscillator

$S_{k}(x_{f},t|x_{i},0)=\int_{0}^{t}\mathrm{d}\tau\left[\frac{m}{2}(\dot{x}_{k})^{2}-\frac{1}{2}m\omega^{2}\,x_{k}^{2}-\frac{\lambda}{4}x_{k}^{4}\right],$
(47)

evaluated on $x_{k}$, the $k$-th solution of Eq. (44) with
$\phi_{k}\equiv\bar{\phi}[x_{k}]$, and $f_{k}\equiv f[\phi_{k}]$. In order to
determine the functional Gaussian integral in Eq. (45) it is necessary to
calculate the determinant of the operator $H_{k}$. As shown in Appendix E,
this computation can be done explicitly, leading to the final expression of
the semiclassical propagator as

$G(x_{f},t|x_{i},0)=\sum_{k}\sqrt{\frac{m}{2\pi
i\hbar\,F_{k}(t)}}\,\,e^{iS_{k}(x_{f},t|x_{i},0)/\hbar},$ (48)

where $F_{k}\equiv F[\phi_{k}]$, similarly to $f_{k}$, satisfies the
differential equation

$\ddot{F}_{k}(\tau)+\left[\omega^{2}+3\frac{\lambda}{m}x_{k}^{2}(\tau)\right]F_{k}(\tau)=0,$
(49)

with $F_{k}(0)=0$ and $\dot{F}_{k}(0)=1$, being $F_{k}$ proportional to the
determinant of the $H_{k}$ operator.

In the process of deriving Eq. (47) within the present approach, we relate the
operator $H_{k}(t_{1},t_{2})$ to the functional derivative of (38) according
to Eq. (110) in Appendix E, leading to $\det(H_{k})=\sqrt{f_{k}(t)/F_{k}(t)}$.
In summary, we have derived an alternative representation of the propagator
$G(x_{f},t|x_{i},0)$ of the quartic oscillator, expressed as a weighted
collection of the propagators of effective harmonic oscillators $G_{\rm
HO}(x_{f},t|x_{i},0)$. Moreover, we provided a parametrization of the time-
evolution of these harmonic oscillators in terms of the stochastic variables.
Finally, we have proven that the semiclassical approximation of
$G(x_{f},t|x_{i},0)$ relies on the calculation of the determinant of the
second variation of the effective action $S_{\rm HO}$, and how this
determinant is linked to the density of paths along the classical trajectory
$\bar{x}_{k}$ of the quartic oscillator.

###  6.2 Partition Function

In this section we show that our formulation is not only restricted to non-
equilibrium problems, but it can be used to extract finite-temperature [15,
16] or ground state [21] properties by Wick-rotating to imaginary time. Here
we provide an additional confirmation of the validity of the stochastic
representation of the quantum quartic oscillator by obtaining the
semiclassical limit ($\hbar\rightarrow 0$) of its partition function
$Z(\beta)$, where $\beta$ is the inverse temperature. The partition function
$Z(\beta)$ is obtained from the propagator $G(x_{f},t|x_{i},0)$ in Eq. (34)
according to [62]

$\displaystyle Z(\beta)$
$\displaystyle=\int_{-\infty}^{+\infty}\mathrm{d}x\,G(x,t=-i\hbar\beta|x,0)=\int_{-\infty}^{+\infty}\mathrm{d}x\,\Biggl{\langle}\frac{1}{\sqrt{-2\pi\hbar\xi^{-}}}\exp\left[\frac{\xi^{z}}{4}+\frac{x^{2}}{2\hbar}\left(\xi^{+}+\frac{(e^{\xi^{z}/2}-1)^{2}}{\xi^{-}}\right)\right]\Biggl{\rangle}_{\phi}.$
(50)

The associated evolution of the stochastic variables as functions of $\beta$
is determined by the set of differential equations

$\frac{1}{\hbar}\frac{{\rm d}\xi^{+}}{{\rm
d}\beta}-\frac{(\xi^{+})^{2}}{m}+m\omega^{2}+2\phi=0,\qquad\frac{1}{\hbar}\frac{{\rm
d}\xi^{z}}{{\rm d}\beta}-\frac{2}{m}\xi^{+}=0,\qquad\frac{1}{\hbar}\frac{{\rm
d}\xi^{-}}{{\rm d}\beta}+\frac{e^{\xi^{z}}}{m}=0,$ (51)

with the usual initial conditions $\xi^{+,-,z}(0)=0$. Equations (51) suggest
that we may retrieve the semiclassical limit by retaining the leading-order
contributions of the series expansion of $\xi^{+,-,z}$ in integer powers of
$\hbar$, i.e.,

$\xi^{+}(\beta)=\sum_{k=1}^{\infty}\hbar^{k}f_{k}(\beta),\qquad\xi^{z}(\beta)=\sum_{k=1}^{\infty}\hbar^{k}g_{k}(\beta),\qquad\xi^{-}(\beta)=\sum_{k=1}^{\infty}\hbar^{k}l_{k}(\beta).$
(52)

By inserting these expansions in Eq. (51) and by retaining terms up to order
$O(\hbar^{3})$, we get closed-form expressions for the first coefficients. We
find that $f_{2}=g_{1}=g_{3}=l_{2}=0$ vanish and the non-vanishing
contributions read

$\begin{split}f_{1}&=-m\omega^{2}\beta-2\int_{0}^{\beta}\mathrm{d}\tau\,\phi(\tau),\\\
f_{3}&=\frac{1}{m}\int_{0}^{\beta}\mathrm{d}\tau(f_{1}(\tau))^{2},\\\
g_{2}&=\frac{2}{m}\int_{0}^{\beta}\mathrm{d}\tau f_{1}(\tau),\\\
\end{split}\qquad\qquad\begin{split}l_{1}&=-\frac{\beta}{m},\\\
l_{3}&=-\frac{2}{m^{2}}\int_{0}^{\beta}\mathrm{d}\tau\int_{0}^{\tau}\mathrm{d}sf_{1}(s).\end{split}$
(53)

It follows that all the coefficients in Eq. (53) are expressed in terms of
$f_{1}(\beta)$, so that the calculation reduces to the evaluation of moments
of $f_{1}$. We now consider the leading contribution up to order
$O(\hbar^{2})$ of Eq. (50) by explicitly inserting the expansions in Eq. (52):

$\displaystyle Z(\beta)=$
$\displaystyle\int_{-\infty}^{+\infty}\mathrm{d}x\,\Biggl{\langle}\frac{1}{\hbar}\sqrt{\frac{m}{2\pi\beta}}\exp\left(f_{1}\frac{x^{2}}{2}\right)\left\\{1+\frac{\hbar^{2}}{2}\left[\frac{g_{2}}{2}-\frac{l_{3}}{l_{1}}+\frac{x^{2}}{2}\left(f_{3}+\frac{g_{2}^{2}}{4l_{1}}\right)\right]+o(\hbar^{2})\right\\}\Biggl{\rangle}_{\phi}.$
(54)

The expectation value with respect to the Gaussian field $\phi$ is then easily
calculated and it is given by the sum of the following expressions

$\displaystyle\Big{\langle}\exp\left(f_{1}\frac{x^{2}}{2}\right)\frac{g_{2}}{2}\Big{\rangle}_{\phi}=-e^{-\beta
V(x)}\frac{\beta^{2}(m\omega^{2}+\lambda x^{2})}{2m},$ (55)
$\displaystyle\Big{\langle}\exp\left(f_{1}\frac{x^{2}}{2}\right)\frac{l_{3}}{l_{1}}\Big{\rangle}_{\phi}=-e^{-\beta
V(x)}\frac{\beta^{2}(m\omega^{2}+\lambda x^{2})}{3m},$
$\displaystyle\Big{\langle}\exp\left(f_{1}\frac{x^{2}}{2}\right)f_{3}\Big{\rangle}_{\phi}=\beta^{2}\frac{e^{-\beta
V(x)}}{m}\left[\frac{\beta(m\omega^{2}+\lambda x^{2})^{2}}{3}-\lambda\right],$
$\displaystyle\Big{\langle}\exp\left(f_{1}\frac{x^{2}}{2}\right)\frac{g_{2}^{2}}{4l_{1}}\Big{\rangle}_{\phi}=\beta^{2}\frac{e^{-\beta
V(x)}}{m}\left[-\frac{\beta(m\omega^{2}+\lambda
x^{2})^{2}}{4}+\frac{2}{3}\lambda\right],$

where $V(x)\equiv\frac{1}{2}m\omega^{2}x^{2}+\frac{\lambda}{4}x^{4}$. Given
that the average is computed with respect to a quadratic measure with zero
average, the odd moments vanish and the contributions in Eq. (55) are real-
valued. By collecting all of the above terms, we finally obtain the
semiclassical expansion of the partition function $Z(\beta)$:

$\displaystyle Z(\beta)$
$\displaystyle=\int_{-\infty}^{+\infty}\mathrm{d}x\,e^{-\beta
V(x)}\left\\{1+\frac{\hbar^{2}\beta^{2}}{12m}\left[x^{2}\beta\frac{(m\omega^{2}+\lambda
x^{2})^{2}}{2}-(m\omega^{2}+3\lambda x^{2})\right]+o(\hbar^{2})\right\\}$ (56)
$\displaystyle=\int_{-\infty}^{+\infty}\mathrm{d}x\,e^{-\beta
V(x)}\left\\{1+\frac{\hbar^{2}\beta^{2}}{12m}\left[\beta\frac{(V^{\prime}(x))^{2}}{2}-V^{\prime\prime}(x)\right]+o(\hbar^{2})\right\\},$

which matches the expression reported in the literature, see, e.g., Ref. [62].
As in the case of the propagator $G(x_{f},t|x_{i},0)$ discussed in the
previous section, we have shown how it is possible to represent the partition
function $Z(\beta)$ of the quartic oscillator in terms of the imaginary-time
version of the stochastic variables.

##  7 Summary and outlook

In this work we generalized the stochastic formalism recently introduced for
quantum spin systems [15, 16, 17, 18] to the case of non-linear bosonic
systems, explicitly considering the quantum quartic oscillator. We derived the
exact disentangled representation of the time-evolution operator of the
quartic oscillator in Eq. (14) and provided exact formulas for the time
evolution of Gaussian wave packets. In particular, we considered the time
evolution of the expectation values of the position and of the momentum
operator and their corresponding higher moments. We benchmarked our approach
(i) in the harmonic and the commuting limit by comparison with the respective
analytic solutions and (ii) for a quartic anharmonicity the comparison was
done numerically by using the stochastic interpretation of the formalism. We
further use the stochastic formalism to derive a perturbative expansion of the
time-evolution operator in powers of the quartic term. We recover the usual
Dyson series for the quantum quartic oscillator which thus implies that our
formalism is viable for evaluating perturbative expansions of observables.
Finally, we provided a semiclassical expansion of the propagator and the
partition function. Our results agree with known expressions, proving the
validity of this alternative formulation.

The stochastic approach presented in this work provides a novel theoretical
formulation of the quantum quartic oscillator as a paradigm of non-linear
bosonic systems. We described the quartic oscillator by an ensemble of
harmonic oscillators under the influence of classical (stochastic) fields.
This exact representation is a suitable starting point for developing a range
of approximations, which we illustrated in sections 5 and 6. Furthermore, the
realizations of the classical stochastic fields fully encode the underlying
quantum problem. Hence, their study should provide information about the
dynamical properties of the corresponding quantum system, as recently found
for spin systems [18]. Finally, the stochastic approach establishes a
connection between bosonic quantum systems and the theory of classical
stochastic processes. In particular, this connection allows us to evaluate the
evolution of physical observables numerically in a novel fashion, e.g., our
method does not require a truncation of the Hilbert space dimension. However,
despite this numerical strategy being quite intuitive and simple to implement,
the non-unitarity of the effective time-evolution operator leads to artificial
divergences such that simulations break down after a finite time which depends
on the strength of the quartic coupling.

Possible further directions include the generalization of our approach to
coupled oscillators [63] and bosonic lattice systems. This could be done by
decoupling interactions between different sites by means of additional
Hubbard-Stratonovich fields, as is done for quantum spin systems [15, 16, 17,
18]. The disentanglement approach could then provide a numerical technique to
simulate bosonic dynamics as well as an analytical framework based on which
further approximations can be developed.

## Acknowledgments

##### Funding information

S. De Nicola acknowledges funding from the Institute of Science and Technology
Austria (ISTA), and from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie Grant Agreement No.
754411. S. De Nicola also acknowledges funding from the EPSRC Center for
Doctoral Training in Cross-Disciplinary Approaches to NonEquilibrium Systems
(CANES) under Grant EP/L015854/1.

##  Appendix A Stochastic integral and time-dependent quartic Hamiltonian

###  A.1 Gaussian integral

In this section, we derive Eq. (4), which is fundamental for the construction
of the stochastic description of the quartic potential. First, we evaluate the
integral

$\oint_{C}\mathrm{d}z\,e^{iaz^{2}},$ (57)

where $a>0$ and the contour $C$ is displayed in Fig. 5, see, e.g., [64].

Figure 5: Integration contour $C$ of Eq. (57). $C$ consists of a circular
section of angle $\pi/4$ and radius $R$, including the arc $S$ and the segment
$D$, in addition to the segment $[0,R]$ on the real line.

Because of the absence of singularities inside the contour $C$, the residue
theorem immediately implies

$\oint_{C}\mathrm{d}z\,e^{iaz^{2}}=0,$ (58)

which results from the sum of the following contributions:

$\oint_{C}\mathrm{d}z\,e^{iaz^{2}}=\int_{0}^{R}\mathrm{d}x\,e^{iax^{2}}+\int_{S}\mathrm{d}z\,e^{iaz^{2}}+\int_{D}\mathrm{d}z\,e^{iaz^{2}}=0,$
(59)

where $S$ represents the circular arc of radius $R$, parametrized as
$z=R\,e^{i\phi}$ with $\phi\in(0,\pi/4)$, and $D$ the radial contribution
corresponding to $z=r\,e^{i\pi/4}$ with $r\in(0,R)$, in the direction shown in
Fig. 5. Finally, we are interested in the limit $R\rightarrow\infty.$ The
integral along $D$ can be expressed as

$\int_{D}\mathrm{d}z\,e^{iaz^{2}}=e^{i\pi/4}\int_{R}^{0}\mathrm{d}r\,\exp\left\\{ia\left(r\,e^{i\pi/4}\right)^{2}\right\\}=-\,e^{i\pi/4}\int_{0}^{R}\mathrm{d}r\,e^{-ar^{2}},$
(60)

coming from the change of variables $z=r\,e^{i\pi/4}.$ Similarly, the integral
along $S$ is given by

$\int_{S}\mathrm{d}z\,e^{iaz^{2}}=iR\int_{0}^{\pi/4}\mathrm{d}\phi\,\exp\left\\{iaR^{2}\left[\cos(2\phi)+i\sin(2\phi)\right]\right\\},$
(61)

and it can be shown to vanish in the limit $R\rightarrow\infty$. In fact, we
start from the inequality

$\Bigg{|}\int_{S}\mathrm{d}z\,e^{iaz^{2}}\Bigg{|}\leq
R\int_{0}^{\pi/4}\mathrm{d}\phi\,e^{-aR^{2}\sin(2\phi)}=\frac{R}{2}\int_{0}^{\pi/2}\mathrm{d}\theta\,e^{-aR^{2}\sin\theta},$
(62)

where last equality follows from introducing $\theta\equiv 2\phi.$ Moreover,
for $\theta\in(0,\pi/2)$, we have that $2\theta/\pi\leq\sin\theta\leq\theta$,
yielding

$\Bigg{|}\int_{S}\mathrm{d}z\,e^{iaz^{2}}\Bigg{|}\leq\frac{R}{2}\int_{0}^{\pi/2}\mathrm{d}\theta\,e^{-aR^{2}\sin\theta}\leq\frac{R}{2}\int_{0}^{\pi/2}\mathrm{d}\theta\,e^{-2aR^{2}\theta/\pi}=\frac{\pi}{4aR}\left(1-e^{-aR^{2}}\right),$
(63)

which vanishes for $R\rightarrow\infty.$ Finally, we get

$\int_{0}^{\infty}\mathrm{d}x\,e^{iax^{2}}=e^{i\pi/4}\int_{0}^{\infty}\mathrm{d}r\,e^{-ar^{2}}=e^{i\pi/4}\sqrt{\frac{\pi}{4a}}=\sqrt{\frac{i\pi}{4a}},$
(64)

where we have fixed $\sqrt{i}=e^{i\pi/4}.$ This last result can be generalized
to

$\int_{-\infty}^{+\infty}\mathrm{d}x\,e^{i(ax^{2}+bx)}=\sqrt{\frac{i\pi}{a}}e^{-ib^{2}/(4a)}.$
(65)

###  A.2 Time-dependent quartic Hamiltonian

In the framework discussed in section 2, it is natural to generalise the
expression of Eq. (14) to the case of the time-dependent Hamiltonian

$\hat{H}(t)=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}(t)\hat{x}^{2}+\frac{\lambda(t)}{4}\hat{x}^{4},$
(66)

with a non-negative quartic coupling, $\lambda(t)\geq 0$. The procedure
follows the same steps as in the time-independent case, i.e.,

  * (i)

we perform a Trotter-Suzuki splitting of the time-evolution operator, yielding

$\hat{U}(t)=\lim_{n\to\infty}\left\\{\exp\left[-\frac{i\tau_{n}}{\hbar}\frac{\hat{p}^{2}}{2m}\right]\exp\left[-\frac{i\tau_{n}}{\hbar}\left(\frac{m}{2}\omega_{n}^{2}\hat{x}^{2}+\frac{\lambda_{n}}{4}\hat{x}^{4}\right)\right]\right\\}^{n},$
(67)

with $\omega_{n}^{2}\equiv\omega^{2}(n\tau_{n})$ and
$\lambda_{n}\equiv\lambda(n\tau_{n})$;

  * (ii)

this is followed by the Hubbard-Stratonovich transformation, performed through
an integral of the type

$\exp\left(-i\frac{\tau_{n}}{\hbar}\frac{\lambda_{n}}{4}\hat{x}^{4}\right)=\sqrt{\frac{\tau_{n}}{i\hbar\pi}}\int_{-\infty}^{\infty}\mathrm{d}\phi\,\exp\left[\frac{i\tau_{n}}{\hbar}\left(\phi^{2}-\sqrt{\lambda_{n}}\,\hat{x}^{2}\,\phi\right)\right];$
(68)

  * (iii)

Finally, combining all previous steps, we retrieve Eq. (14) with
$\Omega^{2}(t)\equiv\omega^{2}(t)+2\sqrt{\lambda(t)}\,\phi(t)/m$ and
$S_{0}[\phi]=\hbar^{-1}\int_{0}^{t}\mathrm{d}\tau\,\phi^{2}(\tau).$

Once again, this allows one to study the dynamics of a quantum problem by
studying a set of classical differential equations.

##  Appendix B Action of $\exp(w\,\hat{S}^{\pm,z})$ on a Gaussian wave packet

We now investigate the action of the operators appearing in Eq. (12), i.e.,

$\hat{U}^{\alpha}\equiv\exp(w\,\hat{S}^{\alpha})$ (69)

with complex $w$ and $\alpha\in\\{+,-,z\\}$, on the Gaussian wave packet
$\ket{\psi}$, reported in Eq. (16). For a general complex-valued $w$, the
exponential operators are not unitary and do not conserve the normalization of
the state. For simplicity, we consider here values of $w$ for which the
corresponding $\hat{U}^{\alpha}$ conserves the state normalization. Note that
the $\hat{S}^{\alpha}$ operators in Eqs. (9) are, at most, of quadratic order
with respect to the operators $\hat{x}$ and $\hat{p}$. In what follows, we
assume $\hbar=1$ and a real-valued $w$. In order to analyse the action of
$\hat{U}^{\alpha}$ in Eq. (69), we have introduced in Eq. (17) the Wigner
function $W(x,p)$ for the generic state $\ket{\psi}$.

$W(x,p)$ provides a phase space description of the state, and allows us to
compute expectation values of operators of the type $O(\hat{x},\hat{p})$ as
$\int\mathrm{d}x\int\mathrm{d}p\,O(x,p)\,W(x,p)$ [53, 54, 55]. The evaluation
of the Wigner function in Eq. (18) for the Gaussian wave packet $\ket{\psi}$
in Eq. (16) is obtained by substitution of Eq. (16) into Eq. (17) and
integrating with respect to $y$. First, we consider the action of the operator
$\hat{U}^{+}=\exp(iw\hat{S}^{+})$ on $\ket{\psi}$, that we denote as
$\ket{\psi_{+}}\equiv U^{+}\ket{\psi}$. The state $\ket{\psi_{+}}$ is simply
given by

$\ket{\psi_{+}}=\int\frac{\mathrm{d}x}{\sqrt[4]{\pi\sigma^{2}}}\exp\left[-\frac{(x-a)^{2}}{2\sigma^{2}}+iw\frac{x^{2}}{2}+ik(x-a)\right]\ket{x},$
(70)

which follows from the fact that $\hat{S}^{+}$ acts trivially on its
eigenstate $\ket{x}$.

By substituting Eq. (70) into (17) we get the Wigner function for
$\ket{\psi_{+}}$

$W_{+}(x,p)=\frac{1}{\pi}\exp\left[-\frac{(x-a)^{2}}{\sigma^{2}}-\sigma^{2}(p-k-
wx)^{2}\right],$ (71)

which is equal to $W$ up to position-dependent shift in the momentum.
Accordingly, the expectation value of operators of the form $O(\hat{x})$ is
unaffected by the $\hat{U}^{+}$ transformation, i.e.,
$\bra{\psi}O(\hat{x})\ket{\psi}=\bra{\psi_{+}}O(\hat{x})\ket{\psi_{+}}$. On
the other hand, the expectation value of a momentum-dependent operator
$O(\hat{p})$ can be expressed as

$\displaystyle\bra{\psi_{+}}O(\hat{p})\ket{\psi_{+}}$
$\displaystyle=\int\mathrm{d}x\int\mathrm{d}p\,O(p)\,W_{+}(x,p)$ (72)
$\displaystyle=\sqrt{\frac{\sigma^{2}}{\pi(1+\sigma^{4}w^{2})}}\int\mathrm{d}p\,O(p)\exp\left[-\sigma^{2}\frac{(p-k-
wa)^{2}}{1+\sigma^{4}w^{2}}\right]$
$\displaystyle=\sqrt{\frac{\sigma^{2}}{\pi}}\int\mathrm{d}q\,O\left(\sqrt{1+\sigma^{4}w^{2}}\,(q-k)+wa+k\right)\,e^{-\sigma^{2}(q-k)^{2}}$
$\displaystyle=\bra{\psi}O\left(\sqrt{1+\sigma^{4}w^{2}}\,(\hat{p}-k)+wa+k\right)\ket{\psi},$

where the second line comes from direct substitution of Eq. (71) and the final
result from the change of variable $p-wa-k=(q-k)\,\sqrt{1+\sigma^{4}w^{2}}$.
Equation (72) tells us that expectation values with respect to the state
$\ket{\psi_{+}}$ of operators depending only on $p$ are equivalent to
expectation values with respect to the Gaussian wave packet $\ket{\psi}$ with
the rescaled and shifted momentum operator
$\sqrt{1+\sigma^{4}w^{2}}\,(\hat{p}-k)+wa+k$. In particular, for the mean and
the variance of the momentum operator we can immediately read off from Eq.
(72) that

$\displaystyle\langle p\rangle$
$\displaystyle\equiv\bra{\psi_{+}}\hat{p}\ket{\psi_{+}}=k+wa,$ (73)
$\displaystyle\langle p^{2}\rangle_{c}$
$\displaystyle\equiv\bra{\psi_{+}}\hat{p}^{2}\ket{\psi_{+}}-\bra{\psi_{+}}\hat{p}\ket{\psi_{+}}^{2}=\frac{1+\sigma^{4}w^{2}}{2\sigma^{2}}.$

These parameters, together with the unaltered cumulants of the position
operator, allow us to fully characterize the state $\ket{\psi_{+}}$. As an
explicit time-dependent example, we consider the evolution of the wave packet
under the action of the harmonic oscillator Hamiltonian. Referring to Eqs.
(25), we find that

$\displaystyle w(t)$ $\displaystyle=-m\omega\tan(\omega t),$ (74)
$\displaystyle\langle p(t)\rangle$ $\displaystyle=k-am\omega\tan(\omega t),$
$\displaystyle\langle p^{2}(t)\rangle_{c}$
$\displaystyle=\frac{1+(m\omega\sigma^{2})^{2}\tan^{2}(\omega
t)}{2\sigma^{2}}.$

Next, we consider the case of the operator
$\hat{U}^{z}=\exp\left(iw\\{\hat{x},\hat{p}\\}/4\right)$ whose action on the
Gaussian wave packet $\ket{\psi}$, which we denote by
$\ket{\psi_{z}}\equiv\hat{U}^{z}\ket{\psi}$, reads

$\ket{\psi_{z}}=\int\frac{\mathrm{d}x}{\sqrt[4]{\pi\sigma^{2}}}\exp\left[\frac{w}{4}-\frac{(x\,e^{w/2}-a)^{2}}{2\sigma^{2}}+ik(x\,e^{w/2}-a)\right]\ket{x}.$
(75)

This is computed considering the direct action of $\hat{S}^{z}$ on $\ket{x}$
according the property of the dilation operator
$e^{by\frac{d}{dy}}f(y)=f(e^{b}y)$, where $f$ is any sufficiently smooth
function, similarly to what has been done for Eq. (86). The Wigner function
$W_{z}(x,p)$ of the state $\ket{\psi_{z}}$ can be directly evaluated as

$W_{z}(x,p)=\frac{1}{\pi}\exp\left[-\frac{(x\,e^{w/2}-a)^{2}}{\sigma^{2}}-\sigma^{2}(p\,e^{-w/2}-k)^{2}\right].$
(76)

Equation (75) shows that the action of $\hat{U}^{z}$ consists in a uniform
rescaling all the $x$ variables by a factor $e^{w/2}$.

As for $\ket{\psi_{z}}$, the Wigner function $W_{z}$ is equivalent to $W$ up
to a rescaling of the variables. It follows that the expectation value of an
operator $O(\hat{x})$, depending only on $x$, is given by

$\displaystyle\bra{\psi_{z}}O(\hat{x})\ket{\psi_{z}}$
$\displaystyle=\int\mathrm{d}x\int\mathrm{d}p\,O(x)W_{z}(x,p)$ (77)
$\displaystyle=\frac{e^{w/2}}{\sqrt{\pi\sigma^{2}}}\int\mathrm{d}x\,O(x)\exp\left[-\frac{(x\,e^{w/2}-a)^{2}}{\sigma^{2}}\right]$
$\displaystyle=\bra{\psi}O(\hat{x}\,e^{-w/2})\ket{\psi},$

and, analogously, for a $p-$dependent operator $O(\hat{p})$, we get

$\bra{\psi_{z}}O(\hat{p})\ket{\psi_{z}}=\bra{\psi}O(\hat{p}\,e^{w/2})\ket{\psi},$
(78)

which reflects the rescaling action of $\hat{U}^{z}$. It follows that the
first connected moments of $\hat{x}$ and $\hat{p}$ on $\ket{\psi_{z}}$ are
given by

$\displaystyle\langle x\rangle$
$\displaystyle\equiv\bra{\psi_{z}}\hat{x}\ket{\psi_{z}}=a\,e^{-w/2},$ (79)
$\displaystyle\langle x^{2}\rangle_{c}$
$\displaystyle\equiv\bra{\psi_{z}}\hat{x}^{2}\ket{\psi_{z}}-\bra{\psi_{z}}\hat{x}\ket{\psi_{z}}^{2}=\frac{\sigma^{2}}{2}\,e^{-w},$
$\displaystyle\langle p\rangle$
$\displaystyle\equiv\bra{\psi_{z}}\hat{p}\ket{\psi_{z}}=k\,e^{w/2},$
$\displaystyle\langle p^{2}\rangle_{c}$
$\displaystyle\equiv\bra{\psi_{z}}\hat{p}^{2}\ket{\psi_{z}}-\bra{\psi_{z}}\hat{p}\ket{\psi_{z}}^{2}=\frac{e^{w}}{2\sigma^{2}}.$

These quantities fully characterize the state $\ket{\psi_{z}}$. According to
Eqs. (25), for a $\ket{\psi}$ evolving under the effect of an harmonic
oscillator Hamiltonian, we have

$\displaystyle w(t)$ $\displaystyle=-\log\cos^{2}(\omega t),$ (80)
$\displaystyle\langle x(t)\rangle$ $\displaystyle=a\cos(\omega
t),\,\,\,\,\langle p\rangle(t)=\frac{k}{\cos(\omega t)},$
$\displaystyle\langle x^{2}(t)\rangle_{c}$
$\displaystyle=\frac{\sigma^{2}}{2}\,\cos^{2}(\omega t),\,\,\,\,\langle
p^{2}(t)\rangle_{c}=\frac{1}{2\sigma^{2}\cos^{2}(\omega t)}.$

Finally, we consider the action of the operator
$\hat{U}^{-}=\exp\left(iw\hat{p}^{2}/2\right)$ on the wave packet
$\ket{\psi}$; the resulting state $\ket{\psi_{-}}\equiv\hat{U}^{-}\ket{\psi}$
is found to be

$\ket{\psi_{-}}=\sqrt{\frac{\sigma^{2}}{\sigma^{2}-iw}}\int\frac{\mathrm{d}x}{(\pi\sigma^{2})^{1/4}}\,\exp\left[-\frac{(x-a-
ik\sigma^{2})^{2}}{2(\sigma^{2}-iw)}-\frac{\sigma^{2}k^{2}}{2}\right]\ket{x}.$
(81)

The associated Wigner function $W_{-}(x,p)$ reads

$W_{-}(x,p)=\frac{1}{\pi}\exp\left[-\sigma^{2}(p-k)^{2}-\frac{(x-a+pw)^{2}}{\sigma^{2}}\right],$
(82)

that is equivalent to $W$ up to a $p-$dependent rescaling of the $x$ variable.
In this case expectation values of $p-$dependent operators are invariant under
the action of $\hat{U}^{-}$, i.e.,
$\bra{\psi_{-}}O(\hat{p})\ket{\psi_{-}}=\bra{\psi}O(\hat{p})\ket{\psi}$, while
the expectation value of an $x-$dependent operator $O(\hat{x})$ transforms as

$\displaystyle\bra{\psi_{-}}O(\hat{x})\ket{\psi_{-}}$
$\displaystyle=\int\mathrm{d}x\int\mathrm{d}p\,O(x)\,W_{-}(x,p)$ (83)
$\displaystyle=\sqrt{\frac{\sigma^{2}}{\pi(\sigma^{4}+w^{2})}}\int\mathrm{d}x\,O(x)\exp\left[-\sigma^{2}\frac{(x-a+wk)^{2}}{\sigma^{4}+w^{2}}\right]$
$\displaystyle=\frac{1}{\sqrt{\sigma^{2}\pi}}\int\mathrm{d}y\,O\left(\frac{y-a}{\sigma^{2}}\,\sqrt{\sigma^{4}+w^{2}}+a-wk\right)\,e^{-(y-a)^{2}/\sigma^{2}}$
$\displaystyle=\bra{\psi}O\left(\frac{\hat{x}-a}{\sigma^{2}}\,\sqrt{\sigma^{4}+w^{2}}+a-wk\right)\ket{\psi},$

where the second line is found by integrating $W_{-}(x,p)$ in Eq. (82) with
respect to $p$, and the last two lines are obtained by performing the change
of variable $x=(y-a)\sigma^{-2}\,\sqrt{\sigma^{4}+w^{2}}+a-wk$. We deduce
that, in case of $x-$dependent operators, the expectation value with respect
to $\ket{\psi_{-}}$ is equivalent to the expectation value with respect to
$\ket{\psi}$ where the position operator has been rescaled and shifted
according to the final line of Eq. (83). In particular, the first two
cumulants of the position operator $\hat{x}$ read

$\displaystyle\langle x\rangle$
$\displaystyle\equiv\bra{\psi_{-}}\hat{x}\ket{\psi_{-}}=a-wk,$ (84)
$\displaystyle\langle x^{2}\rangle_{c}$
$\displaystyle\equiv\bra{\psi_{-}}\hat{x}\ket{\psi_{-}}-\bra{\psi_{-}}\hat{x}^{2}\ket{\psi_{-}}^{2}=\frac{\sigma^{4}+w^{2}}{2\sigma^{2}}.$

The evolution of the state $\ket{\psi_{-}}$ under the harmonic oscillator
dynamics can be explicitly determined from Eq. (25):

$\displaystyle w(t)$ $\displaystyle=-\frac{\tan(\omega t)}{m\omega},$ (85)
$\displaystyle\langle x(t)\rangle$ $\displaystyle=a+k\frac{\tan(\omega
t)}{m\omega},$ $\displaystyle\langle x^{2}(t)\rangle_{c}$
$\displaystyle=\frac{(m\omega\sigma^{2})^{2}+\tan^{2}(\omega
t)}{2(m\omega\sigma)^{2}}.$

As a last remark, we note that by combining the results in Eqs. (74), (80),
and (85) into the factorised expression for the time evolution of the harmonic
oscillator, Eq. (12), we can construct the time evolution of observable, e.g.,
the position and momentum moments in Eqs. (26).

##  Appendix C Time evolution of a Gaussian wave packet

Here we report the detailed calculations of the expectation values of the
moments of the position and momentum operators on the Gaussian wave packet in
Eq. (16). As a preliminary step to the calculation of Eq. (19), we consider
the action of the operator $\hat{U}(t)$ on an eigenstate $\ket{x}$ of the
position operator, given by

$\displaystyle\hat{U}(t)\ket{x}$
$\displaystyle=\Biggl{\langle}\int\frac{\mathrm{d}p}{\sqrt{2\pi}}\,\,e^{\xi^{+}(t)\hat{x}^{2}/2}\,e^{i\xi^{z}(t)\\{\hat{x},\hat{p}\\}/4}\,e^{\xi^{-}(t)p^{2}/2-ipx}\ket{p}\Biggl{\rangle}_{\phi}$
(86)
$\displaystyle=\Biggl{\langle}e^{-\xi^{z}/4}\int\frac{\mathrm{d}p}{\sqrt{2\pi}}\,\,e^{\xi^{+}(t)\hat{x}^{2}/2}\,e^{-(\xi^{z}(t)/2)p\frac{\partial}{\partial
p}}\,e^{\xi^{-}(t)\,p^{2}/2-ipx}\ket{p}\Biggl{\rangle}_{\phi}$
$\displaystyle=\Biggl{\langle}e^{-\xi^{z}/4}\int\frac{\mathrm{d}p}{\sqrt{2\pi}}\,\,e^{\xi^{+}(t)\hat{x}^{2}/2}\,\exp{\left[\xi^{-}(t)\frac{p^{2}}{2}e^{-\xi^{z}}-ipxe^{-\xi^{z}/2}\right]}\ket{p}\Biggl{\rangle}_{\phi}$
$\displaystyle=\Biggl{\langle}e^{-\xi^{z}/4}\int\frac{\mathrm{d}p}{\sqrt{2\pi}}\,\int\mathrm{d}y\,\,e^{\xi^{+}(t)y^{2}/2}\ket{y}\bra{y}\exp{\left[\xi^{-}(t)\frac{p^{2}}{2}e^{-\xi^{z}}-ipxe^{-\xi^{z}/2}\right]}\ket{p}\Biggl{\rangle}_{\phi}$
$\displaystyle=\Biggl{\langle}e^{-\xi^{z}/4}\int\frac{\mathrm{d}y}{\sqrt{2\pi}}\,\,e^{\xi^{+}(t)y^{2}/2}\int\frac{\mathrm{d}p}{\sqrt{2\pi}}\,\exp{\left[\xi^{-}(t)\frac{p^{2}}{2}e^{-\xi^{z}}-ipxe^{-\xi^{z}/2}+ipy\right]}\ket{y}\Biggl{\rangle}_{\phi}$
$\displaystyle=\Biggl{\langle}\frac{\exp\left[\xi^{z}/4+x^{2}/(2\xi^{-})\right]}{\sqrt{-2\pi\xi^{-}}}\int\mathrm{d}y\,\,\exp\left[\frac{y^{2}}{2}\left(\frac{e^{\xi^{z}}}{\xi^{-}}+\xi^{+}(t)\right)-y\frac{xe^{\xi^{z}/2}}{\xi^{-}}\right]\ket{y}\Biggl{\rangle}_{\phi}$

where $\langle\dots\rangle_{\phi}$ denotes the expectation value with respect
to the Gaussian action $S_{0}$.

In the first line a completeness relation for the momentum basis, $\int
dp\ket{p}\bra{p}=\mathbb{I}$, was inserted between the last exponential
operator and the position eigenket, leading to the appearance of the plane
wave $\braket{p}{x}=e^{-ipx}/\sqrt{2\pi}$, where $\hbar=1$. In the second
line, we substituted $\hat{x}\ket{p}=i\frac{\partial}{\partial p}\ket{p}$ and
the consequent action of the dilation operator was written explicitly, i.e.,
$e^{by\frac{d}{dy}}f(y)=f(e^{b}y)$. Finally, a further position completeness
relation insertion and a Gaussian integration was performed. The convergence
of the Gaussian integral is ensured by the fact that the argument of the
exponential is purely imaginary. Analogously, the corresponding dual vector
evolves according to

$\bra{x}U^{\dagger}(t)=\Biggl{\langle}\frac{\exp\left[\bar{\xi}^{z}/4+x^{2}/(2\bar{\xi}^{-})\right]}{\sqrt{-2\pi\bar{\xi}^{-}}}\int\mathrm{d}z\,\exp\left[\frac{z^{2}}{2}\left(\frac{e^{\bar{\xi^{z}}}}{\bar{\xi}^{-}}+\bar{\xi}^{+}(t)\right)-z\frac{xe^{\bar{\xi}^{z}/2}}{\bar{\xi}^{-}}\right]\bra{z}\Biggl{\rangle}_{\bar{\phi}},$
(87)

with $\langle\dots\rangle_{\bar{\phi}}$ denoting the expectation value with
respect to the Gaussian action $S_{0}[\bar{\phi}]$, and
$\bar{\xi}\equiv[\xi(\bar{\phi}(t))]^{*}$ the complex conjugate of
$\xi^{+,-,z}$. Finally, Eq. (19) follows by plugging Eq. (86) into Eq. (16)
and integrating the Gaussian integral with respect to $x$.

Finally, the evolution of the wave packet, reported in Eq. (19), is eventually
computed by integrating the expression of $\hat{U}(t)\ket{x}$ over the
variable $x$ with respect to the Gaussian measure

$\exp\left[-(x-a)^{2}/(2\sigma^{2})^{2}+i(x-a)k\right](\pi\sigma^{2})^{-1/4}.$
(88)

The convergence is ensured by requiring that ${\operatorname{Re}(\gamma)>0}$.
In our description we have $\xi^{\pm}\in i\mathbb{R}$ and real $\xi^{z}$, so
that it is useful to define $\xi^{\pm}\equiv i\xi^{\pm}_{i}$ with
$\xi^{\pm}_{i}\in\mathbb{R}$, which, together with Eq. (20), leads to

$\operatorname{Re}(\gamma)=\frac{\sigma^{2}e^{\xi^{z}}}{\sigma^{4}+(\xi^{-}_{i})^{2}}\geq
0.$ (89)

A better understanding of the behavior of $\operatorname{Re}(\gamma)$ can be
achieved by considering the following real-valued auxiliary variables:

$\displaystyle X$ $\displaystyle\equiv\xi^{-}_{i}e^{-\xi^{z}/2},$ (90)
$\displaystyle Y$ $\displaystyle\equiv e^{-\xi^{z}/2}.$

These variables evolve according to the harmonic equations with time-dependent
frequency given in Eq. (6),

$\displaystyle\ddot{X}(t)+\Omega^{2}(t)X(t)$ $\displaystyle=0,$ (91)
$\displaystyle\ddot{Y}(t)+\Omega^{2}(t)Y(t)$ $\displaystyle=0,$

with initial conditions $X(0)=0$, $\dot{X}(0)=-m^{-1}$, $Y(0)=1$ and
$\dot{Y}(0)=0.$ These new variables allow one to write
$\operatorname{Re}(\gamma)=\sigma^{2}(\sigma^{4}Y^{2}+X^{2})^{-1}$, making it
apparent that $\operatorname{Re}(\gamma(t))=0$ if $X(t)$ or $Y(t)$ are
infinite. In either case, $\beta(t)=[\sigma^{2}Y(t)-iX(t)]^{-1}$ which
multiplies Eq. (19), vanishes, i.e., $\ket{\psi(t)}=0$. Accordingly, the
convergence of the Gaussian integral in Eq. (19) is guaranteed by the fact
that $\operatorname{Re}(\gamma)\geq 0$ and that whenever
$\operatorname{Re}(\gamma)=0$ the whole $\ket{\psi}$ vanishes.

As a last remark, we point out that the introduction of the variables $X$ and
$Y$ in Eq. (90) explains why in the case of the harmonic oscillator, in which
the $\xi^{+,-,z}$ are found to be periodically divergent according to the Eqs.
(25), there are no divergences in the expectation values of Eqs. (26). In
fact, these values depend on a well-behaved combination of the $\xi^{+,-,z}$,
satisfying an harmonic equation with constant frequency but different initial
conditions.

##  Appendix D Derivations of the Dyson Series

In order to prove the equivalence of the Dyson series for the quantum quartic
oscillator and the asymptotic expansion of $\hat{U}(t)$ around the harmonic
case according to Eq. (32), we begin by calculating the second variation of
the time-evolution operator of the harmonic oscillator, which will be useful
to determine the first-order correction according to Eq. (32), namely

$\displaystyle\hat{U}^{(1)}(t)$ $\displaystyle\equiv
i\frac{\lambda}{4}\int_{0}^{t}\mathrm{d}s\frac{\delta^{2}\hat{U}_{S}[\phi]}{\delta\phi(s_{1})\delta\phi(s_{2})}\Bigg{|}_{\begin{subarray}{c}\phi=0\\\
s_{1}=s_{2}=s\end{subarray}}.$ (92)

This expression involves the first functional derivative of $\hat{U}_{S}$,
given by

$\displaystyle\frac{\delta\hat{U}_{S}[\phi]}{\delta\phi(s)}\equiv
G^{(1)}(s|t)\hat{U}_{S}[\phi],$ (93)

where $G^{(1)}$ is explicitly computed by exploiting the $SU(2)$ commutation
relations (10) of the $S$ operators and $\hat{U}_{S}$, leading to

$G^{(1)}(s)\equiv\left[\frac{\delta\xi^{+}}{\delta\phi\,\,}-\xi^{+}\frac{\delta\xi^{z}}{\delta\phi\,\,}-(\xi^{+})^{2}e^{-\xi^{z}}\frac{\delta\xi^{-}}{\delta\phi\,\,}\right]\hat{S}^{+}+\left(\frac{\delta\xi^{z}}{\delta\phi\,\,}+2\xi^{+}e^{-\xi^{z}}\frac{\delta\xi^{-}}{\delta\phi\,\,}\right)\hat{S}^{z}+\frac{\delta\xi^{-}}{\delta\phi\,\,}e^{-\xi^{z}}\hat{S}^{-},$
(94)

where the parametric dependence on the final time $t$ is understood.

Similarly, the second order functional derivative
$\delta^{2}\hat{U}_{S}[\phi]/\delta\phi(s_{1})\delta\phi(s_{2})$, required to
be symmetric under the exchange $s_{1}\leftrightarrow s_{2}$, can be expressed
as

$\displaystyle\frac{\delta^{2}\hat{U}_{S}[\phi]}{\delta\phi(s_{1})\delta\phi(s_{2})}=\left[G^{(1)}(s_{1})G^{(1)}(s_{2})+G^{(2)}(s_{1},s_{2})\right]\hat{U}_{S}[\phi],$
(95)

where $G^{(2)}(s_{1},s_{2})$ is found to be

$\displaystyle G^{(2)}\equiv$
$\displaystyle\left\\{\xi^{+}_{1,2}-\xi^{+}\xi^{z}_{1,2}-\frac{1}{2}\left(\xi^{z}_{1}\xi^{+}_{2}+\xi^{z}_{2}\xi^{+}_{1}\right)-\xi^{+}e^{-\xi^{z}}\left[\xi^{-}_{1}\xi^{+}_{2}+\xi^{-}_{2}\xi^{+}_{1}-\frac{\xi^{+}}{2}\left(\xi^{z}_{1}\xi^{-}_{2}+\xi^{z}_{2}\xi^{-}_{1}\right)+\xi^{+}\xi^{-}_{1,2}\right]\right\\}\hat{S}^{+}$
(96)
$\displaystyle+e^{-\xi^{z}}\left[\xi^{-}_{1,2}-\frac{1}{2}\left(\xi^{z}_{1}\xi^{-}_{2}+\xi^{z}_{2}\xi^{-}_{1}\right)\right]\hat{S}^{-}$
$\displaystyle+\left\\{\xi^{z}_{1,2}+e^{-\xi^{z}}\left[\xi^{-}_{1}\xi^{+}_{2}+\xi^{-}_{2}\xi^{+}_{1}-\xi^{+}\left(\xi^{z}_{1}\xi^{-}_{2}+\xi^{z}_{2}\xi^{-}_{1}\right)+2\xi^{+}\xi^{-}_{1,2}\right]\right\\}\hat{S}^{z};$

in order to streamline the formulas, the subscripts $\\{1,2\\}$ above are used
to denote the functional differentiation with respect to $\phi(s_{1})$ and
$\phi(s_{2})$, i.e.,
$\delta\xi(t)/\delta\phi(s_{1})\equiv\xi(s_{1}|t)=\xi_{1}$. By taking the
functional derivative of Eqs. (13) we obtain a system of differential
equations for the first functional derivatives $\xi_{1,2}$, namely

$\displaystyle
i\frac{\mathrm{d}}{\mathrm{d}t}\xi^{+}(s|t)+2\frac{\xi^{+}}{m}\xi^{+}(s|t)$
$\displaystyle=2\delta(t-s),$ (97) $\displaystyle
i\frac{\mathrm{d}}{\mathrm{d}t}\xi^{z}(s|t)+\frac{2}{m}\xi^{+}(s|t)$
$\displaystyle=0,$ $\displaystyle
i\frac{\mathrm{d}}{\mathrm{d}t}\xi^{-}(s|t)-\frac{e^{\xi^{z}}}{m}\xi^{z}(s|t)$
$\displaystyle=0,$

with initial conditions $\xi(s|s)=0$ for $t\leq s$, reflecting the fact that
we assume an Itô-like discretization in deriving Eqs. (13)[20]. The solution
to these equations reads

$\displaystyle\xi^{+}(s|t)$
$\displaystyle=-\theta(t-s)\,2i\,\exp\left\\{\frac{2i}{m}\int_{s}^{t}\,\mathrm{d}\tau\,\xi^{+}(\tau)\right\\},$
(98) $\displaystyle\xi^{z}(s|t)$
$\displaystyle=\theta(t-s)\frac{2i}{m}\int_{s}^{t}\,\mathrm{d}\tau\,\xi^{+}(s|\tau),$
$\displaystyle\xi^{-}(s|t)$
$\displaystyle=-\theta(t-s)\frac{i}{m}\int_{s}^{t}\,\mathrm{d}\tau\,\xi^{z}(s|\tau)\,e^{\xi^{z}(\tau)}.$

For $\phi=0$ they reduce to

$\displaystyle\frac{\delta\xi^{+}(t)}{\delta\phi(s)}\Bigg{|}_{\phi=0}$
$\displaystyle=-i\theta(t-s)\,2\frac{\cos^{2}(\omega s)}{\cos^{2}(\omega t)},$
(99) $\displaystyle\frac{\delta\xi^{z}(t)}{\delta\phi(s)}\Bigg{|}_{\phi=0}$
$\displaystyle=\theta(t-s)\frac{4}{m\omega}\cos^{2}(\omega s)\left[\tan(\omega
t)-\tan(\omega s)\right],$
$\displaystyle\frac{\delta\xi^{-}(t)}{\delta\phi(s)}\Bigg{|}_{\phi=0}$
$\displaystyle=-i\theta(t-s)\frac{2}{m^{2}\omega^{2}}\cos^{2}(\omega
s)\left[\tan(\omega t)-\tan(\omega s)\right]^{2}.$

As expected, the functional derivative $\xi^{+}(s|t)$ vanishes for $t<s$, as a
consequence of the fact that the differential equation at time $t$ does not
depend on the realizations of $\phi$ at later times, reflecting the causality
of the problem. Following the same line of reasoning as before, the second
functional derivatives can be computed directly from their differential
equations and can be expressed in terms of first functional derivative
according to

$\displaystyle\xi^{+}_{1,2}$
$\displaystyle=\frac{\xi^{z}_{1}\xi^{+}_{2}+\xi^{z}_{2}\xi^{+}_{1}}{2},$ (100)
$\displaystyle\xi^{z}_{1,2}$
$\displaystyle=-e^{-\xi^{z}(t)}\left(\xi^{-}_{1}\xi^{+}_{2}+\xi^{-}_{2}\xi^{+}_{1}\right),$
$\displaystyle\xi^{-}_{1,2}$
$\displaystyle=\frac{\xi^{z}_{1}\xi^{-}_{2}+\xi^{z}_{2}\xi^{-}_{1}}{2},$

and they are non-zero only if $t>\max(s_{1},s_{2}).$ Moreover, by plugging
Eqs. (100) into Eq. (96), we get $G^{(2)}(s_{1},s_{2})=0$, such that the only
contributing term in the functional derivative in Eq. (95) is

$\displaystyle G^{(1)}(s)|_{\phi=0}$
$\displaystyle=-i\theta(t-s)\,x_{0}^{2}\,\left[\cos(\omega(t-s))\frac{\hat{x}}{x_{0}}-\sin(\omega(t-s))x_{0}\,\hat{p}\right]^{2},$
(101)

which finally yields

$\displaystyle\frac{\delta^{2}\hat{U}_{S}[\phi]}{\delta\phi(s_{1})\delta\phi(s_{2})}\Bigg{|}_{\begin{subarray}{c}\phi=0,\\\
s_{1}=s_{2}=s\end{subarray}}$
$\displaystyle=-\theta(t-s)\,x_{0}^{4}\,\left[\cos(\omega(t-s))\frac{\hat{x}}{x_{0}}-\sin(\omega(t-s))x_{0}\,\hat{p}\right]^{4}\hat{U}_{0}(t)$
(102)
$\displaystyle=-\theta(t-s)\hat{U}_{0}(t-s)\,\hat{x}^{4}\,\hat{U}_{0}(s),$

where the time ordering $t>s$ arises naturally from the fact that the equation
for $\xi^{+}$ depends linearly on $\phi$. Collecting the above results of Eqs.
(92) and (102), the first-order correction to $\hat{U}(t)$ reads

$\hat{U}^{(1)}=-i\frac{\lambda}{4}\,\hat{U}_{0}(t)\int_{0}^{t}\mathrm{d}s\,\hat{U}_{0}^{\dagger}(s)\,\hat{x}^{4}\,\hat{U}_{0}(s),$
(103)

which is nothing but the first order term in the Dyson series [56]. This can
be seen by noticing that the time-evolution operator in the Schrödinger
picture $\hat{U}(t)$ can be written in terms of the interaction time-evolution
operator in the interaction picture $\hat{U}_{I}$ as

$\hat{U}_{I}(t)\equiv\hat{U}_{0}^{\dagger}(t)\,\hat{U}(t)\,\hat{U}_{0}(0),$
(104)

so that, since $\hat{U}_{0}(0)=\mathbb{I}$,

$\hat{U}(t)=\hat{U}_{0}(t)\,\hat{U}_{I}(t).$ (105)

The fact that $G^{(2)}=0$ makes all functional derivatives of order larger
than one to depend only on $G^{(1)}$. This allows one to easily generalize the
result above to an arbitrary order $n$, leading to

$\frac{\delta^{2n}\hat{U}_{S}[\phi]}{\delta\phi(s_{1})\cdots\delta\phi(s_{2n})}\Bigg{|}_{\begin{subarray}{c}\phi=0\\\
s_{1}=s_{2}=t_{1}\\\ \cdots\\\
s_{2n-1}=s_{2n}=t_{n}\end{subarray}}=[G^{(1)}(t_{1})]^{2}[G^{(1)}(t_{2})]^{2}\cdots[G^{(1)}(t_{n})]^{2}\hat{U}_{0}(t),$
(106)

with $t>t_{n}>t_{n-1}>\dots>0$. Using the expression for $G^{(n)}$, this
readily yields

$\hat{U}^{(n)}(t)=\left(-i\frac{\lambda}{4}\right)^{n}\hat{U}_{0}(t)\int_{0}^{t}\,\mathrm{d}t_{n}\int_{0}^{t_{n}}\mathrm{d}t_{n-1}\cdots\int_{0}^{t_{2}}\mathrm{d}t_{1}\,\hat{x}^{4}(t_{n})\cdots\hat{x}^{4}(t_{1}),$
(107)

where $\hat{x}^{4}(t)=\hat{U}_{0}^{\dagger}(t)\,\hat{x}^{4}\,\hat{U}_{0}(t)$.
Equation (107) is precisely the $n$-th order contribution to the Dyson series
in the Schrödinger picture.

##  Appendix E Semiclassical limit

In this section we provide details of the computation of the semiclassical
approximation for the propagator. The stationary path in Eq. (43) is computed
considering the first functional derivative of Eq. (36), namely

$\displaystyle\frac{\delta S_{\rm HO}(y_{f},t|y_{i},0)}{\delta\phi(\tau)}$
$\displaystyle=\frac{\delta\,\,\,\,\,\,\,\,\,\,\,}{\delta\phi(\tau)}\frac{m}{2}\int_{0}^{t}\mathrm{d}s\,\left[\dot{y}^{2}(s)-\Omega^{2}(s)y^{2}(s)\right]$
(108)
$\displaystyle=-y^{2}(\tau)+m\int_{0}^{t}\mathrm{d}s\,\left[\dot{y}(s)\dot{y}_{1}(\tau|s)-\Omega^{2}(s)y(s)y_{1}(\tau|s)\right]$
$\displaystyle=-y^{2}(\tau)+m\left[\dot{y}(t)y_{1}(\tau|t)-\dot{y}(0)y_{1}(\tau|0)\right]-m\int_{0}^{t}\mathrm{d}s\,y_{1}(\tau|s)\left[\ddot{y}(s)+\Omega^{2}(s)y(s)\right]$
$\displaystyle=-y^{2}(\tau),$

where $y_{1}(\tau|s)\equiv\delta y(t)/\delta\phi(\tau)$,
$y_{f}=\sqrt{\lambda}\,x_{f}$, $y_{i}=\sqrt{\lambda}\,x_{i}$, and the $t$
dependence is understood. In the third line of Eq. (108) we integrated by
parts the right-hand side and finally we exploited Eq. (38) and the explicit
expression of $y_{1}(\tau|s)$

$\displaystyle y_{1}(\tau|s)=$
$\displaystyle-\left[\frac{\xi^{-}_{1}(\tau|t)}{\xi^{-}(t)}+\frac{\xi^{z}_{1}(\tau|s)}{2}\right]y(\tau|t)$
(109)
$\displaystyle+\frac{e^{-\xi^{z}(s)/2}}{\xi^{-}(t)}\left[y_{f}\,e^{-\xi^{z}(t)/2}\left(\frac{\xi^{-}(\tau)}{2}\xi^{z}_{1}(\tau|t)+\xi^{-}_{1}(\tau|s)\right)+y_{i}\,\left(\xi^{-}_{1}(\tau|t)-\xi^{-}_{1}(\tau|s)\right)\right],$

where $\xi_{1}(t_{1}|t_{2})$ is null if $t_{1}\geq t_{2}$, so that
$y_{1}(\tau|t)=y_{1}(\tau|0)=0.$ Note that by taking the functional derivative
of Eq. (38) computed along the stationary solution $\bar{\phi}$ in Eq. (43)
one has

$\ddot{\bar{y}}_{1}(t_{2}|t_{1})+\left[\omega^{2}+\frac{\bar{y}^{2}(t_{1})}{m}\right]\bar{y}_{1}(t_{2}|t_{1})+\frac{2}{m}\bar{y}(t_{1})\delta(t_{1}-t_{2})=0,$
(110)

with boundary conditions $\bar{y}_{1}(t_{2}|0)=\bar{y}_{1}(t_{2}|0)=0.$ The
path integral in the second line of Eq. (45) can be expressed in terms of the
determinant of the operator $H(t_{1},t_{2})$ , defined in Eq. (46) as

$H(t_{1},t_{2})\equiv\delta(t_{1}-t_{2})+\frac{1}{2}\frac{\delta^{2}\,S_{\rm
HO}(y_{f},t|y_{i},0)}{\delta\phi(t_{1})\delta\phi(t_{2})}\Big{|}_{\bar{\phi}_{k}}=\delta(t_{1}-t_{2})-\bar{y}(t_{1})\bar{y}_{1}(t_{2}|t_{1}),$
(111)

where last equality follows from direct functional derivation of Eq. (108).
The determinant of $H$ can be evaluated by relying on the fact that this
operator can be recast as the product of two operators whose determinant can
be computed exactly. We start by defining the operators

$\displaystyle O_{1}(t_{1},t_{2})$
$\displaystyle\equiv\delta(t_{1}-t_{2})\frac{1}{\bar{y}(t_{1})}\left[\frac{d^{2}}{dt_{1}^{2}}+\omega^{2}+\frac{\bar{y}^{2}(t_{1})}{m}\right],$
(112) $\displaystyle O_{2}(t_{1},t_{2})$
$\displaystyle\equiv\delta(t_{1}-t_{2})\frac{1}{\bar{y}(t_{1})}\left[\frac{d^{2}}{dt_{1}^{2}}+\omega^{2}+3\frac{\bar{y}^{2}(t_{1})}{m}\right]=O_{1}(t_{1},t_{2})+2\,\delta(t_{1}-t_{2})\frac{\bar{y}(t_{1})}{m}.$

It follows that, given Eqs. (112), one can recast $H(t_{1},t_{2})$ in Eq.
(111) as

$H(t_{1},t_{2})=\int_{0}^{t}\mathrm{d}\tau\,O_{2}(t_{1},\tau)O_{1}^{-1}(\tau,t_{2})=\delta(t_{1}-t_{2})+2\frac{\bar{y}(t_{1})}{m}O_{1}^{-1}(t_{1},t_{2}),$
(113)

where the inverse operator satisfies the following relation

$\int_{0}^{t}\mathrm{d}\tau\,O_{1}(t_{1},\tau)\,O_{1}^{-1}(\tau,t_{2})=\delta(t_{1}-t_{2}).$
(114)

By plugging Eq. (110) in the inverse operator definition in Eq. (114) we
identify $O_{1}^{-1}(t_{1},t_{2})$ as

$O_{1}^{-1}(t_{1},t_{2})=-\frac{m}{2}\bar{y}_{1}(t_{2}|t_{1}).$ (115)

Hence, we have proved Eq. (111) to be true. It then follows that the path
integral evaluates to

$\displaystyle\int\mathcal{D}\varphi\exp\left\\{-\int_{0}^{t}\mathrm{d}t_{1}\int_{0}^{t}\mathrm{d}t_{2}\,\varphi(t_{1})H(t_{1},t_{2})\varphi(t_{2})\right\\}$
(116)
$\displaystyle=\left[\mathrm{det}\left(O_{2}\,O_{1}^{-1}\right)\right]^{-1/2}=\sqrt{\frac{\mathrm{det}\,O_{1}}{\mathrm{det}\,O_{2}}}=\sqrt{\frac{f(t)}{F(t)}},$

where in the last relation we exploited the fact that $f(t)$ and $F(t)$ are
respectively proportional to the determinant of $O_{1}$ and $O_{2}$ with the
same proportionality constant, according to Eq. (39) and the fact that the
determinant of a product of operators is given by the product of the
determinants of the individual operators.

## References

  * [1] T. Kinoshita, T. Wenger and D. S. Weiss,  _A quantum Newton’s cradle_ ,  Nature 440, 900 (2006),  10.1038/nature04693.
  * [2] C. N. Weiler, T. W. Neely, D. R. Scherer, A. S. Bradley, M. J. Davis and B. P. Anderson,  _Spontaneous vortices in the formation of Bose–Einstein condensates_ ,  Nature 455, 948 (2008),  10.1038/nature07334.
  * [3] M. Johanning, A. Varón and C. Wunderlich,  _Quantum simulations with cold trapped ions_ ,  J. Phys. B 42, 154009 (2009),  10.1088/0953-4075/42/15/154009.
  * [4] W. Bakr, J. Gillen, A. Peng, S. Fölling and M. Greiner,  _A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice_ ,  Nature 462, 74 (2009),  10.1038/nature08482.
  * [5] C. Weitenberg, M. Endres, J. F. Sherson, M. Cheneau, P. Schauß, T. Fukuhara, I. Bloch and S. Kuhr,  _Single-spin addressing in an atomic Mott insulator_ ,  Nature 471, 319 (2011),  10.1038/nature09827.
  * [6] R. Blatt and C. Roos,  _Quantum simulations with trapped ions_ ,  Nat. Phys. 8, 277 (2012),  10.1038/nphys2252.
  * [7] I. Bloch, J. Dalibard and W. Zwerger,  _Many-body physics with ultracold gases_ ,  Rev. Mod. Phys. 80, 885 (2008),  10.1103/RevModPhys.80.885.
  * [8] T. Langen, R. Geiger and J. Schmiedmayer,  _Ultracold atoms out of equilibrium_ ,  Annu. Rev. Condens. Matter Phys. 6, 201 (2015),  10.1146/annurev-conmatphys-031214-014548.
  * [9] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo and G. Ferrari,  _Spontaneous creation of Kibble–Zurek solitons in a Bose–Einstein condensate_ ,  Nat. Phys. 9, 656 (2013),  10.1038/nphys2734.
  * [10] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Foelling, L. Pollet and M. Greiner,  _Probing the superfluid–to–Mott insulator transition at the single-atom level_ ,  Science 329, 547 (2010),  10.1126/science.1192368.
  * [11] P. Calabrese, F. Essler and G. Mussardo,  _Introduction to ‘quantum integrability in out of equilibrium systems’_ ,  J. Stat. Mech.: Theory Exp. 2016, 064001 (2016),  10.1088/1742-5468/2016/06/064001.
  * [12] J. C. Bridgeman and C. T. Chubb,  _Hand-waving and interpretive dance: an introductory course on tensor networks_ ,  Journal of Physics A: Mathematical and Theoretical 50(22), 223001 (2017),  10.1088/1751-8121/aa6dc3.
  * [13] U. Schollwöck,  _The density-matrix renormalization group_ ,  Rev. Mod. Phys. 77, 259 (2005),  10.1103/RevModPhys.77.259.
  * [14] J. Schachenmayer, A. Pikovski and A. M. Rey,  _Many-body quantum spin dynamics with monte carlo trajectories on a discrete phase space_ ,  Phys. Rev. X 5, 011022 (2015),  10.1103/PhysRevX.5.011022.
  * [15] P. M. Hogan and J. T. Chalker,  _Path integrals, diffusion on SU (2) and the fully frustrated antiferromagnetic spin cluster_ ,  J. Phys. A: Math. Gen. 37, 11751 (2004),  10.1088/0305-4470/37/49/002.
  * [16] V. Galitski,  _Quantum-to-classical correspondence and Hubbard-Stratonovich dynamical systems: A Lie-algebraic approach_ ,  Phys. Rev. A 84, 012118 (2011),  10.1103/PhysRevA.84.012118.
  * [17] M. Ringel and V. Gritsev,  _Dynamical symmetry approach to path integrals of quantum spin systems_ ,  Phys. Rev. A 88, 062105 (2013),  10.1103/PhysRevA.88.062105.
  * [18] S. De Nicola, B. Doyon and M. J. Bhaseen,  _Stochastic approach to non-equilibrium quantum spin systems_ ,  J. Phys. A: Math. Theor. 52, 05LT02 (2019),  10.1088/1751-8121/aaf9be.
  * [19] C. M. Cheng and P. C. W. Fung,  _The evolution operator technique in solving the Schrödinger equation, and its application to disentangling exponential operators and solving the problem of a mass-varying harmonic oscillator_ ,  J. Phys. A: Math. Gen. 21, 4115 (1988),  10.1088/0305-4470/21/22/015.
  * [20] E. Kloeden, E. Platen and H. Schurz,  _Numerical solution of SDE through computer experiments_ ,  Springer Science & Business Media (2012).
  * [21] S. De Nicola,  _Disentanglement approach to quantum spin ground states: field theory and stochastic simulation_ ,  J. Stat. Mech.: Theory Exp. 2021, 013101 (2021),  10.1088/1742-5468/abc7c7.
  * [22] S. De Nicola, B. Doyon and M. J. Bhaseen,  _Non-equilibrium quantum spin dynamics from classical stochastic processes_ ,  J. Stat. Mech.: Theory Exp. 2020, 013106 (2020),  10.1088/1742-5468/ab6093.
  * [23] S. E. Begg, A. G. Green and M. J. Bhaseen,  _Fluctuations and non-Hermiticity in the stochastic approach to quantum spins_ ,  J. Phys. A 53, 50LT02 (2020),  10.1088/1751-8121/abbf87.
  * [24] R. H. Dicke,  _Coherence in spontaneous radiation processes_ ,  Phys. Rev. 93, 99 (1954),  10.1103/PhysRev.93.99.
  * [25] V. Giovannetti, S. Lloyd and L. Maccone,  _Capacity of nonlinear bosonic systems_ ,  Phys. Rev. A 70, 012307 (2004),  10.1103/PhysRevA.70.012307.
  * [26] S. Wald and M. Henkel,  _Lindblad dynamics of a quantum spherical spin_ ,  J. Phys. A: Math. Theor. 49, 125001 (2016),  10.1088/1751-8113/49/12/125001.
  * [27] M. Timpanaro, S. Wald, F. Semião and G. T. Landi,  _Dynamical chaotic phases and constrained quantum dynamics_ ,  Phys. Rev. A 100, 012117 (2019),  10.1103/PhysRevA.100.012117.
  * [28] S. Wald, S. V. Moreira and F. L. Semião,  _In- and out-of-equilibrium quantum metrology with mean-field quantum criticality_ ,  Phys. Rev. E 101, 052107 (2020),  10.1103/PhysRevE.101.052107.
  * [29] P. Romatschke,  _Quantum mechanical out-of-time-ordered-correlators for the anharmonic (quartic) oscillator_ ,  J. High Energy Phys. 2021, 30 (2021),  10.1007/JHEP01(2021)030.
  * [30] B. Simon and A. Dicke,  _Coupling constant analyticity for the anharmonic oscillator_ ,  Ann. Phys. 58, 76 (1970),  10.1016/0003-4916(70)90240-X.
  * [31] C. M. Bender and T. T. Wu,  _Anharmonic oscillator_ ,  Phys. Rev. 184, 1231 (1969),  10.1103/PhysRev.184.1231.
  * [32] L. Schulman,  _Techniques and applications of path integration_ ,  Courier Corporation (2012).
  * [33] M. F. Girard,  _On the semiclassical propagator for the anharmonic oscillator_ ,  J. Phys. A: Math. Gen. 25, 4199 (1992),  10.1088/0305-4470/25/15/027.
  * [34] G. Auberson and M. C. Peyranère,  _Quantum anharmonic oscillator in the Heisenberg picture and multiple scale techniques_ ,  Phys. Rev. A 65, 032120 (2002),  10.1103/PhysRevA.65.032120.
  * [35] F. R. Halpern,  _Quantum mechanics of the anharmonic oscillator_ ,  J. Math. Phys. 14, 219 (1973),  10.1063/1.1666299.
  * [36] C. M. Bender and T. T. Wu,  _Anharmonic oscillator. II. A study of perturbation theory in large order_ ,  Phys. Rev. D 7, 1620 (1973),  10.1103/PhysRevD.7.1620.
  * [37] M. M. Mizrahi,  _The semiclassical expansion of the anharmonic-oscillator propagator_ ,  J. Math. Phys. 20, 844 (1979),  10.1063/1.524157.
  * [38] I. G. Halliday and P. Suranyi,  _Anharmonic oscillator: A new approach_ ,  Phys. Rev. D 21, 1529 (1980),  10.1103/PhysRevD.21.1529.
  * [39] Y. Shalibo, Y. Rofe, I. Barth, L. Friedland, R. Bialczack, J. M. Martinis and N. Katz,  _Quantum and classical chirps in an anharmonic oscillator_ ,  Phys. Rev. Lett. 108, 037701 (2012),  10.1103/PhysRevLett.108.037701.
  * [40] Y. Shalibo, R. Resh, O. Fogel, D. Shwa, R. Bialczak, J. M. Martinis and N. Katz,  _Direct Wigner tomography of a superconducting anharmonic oscillator_ ,  Phys. Rev. Lett. 110, 100404 (2013),  10.1103/PhysRevLett.110.100404.
  * [41] K. Murch, R. Vijay, I. Barth, O. Naaman, J. Aumentado, L. Friedland and I. Siddiqi,  _Quantum fluctuations in the chirped pendulum_ ,  Nat. Phys. 7, 105 (2011),  10.1038/nphys1867.
  * [42] C. Lam,  _Stationary phase approximation of Feynman path integrals_ ,  Il Nuovo Cimento A 47, 451 (1967),  10.1007/BF02738742.
  * [43] D. Brizuela,  _Classical and quantum behavior of the harmonic and the quartic oscillators_ ,  Phys. Rev. D 90, 125018 (2014),  10.1103/PhysRevD.90.125018.
  * [44] Q. Dong, G.-H. Sun, M. A. Aoki, C.-Y. Chen and S.-H. Dong,  _Exact solutions of a quartic potential_ ,  Mod. Phys. Lett. A 34, 1950208 (2019),  10.1142/S0217732319502080.
  * [45] S. R.L.,  _On a method of calculating quantum distribution functions_ ,  Sov. Phys. Dokl. 2, 416 (1957).
  * [46] J. Hubbard,  _Calculation of partition functions_ ,  Phys. Rev. Lett. 3, 77 (1959),  10.1103/PhysRevLett.3.77.
  * [47] J. Wei and E. Norman,  _Lie algebraic solution of linear differential equations_ ,  J. Math. Phys. 4, 575 (1963),  10.1063/1.1703993.
  * [48] I. Kolokolov,  _Functional representation for the partition function of the quantum Heidenberg ferromagnet_ ,  Phys. Lett. A 114, 99 (1986),  10.1016/0375-9601(86)90488-3.
  * [49] M. Suzuki,  _Improved Trotter-like formula_ ,  Phys. Lett. A 180, 232 (1993),  10.1016/0375-9601(93)90701-Z.
  * [50] R. Shankar,  _Principles of quantum mechanics_ ,  Springer Science & Business Media (2012).
  * [51] F. Nicacio and F. L. Semião,  _Coupled harmonic systems as quantum buses in thermal environments_ ,  J. Phys. A: Math. Gen. 49, 375303 (2016),  10.1088/1751-8113/49/37/375303.
  * [52] F. Nicacio, A. Ferraro, A. Imparato, M. Paternostro and F. L. Semião,  _Thermal transport in out-of-equilibrium quantum harmonic chains_ ,  Phys. Rev. E 91, 042116 (2015),  10.1103/PhysRevE.91.042116.
  * [53] W. Case,  _Wigner functions and Weyl transforms for pedestrians_ ,  Am. J. Phys. 76, 937 (2008),  10.1119/1.2957889.
  * [54] M. Belloni, M. A. Doncheski and R. W. Robinett,  _Wigner quasi-probability distribution for the infinite square well: Energy eigenstates and time-dependent wave packets_ ,  Am. J. Phys. 72, 1183 (2004),  10.1119/1.1767100.
  * [55] J. Snygg,  _Use of operator wave functions to construct a refined correspondence principle via the quantum mechanics of Wigner and Moyal_ ,  Am. J. Phys. 48, 964 (1980),  10.1119/1.12197.
  * [56] J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics. 2-nd edition_ ,  Person New International edition (2014).
  * [57] J. Schwinger,  _Brownian motion of a quantum oscillator_ ,  J. Math. Phys. 2(3), 407 (1961).
  * [58] L. V. Keldysh _et al._ ,  _Diagram technique for nonequilibrium processes_ ,  Sov. Phys. JETP 20(4), 1018 (1965).
  * [59] A. Kamenev,  _Field theory of non-equilibrium systems_ ,  Cambridge University Press (2011).
  * [60] C. Moyer,  _Numerov extension of transparent boundary conditions for the Schrödinger equation in one dimension_ ,  Am. J. Phys. 72, 351 (2004),  10.1119/1.1619141.
  * [61] I. M. Gel’fand and A. Yaglom,  _Integration in functional spaces and its applications in quantum physics_ ,  J. Math. Phys. 1, 48 (1960),  10.1063/1.1703636.
  * [62] J. Zinn-Justin,  _Quantum field theory and critical phenomena_ , vol. 171,  Oxford University Press (2021).
  * [63] C. M. Bender, A. Felski, N. Hassanpour, S. P. Klevansky and A. Beygi,  _Analytic structure of eigenvalues of coupled quantum systems_ ,  Phys. Scr. 92, 015201 (2016),  10.1088/0031-8949/92/1/015201.
  * [64] A. Zee,  _Quantum field theory in a nutshell_ , vol. 7,  Princeton university press (2010).