Continuum of classical-field ensembles from canonical to grand canonical and the onset of their equivalence

TL;DR

This paper investigates how interactions in ultracold Bose gases lead to the convergence of canonical and grand-canonical ensembles, revealing robust features and ensemble equivalence at strong interactions in mesoscopic 1D systems.

Contribution

It introduces a method using a modified stochastic Gross-Pitaevskii equation to generate canonical and intermediate classical field ensembles, exploring ensemble equivalence with interactions.

Findings

Ensemble distributions become identical at sufficient interaction strength.

Fragile features in ideal gases become robust with strong interactions.

The modified SGPE effectively generates various classical field ensembles.

Abstract

The canonical and grand-canonical ensembles are two usual marginal cases for ultracold Bose gases, but real collections of experimental runs commonly have intermediate properties. Here we study the continuum of intermediate cases, and look into the appearance of ensemble equivalence as interaction rises for mesoscopic 1d systems. We demonstrate how at sufficient interaction strength the distributions of condensate and excited atoms become practically identical regardless of the ensemble used. Importantly, we find that features that are fragile in the ideal gas and appear only in a strict canonical ensemble can become robust in all ensembles when interactions become strong. As evidence, the steep cliff in the distribution of the number of excited atoms is preserved. To make this study, a straightforward approach for generating canonical and intermediate classical field ensembles using a modified stochastic Gross-Pitaevskii equation (SGPE) is developed.

Equations

$$\hat{H}=\int d^d\mathbf{x}{\hat{\Psi }}^{†}(\mathbf{x})\left[H^{\mathrm{sp}}+\frac{g}{2}{\hat{\Psi }}^{†}(\mathbf{x})\hat{\Psi }(\mathbf{x})\right]\hat{\Psi }(\mathbf{x})$$

$$H^{\mathrm{sp}}=-\frac{{\hslash }^2}{2m}{\nabla }^2+V(\mathbf{x}).$$

$$\hat{H}\to dV\sum _{\mathbf{x}}{\hat{\Psi }}_{\mathbf{x}}^{†}\left[\sum _{\mathbf{y}}{H}_{\mathbf{x}\mathbf{y}}^{\mathrm{sp}}{\hat{\Psi }}_{\mathbf{y}}+\frac{g}{2}{\hat{\Psi }}_{\mathbf{x}}^{†}{\hat{\Psi }}_{\mathbf{x}}{\hat{\Psi }}_{\mathbf{x}}\right].$$

$${\left(H_{\mathbf{x}\mathbf{y}}^{\mathrm{sp}}\right)}^{\ast }=H_{\mathbf{y}\mathbf{x}}^{\mathrm{sp}}.$$

$$\hat{\Psi }(\mathbf{x})=\sum _j{\varphi }_j(\mathbf{x}){\hat{a}}_j.$$

$$\hat{\Psi }(\mathbf{x})\to \left\{\psi (\mathbf{x})=\sum _{j\in \mathcal{C}}{\varphi }_j(\mathbf{x}){\alpha }_j\right\},$$

$$E(\psi )=\int d^d\mathbf{x}\psi (\mathbf{x}{)}^{\ast }\left[H^{\mathrm{sp}}+\frac{g}{2}|\psi (\mathbf{x}){|}^2\right]\psi (\mathbf{x})$$

$$N(\psi )=\int d^d\mathbf{x}|\psi (\mathbf{x}){|}^2.$$

$$P(\psi)\propto\left\{\begin{array}[]{c}\exp\left(-\,\frac{E(\psi)-\mu N(\psi)}{k_{B}T}\right)\quad\text{in the GCE}\\ \exp\left(-\,\frac{E(\psi)}{k_{B}T}\right)\quad\text{in the CE.}\end{array}\right.$$

$$\mathcal{P}f(\mathbf{x})=\sum _{\mathbf{y}}{\mathcal{P}}_{\mathbf{x}\mathbf{y}}{f}_{\mathbf{y}}.$$

$${\mathcal{P}}_{\mathbf{x}\mathbf{y}}=dV\sum _{j\in \mathcal{C}}{\varphi }_j(\mathbf{x}){\varphi }_j^{\ast }(\mathbf{y}).$$

$$\mathcal{P}\mathcal{P}=\mathcal{P},{\mathcal{P}}^{†}={\left({\mathcal{P}}^T\right)}^{\ast }=\mathcal{P},{\mathcal{P}}^T={\mathcal{P}}^{\ast }.$$

$$\frac{\partial \psi (\mathbf{x})}{\partial t}$$

$$⟨\eta (\mathbf{x},t{)}^{\ast }\eta ({\mathbf{x}}^{\prime },t^{\prime })⟩={\delta }^d(\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime }).$$

$$(\mathcal{P}-1){\psi }_0(\mathbf{x})=0.$$

$$E(\alpha )=\omega |\alpha {|}^2+gc|\alpha {|}^4,$$

$$P_{GCE}(N)\propto e^{\frac{(\mu -\omega )N}{T}}.$$

$$P_{GCE}(N)\propto \exp \left[-\frac{{\left(N-N_{\mathrm{mid}}\right)}^2}{2{\sigma }_{1mode}^2}\right]$$

$$N_{\mathrm{mid}}=\frac{\mu -\omega }{2gc}$$

$${\sigma }_{1\mathrm{mode}}=\sqrt{\frac{T}{2gc}}.$$

$$g\gtrsim \frac{(\mu -\omega {)}^2}{2cT}.$$

$$P_{\sigma }(\psi )\propto \exp \left\{-\frac{\left[E(\psi )-\mu N(\psi )\right]}{T}-\frac{{\left(N(\psi )-\overline{N}\right)}^2}{2{\sigma }^2}\right\}.$$

$$\frac{\partial \psi (\mathbf{x})}{\partial t}$$

$$\frac{dv}{dt}=\sum _v{A}_v(\vec{v})+\sum _{uv}{B}_{uv}(\vec{v}){\xi }_v(t),$$

$$⟨{\xi }_u(t){\xi }_v(t^{\prime })⟩={\delta }_{uv}\delta (t-t^{\prime }),$$

$$\frac{\partial P(\vec{v})}{\partial t}=\left[-\sum _v\frac{\partial }{\partial v}{A}_v(\vec{v})+\frac{1}{2}\sum _{uv}\frac{{\partial }^2}{\partial u\partial v}{D}_{uv}(\vec{v})\right]P(\vec{v}).$$

$$D_{uv}=\sum _{v^{\prime }}{B}_{u{v}^{\prime }}{B}_{v{v}^{\prime }}.$$

$$\frac{\partial P(\vec{v})}{\partial t}=0$$

$$P(\vec{v})=P_{\sigma }(\psi )$$

$$A_{{\psi }_{\mathbf{x}}^R}$$

$$D_{{\psi }_{\mathbf{x}}^R,{\psi }_{\mathbf{y}}^R}=$$

Full text

I Introduction

While the canonical and grand canonical ensembles are two dominant ways to describe thermal ultracold Bose gases, ensembles with intermediate fluctuations of particle number are more typical in practice. The difference can matter a lot for ultracold experiments because they take place in a mesoscopic regime where such fluctuations are well resolved.

Intermediate ensembles are also theoretically interesting in their own right. In typical thermodynamic systems without excessively long-range interactions or correlations, the different statistical ensembles are known to give the same result for intensive thermodynamic quantities in the limit of a large system — ensemble equivalence Kocharovsky et al. (2016); Touchette (2015, 2011); Squartini et al. (2015). However, in a flurry of activity some years ago Navez et al. (1997); Weiss and Wilkens (1997); Grossmann and Holthaus (1996); Wilkens and Weiss (1997); Politzer (1996); Gajda and Rzążewski (1997); Grossmann and Holthaus (1997); Giorgini et al. (1998); Kocharovsky et al. (2000a, b); Xiong et al. (2001, 2002); Yukalov (2005a, b) it was found that the ideal Bose gas at ultracold temperatures, does not behave this way. Not only do its fluctuations become extremely large in the vicinity of the critical temperature, but the result depends on the ensemble that is used (canonical/grand canonical/microcanonical), even in the thermodynamic limit $N\to\infty$. Ref. Kocharovsky et al. (2006) gives an extensive review.

Now, with interactions, often even weak ones, it is thought that equivalence between ensembles is restored because the interactions energetically suppress any excessive number fluctuations Wilkens and Weiss (1997). This matter has been, and continues to be, widely debated Wilkens and Weiss (1997); Giorgini et al. (1998); Idziaszek et al. (1999); Kocharovsky et al. (2000a); Yukalov (2005a, b); Kocharovsky et al. (2006); Svidzinsky and Scully (2006); Heller and Strunz (2013); Cockburn et al. (2011a); Touchette (2011); Squartini et al. (2015); Tarasov et al. (2015). Refs. Kocharovsky et al. (2016) and Touchette (2015) explain the current understanding. Ensemble equivalence or the thermodynamic limit is often invoked to justify the use of the most convenient ensemble in calculations. Hence, the details of how ensemble equivalence is imposed by interactions and what happens in mesoscopic systems are of much interest for practical applications as well as for a theoretical understanding.

In experiment, the situation is that repeated runs usually produce a set of states that correspond to something intermediate between a CE and a GCE. The evolution of a single realization conserves particle number when particle loss is neglected, motivating many CE theoretical treatments. However, quite strong fluctuations in total atom number between runs are the norm. For example, Jaskula (2010) reports standard deviations $\delta N/N$ of about 20% or even 35% with a less optimized system, and a recent study Kristensen et al. (2017) about 10%. An ensemble with controlled atom number fluctuations is what best describes an actual set of experimental runs. Preferably, one would like an external parameter to match the degree of fluctuation to empirical observations. A further recent development in this regard are photonic BECs, because the size of the particle reservoir with which they are in contact can be varied to experimentally study the crossover between the GCE and CE in a controlled way Weiss and Tempere (2016). Experiments with photonic BECs have also been able to measure the distribution $P(N_{0})$ directly Schmitt et al. (2014).

In this paper we demonstrate what happens between the CE and GCE in the uniform one-dimensional gas and develop a convenient method to generate intermediate ensembles. The most adaptable technique to describe degenerate thermal interacting gases are ensembles of classical wave fields (“c-fields”). They are often the only way to gain quantitative access to many quantities in the regime with non-perturbative fluctuations Brewczyk et al. (2007); Blakie et al. (2008); Proukakis et al. (2013). The standard methods developed to date generate only a grand canonical ensemble (GCE) Stoof (1999); Duine and Stoof (2001); Gardiner et al. (2002); Gardiner and Davis (2003); Proukakis and Jackson (2008); Cockburn and Proukakis (2013); Bradley and Blakie (2014); Pietraszewicz and Deuar (2015), canonical ensemble (CE) Witkowska et al. (2010); Sinatra and Castin (2008); Heller and Strunz (2009); Rooney et al. (2012) or microcanonical ensemble (MCE) Góral et al. (2001); Davis et al. (2001); Berloff and Svistunov (2002); Davis and Morgan (2003); Davis et al. (2013); Brewczyk et al. (2007); Blakie et al. (2008); Karpiuk et al. (2010); Brewczyk et al. (2013) of classical fields! We develop an approach based on the stochastic projected Gross-Pitaevskii equation (SPGPE) Stoof (1999); Duine and Stoof (2001); Gardiner et al. (2002); Gardiner and Davis (2003); Proukakis and Jackson (2008); Cockburn and Proukakis (2013); Bradley and Blakie (2014) that readily generates ensembles across the entire continuum from CE to GCE. These transitional ensembles are parametrized by $\sigma$, which determines the standard deviation of the total atom number $N$. Additionally, our method gives more convenient access to the canonical ensemble than methods that hardwire exact number conservation into the system such as Góral et al. (2001); Davis et al. (2001); Berloff and Svistunov (2002); Davis and Morgan (2003); Davis et al. (2013); Brewczyk et al. (2007); Blakie et al. (2008); Karpiuk et al. (2010); Brewczyk et al. (2013); Witkowska et al. (2010).

We will also pay attention to an interesting phenomenon in the CE that has not been extensively investigated. Namely, the appearance of a “cliff” in the distribution of the number of excited particles. This occurs at relatively high temperatures $T\sim T_{c}$, when the constraint on $N$ is lower than the number of excited particles suggested by the Bose-Einstein distribution for each mode. Evidence of this feature has been seen in both the ideal Grossmann and Holthaus (1996); Wilkens and Weiss (1997); Weiss and Wilkens (1997); Svidzinsky and Scully (2006); Witkowska et al. (2009, 2010); Tarasov et al. (2015) and 1d interacting gas in the CE Carusotto and Castin (2003); Bezett and Blakie (2009); Bisset and Blakie (2009); Cockburn et al. (2011a); Bienias et al. (2011a). No investigation has been made in the GCE with interactions. It is interesting to find out how robust this phenomenon is to a breaking of the extreme constraint on $N$ that occurs in the CE.

Prior to that, we derive and describe the SPGPE method for transitional and canonical ensembles in Sec. II. We benchmark it on the ideal gas in Sec. III, calculate the distributions and fluctuations in the interacting gas in Sec. IV and look into ensemble equivalence and the “cliff” in Sec. V.

II Stochastic method for canonical and transitional ensembles

II.1 The system and its c-field description

We will consider a single-species gas of contact-interacting bosons. With the Bose field $\widehat{\Psi}({\mathbf{x}})$, the Hamiltonian is written

[TABLE]

in $d$ dimensions. The contact interaction strength is $g$, and $H^{\rm sp}$ is the single particle energy:

[TABLE]

On a discretized spatial lattice with small volume per lattice point $dV$, as often used for calculations, (1) becomes

[TABLE]

The Hermitian nature of $H^{\rm sp}$ implies

[TABLE]

We will use this discretized representation interchangeably with the continuous one, according to convenience.

In a minimalist view: the c-field (classical wave field) description boils down largely to an assumption that the relevant behavior of the system is captured by the highly occupied modes, while those with ${\mathcal{O}}(1)$ occupation or less can be neglected. Two complementary reviews of the c-field approach are Brewczyk et al. (2007) and Blakie et al. (2008). We can write the quantum Bose field in terms of orthogonal modes labeled $j$ with mode functions $\phi_{j}({\mathbf{x}})$ normalized to unity and annihilation operators $\widehat{a}_{j}$:

[TABLE]

Then, the c-field approximation corresponds to

[TABLE]

where ${\mathcal{C}}$ is the subspace of high-occupied modes and $\alpha_{j}$ are complex values that approximate the $\widehat{a}_{j}$. The $\{\ \dots\ \}$ indicates that the quantum operator $\widehat{\Psi}$ is in general going to be described by an ensemble. The numbers $\alpha_{j}$ will differ among different elements of the ensemble. The subspace ${\mathcal{C}}$ is generally chosen a priori and specified by an energy cutoff $E_{\rm cut}$, such that all single particle modes with energies below this cutoff are included in ${\mathcal{C}}$ and all above excluded. This is the most consistent choice for systems that lie close to thermal equilibrium, since occupations will decrease monotonically with energy. A recent detailed study of a broadly applicable cutoff choice for ultracold interacting gases is Pietraszewicz and Deuar (2017).

The c-field Hamiltonian for the low-energy part of the system takes the form:

[TABLE]

and the number of particles is

[TABLE]

The distribution of $\psi({\mathbf{x}})$ is then written as

[TABLE]

We will use $\hbar=m=k_{B}=1$ units in what follows.

II.2 Generation of ensembles

The two most widespread approaches to produce a c-field ensemble for interacting particles involve generating samples from the ergodic time-evolution of an initial state using the Gross-Pitaevskii equation (GPE) and its variants. They are:

(1) Evolution of an initial state using the deterministic but ergodic GPE Góral et al. (2001); Berloff and Svistunov (2002); Brewczyk et al. (2007); Karpiuk et al. (2010); Brewczyk et al. (2013) or its projected version (PGPE) Davis et al. (2001); Davis and Morgan (2003); Blakie et al. (2008); Davis et al. (2013). This corresponds to isolated Hamiltonian evolution of the classical wave field and produces a microcanonical ensemble (MCE) with number and energy set by the initial state.

(2) Evolving a stochastic Gross-Pitaevskii equation (SGPE) Stoof (1999); Duine and Stoof (2001); Gardiner et al. (2002); Proukakis and Jackson (2008); Cockburn and Proukakis (2013) or its more general projected version (SPGPE) Gardiner and Davis (2003); Bradley and Blakie (2014), which corresponds to a model where the above-cutoff modes that are excluded from ${\mathcal{C}}$ are approximated as a particle and energy bath. This produces a grand canonical ensemble (GCE) with chemical potential $\mu$ and temperature $T$ set externally.

Regarding the canonical ensemble, several methods to generate an interacting classical wave field ensemble have been proposed:

A Metropolis algorithm for generating samples with a CE probability. It involves a discrete random walk taken with steps that conserve particle number Witkowska et al. (2010). It has been used in several studies since Witkowska et al. (2011); Karpiuk et al. (2012); Bienias et al. (2011a); Pawłowski et al. (2013); Lang and Witkowska (2014); Karpiuk et al. (2015); Gawryluk et al. (2015); Pawłowski and Rzążewski (2015); Pietraszewicz and Deuar (2015); Gawryluk et al. (2017), and is also easily adapted to the GCE Pietraszewicz and Deuar (2015, 2017). 2. 2.

A particle number filter applied to grand canonical ensembles to obtain a CE, though this is a wasteful procedure. 3. 3.

The noise modifications and projections laid out in Heller and Strunz (2013) constitute another method. 4. 4.

Rooney et al. found that an SPGPE with no particle exchange terms produces a canonical ensemble if the exotic scattering terms are included Rooney et al. (2012). 5. 5.

Another approximate approach that can be very accurate under the right conditions used a Bogoliubov description for excited atoms supplanted with as many condensate atoms as required to match the total assumed atom number Sinatra et al. (2007).

The above canonical ensemble approaches are not always easily done, especially when one wants to have e.g. a set magnetization in spinor or multi-component condensates. For example, we have found that ensuring the right conservation law while preserving detailed balance becomes very tricky with Metropolis for multicomponent gases. It is known also that the implementation of the scattering-only SGPE is a nontrivial endeavor even for single component gases Blakie et al. (2008); Rooney et al. (2012).

In this paper we derive an alternative approach that extends the SPGPE to incorporate a controllable number filter. It restricts the evolution to the vicinity of the CE and avoids the waste of discarding realizations. Moreover, it gives access to natural intermediate ensembles between the marginal CE and GCE.

II.3 The SPGPE

The SPGPE is a flexible way to generate the grand canonical ensembles of classical wave fields given by (9) Blakie et al. (2008); Staliunas (2006). It has been described in detail in Duine and Stoof (2001); Gardiner and Davis (2003); Proukakis and Jackson (2008) and benchmarked in Cockburn et al. (2011a, b).

In general, one works in a projected subspace ${\mathcal{C}}$. It is imposed by acting with a projector ${\mathcal{P}}$ onto spatially dependent fields $f({\mathbf{x}})$ such that ${\mathcal{P}}f({\mathbf{x}})$ lies wholly within ${\mathcal{C}}$. This allows one to work on a simple spatial grid ${\mathbf{x}}$ while restricting the basis in any desired way. On the spatial lattice, this becomes

[TABLE]

An explicit form of the matrix elements is

[TABLE]

The projector ${\mathcal{P}}$ fulfills the usual properties:

[TABLE]

The simplified case of a plane wave basis with cutoff set implicitly by the lattice corresponds to setting ${\mathcal{P}}\to 1$.

The time evolution of $\psi({\mathbf{x}})$ is governed by the SPGPE

[TABLE]

It corresponds to coupling the c-field to a thermal bath at temperature $T$ and chemical potential $\mu$. The dimensionless positive coupling strength is $\gamma({\mathbf{x}})$, which is commonly taken to be constant in space. Such an assumption is certainly a convenience if one is primarily interested in the long-time ensemble, rather than the transient dynamics. Physically justified values are usually small ($\gamma\ll 1$). $\eta$(${\mathbf{x}},t$) is a complex white noise field independent at each spatial position and time, with zero mean, and variance:

[TABLE]

The equation (13) changes only that part of the field $\psi$ that has support in the c-field subspace ${\mathcal{C}}$. To be self-consistent and physically sensible we need the initial state $\psi_{0}({\mathbf{x}})$ to be fully in this subspace, i.e.

[TABLE]

To obtain a GCE one evolves the equation until transients related to the initial state have died off, and only thermally activated fluctuations remain. Let us call this time $t_{*}$. The choice of initial state is, in principle, irrelevant, although it may affect the length of time needed to reach the thermally activated regime. Starting from vacuum $\psi({\mathbf{x}},0)=0$ is a common choice. Independent samples of the distribution can then be obtained from values of the field $\psi({\mathbf{x}},t)$ sufficiently well spaced in time after $t_{*}$. This is reminiscent in many ways of the procedure with the Metropolis method, except that all updates are accepted, and given explicitly by the noise term. Alternatively, one can simply evolve from the same initial state to $t=t_{*}$ but using a different noise realization each time, and the fields $\psi({\mathbf{x}},t_{*})$ will be the independent samples of the GCE. The latter approach removes the need to investigate time correlations. It also simplifies the determination of $t_{*}$ because ensemble-averaged quantities can be tracked for a number of times leading up to $t\gtrsim t_{*}$, to verify when the stationary ergodic ensemble has been reached.

The system evolves to the GCE distribution (9) regardless of the details of $\gamma({\mathbf{x}})>0$, which only affects the time $t_{*}$.

II.4 One mode and ensemble equivalence

It is instructive in the beginning to look at the behavior of ensembles in a single mode, $j=0$ say. We will see that this example encapsulates both the basic physics of how ensemble equivalence is restored by interactions, and suggests a naturally occurring form for the manifold of intermediate ensembles.

The Hamiltonian (7) in the c-field description can be written as an energy

[TABLE]

that depends on the amplitude $\alpha$. The coefficients depend on the shape of the mode function according to $\omega=\int d^{d}{\mathbf{x}}\,\phi_{0}({\mathbf{x}})^{*}H^{\rm sp}\phi_{0}({\mathbf{x}})$ and $c=\tfrac{1}{2}\int d^{d}{\mathbf{x}}\,|\phi_{0}({\mathbf{x}})|^{4}$, while the number of atoms is $N(\alpha)=|\alpha|^{2}$. According to (9) the distribution of the states in the GCE is $P_{GCE}(\alpha)\propto\ \exp\left\{-\,\frac{(\omega-\mu)|\alpha|^{2}+gc|\alpha|^{4}}{T}\right\}$ with all values of $\alpha$ represented.

For an ideal gas, the above exponent produces a very broad distribution of particle number

[TABLE]

This is the most trivial case of the GCE fluctuation catastrophe and inequivalence of ensembles, since the fluctuations of $N$ scale as $N$. In fact $\delta N=\sqrt{\langle N^{2}\rangle-\langle N\rangle^{2}}=\frac{\omega-\mu}{T}=N$, so $P_{GCE}(N)$ never approaches the CE behavior of $\delta N=0$, even as $N\to\infty$.

Interactions, however, make the distribution Gaussian:

[TABLE]

Now we can see that, with the help of the interaction $g$, one can drive the standard deviation of the Gaussian, $\sigma_{1\rm mode}$, to smaller values. Eventually, the fluctuations in particle number, $\delta N$, become $\sigma_{\rm 1mode}\propto 1/\sqrt{g}$. Thus, for large $N$ and $\sigma_{\rm 1mode}\ll N_{\rm mid}$, ensemble equivalence is restored in the thermodynamic limit $N\to\infty$ because $\delta N/N\to 0$ like in the CE. In terms of $g$, this happens for

[TABLE]

This example shows the essence of how ensemble equivalence is restored by interactions.

II.5 Restricting the atom number in the SPGPE via an additional term

It would be convenient to have an equation that explicitly conserves $N(\psi)$ to a set value ${\overline{N}}$ but keeps a similar form as the SPGPE (13). And indeed – the one mode toy problem of Sec. II.4 suggests a way: Terms of a similar form to the interaction term should be capable of imposing a Gaussian distribution of $N$ with a width of our choice, while leaving the rest of the system evolution largely unchanged. In the limit of a narrow Gaussian distribution around the desired value, we would have effectively a CE distribution.

Consider an additional Gaussian factor to (9) thus:

[TABLE]

When $\sigma$ becomes smaller than other widths, only fields $\psi({\mathbf{x}})$ with a number of particles ${\overline{N}}\pm\sigma$ occur with non-negligible probability. In the limit of small $\sigma$ this becomes effectively a CE with ${\overline{N}}$ particles. What terms should be added to the SPGPE to attain this modification?

First note that the exponent of the probability distributions of c-field states contain all the terms of the Hamiltonian, converted to a classical field and scaled. Secondly, the deterministic parts of the SPGPE correspond to the classical field simplification of the Heisenberg equations of motion for the quantum field $\widehat{\Psi}$. i.e. of $d\widehat{\Psi}({\mathbf{x}})/dt=-i\left[\widehat{\Psi}({\mathbf{x}}),\widehat{H}\right]$. Hence, each term $\widehat{H}_{j}$ in the Hamiltonian leads to the c-field version of $-{\mathcal{P}}\{(i+\gamma)\left[\widehat{\Psi}({\mathbf{x}}),\widehat{H}\right]\}$ in the SPGPE. Taken together, these two points suggest that the new term in (20) proportional to the c-field version of $(\widehat{N}-{\overline{N}})^{2}$ will generate a term in the stochastic equation proportional to the c-field version of $[\widehat{\Psi}({\mathbf{x}}),(\widehat{N}-{\overline{N}})^{2}]$. That is, one may expect terms proportional to $\psi({\mathbf{x}})(N(\psi)-{\overline{N}})$. Let us postulate, then, a modified SPGPE:

[TABLE]

with a constant $K$ (possibly space-dependent) to be determined. We will see if and for what value of $K$ the stationary distribution is equal to the desired (20). Note that we have placed the term inside the projection ${\mathcal{P}}$ because the equation should always preserve the property that $\psi({\mathbf{x}})$ has support only in the ${\mathcal{C}}$ subspace to have a consistent c-fields description.

The correspondence between stochastic equations, Fokker-Planck equations for the distribution, and stationary states is well known. A detailed explanation can be found e.g. in Gardiner (2009). When one has a set of real variables $\vec{v}=\{v\}$ governed by Langevin stochastic equations of the form

[TABLE]

with real noises of zero mean and correlations

[TABLE]

then it is a realization of the following Fokker-Planck equation (FPE) for the probability distribution $P(\vec{v})$ of the variables:

[TABLE]

Summing is over all variables in $\vec{v}$, and derivatives act on all factors to the right. The diffusion matrix $D$ is given by elements

[TABLE]

The desired stationary distribution is (20), so that we want to impose

[TABLE]

when the substitution

[TABLE]

is made in the FPE (24).

Consider the system on a numerical lattice as in (3), so that the set of variables in $\vec{v}$ consists of the real and imaginary components of $\psi$ at each point ${\mathbf{x}}$, i.e. $\psi^{R}_{{\mathbf{x}}}={\rm Re}[\psi_{{\mathbf{x}}}]$ and $\psi^{I}_{{\mathbf{x}}}={\rm Im}[\psi_{{\mathbf{x}}}]$, respectively. The equation (21) can be rewritten in the general $A,B,D$ notation of (22)-(25) using the coefficients

[TABLE]

and

[TABLE]

In the interest of clarity we will assume for now a constant $\gamma$ and $K$:

[TABLE]

and report the more general result from Appendix A at the end of this section.

Substituting (27)-(30) into (24) leads (after much algebra) to the FPE:

[TABLE]

The first line of (31) can easily be made zero with an appropriate choice of $K$, but even then, the second line still remains potentially troublesome. However, note that the properties of the equation (21) and initial state (15) ensure that the field $\psi$ stays in the c-field subspace in the overall model. Then,

[TABLE]

Using this, (31) becomes

[TABLE]

The following choice of prefactor on the first line:

[TABLE]

makes the distribution (20) stationary. This is exactly what we required. The final equation to simulate is then simply:

[TABLE]

Appendix A explains what happens when we relax condition (30) allowing the coupling strength $\gamma({\mathbf{x}})$ to be spatially dependent.

The condition needed to obtain a well behaved equation is that $\gamma({\mathbf{x}})$ varies slowly in the region of space around ${\mathbf{x}}$ compared to ${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}$. i.e.

[TABLE]

It is met in most realistic cases. Then, the choice (34) turns out to generalize to

[TABLE]

and the modified SPGPE equation that keeps (20) stationary is

[TABLE]

This is just the usual SPGPE with one extra term.

The most likely place for a breakdown of the condition (36) in typical problems is in the low density tails of a trapped system that is described using harmonic oscillator modes. Away from the main cloud, near the energy cutoff, only very low wavelength parts of modes $\phi_{j}({\mathbf{x}})$ are present, and then ${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}$ may vary on comparable scales to $\gamma({\mathbf{x}})$.

The special but common case of an unprojected “plain” SGPE, where the only projection is an implicit one imposed by the numerical lattice (${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}=\delta_{{\mathbf{x}}{\mathbf{y}}}$) also uses (38) with ${\mathcal{P}}\to 1$, and without the need for the slowly varying condition on $\gamma$ that is (36).

To conclude this derivation, one can safely say that the canonical ensemble has been achieved for small $\sigma$ when all relevant observable quantities cease to depend on $\sigma$ in any significant way.

The equations (35) and (38) are a convenient way by which one can generate the CE. Both equations are stable, straightforward to integrate, and require fewer numerical tweaks than a Metropolis algorithm. We only need to set $\gamma$, which can be chosen over a wide range without ill effect, when the purpose is to generate a stationary ensemble. Furthermore there is no wastage due to particle number filtering.

The equations (35) and (38), can also be used to produce the dynamics of a canonical ensemble, but then one should determine a correct value and spatial dependence of the reservoir coupling $\gamma({\mathbf{x}})$. The question of how realistic the physical model is remains somewhat open, since the system corresponds to having a low-energy part of the field that exchanges only energy but not particles with the high energy components that are treated as a bath. Nevertheless, such a model has been discussed in some detail in the context of a scattering-only SPGPE Rooney et al. (2012), and may be useful in various situations.

Overall, the computational cost scales the same way as in other treatments based around the SGPE. That is, a very lenient $M\log M$ scaling with the number of points on the computational lattice $M$, regardless of the dimensionality. This makes it convenient for 2d and 3d systems. The usual limiting factor is the efficiency of a Fourier transform used to evaluate kinetic energy. An issue to keep in mind is that very small values of $\sigma$ will shorten the required timestep by virtue of introducing a large gradient. It may also tend to increase the time $t_{*}$ required to obtain stationarity. Some precursors of this were seen at the lowest values of $\sigma$ in our 1d calculations.

II.6 Transitional distributions between CE and GCE

Equation (38) generates the family of ensembles (20) as its long time stationary distribution. These span the whole continuum between CE and GCE for interacting systems, with the location on the continuum given by $\sigma$.

A convenient way to specify distributions intermediate between CE and GCE is through the standard deviation of the atom number fluctuations $\delta N$. This captures the foremost difference between the CE and GCE, and can be readily matched to experimental data such as in Jaskula (2010).

There are two contributions to $\delta N$: First, the “natural” one ($\delta^{GCE}N$) that arises as a result of the interplay of the interaction strength $g$ and the particle bath described by the chemical potential $\mu$. For the single mode this is (18c). Then there is also the externally steered fluctuation $\sigma$. It will not increase fluctuations beyond the natural level, but can decrease them. Hence, we expect that

[TABLE]

The largest values of $\sigma$ do not affect the GCE much. Then, when $\sigma$ becomes small enough to limit the natural fluctuation width, it begins to meaningfully steer the distribution. Finally, when $\sigma$ becomes small enough that observable quantities cease to change, we have reached the CE. This changeover will be seen later in Fig. 5.

In the ideal gas, for small enough $\sigma$, the center of the Gaussian-like distribution for $N$ can be quite well estimated by

[TABLE]

This comes from inspection of (20) while omitting the $H^{\rm sp}$ contribution. ${\overline{N}}_{\sigma}$ converges to ${\overline{N}}$ in the CE limit. For large $\sigma$, (40) becomes inaccurate because other factors come into play, such as a nontrivial $H^{\rm sp}$ contribution and the fact that the distribution of $N$ is nonzero only for $N>0$.

In a uniform interacting gas in volume $V$, the properties of the Gaussian can also be estimated. The energy functional is $E(\psi)=\varepsilon_{\rm sp}N(\psi)+\frac{g}{2V}\,g^{(2)}(0)N(\psi)^{2}$, where $\varepsilon_{\rm sp}$ is the mean energy per particle from the single-particle Hamiltonian $H_{\rm sp}$, and $g^{(2)}(0)$ is the density-density correlation function. $g^{(2)}(0)$ lies between 1 and 2 in an equilibrium ensemble. Looking first at the natural GCE in (9), the Gaussian distribution for $N(\psi)$ is centered at

[TABLE]

with a standard deviation

[TABLE]

For the transition distributions $P_{\sigma}$ of (20), the center of the Gaussian for $N(\psi)$ shifts to

[TABLE]

and the standard deviation becomes

[TABLE]

One can see that indeed in the $\sigma\ll s_{GCE}$ limit, both quantities converge to the externally set values of ${\overline{N}}$ and $\sigma$, while in the opposite $\sigma\gg s_{GCE}$ limit, the natural GCE behavior reasserts itself.

Unplanned behavior can occur if the difference between ${\overline{N}}$ and the natural ${\overline{N}}_{GCE}$ is much greater than $\sigma$. In that case, the external constraint ${\overline{N}}$ and the internal chemical potential $\mu$ work against each other. The result is a relatively narrow distribution that is not centered near ${\overline{N}}$ but at a weighted average of ${\overline{N}}$ and ${\overline{N}}_{GCE}$ given by (43). The upshot of this for generation of canonical ensembles in general cases is that one should check the actual resulting mean particle number. If it does not closely match ${\overline{N}}$, then $\mu$ should be modified to bring ${\overline{N}}_{GCE}$ close to ${\overline{N}}$.

III Ideal gas

Let us first check the method on the ideal gas, where exact results are available. The typical observables studied in the context of comparing ensembles are the distributions $P(N_{0})$ and $P(N_{\rm ex})$ of the number of atoms in the ground or excited states, as well as their moments. An experimental method for measuring fluctuations in the condensate occupation $N_{0}$ has been proposed in Bezett and Blakie (2009).

III.1 Procedure

We treat here a 1d uniform gas, and the procedure outlined below was applied for both ideal and interacting gases. The chosen basis consists of plane waves $\phi_{j}(x)=e^{ik_{j}x}/\sqrt{L}$ defined in a box of length $L$ with periodic boundary conditions. Wave vectors are $k_{j}=2\pi j/L=j\Delta k$ with $j=\{0,\pm 1,\pm 2,\dots\}$. We take $L$ to be the computational unit of length in what follows, and only write it out explicitly in a few cases to show scaling. The c-field subspace ${\mathcal{C}}$ is implemented using a maximum kinetic energy cutoff for the plane waves $E_{c}=\hbar^{2}k_{c}^{2}/2m$.

We revisit the regimes that were investigated in the past work of Witkowska et al. (2010) (Fig. 1). Namely, we study a similar condensate fraction $n_{0}=N_{0}/N$ and distribution of excited atoms. We fixed the target total atom number in the CE at the higher value of ${\overline{N}}=500$.

It is convenient to give the temperature scaled with respect to a finite-size characteristic temperature for condensation. In the ideal gas canonical ensemble, the occupation of excited modes is given by $N_{k}\approx N_{k}^{\rm BE}=[e^{k^{2}/2T}-1]^{-1}$, provided the total number of excited atoms $N_{\rm ex}=\sum_{j\neq 0}N_{k_{j}}$ does not reach ${\overline{N}}$. Otherwise, it invokes the constraint and mode occupations reduce below $N^{\rm BE}_{k}$. To estimate the temperature $T_{*}$ below which a significant condensate will appear, one can evaluate the simple condition $N_{\rm ex}(T_{*})={\overline{N}}$ using the estimates $N_{k}^{\rm BE}$. In our particular case of ${\overline{N}}=500$, we find $T_{*}=3195/L^{2}.$ The simplest general estimate comes from considering only the two lowest lying excited states, in which case $T_{*}=\pi^{2}{\overline{N}}/L^{2}$.

The cutoff used for calculations was the recommended value for matching the condensate fraction and $P(N_{\rm ex})$ in a 1d ideal gas in a box in the CE Witkowska et al. (2009) 111This corresponds to $f_{c}=1.9023$ in the global optimized cutoff notation of Pietraszewicz and Deuar (2015, 2017), where $k_{c}=f_{c}\sqrt{2\pi T}$., i.e. $k_{c}^{2}=0.58T$. This leads to a cutoff of $4\Delta k$ for the low temperature $T=0.341\,T_{*}$ of Sec. III.2 and $k_{c}=8\Delta k$ for the high temperature $T=1.365\,T_{*}$ of Sec. III.3, like in the work of Witkowska et al. (2010).

The generation of each member of the ensemble proceeds by starting with the vacuum state $\psi_{{\mathbf{x}}}(0)=0$ on a numerical lattice with spacing $\Delta x=L/2^{7}$. Note that the maximum allowable wavevector on this lattice, $k_{\rm max}=\pi/\Delta x$ is much greater than the cutoff $k_{c}$. This allows us to accurately calculate the interacting evolution, which would otherwise suffer from some small but spurious aliasing and umklapp processes on a lattice with $k_{\rm max}=k_{c}$. The state $\psi_{{\mathbf{x}}}(t)$ is then evolved using (35) with a constant value of $\gamma$ until a stably randomly fluctuating solution is reached above $t_{*}$. We used values of $\gamma$ in the range 0.01 to 0.1. This is repeated for each sample, using a new set of noises $\eta_{{\mathbf{x}}}(t)$. The stationarity of the ensemble is checked by tracking ensemble averages of various observables, and this allows us to determine appropriate $t_{*}$. These times were $t_{*}={\mathcal{O}}(10/\gamma)$, with some variation depending on parameters and $\sigma$.

When changing $\sigma$ to move between the CE and GCE, we keep the chemical potential $\mu$ constant for each temperature and interaction strength. This assumption aids in obtaining a sequence of physically related intermediate ensembles. The value of $\mu$ is chosen so that the mean number of atoms $\langle N\rangle$ in the GCE matches the CE value of ${\overline{N}}$. This helps to avoid the possible competition between ${\overline{N}}$ and $\mu$ that was discussed in Sec. II.6.

III.2 Low temperature case

Let us consider first a low temperature case in which the majority of atoms are in the condensate. This is the regime in which $P(N_{0})$ or $P(N_{\rm ex})$ distributions have most commonly been described. For example Grossmann and Holthaus (1996); Weiss and Wilkens (1997); Wilkens and Weiss (1997); Sinatra et al. (2002); Svidzinsky and Scully (2006); Witkowska et al. (2009, 2010); Bienias et al. (2011a); Heller and Strunz (2013); Tarasov et al. (2015) in the ideal gas, and Sinatra et al. (2002); Carusotto and Castin (2003); Bezett and Blakie (2009); Idziaszek et al. (2009); Bisset and Blakie (2009); Witkowska et al. (2010); Cockburn et al. (2011a); Bienias et al. (2011a); van der Wurff et al. (2014) in the interacting. One reason for its popularity is that it is accessible by the Bogoliubov approximation.

Fig. 1 shows the distribution of the number of excited atoms $N_{\rm ex}$ at $T=0.341\,T_{*}$ in the CE and GCE. The CE has ${\overline{N}}=500$ and the ensemble is obtained using (35) with $\sigma=1$. In the c-field description, $N_{\rm ex}=\sum_{k\neq 0}|{\widetilde{\psi}}_{k}|^{2}\Delta k$, and ${\widetilde{\psi}}$ is the Fourier transformed field

[TABLE]

The GCE is obtained using a simple SPGPE (13), in the limit $\sigma\to\infty$. As explained in Sec. III.1, the chemical potential is chosen so that the mean number of atoms $\langle N\rangle$ in the GCE matches the CE value of ${\overline{N}}=500$. This is $\mu=-2.135/L^{2}$ here. The numerical histograms are compared to exact results which are obtained in Appendix B. We see that despite the not so large shift from CE to GCE in this regime, the distribution tracks it in detail. The histogram is from ${\mathcal{S}}=2.5\times 10^{4}$ samples of the ensemble.

Fig. 2 shows the behavior of the probability distribution of the total atom number, $P(N)$, as $\sigma$ is varied through the transition ensembles. This is for the same low temperature case. As expected, we move from an extremely broad distribution in the GCE, through a Gaussian (initially broad, later narrow) which converges to an extremely narrow distribution around $N={\overline{N}}$.

III.3 High temperature case

Distributions for the high temperature case have been reported for the ideal gas Grossmann and Holthaus (1996); Wilkens and Weiss (1997); Weiss and Wilkens (1997); Svidzinsky and Scully (2006); Witkowska et al. (2009, 2010); Tarasov et al. (2015) and interacting gas Carusotto and Castin (2003); Bezett and Blakie (2009); Bisset and Blakie (2009); Cockburn et al. (2011a); Bienias et al. (2011a) in this regime. They behave very differently, though this has not been analyzed as much in the literature.

Fig. 3 shows the CE, GCE, and two intermediate ensembles for $T=1.365\,T_{*}$. We set ${\overline{N}}=500$ and use $\mu=-32.789/L^{2}$ to have matching $\langle N\rangle=500$ in the GCE. The match to exact CE and GCE results is ideal. Particularly notable is the reconstruction of the CE “cliff” in $P(N_{\rm ex})$ despite the total atom number not being hardwired into the simulation, and all values being at least in principle allowed.

This is an unusual regime in a number of aspects. Apart from the presence of the sharp cliff, another interesting feature appears. Namely, the most commonly occurring values of the number of excited atoms are larger in the CE (and the $\sigma=30$ case) than in the GCE. They are around 400 versus 300, respectively. This is rather counterintuitive compared to the usual impression that the GCE in the ideal gas allows much larger numbers of excited particles. What we observe here is a consequence of the strong restriction on allowable states that the CE (or low $\sigma$) condition imposes.

Further ideal gas results will appear as limiting cases in the later discussion and Figs. 4 to 9, such as $P(N_{0})$, mean values of $\langle N\rangle$, and fluctuations of $N$ and $N_{0}$.

IV Interacting gas and transitional ensembles

Having verified that the method reproduces the expected ideal gas distributions exactly, we now turn to a more detailed analysis of the effect of nonzero interactions on the transitional ensembles.

Fig. 4 shows how the relative fluctuation of the number of condensate atoms

[TABLE]

changes with $\sigma$ and $g$ in the low temperature case $T=0.341\,T_{*}$. The behavior of this quantity when temperature, $N$ or interaction are changed has been studied extensively in the standard ensembles (CE,GCE,MCE) Grossmann and Holthaus (1996); Weiss and Wilkens (1997); Navez et al. (1997); Grossmann and Holthaus (1997); Wilkens and Weiss (1997); Sinatra et al. (2002); Svidzinsky and Scully (2006); Idziaszek et al. (2009); Witkowska et al. (2010); Cockburn et al. (2011a); Bienias et al. (2011a, b); Heller and Strunz (2013); Weiss and Tempere (2016), but not the transition between ensembles or experimentally relevant intermediate cases.

On the figure, the ideal gas case appears as hollow symbols, and unsurprisingly has the highest relative fluctuations. We see two plateau regions. The first, for $\sigma\lesssim 20$, in which there is no discernible difference from the CE. The size of this range is related to the width of typical features in the distribution of $P(N_{0})$. When the allowable fluctuation in $N$ (which is $\sim\sigma$) becomes several times smaller, it will cease to visibly affect $P(N_{0})$. For example, in Fig. 1 features in the distribution of $P(N_{\rm ex})$ have a width of ${\mathcal{O}}(50)$ atoms, and the same applies for $P(N_{0})$ in the CE.

The second plateau area for $\sigma\gtrsim 1000$ displays the same magnitude of fluctuations of $N_{0}$ as in the GCE. Note that this is a point where $\sigma\sim{\overline{N}}$, and indeed we would not expect the Gaussian narrowing caused by $\sigma$ to affect much if it is significantly broader than the natural size of $N_{0}$ fluctuations in the GCE.

As interaction is raised, initially only the fluctuations in the GCE are affected because they are large. This starts for quite small interaction strengths. As interaction grows, the GCE-like region expands somewhat to lower values of $\sigma$. Eventually, though, for strong-enough interaction, the fluctuations of $N_{0}$ begin to reduce also in the CE, somewhat unexpectedly.

Figure 5 shows the relative fluctuation of the total number of atoms, $N$. This has not been studied so much in the standard ensembles, primarily because not much happens in those cases (e.g. in the CE or MCE). For intermediate ensembles, we also see the two plateau regions in the limits of $\sigma$ and strong reductions in fluctuation with increasing $g$.

A comparison of $\delta N/\langle N\rangle$ with the naively expected effect of only the Gaussian narrowing (which is $\sigma/{\overline{N}}$) in Fig. 5 shows that the total particle fluctuations track this estimate faithfully from small $\sigma$ up to $\sigma\approx 150$. This agrees with (39). Note that $\delta^{GCE}N\approx 400$, and the center of the transition between $\sigma$-limited behavior and natural GCE behavior in Fig. 5 is also around this value of $\sigma$.

To see in more detail what goes on in the transitional ensembles, Fig. 6 shows the mean atom number as a function of $\sigma$. There is some (mostly minor) variation despite $\mu$ being chosen to match CE and GCE mean atom numbers in the two limiting cases. At low values of $\sigma$, the ideal gas behaves as predicted by (40) (red line). Overall, there is a dip at intermediate $\sigma$. This comes about because the GCE distribution of $N$ has a positive skewness (long tail at high $N$). The tail is more strongly suppressed by the Gaussian multiplier in (20) than the low $N$ part of the distribution, because the latter is closer to the mean.

Finally, Fig. 7 shows the relative fluctuations of the ground state (“condensate”) occupation for the high temperature case $T=1.365\,T_{*}$. The relative fluctuations are large even in the CE. Note that this occupation is still appreciable (in the range 0-300) despite $T>T_{*}$, since we are considering a mesoscopic system, not one in the limit of $N\to\infty$. The usual plateau behaviors seen before in Fig. 4 are also present. Differently from low $T$, it takes a rather strong interaction $g\geq 0.01$ (units of $1/L$) to invoke a response in the relative fluctuations. Moreover it is difficult to bring the fluctuations down to the CE level by interactions alone, despite the CE value of $\delta N_{0}/N=0.66$ being very high. e.g. $g=2$ is still not fully sufficient.

V Equivalence of ensembles

Let us now look in some more detail at the matter of the equivalence of ensembles as interaction is increased. In most studies, the relative fluctuation of $N_{0}$ is the quantity that has been investigated in this context. For equivalence, we expect the CE and GCE values to be equal (or ideally, to see a horizontal line across all $\sigma$ values). This is indeed what is seen in Fig. 4 for the highest values of $g$ (e.g. it is close for $g=0.3$ and truly equal at $g=1.5$). A very good match for our parameters occurs only once the CE value (and the rest) have fallen below the ideal gas value due to interactions. At this stage we have insufficient information to state whether this is a general feature.

The detailed behavior of the distributions is shown for this low $T$ case in Fig. 8, in which each panel compares the CE and GCE distributions. The lower panels show the ideal gas, and apart from the huge discrepancy in $P(N_{0})$, we see that $P(N_{\rm ex})$ also differs. At the strong interaction of $g=2$, however, both distributions have become very close. A small difference in $P(N_{0})$ remains, though it is of a size that would often be inconsequential operationally.

The above plots quantitatively validate many existing intuitions about ensemble equivalence.

The system is somewhat more resistant to ensemble equivalence at the higher temperatures above $T_{*}$ in Fig. 7. The behavior of the distributions (in which we expect the appearance of the cliff) is shown for this case in Fig. 9. In the ideal gas the CE and GCE distributions are dissimilar for both $N_{\rm ex}$ and $N_{0}$, though the difference for $N_{0}$ is nowhere near as great as at lower temperatures. In the end though, at the high interaction value of $g=8$ shown (larger than in Fig. 7), the shape of the CE and GCE distributions in the upper panels have become very close. The conclusion is that ensemble equivalence has been restored by interactions also here. This is despite the complicated form of the CE/GCE distribution itself.

Of particular note is the fact that the “cliff” near $N_{\rm ex}=500$ is also present in the GCE! In the ideal gas, this feature was due exclusively to the property $P(N_{0})=P({\overline{N}}-N_{\rm,ex})$ caused by the hard-wired external constraint ${\overline{N}}=500$ in the CE, and was rather fragile with respect to a change of the ensemble. For example at $\sigma=30$, the cliff has already practically disappeared in Fig. 1, while other quantities such as $\delta N_{0}/\langle N_{0}\rangle$ are more robust, and have hardly budged from their CE value. However, in Fig. 9 interactions impose the cliff again in the GCE where there is absolutely no explicit constraint on particle number. This quantitatively validates ensemble equivalence at high temperatures, including the preservation of features that appear fragile in the ideal gas.

VI Conclusions

We have derived extended SPGPE-like equations that generate a canonical ensemble (rather than grand canonical) as the stationary state. This occurs in the limit of small $\sigma$, which is a control parameter for the allowed fluctuations of the total atom number

[TABLE]

The equations are easily adaptable to arbitrary external potentials and nonlocal interactions. A major added benefit is the possibility to readily generate a whole range of intermediate ensembles between canonical and grand canonical. The relationship (47) makes it quite simple to use it to match the true experimental variability of atom number in a run with many experimental realizations.

We have tested the ensembles produced (Sec. III), and also shown their utility for studying the transition to the canonical ensemble and the onset of ensemble equivalence as interaction grows. We have also drawn attention to the unusual behavior of canonical ensembles with low numbers of atoms: the appearance of a cliff in the distribution of atoms in excited modes, and its retention also in the grand canonical ensemble when interactions are sufficiently high.

It is hoped that the method will be useful for the study of canonical ensembles, other experimentally obtained ensembles that do not fit neatly into the CE/GCE categorization, as well as for the study of ensemble equivalence and other related phenomena. Here, we have quantitatively validated the ensemble equivalence scenario and shown the details of how it gradually appears with growing $g$ (Figs. 4,5,7-9). Interestingly and importantly, we find that even canonical ensemble features such as the “cliff” that are fragile to a weakening of the canonical ensemble constraint in the ideal gas can nevertheless be robustly reproduced when interactions become strong.

Importantly, the equations should be readily applicable to multiple components. One $-\gamma(N_{j}-{\overline{N}}_{j})\psi_{j}({\mathbf{x}})/\sigma_{j}^{2}$ term can be added to the evolution equation $d\psi_{j}({\mathbf{x}})/dt$ for each component $j$. This is much preferable to trying to set the magnetization or relative populations of different components using Lagrange multipliers $\mu_{j}$ because the latter only set the mean particle number and may allow very large relative fluctuations from shot to shot. This is a work in progress.

Looking further ahead, the equations presented here should be capable of producing ensembles of attractive gases with $g<0$, something that absolutely cannot be stably treated using the standard SGPE. Moreover, they should also be easily extensible to the case of long range interactions, another situation when ensemble equivalence is known to be broken Touchette (2015).

Recent experiments have demonstrated atom number measurements well below atomic shot noise using dispersive imaging Gajdacz et al. (2016); Kristensen et al. (2017), which suggests that high precision studies of ensemble equivalence could be carried out with present technology.

Acknowledgements.

This work was supported by the National Science Centre grant No. 2012/07/E/ST2/01389.

Appendix A The case of nonuniform $\gamma({\mathbf{x}})$

If we do not assume (30), then instead of (31) one obtains the following even more cumbersome form of the FPE:

[TABLE]

We omitted the $\psi$ dependence of $N(\psi)$ for a minor improvement in brevity. Due to the fact that in many terms there is no summation index that involves only ${\mathcal{P}}$ and $\psi$, one cannot apply (32) in all necessary cases as was done for a constant $\gamma$.

A special but very common case is when the projection is made implicitly by the numerical lattice as done in the plain SGPE approach (rather than the SPGPE). Then, ${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}=\delta_{{\mathbf{x}}{\mathbf{y}}}$ and all the inconvenient features of (A) abate. One obtains

[TABLE]

which is similar in complexity to (II.5). This still complicated expression can be made zero with the simple choice

[TABLE]

in full space-dependent analogy to (34).

Now, if we return to the general projected case, the typical situation is that $\gamma({\mathbf{x}})$ is slowly varying in space compared to the highest energy modes in ${\mathcal{C}}$. These produce features of length $\lambda_{E}$. On the other hand, ${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}$ is typically close to diagonal with a width given by the length scale of the highest-energy components. This means that it decays to zero on length scales of the order of $\lambda_{E}$. Thus, as long as $\gamma({\mathbf{x}})$ varies slowly in the region of space around ${\mathbf{x}}$ compared to ${\mathcal{P}}_{{\mathbf{x}}{\mathbf{y}}}$, one will have

[TABLE]

This condition allows us to put the troublesome terms in (A) involving $\gamma$ into a form in which the projection property of the field, (32), can be applied. For example, one has

[TABLE]

However, there are some remaining (also troublesome) terms in (A) which involve $K_{{\mathbf{x}}}$ not $\gamma_{{\mathbf{x}}}$. Judging by the earlier result (50) for a special case, the rate at which $K$ will vary spatially is similar to that of $\gamma$. So, let us also provisionally assume the same slowly varying property for $K$, and check its consistency later. This assumption lets us apply

[TABLE]

Conditions (51)-(53) lead to much simplification in (A):

[TABLE]

The above equation becomes stationary using the same simple expression (50) as for the unprojected case. This confirms the validity of the condition (53) once (51) is assumed.

Substituting (50) into the general postulated stochastic equation (21) gives us the most general transition SPGPE (38).

Appendix B Some exact results for the ideal gas

We follow the same approach as Witkowska et al. (2009) used for the 1d trapped gas, but adapt the procedure for the doubly degenerate levels that occur in the uniform gas.

We have plane wave modes $k_{j}$ with energies

[TABLE]

occupied by $n_{j}$ bosons, and temperature set by $\beta=1/k_{B}T$. The total energy of a state is $E=\sum_{j}\varepsilon_{j}n_{j}$, and the number of atoms is $N=\sum_{j}n_{j}$.

B.1 Canonical ensemble

The fuller version of (1) from Witkowska et al. (2009) that also includes the canonical ensemble constraint $N={\overline{N}}$ is

[TABLE]

$j$ enumerates modes over all integers, including the ground state $j=0$ which is the only nondegenerate mode and sets the energy zero: $\varepsilon_{0}=0$. Normalization of $P$ is ignored. Combining the two deltas immediately implies the obvious $n_{0}={\overline{N}}-N_{\rm ex}$, and the necessity of the “cliff”, i.e. $P(N_{\rm ex}>{\overline{N}})=0$.

Substituting

[TABLE]

into (56) we get

[TABLE]

The sum has the form $\sum_{n=0}^{\infty}f^{n}=\frac{1}{1-f}$ where $f=e^{-\beta\varepsilon_{j}-ix}$. So

[TABLE]

We now have a right-hand contour on the unit circle in the variable $z=e^{ix}$. In particular, $dz=izdx$, so

[TABLE]

The poles (doubly degenerate) are at locations

[TABLE]

an all within the contour because $\beta\varepsilon_{k}>0$. Note that $m\geq 0$ now counts only energy levels, not modes, and we denote $m=0$ to be the ground state. With the help of the Cauchy residue theorem, the result is

[TABLE]

Upon substitution, the sum is now over excited energy levels, and the end result looks like this:

[TABLE]

B.2 Classical field expression

In the c-field approximation, again analogously to Witkowska et al. (2009), the expression corresponding to (56) is

[TABLE]

with Dirac deltas and mode amplitudes $\alpha_{j}$. Here also the two deltas give a deterministic condition $\delta\left(|\alpha_{0}|^{2}-({\overline{N}}-N_{\rm ex})\,\right)$ on the ground state amplitude, and the “cliff” is present as well. Moving on to the evaluation of this expression for degenerate states, we use (57) again and find

[TABLE]

The integrals are easily done, giving

[TABLE]

Changing to contour variable $z$,

[TABLE]

The poles are at the same locations (61) as before in the quantum case, with the same degeneracy, so that evaluation of the residues leads to

[TABLE]

This sum is overall similar to the quantum one (63), replacing the Bose-Einstein occupations with Rayleigh-Jeans, working with only the states in the subspace ${\mathcal{C}}$, and without the extra $+1$ term in the prefactor.

B.3 Grand canonical ensemble

Consider now the GCE in the quantum case. Here, weights for states are

[TABLE]

with $\mu<0$. We have, rather similarly to (56), that

[TABLE]

with no constraint on $N$. In principle the $j=0$ state enters above, but only as a prefactor $\sum_{n_{0}=0}^{\infty}e^{\beta\mu n_{0}}=1/(1-e^{\beta\mu})$ that can be incorporated into the normalization. Hence $\prod_{j}\to\prod_{j\neq 0}$ and the expression for $P(N_{\rm ex})$ differs from the expression (56) for the canonical ensemble only by the replacement $\varepsilon_{j}\to(\varepsilon_{j}-\mu)$. Proceeding as before, one obtains

[TABLE]

and for the c-field case:

[TABLE]