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 [2,9,14] and optical lattice geometries [5,21,26], including gases at temperatures much lower than the degeneracy temperature.
The effective field theory of a cold gas of neutral bosonic atoms with shortrange 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 [13]) whereĤ 0 andĤ int are respectively the free-particle and interaction Hamiltonians,ψ(z) is the field operator, which annihilates a particle at position z, m is the atomic mass, V (z) is the external trap potential and g is the effective interaction strength, given in the experimentally relevant case of a harmonic transversal confinement with trapping frequency ω r by g = 2ω r a s , with a s being the s-wave scattering length. The usual experimental setups deal with thousands of atoms [2], 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ψ(z) by a classical field ψ(z) [4]. This approach is valid for low temperatures, where we have a range of macroscopically occupied modes |ψ n | 2 ≫ 1; the modesψ n are taken to be the eigenfunctions of the one-body non-interacting HamiltonianĤ 0 . The evolution of this redefined classical order parameter ψ(z) is then governed by the celebrated Gross-Pitaevskii equation (GPE) [8].
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 ψ 0 (z) 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 [20].
Conventional methods of initial state sampling include analytical ones [23,12], where the gas is initialized with a Bose-Einstein distribution of Bogoliubov quasiparticles with random phases, as well as stochastic ones [7,10], 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 [16]. For a broad overview of quantum and classical Monte Carlo methods, including the Metropolis-Hastings algorithm, see [18,1] 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 rel-evant 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 [3], 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: 1. Initialization: (a) Choose an initial order parameter ψ 0 (z). Specific choices of ψ 0 (z) will be discussed in the following section.
where β is the inverse temperature, µ is the chemical potential (both β and µ are fixed external parameters), andN = ψ † (z)ψ(z) dz is the particle number operator. Note that the free energy does not enter the expression for S 0 , 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 its absolute value.

For each iteration
(a) Generate a candidate field ψ N (z) by varying the energy. This variation of energy can be achieved by adding either a density phonon or a phase phonon 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.
(b) Vary the particle number (c) Calculate the reduced entropy of the candidate field Z e SN is the Boltzmann probability to find the field in the state ψ N (z). 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 Z, which is practically impossible for interacting systems with many degrees of freedom. Then we check the value of a: Real random numbers, distributed normally with zero mean and unit variance.
Numerical constants governing the rate of convergence to the equilibrium state. In the presented results they have been empirically chosen to be c 1 = 4(n 0 ) −1/2 , 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/ , wheren 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.
i. If a ≥ 1, then the candidate state is more probable than the seeding state, so we keep the former. ii. If a < 1, we pick a uniform random number r ∈ [0, 1]. If r ≤ a, the candidate state is accepted; but if r > a, the candidate state is discarded and the seeding state is kept for the next iteration ψ N (z) := ψ N −1 (z).
(e) Proceed to the next iteration.
As a result we have a Markov chain of states ψ N (z), N ∈ [0, N max ], 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. Neighbouring states ψ N and ψ N +1 are usually highly correlated (as they differ by only one phonon and a small particle number variation), so it is necessary to throw away majority of the results, picking only one state out of N c , where N c 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: 1. Arbitrary trap geometry, as we can freely modify the trapping potential V in the total HamiltonianĤ. In general the phonons in Eqs. 5 and 6 can be modified to be the eigenfunctions of the trapping potential (e.g. in the case of harmonic confinement V (z) ∝ z 2 we can take Hermite functions 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.
2. Any number of dimensions. This requires representing the order parameter as a scalar field on many-dimensional space ψ( z), the phonons (Eqs. 5 and 6) being modified accordingly as sin( k r z + φ r ).
3. Canonical state sampling. Reduced entropy becomes S N = −β ψ N |Ĥ|ψ N , and we have to omit the 2b stage of the algorithm to make sure the particle number does not change.
4. Generalized Gibbs ensemble sampling. Reduced entropy now reads whereÎ i are the local conserved charges (integrals of motion) of the system, in addition to the energy Ĥ and the particle number N , and µ i 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 [25]. 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Ĥ ′ =Ĥ − 1 β i µ iÎi , 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 87 Rb atoms and an experimentally relevant case of the same gas in a harmonic confinement. The parameters of the simulation are summarized in table 2. Typical examples of the grand canonical thermal state of the 1D Bose-Einstein quasicondensate after N max = 10 5 Metropolis-Hastings iterations are presented in figure 1.
The initial state for all the presented results was taken to be the ground state of the non-interacting gas (n 0 (z) = n 0 = const, φ 0 (z) = 0) 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 g 1 autocorrelation function  Figure 1: Typical examples of the grand canonical thermal state with the temperature T = 10 nK of the interacting 1D BEC, achieved after N max = 10 5 Metropolis-Hastings iterations in the untrapped system with periodic boundary conditions (a, b; thin zigzag line) and harmonically trapped case (c, d; thin zigzag line). Quasicondensate local densities n(z) (a, c), measured in atoms per micrometer, and phases φ(z) (b, d), measured in units of π, as a function of the longitudinal direction z in micrometers. Thick horizontal (a, b, d) and parabolic (c) lines represent the initial conditions, which in the case of the untrapped system (a, b) were taken to be the ground state of the non-interacting gas (n 0 (z) = n 0 , φ 0 (z) = 0), and in the case of the harmonic confinement (c, d) as a Thomas-Fermi parabolic density profile with constant zero phase. Extensive fluctuations of the phase at the edges of (d) 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.  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 N = 2 − 6 · 10 4 iterations.  . 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. 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. (b) Natural logarithm of the g 1 correlation function in the homogeneous case for the temperatures T = 10, 60 and 120 nK (from top to bottom) at the last iteration N max of the algorithm, averaged over the ensemble of 70 realizations. These g 1 functions were used to calculate averaged temperatures presented in figure 2(a). The linear region of the logarithm spans from 0 till ≈ 15 µm, 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 L = 200 µm) and are to be discarded.
In thermodynamic equilibrium at positive temperatures in 1D g 1 is exponentially decaying with ∆z, confirming the fact that there can be no true Bose-Einstein condensate in this case where λ T is thermal coherence length with k B being the Boltzman's constant andn = 1 L ′ |ψ(z)| 2 dz the mean density of the cloud. L ′ is the averaging length, which is the length of the integration region in Equation 9 as well. In case of untrapped gas gas L ′ = L is the total simulation region, but in case of harmonic confinement the integration region contains only the points where the local density n(z) 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 figure 2, with one particular example of the g 1 function in figure 3(b). It is evident that the thermal equilibrium is achieved after N = (2 − 6) · 10 4 iterations.
Another independent test whether the achieved state is thermal is the realtime development of the state using Gross-Pitaevskii equation, as by definition the thermal state should remain thermal during such evolution. The results for the untrapped gas, presented in the inset to the figure 3(a), 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 (there exist efficient algorithms for solving real-time GPE, see e.g. [24,17]).
As in all realization of Metropolis-Hastings algorithm a 'good guess' of the initial state is essential for the fast convergence. In figure 3(a) 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 [12,23]. 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) and doesn't immediately lead to the desired thermal equilibrium. This is one particular example where numerical methods successfully compete with 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 figure 4, where C ψ is the two-point correlation function of the last sample ψ Nmax (z) and C n is the density fluctuation correlation function of the last sample where δn(z) = n(z) −n, n(z) = |ψ(z)| 2 , andn = 1 L ′ n(z) dz for the uniform gas.
It is evident that the order parameters still remain phase-correlated after 10 5 iterations, which is consequence of the fact that we observe the system below the thermal gas to quasicondensate crossover temperature [6]: 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 C n 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 N c = (2 − 8) · 10 4 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'.

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 [19], prethermalization of a split quasicondensate [11]), it can be added to the thermal state using conventional methods [4,20].
The proposed algorithm can be generalized to higher dimensions and arbi-trary trap geometry. We see the main advantage of the proposed 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 [22,15].