Quantum dynamics of long-range interacting systems using the positive-P and gauge-P representations

TL;DR

This paper develops and tests stochastic positive-P and gauge-P simulation methods for long-range interacting bosonic systems, showing they can effectively extend simulation times and handle large systems without truncating the quantum state.

Contribution

It introduces optimized stochastic gauges and analytical estimates to improve the efficiency and applicability of positive-P and gauge-P methods for long-range quantum systems.

Findings

Long-range interactions do not significantly limit positive-P and gauge-P methods.

Stochastic gauges can extend simulation times effectively.

Different gauges are optimal depending on system size.

Abstract

We provide the necessary framework for carrying out stochastic positive-P and gauge-P simulations of bosonic systems with long range interactions. In these approaches, the quantum evolution is sampled by trajectories in phase space, allowing calculation of correlations without truncation of the Hilbert space or other approximations to the quantum state. The main drawback is that the simulation time is limited by noise arising from interactions. We show that the long-range character of these interactions does not further increase the limitations of these methods, in contrast to the situation for alternatives such as the density matrix renormalisation group. Furthermore, stochastic gauge techniques can also successfully extend simulation times in the long-range-interaction case, by making using of parameters that affect the noise properties of trajectories, without affecting physical observables. We derive essential results that significantly aid the use of these methods: estimates of the available simulation time, optimized stochastic gauges, a general form of the characteristic stochastic variance and adaptations for very large systems. Testing the performance of particular drift and diffusion gauges for nonlocal interactions, we find that, for small to medium systems, drift gauges are beneficial, whereas for sufficiently large systems, it is optimal to use only a diffusion gauge. The methods are illustrated with direct numerical simulations of interaction quenches in extended Bose-Hubbard lattice systems and the excitation of Rydberg states in a Bose-Einstein condensate, also without the need for the typical frozen gas approximation. We demonstrate that gauges can indeed lengthen the useful simulation time.

Tables

Table 1. Table 1: Summary of main analytic results.

. Result	positive-P	gauge-P
Evolution equations	basic (13)	general gauge-P (34)
for nonlocal interactions:	large system (79)	standard drift gauge (35)
Characteristic variances:	(59)	all gauges (37)
		diffusion gauges only (57)
Optimized global gauges:		all gauges (55)
		diffusion gauge only (81)
Simulation time estimate:	(64)	all gauges (66)
		diffusion gauge only (82)
		gauge comparison (84)

Equations

$$\hat{H}=\sum _{nm}\left[{\hat{a}}_n^{†}{\omega }_{nm}{\hat{a}}_m+\frac{1}{2}{\hat{a}}_n^{†}{\hat{a}}_m^{†}{W}_{nm}{\hat{a}}_n{\hat{a}}_m\right].$$

$$W(|\mathbf{x}-\mathbf{y}|)=-\frac{C_6}{(|\mathbf{x}-\mathbf{y}{|}^a+{ϵ}^a{)}^b}.$$

$$W_{nn}=W_0.$$

$$\hat{H}=\frac{\kappa }{2}\int d\mathbf{x}\left({\hat{\Psi }}_e^{†}(\mathbf{x}){\hat{\Psi }}_g(\mathbf{x})+{\hat{\Psi }}_g^{†}(\mathbf{x}){\hat{\Psi }}_e(\mathbf{x})\right)$$

$$+\frac{1}{2}\int d\mathbf{x}\int d\mathbf{y}{\hat{\Psi }}_e^{†}(\mathbf{x}){\Psi }_e^{†}(\mathbf{y})W(\mathbf{x}-\mathbf{y}){\hat{\Psi }}_e(\mathbf{x}){\Psi }_e(\mathbf{y}).$$

$${\hat{H}}_{\mathrm{kin}}=-\int dx\left({\hat{\Psi }}_e^{†}(\mathbf{x})\frac{{\nabla }^2}{2m}{\hat{\Psi }}_e(\mathbf{x})+{\hat{\Psi }}_g^{†}(\mathbf{x})\frac{{\nabla }^2}{2m}{\hat{\Psi }}_g(\mathbf{x})\right)$$

$$\hat{\rho }(t)=\int d^{2M}\mathbf{\alpha }\int d^{2M}\mathbf{\beta }\int d^2\Omega P_G(\mathbf{\alpha },\mathbf{\beta },\Omega ,t)\left[\Omega \frac{|\mathbf{\alpha }⟩⟨{\mathbf{\beta }}^{\ast }|}{⟨{\mathbf{\beta }}^{\ast }|\mathbf{\alpha }⟩}\right].$$

$$\frac{d}{dt}\hat{\rho }=-\frac{i}{\hslash }[\hat{H},\hat{\rho }].$$

$$\frac{\partial P_G}{\partial t}=-\sum _{\mu }\frac{\partial }{\partial v_{\mu }}{A}_{\mu }{P}_G+\frac{1}{2}\sum _{\mu \nu }\frac{\partial }{\partial v_{\mu }}\frac{\partial }{\partial v_{\nu }}{D}_{\mu \nu }{P}_G$$

$${\hat{a}}_n\hat{\rho }\leftrightarrow {\alpha }_n{P}_G,{\hat{a}}_n^{†}\hat{\rho }\leftrightarrow \left({\beta }_n-\frac{\partial }{\partial {\alpha }_n}\right)P_G,$$

$$\hat{\rho }{\hat{a}}_n\leftrightarrow \left({\alpha }_n-\frac{\partial }{\partial {\beta }_n}\right)P_G,\hat{\rho }{\hat{a}}_n^{†}\leftrightarrow {\beta }_n{P}_G,$$

$$0\leftrightarrow P_G+\frac{\partial }{\partial \Omega }\Omega P_G.$$

$$\frac{d{v}_{\mu }}{dt}=A_{\mu }+\sum _{\nu }{B}_{\mu \nu }{\xi }_{\nu }(t),$$

$$⟨{\hat{a}}_n^{†}{\hat{a}}_m^{†}{\hat{a}}_l^{†}\dots {\hat{a}}_p{\hat{a}}_q{\hat{a}}_r\dots ⟩=\frac{\overline{\Omega {\beta }_n{\beta }_m{\beta }_l\dots {\alpha }_p{\alpha }_q{\alpha }_r\cdots +{\Omega }^{\ast }{\alpha }_n^{\ast }{\alpha }_m^{\ast }{\alpha }_l^{\ast }\dots {\beta }_p^{\ast }{\beta }_q^{\ast }{\beta }_r^{\ast }\dots }}{\overline{\Omega +{\Omega }^{\ast }}}.$$

$$\underline{\underline{D}}\,=\left(\begin{array}[]{cc}D^{(\alpha)}&0\\ 0&D^{(\beta)}\end{array}\right),\qquad D^{(\alpha)}_{nm}=-iW_{nm}\alpha_{n}\alpha_{m},\qquad D^{(\beta)}_{nm}=iW_{nm}\beta_{n}\beta_{m}.\qquad$$

$$\frac{d{\alpha }_n}{dt}$$

$$\frac{d{\beta }_n}{dt}$$

$$\displaystyle\underline{\bm{\gamma}}=\left(\begin{array}[]{c}\bm{\alpha}\\ \bm{\beta}\end{array}\right),\>\>\>\underline{\underline{\omega}}\,=\left(\begin{array}[]{cc}-\omega&0\\ 0&\omega\end{array}\right),\>\>\>\underline{\underline{W}}\,=\left(\begin{array}[]{cc}-W&0\\ 0&W\end{array}\right),\>\>\>$$

$$\displaystyle\underline{\bm{n}}=\left(\begin{array}[]{c}\bm{n}\\ \bm{n}\end{array}\right)=\left(\begin{array}[]{c}(\dots,\alpha_{j}\beta_{j},\dots)^{T}\\ (\dots,\alpha_{j}\beta_{j},\dots)^{T}\end{array}\right),\>\>\>\underline{\bm{A}}=i\underline{\underline{\omega}}\,\underline{\bm{\gamma}}+i\underline{\underline{\Gamma}}\,\underline{\underline{W}}\,\underline{\bm{n}},$$

$$\displaystyle\underline{\bm{\xi}}=\left(\begin{array}[]{c}\bm{\xi}^{(1)}\\ \bm{\xi}^{(2)}\end{array}\right),\>\>\>\underline{\underline{\Gamma}}\,=\mbox{diag}(\underline{\bm{\gamma}}),\>\>\>\underline{\underline{B}}\,=\sqrt{i}\underline{\underline{\Gamma}}\,\underline{\underline{S}}\,,\>\>\>\underline{\underline{S}}\,=\left(\begin{array}[]{cc}-i\sqrt{W}&0\\ 0&\sqrt{W}\end{array}\right)\!.$$

$$n_m=n_m^{\prime }+i{n}_m^{\prime \prime }={\alpha }_m{\beta }_m.$$

$$\frac{\partial }{\partial t}\underline{\mathbf{\gamma }}=\underline{\mathbf{A}}+\underline{\underline{B^{\prime }}}(\underline{\mathbf{\xi }}-\underline{\mathbf{g}}),\frac{\partial }{\partial t}\Omega =\Omega \underline{\mathbf{g}}\cdot \underline{\mathbf{\xi }},$$

$$\displaystyle\underline{\underline{B^{\prime}}}\,=\underline{\underline{B}}\,\underline{\underline{O}}\,=\sqrt{i}\underline{\underline{\Gamma}}\,\underline{\underline{S}}\,\underline{\underline{O}}\,,\>\>\>\>\>\>\underline{\bm{g}}=\left(\begin{array}[]{c}\bm{g}^{(1)}\\ \bm{g}^{(2)}\end{array}\right),$$

$$g_{\mu }=\sqrt{i}\sum _{\nu \sigma }{O}_{\nu \mu }{S}_{\sigma \nu }{f}_{\sigma },$$

$$\frac{\partial }{\partial t}{\gamma }_{\mu }=i\sum _{\nu }{\omega }_{\mu \nu }{\gamma }_{\nu }+i\sum _{\nu }{\gamma }_{\mu }{W}_{\mu \nu }(n_{\nu }-f_{\nu })+\sqrt{i}\sum _{\nu \sigma }{\gamma }_{\mu }{S}_{\mu \nu }{O}_{\nu \sigma }{\xi }_{\sigma },$$

$$\frac{\partial }{\partial t}\Omega =\sqrt{i}\Omega \sum _{\nu \sigma \lambda }{f}_{\lambda }{O}_{\sigma \nu }{S}_{\lambda \sigma }{\xi }_{\nu }=\sqrt{i}\Omega {\underline{\mathbf{f}}}^T\underline{\underline{S}}\underline{\underline{O}}\underline{\mathbf{\xi }}.$$

$$\underline{\mathbf{f}}=i{\underline{\mathbf{n}}}^{\prime \prime }$$

$$\mathcal{V}=\frac{1}{2M}\sum _{\mu }^{2M}\text{mbox}var\left[\log |\Omega {\gamma }_{\mu }|\right].$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\cal V}(t)={\cal V}(0)+\frac{1}{2M}\Bigg{\{}\frac{t}{2}\mbox{Tr}\left[\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger}\right]+t\,\mbox{Tr}\left[{\rm Im}\left\{\left[\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger}\right]\right\}\underline{\underline{N}}\,^{\prime\prime}\right]$$

$$\displaystyle+M\mbox{Tr}\bigg{[}{\rm Re}\left\{\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger}\right\}\underline{\underline{Q}}\,\bigg{]}\Bigg{\}},$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!Q_{\mu\nu}=\frac{1}{2}{\rm Re}\left[\overline{n_{\mu}(0)n^{*}_{\nu}(0)}\Big{(}\exp\left[t\,P_{\mu\nu}\right]-1-t\,P_{\mu\nu}\Big{)}/P_{\mu\nu}\right]+t\,[\overline{n^{\prime\prime}_{\mu}(0)n^{\prime\prime}_{\nu}(0)}]$$

$$P_{\mu \nu }={\left[\underline{\underline{F}}\underline{\underline{S}}\underline{\underline{O}}\underline{\underline{O}}{}^{†}\underline{\underline{S}}{}^{†}\underline{\underline{F}}\right]}_{\mu \nu }.$$

$$\underline{\underline{F}}\,=\left(\begin{array}[]{cc}\mathds{1}&\mathds{1}\\ \mathds{1}&\mathds{1}\end{array}\right),$$

$$\underline{\underline{N}}=\underline{\underline{N}}{}^{\prime }+i\underline{\underline{N}}{}^{\prime \prime }=\text{mbox}diag[\overline{\underline{\mathbf{n}}(0)}].$$

$$\displaystyle\underline{\underline{O}}\,=\left(\begin{array}[]{cc}\cosh a\mathds{1}&-i\sinh a\mathds{1}\\ i\sinh a\mathds{1}&\>\>\cosh a\mathds{1}\end{array}\right)=\exp(g),\>\>\>\>\>\>\ g=\left(\begin{array}[]{cc}0&-ia\mathds{1}\\ ia\mathds{1}&0\end{array}\right).$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}=\left(\begin{array}[]{cc}\>\mathds{1}\cosh 2a&-i\mathds{1}\sinh 2a\\ i\mathds{1}\sinh 2a&\>\>\mathds{1}\cosh 2a\end{array}\right),\qquad\underline{\underline{P}}\,=2e^{-2a}\left(\begin{array}[]{cc}U&U\\ U&U\end{array}\right).$$

Full text

We derive essential results that significantly aid the use of these methods: estimates of the available simulation time, optimized stochastic gauges, a general form of the characteristic stochastic variance and adaptations for very large systems. Testing the performance of particular drift and diffusion gauges for nonlocal interactions, we find that, for small to medium systems, drift gauges are beneficial, whereas for sufficiently large systems, it is optimal to use only a diffusion gauge.

The methods are illustrated with direct numerical simulations of interaction quenches in extended Bose-Hubbard lattice systems and the excitation of Rydberg states in a Bose-Einstein condensate, also without the need for the typical frozen gas approximation. We demonstrate that gauges can indeed lengthen the useful simulation time.

I Introduction

First-principles treatments of many-body quantum problems are notoriously difficult due to the exponential increase of Hilbert space dimension with the number of system components. Tackling this complexity is an outstanding goal of theoretical physics. Long-range interactions usually amplify the difficulties, because they can break symmetries or frustrate many of the tricks used to reduce the Hilbert space to manageable sizes. For example, the entanglement area laws can be broken Eisert et al. (2010).

If limited precision is acceptable, stochastic phase-space methods are a promising contender compared to other standard methods such as Monte Carlo or entangled pair states Steel et al. (1998); Gardiner and Zoller (2004); Drummond et al. (2007), especially for higher-dimensional systems. The computational complexity of phase-space methods tends to scale only linearly or quadratically in system size and to be largely independent of dimensionality. Originating from quantum optics Drummond and Gardiner (1980), these methods have also been successful for cold degenerate gases Kheruntsyan et al. (2005); Savage and Kheruntsyan (2007); Carusotto et al. (2001); Carusotto and Castin (2001); Deuar and Drummond (2007); Drummond et al. (2004); Polkovnikov (2010); Deuar et al. (2011, 2013); Kheruntsyan et al. (2012); Lewis-Swan and Kheruntsyan (2014, 2015) including fermions Corney and Drummond (2004a, 2006) and spin systems Barry and Drummond (2008); Ng et al. (2013); Schachenmayer et al. (2015). They are particularly suited to relatively dilute boson systems deep in the quantum regime.

Apart from a few small forays Drummond et al. (2007); Deuar (2005a), previous work with phase-space methods has been limited to systems with contact interactions. However, long-range interactions have recently become important in the dilute quantum regime due to advances in the production of ultracold Rydberg Gallagher (1994); Löw et al. (2012), dipolar atomic Griesmaier et al. (2005); Lahaye et al. (2009) and molecular Carr et al. (2009) gases, whose physics can not be captured with a contact interaction alone. Motivated by this experimental progress, we develop here tools for the application of phase-space methods to systems with long range interactions. These tools could then go substantially beyond the initial exploration of the positive-P method to a Rydberg system, reported in Wu et al. (2016), which required a simplified interaction model.

The price paid for use of phase-space methods is that the dynamics is tractable usually only over short time scales. This limitation arises from the nonlinear amplification of stochastic fluctuations in individual trajectories, leading to a phase-space distribution that is too broad Gilchrist et al. (1997); Deuar and Drummond (2002); Deuar and Drummond (2006a); Carusotto and Castin (2003). Hence, a central requirement in practice is to estimate the time limits and if possible to extend them. Fortunately, one can exploit the over-completeness of the basis used to keep the distribution as compact as possible and to stabilize trajectories. Simulation schemes exploiting this feature have been termed Gauge-P methods Plimak et al. (2001); Deuar and Drummond (2002); Dowling et al. (2005); Deuar and Drummond (2003); Deuar and Drummond (2006b); Deuar and Drummond (2001); Deuar (2004); Drummond et al. (2004); Deuar and Drummond (2007). Two different types of gauges are available. Drift-gauges modify the drift term of the underlying Fokker-Planck equation for the phase-space distribution, allowing e.g. a modification of nonlinearities in the evolution equation at the expense of introducing a fluctuating trajectory weight. Diffusion gauges alter the noise terms in the stochastic evolution equations and can be used to vary the relative fraction of noise affecting the phase or amplitude of complex variables.

These stochastic gauges can become particularly powerful when it is possible to determine the best choice of gauge for a given simulation. For problems that rely essentially on single bosonic modes, these best choices have been found Deuar (2004); Deuar and Drummond (2006b). For the many-mode case to date, optimization has only been considered insofar as a collection of single modes can offer some guidance, in situations where the interactions are weak or short range.

In this work, we adapt the prescriptions and approaches that are successful for short range interactions to long-range interactions, and we show that they still give useful extensions of simulation time. To this end, we calculate the evolution of a characteristic variance that can serve to estimate the effect of a specific gauge, analogous to the approach of Deuar and Drummond (2006b) for the single-mode case. Though quite an involved expression, it provides the essential basis for the application of Gauge-P techniques to problems with long-range interactions. Estimates of the available simulation time and optimal adaptive gauge choices then follow in a relatively straightforward way.

The article is organized as follows: In section II we describe the class of many-body quantum mechanical models that are being considered here. The general idea behind the positive-P and Gauge-P methods is then briefly reviewed in section III, and the evolution equations for nonlocal interactions derived. In section IV, which contains the core of our analysis, we derive and analyze the characteristic noise variances and optimize gauges. Based on this, section V presents estimates of the available simulation time.

All these results are then applied to two test cases in section VI: an interaction quench in an extended (long-range) Bose-Hubbard model and the excitation dynamics of mobile Rydberg atoms within a BEC. These examples are oriented towards long-range interacting ultra-cold atomic physics, such as Bose-Einstein condensates (BEC) of dipolar atoms Griesmaier et al. (2005), polar molecules Carr et al. (2009) or Rydberg dressed condensates Santos et al. (2000); Henkel et al. (2010); Pupillo et al. (2010); Honer et al. (2010); Maucher et al. (2011); Wüster et al. (2011a); Johnson and Rolston (2010); Balewski et al. (2014); Jau et al. (2016). Following this, section VII presents some further results that facilitate simulation of very large systems: a more efficient representation of the noise, and an optimization of pure diffusion gauges. In section VIII, we report also on some initial explorations of more general nonlocal gauges as a starting-point for future work. Appendices contain various technical details, such as the Stratonovich corrections (A), the complete derivation of the characteristic variances (B and C), and supporting material on non-local noises (D) and gauges (E).

II Long-range interacting systems

The class of Hamiltonian operators with which we work in this article is represented by

[TABLE]

The indices label elements of the chosen single particle basis $|m\rangle$ and the operator $\hat{a}_{m}$ creates a particle in that basis state. Single particle energies and linear mode couplings are included in $\omega_{nm}$, while $W_{nm}$ parametrizes generic two-body interactions. This potential $W$ is assumed to be real here. For Hermiticity, $\omega_{nm}=\omega_{mn}^{*}$.

We now introduce two physical model systems the Hamiltonian of which has the form (1). For these examples, we show simulation results in section VI.

II.1 Extended Bose-Hubbard model

One Hamiltonian of the form Eq. (1) is that of an extended (long-range interacting) Bose-Hubbard model. In a cold atom context, $|m\rangle$ is then the localized, lowest-band, Wannier function at a single site of an optical lattice. The coupling $\omega_{nm}=J(\delta_{n+1,m}+\delta_{n-1,m})$ describes particle tunneling between neighboring lattice sites with amplitude $J$. Many cases of long-range interactions can be captured by the potential

[TABLE]

Possible experimental realizations include optical lattices filled with polar atoms Griesmaier et al. (2005) or Rydberg dressed atoms Gil et al. (2014); Macrì and Pohl (2014); Li et al. (2013); Zeiher et al. (2016). The $W_{nm}$ in (1) can then be obtained from Eq. (2) by use of overlap integrals involving Wannier basis functions. However, here we simply use $W_{nm}=W(r)$, where $r=|\mathbf{x}-\mathbf{y}|$ and $\mathbf{x}_{n}$, $\mathbf{x}_{m}$ are discrete lattice points. We will later present simulation results for $a=2$ and $b=3$, corresponding to a softened van-der-Waals potential. Alternatively we could have used $a=6$ and $b=1$ to represent Rydberg-dressed potentials that include the cut-off due to blockade effects Henkel et al. (2010). The contact part of the interaction plays a prominent mathematical role in the following, despite the presence of long range forces. For a homogenous interaction, we abbreviate it with

[TABLE]

With the many-body model above, we will consider controlled interaction quenches to illustrate our results in section VI.1. These are readily realized with e.g. Rydberg atom interactions due to the possibility of time-dependent external control.

II.2 Rydberg excitations in a Bose-gas

As a second example, we will consider a Bose-gas where atoms are continuously coherently excited from a groundstate $|g\rangle$, in which interactions are neglected, to a Rydberg state $|e\rangle$ with strong long-range interactions between the atoms Singer et al. (2004); Tong et al. (2004); Heidemann et al. (2007); Johnson et al. (2008); Tanner et al. (2008); Raitzsch et al. (2008). Initially all atoms are in the ground-state. With more atoms getting excited, the so-called dipole blockade develops Singer et al. (2004); Tong et al. (2004); Heidemann et al. (2007), where the simultaneous excitation of nearby atoms is strongly suppressed.

For this system we use a continuum formulation of the Hamiltonian (1). To this end, we define field operators via $\hat{a}_{n}=\hat{\Psi}(\mathbf{x}_{n})\sqrt{dV}$, $\hat{a}_{m}=\hat{\Psi}(\mathbf{y}_{m})\sqrt{dV}$, $W_{nm}=W(\mathbf{x}_{n}-\mathbf{y}_{m})$, and $\omega_{nm}=-(\delta_{n+1,m}+\delta_{n-1,m}-2\delta_{n,m})/(2m\>dV^{2})$. Here, $\hat{\Psi}(\mathbf{x})$ is a field operator creating an atom at a location $\mathbf{x}$, $dV$ is the (physically small) volume corresponding to one discrete location $\mathbf{x}_{n}$ on the numerical lattice. We set $\hbar=m=1$ for simplicity, unless stated otherwise.

We then extend our model to include two fields describing atoms in components $|g\rangle$ and $|e\rangle$. Only the latter shall experience the long-range interactions. Rydberg excitation and dynamics often happen much faster than atomic motion, which justifies the so-called “frozen gas” approximation where atomic motion is neglected. With the coherent coupling amplitude $\kappa$ between the components we then have:

[TABLE]

The crucial difference to the extended Bose-Hubbard model is that the number of atoms in the component experiencing the interactions is no longer conserved. We will assume the same interaction potential $W$ as in the previous model (2).

For short times, Eq. (4) represents an adequate model (see e.g. Möbius et al. (2013)). However corrections due to atomic motion can easily become relevant Wüster et al. (2010a, 2011b); K Leonhardt and S Wüster and J.-M. Rost (2014); Leonhardt et al. (2017); Thaicharoen et al. (2015); Thaicharoen et al. (2016a); Celistrino Teixeira et al. (2015), in which case the kinetic-energy term

[TABLE]

must be reinstated to Eq. (4), precluding the use of a spin model. This reduces the range of methods available, but does not constitute a problem for phase-space methods.

With or without motion, we will later be interested in correlation function dynamics at the onset of a dipole-blockade: Assume all atoms are initially condensed in the groundstate. The coupling $\kappa$ then acts for a given time, transferring population into the excited state, which has interactions between atoms. The interaction energy of atoms in the excited state can become large enough to suppress transition to that state after an initial one. As a result, there is only one atom excited per “blockade volume”. The Rydberg atom autocorrelation function $g^{(2)}$ has a pronounced dip for radii less than the blockade radius $r_{b}$ and the frequency of incomplete Rabi oscillations increases to $\kappa_{mb}=\sqrt{N_{b}}\kappa$, where $N_{b}$ is the number of ground state atoms within the blockade volume and $\kappa$ is the single atom Rabi frequency Lukin et al. (2001).

III Gauge-P and positive-P descriptions

If we consider $M$ single particle modes $|\,{m}\,\rangle$, populated with up to $N$ particles, the number of many-body basis states of our problem scales like $M^{N}$, an astronomical number for all but the smallest systems. The first-principles simulation of the Hamiltonian (1) is thus a severe challenge. We now briefly review how the dynamics can be tackled with the gauge-P distribution, by stochastically sampling the many-body state, rather than representing it exactly. For an in-depth explanation we refer to the extensive literature on the Positive-P Drummond and Gardiner (1980); Steel et al. (1998); Gardiner and Zoller (2004); Drummond et al. (2007); Drummond and Corney (1999); Savage and Kheruntsyan (2007); Deuar and Drummond (2007); Polkovnikov (2010); Deuar et al. (2011); Ng and Sørensen (2011); Gilchrist et al. (1997); Carter et al. (1987); Deuar and Drummond (2006a); Deuar (2009); Świsłocki and Deuar (2016) and Gauge-P methods Plimak et al. (2001); Deuar and Drummond (2002); Carusotto et al. (2001); Carusotto and Castin (2001); Deuar and Drummond (2003); Deuar and Drummond (2006b); Deuar and Drummond (2001); Deuar (2004); Aimi and Imada (2007); Drummond et al. (2004); Deuar et al. (2009); Dowling et al. (2007); Ng et al. (2013).

III.1 Phase-space representation

First, we expand the system’s density matrix $\hat{\rho}$ in an over-complete operator basis of coherent states:

[TABLE]

The operator kernel in this expression is $\hat{\Lambda}=\Omega|\bm{\alpha}\rangle\langle\bm{\beta}^{*}|/\langle\bm{\beta}^{*}|\bm{\alpha}\rangle$, where we have made use of many-mode coherent states111The many-mode coherent states are a separable tensor product of coherent states $\otimes_{n=1}^{M}|\alpha_{n}\rangle$ with complex amplitude $\alpha_{n}$ for each of the single particle modes. $|\bm{\alpha}\rangle$ defined by $\hat{a}_{n}|\bm{\alpha}\rangle=\alpha_{n}|\bm{\alpha}\rangle$, with $\alpha_{n}\in{\mathbb{C}}$. The function $P_{G}$ is the “Gauge-P” distribution, which can be chosen positive and real for any density matrix $\hat{\rho}$, similarly to the positive-P distribution $P_{+}(\bm{\alpha},\bm{\beta},t)$, which is the special case that occurs when the global weight $\Omega$ is chosen to be constant.

For unitary dynamics, on which we will focus here, the density matrix obeys the von Neumann equation

[TABLE]

Dissipative environments that would add, for example, Lindblad terms to Eq. (7) and produce a master equation, can be straightforwardly included.

Due to the doubling of the phase space dimension for quantum systems seen above, it will be convenient to introduce the following notation: We collect all $2M+1$ stochastic fields into a vector $\underline{\bm{v}}=(\bm{\alpha},\bm{\beta},\Omega)$ whose components $v_{\mu}$ are labeled by greek indices $\mu,\nu,\dots$, as opposed to the $M$ modes which are labeled by roman indices $n,m,j,\dots$. Continuing this distinction for vectors in general, bold-face, underlined vectors will have $2M$ components222or $2M+1$ components if the weight $\Omega$ is included., without the underline they have $M$ components. Similarly a double underline is used to distinguish $2M\times 2M$ matrices from the $M\times M$ matrices that they are usually composed of. Components of extended matrices have greek subscripts that run over $2M$ values.

For a sufficiently well bounded distribution $P_{G}(\bm{\alpha},\bm{\beta},\Omega)$, the master equation can be converted into a Fokker-Planck equation (FPE) of the generic form

[TABLE]

by virtue of the operator correspondences Hope and Olsen (2001); Deuar (2004):

[TABLE]

Crucially, Eq. (8) is equivalent to the set of coupled stochastic differential equations (SDEs)

[TABLE]

if the diffusion matrix $\underline{\underline{D}}\,$ is positive semi-definite 333A positive semi-definite diffusion matrix has all eigenvalues non-negative, and such a form for $D$ can always be constructed here due to $\hat{\Lambda}$ being an analytic function of complex variables Drummond and Gardiner (1980).  Gardiner and Zoller (2004); Gardiner (2003). These SDEs have to be interpreted in the Itô-form of stochastic calculus Gardiner (2003). The $\xi_{\nu}(t)$ are independent real noise fields, defined by a zero mean $\overline{\xi_{\nu}(t)}$ and correlations $\overline{\xi_{\mu}(t)\xi_{\nu}(t^{\prime})}=\delta_{\mu,\nu}\delta(t-t^{\prime})$, where $\overline{f}$ denotes the stochastic average of $f$. They are usually implemented with independent Gaussian noises at each time step $\Delta t$ with a variance of $1/(\Delta t\>dV)$. The “noise matrix” $\underline{\underline{B}}\,$ must satisfy $\underline{\underline{D}}\,=\underline{\underline{B}}\,\underline{\underline{B}}\,^{T}$. The matrix square root $\underline{\underline{D}}\,^{1/2}$ is a viable noise matrix in most cases of interest, provided that it is self-transposed ($\underline{\underline{D}}\,^{1/2}=[\underline{\underline{D}}\,^{1/2}]^{T}$).

All quantum expectation values (observables) can be found in principle by averaging 444Pathological cases can occur if the distribution $P_{G}$ does not decay sufficiently rapidly for large $|\alpha|$, $|\beta|$. This is discussed e.g. in Gilchrist et al. (1997); Deuar and Drummond (2002); Deuar (2005b). over the set of independent stochastic trajectories. Expectation values of normal ordered operator products take the form Deuar and Drummond (2002)

[TABLE]

The discussion above also applies to the Positive-P distribution, in which case we can set $\Omega\to 1$ and simplify the calculation considerably Deuar and Drummond (2006a).

Upon reaching Eq. (10), we have reduced the generally intractable quantum-many body problem in an $N^{M}$ dimensional Hilbert-space to the seemingly much easier solution of a system of $2M$ coupled complex SDEs. Whether the problem in this new form is tractable depends on the number of stochastic realizations required to obtain converged results for Eq. (11). Therein lies the fundamental limitation of stochastic phase-space methods: After a certain simulated time $t_{\rm sim}$, the noise in the simulation can become nonlinearly amplified which prevents convergence of the averages. As a consequence, these stochastic methods are only useful for understanding properties of the system that manifest themselves before $t_{\rm sim}$.

Fortunately, the above derivation of the equations of motion contains several mathematical degrees of freedom, that do not affect the simulated physics. These degrees of freedom are termed “gauges”. Under favorable conditions, they can be exploited to limit the growth of noise and obtain simulations that give results for longer times $t_{\rm sim}$. The following options for gauges are available: Firstly, we can use the replacement rule (9c) to add functions of our choice to the drift vector $A_{\mu}$ for the evolution of the $\bm{\alpha}$ and $\bm{\beta}$. These functions are the drift gauges. As a trade-off, their use induces a stochastic evolution of the global trajectory weight $\Omega$. Secondly, we can transform the noise matrix $\underline{\underline{B}}\,$ via $\underline{\underline{B}}\,^{\prime}=\underline{\underline{B}}\,\underline{\underline{O}}\,$ with an orthogonal matrix $\underline{\underline{O}}\,$ such that $\underline{\underline{O}}\,\underline{\underline{O}}\,^{T}=I$. If $\underline{\underline{B}}\,$ had fulfilled the diffusion condition $\underline{\underline{D}}\,=\underline{\underline{B}}\,\underline{\underline{B}}\,^{T}$, then clearly so does $\underline{\underline{B}}\,^{\prime}$. An orthogonal matrix $\underline{\underline{O}}\,$ can be written as $\underline{\underline{O}}\,=\exp{\underline{\underline{G}}\,}$, where $\underline{\underline{G}}\,$ is anti-symmetric ($\underline{\underline{G}}\,^{T}=-\underline{\underline{G}}\,$). It has been shown that the real part of $\underline{\underline{G}}\,$ merely induces inconsequential re-summations of the noise, whereas the imaginary part redistributes the noise between amplitude and phase components of the stochastic fields and can have strong influence on the available simulation times Deuar and Drummond (2002, 2003); Deuar (2004); Plimak et al. (2001). The matrices $\underline{\underline{G}}\,$ or $\underline{\underline{O}}\,$ are called diffusion gauges.

Different gauges yield different stochastic equations of motion, whose solutions have different convergence properties. Nonetheless, the evolution still corresponds to one and the same physical density matrix, and averages such as (11) with physical meaning remain unique. Thus gauges affect the mathematical but not the physical properties of our problem, just as their namesake in electro-dynamics.

III.2 Direct form of the equations of motion

Now we apply this formalism to the lattice system described by Hamiltonian (1). In the positive-P representation (without gauges), the diffusion matrix is block-diagonal with separate blocks for the $\alpha$ and $\beta$ variables Gardiner and Zoller (2004); Deuar (2004)., i.e.

[TABLE]

Realistic two-body interactions are symmetric, $W_{nm}=W_{mn}$, which leads to symmetric roots $\sqrt{W}=W^{1/2}=[\sqrt{W}]^{T}$. Then, the simplest decomposition of the noise matrix $\underline{\underline{D}}\,$ gives the following equations:

[TABLE]

The noises $\xi^{(1)}_{j}$ and $\xi^{(2)}_{j}$ are independent, and $[\sqrt{W}]_{nj}$ are components of the matrix square root $\sqrt{W}$ of the interaction potential.

For the gauge-P representation, with stochastic gauge freedoms included, it will be convenient to use the matrix notation described at the beginning of section III.1 to conscisely deal with the $2\times 2$ block nature of the problem. Vectors are assumed to be column vectors unless explicitly transposed, and we use the notation $\mbox{diag}[\bm{v}]$ do assemble a square diagonal matrix with vector $\bm{v}$ on the diagonal. Following Corney and Drummond (2003); Corney and Drummond (2004b), we introduce the following:

[TABLE]

A particularly useful quantity that will recur is the (complex) mode occupation

[TABLE]

Using these definitions, we can write the Gauge-P representation equations of motion

[TABLE]

as per the formalism outlined in Plimak et al. (2001); Deuar and Drummond (2002, 2003); Deuar and Drummond (2006b); Deuar and Drummond (2001); Deuar (2004); Drummond et al. (2004); Deuar and Drummond (2007).

The drift-gauge $\underline{\bm{g}}$ is a yet unspecified function of the variables $\underline{\bm{\gamma}}$. We also have already inserted an arbitrary orthogonal matrix $\underline{\underline{O}}\,=\exp[\underline{\underline{G}}\,]$ with $\underline{\underline{O}}\,\underline{\underline{O}}\,^{T}=\underline{\underline{\mathds{1}}}\,$, the diffusion gauge, into the diffusion term.

IV Optimised gauges

For practical applications of the available stochastic gauges, we introduce more specific forms and then optimize them, similarly to how it was done for local interactions in Deuar and Drummond (2006a, b).

IV.1 Intuitive drift gauge form

Here, we introduce a specific form of $\bm{g}$ and describe its effect on the equations of motion. The objective of $\bm{g}$ is to control the tendency of the interaction term to create so-called moving instabilities and boundary term errors in the distribution $P_{G}$ Deuar and Drummond (2002, 2003); Deuar (2004). Let

[TABLE]

where $\bm{f}$ is a yet unspecified vector that becomes the gauge. In matrix notation this reads $\underline{\bm{g}}=\sqrt{i}\,\underline{\underline{O}}\,^{T}\underline{\underline{S}}\,^{T}\underline{\bm{f}}$. Using the properties of $\underline{\underline{O}}\,$ and $\underline{\underline{S}}\,$,such as $\underline{\underline{S}}\,\underline{\underline{S}}\,^{T}=\underline{\underline{S}}\,\underline{\underline{S}}\,=\underline{\underline{W}}\,$, one can see that the new equations of motion are

[TABLE]

Using the following notation for the real and imaginary part of complex numbers: $z=z^{\prime}+iz^{\prime\prime}$, we can now pick for instance $f_{\lambda}=in_{\lambda}^{\prime\prime}$ to remove the imaginary part of the stochastic density $\underline{\bm{n}}$ from the non-linear interaction term of Eq. (32), or $f_{\lambda}=n_{\lambda}-|n_{\lambda}|$ to replace $\underline{\bm{n}}$ by its modulus. Both choices reduce the tendency of the nonlinear term in Eq. (32) to result in exponentially diverging trajectories Deuar and Drummond (2002).

In the following sections, we will make the choice

[TABLE]

for the drift gauge. This is a form that was found to be particularly useful for the dynamics of the contact-interacting gas in previous work, giving the best extensions of simulation time Deuar and Drummond (2006b).

Most of the efficient numerical time-propagation schemes for stochastic differential equations require the use of the Stratonovich form. The necessary corrections to convert Eq. (34) from the Itô to the Stratonovich form are described in A.

IV.2 Ensemble spread and characteristic variance

To use the diffusion gauge beneficially, we have to know how a given choice for the matrix $\underline{\underline{O}}\,$ affects the noise induced spread of trajectories in the ensemble as time progresses. The trajectory spread determines the simulation time available. A characteristic variance ${\cal V}_{1}=\frac{1}{2}(\mbox{var}[\log|\Omega\alpha|]+\mbox{var}[\log|\Omega\beta|])$ was used for the single mode case in the past to obtain quite accurate estimates of the simulation time $t_{\rm sim}$ Deuar and Drummond (2006b). The quantities $\Omega\gamma$ appearing in ${\cal V}_{1}$, are the estimators that appear in Eq. (11) for expectation values of the one-point correlation function. They are fundamental to most quantities that one might want to calculate. The logarithm is used in anticipation of an exponential increase in noise variances throughout the simulation, as is observed in practice. The first task here, and a keystone of subsequent analysis, is to generalize this quantity for our purposes.

A many mode generalization of ${\cal V}_{1}$ is

[TABLE]

The stochastic uncertainty of fields in all modes contributes to ${\cal V}$, hence this simple scalar quantity can signal when the simulation becomes too noisy to extract useful information via Eq. (11). Even if just the variance of $\gamma_{\mu}$ for a specific mode $\mu$ becomes too large, this will soon contaminate all the remaining modes through the linear coupling terms $\propto\omega$ in the equations of motion. Note, that ${\cal V}$ does not have any physical interpretation, and indeed must not, since it should be gauge-dependent.

For the “standard” drift gauge (35), the approximate time evolution of Eq. (36) is calculated from Eq. (34), in B. The calculation neglects the kinetic $\omega_{\mu\nu}$ part of the drift, as was done for the two-mode case in Deuar and Drummond (2006b), expecting the predominant contribution to the increase of ${\cal V}$ to come from the noise terms and non-linearities. The remaining terms conserve the expectation value of the local density $\langle n_{m}\rangle$, so that ${\cal V}(t)$ depends only on their initial values $n_{m}(0)$. The result is

[TABLE]

This is central to most of the rest of the paper. A useful $2M\times 2M$ matrix is

[TABLE]

and $[\cdots]_{\mu\nu}$ denotes the component $\mu\nu$ of the matrix product within square brackets. The auxilliary diagonal matrix $\underline{\underline{N}}\,$ is

[TABLE]

and has the initial density on the diagonal.

The first term of Eq. (37a) is the contribution from noise placed directly into the amplitudes $\underline{\bm{\gamma}}$. The third term, involving $\underline{\underline{Q}}\,$, arises from noise accumulated in the weight $\Omega$. The remaining term, involving $\underline{\underline{N}}\,^{\prime\prime}$, comes from an accumulation of noise in $\Omega$ due to any initial fluctuations in $\bm{n}^{\prime\prime}(0)$ and is less important.. For single-mode repulsive contact interactions we have $W_{nm}=(g/dV)\delta_{nm}>0$, and thus Eq. (37a) reduces to the known results Deuar and Drummond (2006b).

The variance Eq. (37) forms the essential basis for most other results of this article along with Eq. (57), and may additionally contain useful information for future work on stochastic gauges. The calculation is involved but indispensable to (i) determine simulation time limitations that help to assess whether the methods here are applicable to a given physical problem, (ii) assess whether gauges are expected to outperform the direct application of the (ungauged) Positive-P method, (iii) operate diffusion gauges in analogy to the single mode case and (iv) possibly enable future work to fully harness the large number of gauge degrees of freedom available in the long-range interacting many-mode case.

IV.3 Global diffusion gauge

To make use of Eq. (37) to choose a good gauge, we have to specify the form of the diffusion gauge matrix $\underline{\underline{O}}\,$. In the single mode case Plimak et al. (2001); Deuar and Drummond (2006b) available simulation times for the Gauge-P method can be substantially improved using $\underline{\underline{O}}\,=-\cosh(a)\sigma_{3}-i\sinh(a)\sigma_{1}$, where $a$ is an adjustable real gauge parameter and $\sigma_{i}$ are the standard Pauli matrices. The underlying mechanism is a shift of the numerical noise from the amplitude of $n=\alpha\beta$ to its phase, which less directly leads to deleterious noise amplification. The simplest generalisation to many modes is to use one global value of the parameter $a$, giving the global diffusion gauge

[TABLE]

For a nonuniform system (in a trap, etc.), the obvious adaptation is to have a vector $\bm{a}$ of gauge parameters with spatially varying values $a_{n}=a(\mathbf{x}_{n})$. Such a gauge could most easily be implemented piecewise if the length scale on which the density varies is larger than the characteristic lengthscale of the potential $W$. Then each region would correspond to a block like Eq. (44), optimised using densities inside the block. A local gauge with arbitrarily varying $a(\mathbf{x})$ seriously complicates analysis, because there is an interplay between the lengthscales of the density variation and of the interaction. We will comment on several such more involved choices for the gauge in section VIII.

Returning to the global choice Eq. (44), the quantities appearing in Eq. (37) simplify and are

[TABLE]

IV.3.1 Noise minimisation and gauge choice

Assuming the form Eq. (44), an optimal parameter $a$ can be found by minimizing ${\cal V}(t_{\rm opt})$, where $t_{\rm opt}$ is the chosen target time, typically the longest time for which simulation results are required Deuar and Drummond (2006b). The analytical minimization of Eq. (37) is only tractable with several simplifications. We proceed in similar fashion as was done in Deuar and Drummond (2006b) for a single-mode case. We expand the variance Eq. (37) to second order in $t$:

[TABLE]

where, for the continuum version of the integrals, we define the complex density $\rho_{0}(\mathbf{x})=\rho^{\prime}_{0}(\mathbf{x})+i\rho^{\prime\prime}_{0}(\mathbf{x})=n(\mathbf{x}_{n})/dV$. Here, we have also introduced a semi-positive definite “rectified interaction” matrix

[TABLE]

that plays an important role and will often reappear later. In general, its diagonal elements $U_{0}=U_{nn}$ are not equal to those of $W$ which are $W_{0}$, but they can be related, and are always positive real 555Writing the potential $W$ as an eigenvalue decomposition with orthogonal normalized eigenvectors $|j\rangle$ such that $\langle j|k\rangle=\delta_{jk}$ we have that $W=\sum_{j}\lambda_{j}|j\rangle\langle j|$ and hence $U=\sum_{j}|\lambda_{j}|\,|j\rangle\langle j|$. For cases when the potential $W$ has all eigenvalues positive/negative, one thus finds that $U_{0}=|W_{0}|$. Otherwise, $U_{0}>|W_{0}|$, and the relationship is more complicated. . The potential Eq. (2) used for the examples in section VI, usually has all eigenvalues negative though not always, particularly if periodic boundary conditions occur.

It is instructive to consider the limit of contact interactions, defining their strength $g$ via $U_{0}\to|g|/dV$. The integrals $I_{j}$ in (50) then become

[TABLE]

We see that the $I_{j}$ are extensive quantities, growing with system size, whether measured by mode count $M$ or volume $V$.

We now wish to minimize ${\cal V}(t_{\rm opt})$at a given time $t_{\rm opt}$. Following Deuar and Drummond (2006b), we obtain analytically tractable results in two limits:

Weight dominated,

[TABLE]

when the $I_{1}$ term is small, i.e. $n_{0}$ and $t_{\rm opt}$ are large and direct noise dominated,

[TABLE]

when the $I_{2}$ term is small, i.e. $n_{0}$ or $t_{\rm opt}$ is small. Then we follow the same interpolation procedure that was successful in  Deuar and Drummond (2006b) to give an overall best estimate

[TABLE]

IV.3.2 Adaptive diffusion gauge

Single mode work has shown that the performance can sometimes be improved further, by adapting the gauge parameter $a$ in a time-dependent manner during the simulation, leading to an $a=a(t)$. In that case, one determines $a(t)$ at each timestep according to Eq. (55) using $t_{\rm opt}=t_{\rm fin}-t$, where $t_{\rm fin}$ is the desired final time of the simulation and the values $n_{0}(\mathbf{x})$ appearing in Eq. (55) are changed to the time-dependent $n(\mathbf{x},t)$. We then obtain

[TABLE]

The Stratonovich correction Eq. (96) becomes more involved when $a$ depends on the phase-space variables $\gamma_{\mu}$, details of the derivation are given in A.2.

V Available simulation times

Before attempting to simulate physics with the methods presented here, the crucial question is whether one expects the physical effects of interest to take place prior to materialize before the time ${t}_{\mbox{\scriptsize sim}}$ when the noise level overwhelms the simulation. As in earlier work Deuar and Drummond (2006b), we define this moment via ${\cal V}({t}_{\mbox{\scriptsize sim}})=10$. Based on the optimized gauges from the previous section, we will determine ${t}_{\mbox{\scriptsize sim}}$ for the Gauge-P method in this section.

Gauges do not always offer an advantage over the technically much simpler plain Positive-P method. In some other cases the use of diffusion gauges without drift gauge is preferable. To allow an a-priori assessment, we will also derive ${t}_{\mbox{\scriptsize sim}}$ for those methods. This requires a re-derivation of the characteristic variance for those cases, similar to section IV.2, which is presented in the following.

V.1 Positive-P and diffusion gauged characteristic variance

Let us consider the case without drift gauges ($\underline{\bm{g}}=\underline{\bm{f}}=0$), which includes the plain positive-P representation as a special case ($\underline{\underline{O}}\,=\underline{\underline{I}}\,$), and will also be useful in section VII.2.1. The approximate time evolution of $\cal V$ is calculated similarly to Eq. (37) from Eq. (32), taking $\Omega=1$ and $\underline{\bm{f}}=0$. Details can be found in C. We obtain

[TABLE]

In the above very general case, two initial state covariances appear,

[TABLE]

which are, however, zero in most cases.

The first term of Eq. (57a) is again the direct noise contribution placed into the amplitudes $\underline{\bm{\gamma}}$. The next term involving $\widetilde{\underline{\underline{Q}}\,}$ is the main noise amplification, due to the nonlinearity working on the direct noise. The term involving $\underline{\underline{N}}\,^{\prime}$ is a cross correlation between these, while the remaining terms are due to amplification of any initial fluctuations in $\bm{n}^{\prime\prime}(0)$.

For the case of plain positive-P with $\underline{\underline{O}}\,=\underline{\underline{1}}\,$, Eq. (57) reduces to the following:

[TABLE]

We have used $M\times M$ matrices here.

The initial covariances in (59) are

[TABLE]

The known single mode results  Deuar and Drummond (2006a, b) for repulsive interaction are recovered from Eq. (59) with the replacement $M=1$, $W=(g/dV)\delta_{nm}>0$, and var$[n^{\prime\prime}(0)]=0$. 666 The calculations in Deuar and Drummond (2006a, b) had assumed a deterministic initial state such that $\overline{n(0)}=n(0)$, and thus did not contain the covariance contributions found here.

V.2 Useful simulation time (positive-P)

As has been justified in detail in Deuar and Drummond (2006a), a characteristic variance of ${\cal V}\approx 10$ sets the ultimate limit of usefulness in practice. It gives the variance of the logarithm of a typical observable estimator, and exponentials of Gaussian random variables acquire intractably long distribution tails once this variance reaches a value of ${\cal O}(10)$. The observable estimators are largely of such a form.

Analytical results are only easily obtainable with several simplifications. We proceed in a similar fashion as was done in Deuar and Drummond (2006a) for a single-mode case. First, we assume deterministic initial conditions so that the initial covariances $C^{(0)}$ and $\widetilde{C}^{(0)}$ can be neglected. Then, we expand the variance Eq. (59) to third order in $t$. Unlike for (50), here this is the first order that includes the important noise amplification contribution:

[TABLE]

The rectified interaction $U$ and its diagonal value $U_{0}=U_{jj}\geq 0$ appear again.

To get an idea of the order of magnitude of these quantities, note that in the limit of contact interactions (as before, strength $g$, $U_{0}\to|g|/dV$, and density $\rho_{0}$), we have

[TABLE]

The integrals $I_{j}$ are thus intensive quantities, independent of system volume, unlike those encountered in the corresponding Gauge-P derivation (50). However they are dependent on the numerical lattice spacing $dV$, growing as lattice spacing drops. This is a manifestation of the known feature that simulation time improves as the lattice becomes more coarse.

Consider now the useful case of coherent state (non-entangled) initial conditions, with $n^{\prime\prime}(0)=0$ and no starting variance (${\cal V}(0)=0$). We seek to determine when ${\cal V}^{(P)}=10$ is reached.

The most important regime is when the last term in Eq. (61), that comes from noise amplification, dominates the variance. By using the estimates Eq. (62) to make analogy with the contact-interacting case we can get some bearings. This regime is seen to occur when

[TABLE]

for chemical potentials $\mu\approx g\rho$ that are typical for the superfluid regime of ultracold gases. So, requiring ${\cal V}^{(P)}=10$ in this regime, one readily obtains

[TABLE]

This is a remarkably simple expression. It is roughly independent of the system size, and reduces to the known single-mode result $t_{\sim}\approx 2.5\,dV^{1/3}/[g\rho_{0}^{2/3}]$ in that limit Deuar and Drummond (2006a).

The opposite regime of short times or small chemical potential has evolution still dominated by the direct noise. Then the first term in Eq. (61) dominates and we have

[TABLE]

The situation with dominant $I^{(P)}_{1}$ term does not cover a significant range of parameters.

V.3 Useful simulation time (gauge-P)

With an optimal gauge chosen as discussed in section IV, we can proceed to further estimate the useful simulation time $t_{\rm sim}$ for the gauge-P representation. We follow the same approach as in section V.2 and Deuar and Drummond (2006b). We take again coherent state initial conditions ($n^{\prime\prime}(0)=0$), and ${\cal V}(0)=0$. Then, $I_{1}$ from Eq. (50) vanishes, and we substitute the optimal gauges in the two regimes (53) and (54) into Eq. (50) and impose the condition ${\cal V}(t_{\rm sim})=10$.

In the weight-dominated regime, which applies for large systems or not-too-small densities, one obtains the important result that

[TABLE]

(assuming also $\cosh 2a\approx e^{2a}/2$). This reduces to the results for the single-mode contact interacting case. In particular, for $M$ contact-interacting modes,

[TABLE]

In the direct-noise dominated regime that is relevant for short times or small occupations,

[TABLE]

which roughly applies when it is shorter than Eq. (66), i.e. for $I_{2}\lesssim U_{0}^{2}/40$. Notably, it is the same as Eq. (65) for positive-P. This criterion is usually only met for very small occupations $n_{0}$.

VI Example simulations

We now proceed to demonstrate the utility of these methods with some exemplary applications. The physical models to which we will apply the method have been presented in section II. For both examples, we consider for simplicity a one-dimensional (1D) system with periodic boundary conditions.

VI.1 Extended Bose-Hubbard model

As promised in section II.1, we first consider long-range interactions that are suddenly switched on (interaction quench) in an extended Bose-Hubbard model. This is realistic for engineered interactions in cold atom systems, e.g. with Rydberg dressing. For an initial coherent state, interactions will dephase the different Fock-state components, and coherence between neighbouring sites will be lost. At much longer timescales, the coherence may revive due to the limited number of relevant frequencies involved. Such collapse-and-revival sequences are well understood in the case of local interactions Greiner et al. (2002); Kollath et al. (2007).

Since the expression for the characteristic variance (37) was derived for ${\omega}_{nm}=0$, we shall consider $J=0$ first. We later demonstrate that for small nonzero $J$ the tools derived from (37) still offer useful guidance. For $J=0$, a state $|\bm{n}\rangle=|n_{1},n_{2},n_{3},...\rangle$ with exactly $n_{i}$ bosons on site $i$ is an eigenstate of the Hamiltonian and the time-evolution can be found analytically. Let us consider the case where the initial quantum state is a uniform coherent state on each site such as for a Bose-Einstein condensed gas $|\Psi(t=0)\rangle=|\phi,\phi,...\phi\rangle$, where $\phi\in{\mathbb{C}}$ is the BEC order parameter and $|\phi\rangle=\exp{(-|\phi|^{2}/2)}\sum_{n}\phi^{n}|n\rangle/\sqrt{n!}$. Thus $|\Psi(t=0)\rangle=\sum_{\bm{n}}C_{\bm{n}}|\bm{n}\rangle$, with $C_{\bm{n}}=\exp{(-|\bm{\phi}\cdot\bm{\phi}|/2)}\phi^{\sum n_{i}}/\sqrt{\prod n_{i}!}$. The state approximates the true state of a BEC in a lattice deeply in the superfluid regime. For our demonstration we compare Gauge-P simulations for $J=0$ with the analytical results and then consider various values of $J$ in Gauge-P only.

In Fig. 1 we show the evolution of the inter-site coherence, $g^{(1)}(\Delta n)=\sum_{n}\langle a^{\dagger}_{n}a_{n+\Delta n}\rangle/\sqrt{n_{n}n_{n+\Delta n}}$, the on-site mean field $\langle\hat{a}_{n}\rangle$ and the characteristic variance ${\cal V}$ for an interaction quench on $M=6$ sites spread on a periodic 1D interval of length $L=2$. The potential Eq. (2) is described by $C_{6}=-32$, $\epsilon=1$, $a=2$, $b=3$ and the initial state is $|\phi|^{2}$=$\rho\,dV=1.2$. $J$ is varied over $J=0,2,5,10$.

To know in advance up to which times we can simulate, we determine $I_{2}$ from Eq. (50) using the potential Eq. (2)777$I_{1}$ is zero for our real coherent state initial conditions.. Due to the small number of modes in this example, replacing the integral in $I_{2}$ by the corresponding discrete sum is more accurate, yielding $t<{t}_{\mbox{\scriptsize\rm sim}}=0.14$ using Eq. (66). A comparison with Fig. 1 (b) shows that this indeed describes the available time for the observable $\langle a_{n}\rangle$ quite well, for which the stochastic average is structurally similar to the characteristic variance (36) underlying the estimate for ${t}_{\mbox{\scriptsize sim}}$. In contrast, the sampling of the inter-site coherence $g^{(1)}$ becomes intractable some time before ${t}_{\mbox{\scriptsize\rm sim}}$, around the maximal time shown in panel (a). This is common for higher order quantities. Being interested in $g^{(1)}$, we thus have chosen $t_{\rm opt}=0.05$ in the formula Eq. (55), which then gives $a=0.7$.

Finally, panel (c) numerically verifies Eq. (37). It also shows that for the present model with conserved mean occupation of all the sites, the adaptive diffusion gauge $a(t,n)$ yields no significant improvement of the variance compared to a fixed gauge $a$. The probable reason is the constant mean mode occupation for this example. We expect that there exist cases with time-varying occupation where the adaptive gauge outperforms the fixed one.

Using the ungauged positive-P method we can only get results for a much shorter time shown in Fig. 1 (b). The estimate Eq. (64) $t\lesssim 0.07$ is of the right magnitude, but a little high.

VI.2 Coulomb blockade

We now consider the excitation of atoms within a homogenous 1d BEC to long-range interacting Rydberg states, as was outlined in section II.2. Establishing an interaction blockade within the excited state fraction will correspond to nontrivial correlation function dynamics. To see this within the accessible time interval up to $t_{\rm sim}$, let us consider an echo sequence as used in Raitzsch et al. (2008); Younge and Raithel (2009): Atoms are excited from the ground to the Rydberg state with a Rabi-frequency $\kappa$ during a time-interval $\tau/2$, then de-excited again using $\kappa\rightarrow-\kappa$ for a time $\tau/2$. This yields an additional time-scale $\tau$ which can be chosen such that if it is short enough, it can fit within $t_{\rm sim}$. During this time, we are then interested in the time-evolution of the Rydberg-Rydberg correlation function

[TABLE]

with $\rho_{e}(x)=\langle\Psi^{\dagger}_{e}(x)\Psi_{e}(x)\rangle$. Due to our homogeneous initial conditions, we can additionally average the expression (69) over $x$. The evolution equations are Eq. (34) with the drift gauge Eq. (35) and global diffusion gauge Eq. (44), adapted for two fields and the coupling $\kappa$ in Eq. (4). Namely,

[TABLE]

with $\hbar=m=1$ and kinetic motion, if present, given by $\omega=\frac{1}{2m}\nabla^{2}$. The sign $\mp$ refers to $-$ for $\alpha$ and $+$ for $\beta$ variables.

It would also be interesting to follow the entire development of the blockade-dip for larger times, and many-body Rabi oscillations in the saturated regime. However, we found the time-scale for development of a complete blockade out of reach for all physical and gauge parameters tested here.

VI.2.1 Echo sequence

Without interactions, after completing the echo sequence at time $\tau$, all the atoms would have returned to the ground-state. With interactions, in contrast, the induced dephasing makes this reversal incomplete. We model this scenario with parameters $\tau=0.18$, $\kappa=3$, $C_{6}=-5.96\times 10^{7}$, $\epsilon=12.5$ and initially $N=500$ atoms in state $|\,{g}\,\rangle$ between $x=-L$ and $x=L$ for $L=50$. These parameters are chosen for demonstration, to yield an excitation blockade that only allows a few atoms in state $|\,{e}\,\rangle$ within the domain $[-L,L]$. Fig. 2 shows some of the results. Except where indicated, we invoke the frozen gas approximation, in which the $\omega$ terms in Eq. (70) are neglected.

We find that, counterintuitively, the remnant excited state population after a portion of the de-exitation pulse develops bunching correlations as seen in Fig. 2e. We confirmed this behaviour using the omega-expansion Stanojevic and Côté (2010, 2009), which allows a qualitative calculation of the excited state population $N_{e}$ and correlations $g^{(2)}$ as an expansion in $\omega t$. The leading order results, shown in Fig. 2 for comparison, reveal qualitative agreement with the exact gauge-P with some quantitative differences. Preliminary calculations of the kind shown in Fig. 2 were the first sign of bunching correlations after an echo-sequence discovered by us. This prompted our more detailed investigation using exact diagonalisation for small Rydberg ensembles, published in Wüster et al. (2010b) and experimentally confirmed in Thaicharoen et al. (2016b).

A good a priori estimate of gauge parameters is not so straightforward here like in the last section, because the guiding expressions that we have derived depend on the density of the interacting species $\rho_{e}=N_{e}/2L$, which is strongly time-dependent as seen in Fig. 2(a). First, using $N_{e}\approx 8$ from Fig. 2 (a) as a worst-case estimate, we obtain for the Positive-P method $t_{sim}=0.21$ using Eq. (64), while for Gauge-P we find $t_{sim}=0.32$ using Eq. (66). These are roughly consistent with an empirical simulation time of $t_{\rm sim}=0.2$, but the numerics do not bear out the predicted modest advantage of the gauge-P here. We also note that the optimum fixed gauge found empirically, $a=0.55$, is somewhat smaller than that predicted by Eq. (55), which is $a_{\rm approx}=0.66$ for $t_{\rm opt}=0.135$. These differences do not come as a complete surprise, since our calculations neglect the strong time dependence of $n_{e}$.

While the frozen Rydberg gas dynamics presented so far can be alternatively calculated using exact diagonalisation Wüster et al. (2010b), the quantum mechanical inclusion of motion typically renders that approach intractable, and would require the use of classical or quantum-classical methods Ates et al. (2008); Amthor et al. (2007a, b); Wüster et al. (2010a, 2011b); K Leonhardt and S Wüster and J.-M. Rost (2014); Leonhardt et al. (2017); Thaicharoen et al. (2015); Thaicharoen et al. (2016a); Celistrino Teixeira et al. (2015). Motion however poses no additional challenge to the Gauge-P method, which we now illustrate by dropping the frozen gas approximation and considering atomic motion as outlined in section II.2, with the full equations Eq. (70). Simulated echo-sequences including motion are added to Fig. 2. We used an exemplary mass of $m=4\times 10^{-3}$, chosen to produce a visible effect but not too short ${t}_{\mbox{\scriptsize sim}}$. As one can see, the change is significant.

There are many scenarios in which the frozen Rydberg gas approximation is not applicable such asAtes et al. (2008); Amthor et al. (2007a, b); Wüster et al. (2010a, 2011b); K Leonhardt and S Wüster and J.-M. Rost (2014); Leonhardt et al. (2017); Thaicharoen et al. (2015, 2015); Celistrino Teixeira et al. (2015). Other possible realisations of our long-range interacting model with relevant atomic motion could be a two component polar BEC, where the long-range interactions between one of these components are dramatically enhanced with a Feshbach resonance Griesmaier et al. (2005); Lahaye et al. (2007) or a Rydberg dressed BEC Gil et al. (2014); Macrì and Pohl (2014); Li et al. (2013); Zeiher et al. (2016).

VII Methods for very large systems

The sets of stochastic differential equations Eq. (13) and Eq. (32a) are usually quite straightforward to implement for small and medium size systems up to $M\sim\mathcal{O}(10-1000)$. However, for larger system two problems arise: The direct numerical implementation of the equations as written may become too inefficient, and the drift gauge may cease to be beneficial. We discuss both of these issues, and possible solutions, below.

VII.1 A more efficient treatment of the interaction for large systems

For larger systems, the direct application of Eq. (13) and Eq. (32a) by summing over $m=1,\dots,M$ becomes a problem for two reasons: (i) integrating a single time step takes $\mathcal{O}(M^{2})$ operations (a separate summation of $\sum_{m}W_{nm}\alpha_{m}\beta_{m}$ and $\sum_{j}[\sqrt{W}]_{nj}\xi_{j}$ for each lattice site $n$), and (ii) if $M$ gets really large, e.g., of the order of $M\gtrsim\mathcal{O}(10^{5})$ in a three-dimensional system, calculating the $M\times M$ matrix $\sqrt{W}$ at the beginning of a calculation becomes time and memory consuming.

However, the computational work could be substantially reduced by taking advantage of the fact that physical potentials $W(\mathbf{z})$ depend only on the vector difference $\mathbf{z}=\mathbf{x}-\mathbf{y}$. The overwhelming part of information in the two-index matrix $W_{nm}$ is redundant: $W_{nm}=W_{(n-m)\,{\rm mod}\,M,0}=W_{q0}$ with $q=(n-m)\,{\rm mod}\,M$, an index that embodies only the relative displacement on the numerical lattice.

Consider the discrete Fourier transform (DFT) of this translated interaction potential $W_{q0}=W(\mathbf{z}_{q})$ on a numerical lattice with periodic boundary conditions:

[TABLE]

The last expression brings us back to the two-coordinate form $W_{nm}$ for all values of $m$. Here, $\mathbf{k}_{\widetilde{m}}$ are the standard wavevectors on the numerical lattice. Since the potential $W_{nm}$ is even, that is $W_{n,(n+d)\,{\rm mod}\,M}=W_{n,(n-d)\,{\rm mod}\,M}$ for any $d$ due to 1D periodicity (analogously in higher dimensions), the transform $\widetilde{W}$ will be real. The interaction drift term for $d\alpha_{n}/dt$ in Eq. (13) can then be written as an inverse DFT

[TABLE]

where $V=MdV$ is the total volume, tildes will indicate k-space, and

[TABLE]

is the DFT of the density $n_{p}(t)=\alpha_{p}(t)\beta_{p}(t)$ for a given trajectory at a given time. This way, we can see that the deterministic evolution can be carried out by storing the interaction vector $\widetilde{W}_{\widetilde{m}}$ of just size $M$ from the beginning of a simulation, and by carrying out two DFTs per trajectory and time step to obtain $\widetilde{n}$. This has only $M\log M$ computational cost with “fast Fourier transform” methods.

Remarkably, the noise terms can be dealt with in a related, though more involved way, that follows an approach reported in Deuar (2005a). Denoting the noise term for $d\alpha_{n}/dt$ as $X_{n}=\sqrt{-i}\,\alpha_{n}\sum_{j}[\sqrt{W}]_{nj}\xi_{j}^{(1)}(t)$, the diffusion condition Eq. (12) requires simply that

[TABLE]

(equal times will be assumed in this section). Then, if we rewrite the noise $X_{n}$ in terms of new scaled and Fourier-transformed noise quantities $\widetilde{Y}$ as

[TABLE]

it can be shown by simple substitution that the condition

[TABLE]

on the new noises is sufficient to satisfy the required diffusion condition Eq. (74). If we are able to construct a noise field $\chi_{\widetilde{p}}$ with the somewhat atypical correlation properties

[TABLE]

then the choice $\widetilde{Y}_{\widetilde{p}}=\chi_{\widetilde{p}}\,\sqrt{V\,\widetilde{W}_{\widetilde{p}}}$ will satisfy all that is required. Such a noise can indeed be constructed, as is described in D.

With it, the noise term Eq. (75) can be written as an appropriate inverse DFT:

[TABLE]

This is also implementable with computational cost $M\log M$ per time step. The linear term $\sum_{m}\omega_{nm}\alpha_{m}$ does not usually constitute an efficiency problem, so we leave it as is. Applying the above, using independent noise fields $\chi^{(1)}$ and $\chi^{(2)}$, the final more efficient equations in the positive-P representation are

[TABLE]

VII.2 Drift gauges and system size

Looking at the simulation time estimate for the gauge-P distribution in the contact-interacting case Eq. (67), one sees that after intensive quantities like $g$ and $\rho$ have been factored out, a disadvantageous reduction with the system size $M$ remains. This is a feature that has been also noted previously. It arises because the single weight variable collects noise from the entire systemDeuar and Drummond (2006b), while the limit ${\cal V}\lesssim 10$ remains unforgiving as system size grows. In contrast, positive-P or only diffusion gauged calculations do not show this behaviour. For this reason, though they are consistently less advantageous for single modes, they may be better when the system size becomes large.

VII.2.1 Diffusion gauge only

To address these issues, we first evaluate how a global diffusion gauge Eq. (44) performs in long-range interacting systems without any drift gauges. We proceed the same way as in section IV.3.1 and V.3 and Deuar and Drummond (2006b). The variance Eq. (57) is expanded to third order in $t$, with gauge Eq. (44), and deterministic initial conditions $C^{(0)}=\widetilde{C}^{(0)}=0$. The result is:

[TABLE]

with the same integrals as in Eq. (61). The $I_{1}^{(P)}$ term is irrelevant as it does not depend on $a$. The optimum gauge then is readily shown to be simply

[TABLE]

in agreement with known single-mode results. From the estimates Eq. (62), we see that this is largely independent of $M$, as hoped.

The simulation time, requiring ${\cal V}(t_{\rm sim})=10$ can be found in the usual two regimes: The noise amplification dominated regime, when $n$ and $t_{\rm sim}$ are large enough, has

[TABLE]

This is a different and more advantageous scaling as compared to the plain positive-P Eq. (64). The diffusion gauged value of $t_{\rm sim}$ is longer when

[TABLE]

Notably, $I_{2}^{(P)}$ is an integral involving third powers of the potential, so Eq. (83) is primarily a condition on the density, requiring it to be sufficiently high. Substituting the estimates Eq. (62), one has $n_{0}\gtrsim 0.2$ for the lattice occupations. This result does not depend on system size. The direct-noise dominated regime has the same estimate Eq. (65) as in the other cases.

VII.2.2 Diffusion gauges vs drift gauges

With the different $t_{\rm sim}$ of drift and diffusion-only gauged calculations, a pertinent question is when should the drift gauge be added to the diffusion gauge? Comparing Eq. (66) and Eq. (82), diffusion only calculations last longer when

[TABLE]

For contact interactions in a uniform system when substituting Eq. (52) and Eq. (62) this reduces to the surprisingly simple expression

[TABLE]

For long range interactions, condition Eq. (84) should be used instead of Eq. (85). Also there it is clear, however, that for large enough systems, say $M\gg 100$, the diffusion gauge only approach will be better.

The simulations of section VI.2 were an example of a sufficiently large system, with $M=64$, $I_{2}^{(P)}\approx 2500$, $U_{0}=15.624$, and $I_{2}\approx 1500$, where drift gauges cease to offer a huge advantage.

VIII General many-mode diffusion gauges

We have seen so far, that the global gauge ansatz (44) can already provide significant advantages as well as technical challenges in the implementation. However, it only exploits a small fraction of the degrees of freedom for the many-mody case that are inherent in the full $2M\times 2M$ matrix $\underline{\underline{O}}\,$. In this section we present our initial explorations of non-local many-mode Gauges. These are intended as starting point for future research, beyond the scope of this article.

For many-mode problems with contact potentials $W=W_{0}\mathds{1}$ as considered in previous work, a local diffusion gauge

[TABLE]

is a feasible choice that already goes beyond (44). Since $\underline{\underline{O}}\,$ and $\underline{\underline{S}}\,$ in Eq. (34b) commute in this case, each gauge parameter $a_{m}$ can be individually associated with the mode $m$. For sufficiently small inter-mode coupling, one can then use an optimized gauge $a_{m}(t)$ based on single-mode results that depend on the stochastic occupation $n_{m}$ of this particular mode $m$. Often this gives good results Deuar (2004).

However, for a non-local potential $W$, the matrices $\underline{\underline{O}}\,$ and $\underline{\underline{S}}\,$ no longer commute and hence the $a_{n}$ would have to be viewed as gauge parameter of the noise source $\xi_{n}$. In that case, it is not clear to what degree a naive application of single-mode results is still possible. In the remainder of this section we describe two approaches aimed at harnessing more of the gauge freedom of $\underline{\underline{O}}\,$: an analytical one and a numerical one.

VIII.1 Analytical nonlocal diffusion gauges

Finding the matrix $\underline{\underline{O}}\,$ which globally minimizes the variance for a given interaction potential and simulation time is a difficult task, analytically and numerically. Here, we investigate the ansatz

[TABLE]

where the hyperbolic functions are defined in terms of matrix exponentials, and $A$ is a symmetric $M\times M$ matrix. The matrix $\underline{\underline{O}}\,$ is orthogonal for all choices of $A$ and reduces to (44) for $A=a\mathds{1}$, with scalar $a$.

Now let us consider the subset of nonlocal interactions for which $W=U$, i.e. $\sqrt{W}=\sqrt{W}^{\dagger}=\mbox{Re}[\sqrt{W}]$. We take $A$ real and the starting densities deterministic: $\bm{n}(0)=\overline{\bm{n}(0)}$. Then, inserting $\underline{\underline{O}}\,$ into the variance expression Eq. (37) and simplifying, we obtain the following for short times:

[TABLE]

We used the matrices

[TABLE]

A minimum of the scalar quantity ${\cal V}$ in the space of matrices $A$ is given by the solution of $\partial{\cal V}/\partial A_{ij}=0$.

As described in more detail in E, if $\bm{n}$ is homogeneous and real, we can obtain the matrix equation

[TABLE]

Equation Eq. (94) has the same structure as the single mode version in Deuar and Drummond (2006b). Hence, following Deuar and Drummond (2006b), we can solve this equation separately for cases where the $\exp(-4A)$ term or $\exp(-6A)$ term on the right hand side is dominant, and then suggest a simple interpolation equation:

[TABLE]

where $\log$ denotes a matrix logarithm. We present Eq. (95) for completeness, as for the cases explored, it did not prove more useful than the much simpler adaptive global gauge of section IV.3.

We defer an initial test of the performance of (95) to the next section.

VIII.2 Numerical nonlocal diffusion gauges

Another route by which non-local diffusion gauges might offer an advantage over local adaptive ones, is the numerical minimization of Eq. (37) with respect to an arbitrary orthogonal matrix $O$ for given initial mode occupations $n_{m}(0)$, $t$ and ${W}$. We carried out such a minimization for the same potential as in section VI.2.1 with spatial atom-density $n(x)=\bar{n}\exp[-x^{2}/2/\sigma^{2}]$, $\sigma=10$, $\bar{n}=0.5$ and $t_{\rm opt}=0.15$. A small number of modes $M=16$ covering the range $x\in[-L,L]$ with $L=50$ allowed reasonably fast optimization. In practice the degree of freedom for the conjugate gradient optimization routine Blakie (2001) was a real, anti-symmetric matrix $g$ such that $O=\exp(g)$ as in Eq. (44).

We employed a variety of different initial conditions: (i) local global gauge (44) with $a$ guessed, (ii) non-local gauge (91) with $A_{kl}=$const, (iii) $\underline{\underline{O}}\,=\exp(g)$ with random $g_{kl}$, (iv) $g=0$, (v) the analytical solution Eq. (95). For comparison with local diffusion gauges Eq. (88), we extracted from the optimized $\underline{\underline{O}}\,$ a local diffusion gauge parameter $a(x_{\mu})=\sinh^{-1}[\mbox{Im}(\underline{\underline{O}}\,_{M+\mu,\mu})]$ for $1\leq\mu\leq M$, as well as the non-local shape $a(r_{\mu})=\sinh^{-1}[\mbox{Im}(\underline{\underline{O}}\,_{M,\mu})]$, both shown in Fig. 3 for the initial condition (i). For comparison, a pure local diffusion form Eq. (88) would contain nonzero elements only exactly on the dominant diagonal features described by $g(x)$ as seen in Fig. 3 (b–c), while the magenta line for $g(r)$ in Fig. 3 (a) would only have the one nonzero element at $r=0$.

All initial conditions led to quite similar results as shown, with a significant spatial dependence of the gauge parameter $a(x)$ due to the inhomogeneous density $n(x)$, and small non-local tails visible in $a(r)$. The overall decrease in variance of Eq. (37) compared to a local diffusion gauge was by a factor of $2.3$ in the case shown, with weak non-locality. We also used the numerical optimisation scheme to confirm that (95) is indeed an optimal choice if $n(x)$ is uniform and times are short. Even for the inhomogeneous test case described above it already outperforms the global gauge by a factor of $2.1$, close to the performance of the completely optimized result. The expression (95) also roughly captures the off-diagonal shape of gauge parameters found by the numerical optimisation in those cases.

From this reduction of the characteristic variance ${\cal V}$, one can make one important general conclusion – that there is indeed some potential for gains from a nonlocal optimization of the gauge. Since we just skimmed the surface of gauges more general than Eq. (44), this indicates an extensive future route for possible improvement in our simulation methods, beyond the scope of the present article.

IX Summary

We have provided the necessary apparatus for carrying out stochastic positive-P and Gauge-P simulations for systems of long-range interacting bosons. Our focus has been on directions relevant for ultracold gases such as Rydberg and dipolar atoms, as this is presently a burgeoning field of experimental and theoretical investigation. A summary of our main analytical results is given in table 1 for easy reference.

Our most important results are the dynamical equations with a long-range potential, the estimates for the available simulation times $t_{\rm sim}$, and simple formulae for optimal diffusion gauge parameters $a_{\rm opt}$. The simulation times reveal whether a particular problem is amenable to calculation using this method, while the gauge formulae give a straightforward means of improving simulation times beyond those directly available from the positive-P method. The underlying general expressions for the noise evolution (the characteristic variances ${\cal V}$) constitute a good starting point for future investigations, including the development of truly non-local diffusion gauges. The results of Sec. VIII.2 confirm that nonlocal gauges can offer advantages in principle.

Comparing with the previously known results for contact interactions, we can make several statements: stochastic gauges continue to improve simulation times. For this we use a global gauge that is optimised with the help of nonlocal integrals such as $I_{2}$ and $I_{2}^{(P)}$. The resulting expressions for the optimal gauge, simulation time, and characteristic variance that we have derived bear structural similarity to the contact interacting case, and reduce to them in the appropriate limit. For large systems with long-range interaction, special care must be taken to avoid a quadratic scaling with system size and memory problems, but remarkably, this can be achieved using Fourier transformed interactions and noise fields. This technique makes simulations with, say, $M=10^{5}-10^{6}$ possible, which would otherwise be well out of reach.

Furthermore, we have confirmed here a suspicion that arose already for contact-interacting gases, that very large systems with $\gtrsim{\cal O}(100)$ modes have different stochastic gauge properties than small to medium systems. Drift gauges become inefficient for large systems because of the accumulation of noise from all modes into the one global weight $\Omega$. In that case, diffusion gauges should be used, and continue to provide improvement over the plain positive-P treatment for many cases.

Finally, we have illustrated our techniques by modeling an interaction quench in the extended Bose-Hubbard model and the excitation of Rydberg states in a Bose-Einstein condensate, using an echo sequence. Notably, the inclusion of atomic motion in such a calculation introduces no additional significant complications here, in contrast to many other methods used to study Rydberg excitations.

Acknowledgements

It is a pleasure to thank André Ribeiro de Carvalho, Karen Kheruntsyan and Murray Olsen for comments. We also gratefully acknowledge useful discussions with Cenap Ates, Thomas Pohl and Jovica Stanojevic. P.D. acknowledges support from the National Science Centre (Poland) grant No. 2012/07/E/ST2/01389.

Appendix A Stratonovich corrections

The evolution equations Eq. (34a)–Eq. (34b) are to be interpreted in Itô stochastic calculus. Most numerical implementations that go beyond simple and inefficient first-order methods require the use of the Stratonovich form. This includes the adaptive 8th/9th order Runge-Kutta method Wilkie and Çetinbaş (2005) employed for all simulations in this article using the XMDS package used here Collecutt and Drummond (2001) (see XMD ; Dennis et al. (2013) for the most recent version). For equations with variable-dependent noise such as the ones derived here, the Stratonovich correction can be complicated, especially when the gauges themselves are a function of the variables. Hence, it is of utility to provide the end results here.

If one decomposes the variables in $\overline{\bm{v}}=(\bm{\gamma},\Omega)$ into all the independent variables (either ${\rm Re}[v]_{\mu}$ and ${\rm Im}[v]_{\mu}$ or $v$ and $v^{*}$), labeling them them as $u_{j}$ with $j=1,\dots 2(2M+1)$, then one can apply the usual formula Gardiner and Zoller (2004); Gardiner (2003) for the Stratonovich correction:

[TABLE]

for a real random process with Itô equation

[TABLE]

The Stratonovich form of the equation is then $\frac{du_{j}}{dt}=A_{j}(\overline{\bm{v}})+\sum_{\mu}B_{j\mu}(\overline{\bm{v}})\xi_{\mu}+S_{u_{j}}(\overline{\bm{v}})$.

A.1 Uniform diffusion gauge

With a diffusion gauge $\underline{\underline{O}}\,$ independent of the variables, but allowing arbitrary drift gauges $\underline{\bm{f}}$, we obtain:

[TABLE]

The notation $\left[\underline{\bm{a}}\right]_{\nu}$ is used for the $\nu$-component of the vector inside the brackets.

If we specify the drift gauge as $\underline{\bm{f}}=i\underline{\bm{n}}^{\prime\prime}$ from Eq. (35) as we have used throughout this paper, one obtains

[TABLE]

with $\underline{\underline{F}}\,$ defined in Eq. (38). A case relevant for this paper is the uniform diffusion gauge Eq. (88) using a real $\bm{a}$, where one has:

[TABLE]

which reduces to

[TABLE]

for the global gauge Eq. (44). The real rectified potential Eq. (51) appears again. The case of contact interactions is recovered when $W_{nm}=g\delta_{nm}/dV$. An ungauged positive-P simulation, or one with diffusion gauges that do not depend on the variables has only the Eq. (98) corrections.

A.2 Adaptive gauge Stratonovich correction

In the case of an adaptive diffusion gauge the Stratonovich correction becomes significantly more complicated, since the derivatives in Eq. (96) now also act on the gauge parameter $a$. Taking the global form Eq. (44) but with the variable and time dependence of Eq. (56), after much algebra, one obtains:

[TABLE]

[TABLE]

Appendix B Characteristic variance: with drift gauge

B.1 Logarithmic variables

To estimate Eq. (36), a first step is to turn the stochastic equations (differentials) Eq. (34a)-Eq. (34b) into equations for $\log{|\Omega\gamma_{\mu}|}$. We use the usual multivariate Itô formula Gardiner (2003); Arnold (1992) using the independent variables $\overline{\bm{v}}$ like in Eq. (97), to give

[TABLE]

From this point on, we will neglect kinetic terms, setting $\underline{\underline{\omega}}\,=0$, as explained in section V.1. Formally, Eq. (105) can be integrated over time, with the result

[TABLE]

Next, we move from $\log{(\Omega\gamma_{\mu})}$ to $\log{(|\Omega\gamma_{\mu}|)}$. It is helpful to specify the form of $\bm{f}$ at this point, as the usual $\underline{\bm{f}}=i\underline{\bm{n}}^{\prime\prime}$. This has the convenient feature that the deterministic parts of $\log{(\Omega\gamma_{\mu})}$ are imaginary and do not contribute to $\log{|\Omega\gamma_{\mu}|}$ since $\log|z|={\rm Re}[\log z]$. This is likely an important factor why this is a useful gauge. With this choice of $\underline{\bm{f}}$, we arrive at

[TABLE]

As previously, the indices below square brackets denote the elements of the enclosed matrix/vector.

B.2 Variance assembly

The mean of Eq. (107) is simply $\overline{\log{[|\Omega(t)\gamma_{\mu}(t)|]}}=\overline{\log{[|\Omega(0)\gamma_{\mu}(0)|]}}$, since $\overline{\int_{0}^{t}G(s)\xi_{\sigma}(s)}=0$ for any non-anticipating function $G$ Gardiner (2003). The remaining averaging at $t=0$ is over initial state noise, which is independent of any later dynamical noises $\xi$. The mean of the square is

[TABLE]

We consider the three terms in sequence:

Term (108a): From independence of noises, $\overline{\xi_{\nu}(s)\xi_{\nu^{\prime}}(s^{\prime})}=\delta_{\nu\nu^{\prime}}\,\delta(s-s^{\prime})$ so that

[TABLE]

Hence

[TABLE]

For the second line we note that $\underline{\underline{S}}\,\underline{\underline{S}}\,^{T}=\underline{\underline{W}}\,$ and $\underline{\underline{W}}\,$ is real, as well as that $(\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger})^{*}=(\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger})^{T}$. The transpose is then irrelevant for the $[\ ]_{\mu\mu}$ element. Note that while $\underline{\underline{O}}\,\underline{\underline{O}}\,^{T}=\underline{\underline{\mathds{1}}}\,$, this is not generally the case for $\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}$.

Term (108b): The stochastic average is $\overline{\,n^{\prime\prime}_{\lambda}(s)n^{\prime\prime}_{\lambda^{\prime}}(s^{\prime})\xi_{\nu}(s)\xi_{\nu^{\prime}}(s^{\prime})\,}$. Now consider the following: if $s>s^{\prime}$ then then $\xi(s)$ is uncorrelated with all the other terms, because of causality, and for the case of $n^{\prime\prime}(s)$, the non-anticipating nature of an Itô process. Hence, in this case we can write

[TABLE]

since $\overline{\xi_{\nu}(s)}=0$. The same argument can be made for $s^{\prime}>s$, in which case a factor $\overline{\xi_{\nu^{\prime}}(s^{\prime})}=0$ occurs. The remaining situation is $s=s^{\prime}$, in which case the noises $\xi$ are also uncorrelated with the $n^{\prime\prime}$ because of the non-anticipating nature. Hence, we can write

[TABLE]

The shorthand expression

[TABLE]

for the correlation matrix of integrated expectation values will be useful. Inserting Eq. (115) into Eq. (108b) gives

[TABLE]

Term (108c): Proceeding in the same way, one finds

[TABLE]

and

[TABLE]

Gathering now all the terms and inserting them into Eq. (108) and Eq. (36) we obtain:

[TABLE]

B.3 Stochastic density evolution

To proceed further we must evaluate the expectation values involving $\underline{\bm{n}}^{\prime\prime}$. Applying the Itô formula for $n_{m}=\alpha_{m}\beta_{m}$, and again discarding the $\omega$ terms, leads to

[TABLE]

Using the matrix $\underline{\underline{F}}\,$ defined in Eq. (100), this can be written

[TABLE]

The last follows because $\underline{\underline{W}}\,(\underline{\underline{F}}\,-\underline{\underline{I}}\,)$ has all zeros on the diagonal. Moving to logarithmic variables and formally integrating the differential we obtain

[TABLE]

where

[TABLE]

Since it is a sum of Gaussian random variables, $\zeta_{\mu}$ is itself Gaussian distributed. It has correlations

[TABLE]

as found inDeuar (2004).

B.4 Stochastic density means and variances

The expectation value of Eq. (123) can be written in terms of the probability distribution of Gaussian random variables $\zeta_{\nu}$ that obey Eq. (125):

[TABLE]

Here, $\underline{\bm{K}}^{(\mu)}$ are vectors with elements $K_{\nu}^{(\mu)}=\sqrt{i}\left[\underline{\underline{F}}\,\underline{\underline{S}}\,\underline{\underline{O}}\,\right]_{\mu\nu}$, and the probability distribution is

[TABLE]

We can then write the expectation value as an integral over the quadratic form

[TABLE]

where $\underline{\underline{A}}\,=\underline{\underline{\mathds{1}}}\,/2t$, whence det$\underline{\underline{A}}\,=(2t)^{-2M}$. Using the standard result that

[TABLE]

the expectation value finally becomes

[TABLE]

Note that for our form of matrices $W$ and $F$ in Eq. (25) and Eq. (100),

[TABLE]

Hence $\overline{n_{\mu}(t)}=\overline{n_{\mu}(0)}$ and

[TABLE]

We still have to evaluate $Q$ of Eq. (116). Note that for $z,w\in{\mathbb{C}}$ we have $z^{\prime\prime}w^{\prime\prime}={\rm Re}(zw^{*}-zw)/2$. Thus

[TABLE]

Using Eq. (123) we can write:

[TABLE]

Now the symbol $\underline{\bm{K}}^{(\lambda,\lambda^{\prime})}$ is

[TABLE]

To use Eq. (129), we obtain

[TABLE]

We have used $\underline{\underline{F}}\,\underline{\underline{W}}\,\underline{\underline{F}}\,^{T}=0$ and the other obvious properties of $\underline{\underline{F}}\,$, $\underline{\underline{S}}\,$, $\underline{\underline{O}}\,$. Thus

[TABLE]

Using the same methods we obtain

[TABLE]

Now we have to integrate over time. The result is Eq. (37b). We also make use of Eq. (134) and the proprtties of $\underline{\underline{W}}\,$ and $\underline{\bm{n}}$ to find that

[TABLE]

Finally, substituting Eq. (141) and Eq. (37b) into Eq. (120), we obtain the pivotal result Eq. (37).

Appendix C Characteristic variance: without drift gauge

The case when there is no drift gauge ($\underline{\bm{f}}=0$), follows in the same general way as above in appendix B, though with some additional elements. Summarising the procedure, we begin with

[TABLE]

With the usual assumptions, $\underline{\bm{n}}(t)$ continues to be described by Eq. (123), and $\overline{n_{\mu}}=n_{\mu}(0)$. With the use of this, and algebra

[TABLE]

in terms of a quantity

[TABLE]

that is similar to $Q$.

To evaluate Eq. (143) and Eq. (144) the two-time averages $\overline{n^{\prime\prime}_{\lambda}(s)n^{\prime\prime}_{\lambda^{\prime}}(s^{\prime})}$ and $\overline{\zeta_{\lambda}(t)n^{\prime\prime}_{\lambda^{\prime}}(s)}$ are needed. From the construction Eq. (124) we note that when $s^{\prime}>s$, the random variable $\zeta$ can be split into completely independent early and late ($\widehat{\zeta}$) parts

[TABLE]

with $\widehat{\zeta}(v)$ having the same variances as $\zeta(v)$ given inEq. (125), except that it is independent because it consists of different building block noises $\xi(t)$. This way, after using Eq. (123) and Eq. (145) one has

[TABLE]

The factor in the right-hand integral can be identified as a scaled expectation value from Eq. (130), $\overline{n_{\lambda^{\prime}}(s^{\prime}-s)}^{*}/\overline{n^{*}_{\lambda^{\prime}}(0)}$, and substituting $s^{\prime}=v+s$ one has

[TABLE]

Here, Eq. (134) was used in the 2nd line, and Eq. (139) in the third. Similarly,

[TABLE]

and $\int_{0}^{t}ds\,\int_{0}^{t}ds^{\prime}\,\overline{n^{\prime\prime}_{\lambda}(s)}\overline{n^{\prime\prime}_{\lambda^{\prime}}(s^{\prime})}=t^{2}\,\overline{n^{\prime\prime}_{\lambda}(0)}\overline{n^{\prime\prime}_{\lambda^{\prime}}(0)}$ from Eq. (134), so that using a form similar to Eq. (135) one has

[TABLE]

This can be rearranged to the form

[TABLE]

using $\widetilde{Q}$ from Eq. (57b). Consider now the quantity $\overline{\zeta_{\lambda}(t)n^{\prime\prime}_{\lambda^{\prime}}(s)}$ in the last line of Eq. (143). Splitting $\zeta$ into early and late parts as $\zeta(t)=\zeta(s)+\widehat{\zeta}(t-s)$ like in Eq. (145), one an see that $\widehat{\zeta}(t-s)$ must be uncorrelated with the earlier $n^{\prime\prime}(s)$. From this, and Eq. (123) one has

[TABLE]

The right-most expectation value is a Gaussian integral similar to Eq. (128), but the first moment of the random variable $\zeta(s)$. The matrix $\underline{\underline{A}}\,$ is now $\underline{\underline{1}}\,/(2s)$ and we have $\underline{\bm{K}}^{(\nu^{\prime})}_{\sigma}=\sqrt{i}\left[\underline{\underline{F}}\,\underline{\underline{S}}\,\underline{\underline{O}}\,\right]_{\nu^{\prime}\sigma}$. Here, the result from multivariate normal distributions that we need is

[TABLE]

With this one finds

[TABLE]

where $\underline{\underline{N}}\,={\rm diag}\left[\langle\underline{\bm{n}}(0)\rangle\right]$ as defined before. The bottom line in Eq. (143) is

[TABLE]

Finally substituting this and Eq. (150), the expression Eq. (57) is obtained after some manipulation to separate the mostly useless initial covariance $C^{(0)}$ defined in Eq. (58).

Appendix D Noise performance improvement for large systems

As described in section VII.1, a more efficient formulation of the Gauge-P SDEs is possible using a noise field $\chi_{\widetilde{p}}$ with correlation properties

[TABLE]

Such a noise can indeed be constructed, by first dividing the k-space of the problem into three regions: $\mathcal{R}$ contains half of all the k-space vectors $\mathbf{k}$ that possess a matching wavevector $-\mathbf{k}$ reflected across the origin. In 3D on a standard DFT lattice $\mathcal{R}$ can consist for example of all wavevectors in the $k_{x}>0$ half-space, plus all vectors on the $k_{x}=0$ plane with $k_{y}>0$, plus all vectors on the $k_{z}$ axis with $k_{z}>0$. $\mathcal{R}^{\prime}$ contains all the negatives $-\mathbf{k}$ of the wavevectors $\mathbf{k}\in\mathcal{R}$. Finally, $\mathcal{R}_{0}$ contains the residual unpaired wavevectors. We always have $[\mathbf{k}=0]\in\mathcal{R}_{0}$, but if the number of lattice points $M_{d}$ in a particular dimension $d$ is even, it will also contain the unpaired wavevectors possessing maximally negative values $k_{d}=-\pi L_{d}/M_{d}$ that occur in the DFT for box length $L_{d}$. With this division, the prescription for Eq. (155) is

[TABLE]

The real noises $\xi$ have the same properties as before, i.e., they are all independent with variance $\overline{\xi_{(\widetilde{p},j)}(t)\xi_{(\widetilde{q},j^{\prime})}(t^{\prime})}=\delta_{\widetilde{p}\widetilde{q}}\,\delta_{jj^{\prime}}\,\delta(t-t^{\prime})$.

Appendix E Non-local gauge expressions

The purpose of this appendix is to elaborate how the results of section VIII.1 were obtained. We desire a simpler version of Eq. (37) valid for short times. To this end, the exponential within each element $\lambda,\lambda^{\prime}$ of the matrix $\left(e^{t\left[\underline{\underline{F}}\,\underline{\underline{S}}\,\underline{\underline{O}}\,\underline{\underline{O}}\,^{\dagger}\underline{\underline{S}}\,^{\dagger}\underline{\underline{F}}\,\right]_{\lambda\lambda^{\prime}}}-1\right)$ in Eq. (37) is expanded to second order in $t$. One obtains

[TABLE]

as written in Eq. (92), with the diagonal matrix $N$ given by Eq. (93).

The more difficult task, is to the differentiate Eq. (157) with respect to each element $a_{ij}$ of the diffusion gauge matrix $A$. We make use of the expression for derivatives of matrix exponentials:

[TABLE]

where $[\partial A/\partial a_{ij}]_{lm}=\delta_{li}\delta_{mj}$. With this expression we require:

[TABLE]

Let us assume that $A$ will be a function of $\sqrt{W}$ only888Defined by a power series., and the initial mode occupation is constant and real $n_{0}(x)=n_{0}$. Now matrices such as $\exp{[2\tau A]}$, $N^{\prime}$ and $\sqrt{W}$ commute. The $\tau$ integration is then trivial, and we arrive at the matrix equation:

[TABLE]

from which we can directly write Eq. (94) and then Eq. (95).