Two-body and three-body contacts for three bosons in the unitary regime: Analytic expressions and limiting forms

The two and three-body contacts are central to a set of univeral relations between microscopic few-body physics within an ultracold Bose gas and its thermodynamical properties. They may also be defined in trapped few-particle systems, which is the subject of this work. In this work, we focus on the unitary three-body problem in a trap, where interactions are as strong as allowed by quantum mechanics. We derive analytic results for the two and three-body contacts in this regime and compare with existing limiting expressions and previous numerical studies.


I. INTRODUCTION
As the s-wave scattering length approaches unitarity |a| → ∞, properties of both macroscopic and microscopic strongly-interacting ultracold quantum systems can be described effectively by a reduced set of remaining finite quantities. For three trapped bosonic alkali atoms at unitarity, physics in the unitary regime is parametrized solely by the range of the potential, captured by the van der Waals length r vdW [1], and the trap length a ho = ( /mω) 1/2 , with frequency ω and single-particle mass m. Consequently, these scales parametrize the two and three-body contacts, respectively, that are central to many relevant observables in the system including the tail of the single-particle density, short-distance behavior of correlation functions, high-frequency tail of the rf transition rate, virial theorem, and total energy of the system [2,3]. Intuitively, the two and three-body contacts measure the mean number of clustered pairs and triplets of bosons, respectively [4]. The centrality of the extensive two-body contact C 2 to the thermodynamic properties of the two-component unitary Fermi gas was worked out by Shina Tan [5][6][7], and this quantity frames part of our current understanding of experiments in this regime [8]. For strongly-interacting Bose gases, however, the presence of the Efimov effect [9,10] and associated infinite sequence of bound trimers requires the introduction of an additional parameter, the extensive three-body contact C 3 [2,3]. In the strongly-interacting regime, where interaction length scales are diverging, the contacts may also be used to interpolate in between perturbative results in neighboring regimes [11].
The centrality of the contacts means that they may be measured through various means. By measuring the high-momentum tail of the rf transition rate in a weaklyinteracting Bose gas, Wild et al. [12] measured a nonzero C 2 and found that C 3 was consistent with zero. In the non-equilibrium regime of a nondegenerate Bose gas quenched to unitarity, both C 2 and C 3 were measured interferrometrically in Ref. [13], finding a gradual approach of C 3 to equilibrium. It has also been suggested that the * Corresponding author: colussiv@gmail.com high-momentum tail of the degenerate Bose gas quenched to unitarity measured in Ref. [14] is a nonzero measurement of C 3 [15], although there is currently no consensus and a recent experimental study by Eigen et al. [16] found no evidence of the contacts in the single-particle momentum distribution at early-times after the quench. Theoretically, a prediction for the gradual growth of C 3 in a uniform degenerate Bose gas quenched to unitary was made in Ref. [17] where a coherent beating phenomenon at the frequency of Efimov states was observed.
The contacts can also be measured in microscopic fewbody systems for instance by measuring interaction energies or spectroscopically. Theoretically, the contacts corresponding to the lowest few eigenstates were calculated numerically for three and four trapped bosons in Ref. [18] at fixed trap length. In the present work, we provide analytic results for C 2 and C 3 for the trapped unitary three-body problem including analytic expressions and limiting forms. This includes a review of some of previous results derived already in the supplementary materials of Ref. [17] for C 3 , included here in more detail, in addition to new derivations for C 2 .
In this work, we begin in Sec. II by reviewing analytic solutions of the unitary three-boson problem in a trap in the zero-range model first derived in Ref. [19]. In Sec. III, we use a relation between the short-distance behavior of the three-body correlation function and C 3 from Ref. [17] to derive analytic expressions for the three-body contacts for each trapped three-body eigenstate including limiting expressions in the large-energy limit. In Sec. IV, we perform an analogous derivation except for the two-body contacts for each trapped three-body eigenstate. Finally, in Sec. V, we discuss our findings and compare against the available existing results in [18] which evaluated C 2 and C 3 including more realistic finite-range effects at fixed trap length. We find generally excellent agreement with the available results of that work.

II. THREE TRAPPED BOSONS AT UNITARITY
In this section, we outline analytic solutions of the unitary three-body problem in a trap given in Ref. [19]. Consider three non-interacting bosons located at r 1 , r 2 , r 3 in a harmonic trap whose wave function Ψ(r 1 , r 2 , r 3 ) solves the non-interacting Schrdinger equation The three-body wave function can be separated |Ψ = |Ψ CM |ψ into center of mass |Ψ CM and relative |ψ components. The center of mass coordinate is defined as C ≡ (r 1 + r 2 + r 3 )/3, and the relative motion is parametrized by the hyperradius R 2 ≡ (r 2 + ρ 2 )/2 and hyperangles Ω ≡ {ρ,r, α} where r ≡ r 2 − r 1 and ρ = (2r 3 − r 1 − r 2 )/ √ 3 are the Jacobi vectors with spherical anglesρ andr, respectively. The hyperangle α is defined as tan α = r/ρ, and is related to each individual Jacobi vector as In ultracold quantum gases of alkali atoms, it is generally the case that the typical momentum of the problem are such that kr vdW 1, where r vdW is the van der Waals length which furnishes a natural short-range for pairwise interactions [1]. In this limit, the zero-range approximation can be made, and pairwise interactions can be accounted for via the Bethe-Peierls boundary conditions where the unknown function A is determined by solving the three-body Schrdinger equation (Eq. (1)) with the above boundary condition.
In the unitary regime, |a| → ∞ in Eq. (2) and the relative three-body wave function becomes singular as 1/r in the r → 0 limit. Following Ref. [19] and the original approach of Efimov [9], the relative three-body wave function in this regime can be decomposed as ψ(R, Ω) = R −2 N F (R)φ(Ω) in terms of hyperradial F (R) and hyperrangular φ(Ω)) = (1 +Q)ϕ(α)/ sin 2α √ 4π functions with normalization constant N . The operatorQ = P 13 +P 23 is written in terms of the permuation oper-atorP ij that swaps particles i and j. The hyperangular basis functions ϕ are determined from solutions of the equations [19] −ϕ s (α) = s 2 ϕ s (α), which are solved by ϕ s (α) = sin(s(π/2−α)). The channel index s is obtained by solving the transcendental equation Each channel supports an associated set of hyperradial eigenfunctions F (s) with quantum number j = 0, 1, 2, . . . and relative threebody eigenenergy E where L (s) j is a generalized Laguerre polynomial of degree j, and the associated eigenenergies are E (s,j) 3b = ω(s + 2j + 1). We label the channels with s 2 > 0 as universal because they depend solely on the trapping parameters and not on microscopic details of the interaction.
The lone imaginary root of Eq. (6) s = is 0 denotes the Efimov channel, where s 0 ≈ 1.00624 is Efimov's constant [9]. Equation (7) must be supplemented with the additional boundary condition F (s0) j (R) ∝ sin[s 0 ln(R/R t )] to preserve the Hermiticity of the problem [19,20], where R t is a three-body parameter setting the phase of the log-periodic oscillation. The hyperradial solutions in this channel are where W is a Whittaker function [21]. The eigenenergy spectrum in the Efimov channel is obtained by solving which is understood mod π and was first obtained in Ref. [22]. To make the connection between the threebody parameters R t and κ * , we take the j → −∞ limit of Eq. 10, using the relations [21] arg where ψ in the above equation is the digamma function. Inserting this limiting result into Eq. (10) gives the freespace result for the binding energies of Efimov trimers Setting j = 0 in the above equation matches the ground-state trimer energy at unitarity for neutral atoms studied in Ref. [23]. In that work, the universality of κ * with the van der Waals length was established, finding κ * = 0.226/r vdW . Therefore, we make the arbitrary distinction that j = 0 is the "ground-state" Efimov trimer in the sense of Ref. [23] and j < 0 comprises the remainder of the Thomas collapse [24]. Taking the j → −∞ limit of the Whittaker function in Eq. 9 gives where K is0 is a modified Bessel function of imaginary argument, which is the free-space bound-state wave func-tion for an Efimov trimer [9]. The Efimov channel depends on the trapping parameters and also on R t , which is sensitive to the microscopic length scale r vdW . The Efimov channel is therefore not universal. Furthermore, although R t depends on r vdW , the model remains zerorange involving only contact potentials. We contrast this with the approach of Ref. [18], which used more realistic finite-ranged potentials. For completeness, we discuss also the j → +∞ limit of Eq. 10, where Stirling's formula gives Therefore, Eq. 10 in this limit becomes which is in general difficult to solve. However, Eq. 15 is written in a form "f (x) = a + g(f (x))" that can be iterated. In the E (s0,j) 3b → ∞ limit this produces corrections in powers of j as which was first obtained in Ref. [25]. The normalization constants N s,j are defined as where the hyperangular inner product is defined in terms of the hyperangular solid angle as The result in Eq. 19, was first calculated in Ref. [25] using a change of variables trick due to Efimov in Ref. [26]. The hyperradial inner product is defined as where the relevant integrals can be found tabulated in Ref. [27]. The factor

III. THREE-BODY CONTACT C3
Any one of the universal relations involving the contacts can be used as a starting point, however, we choose to obtain them through the limiting behavior of few-body correlation functions [2,3]. We begin by deriving the extensive three-body contacts C 3 for each trapped threebody eigenstate.
The non-normalized triplet correlation function is defined in second-quantization as [28] in terms of the bosonic field operatorsψ(r) = kâ k e ik·r / √ V and system volume V . For three bosons in vacuum, this is equivalent in first quantization to Following [2], we integrate out the center of mass dependence, and take the limit R → 0 at fixed Ω for the j th eigenstate in the Efimov channel. In the above equation, we have indexed the extensive three-body contact specific to the eigenstate (s, j) as C (s,j) . In Eq. (27), ψ sc (R, Ω) is zero-energy three-body scattering state [2] ψ sc (R, Ω) = 1 R 2 sin s 0 ln Integrating both sides of Eq. (27) over dΩ, and taking advantage of the orthogonality of the hyperangular eigenfunctions so that only s = is 0 need be considered yields From the asymptotic behavior of the Whittaker functions [21], the R → 0 limit in the above equation can be taken with the simple result In the above equation, C 3 is made dimensionless by rescaling in powers of the oscillator length. In the following two subsections, we derive first the |E 3b |/ ω → ∞ limiting forms of Eq. (30) and then show some results at intermediate energies.
Case I (E 3b / ω → −∞): In this limit, the digamma function has the asymptotic form [21] Using the property ln z = ln |z| + i arg z, we find We therefore obtain the result for where we have used the result of Eq. (12). The above result matches the operational definition /m in the free-space limit [2,3].

B. Results
Case II (E 3b / ω → +∞): In this limit, the asymptotic expansion of the digamma function (Eq. (31)) should be used carefully because we are taking the limit | arg z| → π. Instead we use the reflection formula for the digamma function [21] to obtain From the discussion of Case I, we know that the final term on the right in the above equation will vanish in the limit E 3b / ω → +∞. Therefore, if the term involving the imaginary part of the cotangent is nonvanishing, it is the only term that needs to be considered. Following this reasoning and using properties of cotangent for imaginary argument [21], we find that which is nonvanishing. We now take the limit (E 3b / ω → +∞) and use Eq. 16 to find a 2 ho C (s0,j) 3 ≈ s 0 π cosh(πs 0 ) + cos(s 0 ln(jR 2 t /a 2 ho ) − 2 arg Γ(1 + is 0 )) sinh(πs 0 ) , ≈ s 0 π 1 + 0.085 cos(s 0 ln(jR 2 t /a 2 ho ) − 2 arg Γ(1 + is 0 )) .
The leading order result a 2 ho C (s0,j) 3 ≈ s 0 /π was obtained previously in Ref. [25] for an analogous expression for the decay rate Γ in this limit. Having outlined the limiting behavior of C (s0,j) 3 in Sec. III A, we now discuss its behavior in general and over a range of a ho κ * as shown in Fig. 1. We find that C (s0,j) 3 matches the limiting forms (Eqs. (33), (36)) by |E (s0,j) 3b |/ ω ∼ 10 2 as shown in Fig. 1(a). At intermediate energies, as a ho decreases Efimov trimers whose free-space energies are comparable to ω are shifted above the zero-point energy of the trap as studied in Ref. [29]. This process repeats logperiodically as a ho κ * → 0. A single log-period is shown in Fig. 1(a), where the shift of the j = 2 Efimov trimer is clearly visible. For |E 3b |/ ω < 1, we find that a 2 ho C 3 as a function of E 3b / ω interpolates smoothly between the limiting behaviors as shown in Fig. 1(b).
In general we find that C (s0,j) 3 for a particular Efimov trimer is increased compared to the free space limit [see Eq. (33)]. We echo the conclusion of Ref. [18] that the restriction of a trimer to a reduced region of space should correspond to an increase in the probability to find three bosons at short distances as encoded in C where the phase factor Λ ≈ 0.697. The factor 2s 0 ln(κ * a ho ) = 2s 0 ln(κ * a ho e π/s0 ) mod 2π, which gives the log-periodic behavior seen in Fig. 1.

IV. TWO-BODY CONTACT C2
In this section, we derive the extensive two-body contacts C 2 for each trapped three-body eigenstate. Analogous to the treatment in Sec. III, we begin from the definition of the non-normalized pair correlation function [28] For three bosons in vaccum, this is equivalent in first quantization to We then integrate over the center of mass coordinate c ≡ (r 1 + r 2 )/2, and do a change of variables {c, r 3 } → {C, r 3,12 } where r 3,12 ≡ r 3 − c. The Jacobian for this transformation is unity, and we obtain d 3 c G (2) (r 1 , r 2 ) = 3! d 3 r 3,12 |ψ(r, r 3,12 )| 2 . (40) In the limit |r| → 0, the left hand side of the above equation is related to C 2 [2], and therefore we find the following relation between the relative three-body wave function and the extensive two-body contact 3! d 3 r 3,12 |ψ(r, r 3,12 ) In order to take the above limit, we rewrite the threebody relative wave function as Ψ(r, r 3,12 ) = χ 0 (r, ρ) rρ + χ 0 (r (2) , ρ (2) ) r (2) where we have implicitly defined r ≡ r (1) , ρ ≡ ρ (1) . With some rearrangement and the help of the kinematic rotations, we take the limit in Eq. (41) to obtain Written in terms of the trapped eigenstates outlined in Sec. II, this becomes where we have indexed the extensive two-body contact specific to the eigenstate (s, j) as C (s,j) 2 . The notation A s,s j,j is shorthand for diagonal components of the array It is possible to obtain analytic expressions for the array A s,s j,j using the trapped hyperradial eigenstates outlined in Sec. II. We are concerned however only with the components with j = j and s = s required to evaluate C (s,j) 2 via Eq. (45). For s = is 0 , Eq. 46 is an integral over Whittaker functions that can be found in a table [27] with result where 3 F 2 is a generalized hypergeometric function. We note that an equivalent expression to Eq. 47 was first obtain in Ref. [25] for the loss-rate Γ. For channels with s 2 > 0, Eq. 46 is an integral over generalized Laguerre polynomials that can be recast using the recurrence relation [21] L (s) where the (·) is the generalized binomial coefficient. Using this relation, the array A s,s j,j can be found tabulated in Ref. [27], and we find that In the following subsection, we derive some limiting cases of A s,s j,j , and then discuss results for C (s,j) 2 in general.

A. Limiting cases
Due to the 3 F 2 function in Eqs. (47), (49) it is unclear how to derive limiting expressions for C (s,j) 2 in the limit j → +∞ [31]. It is however possible to derive expressions in the opposite limits, which includes the j → −∞ limit in the Efimov channel and the case j = 0 for the universal channels.
For s 2 > 0 with j = 0, the array A s,s 0,0 has the form and we obtain the following expression which vanishes in the limit s → +∞.

B. Results
Having outlined some of the limiting behaviors of C (s,j) 2 in Sec. IV A, we now analyze its behavior in general over a range of universal channels and for various values of a ho κ * in the Efimov channel as shown in Fig. 2. First we discuss results for the universal channels shown in Fig. 2(a-b). In Fig. 2(a), we see that our asymptotic prediction for s → +∞ with j = 0 given by Eq. (54) matches our calculation of a ho C 2 as s n 1 (solid black line). At larger j, the behavior of a ho C 2 is shown in Fig. 2(b), where we see a gradual decay. The j → +∞ asymptotic form of a ho C 2 was not derived in Sec. IV A, although we estimate its form by fitting a ho C 2 to a power law j −α . We find that for extremely large energies E 3b / ω > 10 6 this exponent approaches α ≈ 0.5 in both universal and Efimov channels. This power law can be motivated by approximating 1/2 − j − s ≈ −j − s in the arguement of the 3 F 2 function in Eq. 49. The 3 F 2 function then reduces to 2 F 1 [31] and the j → +∞ limit can be taken to obtain a j −1/2 power law scaling for a ho C 2 .
Our results for a ho C 2 in the Efimov channel are shown in Fig. 2(c-d) for both positive and negative energies starting at j = 0. We see that the E 3b / ω → −∞ asymptotic form given in Eq. (51) (dot-dashed black line) matches our calculation of a ho C 2 as shown in Fig. 2(c). At positive energies, a gradual decay of a ho C 2 is evident as in the universal case [ Fig. 2(c)], and the log-periodic oscillations that were prominent in C 3 in Fig. 1 are absent. For varying a ho κ * , our results follow the same trends. As j → +∞ we observe numerically the same j −1/2 power discussed for the universal channels. Additionally, in 2(d) we show the behavior of a ho C s0,j 2 as a function of E 3b / ω in the region |E 3b |/ ω < 1. Here we find, as in Fig. 1(b), a smooth interpolation between limiting behaviors.
In general, we find that over the same energy ranges a ho C 2 is larger in the Efimov channel than in the universal channels. Physically, we understand this as a larger probability for finding correlated pairs of atoms at shortdistances, which is measured by C 2 [4]. This should be larger in the absence of a repulsive hyperradial barrier for s 2 > 0 (see Eq. (7).) As s increases, this barrier becomes more repulsive [see Eq. (7)], and it is therefore less likely to find correlated pairs at short-distances. This is reflected in Fig. 2(a-b). When comparing to the freespace result in Eq. (51) for an Efimov trimer, we find that the trap leads to a relative increase of C 2 . Tighter traps increasingly confine the atoms, and therefore it makes sense that the probability to find correlated pairs should increase accordingly. In this work, analytic solutions of the trapped unitary three-boson problem from Ref. [19] are utilized to derive analytic expressions for the two and three-body contacts for each trapped three-body eigenstate. These results can be used along with the set of universal relations associated with the contacts [2,3] to predict, for instance, the high-frequency tail of the rf transition rate, the virial theorem, loss-rates, and total energy. We derive the large energy asymptotics of the contacts, and compare with existing results in the literature when possible finding good agreement. At intermediate energies, the only results in the literature that we are aware of are from Ref. [18], which were obtained using more realistic Gaussian potentials with a finite-range. We take a wave function based approach to obtain the contacts from the short-distance behavior of the pair and triplet correlation functions. In Ref. [18], however, the contacts were obtained from their operational definitions [2,3] a∂ a E (s,j) 3b by taking the derivatives of three-body eigenenergies with variations of a and κ * . This amounted to changing the underlying parameters of the finite-range Gaussian potentials and approximating the derivatives via finite differencing. In Table I, we compare against results in that work, finding good agreement for the eigenenergies and contacts to within a few percent or less. In Ref. [18], the convergence of eigenenergies was estimated as 10 −3 E ho or better. The eigenenergies calculated in the present work can be calculated from Eq. (10) to arbitrary precision. We attribute the disagreement therefore to the absence of finite-range effects in the zero-range model in the present work.
Finally, we note that the method outlined in these notes to obtain C 2 and C 3 through the limiting behavior of the correlation functions applies also away from unitarity at finite scattering lengths provided the zerorange approximation can be made [33]. Additionally, the derivation in of C 2 from a three-body model presented in Sec. IV may be generalized in the spirit of Ref. [17] to study the dynamics of the two-body contact following an interaction quench. We address this latter problem elsewhere [34] and leave the former as the subject of future study.