Three-body contact for fermions. I. General relations

We consider the resonant Fermi gas, that is, two-component fermions in three dimensions interacting by a short-range potential of large scattering length. We introduce a quantity, the three-body contact, that determines several observables. Within the zero-range model, the number of nearby fermion triplets, the large-momentum tail of the center-of-mass momentum distribution of nearby fermion pairs, as well as the large-momentum tail of the two-particle momentum distribution, are expressed in terms of the three-body contact. For a small finite interaction range, the formation rate of deeply bound dimers by three-body recombination, as well as the three-body contribution to the finite-range correction to the energy, are expressed in terms of the three-body contact and of a three-body parameter. This three-body parameter, which vanishes in the zero-range limit, is defined through the asymptotic behavior of the zero-energy scattering state at distances intermediate between the range and the two-body scattering length. In general, the three-body contact has different contributions labeled by spin and angular momentum indices, and the three-body parameter can depend on those indices. We also include the generalization to unequal masses for $\uparrow$ and $\downarrow$ particles. With respect to the relation between three-body loss rate and number of nearby triplets stated in [Petrov, Salomon and Shlyapnikov, PRL 93, 090404 (2004)], the present work adds a derivation, expresses the proportionality factor in terms of the three-body parameter, and includes the general case where there are several contributions to the three-body contact and several three-body parameters.


Introduction
Over the last twenty years, the two-component Fermi gas with zero-range interactions in three dimensions has become one of the most extensively studied quantum many-body problems.One considers particles with two internal states (denoted ↑ and ↓) and an interaction of vanishing range characterized by its s-wave scattering length a 2 .When 1/a 2 changes from −∞ to +∞, the interaction changes from weakly to strongly attractive, leading to the BCS to BEC crossover.The strongly correlated regime is reached in the central region of the crossover, around the unitary limit 1/a 2 = 0.While the model was historically introduced as a theoretical abstraction [1][2][3], it accurately describes ultracold gases of fermionic atoms in two hyperfine states near a Feshbach resonance, which are the subject of numerous experimental studies, see e.g. .
For bosons in 3D with resonant interactions, in addition to the two-body contact C 2 , measured in [81,82], a three-body contact C 3 appears in several exact relations [46,52,72,83,84].This appearance of C 3 is linked to the Efimov effect.In particular, a three-body parameter has to be included in the definition of the zero-range model, and C 3 is proportional to the derivative of the energy w.r.t. the three-body parameter.C 3 was measured interferometrically in [82] after an interaction quench to unitarity.A three-body contact also appears for single-component fermions with higher-partial-wave resonant interactions, both in 2D [85] (where a super-Efimov effect occurs) and in 1D [86,87].Two-body and three-body contacts were also found to be useful to describe short-distances or large-momenta properties in clusters of helium atoms [88] and in nuclei [89][90][91][92][93][94][95], although the corresponding relations are only approximate because the interaction range is not much smaller than the interparticle distance.
Here we show that a three-body contact C 3 plays an important role for two-component fermions with resonant interactions in 3D, although there is no three-body Efimov effect so that the zero-range model is parameterized by the scattering length a 2 without any three-body parameter [96,97].We work within the zero-range model in Section 2, and we consider models with a small finite interaction range in Section 3. Within the zero-range model, the number of triplets of particles separated by a small distance (Section 2.1), the tail of the center-of-mass momentum distribution of pairs separated by a small distance (Sec.2.4), and the tail of the two-particle momentum distribution (Sec.2.5) are expressed in terms of C 3 , and C 3 is also related to the third order density correlation function (Sec.2.2) and to the behavior of the manybody wavefunction when three particles approach each other (Sec.2.3).When the interaction range b is non-zero but still small compared to the other typical lengthscales, we consider two additional observables, the formation rate of deeply bound dimers by three-body recombination Γ 3 (Sec.3.1), and the three-body contribution to the energy correction induced by the finite interaction range δE 3 (Sec.3.2).We express Γ 3 and δE 3 in terms of C 3 , and of a three-body parameter a 3 (which is small in the zero-range regime).We define a 3 through the asymptotic behavior of the three-body zero-energy scattering state at distances ≫ b and ≪ |a 2 |.
We consider N ↑ fermions of spin ↑ and N ↓ fermions of spin ↓, either confined by a smooth external trapping potential, or in a box with periodic boundary conditions.We consider equal masses for ↑ and ↓ particles, and discuss the unequal-mass case in Appendix D. We consider a stationary state throughout the article, and discuss statistical mixtures and non-stationary states in Appendix E. In parallel with presenting the relations involving the three-body contact, we will recall for comparison the known Tan relations involving the two-body contact.

Relations for the zero-range model
In this Section we work within the zero-range model, where interactions are characterized by a single parameter, the two-body scattering length a 2 .The zero-range model is defined in Eqs.(13,14).The zero-range limit of finite-range models is expected to be universally described by the zero-range model. 1

Number of nearby fermion triplets
If one measures the positions of all particles, the average number of pairs of particles whose separation is smaller than some ϵ is given by where C 2 is the two-body contact [35,36].Similarly, let us consider the number of triplets of fermions separated by small distances.For three particles 1,2,3, let us define the hyperradius where r i j is the distance between particles i and j .If one measures the positions of all particles, the average number of triplets of particles with hyperradius R < ϵ is given by where the prefactor C 3 is what we call the three-body contact, while the exponent s = 1.772724267 . . . is the lowest positive solution different from 1 of The scaling N 3 (ϵ) ∝ ϵ 2s+2 was already obtained in [98,115] (see Section 2.3 for a rederivation).The anomalous exponent 2s + 2 comes from the analytical solution of the unitary three-body problem [96], and is directly linked to a hidden dynamical symmetry and a separability of the three-body problem in hyperspherical coordinates [38,116,117], or in a field theory point of view, to non-relativistic conformal invariance, with s + 5/2 the scaling dimension of a three-fermion operator [118][119][120].
In general there are two different contributions to the three-body contact, coming from ↑↑↓ and ↑↓↓ spin configurations: Denoting by N 2,1 (resp.N 1,2 ) the contributions to N 3 (ϵ) from triplets of particles of spins ↑↑↓ (resp.↑↓↓), we have where C 2,1 (resp.C 1, 2) is what we call the ↑↑↓ (resp.↑↓↓) three-body contact.Clearly, Remarks: • Due to the antibunching effect associated to the Pauli exclusion between fermions with identical spins, the contribution to Eq. ( 1) coming from pairs of particles with identical spins (↑↑ or ↓↓) is negligible in the ϵ → 0 limit (it scales as ϵ 5 ), and N 2 (ϵ) is dominated by the contribution from pairs of particles with opposite spins (↑↓).
Similarly, the contribution to Eq. ( 3) coming from triplets of particles with identical spins (↑↑↑ or ↓↓↓) is negligible in the ϵ → 0 limit (it scales as ϵ 10 ), and N 3 (ϵ) is dominated by the contributions from triplets of particles with non-identical spins, ↑↑↓ or ↑↓↓, in agreement with Eq. ( 7).• For comparison, in the non-interacting case, the number of nearby pairs and triplets scales as (more generally, these scalings also hold with a finite-range interaction that does not diverge too strongly at small distance, so that the wavefunction is bounded).Equation (8) [resp.Equation ( 9)] is dominated by the contribution from pairs (resp.triplets) of particles with non-identical spins.Equation ( 9) includes the antibunching effect due to the Pauli exclusion between the two identical-spin fermions: The wavefunction vanishes linearly with the distance between these fermions, hence an ϵ 2 suppression factor compared to the completely uncorrelated case of non-interacting distinguishable particles • Compared to the non-interacting case Eqs.(8,9), the exponents in Eqs.(1,3) are reduced, i.e. the probability to find particles near to each other is enhanced.This bunching effect is due to the attractive effect of the resonant zero-range interaction, which causes the wavefunction to diverge: When the distance r between two opposite-spin particles vanishes, ψ ∝ 1/r , which yields Eq. (1), while in the limit of vanishing hyperradius R between three particles, ψ ∝ R s−2 , which yields Eq. ( 3).Note that since s < 2, the wavefunction indeed diverges for R → 0, and the exponent in Eq. ( 3) is smaller than in the uncorrelated case Eq. (10), which means that the bunching effect due to the zerorange interactions overcompensates the antibunching effect due to Pauli exclusion. 2

Link with the many-body wavefunction
When three particles approach each other, the many-body wavefunction has a singular asymptotic behavior, with a prefactor related to three-body contact.Let ψ(r 1 , . . ., r N ) be the (orbital) many-body wavefunction.Without loss of generality, we can assume N ↑ ≥ 2, and consider that particles (1, 2, 3) have spins (↑, ↑, ↓).We take this convention throughout the article.This means that ψ is antisymmetric w.r.t.exchange of r 1 and r 2 (ψ is also antisymmetric w.r.t.exchange of any other pair of same-spin particles).We take the normalization d3 r 1 . . .
Within the zero-range model, the stationary Schrödinger equation contains an external trapping potential U but no interaction potential; instead, ψ should satisfy a contact condition in the limit where two particles of opposite spin approach each other: There exists A such that where a 2 is the two-body scattering length, r = ∥r 1 − r 3 ∥ is the distance between the oppositespin particles 1 and 3, and c = (r 1 + r 3 )/2 is their center-of-mass.The limit r → 0 is taken for fixed c and fixed positions of the remaining particles (r 2 , r 4 , . . ., r N ).By antisymmetry, a similar contact condition automatically also holds for all other pairs of opposite-spin particles, and Eqs.(13,14) are sufficient to define the eigenstates ψ and energies E of the zero-range model. 3 When particles 1, 2 and 3 approach each other, the wavefunction of any stationary state has the asymptotic behavior [38,96,98,115] Here, as in Eq. ( 11), R is the hyperradius of particles (1, 2, 3) defined in Eq. ( 2), C is their centerof-mass, and Ω denotes their hyperangles.The limit R → 0 is taken for fixed Ω, C, r 4 , . . ., r N .The unitary hyperangular wavefunctions φ m (Ω) are such that R s−2 φ m (Ω) is a solution of the three-body problem at zero energy and infinite scattering length with total angular momentum quantum numbers l = 1 and m ∈ {−1, 0, 1}, see Appendix A for more details. 4,5  It is known [35] that C 2 is given by the norm of the function A that appears in Eq. ( 14), Similarly, the three-body contact is given by the norm of the function B that appears in Eq. ( 15), The expression (5) for the number of nearby ↑↑↓ fermion triplets, together with Eq. ( 17), simply follow from Eq. ( 15) by integrating |ψ| 2 over the R < ϵ region. 6Similarly, the relation involving g 2,1 , Eq. ( 11), follows immediately from Eqs. (15,17) and from the expression of g 2,1 in first quantization, g 2,1 (r 1 , r 2 , There is a completely analogous relation between C 1,2 and the behavior of the many-body wavefunction when three particles of spins ↑↓↓ approach each other (provided N ↓ ≥ 2).Specifically, considering that particle 4 has spin ↓ (while particles 1, 3 still have spins ↑, ↓) and denoting by R, Ω, C the hyperradius, hyperangles and center-of-mass associated to particles 3, 4, 1 [obtained by replacing (r 1 , r 2 , r 3 ) with (r 3 , r 4 , r 1 ) in Eqs.(88,89,90)] we have 3 Configurations with a vanishing interparticle distance are implicitly excluded in (13).In an equivalent alternative formulation, these configurations are included and regularized delta pseudopotential terms are added [97,116]. 4There is a similarity between the two-body and three-body short-distance asymptotic behaviors Eqs. ( 14) and (15), given that 1/r − 1/a 2 is a solution of the two-body problem at zero energy. 5An asymptotic behavior similar to (15) holds when any three particles with spins ↑↑↓ approach each other, the functions B m corresponding to different triplets of particles being simply related to each other by antisymmtery. 6Here we used the change of integration variables (r 1 , r 2 , r 3 ) −→ (r, ρ, C) with ρ defined in (88), of Jacobian 3 .We also used the property which yields the relations for N 1,2 and g 1,2 , Eqs. (6,12), with Here we assumed N ↓ ≥ 2; if N ↓ = 1, then N 1,2 is obviously zero, and C 1,2 = 0.
Remark: Higher-body contacts can be defined in the same way than the three-body contacts.
When j ↑ particles of spin ↑ and j ↓ particles of spin ↓ approach each other, the N -body wavefunction factorizes into the product of (i) a function of the relative positions of the j = j ↑ + j ↓ nearby particles, given by a zero-energy solution of the j ↑ + j ↓ body problem at a 2 =∞, proportional to the j body hyperradius to some power, and (ii) a function of the center-of-mass of the j nearby particles and of the positions of the (N − j ) other particles [38,115].The L 2 norm of the latter function defines the j ↑ + j ↓ body contact C j ↑ , j ↓ (up to a prefactor which is a matter of definition).

Large-momentum tail of the center-of-mass momentum distribution of nearby pairs
Since C 2 and C 3 determine short-distance singularities, it is natural that they also determine large-momentum tails.C 2 determines the leading tail of the single-particle momentum distribution [36] N with the normalization N σ (k) d 3 k/(2π) 3 = N σ (in the case of periodic boundary conditions, momenta become discrete and momentum integrals should be replaced by sums).C 3 also determines a large-momentum tail.Suppose that one measures, for a pair of particles with opposite spin, both their spatial separation r and their center-of-mass momentum K (this is allowed since the corresponding operators commute).Let N 2 (ϵ, K) be the probability distribution over K conditional to r < ϵ, with the normalization N 2 (ϵ, K) d 3 K /(2π) 3 = N 2 (ϵ).In other words, N 2 (ϵ, K) is the center-of-mass momentum distribution of the pairs of particles separated by a distance < ϵ, normalized to the total number of such pairs. 7We have where N P (K) is what we call the center-of-mass momentum distribution of nearby fermion pairs.With this definition, one simply has the normalization as a consequence of (1). 8The tail of N P (K) is determined by the three-body contact 9 : 7 More formally, N 2 (ϵ, K) is the expectation value of the operator i :↑, j :↓ θ(ϵ − ri j ) (2π) 3 δ 3 ( Ki j − K), where θ is the Heaviside function, while ri j = ∥r j − ri ∥ and Ki j = ki + kj are the operators corresponding to the relative distance and the center-of-mass momentum of particles i and j . 8N P (K) appears naturally in the diagrammatic formalism: In a homogeneous system, N P (K) divided by the volume equals −(m 2 /ħ 2 ) Γ(K, τ = 0 − ) where Γ is the pair-propagator defined e.g. in [66].This can be shown using a lattice model [56], for which 1/(r r ′ ) in ( 27) can be replaced by φ(r)φ(r ′ ) with φ the zero-energy two-body scattering state, and setting r = r ′ = 0. 9 The fact that Γ(K , τ = 0 − ) has a tail ∝ C 3 /K 2s+4 was pointed out to us by Shina Tan (private communication, Aspen, 2011).
with the prefactor whose numerical value is Here NP (K ) stands for the angular average N P (K) d K/(4π) (with K := K/K , and d K the differential solid angle, so that To derive this result, we consider the two-body reduced density matrix Inserting the two-body contact condition ( 14) into (26) yields where We see that g P can be physically interpreted as a coherence function for pairs of nearby fermions.Accordingly, N P is related to g P by Fourier transformation, as can be formally shown using the definition (21) of N P .Hence We then follow a reasoning resembling the one used to derive Eq. (20) in Sec.IV.A of [56].
In the large K limit, the Fourier transform with respect to c in (30) is dominated by the contributions from the singularities of A(c; r 2 , r 4 , . . ., r N ), which occur when c approaches one of the r i (i = 2, 4, . . ., N ).This corresponds to particles 1, 3 and i being close to each other [since A(c; r 2 , . ..) determines the wavefunction ψ when particles 1 and 3 are close to c, according to Eq. ( 14)].For example, for i = 2, the behavior of A(c; r 2 , r 4 , . . ., r N ) in the limit c → r 2 is determined by the asymptotic behavior of ψ when the three particles 1, 2, 3 are close, given by Eq. (15).Therefore, we just need to take the limit where particles 1 and 3 approach each other in (15) to obtain where we used the expression (98,99) of φ m (Ω).There is a similar singularity of A(c; r 2 , r 4 , . . ., r N ) when c approaches r i with 4 ≤ i ≤ N ; when particle i has spin ↓, the function Bm introduced in (18) appears instead of B m .This gives (1 − 2 δ i ,4 ) e −i K•r i Bm (P 4i (r 2 ; r 4 , . . ., r N )) (31) where the first (resp.second) sum over i is taken over particles with spin ↑ (resp.↓), P j i (r 2 ; r 4 , . . ., r N ) is obtained from (r 2 ; r 4 , . . ., r N ) by exchanging r j with r i (P i i is the identity), and I m (K) := s d 3 u e i K•u u s−1 Y m 1 ( û).The latter integral can be evaluated analytically: Using with j 1 (t ) = (sin t −t cos t )/t 2 , and evaluating the remaining integral over u by integrating along a closed contour including the positive real axis and negative imaginary axis, we get Inserting (31,32) into (30), expanding the modulus squared, and neglecting in the large K limit the cross terms coming from two different values of i , 10 we obtain the result (23,24,25), where we used the value of N given in Eq. ( 100) of Appendix A.

Large-momentum tail of the two-particle momentum distribution
The three-body contact also determines the asymptotic behavior at large momenta of the twoparticle opposite-spin momentum distribution function, defined by where Nσ (k) := ĉ † σ (k) ĉσ (k) with ĉσ (k) = d 3 r e −i k•r ψσ (r) the annihilation operator of a particle of spin σ in the state |k〉 defined by 〈r|k〉 = e i k•r .Since Nσ (k) Experimentally, the two-particle momentum distribution can be accessed from the statistics of time-of-flight images, as was recently demonstrated in 2D [122]. 11 Taking the limit of a large relative momentum k → ∞, and integrating over the center-of-mass momentum K, one obtains a tail proportional to C 2 , as pointed out in [128].If we instead send K to infinity, and average over the direction of K, we obtain a tail proportional to C 3 , lim where the constant M was given in (24).This result immediately follows from (23) and from the relation obtained in [128] and rederived in the sequel. 12,13The definition (33) yields 10 By power counting, these cross terms give rise to a ∝ 1/K 2s 4 +4 tail of NP (K ), where s 4 is the smallest scaling exponent of the unitary four-fermion problem; this is indeed negligible compared to the leading 1/K 2s+4 tail of NP (K ), given that s 4 = 2.509(1) is larger than s ≡ s 3 .This value of s 4 follows from the four-body ground-state energy E 4 in an isotropic harmonic trap computed in [121] and the relation [115,117] E N = (s N + 5/2)ħω. 11In 3D, early measurements in the BEC regime were reported in [123]; see also [124] for a recent numerical study.In optical lattices, detailed experimental studies were carried out in recent years using metastable helium atoms [125][126][127].
12 Relation (34) also follows from (36), given (22). 13In Ref. [128], N P (K) was defined by where ρ 2 is the two-body reduced density matrix, which has the diverging behavior (27) when r and r ′ tend to zero.This short-distance divergence leads to a k → ∞ tail of the Fourier transform Eq. ( 37), which can be computed by replacing ρ 2 with its asymptotic expression (27), and using the identity (in the sense of distributions) d 3 r e −i k•r /r = 4π/k 2 .Using (29) then yields (36).

Relations for finite-range models
In this Section, we go beyond the zero-range model, and consider interactions of small but nonzero range.We express two observable in terms of the three-body contact: the rate of three-body recombinations towards deeply bound dimers (Sec.3.1), and the three-body contribution to the energy difference between finite-range and zero-range models (Sec.3.2).

Three-body loss rate
In ultracold atom experiments, three-body losses generically take place, being a manifestation of the fact that the true equilibrium state at such low temperatures is not gaseous (with the exception of polarized hydrogen).In this Section we relate the rate of three-body losses to the three-body contact.

Simple finite-range interaction
To describe three-body losses, we need to go beyond the zero-range model.In this subsection we consider a simple model where fermions of different spin interact through a rotationally invariant potential V 2 (r ), of finite range b. 14 We consider the resonant regime where the twobody scattering length a 2 is large, In this regime, there are two kinds of two-body bound states: • the weakly bound dimer, of binding energy ≈ ħ 2 /(ma 2 2 ), which exists for a 2 > 0 • deeply bound dimer(s), of binding energy ≳ħ 2 /(mb 2 ), which exist if the interaction potential is deep enough (as in generic cold atom experiments).
We consider a stationary solution of the N -body Schrödinger equation in the zero-range regime 1/k typ ≫ b (41) and ( 36) was deduced from the expression which has the asymptotic behavior where 1/k typ is defined as the smallest scale of variation of the stationary wavefunction ψ in the region where all interparticle distances are ≫ b. 15  Let us first consider the case where there are no deeply bound states.For simplicity we also assume that the spectrum is discrete (which is the case in a trapping potential of infinite depth -e.g. a harmonic trap-or in a box with periodic boundary conditions).Then, in the zerorange regime (41), the zero-range model is valid for any stationary state ψ, in the sense that standard observables tend to their respective values within the zero-range model.This includes observables such as the energy, as well as N 2 (ϵ) and N 3 (ϵ) provided ϵ ≫ b [to reach the asymptotic regimes of Eqs.(1,3) one also needs ϵ ≪ 1/k typ ].
We turn to the experimentally relevant case where deeply bound dimers exist.These deeply bound dimers can be formed through recombination processes between three atoms.Let us denote by Γ 3 the number of such events per unit of time.The recombination products (the deeply bound dimer and the third atom) escape from the trapped gas, provided the trapping potential U (r) has a finite depth much smaller than the binding energy (∼ ħ 2 mb 2 ) of deeply bound dimers.In typical experiments, this condition holds, and other loss processes are negligible, which allows one to measure Γ 3 from the decay of the number of trapped atoms, Ṅ = −3 Γ 3 [11,[129][130][131][132].
In the zero-range regime, this decay is slow compared to the other timescales of the problem, as we will see.A standard way to describe such a slowly decaying state in quantum mechanics is to consider a quasi-stationary Gamow state, i.e., a solution of the Schrödinger equation with a complex energy and an outgoing-wave asymptotic behavior [133][134][135][136][137][138].Accordingly, we will consider a solution ψ of (39,40) with a complex E , and an outgoing-wave asymptotic behavior corresponding to the recombination products (a deeply bound dimer + an atom) flying apart towards large distances. 16For such a quasi-stationary state, in the zero-range regime, standard observables again tend to their respective values within the zero-range model. 17,18The threebody loss rate Γ 3 , however, is simply zero within the zero-range model.To compute Γ 3 , one thus needs to go beyond the zero-range model.As we will see, one only needs to do so for the threebody problem, in order to define a three-body parameter a 3 .We then find where C 3 can be evaluated within the zero-range model.Furthermore, breaking up Γ 3 into the sum of the two contributions Γ 2,1 and Γ 1,2 corresponding to ↑↑↓ and ↑↓↓ loss processes, we have We expect these relations to be asymptotically exact in the resonant zero-range regime (38,41). 15For example, for the ground state of the homogeneous unpolarized gas, for a 2 > 0, where k F is the Fermi momentum; for the ground state of a few particles in an isotropic harmonic trap of frequency ω, 1/k typ is ∼ a ho for a 2 < 0 and ∼ Min(a ho , a 2 ) for a 2 > 0, where a ho := ħ/(mω) is the harmonic oscillator length. 16An alternative approach would be to use the Lindblad equation.We expect that this would lead to the same result for the loss rate, as was checked for three-body losses for bosons in [139]. 17An appropriate normalization of the Gamow state will be given below in Eq. ( 58).Similarly, the expectation value of an observable Ô should be defined as 18 Within the zero-range model, it is convenient to add steep infinite walls to the trapping potential at the boundary of R trap , in order to have truly stationary states, thereby neglecting the exponentially suppressed evaporation process discussed in footnote 23.
To define the three-body parameter a 3 , we consider the zero-energy free-space solution of the Schrödinger equation (39,40) for three particles of spins ↑↑↓ and angular-momentum quantum numbers (l = 1, m) whose asymptotic behavior has the form in the region {b ≪ r i j ≪ |a 2 |, ∀i < j } where all interparticle distances are large compared to the range but small compared to the two-body scattering length.Here we have neglected the deep-dimer + atom outgoing wave, since it is proportional to the dimer wavefunction which is exponentially suppressed at distances ≫ b.The fact that a 3 does not depend on the quantum number m follows from rotational invariance of the interaction. Remarks: • Relation ( 42) is reminiscent of the known relation [44,139] between two-body loss rate and two-body contact 19,20 • Im a 3 must be negative, since the loss rate is positive.
• Typically one has the order of magnitude estimate a 3 ∼ b 2s , assuming that there is no extra fine-tuning. 21Therefore Γ 3 ∝ b 2s is small in the zero-range regime, as anticipated.• Based on heuristic arguments, it was already stated in [98] (see also [130,140,149]) that in the zero-range regime, the formation rate of deeply bound dimers is proportional to the probability of finding three particles at distances ≲ b, times ħ/(m b 2 ), that is, , with N 3 (b) evaluated within the zero-range model, and K a dimensionless prefactor that depends on short-range three-body physics.This statement is equivalent to (42), given the relation (3), with K = −8 s (s+1) Im a 3 /b 2s .The novelties of the present work are (i) to provide a derivation of the relation (42), and hence of the above statement from [98], (ii) to introduce the natural single parameter a 3 (instead of the two parameters K and b), and (iii) to generalize the relation to more complex interactions (in the following Section 3.1.2). • Let us consider the case of the homogeneous unpolarized zero-temperature unitary gas, of number density n.Introducing the three-body contact density C 3 := C 3 /V with V 19 Relation (46) concerns the situation where the states ↑ and ↓ populated in the gas are not the two energetically lowest atomic internal states, so that inelastic two-body collisions towards lower lying states are energetically allowed, and a 2 acquires an imaginary part. 20Relation (46) was obtained in [44] using the relation (84) and Im E = −ħ Γ 2 /2.It was rederived in [139] using the Lindblad equation.For completeness, we note that the relation (46) can also be derived by a flux computation.The main steps of this derivation are as follows.We consider a solution ψ of the zero-range model (13,14) with complex values of 1/a 2 and E .The corresponding time-dependent wavefunction is Ψ(t ) = ψ e −i E t /ħ , and we have Using (55,56), we obtain that Γ 2 equals N ↑ N ↓ times the limit when ϵ → 0 of the probability flux entering into the region {(r 1 , . . ., r N )|r < ϵ}, which can be simplified to (14) then yields (46). 21Indeed, the behavior (45) has to be matched at R of order b with the solution inside the potential, which typically imposes that the two terms in (45) are of the same order of magnitude for R ∼ b [97,140].For simplicity, we excluded here the special regime |a 3 | ≫ b 2s that corresponds to the vicinity of a three-body resonance (see [141,142], and [38,117,[143][144][145][146][147][148] for the mass-imbalanced case with 0 ≤ s < 1).Reaching this regime would require a second fine-tuning of the interaction, in addition to the first fine-tuning that causes |a 2 | ≫ b (for cold atoms, it would require a second control parameter of the interaction, in addition to the magnetic field used to tune a 2 to large values).We expect the relations (42,43,44) and the other results of this article to remain applicable in the three-body resonant regime provided the threebody parameter(s) remain(s) small (in modulus) compared to 1/k 2s typ .
the volume, dimensional analysis gives Therefore, as already found in [98], the timescale of three-body losses τ 3 := n/| ṅ| is of order τ F /(k F b) 2s , much larger than the thermalization time τ F ∼ m/(ħ k 2 F ), so that the gas remains at quasi-equilibrium.
• For bosons (and more generally in presence of the Efimov effect) it is commonly accepted that three-body losses can be described by making the three-body parameter complex [150][151][152].We have transposed this to the fermionic case (where the Efimov effect does not occur) by introducing a complex three-body parameter a 3 .The expression (42) of Γ 3 is reminiscent of the relation for bosons expressing Γ 3 in terms of C 3 and the inelasticity parameter (i.e. the phase of the three-body parameter) [72,139].An important difference is that in the fermionic case, in the zero-range regime, the three-body parameter is small, and can be set to zero when evaluating typical observables other than Γ 3 .• The notion of three-body parameter a 3 differs from the three-body scattering hypervolume D which was defined for bosons in 3D [153,154] and for various other cases [155][156][157][158]. Presumably, D could be defined also for the present case (two-component fermions in 3D), and as in [153][154][155][156][157][158] We turn to the derivation of (42,43,44).Our reasoning is similar to the bosonic case treated in App.B of [72], but the present case is significantly more complicated because we cannot work directly within the zero-range model.For R ≫ b, the a 3 /R s term in (45) is negligible compared to the R s term, consistently with the behavior (15) within the zero-range model.However this a 3 /R s term cannot be neglected to describe three-body losses, since the losses come from the non-zero imaginary part of a 3 , as we will see.Accordingly, for the Gamow state ψ, we go beyond the zero-range model and replace (15) with which we expect to hold provided where and whose appropriate choice will be discussed more precisely below.The purpose of ( 51) is to ensure that the interparticle distances within the triplet of particles 1, 2, 3 are much smaller than all other interparticle distances.
We use the shorthand notation X = (r 1 , . . ., r N ).The solution Ψ(X; t ) of the time-dependent Schrödinger equation [i ħ Ψ = H Ψ, with H given by ( 40)] associated with the Gamow state ψ It satisfies the continuity equation in terms of the probability current We express the loss event rate (i.e. the probability for a loss event to occur per unit of time) as where the region R should physically correspond to N atoms in the trap, and where There is some freedom in how to define R. Let R trap denote the "trapping region" of the potential U (r), which can be defined as the set of points r such that the classical trajectory of a particle with initial position r and zero initial velocity remains bounded. 22A possible definition of R would be R = (R trap ) N , but for later convenience, we choose R to be slightly smaller, by excluding configurations with two or three nearby particles: where In other words, (r 1 , . . ., r N ) ∈ R means that • all particle positions r i are inside the trapping region R trap • all distances between pairs of particles with opposite spin are ≥ d 2 • all hyperradii of triplets of particles with non-identical spins are ≥ d 3 .
We take Conditions (53,62) ensure that both d 2 and d 3 are ≪ 1/k typ .Hence, to leading order, the normalization integral (58) is independent of d 2 and d 3 , because R 2 and R 3 are negligibly small subsets of (R trap ) N .Furthermore, to leading order in the zero-range regime, • there exists a stationary state of the zero-range model whose wavefunction (normalized by the usual integral over the entire space) is close to the Gamow state ψ(X) for X ∈ R • in (48), B m (C; r 4 , . . .r N ) can be evaluated within the zero-range model, provided C, r 4 , . . .r N all belong to R trap ; indeed, the ∝ a 3 term is a small correction, and the condition (48) for the Gamow state must match with the condition (15) for the zerorange model.
Equations (54,57,58) directly yield which we will use for a consistency check in Sec.3.2.Here we will evaluate the loss rate by a flux computation.We use the notation for the probability flux through a surface S .In (57), we interchange the time derivative and the integration, and use the continuity equation ( 55) and the divergence theorem, which yields where ∂R is the boundary of R.
Here and in what follows, the differential surface vector d 3N−1 S appearing in ( 64) is oriented towards the exterior of R. Assuming that the trap depth is large enough for evaporation to be negligible 23 , we have where A visual representation is shown in Figure 1.On physical grounds, we identify Φ(S 2 ) as the twobody loss rate Γ 2 , and Φ(S 3 ) as the three-body loss rate Γ 3 , which determine the average number of lost atoms per unit of time: In the considered regime of small range and large trap depth, we expect that Γ 2 is given by the relation (46), and that in the considered case where a 2 is real, Γ 2 is negligible compared to Γ 3 so that Γ ≃ Γ 3 . 24 It remains to determine Γ 3 = Φ(S 3 ).We will use the notation From (60), we have S 3 = i < j :↑,k:↓ or i < j :↓,k:↑ the contribution from ↑↓↓ triplets. 23Evaporation is the process of an individual atom (or a weakly bound dimer for small positive a 2 ) escaping directly from the trapping region (because its energy is large due to a rare fluctuation, and/or it tunnels through the trappingpotential barrier).This process is exponentially suppressed in the limit where the trap depth is large compared to the typical energy per atom (subtracting the dimer binding energy contribution for a 2 > 0).For evaporation of individual atoms, the trap depth can be defined as , the evaporation rate is Φ(S trap ) with , and this rate is exponentially suppressed because J(X) with X ∈ S trap is exponentially suppressed. 24Indeed, in the zero-range limit, the reasoning of footnote 20 yields the expression (46) for Φ(S 2 ) = Γ 2 ; moreover, in the case where a 2 ∈ R, we expect Γ 2 ≪ Γ 3 , because there is no mechanism that would generate a non-negligible flux exiting R through S 2 (hence entering into R 2 ) and propagating in R 2 with an initial wavevector high enough to climb the trapping potential barrier and escape from (R trap ) N .r ρ Figure 1.Geometric illustration of the three-body loss process.The total decay rate is given by the probability flux exiting from the region R (grey shaded area).Neglecting the flux through S 2 (red straight lines) which corresponds to two-body losses, and the flux through S trap (dashed red lines) which corresponds to evaporation, the dominant contribution is the flux through S 3 (green circular arcs) which corresponds to three-body losses.The blue arrows represent the deep-dimer + atom outgoing wave, corresponding to a deeply bound dimer and an atom flying apart with a large relative momentum and escaping from the trap (this wave propagates in the region R 2 ).For the purpose of making this illustrative drawing two-dimensional, we considered N = 3 particles in one space dimension, and fixed the center-of-mass coordinate to C = 0.The positions of the three particles are then determined by the Jacobi coordinates r and ρ.The trapping region was simply assumed to be a symmetric interval around the origin.
By antisymmetry, each ↑↑↓ triplet gives the same contribution, so that (123) .(67) This can be rewritten as25 We then simplify this expression in two steps: (i) replace ψ with its asymptotic behavior (48), with B m evaluated within the zero-range model (ii) replace the integration domain S (123) by the entire region ∂B 123 .
Step (i) is justified since the conditions (49,50,51) hold except in a negligibly small domain of the integration variables C, r 4 , . . ., r N , Ω.
Step (ii) is justified because given (62), the condition X ∉ R 2 only excludes a negligibly small region of hyperangles Ω.We note that the order of these two steps is important. 26Using (17) then yields the result (43).The expression (44) of Γ 1,2 is derived analogously.
Let us now discuss in more detail the appropriate choice of the length d 3 .The condition (52) is actually not sufficient in order to have the behavior (48) of ψ.Indeed, we expect [based on the small-R expansion of the finite-energy solution J s (kR) of the hyperradial Schrödinger equation (93)] that in addition to the term R s in (48), there is a higher-order correction term of order R s+2 k 2 typ , which is negligible compared to a 3 /R s provided we take The condition d 3 ≪ 1/k typ is then automatically satisfied.Compatibility with (52) then requires (k typ b) 1/(s+1) ≪ 1, which is a quite stringent condition given the smallness of the exponent 1/(s + 1) ≃ 0.36.However, while this condition is necessary for (48), we do not expect it to be necessary for the final expression (42) of the three-body loss rate.
To complete our discussion of validity conditions, we now consider the contribution of the angular-momentum sector l = 0. Another condition to fulfill in order for (48) to be valid is that one can neglect the contribution coming from the l = 0 sector of the unitary three-body problem.Let us denote by s ′ := s l =0,n=0 = 2.166221977 . . . the smallest solution different from 2 of s ′ cos(s ′ π/2) + 4 sin(s ′ π/6)/ 3 = 0, and by φ ′ (Ω) := φ (l =0,n=0) (Ω) the corresponding l =0 hyperangular wavefunction.There is a higher-order correction to the r.h.s. of (48) given by R s ′ −2 φ ′ (Ω) times a function B ′ (C; r 4 , . . ., r N ).Requiring this ∝R s ′ correction to be negligible compared to the ∝a 3 /R s term in (48) would yield an additional condition on d 3 , but this is not necessary for the final result (42).Instead, what is truly necessary for (42) is that the contribution from the l =0 ing a minor error in an intermediate step in [72]).We need to evaluate We then perform the change of variables X −→ (C, R) , and rewrite the integrand as  26 If we would start by replacing S (123) with ∂B 123 ∩ (R trap ) N (keeping the original Gamow-state wavefunction ψ), then we would get a vanishing result for the total flux : The flux through ∂B 123 ∩ (R trap ) N \ R 2 (corresponding to the three-atom wave partially reflected from the small-R region) would be compensated by the flux through the complementary surface R 2 ∩ ∂B 123 ∩ (R trap ) N where the deep-dimer + atom outgoing-wave contribution to ψ gives the main contribution to the flux.This follows from the fact that for the three-body zero-energy scattering state Ψ m , the total flux ) vanishes (a related discussion can be found in [154]).

General interactions
In cold atom experiments, interactions are more complex than the simple model of Sec.3.1.1.Not only two-body interactions, but also three-body interaction are present.Moreover, • interactions are not necessarily rotationally invariant around any axis, due to the presence of the external magnetic field • interactions are not necessarily symmetric w.r.t.exchanging the role of ↑ and ↓.
As a minimal model including these features, we consider a two-body interaction potential V 2 (r) which may now depend on the direction of r, and a three-body interaction potential V 2,1 (resp.V 1,2 ) between triplets of particles of spins ↑↑↓ (resp.↑↓↓).The corresponding stationary N -body Schrödinger equation is Here R i j k denotes the Jacobi coordinates associated to particles i , j , k [defined by replacing the indices 1,2,3 by i , j , k in (88,89)].
The two-body interaction V 2 is still assumed to be resonant, |a We find that the ↑↑↓ and ↑↓↓ three-body loss-rates are given by Here, the m-resolved three-body contacts C (m) 2,1 and C (m) 1,2 are defined by where B m is related to the many-body wavefunction (in the zero-range limit) through Eq. ( 15), and similarly where Bm is defined in Eq. ( 18).Note that from Eqs. (17,19) we have 1 m=−1 C (m) 2,1 = C 2,1 and 1 m=−1 C (m) 1,2 = C 1,2 .The three-body parameters a (m) 2,1 are defined as follows.Setting the solution Ψ m of the zero-energy Schrödinger equation in free space for two ↑ and one ↓ particles with angular-momentum quantum numbers (l = 1, m) has the asymptotic behavior in the region where all interparticle distances are ≫ b and ≪ |a 2 |.The three-body parameters a (m) 1,2 are defined similarly in terms of the ↑↓↓ three-body scattering states.The relation ( 71) is derived by considering the probability current, in a completely analogous way to Eqs. (57,65,67,68) above, using the asymptotic behavior of the many-body wavefunction in the region (49,50,51) which is now given by The expression (72) of Γ 1,2 is derived analogously.

Discussion:
The parameters Im a (m) 2,1 and Im a (m) 1,2 are a priori unknown. 27In principle, one may hope to compute them theoretically by solving a sufficiently realistic three-body problem, but this is a difficult task [163][164][165]. 28Instead, one could determine them by measuring the three-body loss rate in situations where the three-body contacts C (m) 2,1 and C (m) 1,2 are known theoretically.A first possibility is to work with a small number N of particles in a microtrap [166][167][168][169] where the N -body wavefunction can be computed numerically with good accuracy [106,170,171] so that the three-body contacts C (m)  i , j could be calculated.Measuring Γ 3 in six different states and inverting (71,72) would allow to determine the six parameters Im a (m)  i , j .The case of three particles in an isotropic harmonic trap is particularly simple: The problem is analytically solvable [115,172], and if one prepares one of the six degenerate ground states corresponding to a given quantum number m ∈ {−1; 0; 1} and is the only non-zero three-body contact, so that Γ 3 is simply proportional to Im a (m) N ↑ ,N ↓ .Explicitly, taking 27 One may expect a small relative difference between a (0) i , j and a (±1) i , j , and an even smaller one between a (1)  i , j and a (−1) i , j , similarly to the m-dependence of the two-body p-wave scattering volume not too close to a p-wave Feshbach resonance [160][161][162]. 28For such a computation of the three-body parameters, one may need to take into account that an atom has more than two relevant internal states |ν〉.However, we expect the general relations (71,72) to remain applicable.Indeed, we expect that in typical experiments, the atoms mainly occupy two internal states, that we can label ν = ↑ and ↓, and if all interatomic distances r i j ≫ b, the many-body wavefunction Φ(r 1 , ν 1 ; . . .; r N , ν N ) is non-negligible only if all ν i belong to {↑; ↓}, in which case Φ is given to good accuracy by antisymmetrizing the wavefunction ψ(r 1 , . . ., r N ) of the zero-range model.
for example N ↑ = 2 and N ↓ = 1, we get 2,1 (mω/ħ) s × 1.266322 (in agreement with the scaling given in [172]). 29 A second possibility is to work with a homogeneous (or locally homogeneous) unpolarized gas at equilibrium, for which the six three-body contacts C (m)  i , j are all equal, as shown in Appendix B. This gives with One could then determine Im ā3 by measuring Γ 3 in a weakly correlated regime, where the asymptotic behavior of the three-body contact density C 3 can be computed exactly.A first option is the non-degenerate regime, where we have computed C 3 for negative or infinite scattering length. 30Other options are the weakly interacting regimes where k typ |a 2 | is small. 31

Three-body contribution to the finite-range correction to the energy
In this Section, we study the corrections to the zero-range model's energy coming from the finite range of the two-body interaction and/or an additional three-body interaction.The stationary Nbody Schrödinger equation is again given by (70).We consider the case where there are no deeply bound dimers, so that the three-body parameters are real.The zero-range model is approached in the zero-range regime where b In particular, each eigenergy E of ( 70) approaches a corresponding eigenenergy E ZRM of the zero-range model.We are interested in the energy difference between the finite-range and zero-range models, We find that in the zero-range regime, δE ≃ δE 2 + δE 3 (80) plus higher-order corrections, where δE 2 is given to leading order by [56] with r e the effective range of the two-body interaction V 2 , while δE 3 is given to leading order by 29 For the excited state whose energy is 2q ħω above the ground state energy, the value of C (m) 2,1 is multiplied by q+s s = Γ(q+s+1) Γ(s+1) q! (which is a growing function of q, meaning that the growing delocalization in the trap is overcompensated by the growing penetration under the s 2 /R 2 barrier). 30For the homogeneous unpolarized unitary gas of density n, we obtain in the nondegenerate limit [X.Leyronas and F. Werner, "Three-body contact for fermions.II.Non-degenerate limit", to be submitted]. 31Although the three-body parameters depend on magnetic field, this dependence is smooth if no three-body resonance is crossed, and may be neglected in the vicinity of a given Feshbach resonance.
in terms of the three-body contacts and three-body parameters.If the two-body and threebody interaction potentials are rotationally invariant and V 2,1 = V 1,2 , then the six three-body parameters are all equal to a single a 3 and the expression simplifies to Remarks: • δE 2 (resp.δE 3 ) comes from configurations where 2 (resp.3) particles are close to each other.• The relations (82,83) are reminiscent of the relation [36] which holds within the zero-range model. 32 • Typically, a 3 is of order b 2s , which gives δE 3 ∼ b 2s .The latter scaling was already stated (for lattice models) in [173,174] and was derived at the level of the third virial coefficient in [175].the three-body correction to the energy is of higher order than the two-body correction. 33,34 • Let us assume that ( 83) can be analytically continued to complex a 3 , with C 3 still evaluated within the zero-range model.We then recover the expression (42) of the three-body loss rate, simply by substituting Im E = Im δE 3 into (63), and using the fact that Γ = Γ 3 and C 3 ∈ R. Similarly, applying this procedure to (82) yields Γ 3 in agreement with the sum of Eqs.(71,72).
We have ∆ 0 = 0, assuming for simplicity that the interaction potentials (V 2 , V 2,1 and V 1,2 ) are 32 A relation which is more directly analogous to (84) can be formulated in the mass-imbalanced case, see (128). 33In the mass-imbalanced case discussed in App.D, the situation is reversed beyond a critical mass ratio, as was already noted in [175]. 34Apart from the leading-order term (81) which is of order b, there are also higher-order contributions to δE 2 , which we expect to contain no term of order b 2s (since the three-body physics does not enter in δE 2 ) but only integer powers of b.Accordingly, the ∝ b 2s contribution to δE should be entirely given by ( 82) or (83).
finite everywhere (i.e.excluding hard-wall potentials) and do not diverge too quickly at short distances. 35Hence Since the potentials are short-ranged, the integral is dominated by the contributions from three regions, outside of which V becomes negligible: • two particles of spins ↑↓ are nearby (at distance ≲ b 2 ) while all other interparticle distances are ≫ b • there is a triplet of particles of spins ↑↑↓ which are nearby (their hyperradius is ≲ b) while all interparticle distances other than the ones within that triplet are ≫ b • there is a triplet of particles of spins ↑↓↓ which are nearby (their hyperradius is ≲ b) while all interparticle distances other than the ones within that triplet are ≫ b.
Denoting the contributions of these regions by ∆ 2 , ∆ 2,1 and ∆ 1,2 respectively, we thus have The two-nearby-particle contribution ∆ 2 is given by the r.h.s. of ( 81), as shown in Appendix C, in agreement with [56].It remains to evaluate the contribution ∆ 2,1 coming from three nearby particles of spins ↑↑↓.We introduce a length d 3 satisfying (52,53,69).∆ 2,1 is then given by the contribution to the integral (85) coming from the region of (r 1 , . . ., r N ) such that there is a triplet of particles (i , j , k) of spins (↑, ↑, ↓) of hyperradius R i j k < d 3 while all interparticle distances other than the ones within the triplet (i , j , k) are ≫ d 3 (the result does not depend on the value of d 3 within the range (52,53), as we will see).In the region where these conditions are met for the triplet (1, 2, 3), i.e. when (49,51) hold, we expect a factorization as well as [cf.(15)] with Ψ (0) m (R) := R s−2 φ m (Ω).Furthermore, we can approximate B (1)  m by B m in (86), given that we are in the zero-range regime and all distances between the points C, r 4 , . . ., r N are ≫ b.Also using the fact that each triplet gives the same contribution by fermionic antisymmetry, we get where W m,m ′ := and V (R) was defined in (75).To evaluate W m,m ′ , we first use the Schrödinger equation (76) to replace V by ħ 2 m ∆ R .Since we also have and applying the divergence theorem yields We can then use the asymptotic behavior (77) of Ψ m (indeed, for R = d 3 ≫ b, all three interparticle distance are ≫ b except in a small region of hyperangles).This gives W m,m ′ = ħ 2 m 2 s δ m,m ′ a (m) 2,1 .Substituting this into (87) and using (17 The same reasoning gives an analogous expression for ∆ 1,2 , which yields the final expression (82) for

Summary and outlook
We have shown that the three-body contact C 3 is a useful concept for the fermionic N -body problem with resonant interactions, in the standard regime of mass ratio where there is no Efimov effect.Within the zero-range model, the three-body contact controls the number of nearby triplets, the third order density correlation function at short distances, the tail of the centerof-mass momentum distribution of nearby pairs, and the tail of the two-particle momentum distribution; C 3 also has a simple expression in terms of the N -body wavefunction in the limit where three particles are nearby.Beyond the zero-range model, for a small finite interaction range, we introduced a small three-body parameter a 3 , and we showed that the formation rate Γ 3 of deeply bound dimers by three-body recombination equals C 3 Im a 3 times an explicit prefactor; we also showed that the finite-range correction to the energy has a contribution equal to − ħ 2 C 3 a 3 times the same prefactor.With respect to the relation between Γ 3 and the number of nearby triplets stated in [98], the present work adds a derivation, and an expression of the prefactor in terms of Im a 3 .
Furthermore, we considered the general case where there are two different contributions C 2,1 and C 1,2 to the three-body contact, corresponding to the spin configurations ↑↑↓ and ↑↓↓ for the associated three-body problem, which can be further broken up into the contributions C (m) 2,1 and are the ones for finite-range interactions, in the general case where interactions are not invariant under rotation and under exchange between ↑ and ↓.In this case, the three-body parameter also depends on the spin and angular momentum indices.Nevertheless, for a homogeneous unpolarized gas, Γ 3 simply equals C 3 Im ā3 times an explicit prefactor, with ā3 the mean three-body parameter.For the unitary gas in the non-degenerate regime, we announced the result of our computation of C 3 (see footnote 21), which would allow one to determine Im ā3 by measuring Γ 3 .Measuring Γ 3 in the low-temperature regime would then allow one to experimentally test the power-law C 3 = ζ 3 n (2s+5)/3 and determine the associated many-body parameter ζ 3 whose computation is an open theoretical challenge.zero center-of-mass momentum, ψ only depends on the relative positions between the three particles, which can be parameterized by the Jacobi coordinates It is convenient to introduce the six-dimensional vector whose norm R = ∥R∥ = r 2 + ρ 2 is the hyperradius, while its direction R/R can be parameterized by five hyperangles denoted collectively by Ω, The three-body Schrödinger equation then writes with the contact condition: there exists A such that (the contact condition between particles 2 and 3 then automatically holds by antisymmetry).
At the unitary limit a 2 = ∞, this contact condition scale invariant.Accordingly, it only acts on the hyperangles; explicitly, it can be expressed as ∂ ∂α sin α ψ α=0 = 0 where α = arctan(r /ρ).As a result, the unitary three-body problem is separable in hyperspherical coordinates: One can look for eigenstates of the factorized form (Ω) (where the factor 1/R 2 is introduced for later convenience), and the three-body problem (91,92) separates into • a hyperangular problem, defined by together with the contact condition and the antisymmetry constraint.
Here T Ω is a differential operator acting on the hyperangles, defined by There is a discrete spectrum of values for s 2 , which are all real and positive in the present case of equal-mass fermions, so that we can take s real and positive; we denote the set of allowed s by {s ν } where ν is a discrete index.The corresponding hyperangular eigenfunctions {φ ν (Ω)} form an orthonormal basis for the hyperangular scalar product where d 5 Ω denotes the differential solid angle in six-dimensional space, d 6 R = d 5 Ω R 5 d R ; this can be deduced from the self-adjointness of the unitary three-body problem in an isotropic harmonic trap and can also be checked by explicit analytical calculations in the l = 0 subspace [178] (for a mathematical of self-adjointness in free space, see Refs.[179][180][181]).
The hyperangular problem (94,95) is analytically solvable, with two types of solutions.The first type are common solutions of the unitary and non-interacting hyperangular problems, whose wavefunction φ(Ω) vanishes when two particles approach each other (i.e. for α → 0); the corresponding s take integer values.The second type are the following truly interacting solutions: • For each l , let us denote by {s l ,n ; n ∈ N} the allowed values of s, in increasing order (s l ,0 < s l ,1 < . . .).The index ν can be identified with the set of quantum numbers (l , m,n).All the eigenfunctions φ (l ,m,n) with −l ≤ m ≤ l correspond to the same s l ,n , so that each s l ,n is 2l + 1 times degenerate.For l = 1, the s l ,n are the real positive solutions different from 1 of the transcendental equation ( 4) [the solution s = 1 should be discarded since it leads to an identically vanishing ϕ(α)].For arbitrary l , the s l ,n also solve transcendental equations (see, e.g., [38] and refs.therein); the smallest of all s l ,n is s l =1,n=0 , which is denoted by s for short throughout the article.
The value of the normalization constant, for l = 1, is We computed this value as follows.We have N = 3/J (s) with We evaluate the sum over m thanks to the identity: for any unit vectors û and û′ .This yields where α ′ and ρ′ are obtained from α and ρ by permutation of particles 1 and 2. The integrand in (101) only depends on α and α ′ , since ρ . Hence the integrand in ( 101) is a function of α and u.To evaluate the integral we use the formula 36 More explicitly, the considered wavefunction is an eigenstate of L 2 with eigenvalue l (l + 1)ħ 2 , and of L z with eigenvalue mħ, where L is the relative angular momentum of the 3 particles, L = −i ħ(r × ∇ ∇ ∇ r + ρ × ∇ ∇ ∇ ρ ).Note that the total angular momentum of the three particles is the sum of their relative and center-of-mass angular momenta: 37 Like the standard hyperspherical harmonics, the φ ν are eigenstates of the Laplacian on the hypersphere, Eq. ( 94); but they also satisfy the unitary-limit contact condition Eq. ( 95) (together with fermionic antisymmetry) which leads to non-integer eigenvalues s 2 .Félix Werner and Xavier Leyronas π/2 0 dα sin 2 (α) cos 2 (α) dr d ρ f (α, r, ρ).Since the integrand is independent of r and of the azimuthal angle of ρ w.r.t.r, we can integrate over them, which just gives a factor 4π × 2π.We are left with the integral over α and u.The change of variable u α ′ then yields where D is the domain α, α We evaluate analytically the integrals over α and α ′ for the first and second term of (102); for the third term, we evaluate analytically the integral over α ′ , and perform numerically the remaining integration over α.

Appendix B. Homogeneous gas
Let us show that for the homogeneous gas at equilibrium, C (m) 2,1 is independent of m.We will use the following lemma: Let Ô be an operator acting on functions of Ω, such that with θ the Heaviside function.This follows from Eq. ( 15) and the fact that the φ m (Ω) are orthonormal.
Injecting the ansatz (106) into the N -body Schrödinger equation ( 70) yields, in the region R 13 , with up to corrections that are negligible in the zero-range regime, as was checked in [56].We omitted the dependence of E on (c; r 2 , r 4 , . . ., r N ) to alleviate notations.We neglected all interaction potentials other than V 2 (r ), because we are far outside their ranges: In the considered region R 13 , all interparticle distances other than r ≡ r 13 , and hence also all triplet hyperradii, are ≫ b.
For the same reason, we will replace V by V 2 in (105).
For the zero-range model, we also expect a factorization of the many-body wavefunction in R 13 , where χ (0) E (r ) is the s-wave two-body scattering state at energy E for the zero-range model, i.e. the solution of with the contact condition: χ (0) E (r ) = 1/r − 1/a 2 + O(r ) for r → 0.Here we normalized χ (0) E in such a way that the function A in ( 108) is the same one than in (14).The solution is where f (0) k = −(1/a 2 +i k) −1 is the scattering amplitude of the zero-range model, and k := mE /ħ, with the determination k = i −mE /ħ if E < 0. We note that for r → 0, from the Taylor expansion of (110), we get a subleading singular contribution to the asymptotic expansion of ψ (0) given by −k 2 r A/2, in agreement with Eqs.(135,136) of [56].From this we can infer that k is typically ≲ k typ .
Scattering theory gives the large-distance behavior of the finite-range scattering state χ E : where f k is the s-wave scattering amplitude associated to the interaction potential V 2 (r ).
Here we normalized χ E in such a way that χ E (r ) ∼ 1/r for E → 0 and r → ∞, in agreement with the assumption that the same function A appears in (108) and (106).
Next, in (105), we can thus replace V by V 2 and substitute (108,106), which yields where It remains to evaluate W .Using (107,109), we rewrite it as This gives, by the divergence theorem, W = 4πħ 2 m W with W := χ (0) we can directly substitute (110,111) and their derivatives, which yields Since kb 2 ≲ k typ b 2 ≪ 1, we can use the low-energy expansion of the scattering amplitude: e /2 + . . . in the limit kb 2 → 0, where r e is by definition the effective range.This gives W ≃ 2πE r e .Substituting this into (112), we conclude that ∆ 2 is indeed given by the r.h.s. of (81).

Appendix D. Mass-imbalanced case
In this Appendix we consider the case where the ↑ and ↓ fermions have different masses, m ↑ > m ↓ for definiteness.Experimentally, this is realized in a mixture of two fermionic species, such as 40 K-6 Li [182], 161 Dy-40 K [183], or 53 Cr- 6 Li [184].In Section D.1, we extend the relations obtained in the main text to the mass-imbalanced case.In Section D.2 we show that, when the mass ratio exceeds a critical value, a conceptual simplification takes place: A generalized zero-range model can be introduced, within which relations involving the three-body parameters can be formulated and derived more directly.The results of this Appendix are valid if m ↑ /m ↓ is smaller than the threshold where the five-body Efimov effect appears [103], which implies that the fourbody Efimov effect [102] and three-body Efimov effect [176] do not take place either; moreover m ↑ /m ↓ should not be too close to the four-body and five-body Efimov thresholds, as discussed in Section D.3.Obviously, in the N -body Schrödinger equation, the mass m is replaced by the mass m i of particle i (m i = m ↑ or m ↓ depending on the spin of particle i ): for the zero-range model, and for the finite-range model.It will prove convenient to introduce the angles ϕ and φ related to the mass ratio by The definitions of the center-of-mass and Jacobi coordinates should be generalized to the unequal-mass case.The center-of-mass of particles 1 and 3, which appears in the two-body contact condition (14), The continuity equation and the probability current are still given by (55,56) provided we define While in the mass-balanced case, there was a single exponent s associated to the unitary three-body problem, in the mass-imbalanced case we have two exponents s and s, associated respectively to the ↑↑↓ and ↑↓↓ unitary three-body [97,176].As a function of m ↑ /m ↓ , s is continuously decreasing and vanishes at the Efimov-effect threshold m ↑ /m ↓ = 13.6069 . . ., while s is continuously increasing and tends to 2 for m ↑ /m ↓ → ∞.

D.1. Extension of the relations from the mass-balanced case
The number of nearby ↑↑↓ triplets is still given by whereas the number of nearby ↑↓↓ triplets is now given by For our convention m ↑ > m ↓ , we have s < s.Hence the total number of nearby triplets N 3 (ϵ) is dominated by the ↑↑↓ contribution, so that The remarks at the end of Sec.2.1 remain valid.In particular, the bunching effect due to the zerorange interactions still overcompensates the antibunching effect due to Pauli exclusion, both for N 2,1 (ϵ) and N 1,2 (ϵ), because s < s < 2.
When three particles of spins ↑, ↑, ↓ approach each other, the asymptotic behavior of the manybody wavefunction is still given by (15).This yields [using the Jacobian |∂(r 1 , r 2 , r 3 )/∂(C, R)| = cos 3 When three particles of spins ↑, ↓, ↓ approach each other, the asymptotic behavior of the manybody wavefunction is given by ( 18) with s replaced by s, i.e.
This yields The triplet correlation functions satisfy The leading large-momentum tail of NP (K ) is given by where M is still given by the expression (24) in terms of N , and N now stands for the normalization constant of the unitary hyperangular wavefunction of the ↑↑↓ problem. 38 The two-particle momentum distribution has the tail lim The expression (71) for the ↑↑↓ three-body loss rate remains valid.In the expression (72) for the ↑↓↓ three-body loss rate, s is replaced by s Since the ↑↓↓ three-body parameters a (m) 1,2 are now of order b 2 s whereas the ↑↑↓ three-body parameters a (m) 2,1 are of order b 2s , the ↑↓↓ three-body loss rate Γ 1,2 ∝ b 2 s is negligible (in the zero-range regime) compared to the ↑↑↓ three-body loss rate Γ 2,1 ∝ b 2s .Hence The expression (81) for the two-body contribution to the energy correction remains valid.For the three-body contribution to the energy correction, we obtain because the leading-order contribution again comes from ↑↑↓ triplets. 39 A peculiar situation takes place for s < 1/2 (i.e. for m ↑ /m ↓ > 12.313 . . .): δE 3 ∝ b 2s dominates over δE 2 ∝ b ; i.e., the finite-range correction mainly comes from configurations with three nearby particles, rather than two nearby particles.This was already pointed out in [175] at the level of the third virial coefficient.Few-body and many-body numerical computations are often performed with finite-range interactions and extrapolated to the zero-range limit, see e.g.[146,185] and [28,62,65,67,111,112,173,174,[186][187][188][189][190][191][192] respectively; to accurately perform such b → 0 extrapolations in the regime s < 1/2, it may be important to include the b 2s scaling. 40,41  38 The contribution from the ↑↓↓ three-body problem gives rise to a higher-order subleading tail of NP (K ), given by C 1,2 K −2 s−4 ×32 π 3 4 s 3 − s−1/2 ( s +1) Γ( s +2) 2 sin 2 ( sπ) Ñ 2 with Ñ the normalization constant of the unitary hyperangular wavefunction of the ↑↓↓ problem. 39Indeed, the contribution to δE 3 from ↑↑↓ triplets of nearby particles, given by the r.h.s. of (126), is ∝ b 2s , whereas 40 Naturally, the ∝ b 2s scaling may be hard to distinguish from the regular ∝ b scaling if s is close to 1/2 (see e.g.Fig. 9(b) of [146]). 41We take this opportunity to recall two other subtleties relevant to zero-range extrapolations (i.e. to continuum extrapolations in the case of lattice models, where the interaction range is set by the lattice spacing).Being due to breaking of Galilean invariance, they arise for any mass ratio.The first subtlety, reported in [38,56,187], is that for lattice models (where the interaction range b is set by the lattice spacing), the ∝ r e term given by the r.h.s. of ( 81) is not the only ∝ b

D.2. Relations within the generalized zero-range model
In this Section, we assume that m ↑ /m ↓ is larger than 8.619 . . ., so that s < 1.This ensures that a three-body wavefunction diverging as R −s−2 at small R remains square integrable.This allows one to define a generalized zero-range model (GZRM) by supplementing the two-body contact condition ( 14) (involving the two-body scattering length) by the following three-body contact condition (involving the three-body parameters): There exists B m such that in the limit R → 0 where particles 1,2,3 approach each other, while keeping fixed their hyperangles Ω, their center-of-mass C, and the positions of the other particles r 4 , . . ., r N .By antisymmetry, Eq. ( 127) imposes a similar condition when any triplet of particles with spins ↑↑↓ approach each other.Apart from the two-body and three-body contact conditions (14,127), the GZRM is defined (like the standard zero-range model) by the Schrödinger equation without any interaction potential (113).Note that for vanishing three-body parameters, the GZRM reduces to the standard zero-range model (ZRM).

Remarks:
• If the underlying microscopic interactions are rotationally-invariant, then a (m) 3 does not depend on m.
• In the present GZRM, in the limit where three particles of spins ↑↓↓ approach each other, the asymptotic behavior of the wavefunction has the same form (120) than for the ZRM.
In other words, the behavior ψ ∝ R −2− s is forbidden.Such a wavefunction would not be a normalizable at small R, since s is always larger than one. 42This fact does not prevent one from rederiving the relations (125,126), which come from configurations with nearby ↑↑↓ triplets.• Extensions of the GZRM to the regime s ≥ 1, which were introduced recently [142], are beyond the scope of this work.• For mathematically rigorous studies of the GZRM in the three-particle case, see [201,202].contribution to the two-body finite-range energy-correction δE 2 : There is a second contribution to δE 2 , proportional to a parameter R e (that parameter R e quantifies the dependence of the two-body vacuum T-matrix on the center-of-mass momentum, which arises from lattice-induced breaking of Galilean invariance).Therefore, if one wishes to cancel the ∝ b term in the b → 0 continuum extrapolation, one needs to tune to zero not only r e (as done in [173,174,193] and for one of the dispersion relations considered in [112]) but also R e (and if the dispersion relation has a cusp at the edge of the Brillouin zone, there is a third ∝ b contribution to δE 2 , in a finite box with periodic boundary conditions [56]).The second subtlety, reported in [56] and further evidenced in [191], arises if one restricts single-particle momenta to a ball of radius Λ (with Λ = π/b for lattice models, with b the lattice spacing), i.e. if one takes a dispersion relation equal to +∞ for momenta larger than Λ, as was done in [188,[194][195][196][197][198][199][200]: The universal zero-range model is not obtained in the limit Λ → ∞, because the hard cutoff induces a dependence of the two-body T-matrix on the center-of-mass momentum, and this dependence surprisingly survives for Λ → ∞.We expect the same problem in [173,174] where such a spherical cutoff was also used. 42Hence there are no ↑↓↓ three-body parameters within the GZRM, which is why we denoted the ↑↑↓ Here the derivative w.r.t. a (m) 3 is taken at fixed value of a 2 and of the other a (m ′ ) 3 with m ′ ̸ = m.
The three-body contacts C (m) 2,1 are still defined by (119).Remark: Within the ZRM, the derivative of the energy w.r.t. a 2 is given by C 2 , see (84).Within the GZRM, this relation also holds, with d E /d (−1/a 2 ) replaced by the partial derivative ∂E /∂(−1/a 2 ) taken at fixed a (m)  3 .This can be justified by using the derivation presented in Sec.IV.C of [56] (at fixed three-body parameters, there is no additional contribution to the energy variation coming from nearby triplets of particles, as we will see below).
On the other hand, from the Schrödinger equations for ψ and ψ, As we will see below, there is a contribution δ 0 to δ coming from the configurations where particles 1, 2 and 3 are close to each other.By symmetry, there is an identical contribution from the configurations where any set of three particles with spins ↑↑↓ are close to each other.Hence δ equals δ 0 times the number of such three-particle sets, N ↑ (N ↑ − 1)N ↓ /2.To evaluate δ 0 , we only need to keep the terms i = 1, 2, 3 in Eq. ( 129), which gives, after the change of coordinates (r 1 , r 2 , r 3 ) −→ (C, R, Ω), where F := R 2 ψ, F := R 2 ψ, and T Ω is defined by (96).Since the two-body scattering length is the same for ψ and ψ, we only need to keep the terms involving derivatives with respect to R in Eq. ( 130). 44Transforming the integral over R into a boundary term at R → 0, we then get The result (128) follows by using the three-body contact conditions [given by Eq. ( 127) for ψ, and the same condition with ã(m)
To derive this relation, we proceed very similarly to the bosonic case of [72].We consider a stationary state of the GZRM, i.e. a solution ψ(r 1 , . . ., r N ) of the stationary Schrödinger equation (113) together with the two-body and three-body contact conditions (14,127).The corresponding solution of the time-dependent Schrödinger equation is Ψ(r 1 , . . ., r N ; t ) = ψ(r 1 , . . ., r N )e −i E t /ħ , and d 3 r 1 . . .d 3 r N |Ψ(r 1 , . . ., r N ; t )| 2 , with ψ normalized to unity.Excluding from the integration domain the regions where R i j k < ϵ and taking the limit ϵ → 0, the continuity equation and the three-body contact conditions leads to the final expression, where all extra mass-ratio dependent factors, arising e.g. from (115), divide out.

D.3. Validity conditions
Let us denote by s j ↑ , j ↓ the scaling exponent of the unitary ( j ↑ + j ↓ ) body problem with j ↑ particles of spin ↑ and j ↓ particles of spin ↓ (so that s 2,1 ≡ s and s 1,2 ≡ s ).We expect the relations (116,118,123,124) to be valid under the condition Min s j ↑ , j ↓ ; j ↑ ≥ 2, j ↓ ≥ 1, j ↑ + j ↓ > 3 > s.Indeed, we expect a contribution from the ( j ↑ + j ↓ ) body problem to N 2,1 (ϵ) at small ϵ [resp.to NP (K ) at large K ] scaling as ϵ 2s j ↑ , j ↓ +2 resp.1/K 2s j ↑ , j ↓ +4 , which dominates over the contribution from the three-body problem when s j ↑ , j ↓ < s.Similarly, we expect relation (117) to be valid under the condition Min s j ↑ , j ↓ ; j ↑ ≥ 1, j ↓ ≥ 2, j ↑ + j ↓ > 3 > s.The former condition breaks down when m ↑ /m ↓ is near the thresholds where the four-body [102] and five-body [103] Efimov effects appear, where s 3,1 and s 4,1 become small.Based on existing data, we expect that both of the above validity conditions are satisfied at least in the range m ↑ /m ↓ ≤ 10.Indeed, in this range, we have s 3,1 > s and s 4,1 > s [103], s 5,1 > s [104], s 2,2 > s [146], s 1,3 > s [103], and the trends of available data suggest that the conditions will also hold for larger values of j ↑ or j ↓ .We conservatively restricted to m ↑ /m ↓ ≤ 10 because s 2,2 was not computed beyond this range, but the conditions s 3,1 > s and s 4,1 > s actually hold up to at least m ↑ /m ↓ = 13.2 [103].

Appendix E. Generalization to statistical mixtures and non-stationary states
Many of the relations derived for stationary states in the main text and in Appendix D are directly generalizable to non-stationary states and statistical mixtures, similarly to the relations involving the two-body contact [35,36,56].Indeed, Eqs. (15,17,18,19) remain valid for any non-pathological linear combinations of stationary states (including solutions of time-dependent problems, where a 2 and the trapping potential U can depend on time); and relations (3,5,6,11,12,23) remain true for arbitrary non-pathological statistical mixtures of such pure states (including the case of thermal equilibrium).Here, non-pathological means that the occupation probabilities of stationary states should decay sufficiently quickly (which includes the simple case where only a finite number of states are populated); more specifically, for a pure state, it is necessary that Eqs.(14,48) still hold, while for a statistical mixture ρ = i c i |ψ i 〉〈ψ i |, it is necessary that the three-body contact of the mixture [defined by Eq. (3)] equals i c i C 3,i with C 3,i the three-body contact of state ψ i .Moreover, at thermal equilibrium, the thermally averaged loss rates are given by (42,43,44) for simple interactions; for more general interactions they are given by (71,72), and by (79) for a homogeneous unpolarized gas.Furthermore, (82,83) remain valid in the canonical ensemble, with δE 3 the energy difference (between finite-range and zero-range models) taken at fixed entropy (which equals the free-energy difference at fixed temperature).
As for the wavevector k typ that appears in the validity conditions, it should be defined by the same procedure as before and then taking the maximum over all significantly populated eigenstates; for example for the balanced unitary gas at thermal equilibrium, k typ = max(k F , k T ) with k T the thermal wavevector, defined by ħ 2 k 2 T /(2m) = k B T .

14
Typically, b is set by the true range b 0 of the potential (i.e. the length such that V 2 (r ) decays quickly for r ≫ b 0 ).More generally, b = Max(b 0 , |r e |) where r e is the effective range.
, D would govern the asymptotic behavior of the three-body zero-energy wavefunction at distances ≫ |a 2 | (while a 3 governs the regime of distances ≪ |a 2 | and ≫ b) and Γ 3 would have a simple expression in terms of Im D in the weakly interacting regime k typ |a 2 | ≪ 1 (whereas (42) remains valid in the strongly correlated regime k typ |a 2 | ≳ 1).On the other hand, a 3 is only defined for resonant interactions (|a 2 | ≫ b) while D also exists for non-resonant interactions.
By the divergence theorem, the integral over C of the ∇ ∇ ∇ C term vanishes, while the integral over R of the ∇ ∇ ∇ R term yields the result Q = − 3 3 4 d 3 5

2 | ≫ b 2 (
with a 2 and b 2 the scattering length and the range of V 2 ).The three-body interaction potential is assumed to have a finite range b 3 , in the sense that it decays quickly at hyperradii larger than b 3 .We consider the zero-range regime whereb := Max(b 2 , b 3 )is much smaller than |a 2 | and 1/k typ , where 1/k typ is defined as the smallest scale of variation of the wavefunction ψ in the region where all interparticle distances are ≫ b.For alkali atoms near an open-channel dominated Feshbach resonance, b 2 is set by the van der Waals length[159], which is ≪ 1/k typ in typical cold-atom experiments.Instead of a single three-body parameter a 3 , there are in general six three-body parameters a (m) 2,1 and a (m) 1,2 , with m ∈ {−1, 0, 1} the angular momentum quantum number.

D. 2 . 1 .
three-body parameters by a (m) Derivative of the energy with respect to the three-body parameters Within the GZRM, the derivatives of the energy w.r.t. the three-body parameters are given by the three-body contacts,