# Phase-space stochastic quantum hydrodynamics for interacting Bose gases

S. A. Simmons  School of Mathematics and Physics, The University of
Queensland, Brisbane, Queensland 4072, Australia  J. C. Pillay  School of
Mathematics and Physics, The University of Queensland, Brisbane, Queensland
4072, Australia  K. V. Kheruntsyan  School of Mathematics and Physics, The
University of Queensland, Brisbane, Queensland 4072, Australia

###### Abstract

Hydrodynamic theories offer successful approaches that are capable of
simulating the otherwise difficult-to-compute dynamics of quantum many-body
systems. In this work we derive, within the positive-$P$ phase-space
formalism, a new stochastic hydrodynamic method for the description of
interacting Bose gases. It goes beyond existing hydrodynamic approaches, such
as superfluid hydrodynamics or generalized hydrodynamics, in its capacity to
simulate the full quantum dynamics of these systems: it possesses the ability
to compute non-equilibrium quantum correlations, even for short-wavelength
phenomena. Using this description, we derive a linearized stochastic
hydrodynamic scheme which is able to simulate such non-equilibrium situations
for longer times than the full positive-$P$ approach, at the expense of
approximating the treatment of quantum fluctuations, and show that this
linearized scheme can be directly connected with existing Bogoliubov
approaches. Furthermore, we go on to demonstrate the usefulness and advantages
of this formalism by exploring the correlations that arise in a quantum shock
wave scenario and comparing its predictions to other established quantum many-
body approaches.

##  I Introduction

Understanding and simulating the far-from-equilibrium behaviour of quantum
many-body systems is a notoriously challenging task. A particularly powerful
method capable of describing such non-equilibrium behaviour is the incredibly
successful hydrodynamic approach. Here our interest lies in the application of
hydrodynamic theories to the out-of-equilibrium dynamics of ultra-cold quantum
fluids, where there has been significant interest in recent years,
predominantly owing to the development of the theory of generalized
hydrodynamics (GHD) [1, 2, 3, 4, 5, 6], as well as the unforeseen success of
existing theories of classical hydrodynamics to situations where their a
priori applicability would have been questionable [7, 8]. Since hydrodynamic
theories are often postulated phenomenologically, we wish to gain a deeper
understanding of their emergence from a microscopic description of the system,
and to provide such a microscopic theory in terms of hydrodynamic variables
that is capable of going beyond existing hydrodynamic approaches to include
quantum correlations and short-wavelength phenomena.

Conventional or classical hydrodynamics (CHD) itself is postulated using the
local density approximation (LDA) where the fluid is treated by breaking it up
into small uniform slices which are then assumed to be in thermal equilibrium
in the local moving frame. Nevertheless, the classical hydrodynamic equations
themselves are in fact nothing more than equations which express the
conservation of mass, momentum, and energy (or entropy) of a system [9, 10].
Hence, their applicability goes beyond that of a classical ideal gas and they
can be used to simulate the large scale dynamics of both weakly and strongly
interacting quantum fluids by employing the appropriate thermodynamic equation
of state [7, 8, 11, 12, 13, 14]. Yet such an approach is limited in its
applicability in that it usually requires that a system possess fast
themalization times, and due to the LDA it is restricted to the description of
long-wavelength excitations.

To go beyond the classical hydrodynamic description, one can employ superfluid
hydrodynamics, often referred to as quantum hydrodynamics, which is a well-
known and fundamental tool used to describe ultra-cold interacting Bose gases
[15, 16]. Through the use of Madelung’s transformation [17], it is usually
derived from the mean-field Gross-Pitaevskii equation (which describes weakly
interacting Bose-Einstein condensates of ultra-cold atomic gases at zero-
temperature) and goes beyond the classical hydrodynamic approach to include
the so-called quantum pressure term, which allows for the accurate description
of short-wavelength phenomena. It is also possible to extend such a
description to finite temperatures by employing a two-fluid model where the
single fluid is separated into two components; the condensed atoms (or the
superfluid component) which are in the ground-state of the system, and the
non-condensed atoms that comprise the thermal cloud (or the normal component)
[18]. Here, the condensed atoms are treated according to superfluid
hydrodynamics and are then coupled to the non-condensed atoms that are treated
classically. Such an approach can be directly connected to Landau’s celebrated
two-fluid hydrodynamic formalism, which itself was postulated
phenomenologically to describe the strongly interacting system of superfluid
Helium [19, 20, 15]. Despite the powerful utility of the superfluid
hydrodynamic approach, it relies fundamentally on the assumption of a mean-
field or existence of long-range order, which neglects any quantum
fluctuations and correlations, and as such, it is not a truly quantum
description.

The theory of GHD exceeds that of classical hydrodynamics in a different
manner to superfluid hydrodynamics, and it proceeds through the identification
of infinitely many conserved quantities that exist for integrable systems.
Including these higher order conserved quantities in its description allows
GHD to simulate scenarios where CHD would otherwise suffer from the notorious
gradient catastrophe (or derivative discontinuity) problem, and provides a
more accurate description at finite-temperatures [3]. Standard GHD is capable
of simulating ultra-cold Bose gases at any interaction strength, and since its
discovery it has been extended to address (among other scenarios) near-
integrable systems [21, 6] and to include quantum fluctuations [22, 23]. Yet,
to date, these descriptions remain restricted to long-wavelength excitations
and are unable to capture quantum interference and other important short-
wavelength phenomena such as the ones present in dispersive quantum shock
waves [24, 25, 26, 27, 28, 29, 30, 14, 31, 13, 32].

It is an interesting question to ask then, whether an exact quantum
description of interacting Bose gases (which is capable of capturing both
quantum correlations and short-wavelength phenomena) is possible in terms of
hydrodynamic variables and what it would look like. In this work we show that
it is possible, in the weakly-interacting regime, and that such a theory can
be formulated in a phase-space representation; we refer to it as stochastic
quantum hydrodynamics (SQHD). By employing the appropriate assumptions (of
small quantum fluctuations), we additionally derive the new approximate
hydrodynamic equations of linearised stochastic quantum hydrodynamics (LSQHD)
and provide connections between these and other existing approaches capable of
treating interacting Bose gases. Furthermore, we apply our novel approach to
the dynamics of dispersive quantum shock waves [32] where we calculate for the
first time the quantum correlations arising in this situation and benchmark
our results against alternate theoretical approaches.

##  II Exact quantum hydrodynamics: stochastic formulation

We consider an effective field theory Hamiltonian for a trapped Bose gas with
repulsive delta-function interactions in second-quantized form

$\hat{H}=\int
d{\mathbf{r}}\left\\{\frac{\hbar^{2}}{2m}|\mathbf{\nabla}\hat{\Psi}|^{2}+V_{\mathrm{ext}}({\mathbf{r}},t)\hat{\Psi}^{{\dagger}}\hat{\Psi}+\frac{g}{2}\hat{\Psi}^{{\dagger}}\hat{\Psi}^{{\dagger}}\hat{\Psi}\hat{\Psi}\right\\},$
(1)

where $\hat{\Psi}({\mathbf{r}},t)$ is the bosonic field operator with the
usual equal-time commutation relation
$[\hat{\Psi}({\mathbf{r}},t),\hat{\Psi}^{\dagger}({\mathbf{r}}^{\prime},t)]=\delta^{(3)}({\mathbf{r}}-{\mathbf{r}}^{\prime})$
and $m$ is the atomic mass. The first term in Eq. (1) is the kinetic energy,
and the last term describes the $s$-wave scattering interactions between the
particles, with the coupling strength $g=4\pi\hbar^{2}a/m>0$ in three
dimensions (3D), where $a$ is the $s$-wave scattering length. The external
trapping potential is given by $V_{\mathrm{ext}}({\mathbf{r}},t)$, which can
in general be time-dependent.

In an elongated (cigar-shaped) trap geometry, this model can easily be reduced
to the trapped gas Lieb-Liniger model [33] in one-dimension (1D) if $a$ is
much smaller than the amplitude of transverse zero point oscillations (or the
harmonic oscillator length) $l_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$, where
$\omega_{\perp}$ is the frequency of transverse (radial) confinement assumed
harmonic and symmetric. In this case, the interaction constant is given by
$g\rightarrow g_{\mathrm{1D}}\simeq 2\hbar\omega_{\perp}a$ away from
confinement induced resonances [34]. Similarly, in a pancake geometry, the
model can be reduced to two dimensions (2D), with $g\rightarrow
g_{\mathrm{2D}}=2\sqrt{2\pi}\hbar^{2}a/ml_{0}$ [35], where $l_{0}$ is the
harmonic oscillator length in the axial dimension.

We further note that the use of the effective delta-function interaction
potential $g\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$ between atom pairs,
leading to Eq. (1), must assume a UV momentum cutoff $k^{\max}$ to prevent the
theory from UV divergencies [36]. In the numerical simulations below, such a
momentum cutoff is imposed explicitly via the finite computational lattice. If
the lattice spacings ($\Delta x$, $\Delta y$, $\Delta z$) in each spatial
dimensions are chosen to be much larger than the $s$-wave scattering length
$a$, then the respective momentum cutoffs satisfy $k_{x,y,z}^{\max}\ll 1/a$.
In this case, the coupling constant $g$ can be approximated by the expression
given above, without the need for explicit renormalization [37].

To model the dynamics of quantum fields governed by the above Hamiltonian, we
use the positive-$P$ phase-space approach [38, 39, 40, 41, 42, 43, 44, 45]. In
this approach, the dynamics of the bosonic field operators
$\hat{\Psi}(\mathbf{r},t)$ and $\hat{\Psi}^{\dagger}(\mathbf{r},t)$ are
equivalent to the evolution of two independent complex $c$-fields,
$\Psi(\mathbf{r},t)$ and $\widetilde{\Psi}(\mathbf{r},t)$, satisfying the Îto
stochastic differential equations (SDEs) [43, 44, 45]:

$\displaystyle\frac{\partial\Psi(\mathbf{r},t)}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{i}{\hbar}\left(\frac{\hbar^{2}}{2m}\laplacian-
V_{\mathrm{ext}}(\mathbf{r},t)-g\widetilde{\Psi}\Psi\right)\Psi$ (2)
$\displaystyle+\sqrt{-i\frac{g}{\hbar}}\Psi\,\xi_{1}(\mathbf{r},t),$
$\displaystyle\frac{\partial\widetilde{\Psi}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\left(\frac{\hbar^{2}}{2m}\laplacian-
V_{\mathrm{ext}}(\mathbf{r},t)-g\widetilde{\Psi}\Psi\right)\widetilde{\Psi}$
(3)
$\displaystyle+\sqrt{i\frac{g}{\hbar}}\widetilde{\Psi}\,\xi_{2}(\mathbf{r},t),$

where $\laplacian=\bm{\nabla}\cdot\bm{\nabla}$ is the Laplacian operator,
$\xi_{j}(\mathbf{r},t)$ ($j=1,2$) are real independent Gaussian noise sources
with zero mean and the following nonzero correlations
$\langle\xi_{j}(\mathbf{r},t)\xi_{k}(\mathbf{r}^{\prime},t^{\prime})\rangle=\delta_{jk}\delta^{(3)}(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$.

The above SDEs are derived from an equivalent Fokker-Planck equation for the
positive-$P$ quasiprobability distribution function in phase space, which
itself is derived from the master equation for the quantum density operator
[38, 39]. The derivation relies on the assumption that the positive-$P$
distribution function is sufficiently bounded so as to lead to vanishing
boundary terms in the integration-by-parts step of the derivation of the
Fokker-Planck equation. For vanishing boundary terms, the positive-$P$
approach and the resulting SDEs are equivalent to simulating the _exact_
quantum dynamics of the system. Whether this assumption is actually satisfied
or not can in practice be checked _a posteriori, i.e.,_ after carrying out the
numerical simulations of the SDEs. Vanishing boundary terms lead to stable
evolution, with sampling errors that can be controlled and reduced by
increasing the number of stochastic realizations required to calculate the
stochastic averages. Non-vanishing boundary terms, on the other hand, lead to
the well-known occurrence of ‘spiking’ trajectories and a development of
power-law tails in the quasiprobability distribution function [46]. This in
turn leads to rapid and uncontrollable growth of sampling errors as the
simulation progresses in time [42], implying that the simulation results can
no longer be used past the ‘spiking’ time.

The complex stochastic fields $\Psi(\mathbf{r},t)$ and
$\widetilde{\Psi}(\mathbf{r},t)$ are independent of each other
[$\widetilde{\Psi}(\mathbf{r},t)\neq\Psi^{\ast}(\mathbf{r},t)$] except in the
mean,
$\langle\widetilde{\Psi}(\mathbf{r},t)\rangle=\langle\Psi^{\ast}(\mathbf{\mathbf{r}},t)\rangle$,
where the brackets $\langle\ldots\rangle$ refer to stochastic averages with
respect to the positive-$P$ distribution function. In numerical realizations,
this is represented by an ensemble average over a large number of stochastic
realizations (trajectories). Observables described by quantum mechanical
ensemble averages over normally-ordered products of the field operators have
an exact correspondence with stochastic averages over products of the
$c$-fields $\Psi(\mathbf{r},t)$ and $\widetilde{\Psi}(\mathbf{r},t)$, such
that

$\langle[\hat{\Psi}^{\dagger}(\mathbf{r},t)]^{m}[\hat{\Psi}(\mathbf{r}^{\prime},t)]^{n}\rangle=\langle[\widetilde{\Psi}(\mathbf{r},t)]^{m}[\Psi(\mathbf{r}^{\prime},t)]^{n}\rangle.$
(4)

Converting the above SDEs from Îto to Stratonovich form, which has the
advantage of obeying the standard rules of calculus, leads only to a
physically irrelevant constant phase shift term which can be ignored. Thus,
the Stratonovich form of the SDEs is the same as the above Îto form. Once the
Stratonovich form of the SDEs is established, we can transform to a new pair
of complex variables – ‘density’ and ‘phase’ – defined via

$\displaystyle\rho(\mathbf{r},t)$
$\displaystyle=\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r},t),$ (5)
$\displaystyle S(\mathbf{r},t)$
$\displaystyle=\frac{1}{2i}\ln\frac{\Psi(\mathbf{r},t)}{\widetilde{\Psi}(\mathbf{r},t)}=\frac{1}{2i}\left(\ln\Psi(\mathbf{r},t)-\ln\widetilde{\Psi}(\mathbf{r},t)\right).$
(6)

These transformations correspond to

$\displaystyle\Psi(\mathbf{r},t)$
$\displaystyle=\sqrt{\rho(\mathbf{r},t)}e^{iS(\mathbf{r},t)},$ (7)
$\displaystyle\widetilde{\Psi}(\mathbf{r},t)$
$\displaystyle=\sqrt{\rho(\mathbf{r},t)}e^{-iS(\mathbf{r},t)}.$ (8)

The new SDEs for the complex density and phase variables can be written as

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\frac{\hbar}{m}\bm{\nabla}\cdot(\rho\bm{\nabla}S)+\sqrt{-i\frac{g}{\hbar}}\rho\left(\xi_{1}+i\xi_{2}\right),$
(9) $\displaystyle\frac{\partial S}{\partial t}$
$\displaystyle=-\frac{\hbar}{2m}(\bm{\nabla}S)^{2}-\frac{V_{\mathrm{ext}}}{\hbar}-\frac{g}{\hbar}\rho+\frac{\hbar}{2m\sqrt{\rho}}\laplacian\sqrt{\rho}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad-\frac{1}{2}\sqrt{i\frac{g}{\hbar}}\left(\xi_{1}-i\xi_{2}\right),$
(10)

where the term involving $\laplacian\sqrt{\rho}$ can be alternatively
rewritten as

$\frac{1}{\sqrt{\rho}}\laplacian\sqrt{\rho}=\frac{1}{2\rho}\laplacian\rho-\frac{1}{4\rho^{2}}(\bm{\nabla}\rho)^{2}=\frac{1}{2}\laplacian\ln\rho+\frac{1}{4}(\bm{\nabla}\ln\rho)^{2}.$
(11)

Note that the noise term in the equation for the density field $\rho$ is
_multiplicative_ (_i.e_., it depends on the stochastic variable itself),
whereas the noise term in the equation for the phase field $S$ is _additive_
and therefore it corresponds to a random stochastic process that has Gaussian
statistics. The Gaussian statistics for the phase variable implies that, for
example, stochastic averages of the form $\left\langle e^{iS}\right\rangle$
can be expressed as

$\langle e^{iS}\rangle=e^{-\frac{1}{2}\langle S^{2}\rangle}.$ (12)

This step, however, is not required from a numerical point of view as the
stochastic average value of $\langle e^{iS}\rangle$ can be evaluated
explicitly just as easily as $\langle S^{2}\rangle$.

The final step in deriving the stochastic hydrodynamic equations involves the
introduction of a (complex) stochastic ‘velocity’ vector field

$\bm{v}(\mathbf{r},t)=\frac{\hbar}{m}\bm{\nabla}S(\mathbf{r},t),$ (13)

which is analogous to the standard velocity vector in conventional
hydrodynamics: after taking the stochastic average, the expectation value of
the stochastic velocity vector becomes real and acquires the same physical
meaning as in conventional hydrodynamics. This leads us to the following SQHD
equations,

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot(\rho\bm{v})+\sqrt{-i\frac{g}{\hbar}}\rho\left(\xi_{1}+i\xi_{2}\right),$
(14) $\displaystyle\frac{\partial\bm{v}}{\partial t}$
$\displaystyle=-\bm{\nabla}\left[\frac{1}{2}\bm{v}^{2}+\frac{V_{\mathrm{ext}}}{m}+\frac{g}{m}\rho-\frac{\hbar^{2}}{2m^{2}\sqrt{\rho}}\laplacian\sqrt{\rho}\right.$
$\displaystyle\left.\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{1}{2m}\sqrt{i\hbar
g}\left(\xi_{1}-i\xi_{2}\right)\right],$ (15)

which are the first key result of this work.

As we will show below, the utility of these equations is three-fold; (_i_)
they provide a theoretically exact description of interacting Bose gases in
terms of hydrodynamic variables; (_ii_) they can be used to derive new
approximate hydrodynamic schemes, where connections to alternative theories
can be directly shown; and (_iii_) such hydrodynamic theories can be
advantageous over alternatives in a number of different physical scenarios.

##  III Approximate schemes in stochastic quantum hydrodynamics

###  III.1 Suppressed density and phase fluctuations: hydrodynamics of pure
condensates

In this subsection we would like to apply the SQHD approach to Bose-Einstein
condensates characterized by suppressed fluctuations of both the density and
phase such that they can be ignored altogether. The goal is to derive the
mean-field superfluid hydrodynamic equations directly from the stochastic
hydrodynamic equations, providing an alternate derivation to the usual
spontaneously broken symmetry approach (which leads to the Gross-Pitaevskii
equation for the order parameter, and equations for the mean-field density and
phase/velocity through the use of Madelung’s transformation) [15, 16, 47].

By taking the expectation values of both sides of Eqs. (14) and (15), and
taking into account that the noise terms $\xi_{i}$ are independent of the
stochastic variables and have zero mean, we see that these equations will
reduce to hydrodynamic equations for the mean density and the mean velocity,

$\displaystyle\frac{\partial\langle\rho\rangle}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot(\langle\rho\rangle\langle\bm{v}\rangle),$
(16) $\displaystyle\frac{\partial\langle\bm{v}\rangle}{\partial t}$
$\displaystyle=-\bm{\nabla}\left[\frac{1}{2}\langle\bm{v}\rangle^{2}+\frac{V_{\mathrm{ext}}}{m}+\frac{g}{m}\langle\rho\rangle\right.$
$\displaystyle\qquad\qquad\qquad\left.-\frac{\hbar^{2}}{2m^{2}\sqrt{\langle\rho\rangle}}\laplacian\sqrt{\langle\rho\rangle}\right],$
(17)

if we assume that:

(_i_) the average of products of stochastic velocities (phases) factorize into
products of average velocities (phases), so that
$\langle\bm{v}^{2}\rangle=\langle\bm{v}\rangle^{2}$. Such a condition implies
that there are no longer any local velocity fluctuations (and therefore any
local phase fluctuations) in the system since
$\sigma_{\bm{v}}=\sqrt{\langle\bm{v}^{2}\rangle-\langle\bm{v}\rangle^{2}}=0$;

(_ii_) the stochastic density and velocity (phase) variables are uncorrelated,
leading to the replacement
$\langle\rho\bm{v}\rangle=\langle\rho\rangle\langle\bm{v}\rangle$. This
follows from the condition above, since there are no velocity (phase)
fluctuations;

(_iii_) the average of products/ratios of the stochastic densities factorize
into respective products/ratios of average densities, so that the term
$\langle\frac{1}{\sqrt{\rho}}\laplacian\sqrt{\rho}\rangle$ can be replaced by
$\frac{1}{\sqrt{\langle\rho\rangle}}\laplacian\sqrt{\langle\rho\rangle}$.

Assumptions (_i_)–(_iii_) imply that all higher order moments of local
operators factorize, leaving us with a coherent state, or equivalently, a
mean-field description [48]. It is satisfying to note then, that this
alternate derivation has led us back to equations of the same form as the
usual mean-field superfluid hydrodynamic equations [15, 16, 47],

$\displaystyle\frac{\partial\rho_{0}}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot(\rho_{0}\bm{v}_{0}),$ (18)
$\displaystyle\frac{\partial\bm{v}_{0}}{\partial t}$
$\displaystyle=-\bm{\nabla}\left(\frac{1}{2}\bm{v}_{0}^{2}+\frac{V_{\mathrm{ext}}}{m}+\frac{g}{m}\rho_{0}-\frac{\hbar^{2}}{2m^{2}\sqrt{\rho_{0}}}\laplacian\sqrt{\rho_{0}}\right),$
(19) $\displaystyle\frac{\partial S_{0}}{\partial t}$
$\displaystyle=-\frac{\hbar}{2m}(\bm{\nabla}S_{0})^{2}-\frac{V_{\mathrm{ext}}}{\hbar}-\frac{g}{\hbar}\rho_{0}+\frac{\hbar}{2m\sqrt{\rho_{0}}}\laplacian\sqrt{\rho_{0}},$
(20)

where one can identify that,

$\langle\rho\rangle=\rho_{0},\quad\langle\bm{v}\rangle=\bm{v}_{0},\quad\langle
S\rangle=S_{0}.$ (21)

Equations (18)-(20) of superfluid hydrodynamics can be easily derived in the
spontaneously broken symmetry approach by applying Madelung’s transformation
$\Psi_{0}=\sqrt{\rho_{0}}e^{iS_{0}}$ to the order parameter
$\Psi_{0}=\langle\hat{\Psi}\rangle=\langle\Psi\rangle$, with
$\Psi_{0}^{*}=\langle\hat{\Psi}^{\dagger}\rangle=\langle\tilde{\Psi}\rangle$,
along with the definition $\bm{v}_{0}=\frac{\hbar}{m}\bm{\nabla}S_{0}$, where
the order parameter itself evolves according to the mean-field Gross-
Pitaevskii equation (GPE),

$\displaystyle\frac{\partial\Psi_{0}}{\partial t}$
$\displaystyle=\frac{i}{\hbar}\left[\frac{\hbar^{2}}{2m}\laplacian-
V_{\text{ext}}-g|\Psi_{0}|^{2}\right]\Psi_{0}.$ (22)

When dealing with the superfluid hydrodynamic equations (18)-(20), the quantum
pressure term $\frac{1}{\sqrt{\rho_{0}}}\laplacian\sqrt{\rho_{0}}$ is often
neglected in the literature since it is responsible for dynamics on short
length scales $\sim l_{\text{h}}=\hbar/\sqrt{mg\rho_{0}}$ (the healing length)
and can be ignored when considering the bulk dynamics of a system at large
scales [49, 26]. Removing this term results in the classical Euler
hydrodynamic equations for the density and velocity,

$\displaystyle\frac{\partial\rho_{0}}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot(\rho_{0}\bm{v}_{0}),$ (23)
$\displaystyle\frac{\partial\bm{v}_{0}}{\partial t}$
$\displaystyle=-\bm{\nabla}\left(\frac{1}{2}\bm{v}_{0}^{2}+\frac{V_{\mathrm{ext}}}{m}\right)-\frac{1}{m\rho_{0}}\bm{\nabla}P.$
(24)

Such equations can be postulated phenomenologically by considering the
dynamics of small uniform slices of the Bose gas and employing the local
density approximation. Here we have written the equations in standard form
where $P=g\rho_{0}^{2}/2$ is the pressure (or equation of state) of a weakly
interacting Bose gas at zero temperature [15], related to the chemical
potential $\mu=g\rho_{0}$ via the Gibbs-Duhem relation $dP=\rho_{0}d\mu$ [16].

###  III.2 Linearized stochastic quantum hydrodynamics

The next typical situation we would like to address via the stochastic quantum
hydrodynamic approach is the description of the dynamics of quantum fluids
that have small density and phase fluctuations. We would like to simplify the
stochastic equations (14) and (15) in such a way that the _multiplicative_
noise terms (which are known to be the source of uncontrollable sampling
errors past the ‘spiking’ time [46, 42]) can be approximated by _additive_
noise. This can be done in a linearized treatment of fluctuations around their
deterministic (mean-field) value by assuming that the fluctuations are small
and truncating the relevant expansion terms in the equations of motion.

Thus, we proceed by decomposing the density, velocity, and phase variables as
a sum of two components: the deterministic (mean-field) components and the
fluctuating stochastic components,

$\displaystyle\rho$
$\displaystyle=\left\langle\rho\right\rangle+\delta\rho=\rho_{0}+\delta\rho,$
(25) $\displaystyle\bm{v}$
$\displaystyle=\left\langle\bm{v}\right\rangle+\delta\bm{v}=\bm{v}_{0}+\delta\bm{v},$
(26) $\displaystyle S$ $\displaystyle=\left\langle S\right\rangle+\delta
S=S_{0}+\delta S.$ (27)

The deterministic components are described by the superfluid hydrodynamic
equations, Eqs. (18)–(20), whereas the complex fluctuating components are
described to second order by the following SDEs:

$\displaystyle\frac{\partial\delta\rho}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot\left[\rho_{0}\delta\bm{v}+(\bm{v}_{0}+\delta\bm{v})\delta\rho\right]$
$\displaystyle\qquad\qquad\qquad+\sqrt{-i\frac{g}{\hbar}}(\rho_{0}+\delta\rho)\,(\xi_{1}+i\xi_{2}),$
(28) $\displaystyle\frac{\partial\delta\bm{v}}{\partial t}$
$\displaystyle=-\bm{\nabla}\left[(\bm{v}_{0}\cdot\delta\bm{v})+\frac{1}{2}\delta\bm{v}^{2}+\frac{g}{m}\delta\rho\right.$
$\displaystyle\quad-\frac{\hbar^{2}}{2m^{2}}\left\\{-\frac{1}{2}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\sqrt{\rho_{0}}+\frac{1}{2}\frac{1}{\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\sqrt{\rho_{0}}}\right)\right.$
$\displaystyle\quad+\frac{3}{8}\left(\frac{\delta\rho}{\rho_{0}^{5/4}}\right)^{2}\laplacian\sqrt{\rho_{0}}-\frac{1}{4}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\left(\frac{\delta\rho}{\rho_{0}^{3/2}}\right)$
$\displaystyle\quad\left.\left.-\frac{1}{8}\frac{1}{\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\rho_{0}^{3/4}}\right)^{2}+\mathcal{O}(\delta\rho^{3})\right\\}\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad\quad+\frac{1}{2m}\sqrt{i\hbar
g}\,(\xi_{1}-i\xi_{2})\right],$ (29)
$\displaystyle\frac{\partial\mathcal{\delta S}}{\partial t}$
$\displaystyle=-\frac{\hbar}{m}(\bm{\nabla}S_{0}\cdot\bm{\nabla}\delta
S)-\frac{\hbar}{2m}(\bm{\nabla}\delta S)^{2}-\frac{g}{\hbar}\delta\rho$
$\displaystyle\quad+\frac{\hbar}{2m}\left\\{-\frac{1}{2}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\sqrt{\rho_{0}}+\frac{1}{2}\frac{1}{\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\sqrt{\rho_{0}}}\right)\right.$
$\displaystyle\quad+\frac{3}{8}\left(\frac{\delta\rho}{\rho_{0}^{5/4}}\right)^{2}\laplacian\sqrt{\rho_{0}}-\frac{1}{4}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\left(\frac{\delta\rho}{\rho_{0}^{3/2}}\right)$
$\displaystyle\quad-\left.\frac{1}{8}\frac{1}{\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\rho_{0}^{3/4}}\right)^{2}+\mathcal{O}(\delta\rho^{3})\right\\}$
$\displaystyle\qquad\qquad\qquad\qquad\quad-\frac{1}{2}\sqrt{i\frac{g}{\hbar}}\,(\xi_{1}-i\xi_{2}),$
(30)

where the relationship between $\delta S$ and $\delta\bm{v}$ is given by
$\delta\bm{v}=\frac{\hbar}{m}\bm{\nabla}\delta S$, similarly to the usual
$\bm{v}_{0}=\frac{\hbar}{m}\bm{\nabla}S_{0}$. We keep the phase equation here
since it can be a more convenient variable to simulate for the computation of
certain observables (see Sec. IV.1), whereas the velocity equation provides
the more natural set of SDEs from a hydrodynamic perspective.

The form of these equations above is particularly helpful in highlighting the
role of the stochastic many-body quantum pressure term
$\frac{1}{\sqrt{\rho}}\laplacian\sqrt{\rho}$ in the full SQHD equations.
Whilst formally similar to the mean-field quantum pressure
$\frac{1}{\sqrt{\rho_{0}}}\laplacian\sqrt{\rho_{0}}$, the many-body quantum
pressure term contains much more information, incorporating high order density
fluctuations relevant to the many-body nature of interacting Bose gases. This
additional information contained within the many-body quantum pressure term
can be seen in equations (29) and (30) where the terms inside
$\left\\{\cdots\right\\}$ describe the remaining higher order density
fluctuations after truncating the infinite series expansion down to second
order. Whilst in this work we are going to further restrict ourselves to first
order in fluctuations, in principle one can incorporate these second and even
higher order expansion terms in the many-body quantum pressure as an avenue
for future developments (see Sec. VI below).

We also note here that, at second order in the fluctuating components, the
many-body quantum pressure is the only term that has been truncated. All other
terms do not contribute anything higher than second order in fluctuations,
i.e. no terms involving the phase/velocity field have been neglected in
equations (28)-(30). As a foreshadow of Section III.3 below, we mention that
keeping these second order terms corresponds to a description which goes
beyond the Bogoliubov treatment of interacting Bose gases. Again, we leave
this open as an interesting direction to pursue in future work.

Proceeding further and truncating down to linear order, we arrive at

$\displaystyle\frac{\partial\delta\rho}{\partial t}$
$\displaystyle=-\bm{\nabla}\cdot\left[\rho_{0}\delta\bm{v}+\bm{v}_{0}\delta\rho\right]+\sqrt{-i\frac{g}{\hbar}}\rho_{0}\,(\xi_{1}+i\xi_{2}),$
(31) $\displaystyle\frac{\partial\delta\bm{v}}{\partial t}$
$\displaystyle=-\bm{\nabla}\left[(\bm{v}_{0}\cdot\delta\bm{v})+\frac{g}{m}\delta\rho-\frac{\hbar^{2}}{4m^{2}\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\sqrt{\rho_{0}}}\right)\right.$
$\displaystyle\quad+\left.\frac{\hbar^{2}}{4m^{2}}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\sqrt{\rho_{0}}+\frac{1}{2m}\sqrt{i\hbar
g}\,(\xi_{1}-i\xi_{2})\right],$ (32)
$\displaystyle\frac{\partial\mathcal{\delta S}}{\partial t}$
$\displaystyle=-\frac{\hbar}{m}(\bm{\nabla}S_{0}\cdot\bm{\nabla}\delta
S)-\frac{g}{\hbar}\delta\rho+\frac{\hbar}{4m\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\sqrt{\rho_{0}}}\right)$
$\displaystyle\quad-\frac{\hbar}{4m}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\sqrt{\rho_{0}}-\frac{1}{2}\sqrt{i\frac{g}{\hbar}}\,(\xi_{1}-i\xi_{2}).$
(33)

As we see, all fluctuating noise terms are now additive as the contribution of
terms involving products of stochastic components has been neglected on the
grounds of being higher-than-linear order in small fluctuations. By small here
we mean that
$\langle\delta\rho(\mathbf{r},t)\rangle\ll\rho_{0}(\mathbf{r},t)$,
$\langle\delta\bm{v}(\mathbf{r},t)\rangle\ll\bm{v}_{0}(\mathbf{r},t)$, and
$\langle\delta S(\mathbf{r},t)\rangle\ll S_{0}(\mathbf{r},t)$. In the case
when $\langle\delta\rho(\mathbf{r},t)\rangle=0$,
$\langle\delta\bm{v}(\mathbf{r},t)\rangle=0$, $\langle\delta
S(\mathbf{r},t)\rangle=0$, then one would consider
$\langle\delta\rho^{2}(\mathbf{r},t)\rangle\ll\rho_{0}^{2}(\mathbf{r},t)$,
$\langle\delta\bm{v}^{2}(\mathbf{r},t)\rangle\ll\bm{v}_{0}^{2}(\mathbf{r},t)$,
and $\langle\delta S^{2}(\mathbf{r},t)\rangle\ll S_{0}^{2}(\mathbf{r},t)$.

The above LSQHD equations (31)–(33), which have to be solved numerically in
conjunction with the mean-field equations (18)-(20), are the second key result
of this work. Although they are approximate, their advantage over the full
SQHD equations lies in the fact that they remain stable and can be simulated
for significantly longer times (and remain valid provided the fluctuations do
not grow too large) due to the _additive_ nature of the noise terms, as
opposed to the _multiplicative_ noise appearing in the full SQHD equations.

###  III.3 Comparison with Bogoliubov approaches

To elucidate the utility of the LSQHD equations (31)–(33), and before turning
to numerical examples, we now explore their comparison with other known
numerical methods in the literature. In particular, since these equations are
derived under a linearized scheme of quantum fluctuations, we focus on
connections to Bogoliubov type theories.

Specifically, under the transformations

$\displaystyle\delta\psi$
$\displaystyle=\Psi_{0}\left(\frac{\delta\rho}{2\rho_{0}}+i\delta
S\right)=\sqrt{\rho_{0}}e^{iS_{0}}\left(\frac{\delta\rho}{2\rho_{0}}+i\delta
S\right),$ (34) $\displaystyle\delta\tilde{\psi}$
$\displaystyle=\Psi_{0}^{*}\left(\frac{\delta\rho}{2\rho_{0}}-i\delta
S\right)=\sqrt{\rho_{0}}e^{-iS_{0}}\left(\frac{\delta\rho}{2\rho_{0}}-i\delta
S\right),$ (35)

which correspond to
$\delta\rho=\Psi_{0}^{*}\delta\psi+\Psi_{0}\delta\tilde{\psi}$ and $\delta
S=\frac{1}{2i}\left[\frac{\delta\psi}{\Psi_{0}}-\frac{\delta\tilde{\psi}}{\Psi_{0}^{*}}\right]$,
the LSQHD equations become

$\displaystyle\frac{\partial\delta\psi(\mathbf{r},t)}{\partial t}$
$\displaystyle=\frac{i}{\hbar}\left[\frac{\hbar^{2}}{2m}\laplacian-
V_{\text{ext}}(\mathbf{r},t)-2g|\Psi_{0}|^{2}\right]\delta\psi(\mathbf{r},t)$
$\displaystyle\quad-\frac{i}{\hbar}g\Psi_{0}^{2}\delta\tilde{\psi}(\mathbf{r},t)+\sqrt{-i\frac{g}{\hbar}}\Psi_{0}\,\xi_{1}(\mathbf{r},t),$
(36) $\displaystyle\frac{\partial\delta\tilde{\psi}(\mathbf{r},t)}{\partial
t}$ $\displaystyle=-\frac{i}{\hbar}\left[\frac{\hbar^{2}}{2m}\laplacian-
V_{\text{ext}}(\mathbf{r},t)-2g|\Psi_{0}|^{2}\right]\delta\tilde{\psi}(\mathbf{r},t)$
$\displaystyle\quad+\frac{i}{\hbar}g(\Psi_{0}^{*})^{2}\delta\psi(\mathbf{r},t)+\sqrt{i\frac{g}{\hbar}}\Psi_{0}^{*}\,\xi_{2}(\mathbf{r},t),$
(37)

where $\Psi_{0}(\mathbf{r},t)$ is the mean-field order parameter that evolves
according to the GPE (22).

Equations (36) and (37) are equivalent to the stochastic Bogoliubov equations
first introduced in Ref. [45] (see also [50]). Such equations have been used
to successfully describe several experimental results on Bose-Einstein
condensate collisions [45, 51, 52, 53, 54, 55, 56], and a more recent
extension of this approach was used to characterize the quantum depletion of
an expanding condensate [57]. We point out that the stochastic Bogoliubov
equations themselves are equivalent to the usual Bogoliubov approach [15, 16,
58, 59] in which the Bose field operator is expanded as
$\hat{\Psi}=\Psi_{0}+\delta\hat{\psi}$ where $\Psi_{0}$ is the mean-field
order parameter obeying equation (22) and $\delta\hat{\psi}$ represents small
fluctuations around the mean-field obeying

$\displaystyle\frac{\partial\delta\hat{\psi}(\mathbf{r},t)}{\partial t}$
$\displaystyle=\frac{i}{\hbar}\left[\frac{\hbar^{2}}{2m}\laplacian-
V_{\text{ext}}(\mathbf{r},t)-2g|\Psi_{0}|^{2}\right]\delta\hat{\psi}(\mathbf{r},t)$
$\displaystyle\qquad\qquad\qquad\quad-\frac{i}{\hbar}g\Psi_{0}^{2}\delta\hat{\psi}^{\dagger}(\mathbf{r},t).$
(38)

Hence, we find that the evolution equations of LSQHD are equivalent to the
stochastic Bogoliubov equations as well as the usual Bogoliubov approach,
which proceeds through the diagonalization of $\delta\hat{\psi}$ via the mode
decomposition
$\delta\hat{\psi}(\mathbf{r},t)=\sum_{i=1}^{\infty}(u_{i}(\mathbf{r},t)\hat{a}_{i}-v_{i}^{*}(\mathbf{r},t)\hat{a}_{i}^{\dagger})$
leading to the well-known time-dependent Bogoliubov–de Gennes (BdG) equations
[60, 61, 62], where $u_{i}$ and $v_{i}$ are the mode amplitudes and
$\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ the quasiparticle creation and
annihilation operators, respectively.

While the evolution equations of the LSQHD, stochastic Bogoliubov, and usual
Bogoliubov approaches are equivalent, the initial state of each requires
consideration. After diagonalization, the usual Bogoliubov approach begins in
an initial Bogoliubov state which evolves according to the usual time-
dependant BdG equations. By contrast, the stochastic schemes are usually
initialised in a coherent state, since a more accurate description of the
initial ground state of a Bose-Einstein condensate does not currently exist in
the literature for zero temperature positive-P approaches. This coherent state
initial condition is achieved by simply not seeding the initial state with any
noise. For the examples considered in this work, one begins with a desired
mean-field initial state described by $\rho_{0}(\mathbf{r},0)$ and
$\bm{v}_{0}(\mathbf{r},0)$ or $S_{0}(\mathbf{r},0)$ (equivalently
$\Psi_{0}(\mathbf{r},0)$), and sets all initial fluctuating fields to zero:
$\delta\rho(\mathbf{r},0)=\delta\bm{v}(\mathbf{r},0)=\delta S(\mathbf{r},0)=0$
or equally $\delta\psi(\mathbf{r},0)=\delta\tilde{\psi}(\mathbf{r},0)=0$. This
is equivalent to initialising the full positive-P equations with
$\rho(\mathbf{r},0)=\rho_{0}(\mathbf{r},0)=|\Psi_{0}(\mathbf{r},0)|^{2}$,
$\bm{v}(\mathbf{r},0)=\bm{v}_{0}(\mathbf{r},0)$,
$S(\mathbf{r},0)=S_{0}(\mathbf{r},0)$ or
$\Psi(\mathbf{r},0)=\Psi_{0}(\mathbf{r},0)$,
$\tilde{\Psi}(\mathbf{r},0)=\Psi_{0}^{*}(\mathbf{r},0)$. Such an
initialization avoids the need for diagonalization entirely. Hence, for
scenarios where exact diagonalization becomes computationally intractable, the
stochastic methods become invaluable. Some situations where this type of
advantage has already been exploited can be found in Refs. [45, 55, 57]. In
these cases, the systems being considered require an extremely large numerical
lattice or an unreasonable number of modes for simulation. Hence, exact
diagonalization in these types of situations is difficult and undesirable. An
equivalent stochastic approach, like the one derived in this work, provides a
means of avoiding this. Whilst we restrict ourselves to zero temperature here,
we note that it is also possible to consider finite temperature situations by
stochastically sampling the appropriate thermal P-distribution in order to
construct the initial state [63, 64].

Additionally, we mention that since the LSQHD equations are Bogoliubov in
nature (i.e. result from a quadratic Hamiltonian), this means that all higher
order correlations factorize into products of second order moments, according
to Wick’s theorem. At the level of field operators, this implies that we know
everything about the system after computing only the normal and anomalous
correlators,
$G^{(1)}_{n}(\mathbf{r},\mathbf{r}^{\prime},t)\equiv\langle\delta\hat{\psi}^{\dagger}(\mathbf{r},t)\delta\hat{\psi}(\mathbf{r}^{\prime},t)\rangle=\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$
and
$G^{(1)}_{a}(\mathbf{r},\mathbf{r}^{\prime},t)\equiv\langle\delta\hat{\psi}(\mathbf{r},t)\delta\hat{\psi}(\mathbf{r}^{\prime},t)\rangle=\langle\delta\psi(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$,
respectively. In the LSQHD formalism this is equivalent to having knowledge of
the correlators
$\langle\delta\rho(\mathbf{r},t)\delta\rho(\mathbf{r}^{\prime},t)\rangle$,
$\langle\delta S(\mathbf{r},t)\delta S(\mathbf{r}^{\prime},t)\rangle$, and
$\langle\delta\rho(\mathbf{r},t)\delta S(\mathbf{r}^{\prime},t)\rangle$. If
one uses the transformations (34) and (35) to convert the normal and anomalous
correlators to hydrodynamic variables then this results in the relations

$\displaystyle\langle$
$\displaystyle\delta\tilde{\psi}(\mathbf{r})\delta\psi(\mathbf{r}^{\prime})\rangle$
$\displaystyle=\sqrt{\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}e^{-i(S_{0}(\mathbf{r})-S_{0}(\mathbf{r}^{\prime}))}$
$\displaystyle\hskip
14.22636pt\times\left[\frac{\langle\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})\rangle}{4\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}+\langle\delta
S(\mathbf{r})\delta S(\mathbf{r}^{\prime})\rangle\right.$ $\displaystyle\hskip
31.29802pt+\left.\frac{i}{2}\left(\frac{\langle\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})\rangle}{\rho_{0}(\mathbf{r})}-\frac{\langle\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})\rangle}{\rho_{0}(\mathbf{r}^{\prime})}\right)\right]$
(39)

and

$\displaystyle\langle$
$\displaystyle\delta\psi(\mathbf{r})\delta\psi(\mathbf{r}^{\prime})\rangle$
$\displaystyle=\sqrt{\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}e^{i(S_{0}(\mathbf{r})+S_{0}(\mathbf{r}^{\prime}))}$
$\displaystyle\hskip
14.22636pt\times\left[\frac{\langle\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})\rangle}{4\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}-\langle\delta
S(\mathbf{r})\delta S(\mathbf{r}^{\prime})\rangle\right.$ $\displaystyle\hskip
31.29802pt+\left.\frac{i}{2}\left(\frac{\langle\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})\rangle}{\rho_{0}(\mathbf{r})}+\frac{\langle\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})\rangle}{\rho_{0}(\mathbf{r}^{\prime})}\right)\right]$
(40)

where the time index has been omitted for brevity.

The diagonals of these correlators are then given by

$\displaystyle\langle\delta\tilde{\psi}(\mathbf{r})$
$\displaystyle\delta\psi(\mathbf{r})\rangle=\rho_{0}(\mathbf{r})\left[\frac{\langle\delta\rho^{2}(\mathbf{r})\rangle}{4\rho_{0}^{2}(\mathbf{r})}+\langle\delta
S^{2}(\mathbf{r})\rangle\right],$ (41)
$\displaystyle\langle\delta\psi(\mathbf{r})$
$\displaystyle\delta\psi(\mathbf{r})\rangle=\rho_{0}(\mathbf{r})e^{2iS_{0}(\mathbf{r})}$
$\displaystyle\times\left[\frac{\langle\delta\rho^{2}(\mathbf{r})\rangle}{4\rho_{0}^{2}(\mathbf{r})}+i\frac{\langle\delta\rho(\mathbf{r})\delta
S(\mathbf{r})\rangle}{\rho_{0}(\mathbf{r})}-\langle\delta
S^{2}(\mathbf{r})\rangle\right].$ (42)

We take the opportunity now to point out a difference that arises between the
Bogoliubov schemes and the LSQHD approach which we explore further in Sec.
IV.2. From an entirely hydrodynamic perspective, one would compute the real-
space density within the LSQHD scheme using the diagonal of the reduced one-
body density matrix
$G^{(1)}(\mathbf{r},\mathbf{r},t)\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)\rangle$
as

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r},t)$
$\displaystyle=\langle\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r},t)\rangle$
$\displaystyle=\left\langle\rho(\mathbf{r},t)\right\rangle=\rho_{0}(\mathbf{r},t)+\left\langle\delta\rho(\mathbf{r},t)\right\rangle.$
(43)

Due to the additive nature of the LSQHD equations, if they are initialized in
a coherent state where $\delta\rho(\mathbf{r},0)=0$ and
$\delta\bm{v}(\mathbf{r},0)=0$, like we do in this work, then the average
density fluctuations will remain
$\left\langle\delta\rho(\mathbf{r},t)\right\rangle=0$. Hence we see that the
LSQHD approach can be considered a number conserving scheme which reproduces
only the mean-field result for the real-space density. Whilst the scheme is
still capable of computing non-trivial results for higher order correlation
functions, one must consider beyond-linear terms in the hydrodynamic equations
of motion to obtain any corrections to the mean-field density itself.

By contrast however, the usual Bogoliubov and stochastic Bogoliubov schemes
presented here do provide corrections to the mean-field density in the sense
that they are not number conserving. While more elaborate number conserving
approaches are possible [65, 62, 66, 57], where the number of particles in the
mean-field or condensate mode $\Psi_{0}$ remains constant, and the number of
particles in the excitations can vary. This can lead to dynamical population
growth in the total particle number density of the system where the initial
state of the fluctuating component is assumed to be a vacuum state. This is
similar to spontaneous parametric down-conversion in quantum optics in the
undepleted pump approximation, used as the simplest model for generation of
the paradigmatic squeezed vacuum state [67] (see also [68, 69], and references
therein).

Such population growth, along with the discrepancy between the Bogoliubov and
hydrodynamic approaches, can be seen by examining the reduced one-body density
matrix from the Bogoliubov perspective. In this case one arrives at

$\displaystyle G^{(1)}(\mathbf{r},$
$\displaystyle\mathbf{r},t)=\langle\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r},t)\rangle$
$\displaystyle\quad\,\,\,\,=|\Psi_{0}(\mathbf{r},t)|^{2}+\Psi_{0}^{*}(\mathbf{r},t)\langle\delta\psi(\mathbf{r},t)\rangle$
$\displaystyle+\Psi_{0}(\mathbf{r},t)\langle\delta\tilde{\psi}(\mathbf{r},t)\rangle+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle.$
(44)

Similarly to the LSQHD equations, since the stochastic Bogoliubov equations
are additive and we initialize them in a coherent state with
$\delta\tilde{\psi}(\mathbf{r},t)=\delta\psi(\mathbf{r},t)=0$, then the
average of these fields remain
$\langle\delta\tilde{\psi}(\mathbf{r},t)\rangle=\langle\delta\psi(\mathbf{r},t)\rangle=0$.
This means one is left with
$G^{(1)}(\mathbf{r},\mathbf{r},t)=|\Psi_{0}(\mathbf{r},t)|^{2}+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle$
for the real-space density, which contrasts
$G^{(1)}(\mathbf{r},\mathbf{r},t)=\rho_{0}(\mathbf{r},t)$ from Eq. (43). The
beyond-mean-field term here,
$\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle$, is
responsible for particle number growth in these approaches. It is clear then,
that inconsistencies will arise if one wants to convert Eq. (44) into
hydrodynamic variables using equation (41) where averages over second order
density and phase variables appear. Such an inconsistency arises due to the
fact that second order terms are often kept when dealing with observables like
the real-space density, in order to obtain leading order corrections to the
mean-field result. Yet the evolution equations themselves (36), (37) and (31),
(32) along with the transformations (34), (35) are truncated only to first
order.

###  III.4 Low-energy excitations and stochastic Luttinger liquid

Finally, the last situation we consider in this section is that of low-energy
excitations on top of an equilibrium density, which will lead to a Luttinger
liquid description of a weakly-interacting Bose gas in 1D. To treat such a
situation we begin with the full SQHD equations and write the stochastic
density variable as $\rho=\rho_{\text{eq}}+\delta\rho$, where
$\rho_{\text{eq}}$ represents the static equilibrium density and $\delta\rho$
any departure from that equilibrium value. We then proceed by linearizing
equations (14) and (15), assuming that the velocity $\bm{v}$ and the non-
equilibrium density $\delta\rho$ are small quantities. This leads to

$\displaystyle\frac{\partial\delta\rho}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\cdot\left[\rho_{\text{eq}}\bm{v}\right]+\sqrt{-i\frac{g}{\hbar}}\rho_{\text{eq}}\,(\xi_{1}+i\xi_{2}),$
(45) $\displaystyle\frac{\partial\bm{v}}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\left[\frac{\delta\tilde{\mu}}{m}+\frac{1}{2m}\sqrt{i\hbar
g}\,(\xi_{1}-i\xi_{2})\right],$ (46)

where

$\displaystyle\delta\tilde{\mu}=g\delta\rho-\frac{\hbar^{2}}{2m}\left[\frac{1}{2\sqrt{\rho_{\text{eq}}}}\laplacian\frac{\delta\rho}{\sqrt{\rho_{\text{eq}}}}-\frac{\delta\rho}{2\rho_{\text{eq}}^{3/2}}\laplacian\sqrt{\rho_{\text{eq}}}\right]$
(47)

comes from linearizing
$\tilde{\mu}=V_{\text{ext}}+g\rho-\frac{\hbar^{2}}{2m\sqrt{\rho}}\laplacian\sqrt{\rho}$.

Taking the expectation value of these equations results in

$\displaystyle\frac{\partial\expectationvalue{\delta\rho}}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\cdot\left[\rho_{\text{eq}}\expectationvalue{\bm{v}}\right],$
(48) $\displaystyle\frac{\partial\expectationvalue{\bm{v}}}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\left[\frac{\expectationvalue{\delta\tilde{\mu}}}{m}\right],$
(49)

which are precisely the equations of motion one obtains when considering low-
energy excitations on top of the superfluid hydrodynamic equations (18)–(20)
[16]. Computing the time derivative of (48) and eliminating the velocity using
(49) gives

$\displaystyle\frac{\partial^{2}\expectationvalue{\delta\rho}}{\partial
t^{2}}$
$\displaystyle\approx\bm{\nabla}\cdot\left[\rho_{\text{eq}}\bm{\nabla}\left(\frac{\expectationvalue{\delta\tilde{\mu}}}{m}\right)\right],$
(50)

which governs the dynamics of elementary excitations of Bose gases in
arbitrary potentials, and can be used to derive the Bogoliubov excitation
spectrum as well as the frequencies of collective oscillations in trapped
gases [16, 12].

Returning now to equations (45) and (46), we consider the case of a uniform
box where $\rho_{\text{eq}}=$ constant, resulting in
$\delta\tilde{\mu}=g\delta\rho-\frac{\hbar^{2}}{4m}\frac{1}{\rho_{\text{eq}}}\laplacian\delta\rho$.
Then, ignoring the effects of quantum pressure, the stochastic equations of
motion are reduced to

$\displaystyle\frac{\partial\delta\rho}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\cdot\left[\rho_{\text{eq}}\bm{v}\right]+\sqrt{-i\frac{g}{\hbar}}\rho_{\text{eq}}\,(\xi_{1}+i\xi_{2}),$
(51) $\displaystyle\frac{\partial\bm{v}}{\partial t}$
$\displaystyle\approx-\bm{\nabla}\left[\frac{g\delta\rho}{m}+\frac{1}{2m}\sqrt{i\hbar
g}\,(\xi_{1}-i\xi_{2})\right].$ (52)

Since it is the quantum pressure term(s) which contribute to the quadratic
part of the Bogoliubov spectrum, having ignored them restricts our
consideration to only the linear part of the spectrum. We point out however
that these equations of motion go beyond the usual superfluid description due
to their stochastic nature, and by virtue of the noise terms they can be used
to compute the non-equilibrium correlation functions of low-energy
excitations.

When considering the one-dimensional regime, equations (51) and (52) give

$\displaystyle\frac{\partial\delta\rho}{\partial t}$
$\displaystyle\approx-\frac{\partial}{\partial
x}\left[\rho_{\text{eq}}v\right]+\sqrt{-i\frac{g}{\hbar}}\rho_{\text{eq}}\,(\xi_{1}+i\xi_{2}),$
(53) $\displaystyle\frac{\partial v}{\partial t}$
$\displaystyle\approx-\frac{\partial}{\partial
x}\left[\frac{g\delta\rho}{m}+\frac{1}{2m}\sqrt{i\hbar
g}\,(\xi_{1}-i\xi_{2})\right],$ (54)

which we call the stochastic Luttinger liquid (SLL) equations for a weakly-
interacting Bose gas since they represent a stochastic version of the usual
Luttinger liquid equations of motion in the respective operator form [70, 49],

$\displaystyle\frac{\partial\delta\hat{\rho}}{\partial t}$
$\displaystyle=-\frac{\partial}{\partial
x}\left[\rho_{\text{eq}}\hat{v}\right],$ (55)
$\displaystyle\frac{\partial\hat{v}}{\partial t}$
$\displaystyle=-\frac{\partial}{\partial
x}\left[\frac{g}{m}\delta\hat{\rho}\right],$ (56)

with the understanding that quantum fluctuations in the SLL equations are
incorporated via the stochastic noise terms.

We mention also that the SLL equations can be derived directly from the LSQHD
equations (31) and (32) by assuming that the deterministic mean-field
variables are static (do not change in time), and then more specifically that
they describe a homogeneous gas, along with ignoring any terms that resulted
from the many-body quantum pressure. This corresponds formally to taking:
(_i_) $\rho_{0}=\rho_{\text{eq}}$, (_ii_) $\bm{v}_{0}=0$ with the
understanding that now $\delta\bm{v}\leftrightarrow\bm{v}$, and (_iii_)
ignoring the terms
$\frac{\hbar^{2}}{4m^{2}\sqrt{\rho_{0}}}\laplacian\left(\frac{\delta\rho}{\sqrt{\rho_{0}}}\right)$
and
$\frac{\hbar^{2}}{4m^{2}}\frac{\delta\rho}{\rho_{0}^{3/2}}\laplacian\sqrt{\rho_{0}}$
from equation (32).

##  IV Physical observables

###  IV.1 In the full stochastic quantum hydrodynamics

Here we consider a few typical physical observables and reformulate them
within the SQHD approach using equations (7) and (8). We recall that
observables described by averages over normally-ordered products of field
operators can be computed from their stochastic counterparts using Eq. (4).

Then, in this formulation, the reduced one-body density matrix and the
particle number density are given by

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r}^{\prime},t)\rangle=\langle\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle=\langle\sqrt{\rho(\mathbf{r},t)\rho(\mathbf{r}^{\prime},t)}\,e^{-i\left(S(\mathbf{r},t)-S(\mathbf{r}^{\prime},t)\right)^{\,}}\rangle$
(57)

and

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)\rangle=\langle\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r},t)\rangle$
$\displaystyle=\left\langle\rho(\mathbf{r},t)\right\rangle,$ (58)

respectively. The reduced one-body density matrix can then be normalized to
obtain

$\displaystyle g^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\frac{G^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)}{\sqrt{G^{(1)}(\mathbf{r},\mathbf{r},t)}\sqrt{G^{(1)}(\mathbf{r}^{\prime},\mathbf{r}^{\prime},t)}}$
$\displaystyle=\frac{\langle\sqrt{\rho(\mathbf{r},t)\rho(\mathbf{r}^{\prime},t)}\,e^{-i\left(S(\mathbf{r},t)-S(\mathbf{r}^{\prime},t)\right)^{\,}}\rangle}{\sqrt{\left\langle\rho(\mathbf{r},t)\right\rangle}\sqrt{\left\langle\rho(\mathbf{r}^{\prime},t)\right\rangle}}.$
(59)

Furthermore, the density-density correlation function (or the reduced two-body
density matrix) in normally-ordered form can be computed using

$\displaystyle G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime},t)\hat{\Psi}(\mathbf{r}^{\prime},t)\hat{\Psi}(\mathbf{r},t)\rangle$
$\displaystyle=\left\langle\rho(\mathbf{r},t)\rho(\mathbf{r}^{\prime},t)\right\rangle,$
(60)

and normalised as

$\displaystyle g^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\frac{G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)}{G^{(1)}(\mathbf{r},\mathbf{r},t)~{}G^{(1)}(\mathbf{r}^{\prime},\mathbf{r}^{\prime},t)}$
$\displaystyle=\frac{\left\langle\rho(\mathbf{r},t)\rho(\mathbf{r}^{\prime},t)\right\rangle}{\left\langle\rho(\mathbf{r},t)\right\rangle\left\langle\rho(\mathbf{r}^{\prime},t)\right\rangle}.$
(61)

Note that, in particular, the reduced one-body density matrix (59) is easily
represented in terms of the density and phase variables (rather than in terms
of the density-velocity pair) and therefore the respective set of SDEs appears
to be a more natural one to solve numerically for this observable.

###  IV.2 In the linearized treatment of quantum fluctuations

In the LSQHD approach, the same quantities as above can be evaluated as
follows.

The reduced one-body density matrix is given by

$\displaystyle G$ ${}^{(1)}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle\approx\sqrt{\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}e^{-i\left[S_{0}(\mathbf{r})-S_{0}(\mathbf{r}^{\prime})\right]}$
$\displaystyle\quad\times\left(1+\expectationvalue{\delta S(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}-\frac{1}{2}\left[\expectationvalue{\delta
S^{2}(\mathbf{r})}+\expectationvalue{\delta
S^{2}(\mathbf{r}^{\prime})}\right]\right.$
$\displaystyle\qquad\quad+\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{4\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}-\frac{1}{8}\left[\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r})}}{\rho_{0}^{2}(\mathbf{r})}+\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r}^{\prime})}}{\rho_{0}^{2}(\mathbf{r}^{\prime})}\right]$
$\displaystyle\qquad\quad+\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}\right]$
$\displaystyle\qquad+\left.\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r}^{\prime})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r})}}{\rho_{0}(\mathbf{r})}\right]\right),$
(62)

where the time index has been dropped for brevity and averages of first order
terms have been set to zero as per the discussion in Sec. III.3.

The particle number density itself is then simply given by

$G^{(1)}(\mathbf{r},\mathbf{r},t)=\rho_{0}(\mathbf{r},t),$ (63)

and normalizing the reduced one-body density matrix results in

$\displaystyle g$ ${}^{(1)}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle\approx
e^{-i\left[S_{0}(\mathbf{r})-S_{0}(\mathbf{r}^{\prime})\right]}$
$\displaystyle\quad\times\left(1+\expectationvalue{\delta S(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}-\frac{1}{2}\left[\expectationvalue{\delta
S^{2}(\mathbf{r})}+\expectationvalue{\delta
S^{2}(\mathbf{r}^{\prime})}\right]\right.$
$\displaystyle\qquad\quad+\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{4\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}-\frac{1}{8}\left[\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r})}}{\rho_{0}^{2}(\mathbf{r})}+\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r}^{\prime})}}{\rho_{0}^{2}(\mathbf{r}^{\prime})}\right]$
$\displaystyle\qquad\quad+\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}\right]$
$\displaystyle\quad\qquad+\left.\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r}^{\prime})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r})}}{\rho_{0}(\mathbf{r})}\right]\right),$
(64)

where we drop the time index again for brevity.

Furthermore, the density-density correlation function in normally-ordered form
can be computed using

$G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)=\rho_{0}(\mathbf{r},t)\rho_{0}(\mathbf{r}^{\prime},t)+\left\langle\delta\rho(\mathbf{r},t)\delta\rho(\mathbf{r}^{\prime},t)\right\rangle,$
(65)

and normalised as

$\displaystyle g^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\frac{G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)}{G^{(1)}(\mathbf{r},\mathbf{r},t)~{}G^{(1)}(\mathbf{r}^{\prime},\mathbf{r}^{\prime},t)}$
$\displaystyle=1+\frac{\left\langle\delta\rho(\mathbf{r},t)\delta\rho(\mathbf{r}^{\prime},t)\right\rangle}{\rho_{0}(\mathbf{r},t)\rho_{0}(\mathbf{r}^{\prime},t)}.$
(66)

In the stochastic Bogoliubov approach, on the other hand, the reduced one-body
density matrix can be computed using

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle=\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)+\Psi_{0}^{*}(\mathbf{r},t)\langle\delta\psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle\qquad+\Psi_{0}(\mathbf{r}^{\prime},t)\langle\delta\tilde{\psi}(\mathbf{r},t)\rangle+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle=\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle,$
(67)

where averages of first order terms have been set to zero as per the
discussion in Sec. III.3.

The particle number density itself is then given by

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r},t)$
$\displaystyle=\left|\Psi_{0}(\mathbf{r},t)\right|^{2}+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle,$
(68)

and the normalized reduced one-body density matrix is (see Appendix D)

$\displaystyle g$ ${}^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\approx\sqrt{\frac{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}{\Psi_{0}(\mathbf{r},t)\Psi_{0}^{*}(\mathbf{r}^{\prime},t)}}\left(1+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}\right.$
$\displaystyle\,\,\,-\left.\frac{1}{2}\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r},t)}-\frac{1}{2}\frac{\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\Psi_{0}(\mathbf{r}^{\prime},t)}\right).$
(69)

As mentioned previously in Sec. III.3 regarding the real-space density, a
similar discrepancy exists between the normalized reduced one-body density
matrices (64) and (69). These expressions are not equivalent under the
transformations (34) and (35) due to the inconsistency of maintaining second
order terms in these expressions (which are needed to obtain leading order
corrections to the mean-field results) while the equations of motion
themselves have been truncated at linear order. A similar type of
inconsistency has been identified in Ref. [71] when developing an extension of
Bogoliubov theory to quasi-condensates. We explore these differences between
the LSQHD and stochastic Bogoliubov predictions in Appendix F, where we show
that differences in the density are noticeable, whereas differences in
$g^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$ are not.

Finally, the normalized density-density correlation function is given by (see
Appendix D)

$\displaystyle g^{(2)}(\mathbf{r},$
$\displaystyle\mathbf{r}^{\prime},t)\approx
1+2\Re{\frac{\langle\delta\psi(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}}$
$\displaystyle+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}+\frac{\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r},t)\rangle}{\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\Psi_{0}(\mathbf{r},t)}.$
(70)

We point out that whilst the expressions for the normalized first order
correlators (64) and (69) are inconsistent, the normalized second order
correlators (66) and (70) are in fact equivalent under the transformations
(34) and (35).

##  V Simulations & Comparisons with existing methods

###  V.1 Density-density correlations in quantum shock waves

As a first application of the SQHD formalism to systems of interacting bosons,
we begin by exploring the correlations arising in a scenario which develops
dispersive quantum shock waves. In particular we focus on the set-up
previously investigated in Ref. [32]: namely, an initial density bump on top
of a uniform non-zero background in a one-dimensional (1D) Bose gas.
Considering two examples with different particle number $N$, yet with
identical 1D interaction strength $g_{\text{1D}}$, we compare the predictions
made by SQHD with other well-known stochastic methods, as well as exact
quantum calculations using matrix product states (MPS).

In each example we consider an initial density profile of the form
$G^{(1)}(x,x,0)=\rho(x,0)=N_{\text{bg}}[1+\beta e^{-x^{2}/2\sigma^{2}}]^{2}/L$
which subsequently evolves in a uniform potential of length $L$ with periodic
boundary conditions. Here $N_{\text{bg}}=\rho_{\text{bg}}L$ specifies the
number of particles in the background, related to the total particle number
via
$N=N_{\text{bg}}(1+\frac{\sqrt{\pi}\beta\sigma}{L}[\beta\erf(\frac{L}{2\sigma})+2\sqrt{2}\erf(\frac{L}{2\sqrt{2}\sigma})])$.
The amplitude of the bump above the background density is determined by
$\beta$ and its width is governed by $\sigma$.

For this scenario one can define a local dimensionless interaction parameter
$\gamma(x)=mg_{\text{1D}}/\hbar^{2}\rho(x)$ (or a local Lieb-Liniger parameter
[33, 72]), and then use its value at the background density
$\gamma_{\text{bg}}=mg_{\text{1D}}/\hbar^{2}\rho_{\text{bg}}$ to characterize
the initial state. Such a parameter then allows for a simple definition of the
characteristic length-scale associated with interactions—the dimensionless
healing length $l_{\text{h}}/L=1/N_{\text{bg}}\sqrt{\gamma_{\text{bg}}}$,
where we use $L$ as the length scale.

####  V.1.1 Example 1: $\sigma/l_{h}\simeq 0.11$
($N=50$, $\gamma_{\text{bg}}=0.1$, $\bar{g}_{\text{1D}}\simeq 4.773$)

In this first example we consider a weakly-interacting 1D Bose gas of $N=50$
particles, with $\gamma_{\text{bg}}=0.1$. The initial bump is characterized by
$\beta=1$ and $\sigma/L=0.007$. This configuration leads to
$N_{\text{bg}}\simeq 47.73$ particles in the background, and therefore a
background healing length of $l_{\text{h}}/L\simeq 0.066$. As such, the bump
width dominates the dynamics as the shortest length-scale in the problem [32],
being only about 10% of the background healing length, $\sigma\simeq
0.11l_{\text{h}}$. For completeness, the dimensionless interaction strength
for this scenario is given by
$\bar{g}_{\text{1D}}=g_{\text{1D}}mL/\hbar^{2}=\gamma_{\text{bg}}N_{\text{bg}}\simeq
4.773$.

In the following figures we show predictions for the dimensionless real-space
density $G^{(1)}(x,x,t)L=\rho(x,t)L$ and the normally-ordered density-density
correlation function $g^{(2)}(x,x^{\prime},t)$ that result from a number of
different quantum many-body approaches. These include; the full stochastic
positive-P equations, the LSQHD equations, and the truncated Wigner approach
(see Appendix A for details), along with an infinite-MPS calculation (see
Appendix C for details). We have also simulated the stochastic Bogoliubov
equations, but have omitted these results from the figures below for clarity
since we recall that the respective density-density correlation function
$g^{(2)}(x,x^{\prime},t)$ is equivalent to that of LSQHD, whereas results for
the predicted real-space density can be found in Appendix F.

Figure 1: The real-space density $\rho(x,t)$ of a weakly interacting 1D Bose
gas in the quantum shock wave scenario of Ref. [32]. Shown are the initial
($\tau=0$) and time-evolved ($\tau=0.001$) dimensionless density profiles
$\rho(x,t)L$, where the dimensionless time $\tau$ is introduced according to
$\tau=t\hbar/mL^{2}$, for $N=50$ particles and a dimensionless background
interaction strength of $\gamma_{\text{bg}}=0.1$. The initial state (dashed
lines) has a width $\sigma/L=0.007$ and amplitude $\beta=1$. This
configuration leads to $N_{\text{bg}}\simeq 47.73$ particles in the
background, and a background healing length of $l_{\text{h}}/L\simeq 0.066$.
The results displayed are for: the full stochastic positive-P equations (Full
+P), the linearized stochastic hydrodynamic equations (LSQHD), the truncated
Wigner approach (TWA), and a matrix product state calculation initialized in a
coherent state (CS-MPS), with all the different curves being essentially
indistinguishable from each other on this scale. Other technical details on
simulations are provided in the caption to Table 1 below.

Each of the approaches used here have been initialized in a coherent state
(see Appendices A, B and C) for a fair comparison with the positive-P schemes,
since a more suitable approximation to the initial many-body ground state
within this formalism does not currently exist in the literature. Out of
interest, in Appendix E we compare these initial coherent state results with
an infinite-MPS calculation initialized in the exact many-body ground state of
a trapping potential which is tailored to closely reproduce the same initial
density $\rho(x,0)=N_{\text{bg}}[1+\beta e^{-x^{2}/2\sigma^{2}}]^{2}/L$ as the
mean-field (coherent state) approach.

Figure 2: The density-density correlation function $g^{(2)}(x,x^{\prime},t)$
from LSQHD at dimensionless time $\tau=0.001$ for the same parameters as in
Fig. 1.

In Figure 1 we show the real-space density at the initial dimensionless time
$\tau=t\hbar/mL^{2}=0$ and at time $\tau=0.001$. When released into the
uniform trap, the initial density bump begins to split into left and right
moving parts, where each quickly develops into a dispersive quantum shock
wave—an oscillatory wave-train that results from quantum mechanical self-
interference [32]. Unlike each of the other methods presented here, we recall
that the LSQHD scheme does not incorporate the effect of quantum fluctuations
into the density, i.e. it reproduces only the mean-field superfluid result for
this observable. Nevertheless the densities predicted by each method in Fig. 1
agree very well, and on the scale plotted here they lie almost directly on top
of each other making them virtually indistinguishable. Hence, the effect of
quantum fluctuations on the real-space density (which is known to reduce the
amplitude of the interference contrast from that of the pure mean-field GPE
prediction [32]) is negligible in this scenario. Their limited effect on the
density is discussed in Appendix E.

Turning now to examine the quantity of most interest in this section, in
Figure 2 we plot the entire normalized density-density correlation function
$g^{(2)}(x,x^{\prime},t)$ at time $\tau=0.001$, as predicted by LSQHD. We then
proceed by exploring and comparing relevant slices of this correlation
function, demonstrating the advantages and usefulness of the LSQHD approach.

Figure 3: The same-point correlation function $g^{(2)}(x,x,t)$ for the
scenario presented in Fig. 1. (a) shows the initial correlations at $\tau=0$
which are precisely $g^{(2)}(x,x,0)=1$ for each method, except for the
truncated Wigner approach (TWA) which contains initial fluctuations around
this value that go to zero in the limit of infinitely many stochastic
trajectories. (b) provides a comparison of $g^{(2)}(x,x,\tau=0.001)$ for each
method (without any uncertainty shown). For clarity, (c) shows the same data
as in (b) spaced upward by 0.15 from the coherent-state MPS (CS-MPS) result.
The shaded regions in (a) and (c) denote one standard error of uncertainty for
the respective stochastic approaches, which for linearized stochastic quantum
hydrodynamics (LSQHD) grows dynamically to about half the size of the
linewidth used here and can be seen between the dashes in (c). Figure 4: The
opposite-point correlation function $g^{(2)}(x,-x,t)$ for the scenario
presented in Fig. 1. To examine the correlations resulting from the density
bump itself we plot between $|x|/L=0.25$. Due to the periodic boundary
conditions of the system, values outside of this domain correspond to probing
the opposite-point correlations of the background. The results are presented
similarly to Fig. 3 where the upward shift in (c) is now only 0.1.

In Figures 3 and 4, for each of the aforementioned approaches, we plot the
same-point $g^{(2)}(x,x,t)$ and opposite-point $g^{(2)}(x,-x,t)$ correlation
functions respectively. These two correlation functions develop oscillations
during evolution which coincide with those in the real-space density, where
$g^{(2)}(x,-x,t)$ oscillates twice for every oscillation in the density.
Furthermore, in Table 1 we show the time taken to simulate this scenario for
each approach, where identical lattice spacing and time stepping was used for
each of the stochastic schemes.

The first advantage of the LSQHD approach here is readily apparent in the
same- and opposite-point correlations: it provides sufficiently smooth results
that lie close to the exact coherent-state MPS (CS-MPS) predictions, and it
does not suffer from excessive noise like the full positive-P and truncated
Wigner approaches. More specifically, the LSQHD approach possesses noise that
is $\sim\\!17$ times smaller in magnitude compared to the full positive-P
approach, and $\sim\\!54$ times smaller than the noise in the truncated Wigner
result. This means that reducing the Full +P and TWA error bars down to the
same level as for LSQHD would require running $\sim\\!300$ and $\sim\\!3000$
times more stochastic trajectories for each respective approach.

The second advantage we mention is the reduced simulation time compared with
CS-MPS. The LSQHD simulation was approximately $3.3$ times faster than the
exact CS-MPS calculation and yet provides essentially the same results for
each of the correlation functions. Moreover, the LSQHD scheme was similarly
faster to simulate than the stochastic Bogoliubov approach. Whilst the
$g^{(2)}(x,x^{\prime},t)$ results of the stochastic Bogoliubov equations lie
directly underneath the LSQHD result (within one standard error), they took
$\sim 1.5$ times longer to simulate. For clarity, we report that no elaborate
code optimization was performed for any of the stochastic simulations reported
in Table 1 and that they were each carried out on an Apple iMac desktop
hosting a 4.2GHz Intel Core i7 (i7-7700K) processor with 4 cores and 64GB of
RAM, whereas the MPS simulations were performed using the Matrix Product
Toolkit of Ref. [73] on an Apple iMac desktop hosting a 3.5GHz Intel Core i7
(i7-4771) processor with 4 cores and 32GB of RAM.

Finally, we examine the applicability of the linearized positive-P schemes and
the truncated Wigner approach for this scenario. A check of the population
growth present in the stochastic Bogoliubov scheme reveals an increase in
particle number of $\sim 2\%$ by the final time of the simulation, confirming
the validity of the undepleted pump approximation and the reasonableness of
using the linearized positive-P schemes here. The usual requirement for
applicability of the truncated Wigner approach, on the other hand, is that the
number of particles $N$ be much larger than the number of modes $M$, i.e. the
lattice mode occupation should be $N/M\gg 1$ [74]. However in this case it is
only $N/M=50/300\simeq 0.17$ and hence the a priori applicability of this
approach is questionable. Nevertheless, we note that it still provides
reasonable predictions for the density and $g^{(2)}(x,x^{\prime},t)$ (even
though it suffers from excessive noise). This supports the notion that the
condition $N/M\gg 1$ can be somewhat relaxed in 1D to $N/M\sim 1$ [75, 76,
77], or even a smaller value as is the case here.

Table 1: Simulation times for the results presented in Figs. 1-4. The stochastic simulations were run using the XMDS software of Ref. [78]. For each stochastic approach the results we present here are an average over 2,000,000 stochastic trajectories which were run in parallel using the “mpi-multi-path” driver, with 300 spacial lattice points and a time step of $\Delta\tau=3.3333\times 10^{-6}$. Simulation outputs were $\rho(x,t)$ and $g^{(2)}(x,x^{\prime},t)$, along with their associated standard errors. The positive-P schemes were integrated using the SI (semi-implicit) algorithm [79, 80] and the truncated Wigner scheme was integrated using the RK4 (fourth-order Runge-Kutta) algorithm. Method |  Time to simulate
---|---
Full $+P$ |  36 hours
LSQHD |  37 hours
Stoch. Bogoliubov |  54 hours
TWA |  39 hours
CS-MPS |  124 hours

####  V.1.2 Example 2: $\sigma/l_{h}\simeq 0.33$
($N=500$, $\gamma_{\text{bg}}=0.01$, $\bar{g}_{\text{1D}}\simeq 4.773$)

In this second example, similarly to the first, we also consider an initial
bump characterized by $\beta=1$ and $\sigma/L=0.007$. Here however, we
increase the total particle number by an order of magnitude to $N=500$, while
maintaining an equal dimensionless interaction strength
$\bar{g}_{\text{1D}}=g_{\text{1D}}mL/\hbar^{2}=\gamma_{\text{bg}}N_{\text{bg}}\simeq
4.773$. Such a configuration has a background interaction parameter of
$\gamma_{\text{bg}}=0.01$ and $N_{\text{bg}}\simeq 477.3$ particles in the
homogeneous background. This leads to a shorter background healing length of
$l_{\text{h}}/L=1/N_{\text{bg}}\sqrt{\gamma_{\text{bg}}}\simeq 0.021$, where
the initial bump width is now about 30% of the background healing length,
$\sigma/l_{\text{h}}\simeq 0.33$. Whilst only a small increase in the
comparability of these two length scales, the increased role of interactions
in this scenario brings about a clear qualitative change in the correlations
of the system, as we will discuss in a moment.

The usefulness of the LSQHD approach is illustrated in this scenario in that
it is the only one of a few quantum many-body approaches that remains
tractable. The MPS approach becomes intractable here due to the large number
of particles $N=500$, and the multiplicative nature of the noise in the full
positive-P treatment causes this approach to break down halfway through the
full simulation time window, at
$\tau_{{break}}=t_{{break}}\hbar/mL^{2}=0.0005$. This leaves only the
linearized stochastic schemes and the truncated Wigner approach, each of which
are initialized in a coherent state with 300 lattice points, and averages
taken over 500,000 trajectories. Here we present a comparison of these methods
at the final time of the simulation $\tau=0.001$.

Figure 5: Density profiles in a quantum shock wave scenario as in Fig. 1, but
for $N=500$ particles and a dimensionless background interaction strength of
$\gamma_{\text{bg}}=0.01$. These parameters lead to $N_{\text{bg}}\simeq
477.3$, and a background healing length of $l_{\text{h}}/L\simeq 0.021$. The
results displayed are for the LSQHD and TWA approaches. The inset shows more
clearly the effect of quantum fluctuations on the first few interference
fringes; we recall that quantum fluctuations are accounted for in the TWA
result, but are absent in the density of LSQHD such that $\rho(x,t)$ is given
simply by the mean-field prediction
$\rho(x,t)=\rho_{0}(x,t)=|\psi_{0}(x,t)|^{2}$ of superfluid hydrodynamics [see
Eq. (63)].

In Figure 5 we show the evolution of the real-space density. The result is
similar as in Example 1: the initial bump splits into left and right moving
parts which quickly develop into dispersive quantum shock waves. Again, the
density predicted by each method agrees well with each other. In contrast to
Example 1 however, by the time $\tau=0.001$, we see that the left and right
moving parts have almost completely separated, and that the effect of quantum
fluctuations on the density is already becoming apparent on the scale plotted
here – the amplitude of oscillations in the truncated Wigner result is
slightly reduced compared to the LSQHD scheme which reproduces the mean-field
density. This reduced interference contrast can be seen more clearly in the
inset of Fig. 5.

Figure 6: The density-density correlation function $g^{(2)}(x,x^{\prime},t)$
from LSQHD at dimensionless time $\tau=0.001$ for the same parameters as in
Figure 5.

In Figure 6, we plot the entire normalized density-density correlation
function $g^{(2)}(x,x^{\prime})$ at time $\tau=0.001$, as predicted by LSQHD.
Similarly, we proceed by examining the same correlation slices as before,
comparing the linearized stochastic schemes with the truncated Wigner approach
and further demonstrating the utility of the LSQHD equations.

Figure 7: The same-point correlation function $g^{(2)}(x,x,t)$ as in Fig. 3,
but for the parameters of Fig. 5. The only computationally tractable
approaches in this regime are the LSQHD and TWA schemes (see text). For
clarity, (c) shows the same data as in (b) spaced upward by 0.1 from the TWA
result. The shaded regions in (a) and (c) denote one standard error of
uncertainty for the respective stochastic approaches, which for LSQHD grows
dynamically to about the size of the linewidth used here and can be seen
between the dashes in (c). Figure 8: The opposite-point correlation function
$g^{(2)}(x,-x,t)$ as in Fig. 4, but for the parameters of Fig. 5. To examine
the correlations resulting from the density bump itself we plot between
$|x|/L=0.25$. Due to the periodic boundary conditions of the system, values
outside of this domain correspond to probing the opposite-point correlations
of the background. The results are presented similarly to Fig. 7 where the
upward shift of the LSQHD curve in (c) is now only 0.05.

In Figures 7 and 8 we plot the same- and opposite-point correlation functions
respectively. Again, the oscillations which develop coincide with those in the
real-space density. However, in this case, the oscillations in
$g^{(2)}(x,x,t)$ sit predominantly above the background level, and now
oscillate with twice the frequency. Here the predictions of these correlations
in the LSQHD (equally, stochastic Bogoliubov) and truncated Wigner approaches
agree well with each other, where the linearized schemes slightly overestimate
the oscillation amplitudes in comparison to the truncated Wigner result.

Table 2: Simulation times for the results presented in Figs. 5-7. For each stochastic approach the results we present here are an average over 500,000 stochastic trajectories; other details are as in Table 1. Method |  Time to simulate
---|---
Full $+P$ |  2.3 hours (until $\tau_{{break}}$)
LSQHD |  6.3 hours
Stoch. Bogoliubov |  9.2 hours
TWA |  3.6 hours

Table 2 shows the time taken to simulate this scenario for each approach.
Since the full positive-P approach breaks down, we show only the simulation
time to reach $\tau_{{break}}$. Here again for this example we find that the
LSQHD equations are $\sim\\!1.5$ faster to simulate than the stochastic
Bogoliubov approach. On the other hand, the truncated Wigner approach took a
considerably shorter amount of time to simulate than LSQHD. However, we point
out that only one complex deterministic equation must be solved in the
truncated Wigner approach whereas the LSQHD approach requires the simultaneous
integration of two real deterministic equations and two complex stochastic
equations. In light of this, it is quite remarkable then that in Table 1 the
LSQHD equations were faster to simulate than the truncated Wigner approach.
Though we must also remark that no effort has been made to optimize the
integration of equations in any of these approaches. Moreover, the truncated
Wigner result presented here still contains more than twice as much noise on
average as the LSQHD result, and as such would require the simulation of at
least four times as many stochastic trajectories to reduce the noise to a
similar level.

Lastly, we check the applicability conditions of the approaches used here. The
population growth present in the stochastic Bogoliubov scheme reveals an
increase in particle number of $\sim 11\%$ by the final time of the
simulation, which still remains a reasonable application of the undepleted
pump approximation. This is supported by the predicted correlation functions
which do not depart far from unity, suggesting that the system remains in the
regime of small fluctuations where both linearized positive-P schemes are
valid. The lattice mode occupation of the truncated Wigner approach here is
$N/M=500/300\simeq 1.7$. Being larger than one, and considering the 1D nature
of the system, we consider this a reasonable a priori application of the
truncated Wigner approach.

##  VI Conclusions

In summary, we have derived a set of SQHD equations, using the positive-P
phase-space approach, that are capable of describing the exact dynamics of
trapped Bose gases, and shown that they reduce to other well known approaches
under the appropriate approximations. Linearizing such equations allows us to
derive the LSQHD scheme which does not suffer from uncontrollable growth of
sampling errors at late times and we provide connections between this scheme
and Bogoliubov theories. By exploring the density-density correlation
functions of quantum shock waves, which develop from a density bump expanding
into a non-zero background, we have demonstrated the usefulness and advantages
of the LSQHD scheme through a thorough comparison with other dynamical quantum
many-body approaches.

One limitation of the LSQHD scheme however, is that it is not suitable for
situations where the mean-field density $\rho_{0}(\bm{r},t)$ goes to zero
because of the division by $\rho_{0}(\bm{r},t)$ in the quantum pressure terms
of Eqs. (19), (20), (32), and (33). We point out that this limitation is not
restricted to the LSQHD scheme, but is inherent in the standard superfluid
hydrodynamic approach itself which suffers from the same problem when the
mean-field quantum pressure term is kept. Nevertheless, other interesting
dynamical scenarios for which the LSQHD scheme would be an ideal and preferred
method include: the evolution of a density dip (which sheds gray solitons as
it expands into the background [25, 81, 82, 29, 83]), a combination of density
bumps and dips, and density profiles involving step-like functions.

Future outlooks of this work involve going beyond the LSQHD scheme to consider
higher-order terms in the expansion like those included in Eqs. (28)–(30) and
the effect of the many-body quantum pressure term. For example, it would be
extremely interesting and insightful to explore possible new constitutive
relations in the gradient expansion scheme of conventional hydrodynamics [84],
that might follow from the higher-order expansion terms. Another avenue is to
find connections with the theory of nonlinear Luttinger liquids [85, 86] or
more generally with the theories of nonlinear fluctuating hydrodynamics in 1D
[87]. In these theories, nonlinear terms are important for describing
anomalous transport in 1D classical fluids and alter the universality class of
hydrodynamics fluctuations [87, 88]. We further emphasise here that while in
phenomenological theories of fluctuating hydrodynamics the noise and
dissipation terms are added to Euler-scale hydrodynamics essentially by hand,
the theory of stochastic quantum hydrodynamics derived here is fully
microscopic—with the noise terms representing the intrinsic _quantum_
fluctuations. Accordingly, exploring its connections to phenomenological
theories may provide new insights into the mechanisms behind the emergence of
anomalous transport in lower dimensions and the universality classes of
associated fluctuations.

###### Acknowledgements.

The authors acknowledge stimulating discussions with M. J. Davis and D. M.
Gangardt. This work was supported by the Australian Research Council Discovery
Project Grants DP170101423 and DP190101515.

##  Appendix A The truncated Wigner approach

Similarly to the positive-P treatment, the truncated Wigner approach is a
stochastic phase-space approach which relies on the conversion of the master
equation for the quantum density operator to an evolution equation for the
Wigner quasiprobability distribution function. For the Wigner distribution
function, however, this evolution equation, with a quartic (two-body)
interaction Hamiltonian (such as the case here), contains third-order
derivative terms. Such higher-than-second-order derivative terms must be
truncated in order to render the evolution equation a true Fokker-Planck
equation, which can then be equivalently formulated and simulated in terms of
stochastic differential equations for complex $c$-fields [89, 40]. As such,
this approach is only approximate, in contrast to the exact formulation of the
positive-P treatment, and care must be taken to ensure that the contribution
of the truncated terms is small (usually this is true when the lattice mode
occupation of the system is large, i.e. the number of particles $N$ is much
larger than the number of required modes $M$ such that $N/M\gg 1$ [74]). On
the other hand, this approach has the advantage that it does not suffer from
dynamical noise growth like the positive-P treatment does.

In the case of the Hamiltonian (1) with s-wave scattering interactions, the
resulting stochastic differential equations actually take the same form as the
standard time-dependent GPE. Meanwhile, beyond-mean-field quantum effects are
incorporated via the addition of noise into the initial state of the system
$\psi_{W}(\mathbf{r},0)$, which is then sampled stochastically, where each
realization evolves in time according to the GPE

$\displaystyle\frac{\partial\psi_{W}}{\partial t}$
$\displaystyle=\frac{i}{\hbar}\left[\frac{\hbar^{2}}{2m}\laplacian-
V_{\text{ext}}(\mathbf{r},t)-g_{\text{1D}}|\psi_{W}|^{2}\right]\psi_{W}.$ (71)

For the scenarios presented in Section V.1, each trajectory is initialized in
a coherent-state using [40]

$\displaystyle\psi_{W}(x,0)$
$\displaystyle=\Psi_{0}+\sum_{j=0}^{M}\eta_{j}\phi_{j}(x)$ (72)

where $\Psi_{0}=\sqrt{N_{\text{bg}}}[1+\beta e^{-x^{2}/2\sigma^{2}}]/\sqrt{L}$
is the mean-field initial state, $\eta_{j}$ is a complex Gaussian noise term
with zero mean $\langle\eta_{j}\rangle_{stoch.}=0$ and non-zero correlations
$\langle\eta_{j}^{*}\eta_{k}\rangle_{stoch.}=\frac{1}{2}\delta_{jk}$. The sum
is taken over $M$ modes of a complete basis $\\{\phi_{j}(x)\\}$, where $M$ is
chosen to be sufficiently large such that the results are independent of $M$.
Here we choose the discreet position basis, with $\phi_{j}(x)=\delta(x-j\Delta
x)/\sqrt{\Delta x}$.

In the truncated Wigner formalism, we note that stochastic averages correspond
to expectation values of symmetrically ordered products of field creation and
annihilation operators. As such, the reduced one-body density matrix (which is
normally-ordered) is given by

$\displaystyle G^{(1)}(x,x^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x^{\prime},t)\rangle$
$\displaystyle=\langle\psi_{W}^{*}(x,t)\psi_{W}(x^{\prime},t)\rangle-\frac{1}{2}\delta_{c}(x,x^{\prime})$
(73)

where
$\frac{1}{2}\int_{-L/2}^{L/2}\delta_{c}(x,x^{\prime})dx=\frac{1}{2}M\delta_{x,x^{\prime}}$
represents the half quantum of vacuum noise per mode $M$ that is included in
the Wigner formalism. Here $\delta_{x,x^{\prime}}$ is the Kronecker delta
function which gives unity for $x=x^{\prime}$ and zero otherwise. On a
computational grid of spacing $\Delta x=L/\mathcal{N}_{x}$, where
$\mathcal{N}_{x}$ is the number of grid points, the projected delta function
$\delta_{c}(x,x^{\prime})$ is given by
$\delta_{c}(x,x^{\prime})=M\delta_{x,x^{\prime}}/(\mathcal{N}_{x}\Delta
x)=M\delta_{x,x^{\prime}}/L$ [89]. This means that the reduced one-body
density matrix can be computed using

$\displaystyle G^{(1)}(x,x^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x^{\prime},t)\rangle$
$\displaystyle=\langle\psi_{W}^{*}(x,t)\psi_{W}(x^{\prime},t)\rangle-\frac{M}{2L}\delta_{x,x^{\prime}},$
(74)

and the particle number density using

$\displaystyle G^{(1)}(x,x,t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}(x,t)\rangle$
$\displaystyle=\langle|\psi_{W}^{*}(x,t)|^{2}\rangle-\frac{M}{2L}.$ (75)

Similarly, the normally ordered density-density correlation function can be
computed as [90]

$\displaystyle G^{(2)}(x,x^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(x,t)\hat{\Psi}^{\dagger}(x^{\prime},t)\hat{\Psi}(x^{\prime},t)\hat{\Psi}(x,t)\rangle$
$\displaystyle=\langle|\psi_{W}(x,t)|^{2}|\psi_{W}(x^{\prime},t)|^{2}\rangle$
$\displaystyle\quad-\frac{M}{2L}(1+\delta_{x,x^{\prime}})\bigg{[}\langle|\psi_{W}(x,t)|^{2}\rangle$
$\displaystyle\qquad\qquad\qquad+\left.\langle|\psi_{W}(x^{\prime},t)|^{2}\rangle-\frac{M}{2L}\right],$
(76)

and then normalised via

$\displaystyle g^{(2)}(x,x^{\prime},t)$
$\displaystyle\equiv\frac{G^{(2)}(x,x^{\prime},t)}{G^{(1)}(x,x,t)~{}G^{(1)}(x^{\prime},x^{\prime},t)}.$
(77)

##  Appendix B Initial states in the positive-P approaches

As mentioned in Section III.3 of the main text, due to the normally-ordered
nature of the positive-P approach, an initial coherent state is described
simply by not seeding the initial state with any noise.

In the full positive-P approach, this amounts to setting
$\rho(\mathbf{r},0)=\rho_{0}(\mathbf{r},0)$,
$\bm{v}(\mathbf{r},0)=\bm{v}_{0}(\mathbf{r},0)$, and
$S(\mathbf{r},0)=S_{0}(\mathbf{r},0)$, or equivalently,
$\Psi(\mathbf{r},0)=\Psi_{0}(\mathbf{r},0)$ and
$\tilde{\Psi}(\mathbf{r},0)=\Psi_{0}^{*}(\mathbf{r},0)$, where
$\rho_{0}(\mathbf{r},0)$ is the desired initial mean-field density profile,
and $\bm{v}_{0}(\mathbf{r},0)$ and $S_{0}(\mathbf{r},0)$ the desired initial
mean-field velocity and phase profiles respectively. For the scenarios
presented in Section V.1 these profiles correspond to
$\rho_{0}(x,0)=G^{(1)}(x,x,0)=|\Psi_{0}(x,0)|^{2}=N_{\text{bg}}[1+\beta
e^{-x^{2}/2\sigma^{2}}]^{2}/L$, whereas $v_{0}(x,0)=S_{0}(x,0)=0$.

We initialize the linearized positive-P approaches in the same way as the full
positive-P approach; that is, by initializing everything in the mean-field
variables as just discussed, and setting
$\delta\rho(\mathbf{r},0)=\delta\bm{v}(\mathbf{r},0)=\delta S(\mathbf{r},0)=0$
(or equivalently,
$\delta\psi(\mathbf{r},0)=\delta\tilde{\psi}(\mathbf{r},0)=0$), i.e., assuming
that the fluctuating components of the field operators are initially in a
vacuum state.

##  Appendix C Matrix product state implementation

To implement matrix product state approaches, the Lieb-Liniger Hamiltonian
[33]

$\displaystyle\hat{H}_{\text{LL}}$ $\displaystyle=\int
dx~{}\hat{\Psi}^{\dagger}\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x^{2}}+V(x)\right)\hat{\Psi}$
$\displaystyle\qquad\qquad\qquad+\frac{g_{\text{1D}}}{2}\int
dx~{}\hat{\Psi}^{\dagger}\hat{\Psi}^{\dagger}\hat{\Psi}\hat{\Psi}$ (78)

is approximated as a Bose-Hubbard model by discretizing in real-space.
Infinite matrix product state (iMPS) methods [91, 92, 93] can then be employed
to determine the ground and coherent states of the Bose-Hubbard model and to
simulate the ensuing dynamics after performing a trap quench.

The real-space discretization replaces the field operator $\hat{\Psi}(x_{j})$
by $\hat{b}_{j}/\sqrt{\Delta x}$. Here, $\hat{b}_{j}$ is the bosonic creation
operator at site $j$ ($j=1,2,\ldots,M,$ for $M=L/\Delta x$ lattice sites),
$\Delta x$ is the lattice spacing and $x_{j}=j\Delta x$ [94, 95, 96]. This
discretization results in the Bose-Hubbard Hamiltonian

$\displaystyle\hat{H}_{\text{BH}}=-J\sum_{j}$
$\displaystyle\left(\hat{b}_{j}^{\dagger}\hat{b}_{j+1}+\hat{b}_{j+1}^{\dagger}\hat{b}_{j}\right)$
$\displaystyle+\frac{U}{2}\sum_{j}\hat{b}_{j}^{\dagger
2}\hat{b}_{j}^{2}+\sum_{j}V(x_{j})\hat{b}_{j}^{\dagger}\hat{b}_{j},$ (79)

where $J=\hbar^{2}/2m\Delta x^{2}$ and $U=g_{\text{1D}}/\Delta x$. This
approximation is valid in the continuum limit $\Delta x\rightarrow 0$ [96, 97,
98] and small average lattice occupancy,
$\langle\hat{b}_{j}^{\dagger}\hat{b}_{j}\rangle\ll 1$.

The $t=0$ ground and coherent states are prepared using the infinite Density
Matrix Renormalization Group (iDMRG) algorithm, with $M=1000-4000$ sites per
unit cell. For the ground state, the iMPS is optimized for the Hamiltonian
(79) using a bond dimension (basis states kept) of $m=70-80$. This results in
an energy variance $\sigma_{E}^{2}/J^{2}$ of $\sim 10^{-9}$ per site. On the
other hand, the coherent state, which is a classical product state, is
prepared with $m=1$ and $U=0$. This produces $\sigma_{E}^{2}/J^{2}\sim
10^{-13}$ per site. In each case the trapping potential used to prepare the
initial state is given by,

$\displaystyle V(x_{j})$
$\displaystyle=\frac{\hbar^{2}\beta}{2m\sigma^{4}}\left(x_{j}^{2}-\sigma^{2}\right)\left[e^{x_{j}^{2}/2\sigma^{2}}+\beta\right]^{-1}$
$\displaystyle\qquad\qquad\qquad-\frac{g_{\text{1D}}N_{\text{bg}}}{L}\left[1+\beta
e^{-x_{j}^{2}/2\sigma^{2}}\right]^{2}$ (80)

which produces exactly the desired density
$\rho_{0}(x_{j},0)=N_{\text{bg}}[1+\beta e^{-x_{j}^{2}/2\sigma^{2}}]^{2}/L$ in
the mean-field approximation, i.e., as the ground state solution to the GPE
(22).

The time-evolution of the ground and coherent states are executed using the
infinite time-evolving block decimation (iTEBD) algorithm [99] with the
optimized 4th-order decomposition [100]. Here, the quench is performed at
$t=0$ by setting $V(x_{j})=0$ and evolving both states with the Hamiltonian
(79). For the coherent state, the value of $U$ during the dynamics is set to
that of the ground state. A time step of 0.2 (in dimensionless units of time
$Jt/\hbar$), and a cutoff density matrix eigenvalue of $10^{-10}-10^{-9}$,
corresponding to a cutoff singular value of $\sim 10^{-5}$, is used throughout
this work. This results in a distance $d\equiv 1-|F|$ with the $t=0$ state,
after backwards evolution, of $d\sim 10^{-7}-10^{-6}$ per site, where $F$ is
the fidelity with the initial state.

##  Appendix D Observables in the stochastic Bogoliubov formalism

In the stochastic Bogoliubov approach, the reduced one-body density matrix is
given by

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle=\langle[\Psi_{0}^{*}(\mathbf{r},t)+\delta\tilde{\psi}(\mathbf{r},t)]$
$\displaystyle\qquad\qquad\times[\Psi_{0}(\mathbf{r}^{\prime},t)+\delta\psi(\mathbf{r}^{\prime},t)]\rangle$
$\displaystyle=\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle,$
(81)

to second order in fluctuating components.

Hence, the particle number density is simply

$\displaystyle G^{(1)}(\mathbf{r},\mathbf{r},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}(\mathbf{r},t)\rangle$
$\displaystyle=\left|\Psi_{0}(\mathbf{r},t)\right|^{2}+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle.$
(82)

Similarly to the LSQHD equations, the stochastic Bogoliubov equations are
initialized in a coherent-state and so the average of the fluctuating fields
will remain
$\langle\delta\psi(\mathbf{r},t)\rangle=\langle\delta\tilde{\psi}(\mathbf{r},t)\rangle=0$.

Normalizing the reduced one-body density matrix and maintaining only second
order terms leads to

$\displaystyle g^{(1)}$ $\displaystyle(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\frac{G^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)}{\sqrt{G^{(1)}(\mathbf{r},\mathbf{r},t)}\sqrt{G^{(1)}(\mathbf{r}^{\prime},\mathbf{r}^{\prime},t)}}$
$\displaystyle\approx\sqrt{\frac{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}{\Psi_{0}(\mathbf{r},t)\Psi_{0}^{*}(\mathbf{r}^{\prime},t)}}\left(1+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}\right.$
$\displaystyle\qquad-\left.\frac{1}{2}\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r},t)}-\frac{1}{2}\frac{\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\Psi_{0}(\mathbf{r}^{\prime},t)}\right).$
(83)

The density-density correlation function in normally ordered form is given by

$\displaystyle G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\langle\hat{\Psi}^{\dagger}(\mathbf{r},t)\hat{\Psi}^{\dagger}(\mathbf{r}^{\prime},t)\hat{\Psi}(\mathbf{r}^{\prime},t)\hat{\Psi}(\mathbf{r},t)\rangle$
$\displaystyle\approx|\Psi_{0}(\mathbf{r},t)|^{2}|\Psi_{0}(\mathbf{r}^{\prime},t)|^{2}$
$\displaystyle\quad+|\Psi_{0}(\mathbf{r},t)|^{2}\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle\quad+|\Psi_{0}(\mathbf{r}^{\prime},t)|^{2}\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle$
$\displaystyle\quad+\Psi_{0}(\mathbf{r},t)\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle\quad+\Psi_{0}(\mathbf{r}^{\prime},t)\Psi_{0}^{*}(\mathbf{r},t)\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r},t)\rangle$
$\displaystyle\quad+\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\langle\delta\psi(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle$
$\displaystyle\quad+\Psi_{0}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\rangle$
(84)

up to second order. Normalizing gives

$\displaystyle g^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)$
$\displaystyle\equiv\frac{G^{(2)}(\mathbf{r},\mathbf{r}^{\prime},t)}{G^{(1)}(\mathbf{r},\mathbf{r},t)~{}G^{(1)}(\mathbf{r}^{\prime},\mathbf{r}^{\prime},t)}$
$\displaystyle\approx
1+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}+\frac{\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r},t)\rangle}{\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\Psi_{0}(\mathbf{r},t)}$
$\displaystyle\quad+\frac{\langle\delta\psi(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}^{*}(\mathbf{r}^{\prime},t)}$
$\displaystyle=1+2\Re{\frac{\langle\delta\psi(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}}$
$\displaystyle\quad+\frac{\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r}^{\prime},t)\rangle}{\Psi_{0}^{*}(\mathbf{r},t)\Psi_{0}(\mathbf{r}^{\prime},t)}+\frac{\langle\delta\tilde{\psi}(\mathbf{r}^{\prime},t)\delta\psi(\mathbf{r},t)\rangle}{\Psi_{0}^{*}(\mathbf{r}^{\prime},t)\Psi_{0}(\mathbf{r},t)}$
(85)

to second order also.

##  Appendix E Comparisons between an initial coherent state and the full
many-body ground state

Here we make comparisons between the results presented in Section V.1.1 of the
main text (where each of the approaches are initialized in a coherent state)
and those obtained from a matrix product state (MPS) calculation initialized
in the full many-body ground state of a trapping potential which reproduces
$\rho(x,0)=N_{\text{bg}}[1+\beta e^{-x^{2}/2\sigma^{2}}]^{2}/L$ in the mean-
field density, Eq. (80). From here on we refer to this ground-state MPS
calculation as GS-MPS.

###  E.1 Density $\rho(x,t)$

On a global scale, the densities predicted by each method from Sec. V.1.1
agree well with the GS-MPS calculation. We do point out however, that the
effect of quantum fluctuations can still be discerned by examining the peaks
and troughs of the oscillatory wave-train in the time-evolved density
$\rho(x,t)$. In Figure 9 we show the details of the first density trough. Here
we see that the full positive-P, truncated Wigner, and CS-MPS results agree
with each other within one standard error and that the quantum fluctuations
present in these approaches slightly diminish the oscillation amplitude as
compared with LSQHD (and stochastic Bogoliubov which is not plotted here),
whose density is equivalent to the mean-field density at first order in
fluctuations. Moreover, the oscillation amplitude is diminished further for
the GS-MPS result, where the initial state includes quantum fluctuations
already and the initial density profile departs slightly from the coherent
state $\rho(x,0)$ since it is found as a many-body ground state of the
trapping potential, Eq. (80), which only produces
$\rho(x,0)=N_{\text{bg}}[1+\beta e^{-x^{2}/2\sigma^{2}}]^{2}/L$ in the mean-
field. This difference in initial density profiles is shown in Fig. 10.

Figure 9: Predictions of the real-space density $\rho(x,t)L$ from different
approaches at dimensionless times $\tau=0$ and $\tau=0.001$. The inset shows
the details of the first density trough. The shaded regions around the
stochastic results denote one standard error of uncertainty, which for LSQHD
is on the order of the linewidth used here. Figure 10: Comparison of the
initial real space density $\rho(x,0)L=N_{\text{bg}}[1+\beta
e^{-x^{2}/2\sigma^{2}}]^{2}$ with the GS-MPS result obtained by determining
the ground state of the system with Eq. (80) as the trapping potential.

###  E.2 Correlation functions $g^{(2)}(x,x,t)$ and $g^{(2)}(x,-x,t)$ from the
MPS approaches

In figures 11 and 12 we compare the differences in predictions of
$g^{(2)}(x,x,t)$ and $g^{(2)}(x,-x,t)$, respectively, between an initial
coherent state and the full many-body ground state. For clarity we compare
only the CS-MPS result with the GS-MPS calculation and we recall the good
agreement of each stochastic approach with the CS-MPS result, especially for
LSQHD.

We first point out the good qualitative agreement between the two different
initializations; the coherent-state simulations correctly capture all of the
broad features present in the exact ground-state calculation for both
$g^{(2)}(x,x,t)$ and $g^{(2)}(x,-x,t)$.

The primary quantitative disagreement for $g^{(2)}(x,x,t)$ is the shift in
background value due to the initial correlations present in the GS-MPS result.
Additionally, the oscillation amplitude is reduced in the CS-MPS result
compared with GS-MPS. For $g^{(2)}(x,-x,t)$ the primary differences between
the coherent state and ground state initializations are the opposite
orientation of the oscillations, and this time their larger amplitude in the
CS-MPS result compared to GS-MPS.

Figure 11: The same-point correlation function $g^{(2)}(x,x,t)$ of Fig. 3, as
predicted by CS-MPS and GS-MPS at dimensionless times $\tau=0$ and
$\tau=0.001$. Figure 12: The opposite-point correlation function
$g^{(2)}(x,-x,t)$ of Fig. 4, as predicted by CS-MPS and GS-MPS at
dimensionless times $\tau=0$ and $\tau=0.001$.

##  Appendix F Comparison of the density and correlation function
$g^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$ between the LSQHD and stochastic
Bogoliubov approaches

###  F.1 Density comparisons

A comparison of the density predictions of the LSQHD and stochastic Bogoliubov
approaches for the two examples given in the main text are provided in Figures
13 and 14. The LSQHD prediction is given by equation (63), which differs from
the stochastic Bogoliubov expression (68) by the term responsible for particle
number growth,
$\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle$. This
particle number growth can be observed in the figures as the discrepancy
between the respective curves at dimensionless time $\tau=0.001$.

Figure 13: The real-space density $\rho(x,t)$ for the scenario presented in
Fig. 1; here, we compare the results from the LSQHD and stochastic Bogoliubov
(Stoch. Bogo.) approaches. Figure 14: The real-space density $\rho(x,t)$ for
the scenario presented in Fig. 5; here we compare the results from the LSQHD
and stochastic Bogoliubov (Stoch. Bogo.) approaches.

###  F.2 Correlation function $g^{(1)}(\mathbf{r},\mathbf{r}^{\prime},t)$

To identify the difference between the LSQHD expression for
$g^{(1)}(\mathbf{r},\mathbf{r}^{\prime})$ (64) and the stochastic Bogoliubov
expression (69), one can apply Madelung’s transformation
$\Psi_{0}(\mathbf{r},t)=\sqrt{\rho_{0}(\mathbf{r},t)}e^{iS_{0}(\mathbf{r},t)}$
and substitute (39) into the stochastic Bogoliubov result (69). Doing so leads
to

$\displaystyle g$ ${}^{(1)}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle\approx
e^{-i\left[S_{0}(\mathbf{r})-S_{0}(\mathbf{r}^{\prime})\right]}$
$\displaystyle\quad\times\left(1+\expectationvalue{\delta S(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}-\frac{1}{2}\left[\expectationvalue{\delta
S^{2}(\mathbf{r})}+\expectationvalue{\delta
S^{2}(\mathbf{r}^{\prime})}\right]\right.$
$\displaystyle\qquad\quad+\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{4\rho_{0}(\mathbf{r})\rho_{0}(\mathbf{r}^{\prime})}-\frac{1}{8}\left[\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r})}}{\rho_{0}^{2}(\mathbf{r})}+\frac{\expectationvalue{\delta\rho^{2}(\mathbf{r}^{\prime})}}{\rho_{0}^{2}(\mathbf{r}^{\prime})}\right]$
$\displaystyle\qquad\qquad\quad+\left.\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}\right]\right),$
(86)

which differs from the LSQHD expression (64) in that the term

$\displaystyle\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r}^{\prime})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r})}}{\rho_{0}(\mathbf{r})}\right]$ (87)

is absent here in Eq. (86).

We explore the effect of this term by numerically computing
$g^{(1)}(x,x^{\prime},t)$ in the LSQHD, stochastic Bogoliubov, and CS-MPS
approaches for the scenario of Fig. 1 given in Sec. V.1.1 of the main text. A
comparison of the slices $g^{(1)}(x,0,t)$ and $g^{(1)}(x,-x,t)$ are shown in
Figs. 15 and 16, respectively. The differences between the predictions of the
LSQHD and stochastic Bogoliubov approaches is on the order of the uncertainty
in these calculations, and hence we conclude that the term (87) responsible
for the difference in expressions (86) and (64) is negligible here.

Figure 15: The correlation function $g^{(1)}(x,0,t)$ for the scenario
presented in Figure 1 of the main text. Here however, the stochastic results
are an average over only 500,000 trajectories and the shaded regions denote
one standard error of uncertainty. Figure 16: The correlation function
$g^{(1)}(x,-x,t)$ for the scenario presented in Figure 1 of the main text.
Here however, the stochastic results are an average over only 500,000
trajectories and the shaded regions denote one standard error of uncertainty.

Further exploration reveals that the other term containing cross products of
$\delta\rho$ and $\delta S$,

$\displaystyle\frac{i}{2}\left[\frac{\expectationvalue{\delta\rho(\mathbf{r})\delta
S(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r})}-\frac{\expectationvalue{\delta
S(\mathbf{r})\delta\rho(\mathbf{r}^{\prime})}}{\rho_{0}(\mathbf{r}^{\prime})}\right],$
(88)

also contributes negligibly to $g^{(1)}(x,x^{\prime},t)$ in this scenario.

Additionally, the LSQHD and stochastic Bogoliubov results lie close to the
numerically exact CS-MPS prediction, confirming the validity and accuracy of
these linearized schemes in computing $g^{(1)}(x,x^{\prime},t)$ for this
scenario.

## References

  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in Out-of-Equilibrium $XXZ$ Chains: Exact Profiles of Charges and Currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Doyon  _et al._ [2017] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura, Large-Scale Description of Interacting One-Dimensional Bose Gases: Generalized Hydrodynamics Supersedes Conventional Hydrodynamics, Phys. Rev. Lett. 119, 195301 (2017).
  * Schemmer  _et al._ [2019] M. Schemmer, I. Bouchoule, B. Doyon, and J. Dubail, Generalized Hydrodynamics on an Atom Chip, Phys. Rev. Lett. 122, 090601 (2019).
  * Malvania  _et al._ [2021] N. Malvania, Y. Zhang, Y. Le, J. Dubail, M. Rigol, and D. S. Weiss, Generalized hydrodynamics in strongly interacting 1D Bose gases, Science 373, 1129 (2021).
  * Bouchoule and Dubail [2022] I. Bouchoule and J. Dubail, Generalized hydrodynamics in the one-dimensional Bose gas: theory and experiments, Journal of Statistical Mechanics: Theory and Experiment , 014003 (2022).
  * Bouchoule  _et al._ [2016] I. Bouchoule, S. S. Szigeti, M. J. Davis, and K. V. Kheruntsyan, Finite-temperature hydrodynamics for one-dimensional Bose gases: Breathing-mode oscillations as a case study, Phys. Rev. A 94, 051602 (2016).
  * Atas  _et al._ [2019] Y. Y. Atas, S. A. Simmons, and K. V. Kheruntsyan, Finite-temperature dynamics of a Tonks-Girardeau gas in a frequency-modulated harmonic trap, Phys. Rev. A 100, 043602 (2019).
  * Euler [1757] L. Euler, Principes generaux du mouvement des fluides, Mémoires de l’académie des sciences de Berlin 11, 274 (1757).
  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz,  _Fluid Mechanics: Course of Theoretical Physics_ , 2nd ed. (Pergamon Press, 1987) p. 539.
  * Minguzzi  _et al._ [2001] A. Minguzzi, P. Vignolo, M. L. Chiofalo, and M. P. Tosi, Hydrodynamic excitations in a spin-polarized Fermi gas under harmonic confinement in one dimension, Phys. Rev. A 64, 033605 (2001).
  * Menotti and Stringari [2002] C. Menotti and S. Stringari, Collective oscillations of a one-dimensional trapped Bose-Einstein gas, Phys. Rev. A 66, 043610 (2002).
  * Peotta and Ventra [2014] S. Peotta and M. D. Ventra, Quantum shock waves and population inversion in collisions of ultracold atomic clouds, Phys. Rev. A 89, 013621 (2014).
  * Kulkarni and Abanov [2012] M. Kulkarni and A. G. Abanov, Hydrodynamics of cold atomic gases in the limit of weak nonlinearity, dispersion, and dissipation, Phys. Rev. A 86, 033614 (2012).
  * Pitaevskii and Stringari [2003] L. P. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation_ (Clarendon Press, 2003) p. 400.
  * Pethick and Smith [2008] C. J. Pethick and H. Smith,  _Bose-Einstein Condensation in Dilute Gases_ , 2nd ed. (Cambridge University Press, 2008) p. 400.
  * Madelung [1927] E. Madelung, Quantentheorie in hydrodynamischer Form, Zeitschrift für Physik 40, 322 (1927).
  * Griffin and Nikuni [2000] A. Griffin and T. Nikuni, Two-Fluid Hydrodynamics in Trapped Bose Gases and in Superfluid Helium, Journal of Low Temperature Physics 121, 247 (2000).
  * Tisza [1938] L. Tisza, Transport Phenomena in Helium II, Nature 141, 913 (1938).
  * Landau [1941] L. Landau, Theory of the Superfluidity of Helium II, Phys. Rev. 60, 356 (1941).
  * Doyon and Yoshimura [2017] B. Doyon and T. Yoshimura, A note on generalized hydrodynamics: inhomogeneous fields and other concepts, SciPost Phys. 2, 014 (2017).
  * Ruggiero  _et al._ [2020] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Quantum Generalized Hydrodynamics, Phys. Rev. Lett. 124, 140603 (2020).
  * Ruggiero  _et al._ [2021] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, Quantum generalized hydrodynamics of the Tonks–Girardeau gas: density fluctuations and entanglement entropy, Journal of Physics A: Mathematical and Theoretical 55, 024003 (2021).
  * El and Hoefer [2016] G. El and M. Hoefer, Dispersive shock waves and modulation theory, Physica D: Nonlinear Phenomena 333, 11 (2016).
  * Dutton  _et al._ [2001] Z. Dutton, M. Budde, C. Slowe, and L. V. Hau, Observation of Quantum Shock Waves Created with Ultra- Compressed Slow Light Pulses in a Bose-Einstein Condensate, Science 293, 663 (2001).
  * Damski [2004] B. Damski, Formation of shock waves in a Bose-Einstein condensate, Phys. Rev. A 69, 043610 (2004).
  * Simula  _et al._ [2005] T. P. Simula, P. Engels, I. Coddington, V. Schweikhard, E. A. Cornell, and R. J. Ballagh, Observations on Sound Propagation in Rapidly Rotating Bose-Einstein Condensates, Phys. Rev. Lett. 94, 080404 (2005).
  * Hoefer  _et al._ [2006] M. A. Hoefer, M. J. Ablowitz, I. Coddington, E. A. Cornell, P. Engels, and V. Schweikhard, Dispersive and classical shock waves in Bose-Einstein condensates and gas dynamics, Phys. Rev. A 74, 023623 (2006).
  * Chang  _et al._ [2008] J. J. Chang, P. Engels, and M. A. Hoefer, Formation of Dispersive Shock Waves by Merging and Splitting Bose-Einstein Condensates, Phys. Rev. Lett. 101, 170404 (2008).
  * Meppelink  _et al._ [2009] R. Meppelink, S. B. Koller, J. M. Vogels, P. van der Straten, E. D. van Ooijen, N. R. Heckenberg, H. Rubinsztein-Dunlop, S. A. Haine, and M. J. Davis, Observation of shock waves in a large Bose-Einstein condensate, Phys. Rev. A 80, 043606 (2009).
  * Bulgac  _et al._ [2012] A. Bulgac, Y.-L. Luo, and K. J. Roche, Quantum Shock Waves and Domain Walls in the Real-Time Dynamics of a Superfluid Unitary Fermi Gas, Phys. Rev. Lett. 108, 150401 (2012).
  * Simmons  _et al._ [2020] S. A. Simmons, F. A. Bayocboc, J. C. Pillay, D. Colas, I. P. McCulloch, and K. V. Kheruntsyan, What is a Quantum Shock Wave?, Phys. Rev. Lett. 125, 180401 (2020).
  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Exact Analysis of an Interacting Bose Gas. I. The General Solution and the Ground State, Phys. Rev. 130, 1605 (1963).
  * Olshanii [1998] M. Olshanii, Atomic Scattering in the Presence of an External Confinement and a Gas of Impenetrable Bosons, Phys. Rev. Lett. 81, 938 (1998).
  * Petrov  _et al._ [2004] D. Petrov, D. M. Gangardt, and G. V. Shlyapnikov, Low-dimensional trapped gases, in  _Journal de Physique IV (Proceedings)_, Vol. 116 (EDP sciences, 2004) pp. 5–44.
  * Abrikosov  _et al._ [1965] A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzialoshinskiĭ,  _Quantum field theoretical methods in statistical physics_ , Vol. 4 (Pergamon, 1965).
  * Morgan [2000] S. A. Morgan, A gapless theory of Bose-Einstein condensation in dilute gases at finite temperature, Journal of Physics B: Atomic, Molecular and Optical Physics 33, 3847 (2000).
  * Drummond and Gardiner [1980] P. D. Drummond and C. W. Gardiner, Generalised P-representations in quantum optics, Journal of Physics A: Mathematical and General 13, 2353 (1980).
  * Gardiner and Crispin [1985] C. W. Gardiner and W. Crispin,  _Handbook of stochastic methods for physics, chemistry, and the natural sciences_ , 2nd ed., Springer series in synergetics; v. 13 (Springer-Verlag, Berlin; New York, 1985).
  * Steel  _et al._ [1998] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collett, D. F. Walls, and R. Graham, Dynamical quantum noise in trapped Bose-Einstein condensates, Phys. Rev. A 58, 4824 (1998).
  * Hope and Olsen [2001] J. J. Hope and M. K. Olsen, Quantum Superchemistry: Dynamical Quantum Effects in Coupled Atomic and Molecular Bose-Einstein Condensates, Phys. Rev. Lett. 86, 3220 (2001).
  * Deuar and Drummond [2006] P. Deuar and P. D. Drummond, First-principles quantum dynamics in interacting Bose gases: I. The positive P representation, Journal of Physics A: Mathematical and General 39, 1163 (2006).
  * Deuar and Drummond [2007] P. Deuar and P. D. Drummond, Correlations in a BEC Collision: First-Principles Quantum Dynamics with 150 000 Atoms, Phys. Rev. Lett. 98, 120402 (2007).
  * Perrin  _et al._ [2008] A. Perrin, C. M. Savage, D. Boiron, V. Krachmalnicoff, C. I. Westbrook, and K. V. Kheruntsyan, Atomic four-wave mixing via condensate collisions, New Journal of Physics 10, 045021 (2008).
  * Krachmalnicoff  _et al._ [2010] V. Krachmalnicoff, J.-C. Jaskula, M. Bonneau, V. Leung, G. B. Partridge, D. Boiron, C. I. Westbrook, P. Deuar, P. Ziń, M. Trippenbach, and K. V. Kheruntsyan, Spontaneous Four-Wave Mixing of de Broglie Waves: Beyond Optics, Phys. Rev. Lett. 104, 150402 (2010).
  * Gilchrist  _et al._ [1997] A. Gilchrist, C. W. Gardiner, and P. D. Drummond, Positive P representation: Application and validity, Phys. Rev. A 55, 3014 (1997).
  * Castin [2001] Y. Castin,  _Coherent atomic matter waves_, edited by R. Kaiser, C. Westbrook, and F. David (Springer-Verlag, 2001) pp. 1–136.
  * [48] This is similar to Glauber’s definition of a coherent state of an optical field as a factorization property of all higher-order correlation functions, equivalent to its definition as an eigenstate of the field annihilation operator [101, 102, 103].
  * Cazalilla  _et al._ [2011] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, One dimensional bosons: From condensed matter systems to ultracold gases, Rev. Mod. Phys. 83, 1405 (2011).
  * [50] Whilst the stochastic Bogoliubov equations are often written in Îto form, we remind the reader that the Îto and Stratonovich forms of equations (2) and (3) are equivalent (up to a physically irrelevant constant phase shift). Hence, the Îto and Stratonovich forms of linearized schemes based on these equations are also equivalent.
  * Jaskula  _et al._ [2010] J.-C. Jaskula, M. Bonneau, G. B. Partridge, V. Krachmalnicoff, P. Deuar, K. V. Kheruntsyan, A. Aspect, D. Boiron, and C. I. Westbrook, Sub-Poissonian Number Differences in Four-Wave Mixing of Matter Waves, Phys. Rev. Lett. 105, 190402 (2010).
  * Deuar  _et al._ [2011] P. Deuar, J. Chwedeńczuk, M. Trippenbach, and P. Ziń, Bogoliubov dynamics of condensate collisions using the positive-$P$ representation, Phys. Rev. A 83, 063625 (2011).
  * Kheruntsyan  _et al._ [2012] K. V. Kheruntsyan, J.-C. Jaskula, P. Deuar, M. Bonneau, G. B. Partridge, J. Ruaudel, R. Lopes, D. Boiron, and C. I. Westbrook, Violation of the Cauchy-Schwarz Inequality with Matter Waves, Phys. Rev. Lett. 108, 260401 (2012).
  * Deuar  _et al._ [2013] P. Deuar, T. Wasak, P. Ziń, J. Chwedeńczuk, and M. Trippenbach, Tradeoffs for number squeezing in collisions of Bose-Einstein condensates, Phys. Rev. A 88, 013617 (2013).
  * Deuar  _et al._ [2014] P. Deuar, J.-C. Jaskula, M. Bonneau, V. Krachmalnicoff, D. Boiron, C. I. Westbrook, and K. V. Kheruntsyan, Anisotropy in $s$-wave Bose-Einstein condensate collisions and its relationship to superradiance, Phys. Rev. A 90, 033613 (2014).
  * Lewis-Swan and Kheruntsyan [2014] R. Lewis-Swan and K. Kheruntsyan, Proposal for demonstrating the Hong–Ou–Mandel effect with matter waves, Nature Communications 5, 1 (2014).
  * Ross  _et al._ [2021] J. A. Ross, P. Deuar, D. K. Shin, K. F. Thomas, B. M. Henson, S. S. Hodgman, and A. G. Truscott, Survival of the quantum depletion of a condensate after release from a harmonic trap in theory and experiment, arXiv:2103.15283 [cond-mat.quant-gas] (2021).
  * Bogoliubov [1947] N. N. Bogoliubov, On the theory of superfluidity, J. Phys. (USSR) 11, 23 (1947).
  * Noziéres and Pines [1990] P. Noziéres and D. Pines,  _The Theory of Quantum Liquids_ , Vol. II (Addison-Wesley, 1990).
  * Fetter [1972] A. L. Fetter, Nonuniform states of an imperfect bose gas, Annals of Physics 70, 67 (1972).
  * Gennes [1999] P.-G. d. Gennes,  _Superconductivity of metals and alloys_ , Advanced book classics (Advanced Book Program, Perseus Books, 1999).
  * Dziarmaga and Sacha [2005] J. Dziarmaga and K. Sacha, Images of a Bose–Einstein condensate: diagonal dynamical Bogoliubov vacuum, Journal of Physics B: Atomic, Molecular and Optical Physics 39, 57 (2005).
  * Olsen and Bradley [2009] M. Olsen and A. Bradley, Numerical representation of quantum states in the positive-P and Wigner representations, Optics Communications 282, 3924 (2009).
  * Lewis-Swan and Kheruntsyan [2020] R. J. Lewis-Swan and K. V. Kheruntsyan, Atomic twin beams and violation of a motional-state Bell inequality from a phase-fluctuating quasicondensate source, Phys. Rev. A 101, 043615 (2020).
  * Gardiner [1997] C. W. Gardiner, Particle-number-conserving Bogoliubov method which demonstrates the validity of the time-dependent Gross-Pitaevskii equation for a highly condensed Bose gas, Phys. Rev. A 56, 1414 (1997).
  * Castin and Dum [1998] Y. Castin and R. Dum, Low-temperature Bose-Einstein condensates in time-dependent traps: Beyond the $U(1)$ symmetry-breaking approach, Phys. Rev. A 57, 3008 (1998).
  * Walls and Milburn [2007] D. F. Walls and G. J. Milburn,  _Quantum optics_ (Springer Science & Business Media, 2007).
  * Ögren and Kheruntsyan [2009] M. Ögren and K. V. Kheruntsyan, Atom-atom correlations in colliding Bose-Einstein condensates, Phys. Rev. A 79, 021606 (2009).
  * Kheruntsyan [2006] K. V. Kheruntsyan, Quantum Atom Optics with Fermions from Molecular Dissociation, Phys. Rev. Lett. 96, 110401 (2006).
  * Cazalilla [2004] M. A. Cazalilla, Bosonizing one-dimensional cold atomic gases, Journal of Physics B: Atomic, Molecular and Optical Physics 37, S1 (2004).
  * Mora and Castin [2003] C. Mora and Y. Castin, Extension of Bogoliubov theory to quasicondensates, Phys. Rev. A 67, 053615 (2003).
  * Kheruntsyan  _et al._ [2005] K. V. Kheruntsyan, D. M. Gangardt, P. D. Drummond, and G. V. Shlyapnikov, Finite-temperature correlations and density profiles of an inhomogeneous interacting one-dimensional Bose gas, Phys. Rev. A 71, 053615 (2005).
  * [73] I. P. McCulloch, The Matrix Product Toolkit, https://people.smp.uq.edu.au/IanMcCulloch/mptoolkit.
  * Sinatra  _et al._ [2002] A. Sinatra, C. Lobo, and Y. Castin, The truncated Wigner method for Bose-condensed gases: limits of validity and applications, Journal of Physics B: Atomic, Molecular and Optical Physics 35, 3599 (2002).
  * Isella and Ruostekoski [2005] L. Isella and J. Ruostekoski, Nonadiabatic dynamics of a Bose-Einstein condensate in an optical lattice, Phys. Rev. A 72, 011601 (2005).
  * Isella and Ruostekoski [2006] L. Isella and J. Ruostekoski, Quantum dynamics in splitting a harmonically trapped Bose-Einstein condensate by an optical lattice: Truncated Wigner approximation, Phys. Rev. A 74, 063625 (2006).
  * Martin and Ruostekoski [2012] A. D. Martin and J. Ruostekoski, Truncated Wigner Method for Bose Gases, arXiv:1009.1073v2 [cond-mat.quant-gas] (2012).
  * Dennis  _et al._ [2013] G. R. Dennis, J. J. Hope, and M. T. Johnsson, XMDS2: Fast, scalable simulation of coupled stochastic partial differential equations, Computer Physics Communications 184, 201 (2013).
  * Drummond [1983] P. Drummond, Central partial difference propagation algorithms, Computer Physics Communications 29, 211 (1983).
  * Werner and Drummond [1997] M. Werner and P. Drummond, Robust Algorithms for Solving Stochastic Partial Differential Equations, Journal of Computational Physics 132, 312 (1997).
  * Kamchatnov  _et al._ [2002] A. M. Kamchatnov, R. A. Kraenkel, and B. A. Umarov, Asymptotic soliton train solutions of the defocusing nonlinear Schrödinger equation, Phys. Rev. E 66, 036609 (2002).
  * Damski [2006] B. Damski, Shock waves in a one-dimensional Bose gas: From a Bose-Einstein condensate to a Tonks gas, Phys. Rev. A 73, 043601 (2006).
  * Hoefer  _et al._ [2009] M. Hoefer, P. Engels, and J. Chang, Matter–wave interference in Bose–Einstein condensates: A dispersive hydrodynamic perspective, Physica D: Nonlinear Phenomena 238, 1311 (2009), Nonlinear Phenomena in Degenerate Quantum Gases.
  * Grozdanov  _et al._ [2019] S. Grozdanov, P. K. Kovtun, A. O. Starinets, and P. Tadić, Convergence of the Gradient Expansion in Hydrodynamics, Phys. Rev. Lett. 122, 251601 (2019).
  * Imambekov and Glazman [2009] A. Imambekov and L. I. Glazman, Universal Theory of Nonlinear Luttinger Liquids, Science 323, 228 (2009).
  * Imambekov  _et al._ [2012] A. Imambekov, T. L. Schmidt, and L. I. Glazman, One-dimensional quantum liquids: Beyond the Luttinger liquid paradigm, Rev. Mod. Phys. 84, 1253 (2012).
  * [87] H. Spohn, Fluctuating hydrodynamics approach to equilibrium time correlations for anharmonic chains, in _Thermal Transport in Low Dimensions_ ; Ed. S. Lepri (Berlin, Springer, 2016).
  * Bulchandani  _et al._ [2021] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, Superdiffusion in spin chains, Journal of Statistical Mechanics: Theory and Experiment 2021, 084001 (2021).
  * Blakie  _et al._ [2008] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner, Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, Advances in Physics 57, 363 (2008).
  * Ng  _et al._ [2019] K. L. Ng, B. Opanchuk, M. D. Reid, and P. D. Drummond, Nonlocal Pair Correlations in a Higher-Order Bose Gas Soliton, Phys. Rev. Lett. 122, 203604 (2019).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * McCulloch [2007] I. P. McCulloch, From density-matrix renormalization group to matrix product states, Journal of Statistical Mechanics: Theory and Experiment 2007, P10014 (2007).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Annals of Physics 326, 96 (2011), January 2011 Special Issue.
  * Cazalilla [2003] M. A. Cazalilla, One-dimensional optical lattices and impenetrable bosons, Phys. Rev. A 67, 053606 (2003).
  * Schmidt  _et al._ [2005] B. Schmidt, L. I. Plimak, and M. Fleischhauer, Stochastic simulation of a finite-temperature one-dimensional Bose gas: From the Bogoliubov to the Tonks-Girardeau regime, Phys. Rev. A 71, 041601 (2005).
  * Schmidt and Fleischhauer [2007] B. Schmidt and M. Fleischhauer, Exact numerical simulations of a one-dimensional trapped Bose gas, Phys. Rev. A 75, 021601 (2007).
  * Muth  _et al._ [2010a] D. Muth, B. Schmidt, and M. Fleischhauer, Fermionization dynamics of a strongly interacting one-dimensional Bose gas after an interaction quench, New Journal of Physics 12, 083065 (2010a).
  * Muth  _et al._ [2010b] D. Muth, M. Fleischhauer, and B. Schmidt, Discretized versus continuous models of $p$-wave interacting fermions in one dimension, Phys. Rev. A 82, 013602 (2010b).
  * Vidal [2007] G. Vidal, Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension, Phys. Rev. Lett. 98, 070201 (2007).
  * Barthel and Zhang [2020] T. Barthel and Y. Zhang, Optimized Lie-Trotter-Suzuki decompositions for two and three non-commuting terms, Annals of Physics 418, 168165 (2020).
  * Glauber [2007] R. J. Glauber,  _Quantum theory of optical coherence: selected papers and lectures_ (John Wiley & Sons, 2007).