Supplementary information to "Efimov-driven phase transitions of the unitary Bose gas"
EEfimov-driven phase transitions of the unitary Bose gas
Supplementary information
Swann Piatecki and Werner Krauth
Laboratoire de Physique Statistique, ´Ecole Normale Sup´erieure, UPMC, Universit´eParis Diderot, CNRS, 24 rue Lhomond, 75005 Paris, France
The pair interaction in the model hamiltonian may be viewed as the zero-range limit of the square-wellinteraction potential shown in Supplementary Fig. 1, whose range r and depth V are simultaneously takento 0 and ∞ so that the scattering length a = r (cid:104) − tan( r (cid:112) mV / (cid:126) ) / ( r (cid:112) mV / (cid:126) ) (cid:105) remains constant. Atzero range, the pair interaction only a ff ects s -waves by the so-called Bethe-Peierls limit condition when theseparation r between two particles goes to zero.Unitarity corresponds to the pair interaction having a bound dimer of zero energy and infinite extension,a limit in which the k -vector of the particle is much smaller than 1 / r , and the scattering cross-section σ saturates the unitary limit σ < π/ k . This follows from the unitarity of the scattering operator.0 r Separation r − V V ( r ) r Separation r − V r Separation r − V R / a ∼ − (left) , in the unitary region R / a = (centre) , and in the bound pair + particle region R / a ∼ (right) .The zero-range potential is fixed- a limit of these potentials when V → ∞ and r →
0. The (possible) pairbound state correspond to cyan levels.While the two-particle properties are universal, the three-boson trimer ground state at unitarity generallydepends on the details of the pair interaction. Excited trimers form a geometric sequence of asymptoticallyuniversal Efimov trimers with a fixed energetic ratio E n / E n + → n →∞ .
03, where E n is the energy of the n -th excited trimer. For the zero-range pair interaction and the three-body hard core, the ground-state trimeris almost identical to universal Efimov trimers[1]. In our Path-Integral Monte Carlo algorithm, statistical weights due to the two-body and three-body inter-actions in the model hamiltonian are computed as follows: • the zero-range interaction is implemented through the pair-product approximation for the densitymatrix, that estimates the statistical weight of a configuration from the isolated two-body wavefunc-tion of nearby particles. For the zero-range interaction, only s -waves di ff er from the isolated system1 a r X i v : . [ c ond - m a t . qu a n t - g a s ] A ug f two non-interacting bosons. The pair-product approximation is valid when the imaginary timediscretization step ∆ τ is small enough. • The hyperradial cuto ff V is enforced by rejecting the configurations where the hyperradius is smallerthan a threshold R , an approximation also valid if ∆ τ is small enough. In practice, for the valuesof ∆ τ retained in our simulations, we need to take in account a finite ∆ τ shift of the input valueof R , which we compute from a fit of the three-particle hyperradial probability distribution to itsanalytically-known value[1], p ( R ) ∝ RK is (cid:16) √ κ R (cid:17) , where K is a Bessel function of imaginary index is ≈ . i , and κ = mE t / (cid:126) ∝ R − .These weights set the probability of acceptance of the following four types of Markov-chain moves,following a Metropolis-Hastings procedure.(1) Rebuild a single-particle path on a fraction of the total imaginary time,(2) Move a whole permutation cycle,(3) Exchange two strands of paths on a fraction of the total imaginary time,(4) Perform a compression-dilation move on one imaginary time slice (see Supplementary Fig. 2). r Positions I m a g i n a r y ti m e τ r ( + h ) Positions
DilationCompression
Supplementary Figure 2: Compression-dilation move (move (4)). At one imaginary time slice, the newproposed positions for two particles (blue paths) are either more distant or closer by a factor 1 + h whilethe centre of mass (white circles) and positions at other time slices are kept in place.Updates (1-3) are commonly used in Path-Integral Monte Carlo simulations[2, 3, 4]. The compression-dilation move (4) specifically addresses the divergence of the pair-distance distribution function (see Fig.1 d ):Even when it rebuilds a path on one single imaginary time slice, the update (1) is linear: if two particles areseparated by a distance r , it will typically propose a new configuration with a separation r + ∆ . The weightof the proposed configuration w ( r + ∆ ) is quite di ff erent from the weight of the original configuration w ( r )when r →
0, which results in a very low acceptance rate. The update (4) proposes a configuration r (1 + h ),whose weight w ( r ) + hrw (cid:48) ( r ) is close to w ( r ): Because w ( r ) diverges as 1 / r , rw (cid:48) ( r ) is of the same order as w ( r ), which generates much higher acceptance rates. In our simulations of unitary bosons, the system is contained in a harmonic trap, which regulates the avail-able configuration space. To reduce entropic e ff ects for the three-body calculations in Fig.1, we regulatethe volume by imposing that they lie on a single permutation cycle: In Fig.1 a-c , at imaginary time β , the2lue, red and green bosons are respectively exchanged with the red, green and blue bosons. This conditiondoes not modify the properties of the fundamental trimers as other permutations could be sampled at nocost at points of close encounters, some of which are highlighted in Fig.1 b .In Fig.1 a-c , four-dimensional co-cyclic path-integral configurations are represented in a three-dimensionalplot. As the motion of the centre of mass is decoupled from the e ff ects of the interactions, its position is setto zero at all τ . The three spatial dimensions are then reduced to two dimensions by rotating the triangleformed by the three particles at each imaginary time to the same plane in a way that does not favour anyof the three spatial dimensions, but that conserves the permutation cycle structure. This transformationconserves the pair distances. The local density approximation consists in assuming that, at each point x = ( x , y , z ) of the trap, with | x | = r , the Bose gas is in equilibrium at a chemical potential µ = µ − ω r /
2, where µ is the chemicalpotential in the trap centre. Within this approximation, thermodynamical identities allow one to relate thedoubly-integrated density profile ¯ n ( x ), and the grand-canonical pressure[5]: P ( µ ( x )) = ω π ¯ n ( x ) , (1)where µ ( x ) = µ − ω x /
2. The chemical potential µ at the centre of the trap can be obtained from a fit inthe outer region of the trap where P ( µ ) is the pressure of a classical ideal gas.The doubly-integrated density profile is obtained by ensemble averaging, and the numerical equationof state P ( µ ) is compared to the expansion of the pressure in powers of the fugacity e βµ (see Fig.3), P = k B T λ (cid:88) l ≥ b l e l βµ . (2)The l -th cluster integral b l corresponds to the l -body e ff ects that cannot be reduced to smaller non-interacting groups of interacting particles[6]. The classical ideal gas corresponds to b =
1, and highercoe ffi cients both describe both energetic and quantum-statistical correlations. The l -th cluster integrals arerelated to the virial coe ffi cients { a l (cid:48) } l (cid:48) ≤ l [6] as b = a = , b = − a , b = − a + a . (3)At unitarity, an analytical expression of the coe ffi cients b and b was obtained recently[7], b = / (4 √ ≈ .
59, and b √ = C + (cid:88) q ≥ (cid:16) e − β(cid:15) q − (cid:17) + | s | π
12 ln ( e γ β E t ) − (cid:88) p ≥ e − p π | s | (cid:60) (cid:104) Γ ( − ip | s | )( β E t ) ip | s | (cid:105) , (4)where C = .
648 and | s | = . γ = .
577 is the Euler constant, (cid:15) q = − E t e − π q / | s | represents the energy levels of the three-body bound states, and − E t is the energy of the fundamentalEfimov trimer, which is related to the average squared hyperradius in the fundamental state[1]: E t = + s )3 (cid:10) R (cid:11) , (5)an expression that may be used to obtain E t ≈ . e − (cid:126) / ( mR ) by numerically integrating Eq. 2.In the temperature range of our simulations, it is su ffi cient to perform the sums in Eq. 4 up to q = Monitoring the phase transition
When the densities of the gas and the superfluid Efimov liquid approach each other, observing directlythe two-dimensional histogram of pair distances and centre-of-mass positions does not allow to distinguishbetween a weakly first-order phase transition and a cross-over (see Supplementary Fig. 3).0 . . . r / R r s e p / R a . . . r / R b . . . r / R c Supplementary Figure 3: Two-dimensional histogram of pair distances r sep and centre-of-mass positions ¯ r for R = .
3, at k B T / (cid:126) ω = . . (a), k B T / (cid:126) ω = . k B T / (cid:126) ω = .
0. The sole observation of thishistogram does not allow here to state that the system is undergoing a phase transition as the densities ofboth phases are too close.In this regime, we monitor the normal-gas-to-superfluid-liquid phase transition more accurately byfollowing the evolution of the first peak of the pair correlation function (obtained by ensemble averaging)with temperature (see Supplementary Fig. 4).1 . . T / T r p ea k / R R / a ω = .
07 1 . . T / T R / a ω = .
13 0 . . T / T . . . R / a ω = . R . For R / a ω = .
07, the large step indicates the first order phase transition (thiscorresponds to Fig.3 a-c ); at R / a ω = .
13, the step persists but is considerably smaller, a sign that thesystem also undergoes a first order phase transition. At R / a ω = .
23, there is no sign of a first order phasetransition.For each value of R , simulations in the trap are run at a discrete set of inverse temperature β . InFig.4 a , the error bars in the normal-gas-to-superfluid-liquid transition line show the interval between thelast temperature at which there is no liquid droplet, and the first temperature at which there is one. Theerrors for the conventional superfluid transition are computed in the same way, with the criterion that the gasis superfluid if particles lie in a permutation cycle of length at least ten with a probability higher than 0 . Approximate semi-analytical phase diagram
When the free energy F of a homogeneous physical system is not a convex function of its volume V , itbecomes more favourable to split the system into two phases than to keep the system homogeneous, asituation from which first order phase transitions originate (see Supplementary Fig. 5), and that yields theequality of the pressures and of the chemical potentials of both phases at coexistence in absence of interfaceenergy[8]. v l ∝ R v g V / Nf l ( v l ) − ε F / N CoexistenceIncompressibleliquid Normal gas(virialexpansion) v l ∝ R v g V / Np vap P ≡− ∂ V F Coexistence
Supplementary Figure 5: Sketch of the behaviour of the free energy and the pressure of unitary bosonsthrough the transition in our incompressible fluid model.In the system of unitary bosons, the virial expansion is an excellent approximation for the normal gasphase far from the superfluid transition. Although it becomes irrelevant in the superfluid gas, its analyticcontinuation conveys important qualitative features, because of the continuous nature of the normal-gas-to-superfluid-gas phase transition, and is therefore a suitable poor man’s approximation to the behaviour ofthe superfluid gas.The simplest theoretical model for the superfluid liquid is that of an incompressible liquid of specificvolume v l and negligible entropic contribution to the free energy. The incompressibility approximation isacceptable for unitary bosons as the first peak of the pair correlation function seems to scale only with R (see Fig.2); simulations at high R yield v l ≈ (6 R ) . The negligibility of the entropic contribution tothe free energy is ensured for non-pathological systems at low temperature. As results of Ref.[9] may beextrapolated to obtain an energy per particle − (cid:15) = − . E t in the liquid phase (see Supplementary Fig. 6),in this approximation, the free energy is F l = − N (cid:15) .In practice, we draw the transition line into the superfluid liquid by finding the smallest chemicalpotential at which the pressures of the incompressible liquid and of the gas phase (approximated by thevirial expansion) coincide: µ − (cid:15) v l = k B T λ (cid:16) e βµ + b e βµ + b e βµ (cid:17) . (6)As 1 / v = ∂ µ P , the crossing to the regime where this equation has no solution corresponds to the criticalpoint, where both densities are equal.To draw the coexistence line for the trap centre with N =
100 particles, the chemical potential µ at thecentre of the trap is found from integrating the density throughout the trap:4 πλ (cid:90) ∞ r dr (cid:16) e βµ ( r ) + b e βµ ( r ) + b e βµ ( r ) (cid:17) = N , (7)where µ ( r ) = µ − ω r /
2. 5 5 10 15 N | E | / E t Supplementary Figure 6: Energy of clusters of unitary bosons at zero temperature (black crosses) as afunction of the number N of bosons in the cluster, computed by J. von Stecher in Ref.[9], in units of thefundamental Efimov trimer energy E t . A linear regression (solid red line) gives | E | / N = . E t for large N values. For the phase diagram of 100 bosons in a harmonic trap (Fig.4 a ), the value of | E | / N is empiricallyadjusted to 8 E t to account for the small number of bosons in the liquid droplet. References [1] Braaten, E. & Hammer, H.-W. Universality in few-body systems with large scattering length.
Phys.Rep. , 259–390 (2006).[2] Pollock, E. L. & Ceperley, D. M. Simulation of quantum many-body systems by path-integral methods.
Phys. Rev. B , 2555–2568 (1984).[3] Krauth, W. Quantum Monte Carlo calculations for a large number of bosons in a harmonic trap. Phys.Rev. Lett. , 3695–3699 (1996).[4] Krauth, W. Statistical Mechanics: Algorithms and Computations (Oxford University Press, Oxford,Great Britain, 2006).[5] Ho, T.-L. & Zhou, Q. Obtaining the phase diagram and thermodynamic quantities of bulk systemsfrom the densities of trapped gases.
Nature Phys. , 131–134 (2010).[6] Huang, K. Statistical mechanics (Wiley, New York, 1987), 2nd edn.[7] Castin, Y. & Werner, F. Le troisi`eme coe ffi cient du viriel du gaz de bose unitaire. Canad. J. Phys. ,382–389 (2013). arxiv:1212/5512(Englishversion) .[8] Landau, L. D. & Lifshitz, L. M. Statistical physics . No. 5 in Course of theoretical physics (Butterworth-Heinemann, 1980), 3rd edn.[9] von Stecher, J. Weakly bound cluster states of Efimov character.
J. Phys. B43