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 behavior of quantum many-body systems is a notoriously challenging task. A particularly powerful method capable of describing such non-equilibrium behavior 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 those 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 linearized 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}\hat{\Psi}^{\dagger}\laplacian\hat{\Psi}+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 use of the contact interaction potential $g\delta({\mathbf{r}}-{\mathbf{r}}^{\prime})$ in the Hamiltonian (1), instead of the true interatomic potential with non-zero range, is restricted to low energies and momenta. Accordingly, it must be used together with a UV momentum cutoff as to prevent the theory from ultraviolet divergences that show up in perturbative schemes [36, 16, 15, 37]. In the numerical simulations below, the required UV momentum cutoff is imposed explicitly via the finite computational lattice. If the lattice spacings ( $\Delta x$ , $\Delta y$ , $\Delta z$ ) in each spatial dimension 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$ and the coupling constant $g$ in the contact potential can indeed be expressed in terms of the 3D scattering length $a$ as $g=4\pi\hbar^{2}a/m$ [37, 38, 39]. Physically, the restriction $k_{x,y,z}^{\max}\ll 1/a$ means that the theory with contact interactions cannot be used to describe phenomena on lengthscales shorter $a$ .

To model the dynamics of quantum fields governed by the above Hamiltonian, we use the positive- $P$ phase-space approach [40, 41, 42, 43, 44, 45, 46, 47]. 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) [45, 46, 47]:

$\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),$	(2)
$\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),$	(3)

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 [40, 41]. 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 [48]. This in turn leads to rapid and uncontrollable growth of sampling errors as the simulation progresses in time [44], 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 section 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 or velocity through the use of Madelung’s transformation) [15, 16, 49].

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 or ratios of the stochastic densities factorize into respective products or 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. 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 [50, 51, 52]. 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, 49],

$\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 are usually 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 [53, 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 [48, 44]) 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. The assumption that density and phase fluctuations are small is valid for 3D weakly-interacting Bose-Einstein condensates where true long-range order exists. However, when considering 1D and 2D systems more care must be taken, since phase-fluctuations persist at low temperatures (even at zero temperature in 1D) and true long-range order is absent for uniform systems in the thermodynamic limit [54, 35, 55, 56]. Nevertheless, for finite-size trapped systems it is still possible to employ a linearized treatment of fluctuations in these cases, provided that the first order correlation function decays slowly enough, or in other words, that the system size $L$ is smaller than the relevant phase-coherence length $l_{\phi}$ so that long-range order is recovered within the system [35, 57, 58, 37, 59, 60, 61].

Thus, we proceed by decomposing the density, velocity, and phase variables as a sum of two components: a real-valued deterministic component that obeys the superfluid hydrodynamic equations, and a complex-valued fluctuating component that describes any perturbations on top of the mean-field description. We denote such a decomposition using,

$\displaystyle\rho(\mathbf{r},t)$	$\displaystyle=\rho_{0}(\mathbf{r},t)+\delta\rho(\mathbf{r},t),$	(25)
$\displaystyle\bm{v}(\mathbf{r},t)$	$\displaystyle=\bm{v}_{0}(\mathbf{r},t)+\delta\bm{v}(\mathbf{r},t),$	(26)
$\displaystyle S(\mathbf{r},t)$	$\displaystyle=S_{0}(\mathbf{r},t)+\delta S(\mathbf{r},t),$	(27)

where $\rho_{0}(\mathbf{r},t)$ , $\bm{v}_{0}(\mathbf{r},t)$ , $S_{0}(\mathbf{r},t)$ obey Eqs. (18)–(20), and the evolution equations for $\delta\rho(\mathbf{r},t)$ , $\delta\bm{v}(\mathbf{r},t)$ , $\delta S(\mathbf{r},t)$ are yet to be determined.

Taking the expectation value of the SQHD equations (14), (15) and (10), one can show that, to linear order in fluctuating quantities, the superfluid hydrodynamic equations exactly capture the mean behavior of such a description, i.e. $\langle\rho(\mathbf{r},t)\rangle=\rho_{0}(\mathbf{r},t)$ , $\langle\bm{v}(\mathbf{r},t)\rangle=\bm{v}_{0}(\mathbf{r},t)$ and $\langle S(\mathbf{r},t)\rangle=S_{0}(\mathbf{r},t)$ . This means that, at linear order, $\langle\delta\rho(\mathbf{r},t)\rangle=\langle\delta\bm{v}(\mathbf{r},t)\rangle=\langle\delta S(\mathbf{r},t)\rangle=0$ . Therefore, while a linearized scheme cannot provide corrections to quantities like the real-space density, it instead finds its primary utility in the ability to calculate quantum correlation functions in dynamical scenarios, since leading order corrections can be kept when considering these observables (see Sec. IV.2). We derive such linearized equations shortly, and calculate correlation functions later in this paper, however we pause briefly now to consider fluctuations up to second order. In this case $\langle\delta\rho(\mathbf{r},t)\rangle,\langle\delta\bm{v}(\mathbf{r},t)\rangle,\langle\delta S(\mathbf{r},t)\rangle\neq 0$ and it is therefore possible to obtain quantum corrections to the mean-field density, velocity and phase. Such a second-order scheme remains valid provided that perturbations from the mean-field result can be considered small, that is, $\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)$ , $\langle\delta S(\mathbf{r},t)\rangle\ll S_{0}(\mathbf{r},t)$ . The equations governing the fluctuating fields can be derived by substituting the decompositions (25)–(27) into the SQHD equations [(14), (15) and (10)] and employing the superfluid hydrodynamic equations (18)–(20) for the time evolution of the deterministic components. Then, truncating to second order in fluctuating quantities leads to 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}{\sqrt{\rho_{0}}}\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}{\sqrt{\rho_{0}}}\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. While 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. While 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 or 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^{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)$ , since $\langle\delta\rho(\mathbf{r},t)\rangle=\langle\delta\bm{v}(\mathbf{r},t)\rangle=\langle\delta S(\mathbf{r},t)\rangle=0$ in this case (as mentioned previously).

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. [47] (see also [62]). Such equations have been used to successfully describe several experimental results on Bose-Einstein condensate collisions [47, 63, 64, 65, 66, 67, 68], and a more recent extension of this approach was used to characterize the quantum depletion of an expanding condensate [69]. We point out that the stochastic Bogoliubov equations themselves are equivalent to the usual Bogoliubov approach [15, 16, 70, 71] 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)=e^{-i\mu t/\hbar}\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 [72, 73, 74], where $\mu$ is the chemical potential of the condensate mode, $u_{i}$ and $v_{i}$ are the Bogoliubov mode amplitudes, and $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$ the respective quasiparticle creation and annihilation operators.

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 initialized 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 initializing 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)$ , $\widetilde{\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. [47, 67, 69]. 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. While 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 [75, 76].

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)=\langle\widetilde{\Psi}(\mathbf{r},t)\Psi(\mathbf{r},t)\rangle=\left\langle\rho(\mathbf{r},t)\right\rangle=\rho_{0}(\mathbf{r},t),

(43)

since the superfluid hydrodynamic density $\rho_{0}(\mathbf{r},t)$ captures exactly the mean of $\rho(\mathbf{r},t)$ such that $\langle\delta\rho(\mathbf{r},t)\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. While the scheme is still capable of computing non-trivial results for higher order correlation functions, we recall that 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 [77, 78, 74, 79, 80, 69], the usual BdG equations and stochastic Bogoliubov scheme above rely on the undepleted pump approximation, 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 [81] (see also [82, 83], 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 expands the stochastic field using $\Psi=\Psi_{0}+\delta\psi$ (or $\widetilde{\Psi}=\Psi_{0}^{*}+\delta\tilde{\psi}$ ) and 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}+\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle.$		(44)

Similarly to the LSQHD equations, since the stochastic Bogoliubov equations possess additive noise, the mean-field GPE describes exactly the average of the full stochastic field $\Psi(\mathbf{r},t)$ such that $\langle\Psi(\mathbf{r},t)\rangle=\Psi_{0}(\mathbf{r},t)$ and $\langle\delta\tilde{\psi}(\mathbf{r},t)\rangle=\langle\delta\psi(\mathbf{r},t)\rangle=0$ . Equation (44) directly contrasts the hydrodynamic result $G^{(1)}(\mathbf{r},\mathbf{r},t)=\rho_{0}(\mathbf{r},t)=|\Psi_{0}(\mathbf{r},t)|^{2}$ of Eq. (43) for which the beyond-mean-field term $\langle\delta\tilde{\psi}(\mathbf{r},t)\delta\psi(\mathbf{r},t)\rangle$ (that is responsible for particle number growth) is missing. 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 or perturbation 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 [55, 53],

	$\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 normalized 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 we recall that averages of first order terms are zero for linearized approaches.

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)

since averages of first order terms are zero.

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. [39] 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, 84]), 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.

Refer to caption — 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 quantum 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.

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.

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 $\approx\!17$ times smaller in magnitude compared to the full positive-P approach, and $\approx\!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. While 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. [85] 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 $\approx 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$ [38]. 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\approx 1$ [86, 87, 88], 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. [89]. 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 spatial 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 [90, 91] 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\simeq 0.33l_{\text{h}}$ . 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$ .

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.

In Figure 6, we plot the entire normalized density-density correlation function $g^{(2)}(x,x^{\prime},t)$ 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.

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, although 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 $\approx 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 numerical implementation of the LSQHD scheme, however, is that it is not suitable for situations where the mean-field density $\rho_{0}(\bm{r},t)$ tends 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 [92]. 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, 93, 59, 29, 94, 95]), 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 [96], that might follow from the higher-order expansion terms. Another avenue is to find connections with the theory of nonlinear Luttinger liquids [97, 98] or more generally with the theories of nonlinear fluctuating hydrodynamics in 1D [99]. In these theories, nonlinear terms are important for describing anomalous transport in 1D classical fluids and alter the universality class of hydrodynamic fluctuations [99, 100]. We further emphasize 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 [101, 42]. 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$ [38]). 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|\psi_{W}|^{2}\right]\psi_{W}.

(71)

For the scenarios presented in Section V.1, which are carried out in 1D with $g\rightarrow g_{\text{1D}}$ , each trajectory is initialized in a coherent-state using [42]

\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, and $\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$ [101]. 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 [102]

$\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 normalized 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)$ are 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 [103, 104, 105] 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$ [106, 107, 108]. 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$ [108, 109, 110] 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}$ on the order of $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 [111] with the optimized 4th-order decomposition [112]. 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 on the order of $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\approx 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)

as given in the main text.

Normalizing the reduced one-body density matrix and maintaining up to 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.

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.

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$ .

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.

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(R) (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 [2016] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation and Superfluidity (Oxford University Press, 2016).
Pethick and Smith [2008] C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases, 2nd ed. (Cambridge University Press, Cambridge, 2008).
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. S. 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).
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).
Mora and Castin [2003] C. Mora and Y. Castin, Extension of Bogoliubov theory to quasicondensates, Phys. Rev. A 67, 053615 (2003).
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.
Glauber [2007] R. J. Glauber, Quantum theory of optical coherence: selected papers and lectures (Wiley-VCH, Wieinheim, 2007).
Glauber [1963a] R. J. Glauber, The Quantum Theory of Optical Coherence, Phys. Rev. 130, 2529 (1963a).
Glauber [1963b] R. J. Glauber, Coherent and Incoherent States of the Radiation Field, Phys. Rev. 131, 2766 (1963b).
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).
Hohenberg [1967] P. C. Hohenberg, Existence of Long-Range Order in One and Two Dimensions, Phys. Rev. 158, 383 (1967).
Cazalilla [2004] M. A. Cazalilla, Bosonizing one-dimensional cold atomic gases, Journal of Physics B: Atomic, Molecular and Optical Physics 37, S1 (2004).
Castin [2004] Y. Castin, Simple theoretical tools for low dimension Bose gases, in Journal de Physique IV (Proceedings), Vol. 116 (EDP sciences, 2004) pp. 89–132.
Petrov et al. [2000] D. S. Petrov, G. V. Shlyapnikov, and J. T. M. Walraven, Regimes of Quantum Degeneracy in Trapped 1D Gases, Phys. Rev. Lett. 85, 3745 (2000).
Bouchoule et al. [2012] I. Bouchoule, M. Arzamasovs, K. V. Kheruntsyan, and D. M. Gangardt, Two-body momentum correlations in a weakly interacting one-dimensional Bose gas, Phys. Rev. A 86, 033626 (2012).
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).
Mathey et al. [2009] L. Mathey, A. Vishwanath, and E. Altman, Noise correlations in low-dimensional systems of ultracold atoms, Phys. Rev. A 79, 013609 (2009).
Imambekov et al. [2009] A. Imambekov, I. E. Mazets, D. S. Petrov, V. Gritsev, S. Manz, S. Hofferberth, T. Schumm, E. Demler, and J. Schmiedmayer, Density ripples in expanding low-dimensional gases as a probe of correlations, Phys. Rev. A 80, 033604 (2009).
[62] 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, 3752 (2014).
Ross et al. [2022] J. A. Ross, P. Deuar, D. K. Shin, K. F. Thomas, B. M. Henson, S. S. Hodgman, and A. G. Truscott, On the survival of the quantum depletion of a condensate after release from a magnetic trap, Scientific Reports 12, 13178 (2022).
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).
Billam and Gardiner [2012] T. P. Billam and S. A. Gardiner, Coherence and instability in a driven Bose–Einstein condensate: a fully dynamical number-conserving approach, New Journal of Physics 14, 013038 (2012).
Billam et al. [2013] T. P. Billam, P. Mason, and S. A. Gardiner, Second-order number-conserving description of nonequilibrium dynamics in finite-temperature Bose-Einstein condensates, Phys. Rev. A 87, 033628 (2013).
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(R) (2009).
Kheruntsyan [2006] K. V. Kheruntsyan, Quantum Atom Optics with Fermions from Molecular Dissociation, Phys. Rev. Lett. 96, 110401 (2006).
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).
[85] I. P. McCulloch, The Matrix Product Toolkit, https://people.smp.uq.edu.au/IanMcCulloch/mptoolkit.
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(R) (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).
Ruostekoski and Martin [2013] J. Ruostekoski and A. D. Martin, The Truncated Wigner Method for Bose Gases, in Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics (Imperial College Press, 2013) Chap. 13, pp. 203–214.
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).
v= [0] Moreover, in situations where the density goes precisely to zero, the velocity field is no longer well-defined, and even the classical hydrodynamic equations (23) and (24) are no longer suitable to describe the system in that case.
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).
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.
Watson et al. [2022] R. S. Watson, S. A. Simmons, and K. V. Kheruntsyan, Benchmarks of Generalized Hydrodynamics for 1D Bose Gases, arXiv:2208.06614 [cond-mat.quant-gas] (2022).
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).
Spohn [2016] H. Spohn, Fluctuating hydrodynamics approach to equilibrium time correlations for anharmonic chains, in Thermal Transport in Low Dimensions, edited by S. Lepri (Springer International Publishing, Cham, 2016) pp. 107–158.
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(R) (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(R) (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).