Resonant Five-body Recombination in an Ultracold Gas of Bosonic Atoms

We combine theory and experiment to investigate five-body recombination in an ultracold gas of atomic cesium at negative scattering length. A refined theoretical model, in combination with extensive laboratory tunability of the interatomic interactions, enables the five-body resonant recombination rate to be calculated and measured. The position of the new observed recombination feature agrees with a recent theoretical prediction and supports the prediction of a family of universal cluster states at negative $a$ that are tied to an Efimov trimer.


Introduction
Few-body physics with ultracold atoms has emerged as a new research field combining concepts from atomic, nuclear and condensed-matter physics. A growing number of experimental and theoretical studies have been focused on both the fundamentals of few-body phenomena [1] and the connections with many-body systems [2,3]. The cornerstone of recent experimental advances is the control of interactions in an ultracold atomic gas, offered by magnetically tuned Feshbach resonances [4]. In particular, the tunable s-wave scattering length a allows to access the regime of resonant twobody interactions. Here, the system is governed by universal behavior, independent of the short-range details of the interaction potential. The paradigm of universality is Efimov's solution to the problem of three resonantly interacting particles [5]. Once the intimate connection between Efimov states and three-body recombination has been established [1,6,7], resonant loss features became the fingerprint of Efimov physics in experimental studies with ultracold atoms [8][9][10][11][12][13][14][15][16][17][18][19][20].
Advances in three-body physics led to intriguing questions on the generalization of Efimov's scenario to more particles and on the existence and observability of universal few-body cluster states [21][22][23][24][25][26]. N-body cluster states known as Brunnian states exist in a range of interaction where no (N − 1)-body weakly bound subsystems are present [27,28]. The general connection, however, between cluster states and the threebody Efimov states has remained an open issue. It was soon realized that no "true" Efimov states with N > 3 exist, because of the quite different scaling and threshold properties of the cluster states [29]. However, other approaches to extend universal Efimov theory to larger systems have been pursued along different lines [30][31][32]. The development of accurate descriptions of four-boson systems [33][34][35] demonstrates the existence of four-body states tied to each three-body Efimov state. In a major extension of Efimov physics, Ref. [36] predicts the existence of a family of cluster states tied to an Efimov trimer. According to Refs. [34][35][36][37], the binding energy of the cluster states follows universal scaling laws, which are directly connected to the Efimov effect. Figure 1 shows the calculated energy spectrum for the ground states of three-, four-and five-body clusters (lower panel). The corresponding experimental observables are loss peaks in the atom number (upper panel), which appear at values of the scattering length a − , a 4,− , and a 5,− , where the three-, four-and five-body states cross the free atom threshold, respectively. The resonant values of the scattering length are predicted to be universally connected by the simple relations a 4,− = 0.44(1) a − and a 5,− = 0.65(1) a 4,− [37]. A more recent study [38] has theoretically explored the four-boson resonance and has determined its position with greater accuracy. For four-body states, the universal relation has been confirmed in experiments [12,19,20,39] and the four-body recombination rate has been measured [39] and calculated within the hyperspherical framework [34]. A straightforward extension to five-body systems is not currently possible for experiment nor theory. The experimental challenge is to discriminate the five-body recombination signal against a strong background resulting from fewer-body processes. The numerical difficulty of the scattering few-body problem grows exponentially with the number of particles making the description of five and larger systems beyond current theoretical capabilities.
This article presents a combined theoretical and experimental study of universal few-body physics up to five-body states. We present strong evidence for the existence of an Efimov-related cluster state of five identical bosons and we provide quantitative results for the corresponding five-body recombination rate. Our results highlight a new level of understanding concerning few-body physics and its experimental manifestations in ultracold atomic quantum gases.
The "true" Efimov effect refers to the appearance of an infinite number of N -body bound states, which have a discrete scale invariance and which exhibit well-defined thresholds given by the (N − 1)body subsystem.

Theoretical approach
The theoretical analysis of N-body recombination processes requires the description of the N-body scattering continuum. The hyperspherical framework has been successfully applied to describe recombination processes for N > 3 [34,40]. In this framework, the Hamiltonian is diagonalized adiabatically as a function of the hyperradius R, which describes the overall size of the system, leading to a set of coupled one-dimensional Schrödinger equations. At ultralow temperatures and large scattering lengths, N-body recombination events are mainly controlled by scattering processes with incoming flux in the lowest N-body scattering channel and outgoing flux in deeper loss channels. The coupling to the deep channels is assumed to remain approximately unaffected as the scattering length and the collision energy are varied throughout our regime of interest. Thus, all the relevant information comes from the analysis of the lowest N-body potential curve corresponding to the incoming scattering channel.
However, the extraction of the hyperspherical potential curves and couplings becomes computationally unfeasible as the number of particles increases. Currently, the numerical technology allows for the calculation of potential curves of three and fourbody systems but no tractable method has yet been implemented that can calculate both the potential curves and couplings for N > 4. However, the trapped energy eigenvalue spectrum of a five-body system can be accurately obtained with current technology. Thus, the present study extracts the relevant recombination information from an analysis of the trapped spectrum.
Our starting point is the hyperspherical description of the ultracold N-body recombination rate [40], where µ N = m/ N−1 √ N with m being the atomic mass, k = (2µ N E/h 2 ) 1/2 is the incoming hyperradial scattering wavenumber and S 0 + 00 is the diagonal element of the S-matrix for the lowest channel (00) in the J Π = 0 + symmetry. For purely elastic scattering |S 0 + 00 | 2 = 1 and the recombination rate is zero. In the limit in which every N-body collision leads to losses, |S 0 + 00 | 2 = 0. Taking the thermal average in the full loss case at a temperature T , one obtains the unitary limit for N-body recombination at low energy: Based on potential curves computed for the three-and four-body cases, we expect that the five-body potential curve should have the topology depicted as a solid curve in Fig. 2, i.e., the lowest potential curve exhibits a barrier that separates the inner region (small R) from the asymptotic scattering region at large R. Note that the effective mock-centrifugal barrier [41] in the lowest N-body continuum channel at large R is guaranteed to have the form (for finite a), Figure 2. Schematic representation of lowest five-body hyperspherical potential curve. The solid curve represents the free-space hyperspherical potential while dashed curves represent hyperspherical potential curves in the presence of a harmonic potential. Solid lines represent the inner and outer WKB phases. Inset: Five-body trapped energy spectrum in the region where the bound five-body state crosses the lowest trapped states.
this potential curve topology, an extension of the semiclassical (WKBJ) treatment of Berry [42] to include the decay to the lowest channels yields [40]: where φ in is the WKBJ phase for the inner allowed region, γ is the WKBJ tunneling integral in the barrier region, η − describes the decay to deeper, non-universal, channels and is treated as a fitting parameter, and A(η − , γ, φ in ) ensures the proper normalization. Thus, the determination of γ and φ in for the relevant range of energy and scattering length gives an approximate description of five-body recombination. The novelty of our approach is the determination of γ and φ in from an analysis of the trapped spectrum. The five-body trapped energy spectrum is obtained using a correlated Gaussian basis set expansion [37]; for details see Appendix A. In the region where the five-body resonance occurs, the spectrum exhibits a series of avoided crossings between the five-body states bound in the inner region and the outer trap region states (see Fig. 2). In the WKBJ approximation, the quantization condition for the trapped potential curve (dashed curve in Fig.2 For collision energies below the barrier local maximum and away from the avoided crossings, the allowed energy eigenvalues occur when φ α (E, a) ≈ π(i + 1/2), where i is an integer and α = in, out. The phase φ in in the four-boson case near the resonance energy is known from our previous work [40] to be well-described by φ in ≈ φ in,0 + b(a/r vdW ) + caE/(r vdW E vdW ), where φ in,0 , b and c are fitting parameters and r vdW and E vdW are the van der Waals radius and energy, respectively, as defined in [4]. Using this simple form in the five-body case and imposing the eigenstate condition for bound states, the values of φ in,0 , b and c for N = 5 are extracted.
Next, an analysis of the spectrum at the avoided crossings (see e.g. Fig. 2) determines γ, as explained in detail in Appendix B. The relevant avoided crossings occur when φ in ≈ φ out ≈ π(i + 1/2). For narrow avoided crossings (∆ ≪ 1), the quantization condition reduces to δφ in δφ out ≈ ∆ 2 where φ α = π(i + 1/2) + δφ α . Right at the avoided crossing, the energy difference between the two states (2∆E) is related to ∆, namely as δφ α ≈ (dφ α /dE)∆E. The quantization condition thus reduces to ∆ ≈ 2∆E 2 (dφ in /dE)(dφ out /dE). Consequently knowledge of φ in , φ out and the energy avoided crossings allows the tunneling γ to be determined. At large R, interactions can be treated perturbatively leading to an hyperspherical potential curve valid at large R that determines φ out .
By changing the trapping confinement we can change φ out and by introducing short range three-body forces we can modify φ in without affecting the barrier. This allow us to explore how γ depends on both E and a. Interestingly, at low energies our numerical results are nicely fitted by the formula: e −2γ ∝ (E/E vdW ) 5 |a/r vdW | 9.6 , which is in good agreement with the predicted threshold behavior (

Experiment
We prepare an optically trapped sample of cesium atoms in the lowest sub-level of the electronic ground state under similar conditions as described in Ref. [43]. The final evaporation process, performed at a magnetic field of 894 G (a = +285 a 0 , where a 0 is the Bohr radius), is stopped before the onset of Bose-Einstein condensation. The sample is then adiabatically recompressed to avoid further evaporation from the trap. At this point, the sample contains about 6 × 10 4 atoms at a temperature T = 78(3) nK and the confining optical potential has a mean trap frequencyω = 2π × 36.2(2) Hz, which results in a peak number density of 4.2 µm −3 and a peak phase space density close to unity. The wide s-wave Feshbach resonance with its pole at 786 G offers ideal tuning properties [43,44], superior to the low-field region investigated in our previous works [8,39]. We measured the decay of the atom number in the region of negative scattering length (from −500 a 0 at 863 G to −200 a 0 at 873 G), where the four-and five-body recombination resonances are expected. After the recompression stage, we tune the scattering length to its target value, and we measure the atom number after a variable hold time by absorption imaging. After 100 ms, the typical loss fraction is around 10% at −300 a 0 and almost 35% at −450 a 0 .
The time evolution of the number N of trapped atoms and temperature T are determined by the different N-body loss processes and can be expressed [45] in terms of a system of coupled differential equations,  (5) incorporates the antievaporation heating [45] that results from the higher recombination rate in the densest part of the cloud (2ǫ N ≡ 1 − 1/N). Under our experimental conditions, the first two terms of the sums can be omitted. One-body losses, resulting from collisions with the background gas, are negligible on our experimental time scale (L 1 = 0), and the use of atoms in their lowest sub-state assures that two-body losses vanish (L 2 = 0). Therefore, the three-body rate coefficient L 3 is the leading term contributing to the atom losses, and its resonant behavior is related to the Efimov trimers. This contribution is well understood and the coefficient L 3 can be described by the well-established result of effective field theory [46]. The rates of recombination events involving more particles are generally smaller than the one related to three-body losses and the contributions are difficult to discriminate because of the very similar behavior. Since the rate of recombination events for typical gas densities decreases rapidly with N, contributions with N > 5 are considered negligible in the following.
A general fit to the experimental decay curves with L 3 , L 4 , and L 5 as free parameters in Eq. (4) and Eq. (5) turns out to be practically impossible. Therefore, we fix L 3 according to effective field theory, with parameters a − = −955 a 0 and η − = 0.08 as determined for three-body recombination (N = 3) in our previous experiment [43]. We can now interpret the additional losses in terms of four-body and five-body decay. In order to avoid any fitting ambiguities, we chose the simple approach to describe these losses either in terms of an effective four-body loss coefficient L 4,eff (setting L 5 = 0) or an effective five-body loss coefficient L 5,eff (setting L 4 = 0). Figure 3(a) shows L 4,eff as extracted from our experimental data in comparison with the theoretical predictions, obtained by numerically evaluating L 4 from Eqs. (1) and (3) for our experimental conditions. Here we have adjusted the decay parameter to η − = 0.33, which as a non-universal parameter depends on N. The comparison shows that the losses observed around −450 a 0 can be fully attributed to the fourbody recombination resonance. The four-body loss peak position a 4,− = −440(10) a 0 corresponds to 0.46(1) a − and is in very good agreement with the theoretical value 0.44(1) a − [34] and previous observations [12,19,20,39]. In contrast the enhancement of losses centered at about −300 a 0 cannot be interpreted in terms of the known universal four-body cluster states ¶, suggesting that a different loss mechanism is present.
An alternative representation of the same data in terms of L 5,eff is shown in Fig. 3(b), together with the results of our theoretical model for L 5 ; here we adjusted the relevant decay parameter to η − = 0.20. The model nicely explains the loss rates in the region where three-and four-body losses cannot account for the experimental observations. Remarkably, the resonance position a 5,− = 0.64(2) a 4,− is in agreement with the theoretical predictions 0.65(1) a 4,− [36,37]. However, quantitatively, the experimental values for L 5 are about 15 times larger than the calculated ones. To account for this, we introduce a corresponding scaling factor. We find that this deviation might derive from non-universal effects that modify the value of the calculated WKB phase γ by about 10%, which remains within the realistic uncertainty range of our theory.
An experimental search for higher-order recombination resonances (N > 5) that would be expected at lower values of the scattering length did not show clear signatures. A general problem arises, namely that the phase-space density cannot be further increased without causing the collapse of a Bose-Einstein condensate at negative values of a. To induce faster losses, adiabatic compression can increase the density, but then the higher temperatures cause increasing problems with the unitarity limit for high N. By decreasing the temperature, constraints by the unitarity limit can be avoided, but then losses for high N get so small that they become practically unobservable.
Based on the above results and parameters, we model the general loss behavior in the region of interest. Figure 4 shows an example for three, four-and five-body recombination by plotting the atom losses, under typical experimental conditions, for a fixed hold time and variable scattering length. The "family portrait" of N-body recombination highlights the different contributions and confirms how the different loss features dominate the losses at the resonant positions. The experimental data plotted in the inset show that the peak positions and the magnitude of losses are in very good agreement with the simulated losses. Note that the somewhat higher experimental losses can be attributed to an additional loss that occurs during the ramp to the target ¶ We cannot rule out the possibility of a non-universal few-body state, which may be associated with higher partial waves. However, such an accidental coincidence appears to be rather unlikely.  magnetic field strength.

Conclusion
By pushing the limits of ultracold few-body physics, we have explored a universal fivebody recombination resonance both experimentally and theoretically. The observed series of recombination features, which we interpret in terms of three-, four-and fivebody recombination resonances, provides crucial evidence for the existence of a family of universal N-body bound states tied to Efimov trimers. The infinite series of N-boson cluster states represents a paradigm for the general implications of Efimov physics for many-body systems. We speculate that similar scenarios also exist for other few-body systems of increasing size, containing fermionic constituents or particles of different masses, with important consequences for the interaction properties of the many-body system.

Appendix A. Extraction of the trapped five-body spectrum
A crucial ingredient of our recombination analysis is the accurate extraction of the five-body trapped spectrum. The trapped few-boson Hamiltonian is given by where ω is the trapping frequency, V 2b (r) = V 2b0 e −r 2 /(2r 2 0 ) is the two-body potential and V 3b (R) = V 3b0 e −R 2 /(2R 2 0 ) is a three-body potential. Here, r 0 and R 0 are the ranges of the two and three-body potentials that are fixed during the calculation, and R ijk = (r 2 ij + r 2 ik + r 2 jk )/3. The two-body potential strength V 2b0 is used to tune the two-body scattering length and V 3b is set to zero in most calculations (see discussion below). To relate the model potential with the experimental we rely on the universal character of the low energy few-boson physics and we relate r 0 with the van der Waals length r vdw so that the three-and four-body recombination peaks coincide with those experimentally observed. In this transformation, we relate the energy scales so that they scale as inverse length squared.
To extract the recombination parameters, the five-body spectrum is analyzed as a function of the scattering length and the trapping confinement. The scattering length is changed in the region −4r 0 < a < 0 and the trapping frequency is changed so that the corresponding trapping length a ho = h/(mω) varies in the range 5r 0 < a ho ≤ 100r 0 . This region of scattering lengths and energies corresponds to the low energy region where the five-body resonance occurs.
The three-body interaction is used to shift the position of the five-body resonance and explore the dependence of the recombination on the scattering length. The range of the three-body potential is taken to be R 0 , which is smaller than the two-body r 0 so that it is mainly relevant at small hyperradii, i.e. so that its main contribution in our recombination formula is to change the inner phase. For the three-body interactions considered, the spectrum follows the linear dependence with the scattering length that is expected at small and negative scattering lengths. This linear dependence arises from only two-body physics which, in the hyperspherical framework, is described by the long-range behavior of the potential curves. Using the zero-range model of the two-body interaction, one can derive a first order correction of the hyperspherical potential curve. Here, we follow a similar procedure to that on Ref. [48], but for a set of coordinates in which the center of mass has been removed. In this approximation, the lowest potential takes the form The first two terms correspond, respectively, to the hyperangular or "mock-centrifugal" kinetic energy and the trapping potential; and the third term represents the interaction corrections. b Figure A1. (a) The lowest two eigenenergies of a trapped five body system are shown as functions of the scattering length for different trapping frequencies. Different colors represent different trapping frequencies. The combination of these states essentially describes the energy of the five-body state in the inner region of the potential E mol (a) (the diagonal curve). Here E sr =h 2 /(mr 2 0 ) and r 0 is the characteristic range of the two-body model potential that can be tuned to obtain the five-body resonance (i.e. r 0 ∼ 1.7r vdw where r vdw is the van der Waals length). (b) The near-threshold behavior of ∆. The fitting of the lowest energy points leads that ∆ ∝ AE b . The lowest three points lead to b ≈ 5.004 as expected from the known threshold behavior [40].
than the interaction range. Using perturbation analysis, one derives the well-known corrected energy of the lowest state, Our numerical calculations, with and without the three-body forces, show the linear behavior described in Eq. A.3 in the region |a| < |a 5,− |. This suggests that the hyperradial potential in the R ≫ r 0 is well described by Eq. A.2 and that, in this region, the potential is independent of the three-body forces. Thus, it is consistent to interpret the main contribution of the three-body interaction as a modification of the short range physics that controls φ in . The calculations are carried out using a correlated Gaussian basis set expansion limited to describe L P = 0 + trapped states since the energetically lowest scattering continuum corresponds to zero angular momentum and positive parity. In this computation, we are only interested in the lowest trapped states which, in the hyperspherical picture, are supported mainly by the lowest potential curve. Therefore, the basis set is designed to accurately describe those states. The calculations include thousands of basis functions which are optimized for different scattering lengths and trap lengths. To verify the convergence of the energies, we carry out several optimization steps. The typical spectrum obtained by this analysis is shown in Fig. A1(a).
To extract the semiclassical recombination parameters, we repeat the same semiclassical analysis used to derive Eqs. (2,3) but for a hyperspherical potential with the trapping potential (see Fig. 2). This basically amounts to solve a double well problem using a semiclassical analysis. For this analysis, we follow closely the prescription in Ref. [42] to determine the quantization condition where φ in is the phase in the inner region, γ is the barrier phase and φ out is the phase in the external trapped region (see Fig. 2). The semiclassical phases are where q(R) = 2µ[E − V (R)] and V (R) is the hyperspherical curve with the Langer correction. Here, the different integration regions are bounded by the classical turning points, i.e. the R positions at which q(R) = 0. In our analysis, φ and γ are assumed to be unaffected by the trapping potential which is expected to be negligible at small hyperradii. After some mathematical manipulation, the quantization condition can be written as 1 tan(φ in ) tan(φ out ) = ∆ 2 (B.5) where ∆ = e −2γ /2. For tunneling energies well below the barrier height ∆ ≪ 1 which implies that the quantization condition is fulfilled when either φ in ≈ π(i in + 1/2) or φ out ≈ π(i out + 1/2). These conditions (φ in ≈ π(i in + 1/2) and φ out ≈ π(i out + 1/2)) can be interpreted as having a bound state either in the inner or outer region. We call E in (a) and E out (a) the energies at which φ in = π(i in + 1/2) and φ out = π(i out + 1/2), respectively. Here i α are integers representing the number of bound states supported by the inner and outer well respectively. Thus, close to those energies, the inner and outer phases take the form φ α = π(i α + 1/2) + φ ′ α (E − E α (a)) where φ ′ α ≡ dφ α /dE. In the vicinity of the avoided crossings, the quantization condition can be approximated to Equation B.6 resembles a determinant of a 2 × 2 Hamiltonian matrix which lends itself to the following physical interpretation. The states supported by the inner and outer regions are coupled to each other and the coupling energy V in,out is related to the tunneling through the barrier by V 2 in,out = ∆ 2φ ′ in φ ′ out . The eigenenergies of this double well problem (E 1 and E 2 ) reach their squared minimal difference at the avoided crossing, equal to: To verify the validity of this equation, the energy dependence of ∆ is explored by varying the strength of the trapping potential. According to the threshold laws [40], for a five-boson problem, ∆ increases with energy as ∆ ∝ E 5 . Figure A1(b) confirms that the scaling of ∆ calculated through the avoided crossing analysis shows excellent agreement with the threshold law prediction.