Metropolis-Hastings thermal state sampling for numerical simulations of Bose-Einstein condensates.

We demonstrate the application of the Metropolis-Hastings algorithm to sampling of classical thermal states of one-dimensional Bose-Einstein quasicondensates in the classical fields approximation, both in untrapped and harmonically trapped case. The presented algorithm can be easily generalized to higher dimensions and arbitrary trap geometry. For truncated Wigner simulations the quantum noise can be added with conventional methods (half a quantum of energy in every mode). The advantage of the presented method over the usual analytical and stochastic ones lies in its ability to sample not only from canonical and grand canonical distributions, but also from the generalized Gibbs ensemble, which can help to shed new light on thermodynamics of integrable systems.

Introduction

The recent advances in experimental methods allowed precise control and manipulation of ultracold atoms in various trap [1-3] and optical lattice geometries  [4-6], including gases at temperatures much lower than the degeneracy temperature.

The effective field theory of a cold gas of neutral bosonic atoms with short-range repulsive interactions is given by the second quantized Hamiltonian (in the following we deal explicitly with a one-dimensional (1D) case, where quasicondensation takes place instead of true condensation  [7]) where and are respectively the free-particle and interaction Hamiltonians, is the field operator, which annihilates a particle at position , is the atomic mass, is the external trap potential and is the effective interaction strength, given in the experimentally relevant case of a harmonic transversal confinement with trapping frequency by , with being the s-wave scattering length.

The usual experimental setups deal with thousands of atoms  [1], so the quantum dynamics can be numerically simulated only using various approximations. The one approximation especially suited for studies of weakly interacting cold atomic gases is the classical field approximation, where we replace the quantum field operator of the effective field theory by a classical field   [8]. This approach is valid for low temperatures, where we have a range of macroscopically occupied modes ; the operators are defined through the normalized eigenfunctions of the one-body non-interacting Hamiltonian . The evolution of this redefined classical order parameter is then governed by the celebrated Gross–Pitaevskii equation (GPE)  [9].

In experiments with cold atomic gases the system is usually prepared in thermal equilibrium, before a quench or another manipulation is applied, therefore the numerical methods for sampling the thermal initial condition are of great importance. The quantum correction for the classical thermal state of a weakly interacting system can be introduced using the so-called truncated Wigner approximation (TWA), where zero-point quantum oscillations are taken into account in the initial state only, but the subsequent evolution is classical  [10].

Conventional methods of initial state sampling include analytical ones  [11,12], where the gas is initialized with a Bose–Einstein distribution of Bogoliubov quasiparticles with random phases, as well as stochastic ones  [13,14], where the thermal state is achieved during imaginary time GPE evolution with Langevin noise.

In the present paper we propose another way of sampling the initial distribution, namely using the Metropolis–Hastings algorithm. We believe that in some cases it might be preferable over the analytical methods, as it does not use Bogoliubov-type approximations, and may be used to sample states out of a generalized Gibbs ensemble, which is impossible with existing stochastic realizations.

Metropolis–Hastings algorithm

The Metropolis–Hastings algorithm is a Markov chain Monte Carlo method for sampling a probability distribution, especially suited for systems with many degrees of freedom  [15]. For a broad overview of quantum and classical Monte Carlo methods, including the Metropolis–Hastings algorithm, see  [16,17] and references therein.

In the present paper we demonstrate the implementation of the Metropolis–Hastings method for 1D Bose–Einstein quasicondensate without confinement (implying periodic boundary conditions) as well as for the experimentally relevant case of a harmonic longitudinal confinement. The method can be easily generalized to higher dimensions and other trap geometries.

This method has been already applied to classical simulations of cold Bose gases  [18], but it has not been explicitly formulated as a step-by-step algorithm. In the present paper we systematically study the convergence properties of this method and outline its application to sampling the generalized Gibbs ensemble (GGE).

In our particular realization the algorithm reads as follows:

Initialization:

Choose an initial order parameter . Specific choices of will be discussed in the following section.

Calculate the reduced entropy , where is the inverse temperature, is the chemical potential (both and are fixed external parameters), and is the particle number operator. Note that the free energy does not enter the expression for , meaning that the zero-level of the latter is not defined. This is justified by the fact that we are interested only in differences of the reduced entropy, and not in its absolute value.

For each iteration :

Generate a candidate field by varying the energy. This variation of energy can be achieved by adding either a density perturbation or a phase perturbation to the field from the previous iteration (‘the seeding field’). Whether to choose the one or the other is decided at random (by a ‘coin toss’). The meaning and values of the parameters are summarized in Table 1.

Table 1

Numerical parameters of the Metropolis–Hastings algorithm.

Parameter	Description
v_r,u_r	Real random numbers, distributed normally with zero mean and unit variance.
c_1, c_2, c_3	Numerical constants governing the rate of convergence to the equilibrium state. In the presented results they have been empirically chosen to be c_1=2(n_0)^−1,c_2=0.1 and c_3=0.001, where n_0=max|ψ_0|^2 is the maximal initial density. This particular choice provided typical values of the acceptance ratio in each iteration a∈[0.4,0.6], which gave the fastest convergence to equilibrium. It was numerically checked that different choices of those constants did not affect the resulting state, only the rate of convergence.
ϕ_r	Random phase ϕ_r∈[0,2π) picked from the uniform distribution.
k_r	Random wave number picked from the set {±δk,±2δk,…,±k_max}, where δk=2π/L, L is the length of the simulated region and k_max is the cutoff wave number. It was numerically checked that the results do not depend on this cutoff as long as it is larger than the inverse healing length ξ^−1=mgn¯/ħ^2, where n¯ is the mean density. So we present results where k_max=N_zδk/2 is the maximal allowed wave number on a lattice of N_z sampling points.

Vary the particle number

Calculate the reduced entropy of the candidate field

Calculate the acceptance ratio , where is the Boltzmann probability to find the field in the state . The main advantage of the Metropolis–Hastings algorithm lies the fact that we have to evaluate only the ratio of probabilities, in this way avoiding to calculate the partition function , which is practically impossible for interacting systems with many degrees of freedom. Then we check the value of :

If , then the candidate state is more probable than the seeding state, so we keep the former.

If , we pick a uniform random number . If , the candidate state is accepted; but if , the candidate state is discarded and the seeding state is kept for the next iteration .

Proceed to the next iteration.

As a result we have a Markov chain of states , , which can be used as thermal initial states for classical fields simulations. It is important to throw away the states obtained at early iterations (so-called ‘burn-in’ period), where the thermal state is not yet achieved. Neighboring states and are usually highly correlated (as they differ by only one elementary perturbation), so it is necessary to throw away majority of the results, picking only one state out of , where is calculated from the iteration-to-iteration autocorrelation length. We will return to these two issues in the results section.

Straightforward generalizations of the algorithm are easily conceivable:

Arbitrary trap geometry, as we can freely modify the trapping potential in the total Hamiltonian . In general the perturbations in Eqs. (5) and (6) can be modified to be the eigenfunctions of the trapping potential (e.g. in the case of harmonic confinement we can take harmonic oscillator eigenfunctions instead of sine-waves). But in practice using potential-specific eigenfunctions instead of plane waves did not give any speed-up to the achievement of the steady state, so the algorithm can be used without this modification.

Any number of dimensions. This requires representing the order parameter as a scalar field on many-dimensional space , the perturbations (Eqs. (5) and (6)) being modified accordingly as .

Canonical state sampling. Reduced entropy becomes , and we have to omit the 2b stage of the algorithm to make sure the particle number does not change.

Generalized Gibbs ensemble sampling. Reduced entropy now reads where are the local conserved charges (integrals of motion) of the system, in addition to the energy and the particle number , and are generalized potentials. For instance, in case of 1D GPE there exists an infinite number of local conserved charges, which can be explicitly calculated using Zakharov–Shabat construction  [19]. We regard this possibility of GGE sampling as the primary advantage of the presented algorithm. In fact, simulation of the GGE requires only redefinition of the Hamiltonian to , to which the previously described algorithm can be applied without further modification. We reserve the detailed analysis of this case for a separate publication.

Results

In the following we demonstrate the application of the algorithm to generate a grand canonical thermal state for an untrapped gas of neutral 87Rb atoms and an experimentally relevant case of the same gas in a harmonic confinement. The parameters of the simulation are summarized in Table 2.

Table 2

Simulation parameters of the systems presented in the results section.

Parameter	Description
m=87⋅1.67⋅10^−24g	Atomic mass of 87Rb atoms
a_s=5.3⋅10^−7cm	s-wave scattering length
T=10,60 or 120 nK	Temperature
ω_r=2π⋅2000s^−1	Transversal trapping frequency
n_0=90atoms/μm	Maximal initial linear atom density of the cloud
g=2ω_ra_sħ	1D interaction strength
μ=gn_0	Chemical potential
ω_l=2π⋅10s^−1	Longitudinal trapping frequency in case of a harmonic confinement
L=200μm	Total length of the simulation region
N_z=512	Number of spatial discretization points, so the state ψ(z) has N_z degrees of freedom
N_max=10^5	Total number of Metropolis–Hastings iterations
n(z)=|ψ(z)|^2	Local density
ϕ(z)=argψ(z)	Local phase

Typical examples of the grand canonical thermal state of the 1D Bose–Einstein quasicondensate after Metropolis–Hastings iterations are presented in Fig. 1.

Fig. 1

Typical examples of the grand canonical thermal state with the temperatures and 60 nK (labels on the panels) of the interacting 1D BEC, achieved after Metropolis–Hastings iterations in the untrapped system with periodic boundary conditions (four top panels) and harmonically trapped case (four bottom panels). Quasicondensate local densities (left), measured in atoms per micrometer, and phases (right), measured in radians, as a function of the longitudinal direction in micrometers. The initial conditions in the case of the untrapped system were taken to be the ground state of the non-interacting gas, and in the case of the harmonic confinement as a Thomas–Fermi parabolic density profile with constant zero phase. Note that achieved thermal state is not dependent on the initial conditions (see discussion in the text). Extensive fluctuations of the phase at the edges of harmonically trapped quasicondensate are due to the fact that the density there is close to zero, and the phase can take arbitrary values. Physical parameters of the simulations are summarized in Table 2.

The initial state for all the presented results was taken to be the ground state of the non-interacting gas (, ) in the untrapped case, and a Thomas–Fermi parabolic density profile with constant zero phase in the case of the harmonic confinement.

The achievement of the steady state is controlled by temperature measurement at the each iteration of the algorithm, calculated from the autocorrelation function

In thermodynamic equilibrium at positive temperatures in 1D is exponentially decaying with , confirming the fact that there can be no true Bose–Einstein condensate in this case where is thermal coherence length with being the Boltzmann constant and the mean density of the cloud. is the averaging length, which is the length of the integration region in Eq. (9) as well. In case of untrapped gas is the total simulation region, but in case of harmonic confinement the integration region contains only the points where the local density is larger than one tenth of the mean density. This helps to get rid of unessential boundary perturbations, probing the temperature of ‘the bulk’ of the condensate.

The Metropolis–Hastings ‘evolution’ of the temperature is presented in Fig. 3, with one particular example of the function in Fig. 2. It is evident that the thermal equilibrium is achieved after iterations.

Fig. 3

Temperatures during the Metropolis–Hastings ‘evolution’ as a function of the iteration number in the case of untrapped (a) and harmonically trapped (b) gas for three equilibrium temperatures (given as external parameters) and 120 nK. These temperatures are represented by three horizontal dashed lines serving as guides for the eye. Dots stand for one particular realization of the algorithm for the three temperatures (respectively, from bottom to top), and the corresponding solid lines show the averaged temperature over an ensemble of 70 realizations, each having the same initial conditions. Large temperature fluctuations in a single realization stem from the finite size of the simulation region, as they should converge to the equilibrium value only in thermodynamic limit. But from the ensemble averages it is evident that the thermal equilibrium is achieved after iterations.

Fig. 2

Natural logarithm of the correlation function in the homogeneous case for the temperatures and 120 nK (from top to bottom) at the last iteration of the algorithm, averaged over the ensemble of 70 realizations. These functions are used to calculate averaged temperatures presented in Fig. 3(a). The linear region of the logarithm spans from 0 till , and it is used in temperature measurement. The bending and fluctuations in the subsequent region are due to the finite size effects (as the total size of the system is ) and are to be discarded.

During the initial phases of this ‘evolution’ the system passes through a non-equilibrium region, so temperature in general sense in not well defined until the final equilibrium is reached. Though even in the non-equilibrium case the system always possesses a well-defined energy, so we can always formally define an emergent ‘non-equilibrium’ temperature as the temperature of a completely thermalized system with the given average energy. This emergent temperature is then measured by . We stress though that before the steady state is achieved the ‘temperature’ can be used as a convergence monitor only. Clearly, there is no real physical process behind the apparent ‘heating’ in Fig. 3 as there is no real time evolution. We also note that the temperature calculated from is independent of in Eq. (7), the latter being an external parameter, fixing our desired temperature.

Another independent test whether the achieved state is thermal is the real-time development of the state, as by definition the thermal state should remain thermal during such evolution. To check this criterion we prepared the thermal state of the untrapped gas with Metropolis algorithm and then propagated it in real time with Gross–Pitaevskii equation (there exist efficient algorithms for solving real-time GPE, see e.g.  [20,21]). The results, presented in the inset to Fig. 4, show that indeed the temperature of the state does not change on average, assuring that the initial state was thermal with respect to the Gross–Pitaevskii Hamiltonian.

Fig. 4

Influence of the initial state on the rate of convergence to the thermal state. Temperatures during the Metropolis–Hastings ‘evolution’ as a function of the iteration number in the case of untrapped gas for , averaged over 70 realizations. Thick line: zero-temperature state of the non-interacting gas, cf. Fig. 3(a). Thin line: thermal gas of Bogoliubov quasiparticles with random phases and constant amplitudes (see explanation in the text). Both choices of initial conditions eventually lead to equilibrium, but in case of the ‘Bogoliubov gas’ the convergence is faster, meaning that it is a better ‘initial guess’ for the thermal state. In this particular realization the temperature is rising during the ‘evolution’, but we note that if we had chosen a higher-than-desired initial temperature, then the temperature would be dropping to the desired value. Inset. Temperature of the state, produced by the real-time GPE evolution starting from the achieved thermal state as a function of time. Dots: one single realization, solid line: average over 70 realizations. The stability of the temperature shows that the initial state was indeed the thermal state of the Gross–Pitaevskii Hamiltonian (see the text for further discussion).

As in all realization of Metropolis–Hastings algorithm a ‘good guess’ of the initial state is essential for the fast convergence. In Fig. 4 we compare the beforementioned zero-temperature initial conditions with the initial state given by the thermal gas of Bogoliubov quasiparticles with random phases and constant amplitudes, given by the equilibrium Bose–Einstein distribution at the desired temperature of 60 nK  [11,12]. This initial condition seems to be a much better ‘initial guess’, leading to faster convergence. Note that the analytical method is only an approximation (implying weak interactions and neglecting the variance of the amplitudes of the quasiparticles), meaning that it gives a non-equilibrium state, which though is expected to be close to the desired thermal equilibrium: for instance, in Fig. 4 we see that the emergent ‘non-equilibrium’ temperature of this initial condition is about 48 nK, which is lower than the expected 60 nK. So this is one particular example where numerical methods are superior to the analytical ones.

It is well known that Markov chain methods give highly correlated samples from one iteration to the other. We present some correlators for the untrapped gas in Fig. 5, where is the two-point correlation function of the last sample and is the density fluctuation correlation function of the last sample where , , and for the uniform gas.

Fig. 5

(a) Order parameter correlation function for the untrapped gas as a function of the iteration number for temperatures 10, 60 and 120 nK (from top to bottom), averaged over 70 realizations. Remaining strong phase coherence after 105 iterations is due to the existence of long-range order in finite-size quasi-BEC. (b) Density fluctuation correlation function for the same realizations as in subfigure (a) for temperatures 10, 60 and 120 nK (from bottom to top).

It is evident that the order parameters still remain phase-correlated after 105 iterations, which is a consequence of the fact that we observe the system below the thermal gas to quasicondensate crossover temperature  [22]: the thermal fluctuations are too weak to randomize the overall phase (note that the effects of phase diffusion are absent as there is no real-time propagation).

This remaining phase correlation has to be taken into account when performing simulations involving two or more independently prepared condensates, where a random constant overall phase difference should be added to the initial conditions at each realization. For one condensate it is not necessary, as only the phase difference is observable, and not the phase itself.

Density fluctuation correlation function gives a better representation of the correlations in Metropolis–Hastings algorithm, and from the numerical simulations it follows that one should pick one state out of iterations (depending on the temperature) to assure statistical independence. It is always a safe choice to pick only one last realization out of the whole Markov chain, reinitializing the simulation for each ‘measurement’.

The proposed Metropolis–Hastings method is generally slower (requiring more CPU time) than the analytical method of thermal Bogoliubov gas  [11,12] or stochastic ones  [13,14]. For instance, propagation of steps to achieve the thermal states in Fig. 1 takes about 100 s on a 3.40 GHz Intel Core workstation, while a stochastic algorithm would require about 10 s and the analytical one would be instantaneous. Nevertheless, the proposed method is believed to be more precise in comparison with the analytical one (see Fig. 4 and discussion in the text) and more versatile in applications to generalized thermal ensembles in comparison to stochastic methods. In addition, it can be used as an independent benchmark for other numerical algorithms.

Conclusion

We developed an application of Metropolis–Hastings algorithm to sampling the classical thermal states of one-dimensional Bose–Einstein quasicondensates in classical field approximation in the case of untrapped gas with periodic boundary conditions and in experimentally relevant case of harmonic confinement. The achieved thermal steady state can be further used as an initial state for truncated Wigner simulations. In case when the quantum noise is important (e.g. collisions of condensates  [23], prethermalization of a split quasicondensate  [24]), it can be added to the thermal state using conventional methods  [8,10].

The proposed algorithm can be generalized to higher dimensions and arbitrary trap geometry. We see the main advantage of the method in its ability to sample not only the conventional thermodynamic ensembles, but also the generalized Gibbs ensemble, which is believed to arise in the integrable one-dimensional bosonic gas  [25,26].