# Three-body contact for fermions. I. General relations

Félix Werner, Xavier Leyronas

► To cite this version:

Félix Werner, Xavier Leyronas. Three-body contact for fermions. I. General relations. Comptes Rendus. Physique, 2024, 25 (G1), pp.179-218. 10.5802/crphys.181. hal-03882472v2

HAL Id: hal-03882472

https://hal.science/hal-03882472v2

Submitted on 23 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Physique

Félix Werner and Xavier Leyronas

Three-body contact for fermions. I. General relations

Volume 25 (2024), p. 179-218

https://doi.org/10.5802/crphys.181

<figure><img src="image_1.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
Comptes Rendus
Physique

Félix Werner and Xavier Leyronas

Three-body contact for fermions. I. General relations

Volume 25 (2024), p. 179-218

https://doi.org/10.5802/crphys.181

<figure><img src="image_2.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
Research article / *Article de recherche*

# Three-body contact for fermions.

## I. General relations

*Contact à trois corps pour les fermions.*

### I. Relations générales

Félix Werner<sup>●,a</sup> and Xavier Leyronas<sup>●,b</sup>

**Abstract.** 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 ↑ and ↓ 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.

**Résumé.** Nous considérons le gaz de Fermi résonnant, à savoir des fermions avec deux états internes à trois dimensions avec des interactions à courte portée de grande longueur de diffusion. Nous introduisons une quantité, le contact à trois corps, qui détermine plusieurs observables. Pour le modèle de portée nulle, le nombre de triplets de fermions proches, la queue de la distribution selon l'impulsion du centre de masse des paire de fermions proches, ainsi que la queue de la distribution en impulsion à deux particules, sont exprimées en termes du contact à trois corps. Pour une portée non nulle, le taux de formation de dimères fortement liés par recombinaison à trois corps, ainsi que la contribution à trois corps à la correction de portée finie à l'énergie, sont exprimées en termes du contact à trois corps et d'un paramètre à trois corps. Ce paramètre à trois corps, qui tend vers zéro dans la limite de portée nulle, est défini *via* le comportement asymptotique de l'état de diffusion d'énergie nulle à des distances intermédiaires entre la portée et la longueur de diffusion à deux corps. En général, le contact à trois corps a différentes contributions repérées par des indices de spin et de moment cinétique, et le paramètre à trois corps peut dépendre de ces indices. Nous incluons aussi la généralisation à des masses différentes pour les particules ↑ et ↓. Par rapport à la relation donnée dans [Petrov, Salomon et Shlyapnikov, PRL 93, 090404 (2004)] entre taux de pertes à trois corps et nombre de triplets de fermions proches, le présent travail ajoute une dérivation, exprime le facteur de proportionnalité en termes du paramètre à trois corps, et inclut le cas général où il y a plusieurs contributions au contact à trois corps et plusieurs paramètres à trois corps.

<footnote><sup>a</sup> Laboratoire Kastler Brossel, Ecole Normale Supérieure - Université PSL, CNRS, Collège de France, Sorbonne Université, 75005 Paris, France</footnote>

<footnote><sup>b</sup> Laboratoire de Physique de l'Ecole Normale Supérieure, ENS - Université PSL, Sorbonne Université, Université Paris Cité, CNRS, 75005 Paris, France</footnote>

<footnote><sup>E-mails:</sup> werner@lkb.ens.fr (F. Werner), xavier.leyronas@phys.ens.fr (X. Leyronas)</footnote>
Research article / *Article de recherche*

# Three-body contact for fermions.

## I. General relations

*Contact à trois corps pour les fermions.*

### I. Relations générales

Félix Werner<sup>●,a</sup> and Xavier Leyronas<sup>●,b</sup>

**Abstract.** 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 ↑ and ↓ 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.

**Résumé.** Nous considérons le gaz de Fermi résonnant, à savoir des fermions avec deux états internes à trois dimensions avec des interactions à courte portée de grande longueur de diffusion. Nous introduisons une quantité, le contact à trois corps, qui détermine plusieurs observables. Pour le modèle de portée nulle, le nombre de triplets de fermions proches, la queue de la distribution selon l'impulsion du centre de masse des paires de fermions proches, ainsi que la queue de la distribution en impulsion à deux particules, sont exprimées en termes du contact à trois corps. Pour une portée non nulle, le taux de formation de dimères fortement liés par recombinaison à trois corps, ainsi que la contribution à trois corps à la correction de portée finie à l'énergie, sont exprimées en termes du contact à trois corps et d'un paramètre à trois corps. Ce paramètre à trois corps, qui tend vers zéro dans la limite de portée nulle, est défini *via* le comportement asymptotique de l'état de diffusion d'énergie nulle à des distances intermédiaires entre la portée et la longueur de diffusion à deux corps. En général, le contact à trois corps a différentes contributions repérées par des indices de spin et de moment cinétique, et le paramètre à trois corps peut dépendre de ces indices. Nous incluons aussi la généralisation à des masses différentes pour les particules ↑ et ↓. Par rapport à la relation donnée dans [Petrov, Salomon et Shlyapnikov, PRL 93, 090404 (2004)] entre taux de pertes à trois corps et nombre de triplets de fermions proches, le présent travail ajoute une dérivation, exprime le facteur de proportionnalité en termes du paramètre à trois corps, et inclut le cas général où il y a plusieurs contributions au contact à trois corps et plusieurs paramètres à trois corps.

<footnote><sup>a</sup> Laboratoire Kastler Brossel, Ecole Normale Supérieure - Université PSL, CNRS, Collège de France, Sorbonne Université, 75005 Paris, France</footnote>

<footnote><sup>b</sup> Laboratoire de Physique de l'Ecole Normale Supérieure, ENS - Université PSL, Sorbonne Université, Université Paris Cité, CNRS, 75005 Paris, France</footnote>

<footnote><sup>E-mails:</sup> werner@lkb.ens.fr (F. Werner), xavier.leyronas@phys.ens.fr (X. Leyronas)</footnote>
**Keywords.** Unitary gas, fermions, cold atoms.

**Mots-clés.** Gaz unitaire, fermions, atomes froids.

**Funding.** F. W. acknowledges the hospitality of the Aspen Center for Physics, and support from ERC (project Critisup2, H2020 Adv-743159) and ANR (project LODIS, ANR-21-CE30-0033).

Manuscript received 17 November 2022, revised 7 December 2023, accepted 26 February 2024.

# 1. 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 $\uparrow$ and $\downarrow$) and an interaction of vanishing range characterized by its $s$-wave scattering length $a_2$. When $1/a_2$ changes from $-\infty$ to $+\infty$, 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–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. [4–34].

From a theoretical viewpoint, this resonant Fermi gas is a difficult problem. As for most strongly correlated many-body problems in dimension $> 1$, numerical methods are generally the only option to make precise predictions. Furthermore the zero-range nature of the interactions typically constitutes an additional difficulty for numerical computations. But zero-range interactions also give rise to specific exact relations, called Tan relations, involving a ubiquitous quantity, the two-body contact $C_2$ [35–56]. In particular, $C_2$ determines the probability to find two particles close to each other. For the resonant Fermi gas, $C_2$ was measured and computed in numerous studies, see e.g. [15, 25–34, 49] and [28, 57–67] respectively.

Similarly, universal relations involving a two-body contact $C_2$ hold for two-component fermions in 2D [46, 51–53, 55, 56, 68, 69] and 1D [46, 70], and for single-component bosons in 1D [46, 71] and 2D [46, 72]. For a 2D Bose gas, $C_2$ was recently measured interferometrically [73]. Several two-body contacts appear in the general relations for single-component fermions with $p$-wave [74–78] or higher partial-wave [79] short-range interactions in 3D, and dipolar plus short-range interactions in 2D [80].

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–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 Sec. 2, and we consider models with a small finite interaction range in Sec. 3. Within the zero-range model, the number of triplets
**Keywords.** Unitary gas, fermions, cold atoms.

**Mots-clés.** Gaz unitaire, fermions, atomes froids.

**Funding.** F. W. acknowledges the hospitality of the Aspen Center for Physics, and support from ERC (project Critisup2, H2020 Adv-743159) and ANR (project LODIS, ANR-21-CE30-0033).

Manuscript received 17 November 2022, revised 7 December 2023, accepted 26 February 2024.

# 1. 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 $\uparrow$ and $\downarrow$) and an interaction of vanishing range characterized by its $s$-wave scattering length $a_2$. When $1/a_2$ changes from $-\infty$ to $+\infty$, 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–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. [4–34].

From a theoretical viewpoint, this resonant Fermi gas is a difficult problem. As for most strongly correlated many-body problems in dimension $> 1$, numerical methods are generally the only option to make precise predictions. Furthermore the zero-range nature of the interactions typically constitutes an additional difficulty for numerical computations. But zero-range interactions also give rise to specific exact relations, called Tan relations, involving a ubiquitous quantity, the two-body contact $C_2$ [35–56]. In particular, $C_2$ determines the probability to find two particles close to each other. For the resonant Fermi gas, $C_2$ was measured and computed in numerous studies, see e.g. [15, 25–34, 49] and [28, 57–67] respectively.

Similarly, universal relations involving a two-body contact $C_2$ hold for two-component fermions in 2D [46, 51–53, 55, 56, 68, 69] and 1D [46, 70], and for single-component bosons in 1D [46, 71] and 2D [46, 72]. For a 2D Bose gas, $C_2$ was recently measured interferometrically [73]. Several two-body contacts appear in the general relations for single-component fermions with $p$-wave [74–78] or higher partial-wave [79] short-range interactions in 3D, and dipolar plus short-range interactions in 2D [80].

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–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 Sec. 2, and we consider models with a small finite interaction range in Sec. 3. Within the zero-range model, the number of triplets
of particles separated by a small distance (Sec. 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 many-body 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 $\Gamma_3$ (Sec. 3.1), and the three-body contribution to the energy correction induced by the finite interaction range $\delta E_3$ (Sec. 3.2). We express $\Gamma_3$ and $\delta 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 $\gg b$ and $\ll |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 App. D. We consider a stationary state throughout the article, and discuss statistical mixtures and non-stationary states in App. 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.

## 2. 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.<fnref n="1" />

### 2.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 $\epsilon$ is given by

$$N_2(\epsilon) \sim_{\epsilon \to 0} C_2 \frac{\epsilon}{4\pi} \qquad (1)$$

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

$$R = \sqrt{\frac{2}{3}(r_{12}^2 + r_{13}^2 + r_{23}^2)} \qquad (2)$$

where $r_{ij}$ 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 < \epsilon$ is given by

$$N_3(\epsilon) \sim_{\epsilon \to 0} C_3 \epsilon^{2s+2} \qquad (3)$$

<footnote n="1">¹The absence of $(N_↑ + N_↓)$-body Efimov effect (in the equal-mass case) was shown for $(N_↑, N_↓) = (2,1)$ [96, 97], $(2,2)$ [98–101], $(3,1)$ [102], $(4,1)$ [103], $(5,1)$ [104], and $(N,1)$ for any $N$ [105]. The zero-range model was proven to be self-adjoint in the $(2,2)$ [101] and $(N,1)$ [105] cases. In the latter case, some rigorous results about the Tan relations were also obtained [105]. The convergence of finite-range models towards the zero-range model in the zero-range limit was confirmed by various theoretical studies, see e.g. [106, 107] for single-channel models and [48–50, 108, 109] for two-channel models, and by numerous theory-experiment comparisons for the many-body problem, e.g. [12–16, 28, 32, 33, 66, 67, 110–114].</footnote>
of particles separated by a small distance (Sec. 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 many-body 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 $\Gamma_3$ (Sec. 3.1), and the three-body contribution to the energy correction induced by the finite interaction range $\delta E_3$ (Sec. 3.2). We express $\Gamma_3$ and $\delta 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 $\gg b$ and $\ll |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 App. D. We consider a stationary state throughout the article, and discuss statistical mixtures and non-stationary states in App. 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.

## 2. 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.<fnref n="1" />

### 2.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 $\epsilon$ is given by

$$N_2(\epsilon) \sim_{\epsilon \to 0} C_2 \frac{\epsilon}{4\pi} \quad (1)$$

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

$$R = \sqrt{\frac{2}{3} (r_{12}^2 + r_{13}^2 + r_{23}^2)} \quad (2)$$

where $r_{ij}$ 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 < \epsilon$ is given by

$$N_3(\epsilon) \sim_{\epsilon \to 0} C_3 \epsilon^{2s+2} \quad (3)$$

<footnote n="1">¹The absence of $(N_↑ + N_↓)$-body Efimov effect (in the equal-mass case) was shown for $(N_↑, N_↓) = (2,1)$ [96, 97], $(2,2)$ [98–101], $(3,1)$ [102], $(4,1)$ [103], $(5,1)$ [104], and $(N,1)$ for any $N$ [105]. The zero-range model was proven to be self-adjoint in the $(2,2)$ [101] and $(N,1)$ [105] cases. In the latter case, some rigorous results about the Tan relations were also obtained [105]. The convergence of finite-range models towards the zero-range model in the zero-range limit was confirmed by various theoretical studies, see e.g. [106, 107] for single-channel models and [48–50, 108, 109] for two-channel models, and by numerous theory-experiment comparisons for the many-body problem, e.g. [12–16, 28, 32, 33, 66, 67, 110–114].</footnote>
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

$$ (1 - s^2) \sin\left(\frac{s\pi}{2}\right) - \frac{4}{\sqrt{3}} s \cos\left(\frac{s\pi}{6}\right) + 4 \sin\left(\frac{s\pi}{6}\right) = 0. \quad (4) $$

The scaling $N_3(\epsilon) \propto \epsilon^{2s+2}$ was already obtained in [98, 115] (see Sec. 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–120].

In general there are two different contributions to the three-body contact, coming from $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ spin configurations: Denoting by $N_{2,1}$ (resp. $N_{1,2}$) the contributions to $N_3(\epsilon)$ from triplets of particles of spins $\uparrow\uparrow$ (resp. $\uparrow\downarrow\downarrow$), we have

$$ N_{2,1}(\epsilon) \sim_{\epsilon \to 0} C_{2,1} \epsilon^{2s+2} \quad (5) $$

$$ N_{1,2}(\epsilon) \sim_{\epsilon \to 0} C_{1,2} \epsilon^{2s+2} \quad (6) $$

where $C_{2,1}$ (resp. $C_{1,2}$) is what we call the $\uparrow\uparrow\downarrow$ (resp. $\uparrow\downarrow\downarrow$) three-body contact. Clearly,

$$ C_3 = C_{2,1} + C_{1,2}. \quad (7) $$

### Remarks 1.

* 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 ($\uparrow\uparrow$ or $\downarrow\downarrow$) is negligible in the $\epsilon \to 0$ limit (it scales as $\epsilon^5$), and $N_2(\epsilon)$ is dominated by the contribution from pairs of particles with opposite spins ($\uparrow\downarrow$).

Similarly, the contribution to Eq. (3) coming from triplets of particles with identical spins ($\uparrow\uparrow\uparrow$ or $\downarrow\downarrow\downarrow$) is negligible in the $\epsilon \to 0$ limit (it scales as $\epsilon^{10}$), and $N_3(\epsilon)$ is dominated by the contributions from triplets of particles with non-identical spins, $\uparrow\uparrow\downarrow$ or $\uparrow\downarrow\downarrow$, in agreement with Eq. (7).

* For comparison, in the non-interacting case, the number of nearby pairs and triplets scales as

$$ N_2^{(0)}(\epsilon) \propto \epsilon^3 \quad (8) $$

$$ N_3^{(0)}(\epsilon) \propto \epsilon^8 \quad (9) $$

(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).

Eq. (8) [resp. Eq. (9)] is dominated by the contribution from pairs (resp. triplets) of particles with non-identical spins.

Eq. (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 $\epsilon^2$ suppression factor compared to the completely uncorrelated case of non-interacting distinguishable particles

$$ N_{3, \text{ distinguishable}}^{(0)}(\epsilon) \propto \epsilon^6. \quad (10) $$

* 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,
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

$$ (1 - s^2) \sin\left(\frac{s\pi}{2}\right) - \frac{4}{\sqrt{3}} s \cos\left(\frac{s\pi}{6}\right) + 4 \sin\left(\frac{s\pi}{6}\right) = 0. \quad (4) $$

The scaling $N_3(\epsilon) \propto \epsilon^{2s+2}$ was already obtained in [98, 115] (see Sec. 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–120].

In general there are two different contributions to the three-body contact, coming from $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ spin configurations: Denoting by $N_{2,1}$ (resp. $N_{1,2}$) the contributions to $N_3(\epsilon)$ from triplets of particles of spins $\uparrow\uparrow$ (resp. $\uparrow\downarrow\downarrow$), we have

$$ N_{2,1}(\epsilon) \sim_{\epsilon \to 0} C_{2,1} \epsilon^{2s+2} \quad (5) $$

$$ N_{1,2}(\epsilon) \sim_{\epsilon \to 0} C_{1,2} \epsilon^{2s+2} \quad (6) $$

where $C_{2,1}$ (resp. $C_{1,2}$) is what we call the $\uparrow\uparrow\downarrow$ (resp. $\uparrow\downarrow\downarrow$) three-body contact. Clearly,

$$ C_3 = C_{2,1} + C_{1,2}. \quad (7) $$

## Remarks 1.

* 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 ($\uparrow\uparrow$ or $\downarrow\downarrow$) is negligible in the $\epsilon \to 0$ limit (it scales as $\epsilon^5$), and $N_2(\epsilon)$ is dominated by the contribution from pairs of particles with opposite spins ($\uparrow\downarrow$).

Similarly, the contribution to Eq. (3) coming from triplets of particles with identical spins ($\uparrow\uparrow\uparrow$ or $\downarrow\downarrow\downarrow$) is negligible in the $\epsilon \to 0$ limit (it scales as $\epsilon^{10}$), and $N_3(\epsilon)$ is dominated by the contributions from triplets of particles with non-identical spins, $\uparrow\uparrow\downarrow$ or $\uparrow\downarrow\downarrow$, in agreement with Eq. (7).

* For comparison, in the non-interacting case, the number of nearby pairs and triplets scales as

$$ N_2^{(0)}(\epsilon) \propto \epsilon^3 \quad (8) $$

$$ N_3^{(0)}(\epsilon) \propto \epsilon^8 \quad (9) $$

(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).

Eq. (8) [resp. Eq. (9)] is dominated by the contribution from pairs (resp. triplets) of particles with non-identical spins.

Eq. (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 $\epsilon^2$ suppression factor compared to the completely uncorrelated case of non-interacting distinguishable particles

$$ N_{3, \text{ distinguishable}}^{(0)}(\epsilon) \propto \epsilon^6. \quad (10) $$

* 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, $\psi \propto 1/r$, which yields Eq. (1), while in the limit of vanishing hyperradius $R$ between three particles, $\psi \propto R^{s-2}$, which yields Eq. (3). Note that since $s < 2$, the wavefunction indeed diverges for $R \to 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 zero-range interactions overcompensates the antibunching effect due to Pauli exclusion.<fnref n="2" />

## 2.2. Density correlation functions

The probability density of finding a spin-↑ particle at $\mathbf{r}_1$ and a spin-↓ particle at $\mathbf{r}_2$ is given by the pair correlation function $g_2(\mathbf{r}_1, \mathbf{r}_2) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle = \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\downarrow(\mathbf{r}_2) \rangle$. Similarly, the probability density of finding a spin-↑ particle at $\mathbf{r}_1$, a second spin-↑ particle at $\mathbf{r}_2$, and a spin-↓ particle at $\mathbf{r}_3$ is given by triplet correlation function

$$g_{2,1}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\uparrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_3) \hat{\psi}_\uparrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle \\ = \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\uparrow(\mathbf{r}_2) \hat{n}_\downarrow(\mathbf{r}_3) \rangle - \delta^3(\mathbf{r}_1 - \mathbf{r}_2) \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\downarrow(\mathbf{r}_3) \rangle.$$

The second-order density correlation function has the short-distance asymptotic behavior [35]

$$\int d^3c \ g_2(\mathbf{c} + \frac{\mathbf{r}}{2}, \mathbf{c} - \frac{\mathbf{r}}{2}) \underset{r \to 0}{\sim} \frac{1}{(4\pi)^2} \frac{C_2}{r^2}.$$

For the third-order density correlation function, we find the short-distance asymptotic behavior

$$\int d^3C \ d^5\Omega \ g_{2,1}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) \underset{R \to 0}{\sim} C_{2,1} R^{2s-4} \frac{32(s+1)}{3\sqrt{3}} \quad (11)$$

where the limit $R \to 0$ is taken for fixed $C$ and $\Omega$. Here a change of coordinates is implied between $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ and $(C, R, \Omega)$, with $R$ the hyperradius defined in Eq. (2), $C = (\mathbf{r}_1 + \mathbf{r}_2 + \mathbf{r}_3)/3$ the center-of-mass, and $\Omega$ the hyperangles which are five dimensionless coordinates that remain to determine the positions of the three particles, see Eq. (91). Similarly, the ↑↑↑ triplet correlation function $g_{1,2}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle$ satisfies

$$\int d^3C \ d^5\Omega \ g_{1,2}(\mathbf{r}_3, \mathbf{r}_1, \mathbf{r}_2) \underset{R \to 0}{\sim} C_{1,2} R^{2s-4} \frac{32(s+1)}{3\sqrt{3}}. \quad (12)$$

## 2.3. 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 $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ be the (orbital) many-body wavefunction. Without loss of generality, we can assume $N_\uparrow \ge 2$, and consider that particles $(1, 2, 3)$ have spins ($\uparrow, \downarrow, \downarrow$). We take this convention throughout the article. This means that $\psi$ is antisymmetric w.r.t. exchange of $\mathbf{r}_1$ and $\mathbf{r}_2$ ($\psi$ is also antisymmetric w.r.t. exchange of any other pair of same-spin particles). We take the normalization $\int d^3\mathbf{r}_1 \dots d^3\mathbf{r}_N |\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)|^2 = 1$.

Within the zero-range model, the stationary Schrödinger equation

$$\sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \psi = E \psi \quad (13)$$

<footnote n="2">²On the other hand, the Pauli exclusion effect is not overcompensated by too much: Due to the repulsive effective three-body potential $s^2/R^2$ (see App. A), we still have $N_3(\epsilon)/\epsilon^2 \propto \epsilon^{2s} \underset{\epsilon \to 0}{\to} 0$, which implies that the three-body loss rate divided by the thermalization rate vanishes in the zero-density limit, a crucial ingredient for the zero-range model to be an accurate description of ultracold-atom experiments (see Sec. 3.1).</footnote>
which causes the wavefunction to diverge: When the distance $r$ between two opposite-spin particles vanishes, $\psi \propto 1/r$, which yields Eq. (1), while in the limit of vanishing hyperradius $R$ between three particles, $\psi \propto R^{s-2}$, which yields Eq. (3). Note that since $s < 2$, the wavefunction indeed diverges for $R \to 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 zero-range interactions overcompensates the antibunching effect due to Pauli exclusion.<fnref n="2" />

## 2.2. Density correlation functions

The probability density of finding a spin-↑ particle at $\mathbf{r}_1$ and a spin-↓ particle at $\mathbf{r}_2$ is given by the pair correlation function $g_2(\mathbf{r}_1, \mathbf{r}_2) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle = \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\downarrow(\mathbf{r}_2) \rangle$. Similarly, the probability density of finding a spin-↑ particle at $\mathbf{r}_1$, a second spin-↑ particle at $\mathbf{r}_2$, and a spin-↓ particle at $\mathbf{r}_3$ is given by triplet correlation function

$$g_{2,1}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_3) \hat{\psi}_\uparrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle \\ = \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\downarrow(\mathbf{r}_2) \hat{n}_\downarrow(\mathbf{r}_3) \rangle - \delta^3(\mathbf{r}_1 - \mathbf{r}_2) \langle \hat{n}_\uparrow(\mathbf{r}_1) \hat{n}_\downarrow(\mathbf{r}_3) \rangle.$$

The second-order density correlation function has the short-distance asymptotic behavior [35]

$$\int d^3c \ g_2(\mathbf{c} + \frac{\mathbf{r}}{2}, \mathbf{c} - \frac{\mathbf{r}}{2}) \underset{r \to 0}{\sim} \frac{1}{(4\pi)^2} \frac{C_2}{r^2}.$$

For the third-order density correlation function, we find the short-distance asymptotic behavior

$$\int d^3C \ d^5\Omega \ g_{2,1}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) \underset{R \to 0}{\sim} C_{2,1} R^{2s-4} \frac{32(s+1)}{3\sqrt{3}} \quad (11)$$

where the limit $R \to 0$ is taken for fixed $C$ and $\Omega$. Here a change of coordinates is implied between $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ and $(C, R, \Omega)$, with $R$ the hyperradius defined in Eq. (2), $C = (\mathbf{r}_1 + \mathbf{r}_2 + \mathbf{r}_3)/3$ the center-of-mass, and $\Omega$ the hyperangles which are five dimensionless coordinates that remain to determine the positions of the three particles, see Eq. (91). Similarly, the ↑↓↓ triplet correlation function $g_{1,2}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_2) \hat{\psi}_\downarrow^\dagger(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_3) \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle$ satisfies

$$\int d^3C \ d^5\Omega \ g_{1,2}(\mathbf{r}_3, \mathbf{r}_1, \mathbf{r}_2) \underset{R \to 0}{\sim} C_{1,2} R^{2s-4} \frac{32(s+1)}{3\sqrt{3}}. \quad (12)$$

## 2.3. 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 $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ be the (orbital) many-body wavefunction. Without loss of generality, we can assume $N_\uparrow \ge 2$, and consider that particles $(1, 2, 3)$ have spins ($\uparrow, \downarrow, \downarrow$). We take this convention throughout the article. This means that $\psi$ is antisymmetric w.r.t. exchange of $\mathbf{r}_1$ and $\mathbf{r}_2$ ($\psi$ is also antisymmetric w.r.t. exchange of any other pair of same-spin particles). We take the normalization $\int d^3\mathbf{r}_1 \dots d^3\mathbf{r}_N |\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)|^2 = 1$. Within the zero-range model, the stationary Schrödinger equation

$$\sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \psi = E \psi \quad (13)$$

<footnote n="2">²On the other hand, the Pauli exclusion effect is not overcompensated by too much: Due to the repulsive effective three-body potential $s^2/R^2$ (see App. A), we still have $N_3(\epsilon)/\epsilon^2 \propto \epsilon^{2s} \underset{\epsilon \to 0}{\to} 0$, which implies that the three-body loss rate divided by the thermalization rate vanishes in the zero-density limit, a crucial ingredient for the zero-range model to be an accurate description of ultracold-atom experiments (see Sec. 3.1).</footnote>
contains an external trapping potential $U$ but no interaction potential; instead, $\psi$ should satisfy a contact condition in the limit where two particles of opposite spin approach each other: There exists $A$ such that

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) = \int_{r \to 0} \left( \frac{1}{r} - \frac{1}{a_2} \right) A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) + O(r) \quad (14) $$

where $a_2$ is the two-body scattering length, $r = \|\mathbf{r}_1 - \mathbf{r}_3\|$ is the distance between the opposite-spin particles 1 and 3, and $\mathbf{c} = (\mathbf{r}_1 + \mathbf{r}_3)/2$ is their center-of-mass. The limit $r \to 0$ is taken for fixed $\mathbf{c}$ and fixed positions of the remaining particles $(\mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{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 $\psi$ and energies $E$ of the zero-range model.<fnref n="3" />

When particles 1, 2 and 3 approach each other, the wavefunction of any stationary state has the asymptotic behavior [38, 96, 98, 115]

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \sim R^{s-2} \sum_{m=-1}^{+1} \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (15) $$

Here, as in Eq. (11), $R$ is the hyperradius of particles (1, 2, 3) defined in Eq. (2), $\mathbf{C}$ is their center-of-mass, and $\Omega$ denotes their hyperangles. The limit $R \to 0$ is taken for fixed $\Omega, \mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$. The unitary hyperangular wavefunctions $\phi_m(\Omega)$ are such that $R^{s-2}\phi_m(\Omega)$ 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 \in \{-1, 0, 1\}$, see App. A for more details.<fnref n="4" />,<fnref n="5" />

It is known [35] that $C_2$ is given by the norm of the function $A$ that appears in Eq. (14),

$$ C_2 = (4\pi)^2 N_\uparrow N_\downarrow \int |A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3c d^3r_2 d^3r_4 \dots d^3r_N. \quad (16) $$

Similarly, the three-body contact is given by the norm of the function $B$ that appears in Eq. (15),

$$ C_{2,1} = N_\uparrow (N_\downarrow - 1) N_\downarrow \frac{3\sqrt{3}}{32(s+1)} \sum_{m=-1}^{+1} \int |B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3C d^3r_4 \dots d^3r_N. \quad (17) $$

The expression (5) for the number of nearby $\uparrow\uparrow\downarrow\downarrow$ fermion triplets, together with Eq. (17), simply follow from Eq. (15) by integrating $|\psi|^2$ over the $R < \epsilon$ region.<fnref n="6" /> Similarly, 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}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = N_\uparrow (N_\downarrow - 1) N_\downarrow \int d^3r_4 \dots d^3r_N |\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)|^2$.

There is a completely analogous relation between $C_{1,2}$ and the behavior of the many-body wavefunction when three particles of spins $\uparrow\downarrow\downarrow\downarrow$ approach each other (provided $N_\uparrow \ge 2$). Specifically, considering that particle 4 has spin $\downarrow$ (while particles 1, 3 still have spins $\uparrow, \downarrow$) and denoting by $\tilde{R}, \tilde{\Omega}, \tilde{C}$ the hyperradius, hyperangles and center-of-mass associated to particles 3, 4, 1 [obtained by replacing $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ with $(\mathbf{r}_3, \mathbf{r}_4, \mathbf{r}_1)$ in Eqs. (89), (90), (91)] we have

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \sim \tilde{R}^{s-2} \sum_{\tilde{R} \to 0} \phi_m(\tilde{\Omega}) \tilde{B}_m(\tilde{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N) \quad (18) $$

<footnote n="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].</footnote>

<footnote n="4">⁴There 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.</footnote>

<footnote n="5">⁵An asymptotic behavior similar to (15) holds when any three particles with spins $\uparrow\downarrow\downarrow$ approach each other, the functions $B_m$ corresponding to different triplets of particles being simply related to each other by antisymmetry.</footnote>

<footnote n="6">⁶Here we used the change of integration variables $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) \to (\mathbf{r}, \mathbf{\rho}, \mathbf{C})$ with $\mathbf{\rho}$ defined in (89), of Jacobian $|\partial(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)/\partial(\mathbf{r}, \mathbf{\rho}, \mathbf{C})| = (\sqrt{3}/2)^3$. We also used the property $\int d^5\Omega \phi_m(\Omega)^* \phi_{m'}(\Omega) = \delta_{m,m'}$.</footnote>
contains an external trapping potential $U$ but no interaction potential; instead, $\psi$ should satisfy a contact condition in the limit where two particles of opposite spin approach each other: There exists $A$ such that

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) = \int_{r \to 0} \left( \frac{1}{r} - \frac{1}{a_2} \right) A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) + O(r) \quad (14) $$

where $a_2$ is the two-body scattering length, $r = \|\mathbf{r}_1 - \mathbf{r}_3\|$ is the distance between the opposite-spin particles 1 and 3, and $\mathbf{c} = (\mathbf{r}_1 + \mathbf{r}_3)/2$ is their center-of-mass. The limit $r \to 0$ is taken for fixed $\mathbf{c}$ and fixed positions of the remaining particles $(\mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{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 $\psi$ and energies $E$ of the zero-range model.<fnref n="3" />

When particles 1, 2 and 3 approach each other, the wavefunction of any stationary state has the asymptotic behavior [38, 96, 98, 115]

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \sim R^{s-2} \sum_{m=-1}^{+1} \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (15) $$

Here, as in Eq. (11), $R$ is the hyperradius of particles (1, 2, 3) defined in Eq. (2), $\mathbf{C}$ is their center-of-mass, and $\Omega$ denotes their hyperangles. The limit $R \to 0$ is taken for fixed $\Omega, \mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$. The unitary hyperangular wavefunctions $\phi_m(\Omega)$ are such that $R^{s-2}\phi_m(\Omega)$ 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 \in \{-1, 0, 1\}$, see App. A for more details.<fnref n="4" />,<fnref n="5" />

It is known [35] that $C_2$ is given by the norm of the function $A$ that appears in Eq. (14),

$$ C_2 = (4\pi)^2 N_\uparrow N_\downarrow \int |A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3\mathbf{c} d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N. \quad (16) $$

Similarly, the three-body contact is given by the norm of the function $B$ that appears in Eq. (15),

$$ C_{2,1} = N_\uparrow (N_\downarrow - 1) N_\downarrow \frac{3\sqrt{3}}{32(s+1)} \sum_{m=-1}^{+1} \int |B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3\mathbf{C} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N. \quad (17) $$

The expression (5) for the number of nearby $\uparrow\uparrow\downarrow\downarrow$ fermion triplets, together with Eq. (17), simply follow from Eq. (15) by integrating $|\psi|^2$ over the $R < \epsilon$ region.<fnref n="6" /> Similarly, 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}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = N_\uparrow (N_\downarrow - 1) N_\downarrow \int d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N |\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)|^2$.

There is a completely analogous relation between $C_{1,2}$ and the behavior of the many-body wavefunction when three particles of spins $\uparrow\downarrow\downarrow\downarrow$ approach each other (provided $N_\uparrow \ge 2$). Specifically, considering that particle 4 has spin $\downarrow$ (while particles 1, 3 still have spins $\uparrow, \downarrow$) and denoting by $\tilde{R}, \tilde{\Omega}, \tilde{C}$ the hyperradius, hyperangles and center-of-mass associated to particles 3, 4, 1 [obtained by replacing $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ with $(\mathbf{r}_3, \mathbf{r}_4, \mathbf{r}_1)$ in Eqs. (89), (90), (91)] we have

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \sim \tilde{R}^{s-2} \sum_{\tilde{R} \to 0} \phi_m(\tilde{\Omega}) \tilde{B}_m(\tilde{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N) \quad (18) $$

<footnote n="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].</footnote>

<footnote n="4">⁴There 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.</footnote>

<footnote n="5">⁵An asymptotic behavior similar to (15) holds when any three particles with spins $\uparrow\downarrow\downarrow$ approach each other, the functions $B_m$ corresponding to different triplets of particles being simply related to each other by antisymmetry.</footnote>

<footnote n="6">⁶Here we used the change of integration variables $(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) \to (\mathbf{r}, \mathbf{\rho}, \mathbf{C})$ with $\mathbf{\rho}$ defined in (89), of Jacobian $|\partial(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)/\partial(\mathbf{r}, \mathbf{\rho}, \mathbf{C})| = (\sqrt{3}/2)^3$. We also used the property $\int d^5\Omega \phi_m(\Omega)^* \phi_{m'}(\Omega) = \delta_{m,m'}$.</footnote>
which yields the relations for $N_{1,2}$ and $g_{1,2}$, Eqs. (6), (12), with

$$C_{1,2} = N_{\downarrow}(N_{\downarrow}-1)N_{\uparrow} \frac{3\sqrt{3}}{32(s+1)} \sum_{m=-1}^{+1} \int |\tilde{B}_{m}(\mathbf{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N)|^2 d^3C d^3r_2 d^3r_5 \dots d^3r_N. \quad (19)$$

Here we assumed $N_{\downarrow} \ge 2$; if $N_{\downarrow} = 1$, then $N_{1,2}$ is obviously zero, and $C_{1,2} = 0$.

**Remark 2.** Higher-body contacts can be defined in the same way than the three-body contacts. When $j_{\uparrow}$ particles of spin $\uparrow$ and $j_{\downarrow}$ particles of spin $\downarrow$ approach each other, the N-body wavefunction factorizes into the product of (i) a function of the relative positions of the $j = j_{\uparrow} + j_{\downarrow}$ nearby particles, given by a zero-energy solution of the $j_{\uparrow} + j_{\downarrow}$ body problem at $a_2=\infty$, 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_{\uparrow} + j_{\downarrow}$ body contact $C_{j_{\uparrow}, j_{\downarrow}}$ (up to a prefactor which is a matter of definition).

## 2.4. 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_{\sigma}(\mathbf{k}) \underset{k \to \infty}{\sim} \frac{C_2}{k^4} \quad (20)$$

with the normalization $\int N_{\sigma}(\mathbf{k}) d^3\mathbf{k}/(2\pi)^3 = N_{\sigma}$ (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 $\mathbf{K}$ (this is allowed since the corresponding operators commute). Let $N_2(\epsilon, \mathbf{K})$ be the probability distribution over $\mathbf{K}$ conditional to $r < \epsilon$, with the normalization $\int N_2(\epsilon, \mathbf{K}) d^3\mathbf{K}/(2\pi)^3 = N_2(\epsilon)$. In other words, $N_2(\epsilon, \mathbf{K})$ is the center-of-mass momentum distribution of the pairs of particles separated by a distance $< \epsilon$, normalized to the total number of such pairs.<fnref n="7" /> We have

$$N_2(\epsilon, \mathbf{K}) \underset{\epsilon \to 0}{\sim} N_P(\mathbf{K}) \frac{\epsilon}{4\pi} \quad (21)$$

where $N_P(\mathbf{K})$ is what we call the center-of-mass momentum distribution of nearby fermion pairs. With this definition, one simply has the normalization

$$\int N_P(\mathbf{K}) \frac{d^3\mathbf{K}}{(2\pi)^3} = C_2 \quad (22)$$

as a consequence of (1).⁸ The tail of $N_P(\mathbf{K})$ is determined by the three-body contact⁹:

$$\tilde{N}_P(\mathbf{K}) \underset{\mathbf{K} \to \infty}{\sim} \mathcal{M} \frac{C_3}{\mathbf{K}^{2s+4}} \quad (23)$$

$$\sum_{i:\uparrow, j:\downarrow} \theta(\epsilon - \hat{r}_{ij}) (2\pi)^3 \delta^3(\hat{\mathbf{k}}_{ij} - \mathbf{K}),$$

where $\theta$ is the Heaviside function, while $\hat{r}_{ij} = \|\hat{\mathbf{r}}_j - \hat{\mathbf{r}}_i\|$ and $\hat{\mathbf{k}}_{ij} = \hat{\mathbf{k}}_i + \hat{\mathbf{k}}_j$ are the operators corresponding to the relative distance and the center-of-mass momentum of particles i and j.

⁸$N_P(\mathbf{K})$ appears naturally in the diagrammatic formalism: In a homogeneous system, $N_P(\mathbf{K})$ divided by the volume equals $-(m^2/\hbar^2)\Gamma(\mathbf{K}, \tau = 0^{-})$ where $\Gamma$ 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 $\phi(\mathbf{r})\phi(\mathbf{r}')$ with $\phi$ the zero-energy two-body scattering state, and setting $\mathbf{r} = \mathbf{r}' = 0$.

⁹The fact that $\Gamma(\mathbf{K}, \tau = 0^{-})$ has a tail $\propto C_3/K^{2s+4}$ was pointed out to us by Shina Tan (private communication, Aspen, 2011).

<footnote n="7">⁷More formally, $N_2(\epsilon, \mathbf{K})$ is the expectation value of the operator</footnote>
which yields the relations for $N_{1,2}$ and $g_{1,2}$, Eqs. (6), (12), with

$$C_{1,2} = N_{\uparrow}(N_{\downarrow}-1)N_{\uparrow} \frac{3\sqrt{3}}{32(s+1)} \sum_{m=-1}^{+1} \int |\tilde{B}_{m}(\mathbf{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N)|^2 d^3C d^3r_2 d^3r_5 \dots d^3r_N. \quad (19)$$

Here we assumed $N_{\downarrow} \ge 2$; if $N_{\downarrow} = 1$, then $N_{1,2}$ is obviously zero, and $C_{1,2} = 0$.

**Remark 2.** Higher-body contacts can be defined in the same way than the three-body contacts. When $j_\uparrow$ particles of spin $\uparrow$ and $j_\downarrow$ particles of spin $\downarrow$ approach each other, the N-body wavefunction factorizes into the product of (i) a function of the relative positions of the $j = j_\uparrow + j_\downarrow$ nearby particles, given by a zero-energy solution of the $j_\uparrow + j_\downarrow$ body problem at $a_2=\infty$, 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_\uparrow + j_\downarrow$ body contact $C_{j_\uparrow, j_\downarrow}$ (up to a prefactor which is a matter of definition).

## 2.4. 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_{\sigma}(\mathbf{k}) \underset{k \to \infty}{\sim} \frac{C_2}{k^4} \quad (20)$$

with the normalization $\int N_{\sigma}(\mathbf{k}) d^3\mathbf{k} / (2\pi)^3 = N_{\sigma}$ (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 $\mathbf{K}$ (this is allowed since the corresponding operators commute). Let $N_2(\epsilon, \mathbf{K})$ be the probability distribution over $\mathbf{K}$ conditional to $r < \epsilon$, with the normalization $\int N_2(\epsilon, \mathbf{K}) d^3\mathbf{K} / (2\pi)^3 = N_2(\epsilon)$. In other words, $N_2(\epsilon, \mathbf{K})$ is the center-of-mass momentum distribution of the pairs of particles separated by a distance $< \epsilon$, normalized to the total number of such pairs.<fnref n="7" /> We have

$$N_2(\epsilon, \mathbf{K}) \underset{\epsilon \to 0}{\sim} N_P(\mathbf{K}) \frac{\epsilon}{4\pi} \quad (21)$$

where $N_P(\mathbf{K})$ is what we call the center-of-mass momentum distribution of nearby fermion pairs. With this definition, one simply has the normalization

$$\int N_P(\mathbf{K}) \frac{d^3\mathbf{K}}{(2\pi)^3} = C_2 \quad (22)$$

as a consequence of (1).⁸ The tail of $N_P(\mathbf{K})$ is determined by the three-body contact⁹:

$$\tilde{N}_P(\mathbf{K}) \underset{K \to \infty}{\sim} \mathcal{M} \frac{C_3}{K^{2s+4}} \quad (23)$$

$$\sum_{i:\uparrow, j:\downarrow} \theta(\epsilon - \hat{r}_{ij}) (2\pi)^3 \delta^3(\hat{\mathbf{k}}_{ij} - \mathbf{K}),$$

where $\theta$ is the Heaviside function, while $\hat{r}_{ij} = \|\hat{\mathbf{r}}_j - \hat{\mathbf{r}}_i\|$ and $\hat{\mathbf{k}}_{ij} = \hat{\mathbf{k}}_i + \hat{\mathbf{k}}_j$ are the operators corresponding to the relative distance and the center-of-mass momentum of particles i and j.

⁸$N_P(\mathbf{K})$ appears naturally in the diagrammatic formalism: In a homogeneous system, $N_P(\mathbf{K})$ divided by the volume equals $-(m^2/\hbar^2)\Gamma(\mathbf{K}, \tau = 0^{-})$ where $\Gamma$ 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 $\phi(\mathbf{r})\phi(\mathbf{r}')$ with $\phi$ the zero-energy two-body scattering state, and setting $\mathbf{r} = \mathbf{r}' = 0$.

⁹The fact that $\Gamma(\mathbf{K}, \tau = 0^{-})$ has a tail $\propto C_3/K^{2s+4}$ was pointed out to us by Shina Tan (private communication, Aspen, 2011).

<footnote n="7">⁷More formally, $N_2(\epsilon, \mathbf{K})$ is the expectation value of the operator</footnote>
with the prefactor

$$ \mathcal{M} = 32\pi^3 \frac{4^s}{3^{s+1/2}} (s+1) \Gamma(s+2)^2 \sin^2(s\pi) \mathcal{N}^2 \quad (24) $$

whose numerical value is

$$ \mathcal{M} \approx 2272. \quad (25) $$

Here $\tilde{N}_P(K)$ stands for the angular average $\int N_P(K) d\hat{K}/(4\pi)$ (with $\hat{K} := K/K$, and $d\hat{K}$ the differential solid angle, so that $d^3K = d\hat{K}K^2dK$).

To derive this result, we consider the two-body reduced density matrix $\rho_2(\mathbf{r}_\uparrow, \mathbf{r}_\downarrow; \mathbf{r}'_\uparrow, \mathbf{r}'_\downarrow) := \langle \tilde{\psi}_\uparrow^\dagger(\mathbf{r}'_\uparrow) \tilde{\psi}_\downarrow(\mathbf{r}'_\downarrow) \tilde{\psi}_\uparrow(\mathbf{r}_\uparrow) \tilde{\psi}_\downarrow(\mathbf{r}_\downarrow) \rangle$ or equivalently in first quantization

$$ \rho_2(\mathbf{r}_\uparrow, \mathbf{r}_\downarrow; \mathbf{r}'_\uparrow, \mathbf{r}'_\downarrow) = N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \psi^*(\mathbf{r}'_\uparrow, \mathbf{r}_2, \mathbf{r}'_\downarrow, \mathbf{r}_4, \dots, \mathbf{r}_N) \psi(\mathbf{r}_\uparrow, \mathbf{r}_2, \mathbf{r}'_\downarrow, \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (26) $$

Inserting the two-body contact condition (14) into (26) yields

$$ \rho_2 \left( \mathbf{c} - \frac{\mathbf{r}}{2}, \mathbf{c} + \frac{\mathbf{r}}{2}; \mathbf{c}' - \frac{\mathbf{r}'}{2}, \mathbf{c}' + \frac{\mathbf{r}'}{2} \right) \underset{r,r' \to 0}{\sim} \frac{g_P(\mathbf{c}, \mathbf{c}')}{(4\pi)^2 r r'} \quad (27) $$

where

$$ g_P(\mathbf{c}, \mathbf{c}') = (4\pi)^2 N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N A^*(\mathbf{c}'; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (28) $$

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,

$$ N_P(\mathbf{K}) = \int d^3\mathbf{c} d^3\mathbf{c}' e^{-i\mathbf{K} \cdot (\mathbf{c}-\mathbf{c}')} g_P(\mathbf{c}, \mathbf{c}') \quad (29) $$

as can be formally shown using the definition (21) of $N_P$. Hence

$$ N_P(\mathbf{K}) = (4\pi)^2 N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \left| \int d^3\mathbf{c} e^{-i\mathbf{K}\cdot\mathbf{c}} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \right|^2. \quad (30) $$

We then follow a reasoning resembling the one used to derive Eq. (20) in [56, Sec. IV A].

In the large $K$ limit, the Fourier transform with respect to $\mathbf{c}$ in (30) is dominated by the contributions from the singularities of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$, which occur when $\mathbf{c}$ approaches one of the $\mathbf{r}_i$ ($i=2,4,\dots, N$). This corresponds to particles 1, 3 and $i$ being close to each other [since $A(\mathbf{c}; \mathbf{r}_2, \dots)$ determines the wavefunction $\psi$ when particles 1 and 3 are close to $\mathbf{c}$, according to Eq. (14)]. For example, for $i=2$, the behavior of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$ in the limit $\mathbf{c} \to \mathbf{r}_2$ is determined by the asymptotic behavior of $\psi$ 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

$$ A(\mathbf{r}_2 - \mathbf{u}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \underset{u \to 0}{\sim} \mathcal{N} \frac{s}{2} \cos\left(\frac{s\pi}{2}\right) \left(\frac{2u}{\sqrt{3}}\right)^{s-1} \sum_{m=-1}^{1} Y_1^m(\hat{\mathbf{u}}) B_m(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N) $$

where we used the expression ((99), (100)) of $\phi_m(\Omega)$. There is a similar singularity of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$ when $\mathbf{c}$ approaches $\mathbf{r}_i$ with $4 \le i \le N$; when particle $i$ has spin $\downarrow$, the function $\tilde{B}_m$ introduced in (18) appears instead of $B_m$. This gives

$$ \int d^3\mathbf{c} e^{-i\mathbf{K}\cdot\mathbf{c}} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \underset{K \to \infty}{\sim} \frac{2^{s-2}}{3^{(s-1)/2}} \mathcal{N} \cos\left(\frac{s\pi}{2}\right) \sum_{m=-1}^{1} I_m(\mathbf{K}) \times \left[ \sum_{i:\uparrow, i \neq \uparrow} (1-2\delta_{i,2}) e^{-i\mathbf{K}\cdot\mathbf{r}_i} B_m(P_{2i}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)) + \sum_{i:\downarrow, i \neq \downarrow} (1-2\delta_{i,4}) e^{-i\mathbf{K}\cdot\mathbf{r}_i} \tilde{B}_m(P_{4i}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)) \right] \quad (31) $$

where the first (resp. second) sum over $i$ is taken over particles with spin $\uparrow$ (resp. $\downarrow$), $P_{ji}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)$ is obtained from $(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)$ by exchanging $\mathbf{r}_j$ with $\mathbf{r}_i$ ($P_{ii}$ is the identity),
with the prefactor

$$ \mathcal{M} = 32\pi^3 \frac{4^s}{3^{s+1/2}} (s+1) \Gamma(s+2)^2 \sin^2(s\pi) \mathcal{N}^2 \quad (24) $$

whose numerical value is

$$ \mathcal{M} \approx 2272. \quad (25) $$

Here $\tilde{N}_P(K)$ stands for the angular average $\int N_P(K) d\hat{K}/(4\pi)$ (with $\hat{K} := K/K$, and $d\hat{K}$ the differential solid angle, so that $d^3K = d\hat{K}K^2 dK$).

To derive this result, we consider the two-body reduced density matrix $\rho_2(\mathbf{r}_\uparrow, \mathbf{r}_\downarrow; \mathbf{r}'_\uparrow, \mathbf{r}'_\downarrow) := \langle \hat{\psi}_\uparrow^\dagger(\mathbf{r}'_\uparrow) \hat{\psi}_\downarrow(\mathbf{r}'_\downarrow) \hat{\psi}_\uparrow(\mathbf{r}_\uparrow) \hat{\psi}_\downarrow(\mathbf{r}_\downarrow) \rangle$ or equivalently in first quantization

$$ \rho_2(\mathbf{r}_\uparrow, \mathbf{r}_\downarrow; \mathbf{r}'_\uparrow, \mathbf{r}'_\downarrow) = N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \psi^*(\mathbf{r}'_\uparrow, \mathbf{r}_2, \mathbf{r}'_\downarrow, \mathbf{r}_4, \dots, \mathbf{r}_N) \psi(\mathbf{r}_\uparrow, \mathbf{r}_2, \mathbf{r}'_\downarrow, \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (26) $$

Inserting the two-body contact condition (14) into (26) yields

$$ \rho_2 \left( \mathbf{c} - \frac{\mathbf{r}}{2}, \mathbf{c} + \frac{\mathbf{r}}{2}; \mathbf{c}' - \frac{\mathbf{r}'}{2}, \mathbf{c}' + \frac{\mathbf{r}'}{2} \right) \underset{r,r' \to 0}{\sim} \frac{g_P(\mathbf{c}, \mathbf{c}')}{(4\pi)^2 r r'} \quad (27) $$

where

$$ g_P(\mathbf{c}, \mathbf{c}') = (4\pi)^2 N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N A^*(\mathbf{c}'; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (28) $$

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,

$$ N_P(\mathbf{K}) = \int d^3\mathbf{c} d^3\mathbf{c}' e^{-i\mathbf{K} \cdot (\mathbf{c}-\mathbf{c}')} g_P(\mathbf{c}, \mathbf{c}') \quad (29) $$

as can be formally shown using the definition (21) of $N_P$. Hence

$$ N_P(\mathbf{K}) = (4\pi)^2 N_\uparrow N_\downarrow \int d^3\mathbf{r}_2 d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \left| \int d^3\mathbf{c} e^{-i\mathbf{K}\cdot\mathbf{c}} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \right|^2. \quad (30) $$

We then follow a reasoning resembling the one used to derive Eq. (20) in [56, Sec. IV A].

In the large $K$ limit, the Fourier transform with respect to $\mathbf{c}$ in (30) is dominated by the contributions from the singularities of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$, which occur when $\mathbf{c}$ approaches one of the $\mathbf{r}_i$ ($i=2,4,\dots, N$). This corresponds to particles 1, 3 and $i$ being close to each other [since $A(\mathbf{c}; \mathbf{r}_2, \dots)$ determines the wavefunction $\psi$ when particles 1 and 3 are close to $\mathbf{c}$, according to Eq. (14)]. For example, for $i=2$, the behavior of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$ in the limit $\mathbf{c} \to \mathbf{r}_2$ is determined by the asymptotic behavior of $\psi$ 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

$$ A(\mathbf{r}_2 - \mathbf{u}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \underset{u \to 0}{\sim} \mathcal{N} \frac{s}{2} \cos\left(\frac{s\pi}{2}\right) \left(\frac{2u}{\sqrt{3}}\right)^{s-1} \sum_{m=-1}^{1} Y_1^m(\hat{\mathbf{u}}) B_m(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N) $$

where we used the expression ((99), (100)) of $\phi_m(\Omega)$. There is a similar singularity of $A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N)$ when $\mathbf{c}$ approaches $\mathbf{r}_i$ with $4 \le i \le N$; when particle $i$ has spin $\downarrow$, the function $\tilde{B}_m$ introduced in (18) appears instead of $B_m$. This gives

$$ \int d^3\mathbf{c} e^{-i\mathbf{K}\cdot\mathbf{c}} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N) \underset{K \to \infty}{\sim} \frac{2^{s-2}}{3^{(s-1)/2}} \mathcal{N} \cos\left(\frac{s\pi}{2}\right) \sum_{m=-1}^{1} I_m(\mathbf{K}) \times \left[ \sum_{i:\uparrow, i \neq \uparrow} (1-2\delta_{i,2}) e^{-i\mathbf{K}\cdot\mathbf{r}_i} B_m(P_{2i}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)) + \sum_{i:\downarrow, i \neq \downarrow} (1-2\delta_{i,4}) e^{-i\mathbf{K}\cdot\mathbf{r}_i} \tilde{B}_m(P_{4i}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)) \right] \quad (31) $$

where the first (resp. second) sum over $i$ is taken over particles with spin $\uparrow$ (resp. $\downarrow$), $P_{ji}(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)$ is obtained from $(\mathbf{r}_2; \mathbf{r}_4, \dots, \mathbf{r}_N)$ by exchanging $\mathbf{r}_j$ with $\mathbf{r}_i$ ($P_{ii}$ is the identity),
and $I_m(\mathbf{K}) := s \int d^3\mathbf{u} e^{i\mathbf{K}\cdot\mathbf{u}} u^{s-1} Y_1^m(\hat{\mathbf{u}})$. The latter integral can be evaluated analytically: Using $\int d\hat{\mathbf{u}} e^{i\mathbf{K}\cdot\mathbf{u}} Y_1^m(\hat{\mathbf{u}}) = 4\pi i Y_1^m(\hat{\mathbf{K}}) j_1(Ku)$ 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

$$I_m(\mathbf{K}) = 4\pi i \Gamma(s+2) \sin\left(\frac{s\pi}{2}\right) \frac{Y_1^m(\hat{\mathbf{K}})}{K^{s+2}}. \qquad (32)$$

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$,<fnref n="10" /> we obtain the result (23), (24), (25), where we used the value of $\mathcal{N}$ given in Eq. (101) of App. A.

## 2.5. Large-momentum tail of the two-particle momentum distribution

The three-body contact also determines the asymptotic behavior at large momenta of the two-particle opposite-spin momentum distribution function, defined by

$$N(\mathbf{k}_1, \mathbf{k}_2) = \langle \hat{N}_1(\mathbf{k}_1) \hat{N}_1(\mathbf{k}_2) \rangle \qquad (33)$$

where $\hat{N}_\sigma(\mathbf{k}) := \hat{c}_\sigma^\dagger(\mathbf{k}) \hat{c}_\sigma(\mathbf{k})$ with $\hat{c}_\sigma(\mathbf{k}) = \int d^3\mathbf{r} e^{-i\mathbf{k}\cdot\mathbf{r}} \hat{\psi}_\sigma(\mathbf{r})$ the annihilation operator of a particle of spin $\sigma$ in the state $|\mathbf{k}\rangle$ defined by $\langle\mathbf{r}|\mathbf{k}\rangle = e^{i\mathbf{k}\cdot\mathbf{r}}$.

Since $\int \hat{N}_\sigma(\mathbf{k}) \frac{d^3\mathbf{k}}{(2\pi)^3} = N_\sigma$, we have the normalization $\int N(\mathbf{k}_1, \mathbf{k}_2) \frac{d^3\mathbf{k}_1}{(2\pi)^3} \frac{d^3\mathbf{k}_2}{(2\pi)^3} = N_1 N_2$.

Experimentally, the two-particle momentum distribution can be accessed from the statistics of time-of-flight images, as was recently demonstrated in 2D [122].<fnref n="11" />

Taking the limit of a large relative momentum $k \to \infty$, and integrating over the center-of-mass momentum $\mathbf{K}$, one obtains a tail proportional to $C_2$,

$$\int N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) \frac{d^3\mathbf{K}}{(2\pi)^3} \stackrel{k \to \infty}{\sim} \frac{C_2}{k^4} \qquad (34)$$

as pointed out in [128]. If we instead send $K$ to infinity, and average over the direction of $\mathbf{K}$, we obtain a tail proportional to $C_3$,

$$\lim_{K \to \infty} \lim_{k \to \infty} K^{2s+4} k^4 \frac{1}{4\pi} \int d\hat{\mathbf{K}} N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) = \mathcal{M} C_3 \qquad (35)$$

where the constant $\mathcal{M}$ was given in (24).

This result immediately follows from (23) and from the relation

$$N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) \stackrel{k \to \infty}{\sim} \frac{N_P(\mathbf{K})}{k^4}, \qquad (36)$$

obtained in [128] and rederived in the sequel.<fnref n="12" />,<fnref n="13" /> The definition (33) yields

$$N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) = \int d^3\mathbf{c} \, d^3\mathbf{c}' \, d^3\mathbf{r} \, d^3\mathbf{r}' e^{i\mathbf{K}\cdot(\mathbf{c}' - \mathbf{c})} e^{i\mathbf{k}\cdot(\mathbf{r}' - \mathbf{r})} \rho_2\left(\mathbf{c} - \frac{\mathbf{r}}{2}, \mathbf{c} + \frac{\mathbf{r}}{2}; \mathbf{c}' - \frac{\mathbf{r}'}{2}, \mathbf{c}' + \frac{\mathbf{r}'}{2}\right) \qquad (38)$$

$$N_P(\mathbf{K}) = (4\pi)^2 \int \frac{d^3 k_2}{(2\pi)^3} \frac{d^3 k_4}{(2\pi)^3} \cdots \frac{d^3 k_N}{(2\pi)^3} |\tilde{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N)|^2 \qquad (37)$$

with $\tilde{A}$ the Fourier transform of $A$

$$\tilde{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N) := \int d^3\mathbf{c} \, d^3\mathbf{r}_2 \, d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N e^{-i(\mathbf{K}\cdot\mathbf{c} + \mathbf{k}_2\cdot\mathbf{r}_2 + \mathbf{k}_4\cdot\mathbf{r}_4 + \dots + \mathbf{k}_N\cdot\mathbf{r}_N)} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N),$$

<footnote n="10">¹⁰By power counting, these cross terms give rise to a $\propto 1/K^{2s_4+4}$ tail of $\tilde{N}_P(\mathbf{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 $\tilde{N}_P(\mathbf{K})$, given that $s_4 = 2.509(1)$ is larger than $s \equiv 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)\hbar\omega$.</footnote>

<footnote n="11">¹¹In 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–127].</footnote>

<footnote n="12">¹²Relation (34) also follows from (36), given (22).</footnote>

<footnote n="13">¹³In Ref. [128], $N_P(\mathbf{K})$ was defined by</footnote>
and $I_m(\mathbf{K}) := s \int d^3\mathbf{u} e^{i\mathbf{K}\cdot\mathbf{u}} u^{s-1} Y_1^m(\hat{\mathbf{u}})$. The latter integral can be evaluated analytically: Using $\int d\hat{\mathbf{u}} e^{i\mathbf{K}\cdot\mathbf{u}} Y_1^m(\hat{\mathbf{u}}) = 4\pi i Y_1^m(\hat{\mathbf{K}}) j_1(Ku)$ 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

$$I_m(\mathbf{K}) = 4\pi i \Gamma(s+2) \sin\left(\frac{s\pi}{2}\right) \frac{Y_1^m(\hat{\mathbf{K}})}{K^{s+2}}. \qquad (32)$$

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$,<fnref n="10" /> we obtain the result (23), (24), (25), where we used the value of $\mathcal{N}$ given in Eq. (101) of App. A.

## 2.5. Large-momentum tail of the two-particle momentum distribution

The three-body contact also determines the asymptotic behavior at large momenta of the two-particle opposite-spin momentum distribution function, defined by

$$N(\mathbf{k}_1, \mathbf{k}_2) = \langle \hat{N}_1(\mathbf{k}_1) \hat{N}_1(\mathbf{k}_2) \rangle \qquad (33)$$

where $\hat{N}_\sigma(\mathbf{k}) := \hat{c}_\sigma^\dagger(\mathbf{k}) \hat{c}_\sigma(\mathbf{k})$ with $\hat{c}_\sigma(\mathbf{k}) = \int d^3\mathbf{r} e^{-i\mathbf{k}\cdot\mathbf{r}} \hat{\psi}_\sigma(\mathbf{r})$ the annihilation operator of a particle of spin $\sigma$ in the state $|\mathbf{k}\rangle$ defined by $\langle\mathbf{r}|\mathbf{k}\rangle = e^{i\mathbf{k}\cdot\mathbf{r}}$.

Since $\int \hat{N}_\sigma(\mathbf{k}) \frac{d^3\mathbf{k}}{(2\pi)^3} = N_\sigma$, we have the normalization $\int N(\mathbf{k}_1, \mathbf{k}_2) \frac{d^3\mathbf{k}_1}{(2\pi)^3} \frac{d^3\mathbf{k}_2}{(2\pi)^3} = N_1 N_2$.

Experimentally, the two-particle momentum distribution can be accessed from the statistics of time-of-flight images, as was recently demonstrated in 2D [122].<fnref n="11" />

Taking the limit of a large relative momentum $k \to \infty$, and integrating over the center-of-mass momentum $\mathbf{K}$, one obtains a tail proportional to $C_2$,

$$\int N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) \frac{d^3\mathbf{K}}{(2\pi)^3} \stackrel{k \to \infty}{\sim} \frac{C_2}{k^4} \qquad (34)$$

as pointed out in [128]. If we instead send $K$ to infinity, and average over the direction of $\mathbf{K}$, we obtain a tail proportional to $C_3$,

$$\lim_{K \to \infty} \lim_{k \to \infty} K^{2s+4} k^4 \frac{1}{4\pi} \int d\hat{\mathbf{K}} N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) = \mathcal{M} C_3 \qquad (35)$$

where the constant $\mathcal{M}$ was given in (24).

This result immediately follows from (23) and from the relation

$$N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) \stackrel{k \to \infty}{\sim} \frac{N_P(\mathbf{K})}{k^4}, \qquad (36)$$

obtained in [128] and rederived in the sequel.<fnref n="12" />,<fnref n="13" /> The definition (33) yields

$$N\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \frac{\mathbf{K}}{2} + \mathbf{k}\right) = \int d^3\mathbf{c} d^3\mathbf{c}' d^3\mathbf{r} d^3\mathbf{r}' e^{i\mathbf{K}\cdot(\mathbf{c}' - \mathbf{c})} e^{i\mathbf{k}\cdot(\mathbf{r}' - \mathbf{r})} \rho_2\left(\mathbf{c} - \frac{\mathbf{r}}{2}, \mathbf{c} + \frac{\mathbf{r}}{2}; \mathbf{c}' - \frac{\mathbf{r}'}{2}, \mathbf{c}' + \frac{\mathbf{r}'}{2}\right) \qquad (38)$$

$$N_P(\mathbf{K}) = (4\pi)^2 \int \frac{d^3 k_2}{(2\pi)^3} \frac{d^3 k_4}{(2\pi)^3} \cdots \frac{d^3 k_N}{(2\pi)^3} |\tilde{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N)|^2 \qquad (37)$$

with $\tilde{A}$ the Fourier transform of $A$

$$\tilde{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N) := \int d^3\mathbf{c} \, d^3\mathbf{r}_2 \, d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N e^{-i(\mathbf{K}\cdot\mathbf{c} + \mathbf{k}_2\cdot\mathbf{r}_2 + \mathbf{k}_4\cdot\mathbf{r}_4 + \dots + \mathbf{k}_N\cdot\mathbf{r}_N)} A(\mathbf{c}; \mathbf{r}_2, \mathbf{r}_4, \dots, \mathbf{r}_N),$$

<footnote n="10">¹⁰By power counting, these cross terms give rise to a $\propto 1/K^{2s_4+4}$ tail of $\tilde{N}_P(\mathbf{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 $\tilde{N}_P(\mathbf{K})$, given that $s_4 = 2.509(1)$ is larger than $s \equiv 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)\hbar\omega$.</footnote>

<footnote n="11">¹¹In 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–127].</footnote>

<footnote n="12">¹²Relation (34) also follows from (36), given (22).</footnote>

<footnote n="13">¹³In Ref. [128], $N_P(\mathbf{K})$ was defined by</footnote>
where $\rho_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 \to \infty$ tail of the Fourier transform Eq. (38), which can be computed by replacing $\rho_2$ with its asymptotic expression (27), and using the identity (in the sense of distributions) $\int d^3r e^{-ik \cdot r}/r = 4\pi/k^2$. Using (29) then yields (36).

### 3. Relations for finite-range models

In this section, we go beyond the zero-range model, and consider interactions of small but non-zero 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).

#### 3.1. 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.

##### 3.1.1. 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$.<fnref n="14" /> We consider the resonant regime where the two-body scattering length $a_2$ is large,

$$|a_2| \gg b. \qquad (39)$$

In this regime, there are two kinds of two-body bound states:

* the weakly bound dimer, of binding energy $\approx \hbar^2/(ma_2^2)$, which exists for $a_2 > 0$

* deeply bound dimer(s), of binding energy $\gtrsim \hbar^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

$$H \psi = E \psi \qquad (40)$$

$$H = \sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] + \sum_{i:\uparrow, j:\downarrow} V_2(r_{ij}) \qquad (41)$$

in the zero-range regime

$$1/k_{\text{typ}} \gg b \qquad (42)$$

and (36) was deduced from the expression

$$N(\mathbf{k}_1, \mathbf{k}_3) = N_1 N_3 \int \frac{d^3 k_2}{(2\pi)^3} \frac{d^3 k_4}{(2\pi)^3} \cdots \frac{d^3 k_N}{(2\pi)^3} |\bar{\psi}(\mathbf{k}_1, \dots, \mathbf{k}_N)|^2$$

in terms of the momentum-space wavefunction

$$\bar{\psi}(\mathbf{k}_1, \dots, \mathbf{k}_N) := \int d^3\mathbf{r}_1 \dots d^3\mathbf{r}_N e^{-i(\mathbf{k}_1 \cdot \mathbf{r}_1 + \dots + \mathbf{k}_N \cdot \mathbf{r}_N)} \psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$$

which has the asymptotic behavior

$$\bar{\psi}\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \mathbf{k}_2, \frac{\mathbf{K}}{2} + \mathbf{k}, \mathbf{k}_4, \dots, \mathbf{k}_N\right) \underset{k \to \infty}{\sim} \frac{4\pi}{k^2} \bar{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N).$$

<footnote n="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 \gg b_0$). More generally, $b = \text{Max}(b_0, |r_e|)$ where $r_e$ is the effective range.</footnote>
where $\rho_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 \to \infty$ tail of the Fourier transform Eq. (38), which can be computed by replacing $\rho_2$ with its asymptotic expression (27), and using the identity (in the sense of distributions) $\int d^3r e^{-ik \cdot r}/r = 4\pi/k^2$. Using (29) then yields (36).

### 3. Relations for finite-range models

In this section, we go beyond the zero-range model, and consider interactions of small but non-zero 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).

#### 3.1. 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.

##### 3.1.1. 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$.<fnref n="14" /> We consider the resonant regime where the two-body scattering length $a_2$ is large,

$$|a_2| \gg b. \qquad (39)$$

In this regime, there are two kinds of two-body bound states:

* the weakly bound dimer, of binding energy $\approx \hbar^2/(ma_2^2)$, which exists for $a_2 > 0$

* deeply bound dimer(s), of binding energy $\gtrsim \hbar^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

$$H \psi = E \psi \qquad (40)$$

$$H = \sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] + \sum_{i:\uparrow, j:\downarrow} V_2(r_{ij}) \qquad (41)$$

in the zero-range regime

$$1/k_{\text{typ}} \gg b \qquad (42)$$

and (36) was deduced from the expression

$$N(\mathbf{k}_1, \mathbf{k}_3) = N_1 N_3 \int \frac{d^3 k_2}{(2\pi)^3} \frac{d^3 k_4}{(2\pi)^3} \cdots \frac{d^3 k_N}{(2\pi)^3} |\bar{\psi}(\mathbf{k}_1, \dots, \mathbf{k}_N)|^2$$

in terms of the momentum-space wavefunction

$$\bar{\psi}(\mathbf{k}_1, \dots, \mathbf{k}_N) := \int d^3\mathbf{r}_1 \dots d^3\mathbf{r}_N e^{-i(\mathbf{k}_1 \cdot \mathbf{r}_1 + \dots + \mathbf{k}_N \cdot \mathbf{r}_N)} \psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$$

which has the asymptotic behavior

$$\bar{\psi}\left(\frac{\mathbf{K}}{2} - \mathbf{k}, \mathbf{k}_2, \frac{\mathbf{K}}{2} + \mathbf{k}, \mathbf{k}_4, \dots, \mathbf{k}_N\right) \underset{k \to \infty}{\sim} \frac{4\pi}{k^2} \bar{A}(\mathbf{K}; \mathbf{k}_2, \mathbf{k}_4, \dots, \mathbf{k}_N).$$

<footnote n="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 \gg b_0$). More generally, $b = \text{Max}(b_0, |r_e|)$ where $r_e$ is the effective range.</footnote>
where $1/k_{typ}$ is defined as the smallest scale of variation of the stationary wavefunction $\psi$ in the region where all interparticle distances are $\gg b$.<fnref n="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 zero-range regime (42), the zero-range model is valid for any stationary state $\psi$, 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(\epsilon)$ and $N_3(\epsilon)$ provided $\epsilon \gg b$ [to reach the asymptotic regimes of Eqs. (1),(3) one also needs $\epsilon \ll 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 $\Gamma_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 ($\sim \frac{\hbar^2}{mb^2}$) of deeply bound dimers. In typical experiments, this condition holds, and other loss processes are negligible, which allows one to measure $\Gamma_3$ from the decay of the number of trapped atoms, $\dot{N} = -3\Gamma_3$ [11, 129–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–138]. Accordingly, we will consider a solution $\psi$ of (40), (41) 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.<fnref n="16" /> For such a quasi-stationary state, in the zero-range regime, standard observables again tend to their respective values within the zero-range model.<fnref n="17" />,<fnref n="18" /> The three-body loss rate $\Gamma_3$, however, is simply zero within the zero-range model. To compute $\Gamma_3$, one thus needs to go beyond the zero-range model. As we will see, one only needs to do so for the three-body problem, in order to define a three-body parameter $a_3$. We then find

$$ \Gamma_3 \approx -\frac{\hbar}{m} 8s(s+1) C_3 \operatorname{Im} a_3 \quad (43) $$

where $C_3$ can be evaluated within the zero-range model. Furthermore, breaking up $\Gamma_3$ into the sum of the two contributions $\Gamma_{2,1}$ and $\Gamma_{1,2}$ corresponding to $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ loss processes, we have

$$ \Gamma_{2,1} \approx -\frac{\hbar}{m} 8s(s+1) C_{2,1} \operatorname{Im} a_3 \quad (44) $$

$$ \Gamma_{1,2} \approx -\frac{\hbar}{m} 8s(s+1) C_{1,2} \operatorname{Im} a_3. \quad (45) $$

We expect these relations to be asymptotically exact in the resonant zero-range regime (39),(42).

<footnote n="15">¹⁵For example, for the ground state of the homogeneous unpolarized gas, $k_{typ}$ is $\sim k_F$ for $a_2 < 0$ and $\sim \operatorname{Max}(k_F, 1/a_2)$ 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 $\omega$, $1/k_{typ}$ is $\sim a_{ho}$ for $a_2 < 0$ and $\sim \operatorname{Min}(a_{ho}, a_2)$ for $a_2 > 0$, where $a_{ho} := \sqrt{\hbar/(m\omega)}$ is the harmonic oscillator length.</footnote>

<footnote n="16">¹⁶An 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].</footnote>

<footnote n="17">¹⁷An appropriate normalization of the Gamow state will be given below in Eq. (59). Similarly, the expectation value of an observable $\hat{O}$ should be defined as $\int_{\mathcal{R}} d^{3N}X \int_{\mathcal{R}} d^{3N}X' \psi^*(X) \langle X|\hat{O}|X'\rangle \psi(X')$.</footnote>

<footnote n="18">¹⁸Within the zero-range model, it is convenient to add steep infinite walls to the trapping potential at the boundary of $\mathcal{R}_{trap}$, in order to have truly stationary states, thereby neglecting the exponentially suppressed evaporation process discussed in Fn. 23.</footnote>
where $1/k_{typ}$ is defined as the smallest scale of variation of the stationary wavefunction $\psi$ in the region where all interparticle distances are $\gg b$.<fnref n="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 zero-range regime (42), the zero-range model is valid for any stationary state $\psi$, 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(\epsilon)$ and $N_3(\epsilon)$ provided $\epsilon \gg b$ [to reach the asymptotic regimes of Eqs. (1), (3) one also needs $\epsilon \ll 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 $\Gamma_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 ($\sim \frac{\hbar^2}{mb^2}$) of deeply bound dimers. In typical experiments, this condition holds, and other loss processes are negligible, which allows one to measure $\Gamma_3$ from the decay of the number of trapped atoms, $\dot{N} = -3\Gamma_3$ [11, 129–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–138]. Accordingly, we will consider a solution $\psi$ of (40), (41) 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.<fnref n="16" /> For such a quasi-stationary state, in the zero-range regime, standard observables again tend to their respective values within the zero-range model.<fnref n="17" />,<fnref n="18" /> The three-body loss rate $\Gamma_3$, however, is simply zero within the zero-range model. To compute $\Gamma_3$, one thus needs to go beyond the zero-range model. As we will see, one only needs to do so for the three-body problem, in order to define a three-body parameter $a_3$. We then find

$$ \Gamma_3 \approx -\frac{\hbar}{m} 8s(s+1) C_3 \operatorname{Im} a_3 \quad (43) $$

where $C_3$ can be evaluated within the zero-range model. Furthermore, breaking up $\Gamma_3$ into the sum of the two contributions $\Gamma_{2,1}$ and $\Gamma_{1,2}$ corresponding to $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ loss processes, we have

$$ \Gamma_{2,1} \approx -\frac{\hbar}{m} 8s(s+1) C_{2,1} \operatorname{Im} a_3 \quad (44) $$

$$ \Gamma_{1,2} \approx -\frac{\hbar}{m} 8s(s+1) C_{1,2} \operatorname{Im} a_3. \quad (45) $$

We expect these relations to be asymptotically exact in the resonant zero-range regime (39), (42).

<footnote n="15">¹⁵For example, for the ground state of the homogeneous unpolarized gas, $k_{typ}$ is $\sim k_F$ for $a_2 < 0$ and $\sim \operatorname{Max}(k_F, 1/a_2)$ 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 $\omega$, $1/k_{typ}$ is $\sim a_{ho}$ for $a_2 < 0$ and $\sim \operatorname{Min}(a_{ho}, a_2)$ for $a_2 > 0$, where $a_{ho} := \sqrt{\hbar/(m\omega)}$ is the harmonic oscillator length.</footnote>

<footnote n="16">¹⁶An 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].</footnote>

<footnote n="17">¹⁷An appropriate normalization of the Gamow state will be given below in Eq. (59). Similarly, the expectation value of an observable $\hat{O}$ should be defined as $\int_{\mathcal{R}} d^{3N}X \int_{\mathcal{R}} d^{3N}X' \psi^*(\mathbf{X}) \langle \mathbf{X}|\hat{O}|\mathbf{X}' \rangle \psi(\mathbf{X}')$.</footnote>

<footnote n="18">¹⁸Within the zero-range model, it is convenient to add steep infinite walls to the trapping potential at the boundary of $\mathcal{R}_{trap}$, in order to have truly stationary states, thereby neglecting the exponentially suppressed evaporation process discussed in Fn. 23.</footnote>
To define the three-body parameter $a_3$, we consider the zero-energy free-space solution of the Schrödinger equation (40), (41) for three particles of spins ↑↑↓ and angular-momentum quantum numbers ($l = 1, m$) whose asymptotic behavior has the form

$$ \Psi_m(\mathbf{R}) \approx \left(R^s - \frac{a_3}{R^s}\right) \frac{1}{R^2} \phi_m(\mathbf{\Omega}) \quad (46) $$

in the region {$b \ll r_{ij} \ll |a_2|$, $\forall 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 $\gg b$. The fact that $a_3$ does not depend on the quantum number $m$ follows from rotational invariance of the interaction.

### Remarks 3.

* Relation (43) is reminiscent of the known relation [44, 139] between two-body loss rate and two-body contact 19,20

$$ \Gamma_2 = \frac{\hbar}{2\pi m} C_2 \operatorname{Im}\left(\frac{1}{a_2}\right). \quad (47) $$

* $\operatorname{Im} a_3$ must be negative, since the loss rate is positive.

* Typically one has the order of magnitude estimate $a_3 \sim b^{2s}$, assuming that there is no extra fine-tuning.21 Therefore $\Gamma_3 \propto 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 $\lesssim b$, times $\hbar/(mb^2)$, that is, $\Gamma_3 = \mathcal{K} N_3(b) \hbar/(mb^2)$, with $N_3(b)$ evaluated within the zero-range model, and $\mathcal{K}$ a dimensionless prefactor that depends on short-range three-body physics. This statement is equivalent to (43), given the relation (3), with $\mathcal{K} = -8s(s+1) \operatorname{Im} a_3/b^{2s}$. The novelties of the present work are (i) to provide a derivation of the relation (43), and hence of the above statement from [98], (ii) to introduce the natural single parameter $a_3$ (instead of the two parameters $\mathcal{K}$ and $b$), and (iii) to generalize the relation to more complex interactions (in the following Sec. 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$

<footnote n="19">¹⁹Relation (47) 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.</footnote>

<footnote n="20">²⁰ Relation (47) was obtained in [44] using the relation (85) and $\operatorname{Im} E = -\hbar\Gamma_2/2$. It was rederived in [139] using the Lindblad equation. For completeness, we note that the relation (47) can also be derived by a flux computation. The main steps of this derivation are as follows. We consider a solution $\psi$ of the zero-range model (13), (14) with complex values of $1/a_2$ and $E$. The corresponding time-dependent wavefunction is $\Psi(t) = \psi e^{-iEt/\hbar}$, and we have $\Gamma_2 = -\partial_t \langle \Psi(t) | \Psi(t) \rangle|_{t=0}$. Using (56), (57), we obtain that $\Gamma_2$ equals $N_1 N_1$ times the limit when $\epsilon \to 0$ of the probability flux entering into the region $\{ (\mathbf{r}_1, \dots, \mathbf{r}_N) | r < \epsilon \}$, which can be simplified to $\Gamma_2 = -2N_1 N_1 \frac{\hbar}{m} \lim_{\epsilon \to 0} \int d^3c \, d^3r_2 \, d^3r_4 \dots d^3r_N e^2 \int d\mathbf{r} \, \operatorname{Im}(\psi^* \partial_r \psi)|_{r=\epsilon}$. Using (14) then yields (47).</footnote>

<footnote n="21">²¹Indeed, the behavior (46) has to be matched at $R$ of order $b$ with the solution inside the potential, which typically imposes that the two terms in (46) are of the same order of magnitude for $R \sim b$ [97, 140]. For simplicity, we excluded here the special regime $|a_3| \gg b^{2s}$ that corresponds to the vicinity of a three-body resonance (see [141, 142], and [38, 117, 143–148] for the mass-imbalanced case with $0 \le 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| \gg 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 (43), (44), (45) and the other results of this article to remain applicable in the three-body resonant regime provided the three-body parameter(s) remain(s) small (in modulus) compared to $1/k_{\text{typ}}^{2s}$.</footnote>
To define the three-body parameter $a_3$, we consider the zero-energy free-space solution of the Schrödinger equation (40), (41) for three particles of spins ↑↑↓ and angular-momentum quantum numbers ($l = 1, m$) whose asymptotic behavior has the form

$$ \Psi_m(\mathbf{R}) \approx \left(R^s - \frac{a_3}{R^s}\right) \frac{1}{R^2} \phi_m(\mathbf{\Omega}) \qquad (46) $$

in the region {$b \ll r_{ij} \ll |a_2|, \forall 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 $\gg b$. The fact that $a_3$ does not depend on the quantum number $m$ follows from rotational invariance of the interaction.

### Remarks 3.

* Relation (43) is reminiscent of the known relation [44, 139] between two-body loss rate and two-body contact 19,20

$$ \Gamma_2 = \frac{\hbar}{2\pi m} C_2 \operatorname{Im}\left(\frac{1}{a_2}\right). \qquad (47) $$

* Im $a_3$ must be negative, since the loss rate is positive.

* Typically one has the order of magnitude estimate $a_3 \sim b^{2s}$, assuming that there is no extra fine-tuning.21 Therefore $\Gamma_3 \propto 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 $\lesssim b$, times $\hbar/(m b^2)$, that is, $\Gamma_3 = \mathcal{K} N_3(b) \hbar/(m b^2)$, with $N_3(b)$ evaluated within the zero-range model, and $\mathcal{K}$ a dimensionless prefactor that depends on short-range three-body physics. This statement is equivalent to (43), given the relation (3), with $\mathcal{K} = -8s(s+1) \operatorname{Im} a_3/b^{2s}$. The novelties of the present work are (i) to provide a derivation of the relation (43), and hence of the above statement from [98], (ii) to introduce the natural single parameter $a_3$ (instead of the two parameters $\mathcal{K}$ and $b$), and (iii) to generalize the relation to more complex interactions (in the following Sec. 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$

<footnote n="19">¹⁹Relation (47) 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.</footnote>

<footnote n="20">²⁰ Relation (47) was obtained in [44] using the relation (85) and Im $E = -\hbar\Gamma_2/2$. It was rederived in [139] using the Lindblad equation. For completeness, we note that the relation (47) can also be derived by a flux computation. The main steps of this derivation are as follows. We consider a solution $\psi$ of the zero-range model (13), (14) with complex values of $1/a_2$ and $E$. The corresponding time-dependent wavefunction is $\Psi(t) = \psi e^{-iEt/\hbar}$, and we have $\Gamma_2 = -\partial_t \langle \Psi(t) | \Psi(t) \rangle|_{t=0}$. Using (56), (57), we obtain that $\Gamma_2$ equals $N_1 N_1$ times the limit when $\epsilon \to 0$ of the probability flux entering into the region $\{ (\mathbf{r}_1, \dots, \mathbf{r}_N) | r < \epsilon \}$, which can be simplified to $\Gamma_2 = -2N_1 N_1 \frac{\hbar}{m} \lim_{\epsilon \to 0} \int d^3c \, d^3r_2 \, d^3r_4 \dots d^3r_N e^2 \int d\mathbf{r} \, \operatorname{Im}(\psi^* \partial_r \psi)|_{r=\epsilon}$. Using (14) then yields (47).</footnote>

<footnote n="21">²¹Indeed, the behavior (46) has to be matched at $R$ of order $b$ with the solution inside the potential, which typically imposes that the two terms in (46) are of the same order of magnitude for $R \sim b$ [97, 140]. For simplicity, we excluded here the special regime $|a_3| \gg b^{2s}$ that corresponds to the vicinity of a three-body resonance (see [141, 142], and [38, 117, 143–148] for the mass-imbalanced case with $0 \le 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| \gg 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 (43), (44), (45) and the other results of this article to remain applicable in the three-body resonant regime provided the three-body parameter(s) remain(s) small (in modulus) compared to $1/k_{typ}^{2s}$.</footnote>
the volume, dimensional analysis gives $\mathcal{C}_3 = \zeta_3 n^{(2s+5)/3}$ where $\zeta_3$ is a dimensionless constant. Hence

$$ -\dot{n} = 24s(s+1)\frac{\hbar}{m}(-\mathrm{Im}\,a_3)\zeta_3 n^{(2s+5)/3}. \quad (48) $$

Therefore, as already found in [98], the timescale of three-body losses $\tau_3 := n/|\dot{n}|$ is of order $\tau_F/(k_F b)^{2s}$, much larger than the thermalization time $\tau_F \sim m/(\hbar k_F^2)$, 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–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 (43) of $\Gamma_3$ is reminiscent of the relation for bosons expressing $\Gamma_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 $\Gamma_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–158]. Presumably, $D$ could be defined also for the present case (two-component fermions in 3D), and as in [153–158], $D$ would govern the asymptotic behavior of the three-body zero-energy wavefunction at distances $\gg |a_2|$ (while $a_3$ governs the regime of distances $\ll |a_2|$ and $\gg b$) and $\Gamma_3$ would have a simple expression in terms of $\mathrm{Im}\,D$ in the weakly interacting regime $k_{\mathrm{typ}}|a_2| \ll 1$ (whereas (43) remains valid in the strongly correlated regime $k_{\mathrm{typ}}|a_2| \gtrsim 1$). On the other hand, $a_3$ is only defined for resonant interactions ($|a_2| \gg b$) while $D$ also exists for non-resonant interactions.

We turn to the derivation of (43), (44), (45). Our reasoning is similar to the bosonic case treated in [72, App. B], but the present case is significantly more complicated because we cannot work directly within the zero-range model. For $R \gg b$, the $a_3/R^s$ term in (46) 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 $\psi$, we go beyond the zero-range model and replace (15) with

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \left(R^s - \frac{a_3}{R^s}\right) \frac{1}{R^2} \sum_{m=-1}^{+1} \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) \quad (49) $$

which we expect to hold provided

$$ R \le d_3 \quad (50) $$

$$ b \ll r_{ij}, \quad \forall i < j \le 3 \quad (51) $$

$$ d_3 \ll r_{ij}, \quad \forall i < j, j \ge 4 \quad (52) $$

where $d_3$ is a length that satisfies

$$ b \ll d_3 \quad (53) $$

$$ d_3 \ll \mathrm{Min}(|a_2|, 1/k_{\mathrm{typ}}) \quad (54) $$

and whose appropriate choice will be discussed more precisely below. The purpose of (52) is to ensure that the interparticle distances within the triplet of particles 1, 2, 3 are much smaller than all other interparticle distances.
the volume, dimensional analysis gives $\mathcal{C}_3 = \zeta_3 n^{(2s+5)/3}$ where $\zeta_3$ is a dimensionless constant. Hence

$$ -\dot{n} = 24s(s+1) \frac{\hbar}{m} (-\mathrm{Im} a_3) \zeta_3 n^{(2s+5)/3}. \quad (48) $$

Therefore, as already found in [98], the timescale of three-body losses $\tau_3 := n/|\dot{n}|$ is of order $\tau_F/(k_F b)^{2s}$, much larger than the thermalization time $\tau_F \sim m/(\hbar k_F^2)$, 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–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 (43) of $\Gamma_3$ is reminiscent of the relation for bosons expressing $\Gamma_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 $\Gamma_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–158]. Presumably, $D$ could be defined also for the present case (two-component fermions in 3D), and as in [153–158], $D$ would govern the asymptotic behavior of the three-body zero-energy wavefunction at distances $\gg |a_2|$ (while $a_3$ governs the regime of distances $\ll |a_2|$ and $\gg b$) and $\Gamma_3$ would have a simple expression in terms of $\mathrm{Im} D$ in the weakly interacting regime $k_{\mathrm{typ}}|a_2| \ll 1$ (whereas (43) remains valid in the strongly correlated regime $k_{\mathrm{typ}}|a_2| \gtrsim 1$). On the other hand, $a_3$ is only defined for resonant interactions ($|a_2| \gg b$) while $D$ also exists for non-resonant interactions.

We turn to the derivation of (43), (44), (45). Our reasoning is similar to the bosonic case treated in [72, App. B], but the present case is significantly more complicated because we cannot work directly within the zero-range model. For $R \gg b$, the $a_3/R^s$ term in (46) 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 $\psi$, we go beyond the zero-range model and replace (15) with

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \left(R^s - \frac{a_3}{R^s}\right) \frac{1}{R^2} \sum_{m=-1}^{+1} \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) \quad (49) $$

which we expect to hold provided

$$ R \le d_3 \quad (50) $$

$$ b \ll r_{ij}, \quad \forall i < j \le 3 \quad (51) $$

$$ d_3 \ll r_{ij}, \quad \forall i < j, j \ge 4 \quad (52) $$

where $d_3$ is a length that satisfies

$$ b \ll d_3 \quad (53) $$

$$ d_3 \ll \mathrm{Min}(|a_2|, 1/k_{\mathrm{typ}}) \quad (54) $$

and whose appropriate choice will be discussed more precisely below. The purpose of (52) 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 $\mathbf{X} = (\mathbf{r}_1, \dots, \mathbf{r}_N)$. The solution $\Psi(\mathbf{X}; t)$ of the time-dependent Schrödinger equation $[i\hbar\dot{\Psi} = H\Psi$, with $H$ given by (41)] associated with the Gamow state $\psi(\mathbf{X})$ is

$$ \Psi(\mathbf{X}; t) = \psi(\mathbf{X}) e^{-iEt/\hbar}. \qquad (55) $$

It satisfies the continuity equation

$$ \partial_t (|\Psi|^2) = -\nabla_X \cdot \mathbf{J} \qquad (56) $$

in terms of the probability current

$$ \mathbf{J} := \frac{\hbar}{m} \operatorname{Im} (\Psi^* \nabla_X \Psi). \qquad (57) $$

We express the loss event rate (i.e. the probability for a loss event to occur per unit of time) as

$$ \Gamma = - \left. \frac{\partial}{\partial t} \right|_{t=0} \int_{\mathcal{R}} d^{3N}X |\Psi(\mathbf{X}; t)|^2 \qquad (58) $$

where the region $\mathcal{R}$ should physically correspond to $N$ atoms in the trap, and where $\psi(\mathbf{X}) = \Psi(\mathbf{X}, t=0)$ is normalized by the condition

$$ \int_{\mathcal{R}} d^{3N}X |\psi(\mathbf{X})|^2 = 1. \qquad (59) $$

There is some freedom in how to define $\mathcal{R}$. Let $\mathcal{R}_{\text{trap}}$ denote the “trapping region” of the potential $U(\mathbf{r})$, which can be defined as the set of points $\mathbf{r}$ such that the classical trajectory of a particle with initial position $\mathbf{r}$ and zero initial velocity remains bounded.<fnref n="22" /> A possible definition of $\mathcal{R}$ would be $\mathcal{R} = (\mathcal{R}_{\text{trap}})^N$, but for later convenience, we choose $\mathcal{R}$ to be slightly smaller, by excluding configurations with two or three nearby particles:

$$ \mathcal{R} = (\mathcal{R}_{\text{trap}})^N \setminus (\mathcal{R}_2 \cup \mathcal{R}_3) $$

where

$$ \mathcal{R}_2 = \bigcup_{i:\uparrow, j:\downarrow} B_{ij} \quad \text{with} \quad B_{ij} = \{(\mathbf{r}_1, \dots, \mathbf{r}_N) \mid r_{ij} < d_2\} \qquad (60) $$

$$ \mathcal{R}_3 = \bigcup_{i<j:\uparrow, k:\downarrow \text{ or } i<j:\downarrow, k:\uparrow} B_{ijk} \quad \text{with} \quad B_{ijk} = \{(\mathbf{r}_1, \dots, \mathbf{r}_N) \mid R_{ijk} < d_3\}. \qquad (61) $$

In other words, $(\mathbf{r}_1, \dots, \mathbf{r}_N) \in \mathcal{R}$ means that

* all particle positions $\mathbf{r}_i$ are inside the trapping region $\mathcal{R}_{\text{trap}}$

* all distances between pairs of particles with opposite spin are $\ge d_2$

* all hyperradii of triplets of particles with non-identical spins are $\ge d_3$.

We take $d_2$ such that

$$ b \ll d_2 \qquad (62) $$

$$ d_2 \ll d_3. \qquad (63) $$

Conditions (54), (63) ensure that both $d_2$ and $d_3$ are $\ll 1/k_{\text{typ}}$. Hence, to leading order, the normalization integral (59) is independent of $d_2$ and $d_3$, because $\mathcal{R}_2$ and $\mathcal{R}_3$ are negligibly small subsets of $(\mathcal{R}_{\text{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 $\psi(\mathbf{X})$ for $\mathbf{X} \in \mathcal{R}$

* in (49), $B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)$ can be evaluated within the zero-range model, provided $\mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$ all belong to $\mathcal{R}_{\text{trap}}$; indeed, the $\propto a_3$ term is a small correction, and the condition (49) for the Gamow state must match with the condition (15) for the zero-range model.

<footnote n="22">²²For example, if $U(\mathbf{r}) = m\omega^2 r^2 - \lambda r^4$, then $\mathcal{R}_{\text{trap}}$ is the sphere centered at the origin of radius $r_{\text{trap}}$ such that $U(r)$ reaches its maximum at $r = r_{\text{trap}}$, i.e. $r_{\text{trap}} = \sqrt{m/\lambda} \times \omega/2$. This definition of $\mathcal{R}_{\text{trap}}$ is merely one particularly simple and natural choice among a range of possibilities.</footnote>
We use the shorthand notation $\mathbf{X} = (\mathbf{r}_1, \dots, \mathbf{r}_N)$. The solution $\Psi(\mathbf{X}; t)$ of the time-dependent Schrödinger equation $[i\hbar \dot{\Psi} = H\Psi$, with $H$ given by (41)] associated with the Gamow state $\psi(\mathbf{X})$ is

$$ \Psi(\mathbf{X}; t) = \psi(\mathbf{X}) e^{-iEt/\hbar}. \qquad (55) $$

It satisfies the continuity equation

$$ \partial_t (|\Psi|^2) = -\nabla_X \cdot \mathbf{J} \qquad (56) $$

in terms of the probability current

$$ \mathbf{J} := \frac{\hbar}{m} \operatorname{Im} (\Psi^* \nabla_X \Psi). \qquad (57) $$

We express the loss event rate (i.e. the probability for a loss event to occur per unit of time) as

$$ \Gamma = - \left. \frac{\partial}{\partial t} \right|_{t=0} \int_{\mathcal{R}} d^{3N}X |\Psi(\mathbf{X}; t)|^2 \qquad (58) $$

where the region $\mathcal{R}$ should physically correspond to $N$ atoms in the trap, and where $\psi(\mathbf{X}) = \Psi(\mathbf{X}, t=0)$ is normalized by the condition

$$ \int_{\mathcal{R}} d^{3N}X |\psi(\mathbf{X})|^2 = 1. \qquad (59) $$

There is some freedom in how to define $\mathcal{R}$. Let $\mathcal{R}_{\text{trap}}$ denote the “trapping region” of the potential $U(\mathbf{r})$, which can be defined as the set of points $\mathbf{r}$ such that the classical trajectory of a particle with initial position $\mathbf{r}$ and zero initial velocity remains bounded.<fnref n="22" /> A possible definition of $\mathcal{R}$ would be $\mathcal{R} = (\mathcal{R}_{\text{trap}})^N$, but for later convenience, we choose $\mathcal{R}$ to be slightly smaller, by excluding configurations with two or three nearby particles:

$$ \mathcal{R} = (\mathcal{R}_{\text{trap}})^N \setminus (\mathcal{R}_2 \cup \mathcal{R}_3) $$

where

$$ \mathcal{R}_2 = \bigcup_{i:\uparrow, j:\downarrow} B_{ij} \quad \text{with} \quad B_{ij} = \{(\mathbf{r}_1, \dots, \mathbf{r}_N) \mid r_{ij} < d_2\} \qquad (60) $$

$$ \mathcal{R}_3 = \bigcup_{i<j:\uparrow, k:\downarrow \text{ or } i<j:\downarrow, k:\uparrow} B_{ijk} \quad \text{with} \quad B_{ijk} = \{(\mathbf{r}_1, \dots, \mathbf{r}_N) \mid R_{ijk} < d_3\}. \qquad (61) $$

In other words, $(\mathbf{r}_1, \dots, \mathbf{r}_N) \in \mathcal{R}$ means that

* all particle positions $\mathbf{r}_i$ are inside the trapping region $\mathcal{R}_{\text{trap}}$

* all distances between pairs of particles with opposite spin are $\ge d_2$

* all hyperradii of triplets of particles with non-identical spins are $\ge d_3$.

We take $d_2$ such that

$$ b \ll d_2 \qquad (62) $$

$$ d_2 \ll d_3. \qquad (63) $$

Conditions (54), (63) ensure that both $d_2$ and $d_3$ are $\ll 1/k_{\text{typ}}$. Hence, to leading order, the normalization integral (59) is independent of $d_2$ and $d_3$, because $\mathcal{R}_2$ and $\mathcal{R}_3$ are negligibly small subsets of $(\mathcal{R}_{\text{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 $\psi(\mathbf{X})$ for $\mathbf{X} \in \mathcal{R}$

* in (49), $B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)$ can be evaluated within the zero-range model, provided $\mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$ all belong to $\mathcal{R}_{\text{trap}}$; indeed, the $\propto a_3$ term is a small correction, and the condition (49) for the Gamow state must match with the condition (15) for the zero-range model.

<footnote n="22">²²For example, if $U(\mathbf{r}) = m\omega^2 r^2 - \lambda r^4$, then $\mathcal{R}_{\text{trap}}$ is the sphere centered at the origin of radius $r_{\text{trap}}$ such that $U(r)$ reaches its maximum at $r = r_{\text{trap}}$, i.e. $r_{\text{trap}} = \sqrt{m/\lambda} \times \omega/2$. This definition of $\mathcal{R}_{\text{trap}}$ is merely one particularly simple and natural choice among a range of possibilities.</footnote>
Equations (55), (58), (59) directly yield

$$ \Gamma = -\frac{2}{\hbar} \operatorname{Im} E, \qquad (64) $$

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

$$ \Phi(\mathcal{S}) := \int_{\mathcal{S}} \mathbf{d}^{3N-1}\mathbf{S} \cdot \mathbf{J} \qquad (65) $$

for the probability flux through a surface $\mathcal{S}$. In (58), we interchange the time derivative and the integration, and use the continuity equation (56) and the divergence theorem, which yields

$$ \Gamma = \Phi(\partial \mathcal{R}) \qquad (66) $$

where $\partial \mathcal{R}$ is the boundary of $\mathcal{R}$. Here and in what follows, the differential surface vector $\mathbf{d}^{3N-1}\mathbf{S}$ appearing in (65) is oriented towards the exterior of $\mathcal{R}$. Assuming that the trap depth is large enough for evaporation to be negligible$^{23}$, we have

$$ \Gamma = \Phi(\mathcal{S}_2) + \Phi(\mathcal{S}_3) \qquad (67) $$

where

$$ \mathcal{S}_2 = \partial \mathcal{R}_2 \cap \mathcal{R}, \quad \mathcal{S}_3 = \partial \mathcal{R}_3 \cap \mathcal{R}. $$

A visual representation is shown in Fig. 1. On physical grounds, we identify $\Phi(\mathcal{S}_2)$ as the two-body loss rate $\Gamma_2$, and $\Phi(\mathcal{S}_3)$ as the three-body loss rate $\Gamma_3$, which determine the average number of lost atoms per unit of time: $-\dot{N} = 2\Gamma_2 + 3\Gamma_3$. In the considered regime of small range and large trap depth, we expect that $\Gamma_2$ is given by the relation (47), and that in the considered case where $a_2$ is real, $\Gamma_2$ is negligible compared to $\Gamma_3$ so that $\Gamma \approx \Gamma_3$. $^{24}$

It remains to determine $\Gamma_3 = \Phi(\mathcal{S}_3)$. We will use the notation

$$ \mathcal{S}^{(ijk)} := \partial B_{ijk} \cap \mathcal{R} = \{ \mathbf{X} \in \mathcal{R} \mid R_{ijk} = d_3 \}. $$

From (61), we have $\mathcal{S}_3 = \bigcup_{i<j:\uparrow, k:\downarrow \text{ or } i<j:\downarrow, k:\uparrow} \mathcal{S}^{(ijk)}$. Hence $\Phi(\mathcal{S}_3) = \Gamma_{2,1} + \Gamma_{1,2}$ with

* $\Gamma_{2,1} = \sum_{i<j:\uparrow, k:\downarrow} \Phi(\mathcal{S}^{(ijk)})$ the contribution from $\uparrow\uparrow\downarrow$ triplets

* $\Gamma_{1,2} = \sum_{i<j:\uparrow, k:\uparrow} \Phi(\mathcal{S}^{(ijk)})$ the contribution from $\uparrow\downarrow\downarrow$ triplets.

By antisymmetry, each $\uparrow\uparrow\downarrow$ triplet gives the same contribution, so that

$$ \Gamma_{2,1} = \frac{N_{\uparrow}(N_{\downarrow} - 1)N_{\downarrow}}{2} \Phi(\mathcal{S}^{(123)}). \qquad (68) $$

<footnote>$^{23}$ Evaporation 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 trapping-potential 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 $[\operatorname{Min}_{\mathbf{r} \in \partial \mathcal{R}_{\text{trap}}} U(\mathbf{r})] - [\operatorname{Min}_{\mathbf{r} \in \mathcal{R}_{\text{trap}}} U(\mathbf{r})]$, the evaporation rate is $\Phi(\mathcal{S}_{\text{trap}})$ with $\mathcal{S}_{\text{trap}} = \partial((\mathcal{R}_{\text{trap}})^N) \cap \mathcal{R} = \{\mathbf{X} \in \mathcal{R} \mid \exists i, \mathbf{r}_i \in \partial \mathcal{R}_{\text{trap}}\}$, and this rate is exponentially suppressed because $J(\mathbf{X})$ with $\mathbf{X} \in \mathcal{S}_{\text{trap}}$ is exponentially suppressed.</footnote>

<footnote>$^{24}$ Indeed, in the zero-range limit, the reasoning of footnote 20 yields the expression (47) for $\Phi(\mathcal{S}_2) = \Gamma_2$; moreover, in the case where $a_2 \in \mathbb{R}$, we expect $\Gamma_2 \ll \Gamma_3$, because there is no mechanism that would generate a non-negligible flux exiting $\mathcal{R}$ through $\mathcal{S}_2$ (hence entering into $\mathcal{R}_2$) and propagating in $\mathcal{R}_2$ with an initial wavevector high enough to climb the trapping potential barrier and escape from $(\mathcal{R}_{\text{trap}})^N$.</footnote>
Equations (55), (58), (59) directly yield

$$ \Gamma = -\frac{2}{\hbar} \operatorname{Im} E, \qquad (64) $$

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

$$ \Phi(\mathcal{S}) := \int_{\mathcal{S}} \mathbf{d}^{3N-1}\mathbf{S} \cdot \mathbf{J} \qquad (65) $$

for the probability flux through a surface $\mathcal{S}$. In (58), we interchange the time derivative and the integration, and use the continuity equation (56) and the divergence theorem, which yields

$$ \Gamma = \Phi(\partial \mathcal{R}) \qquad (66) $$

where $\partial \mathcal{R}$ is the boundary of $\mathcal{R}$. Here and in what follows, the differential surface vector $\mathbf{d}^{3N-1}\mathbf{S}$ appearing in (65) is oriented towards the exterior of $\mathcal{R}$. Assuming that the trap depth is large enough for evaporation to be negligible$^{23}$, we have

$$ \Gamma = \Phi(\mathcal{S}_2) + \Phi(\mathcal{S}_3) \qquad (67) $$

where

$$ \mathcal{S}_2 = \partial \mathcal{R}_2 \cap \mathcal{R}, \quad \mathcal{S}_3 = \partial \mathcal{R}_3 \cap \mathcal{R}. $$

A visual representation is shown in Fig. 1. On physical grounds, we identify $\Phi(\mathcal{S}_2)$ as the two-body loss rate $\Gamma_2$, and $\Phi(\mathcal{S}_3)$ as the three-body loss rate $\Gamma_3$, which determine the average number of lost atoms per unit of time: $-\dot{N} = 2\Gamma_2 + 3\Gamma_3$. In the considered regime of small range and large trap depth, we expect that $\Gamma_2$ is given by the relation (47), and that in the considered case where $a_2$ is real, $\Gamma_2$ is negligible compared to $\Gamma_3$ so that $\Gamma \approx \Gamma_3$.<fnref n="24" />

It remains to determine $\Gamma_3 = \Phi(\mathcal{S}_3)$. We will use the notation

$$ \mathcal{S}^{(ijk)} := \partial B_{ijk} \cap \mathcal{R} = \{ \mathbf{X} \in \mathcal{R} \mid R_{ijk} = d_3 \}. $$

From (61), we have $\mathcal{S}_3 = \bigcup_{i<j:\uparrow, k:\downarrow \text{ or } i<j:\downarrow, k:\uparrow} \mathcal{S}^{(ijk)}$. Hence $\Phi(\mathcal{S}_3) = \Gamma_{2,1} + \Gamma_{1,2}$ with

* $\Gamma_{2,1} = \sum_{i<j:\uparrow, k:\downarrow} \Phi(\mathcal{S}^{(ijk)})$ the contribution from $\uparrow\uparrow\downarrow$ triplets

* $\Gamma_{1,2} = \sum_{i<j:\uparrow, k:\uparrow} \Phi(\mathcal{S}^{(ijk)})$ the contribution from $\uparrow\downarrow\downarrow$ triplets.

By antisymmetry, each $\uparrow\uparrow\downarrow$ triplet gives the same contribution, so that

$$ \Gamma_{2,1} = \frac{N_{\uparrow}(N_{\downarrow} - 1)N_{\downarrow}}{2} \Phi(\mathcal{S}^{(123)}). \qquad (68) $$

<footnote n="23">²³ Evaporation 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 trapping-potential 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 $[\operatorname{Min}_{\mathbf{r} \in \partial \mathcal{R}_{\text{trap}}} U(\mathbf{r})] - [\operatorname{Min}_{\mathbf{r} \in \mathcal{R}_{\text{trap}}} U(\mathbf{r})]$, the evaporation rate is $\Phi(\mathcal{S}_{\text{trap}})$ with $\mathcal{S}_{\text{trap}} = \partial((\mathcal{R}_{\text{trap}})^N) \cap \mathcal{R} = \{\mathbf{X} \in \mathcal{R} \mid \exists i, \mathbf{r}_i \in \partial \mathcal{R}_{\text{trap}}\}$, and this rate is exponentially suppressed because $J(\mathbf{X})$ with $\mathbf{X} \in \mathcal{S}_{\text{trap}}$ is exponentially suppressed.</footnote>

<footnote n="24">²⁴ Indeed, in the zero-range limit, the reasoning of footnote 20 yields the expression (47) for $\Phi(\mathcal{S}_2) = \Gamma_2$; moreover, in the case where $a_2 \in \mathbb{R}$, we expect $\Gamma_2 \ll \Gamma_3$, because there is no mechanism that would generate a non-negligible flux exiting $\mathcal{R}$ through $\mathcal{S}_2$ (hence entering into $\mathcal{R}_2$) and propagating in $\mathcal{R}_2$ with an initial wavevector high enough to climb the trapping potential barrier and escape from $(\mathcal{R}_{\text{trap}})^N$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Geometric illustration of the three-body loss process. The total decay rate is given by the probability flux exiting from the region $\mathcal{R}$ (grey shaded area). Neglecting the flux through $S_2$ (red straight lines) which corresponds to two-body losses, and the flux through $S_{\text{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 $\mathcal{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 $\mathbf{C}=\mathbf{0}$. The positions of the three particles are then determined by the Jacobi coordinates $\mathbf{r}$ and $\rho$. The trapping region was simply assumed to be a symmetric interval around the origin.</figcaption></figure>

This can be rewritten as<fnref n="25" />

$$ \Gamma_{2,1} = -N_{\uparrow} (N_{\downarrow} - 1) N_{\downarrow} \frac{3\sqrt{3}}{8} d_3^5 \int_{\mathcal{S}^{(123)}} d^3\mathbf{C} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N d^5\Omega \frac{\hbar}{m} \text{Im}[\psi^* \partial_R \psi]. \quad (69) $$

<footnote n="25">²⁵To justify this rewriting, we mostly follow the bosonic case treated in [72, App. B] (we will also correct in passing a minor error in an intermediate step in [72]). We need to evaluate $Q := \int_{\mathcal{S}^{(123)}} (\psi^* \nabla_X \psi - \psi \nabla_X \psi^*) \cdot d^3\mathbf{r}$. The constraint $R \equiv R_{123} > d_3$ which defines the domain $B_{123}$ does not impose any constraint on $(\mathbf{r}_4, \dots, \mathbf{r}_N)$. Therefore the differential surface vector $d^3\mathbf{r}^*$, being normal to $\partial B_{123}$, only has its 9 first coordinates which are non-zero, while its coordinates 10 to 3N are vanishing. More precisely, denoting by $d^3\mathbf{r}^*$ the 9-dimensional vector whose coordinates are equal to the 9 first coordinates of $d^3\mathbf{r}^*$, we have $d^3\mathbf{r}^* = d^8\mathbf{S} \cdot d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N$ where $d^8\mathbf{S}$ is the differential surface vector of the domain $\{(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) / R > d_3\}$. Hence $Q = \int_{\mathcal{S}^{(123)}} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \cdot d^8\mathbf{S} \cdot (\psi^* \nabla_{\mathbf{X}} \psi - \psi \nabla_{\mathbf{X}} \psi^*)$ where $\mathbf{X} := (\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$. Applying the divergence theorem backwards, this can be rewritten $Q = \int_{\mathcal{S}^{(123)}} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \cdot d^9\mathbf{X} \cdot \nabla_{\mathbf{X}} \cdot (\psi^* \nabla_{\mathbf{X}} \psi - \psi \nabla_{\mathbf{X}} \psi^*)$. We then perform the change of variables $\mathbf{X} \to (\mathbf{C}, \mathbf{R})$, and rewrite the integrand as $\psi^* \Delta_{\mathbf{X}} \psi$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geometric illustration of the three-body loss process. The total decay rate is given by the probability flux exiting from the region $\mathcal{R}$ (grey shaded area). Neglecting the flux through $S_2$ (red straight lines) which corresponds to two-body losses, and the flux through $S_{\text{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 $\mathcal{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 $\mathbf{C}=\mathbf{0}$. The positions of the three particles are then determined by the Jacobi coordinates $\mathbf{r}$ and $\rho$. The trapping region was simply assumed to be a symmetric interval around the origin.</figcaption></figure>

This can be rewritten as<fnref n="25" />

$$ \Gamma_{2,1} = -N_{\uparrow} (N_{\downarrow} - 1) N_{\downarrow} \frac{3\sqrt{3}}{8} d_3^5 \int_{\mathcal{S}^{(123)}} d^3\mathbf{C} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N d^5\Omega \frac{\hbar}{m} \text{Im}[\psi^* \partial_R \psi]. \quad (69) $$

<footnote n="25">²⁵To justify this rewriting, we mostly follow the bosonic case treated in [72, App. B] (we will also correct in passing a minor error in an intermediate step in [72]). We need to evaluate $Q := \int_{\mathcal{S}^{(123)}} (\psi^* \nabla_X \psi - \psi \nabla_X \psi^*) \cdot d^3\mathbf{r}$. The constraint $R \equiv R_{123} > d_3$ which defines the domain $B_{123}$ does not impose any constraint on $(\mathbf{r}_4, \dots, \mathbf{r}_N)$. Therefore the differential surface vector $d^3\mathbf{r}^*$, being normal to $\partial B_{123}$, only has its 9 first coordinates which are non-zero, while its coordinates 10 to 3N are vanishing. More precisely, denoting by $d^3\mathbf{r}^*$ the 9-dimensional vector whose coordinates are equal to the 9 first coordinates of $d^3\mathbf{r}^*$, we have $d^3\mathbf{r}^* = d^8\mathbf{S} \cdot d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N$ where $d^8\mathbf{S}$ is the differential surface vector of the domain $\{(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) / R > d_3\}$. Hence $Q = \int_{\mathcal{S}^{(123)}} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \cdot d^8\mathbf{S} \cdot (\psi^* \nabla_{\mathbf{r}^*} \psi - \psi \nabla_{\mathbf{r}^*} \psi^*)$ where $\mathbf{r}^* := (\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$. Applying the divergence theorem backwards, this can be rewritten $Q = \int_{\mathcal{S}^{(123)}} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N \cdot d^9\mathbf{X} \cdot \nabla_{\mathbf{X}} \cdot (\psi^* \nabla_{\mathbf{X}} \psi - \psi \nabla_{\mathbf{X}} \psi^*)$. We then perform the change of variables $\mathbf{X} \to (\mathbf{C}, \mathbf{R})$, and rewrite the integrand as $\psi^* \Delta_{\mathbf{X}} \psi$.</footnote>
We then simplify this expression in two steps:

(i) replace $\psi$ with its asymptotic behavior (49), with $B_m$ evaluated within the zero-range model,

(ii) replace the integration domain $\mathcal{S}^{(123)}$ by the entire region $\partial B_{123}$.

Step (i) is justified since the conditions (50), (51), (52) hold except in a negligibly small domain of the integration variables $C, r_4, \dots, r_N, \Omega$. Step (ii) is justified because given (62), the condition $X \notin \mathcal{R}_2$ only excludes a negligibly small region of hyperangles $\Omega$. We note that the order of these two steps is important.<fnref n="26" /> Using (17) then yields the result (44). The expression (45) of $\Gamma_{1,2}$ is derived analogously.

Let us now discuss in more detail the appropriate choice of the length $d_3$. The condition (53) is actually not sufficient in order to have the behavior (49) of $\psi$. Indeed, we expect [based on the small-$R$ expansion of the finite-energy solution $J_s(kR)$ of the hyperradial Schrödinger equation (94)] that in addition to the term $R^s$ in (49), there is a higher-order correction term of order $R^{s+2} k_{\text{typ}}^2$, which is negligible compared to $a_3/R^s$ provided we take

$$d_3 \ll b/(k_{\text{typ}} b)^{1/(s+1)}. \qquad (70)$$

The condition $d_3 \ll 1/k_{\text{typ}}$ is then automatically satisfied. Compatibility with (53) then requires $(k_{\text{typ}} b)^{1/(s+1)} \ll 1$, which is a quite stringent condition given the smallness of the exponent $1/(s+1) \approx 0.36$. However, while this condition is necessary for (49), we do not expect it to be necessary for the final expression (43) 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 (49) 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'\pi/2) + 4 \sin(s'\pi/6) / \sqrt{3} = 0$, and by $\phi'(\Omega) := \phi(l=0,n=0)(\Omega)$ the corresponding $l=0$ hyperangular wavefunction. There is a higher-order correction to the r.h.s. of (49) given by $R^{s'-2} \phi'(\Omega)$ times a function $B'(C; r_4, \dots, r_N)$. Requiring this $\propto R^{s'}$ correction to be negligible compared to the $\propto a_3/R^s$ term in (49) would yield an additional condition on $d_3$, but this is not necessary for the final result (43). Instead, what is truly necessary for (43) is that the contribution from the $l=0$ angular-momentum sector to $\Gamma_3$ is negligible compared to the $l=1$ contribution. Adding the term $(R^{s'} - a_3'/R^{s'}) R^{-2} \phi'(\Omega) B'(C; r_4, \dots, r_N)$ to the r.h.s. of (49), and still using (69), we obtain the additional term $\Gamma_3' = -(\hbar/m) 4s'(s'+1) C_3' \text{Im} a_3'$ on the r.h.s. of (43). Here $C_3'$ is the “$l=0$ three-body contact”, defined in such a way that there is a higher-order correction $C_3' e^{2s'+2}$ to the r.h.s. of (3). The ratio $\Gamma_3'/\Gamma_3 \propto b^{2(s'-s)}$ is small in the zero-range regime (except if $C_3$ is anomalously small). For example, for the degenerate unpolarized unitary gas, this yields the truly necessary condition $(n^{1/3} b)^{2(s'-s)} \ll 1$ for the validity of (43).

### 3.1.2. 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,

<footnote>$\psi \Delta_X \psi^* = \frac{1}{3}(\psi^* \Delta_C \psi - \psi \Delta_C \psi^*) + 2(\psi^* \Delta_R \psi - \psi \Delta_R \psi^*) = \frac{1}{3} \nabla_C \cdot (\psi^* \nabla_C \psi - \psi \nabla_C \psi^*) + 2 \nabla_R \cdot (\psi^* \nabla_R \psi - \psi \nabla_R \psi^*)$. By the divergence theorem, the integral over $C$ of the $\nabla_C$ term vanishes, while the integral over $R$ of the $\nabla_R$ term yields the result $Q = -\frac{3\sqrt{3}}{4} d_3^5 \int_{\mathcal{S}^{(123)}} d^3C d^3r_4 \dots d^3r_N d^5\Omega (\psi^* \partial_R \psi - \psi \partial_R \psi^*)$.</footnote>

<footnote n="26">²⁶If we would start by replacing $\mathcal{S}^{(123)}$ with $\partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N$ (keeping the original Gamow-state wavefunction $\psi$), then we would get a vanishing result for the total flux : The flux through $\partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N \setminus \mathcal{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 $\mathcal{R}_2 \cap \partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N$ where the deep-dimer + atom outgoing-wave contribution to $\psi$ gives the main contribution to the flux. This follows from the fact that for the three-body zero-energy scattering state $\Psi_m$, the total flux $\int d^5\Omega \text{Im}(\Psi_m^* \partial_R \Psi_m)$ vanishes (a related discussion can be found in [154]).</footnote>
We then simplify this expression in two steps:

(i) replace $\psi$ with its asymptotic behavior (49), with $B_m$ evaluated within the zero-range model,

(ii) replace the integration domain $\mathcal{S}^{(123)}$ by the entire region $\partial B_{123}$.

Step (i) is justified since the conditions (50), (51), (52) hold except in a negligibly small domain of the integration variables $C, r_4, \dots, r_N, \Omega$. Step (ii) is justified because given (62), the condition $X \notin \mathcal{R}_2$ only excludes a negligibly small region of hyperangles $\Omega$. We note that the order of these two steps is important.<fnref n="26" /> Using (17) then yields the result (44). The expression (45) of $\Gamma_{1,2}$ is derived analogously.

Let us now discuss in more detail the appropriate choice of the length $d_3$. The condition (53) is actually not sufficient in order to have the behavior (49) of $\psi$. Indeed, we expect [based on the small-$R$ expansion of the finite-energy solution $J_s(kR)$ of the hyperradial Schrödinger equation (94)] that in addition to the term $R^s$ in (49), there is a higher-order correction term of order $R^{s+2} k_{\text{typ}}^2$, which is negligible compared to $a_3/R^s$ provided we take

$$d_3 \ll b/(k_{\text{typ}} b)^{1/(s+1)}. \quad (70)$$

The condition $d_3 \ll 1/k_{\text{typ}}$ is then automatically satisfied. Compatibility with (53) then requires $(k_{\text{typ}} b)^{1/(s+1)} \ll 1$, which is a quite stringent condition given the smallness of the exponent $1/(s+1) \approx 0.36$. However, while this condition is necessary for (49), we do not expect it to be necessary for the final expression (43) 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 (49) 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'\pi/2) + 4 \sin(s'\pi/6)/\sqrt{3} = 0$, and by $\phi'(\Omega) := \phi(l=0,n=0)(\Omega)$ the corresponding $l=0$ hyperangular wavefunction. There is a higher-order correction to the r.h.s. of (49) given by $R^{s'-2}\phi'(\Omega)$ times a function $B'(C; r_4, \dots, r_N)$. Requiring this $\propto R^{s'}$ correction to be negligible compared to the $\propto a_3/R^s$ term in (49) would yield an additional condition on $d_3$, but this is not necessary for the final result (43). Instead, what is truly necessary for (43) is that the contribution from the $l=0$ angular-momentum sector to $\Gamma_3$ is negligible compared to the $l=1$ contribution. Adding the term $(R^{s'} - a_3'/R^{s'}) R^{-2}\phi'(\Omega) B'(C; r_4, \dots, r_N)$ to the r.h.s. of (49), and still using (69), we obtain the additional term $\Gamma_3' = -(\hbar/m)4s'(s'+1)C_3' \text{Im} a_3'$ on the r.h.s. of (43). Here $C_3'$ is the “$l=0$ three-body contact”, defined in such a way that there is a higher-order correction $C_3' e^{2s'+2}$ to the r.h.s. of (3). The ratio $\Gamma_3'/\Gamma_3 \propto b^{2(s'-s)}$ is small in the zero-range regime (except if $C_3$ is anomalously small). For example, for the degenerate unpolarized unitary gas, this yields the truly necessary condition $(n^{1/3}b)^{2(s'-s)} \ll 1$ for the validity of (43).

### 3.1.2. 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,

<footnote>$\psi \Delta_X \psi^* = \frac{1}{3}(\psi^* \Delta_C \psi - \psi \Delta_C \psi^*) + 2(\psi^* \Delta_R \psi - \psi \Delta_R \psi^*) = \frac{1}{3} \nabla_C \cdot (\psi^* \nabla_C \psi - \psi \nabla_C \psi^*) + 2 \nabla_R \cdot (\psi^* \nabla_R \psi - \psi \nabla_R \psi^*)$. By the divergence theorem, the integral over $C$ of the $\nabla_C$ term vanishes, while the integral over $R$ of the $\nabla_R$ term yields the result $Q = -\frac{3\sqrt{3}}{4} d_3^5 \int_{\mathcal{S}^{(123)}} d^3C d^3r_4 \dots d^3r_N d^5\Omega (\psi^* \partial_R \psi - \psi \partial_R \psi^*)$.</footnote>

<footnote n="26">²⁶If we would start by replacing $\mathcal{S}^{(123)}$ with $\partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N$ (keeping the original Gamow-state wavefunction $\psi$), then we would get a vanishing result for the total flux : The flux through $\partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N \setminus \mathcal{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 $\mathcal{R}_2 \cap \partial B_{123} \cap (\mathcal{R}_{\text{trap}})^N$ where the deep-dimer + atom outgoing-wave contribution to $\psi$ gives the main contribution to the flux. This follows from the fact that for the three-body zero-energy scattering state $\Psi_m$, the total flux $\int d^5\Omega \text{Im}(\Psi_m^* \partial_R \Psi_m)$ vanishes (a related discussion can be found in [154]).</footnote>
* 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(\mathbf{r})$ which may now depend on the direction of $\mathbf{r}$, and a three-body interaction potential $V_{2,1}$ (resp. $V_{1,2}$) between triplets of particles of spins $\uparrow\uparrow\downarrow$ (resp. $\uparrow\downarrow\downarrow$).

The corresponding stationary N-body Schrödinger equation is

$$ \sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \psi + \sum_{i:\uparrow, j:\downarrow} V_2(\mathbf{r}_{ij}) \psi + \sum_{i<j:\uparrow, k:\downarrow} V_{2,1}(\mathbf{R}_{ijk}) \psi + \sum_{i<j:\uparrow, k:\downarrow} V_{1,2}(\mathbf{R}_{ijk}) \psi = E \psi. \quad (71) $$

Here $\mathbf{R}_{ijk}$ denotes the Jacobi coordinates associated to particles $i, j, k$ [defined by replacing the indices 1,2,3 by $i, j, k$ in (89), (90)].

The two-body interaction $V_2$ is still assumed to be resonant, $|a_2| \gg 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 where

$$ b := \text{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 $\psi$ in the region where all interparticle distances are $\gg b$. For alkali atoms near an open-channel dominated Feshbach resonance, $b_2$ is set by the van der Waals length [159], which is $\ll 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_{2,1}^{(m)}$ and $a_{1,2}^{(m)}$, with $m \in \{-1, 0, 1\}$ the angular momentum quantum number.

We find that the $\uparrow\uparrow\downarrow$ and $\uparrow\downarrow\downarrow$ three-body loss-rates are given by

$$ \Gamma_{2,1} \approx -\frac{\hbar}{m} 8s(s+1) \sum_{m=-1}^{1} C_{2,1}^{(m)} \operatorname{Im} a_{2,1}^{(m)} \quad (72) $$

$$ \Gamma_{1,2} \approx -\frac{\hbar}{m} 8s(s+1) \sum_{m=-1}^{1} C_{1,2}^{(m)} \operatorname{Im} a_{1,2}^{(m)}. \quad (73) $$

Here, the m-resolved three-body contacts $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are defined by

$$ C_{2,1}^{(m)} := N_{\uparrow}(N_{\downarrow}-1)N_{\downarrow} \frac{3\sqrt{3}}{32(s+1)} \int |B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3C d^3r_4 \dots d^3r_N \quad (74) $$

where $B_m$ is related to the many-body wavefunction (in the zero-range limit) through Eq. (15), and similarly

$$ C_{1,2}^{(m)} := N_{\downarrow}(N_{\uparrow}-1)N_{\uparrow} \frac{3\sqrt{3}}{32(s+1)} \int |\tilde{B}_m(\mathbf{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N)|^2 d^3C d^3r_2 d^3r_5 \dots d^3r_N \quad (75) $$

where $\tilde{B}_m$ is defined in Eq. (18).

Note that from Eqs. (17), (19) we have $\sum_{m=-1}^{1} C_{2,1}^{(m)} = C_{2,1}$ and $\sum_{m=-1}^{1} C_{1,2}^{(m)} = C_{1,2}$.

The three-body parameters $a_{2,1}^{(m)}$ are defined as follows. Setting

$$ \mathcal{V}(\mathbf{R}) := V_2(\mathbf{r}_{13}) + V_2(\mathbf{r}_{23}) + V_{2,1}(\mathbf{R}), \quad (76) $$

the solution $\Psi_m$ of the zero-energy Schrödinger equation in free space for two $\uparrow$ and one $\downarrow$ particles

$$ -\frac{\hbar^2}{m} \Delta_{\mathbf{R}} \Psi_m + \mathcal{V}(\mathbf{R}) \Psi_m = 0 \quad (77) $$
* 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(\mathbf{r})$ which may now depend on the direction of $\mathbf{r}$, and a three-body interaction potential $V_{2,1}$ (resp. $V_{1,2}$) between triplets of particles of spins $\uparrow\uparrow\downarrow$ (resp. $\uparrow\downarrow\downarrow$).

The corresponding stationary N-body Schrödinger equation is

$$ \sum_{i=1}^{N} \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \psi + \sum_{i:\uparrow, j:\downarrow} V_2(\mathbf{r}_{ij}) \psi + \sum_{i<j:\uparrow, k:\downarrow} V_{2,1}(\mathbf{R}_{ijk}) \psi + \sum_{i<j:\uparrow, k:\downarrow} V_{1,2}(\mathbf{R}_{ijk}) \psi = E \psi. \quad (71) $$

Here $\mathbf{R}_{ijk}$ denotes the Jacobi coordinates associated to particles $i, j, k$ [defined by replacing the indices 1,2,3 by $i, j, k$ in (89), (90)].

The two-body interaction $V_2$ is still assumed to be resonant, $|a_2| \gg 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 where

$$ b := \text{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 $\psi$ in the region where all interparticle distances are $\gg b$. For alkali atoms near an open-channel dominated Feshbach resonance, $b_2$ is set by the van der Waals length [159], which is $\ll 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_{2,1}^{(m)}$ and $a_{1,2}^{(m)}$, with $m \in \{-1, 0, 1\}$ the angular momentum quantum number.

We find that the $\uparrow\uparrow\downarrow$ and $\uparrow\downarrow\downarrow$ three-body loss-rates are given by

$$ \Gamma_{2,1} \approx -\frac{\hbar}{m} 8s(s+1) \sum_{m=-1}^{1} C_{2,1}^{(m)} \text{Im} a_{2,1}^{(m)} \quad (72) $$

$$ \Gamma_{1,2} \approx -\frac{\hbar}{m} 8s(s+1) \sum_{m=-1}^{1} C_{1,2}^{(m)} \text{Im} a_{1,2}^{(m)}. \quad (73) $$

Here, the m-resolved three-body contacts $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are defined by

$$ C_{2,1}^{(m)} := N_{\uparrow}(N_{\downarrow}-1)N_{\downarrow} \frac{3\sqrt{3}}{32(s+1)} \int |B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N)|^2 d^3C d^3r_4 \dots d^3r_N \quad (74) $$

where $B_m$ is related to the many-body wavefunction (in the zero-range limit) through Eq. (15), and similarly

$$ C_{1,2}^{(m)} := N_{\downarrow}(N_{\uparrow}-1)N_{\uparrow} \frac{3\sqrt{3}}{32(s+1)} \int |\tilde{B}_m(\mathbf{C}; \mathbf{r}_2, \mathbf{r}_5, \dots, \mathbf{r}_N)|^2 d^3C d^3r_2 d^3r_5 \dots d^3r_N \quad (75) $$

where $\tilde{B}_m$ is defined in Eq. (18).

Note that from Eqs. (17), (19) we have $\sum_{m=-1}^{1} C_{2,1}^{(m)} = C_{2,1}$ and $\sum_{m=-1}^{1} C_{1,2}^{(m)} = C_{1,2}$.

The three-body parameters $a_{2,1}^{(m)}$ are defined as follows. Setting

$$ \mathcal{V}(\mathbf{R}) := V_2(\mathbf{r}_{13}) + V_2(\mathbf{r}_{23}) + V_{2,1}(\mathbf{R}), \quad (76) $$

the solution $\Psi_m$ of the zero-energy Schrödinger equation in free space for two $\uparrow$ and one $\downarrow$ particles

$$ -\frac{\hbar^2}{m} \Delta_{\mathbf{R}} \Psi_m + \mathcal{V}(\mathbf{R}) \Psi_m = 0 \quad (77) $$
with angular-momentum quantum numbers ($l=1, m$) has the asymptotic behavior

$$ \Psi_m(\mathbf{R}) \approx \left( R^s - \frac{a_{2,1}^{(m)}}{R^s} \right) \frac{1}{R^2} \phi_m(\Omega) \quad (78) $$

in the region where all interparticle distances are $\gg b$ and $\ll |a_2|$. The three-body parameters $a_{1,2}^{(m)}$ are defined similarly in terms of the $\uparrow\downarrow\downarrow$ three-body scattering states.

The relation (72) is derived by considering the probability current, in a completely analogous way to Eqs. (58), (66), (68), (69) above, using the asymptotic behavior of the many-body wavefunction in the region ((50), (51), (52)) which is now given by

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \frac{1}{R^2} \sum_{m=-1}^{+1} \left( R^s - \frac{a_{2,1}^{(m)}}{R^s} \right) \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (79) $$

The expression (73) of $\Gamma_{1,2}$ is derived analogously.

**Discussion:** The parameters $Im\,a_{2,1}^{(m)}$ and $Im\,a_{1,2}^{(m)}$ are a priori unknown.<fnref n="27" /> In principle, one may hope to compute them theoretically by solving a sufficiently realistic three-body problem, but this is a difficult task [163–165].<fnref n="28" /> Instead, one could determine them by measuring the three-body loss rate in situations where the three-body contacts $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are known theoretically.

A first possibility is to work with a small number $N$ of particles in a microtrap [166–169] where the $N$-body wavefunction can be computed numerically with good accuracy [106, 170, 171] so that the three-body contacts $C_{i,j}^{(m)}$ could be calculated. Measuring $\Gamma_3$ in six different states and inverting (72), (73) would allow to determine the six parameters $Im\,a_{i,j}^{(m)}$. 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 \in \{-1; 0; 1\}$ and $(N_\uparrow, N_\downarrow) \in \{(2, 1); (1, 2)\}$, then $C_{N_\uparrow, N_\downarrow}^{(m)}$ is the only non-zero three-body contact, so that $\Gamma_3$ is simply proportional to $Im\,a_{N_\uparrow, N_\downarrow}^{(m)}$. Explicitly, taking for example $N_\uparrow = 2$ and $N_\downarrow = 1$, we get $C_3 = C_{2,1}^{(m)} = (m\omega/\hbar)^{s+1} / [2^{s+1}\Gamma(s+2)]$, hence $\Gamma_3/\omega = -Im\,a_{2,1}^{(m)} (m\omega/\hbar)^s \times 1.266322$ (in agreement with the scaling given in [172]).<fnref n="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_{i,j}^{(m)}$ are all equal, as shown in App. B. This gives

$$ \Gamma_3 \approx -\frac{\hbar}{m} 8s(s+1) C_3 \operatorname{Im} \bar{a}_3 \quad (80) $$

$$ \begin{pmatrix} q+s \\ s \end{pmatrix} = \frac{\Gamma(q+s+1)}{\Gamma(s+1)q!} $$

<footnote n="27">²⁷One may expect a small relative difference between $a_{i,j}^{(0)}$ and $a_{i,j}^{(\pm 1)}$, and an even smaller one between $a_{i,j}^{(1)}$ and $a_{i,j}^{(-1)}$, similarly to the $m$-dependence of the two-body p-wave scattering volume not too close to a p-wave Feshbach resonance [160–162].</footnote>

<footnote n="28">²⁸For 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 $|v\rangle$. However, we expect the general relations (72), (73) to remain applicable. Indeed, we expect that in typical experiments, the atoms mainly occupy two internal states, that we can label $v = \uparrow$ and $\downarrow$, and if all interatomic distances $r_{ij} \gg b$, the many-body wavefunction $\Phi(\mathbf{r}_1, v_1; \dots; \mathbf{r}_N, v_N)$ is non-negligible only if all $v_i$ belong to $\{|\uparrow\rangle; \downarrow\rangle\}$, in which case $\Phi$ is given to good accuracy by antisymmetrizing the wavefunction $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ of the zero-range model.</footnote>

<footnote n="29">²⁹For the excited state whose energy is $2q\hbar\omega$ above the ground state energy, the value of $C_{2,1}^{(m)}$ is multiplied by</footnote>

<footnote>(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).</footnote>
with angular-momentum quantum numbers ($l=1, m$) has the asymptotic behavior

$$ \Psi_m(\mathbf{R}) \approx \left( R^s - \frac{a_{2,1}^{(m)}}{R^s} \right) \frac{1}{R^2} \phi_m(\Omega) \quad (78) $$

in the region where all interparticle distances are $\gg b$ and $\ll |a_2|$. The three-body parameters $a_{1,2}^{(m)}$ are defined similarly in terms of the $\uparrow\downarrow\downarrow$ three-body scattering states.

The relation (72) is derived by considering the probability current, in a completely analogous way to Eqs. (58), (66), (68), (69) above, using the asymptotic behavior of the many-body wavefunction in the region ((50), (51), (52)) which is now given by

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \frac{1}{R^2} \sum_{m=-1}^{+1} \left( R^s - \frac{a_{2,1}^{(m)}}{R^s} \right) \phi_m(\Omega) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N). \quad (79) $$

The expression (73) of $\Gamma_{1,2}$ is derived analogously.

**Discussion:** The parameters $Im\,a_{2,1}^{(m)}$ and $Im\,a_{1,2}^{(m)}$ are a priori unknown.<fnref n="27" /> In principle, one may hope to compute them theoretically by solving a sufficiently realistic three-body problem, but this is a difficult task [163–165].<fnref n="28" /> Instead, one could determine them by measuring the three-body loss rate in situations where the three-body contacts $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are known theoretically.

A first possibility is to work with a small number $N$ of particles in a microtrap [166–169] where the $N$-body wavefunction can be computed numerically with good accuracy [106, 170, 171] so that the three-body contacts $C_{i,j}^{(m)}$ could be calculated. Measuring $\Gamma_3$ in six different states and inverting (72), (73) would allow to determine the six parameters $Im\,a_{i,j}^{(m)}$. 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 \in \{-1; 0; 1\}$ and $(N_\uparrow, N_\downarrow) \in \{(2, 1); (1, 2)\}$, then $C_{N_\uparrow, N_\downarrow}^{(m)}$ is the only non-zero three-body contact, so that $\Gamma_3$ is simply proportional to $Im\,a_{N_\uparrow, N_\downarrow}^{(m)}$. Explicitly, taking for example $N_\uparrow = 2$ and $N_\downarrow = 1$, we get $C_3 = C_{2,1}^{(m)} = (m\omega/\hbar)^{s+1} / [2^{s+1}\Gamma(s+2)]$, hence $\Gamma_3/\omega = -Im\,a_{2,1}^{(m)} (m\omega/\hbar)^s \times 1.266322$ (in agreement with the scaling given in [172]).<fnref n="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_{i,j}^{(m)}$ are all equal, as shown in App. B. This gives

$$ \Gamma_3 \approx -\frac{\hbar}{m} 8s(s+1) C_3 \operatorname{Im} \bar{a}_3 \quad (80) $$

$$ \begin{pmatrix} q+s \\ s \end{pmatrix} = \frac{\Gamma(q+s+1)}{\Gamma(s+1)q!} $$

<footnote n="27">²⁷One may expect a small relative difference between $a_{i,j}^{(0)}$ and $a_{i,j}^{(\pm 1)}$, and an even smaller one between $a_{i,j}^{(1)}$ and $a_{i,j}^{(-1)}$, similarly to the $m$-dependence of the two-body $p$-wave scattering volume not too close to a $p$-wave Feshbach resonance [160–162].</footnote>

<footnote n="28">²⁸For 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 $|v\rangle$. However, we expect the general relations (72), (73) to remain applicable. Indeed, we expect that in typical experiments, the atoms mainly occupy two internal states, that we can label $v = \uparrow$ and $\downarrow$, and if all interatomic distances $r_{ij} \gg b$, the many-body wavefunction $\Phi(\mathbf{r}_1, v_1; \dots; \mathbf{r}_N, v_N)$ is non-negligible only if all $v_i$ belong to $\{|\uparrow\rangle; |\downarrow\rangle\}$, in which case $\Phi$ is given to good accuracy by antisymmetrizing the wavefunction $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ of the zero-range model.</footnote>

<footnote n="29">²⁹For the excited state whose energy is $2q\hbar\omega$ above the ground state energy, the value of $C_{2,1}^{(m)}$ is multiplied by</footnote>

<footnote>(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).</footnote>
with

$$ \bar{a}_3 := \frac{1}{6} \sum_{m=-1}^{1} \left[ a_{2,1}^{(m)} + a_{1,2}^{(m)} \right]. $$

One could then determine Im $\bar{a}_3$ by measuring $\Gamma_3$ in a weakly correlated regime, where the asymptotic behavior of the three-body contact density $\mathscr{C}_3$ can be computed exactly. A first option is the non-degenerate regime, where we have computed $\mathscr{C}_3$ for negative or infinite scattering length.<fnref n="30" /> Other options are the weakly interacting regimes where $k_{typ}|a_2|$ is small.<fnref n="31" />

### 3.2. 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 N-body Schrödinger equation is again given by (71). 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 := \text{Max}(b_2, b_3) \ll |a_2|, 1/k_{typ}$. In particular, each eigenenergy $E$ of (71) 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,

$$ \delta E := E - E_{ZRM}. $$

We find that in the zero-range regime,

$$ \delta E \approx \delta E_2 + \delta E_3 \qquad (81) $$

plus higher-order corrections, where $\delta E_2$ is given to leading order by [56]

$$ \delta E_2 = 2\pi r_e N_1 N_1 \int d^3c \left( \prod_{i \neq 1,3} d^3r_i \right) A^*(c, r_2, r_4, \dots, r_N) \left[ E + \frac{\hbar^2}{4m} \Delta_c - 2U(c) + \sum_{i \neq 1,3} \left( \frac{\hbar^2}{2m} \Delta_{r_i} - U(r_i) \right) \right] A(c, r_2, r_4, \dots, r_N) \quad (82) $$

with $r_e$ the effective range of the two-body interaction $V_2$, while $\delta E_3$ is given to leading order by

$$ \delta E_3 \approx \frac{\hbar^2}{m} 4s(s+1) \sum_{m=-1}^{1} \left[ C_{2,1}^{(m)} a_{2,1}^{(m)} + C_{1,2}^{(m)} a_{1,2}^{(m)} \right] \qquad (83) $$

in terms of the three-body contacts and three-body parameters. If the two-body and three-body 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

$$ \delta E_3 = \frac{\hbar^2}{m} 4s(s+1) C_3 a_3. \qquad (84) $$

### Remarks 4.

* $\delta E_2$ (resp. $\delta E_3$) comes from configurations where 2 (resp. 3) particles are close to each other.

<footnote n="30">³⁰For the homogeneous unpolarized unitary gas of density $n$, we obtain $\mathscr{C}_3 \sim n^3 (\frac{\hbar^2}{mk_B T})^{2-s} \frac{\sqrt{243} \pi^3}{2^{2s+1} \Gamma(s+2)}$ in the non-degenerate limit [X. Leyronas and F. Werner, "Three-body contact for fermions. II. Non-degenerate limit", to be submitted].</footnote>

<footnote n="31">³¹Although 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.</footnote>
with

$$ \bar{a}_3 := \frac{1}{6} \sum_{m=-1}^{1} \left[ a_{2,1}^{(m)} + a_{1,2}^{(m)} \right]. $$

One could then determine Im $\bar{a}_3$ by measuring $\Gamma_3$ in a weakly correlated regime, where the asymptotic behavior of the three-body contact density $\mathscr{C}_3$ can be computed exactly. A first option is the non-degenerate regime, where we have computed $\mathscr{C}_3$ for negative or infinite scattering length.<fnref n="30" /> Other options are the weakly interacting regimes where $k_{typ}|a_2|$ is small.<fnref n="31" />

### 3.2. 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 N-body Schrödinger equation is again given by (71). 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 := \text{Max}(b_2, b_3) \ll |a_2|, 1/k_{typ}$. In particular, each eigenenergy $E$ of (71) 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,

$$ \delta E := E - E_{ZRM}. $$

We find that in the zero-range regime,

$$ \delta E \approx \delta E_2 + \delta E_3 \qquad (81) $$

plus higher-order corrections, where $\delta E_2$ is given to leading order by [56]

$$ \delta E_2 = 2\pi r_e N_1 N_1 \int d^3c \left( \prod_{i \neq 1,3} d^3r_i \right) A^*(c, r_2, r_4, \dots, r_N) \left[ E + \frac{\hbar^2}{4m} \Delta_c - 2U(c) + \sum_{i \neq 1,3} \left( \frac{\hbar^2}{2m} \Delta_{r_i} - U(r_i) \right) \right] A(c, r_2, r_4, \dots, r_N) \quad (82) $$

with $r_e$ the effective range of the two-body interaction $V_2$, while $\delta E_3$ is given to leading order by

$$ \delta E_3 \approx \frac{\hbar^2}{m} 4s(s+1) \sum_{m=-1}^{1} \left[ C_{2,1}^{(m)} a_{2,1}^{(m)} + C_{1,2}^{(m)} a_{1,2}^{(m)} \right] \qquad (83) $$

in terms of the three-body contacts and three-body parameters. If the two-body and three-body 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

$$ \delta E_3 = \frac{\hbar^2}{m} 4s(s+1) C_3 a_3. \qquad (84) $$

### Remarks 4.

* $\delta E_2$ (resp. $\delta E_3$) comes from configurations where 2 (resp. 3) particles are close to each other.

<footnote n="30">³⁰For the homogeneous unpolarized unitary gas of density $n$, we obtain $\mathscr{C}_3 \sim n^3 (\frac{\hbar^2}{mk_B T})^{2-s} \frac{\sqrt{243} \pi^3}{2^{2s+1} \Gamma(s+2)}$ in the non-degenerate limit [X. Leyronas and F. Werner, "Three-body contact for fermions. II. Non-degenerate limit", to be submitted].</footnote>

<footnote n="31">³¹Although 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.</footnote>
• The relations (83), (84) are reminiscent of the relation [36]

$$ \frac{dE}{d(-1/a_2)} = \frac{\hbar^2}{4\pi m} C_2 \qquad (85) $$

which holds within the zero-range model.<fnref n="32" />

• Typically, $a_3$ is of order $b^{2s}$, which gives $\delta E_3 \sim 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].

• If $b_3 \lesssim b_2$, then $b_2 \sim b$, and since $\delta E_2$ is $\propto r_e$, which is typically of order $b_2$, we have $|\delta E_3 / \delta E_2| \propto b^{2s-1}$. Since $s > 1/2$, we get $|\delta E_3 / \delta E_2| \ll 1$ in the zero-range regime, i.e. the three-body correction to the energy is of higher order than the two-body correction.<fnref n="33" />,<fnref n="34" />

• Let us assume that (84) can be analytically continued to complex $a_3$, with $C_3$ still evaluated within the zero-range model. We then recover the expression (43) of the three-body loss rate, simply by substituting $ImE = Im \delta E_3$ into (64), and using the fact that $\Gamma = \Gamma_3$ and $C_3 \in \mathbb{R}$. Similarly, applying this procedure to (83) yields $\Gamma_3$ in agreement with the sum of Eqs. ((72), (73)).

To derive (83), we consider a stationary state $\psi^{(0)}$ and the associated eigenenergy $E^{(0)}$ of the zero-range model, and the stationary state $\psi^{(1)}$ of the finite-range model whose energy $E^{(1)}$ is close to $E^{(0)}$ in the zero-range regime. We are interested in $\delta E = E^{(1)} - E^{(0)}$.

Using the shorthand notations

$$ \begin{aligned} V(\mathbf{r}_1, \dots, \mathbf{r}_N) &:= \sum_{i:\uparrow} V_2(\mathbf{r}_{i\uparrow}) + \sum_{i<j:\uparrow, k:\downarrow} V_{2,1}(\mathbf{r}_{ijk}) + \sum_{i<j:\uparrow, k:\uparrow} V_{1,2}(\mathbf{r}_{ijk}) \\ H_0 &:= \sum_{i=1}^N \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \end{aligned} $$

we have $(H_0 + V)\psi^{(1)} = E^{(1)}\psi^{(1)}$, while $\psi^{(0)}$ satisfies $H_0\psi^{(0)} = E^{(0)}\psi^{(0)}$ together with the two-body contact condition (14). Let us consider

$$ \Delta := \langle \psi^{(0)}, (H_0 + V) \psi^{(1)} \rangle - \langle H_0 \psi^{(0)}, \psi^{(1)} \rangle. $$

We have $\Delta = \delta E \langle \psi^{(0)} | \psi^{(1)} \rangle$, and hence $\delta E / \Delta \to 1$ in the zero-range limit.

To evaluate $\Delta$, we write it as $\Delta = \Delta_0 + \langle \psi^{(0)} | V | \psi^{(1)} \rangle$ where $\Delta_0 := \langle \psi^{(0)}, H_0 \psi^{(1)} \rangle - \langle H_0 \psi^{(0)}, \psi^{(1)} \rangle$. We have $\Delta_0 = 0$, assuming for simplicity that the interaction potentials ($V_{2,1}$ and $V_{1,2}$) are finite everywhere (i.e. excluding hard-wall potentials) and do not diverge too quickly at short distances.<fnref n="35" /> Hence

$$ \Delta \approx \int d^3 r_1 \dots d^3 r_N (\psi^{(0)*} V \psi^{(1)})(\mathbf{r}_1, \dots, \mathbf{r}_N). \qquad (86) $$

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 $\uparrow\downarrow$ are nearby (at distance $\lesssim b_2$) while all other interparticle distances are $\gg b$

• there is a triplet of particles of spins $\uparrow\downarrow\downarrow$ which are nearby (their hyperradius is $\lesssim b$) while all interparticle distances other than the ones within that triplet are $\gg b$

<footnote n="32">³²A relation which is more directly analogous to (85) can be formulated in the mass-imbalanced case, see (129).</footnote>

<footnote n="33">³³In the mass-imbalanced case discussed in App. D, the situation is reversed beyond a critical mass ratio, as was already noted in [175].</footnote>

<footnote n="34">³⁴Apart from the leading-order term (82) which is of order $b$, there are also higher-order contributions to $\delta E_2$, which we expect to contain no term of order $b^{2s}$ (since the three-body physics does not enter in $\delta E_2$) but only integer powers of $b$. Accordingly, the $\propto b^{2s}$ contribution to $\delta E$ should be entirely given by (83) or (84).</footnote>

<footnote n="35">³⁵Specifically, by using the divergence theorem as in Sec. D.2.1, we find the following sufficient conditions: $R^{2+s}\psi^{(1)}$ and $R^{3+s}\partial\psi^{(1)} / \partial R$ tend to zero when $R \to 0$, and similarly, $\bar{R}^{2+s}\psi^{(1)}$ and $\bar{R}^{3+s}\partial\psi^{(1)} / \partial \bar{R}$ tend to zero when $\bar{R} \to 0$.</footnote>
• The relations (83), (84) are reminiscent of the relation [36]

$$ \frac{dE}{d(-1/a_2)} = \frac{\hbar^2}{4\pi m} C_2 \qquad (85) $$

which holds within the zero-range model.<fnref n="32" />

• Typically, $a_3$ is of order $b^{2s}$, which gives $\delta E_3 \sim 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].

• If $b_3 \lesssim b_2$, then $b_2 \sim b$, and since $\delta E_2$ is $\propto r_e$, which is typically of order $b_2$, we have $|\delta E_3 / \delta E_2| \propto b^{2s-1}$. Since $s > 1/2$, we get $|\delta E_3 / \delta E_2| \ll 1$ in the zero-range regime, i.e. the three-body correction to the energy is of higher order than the two-body correction.<fnref n="33" />,<fnref n="34" />

• Let us assume that (84) can be analytically continued to complex $a_3$, with $C_3$ still evaluated within the zero-range model. We then recover the expression (43) of the three-body loss rate, simply by substituting $ImE = Im \delta E_3$ into (64), and using the fact that $\Gamma = \Gamma_3$ and $C_3 \in \mathbb{R}$. Similarly, applying this procedure to (83) yields $\Gamma_3$ in agreement with the sum of Eqs. ((72), (73)).

To derive (83), we consider a stationary state $\psi^{(0)}$ and the associated eigenenergy $E^{(0)}$ of the zero-range model, and the stationary state $\psi^{(1)}$ of the finite-range model whose energy $E^{(1)}$ is close to $E^{(0)}$ in the zero-range regime. We are interested in $\delta E = E^{(1)} - E^{(0)}$.

Using the shorthand notations

$$ \begin{aligned} V(\mathbf{r}_1, \dots, \mathbf{r}_N) &:= \sum_{i:\uparrow} V_2(\mathbf{r}_{i\uparrow}) + \sum_{i<j:\uparrow, k:\downarrow} V_{2,1}(\mathbf{r}_{ijk}) + \sum_{i<j:\uparrow, k:\uparrow} V_{1,2}(\mathbf{r}_{ijk}) \\ H_0 &:= \sum_{i=1}^N \left[ -\frac{\hbar^2}{2m} \Delta_{\mathbf{r}_i} + U(\mathbf{r}_i) \right] \end{aligned} $$

we have $(H_0 + V)\psi^{(1)} = E^{(1)}\psi^{(1)}$, while $\psi^{(0)}$ satisfies $H_0\psi^{(0)} = E^{(0)}\psi^{(0)}$ together with the two-body contact condition (14). Let us consider

$$ \Delta := \langle \psi^{(0)}, (H_0 + V) \psi^{(1)} \rangle - \langle H_0 \psi^{(0)}, \psi^{(1)} \rangle. $$

We have $\Delta = \delta E \langle \psi^{(0)} | \psi^{(1)} \rangle$, and hence $\delta E / \Delta \to 1$ in the zero-range limit.

To evaluate $\Delta$, we write it as $\Delta = \Delta_0 + \langle \psi^{(0)} | V | \psi^{(1)} \rangle$ where $\Delta_0 := \langle \psi^{(0)}, H_0 \psi^{(1)} \rangle - \langle H_0 \psi^{(0)}, \psi^{(1)} \rangle$. We have $\Delta_0 = 0$, assuming for simplicity that the interaction potentials ($V_{2,1}$ and $V_{1,2}$) are finite everywhere (i.e. excluding hard-wall potentials) and do not diverge too quickly at short distances.<fnref n="35" /> Hence

$$ \Delta \approx \int d^3 r_1 \dots d^3 r_N (\psi^{(0)} * V \psi^{(1)})(\mathbf{r}_1, \dots, \mathbf{r}_N). \qquad (86) $$

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 $\uparrow\downarrow$ are nearby (at distance $\lesssim b_2$) while all other interparticle distances are $\gg b$

• there is a triplet of particles of spins $\uparrow\downarrow\downarrow$ which are nearby (their hyperradius is $\lesssim b$) while all interparticle distances other than the ones within that triplet are $\gg b$

<footnote n="32">³²A relation which is more directly analogous to (85) can be formulated in the mass-imbalanced case, see (129).</footnote>

<footnote n="33">³³In the mass-imbalanced case discussed in App. D, the situation is reversed beyond a critical mass ratio, as was already noted in [175].</footnote>

<footnote n="34">³⁴Apart from the leading-order term (82) which is of order $b$, there are also higher-order contributions to $\delta E_2$, which we expect to contain no term of order $b^{2s}$ (since the three-body physics does not enter in $\delta E_2$) but only integer powers of $b$. Accordingly, the $\propto b^{2s}$ contribution to $\delta E$ should be entirely given by (83) or (84).</footnote>

<footnote n="35">³⁵Specifically, by using the divergence theorem as in Sec. D.2.1, we find the following sufficient conditions: $R^{2+s}\psi^{(1)}$ and $R^{3+s}\partial\psi^{(1)} / \partial R$ tend to zero when $R \to 0$, and similarly, $\bar{R}^{2+s}\psi^{(1)}$ and $\bar{R}^{3+s}\partial\psi^{(1)} / \partial \bar{R}$ tend to zero when $\bar{R} \to 0$.</footnote>
• there is a triplet of particles of spins $\uparrow\downarrow\downarrow$ which are nearby (their hyperradius is $\lesssim b$) while all interparticle distances other than the ones within that triplet are $\gg b$.

Denoting the contributions of these regions by $\Delta_2$, $\Delta_{2,1}$ and $\Delta_{1,2}$ respectively, we thus have $\Delta \approx \Delta_2 + \Delta_{2,1} + \Delta_{1,2}$. The two-nearby-particle contribution $\Delta_2$ is given by the r.h.s. of (82), as shown in App. C, in agreement with [56]. It remains to evaluate the contribution $\Delta_{2,1}$ coming from three nearby particles of spins $\uparrow\downarrow\downarrow$. We introduce a length $d_3$ satisfying (53), (54), (70)). $\Delta_{2,1}$ is then given by the contribution to the integral (86) coming from the region of $(\mathbf{r}_1, \dots, \mathbf{r}_N)$ such that there is a triplet of particles $(i, j, k)$ of spins $(\uparrow, \uparrow, \downarrow)$ of hyperradius $R_{ijk} < d_3$ while all interparticle distances other than the ones within the triplet $(i, j, k)$ are $\gg d_3$ (the result does not depend on the value of $d_3$ within the range (53), (54), as we will see). In the region where these conditions are met for the triplet (1,2,3), i.e. when (50), (52) hold, we expect a factorization

$$ \psi^{(1)}(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \sum_{m=-1}^{+1} \Psi_m(\mathbf{R}) B_m^{(1)}(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N), \quad (87) $$

as well as [cf. (15)]

$$ \psi^{(0)}(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \sum_{m=-1}^{+1} \Psi_m^{(0)}(\mathbf{R}) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) $$

with $\Psi_m^{(0)}(\mathbf{R}) := R^{s-2} \phi_m(\Omega)$. Furthermore, we can approximate $B_m^{(1)}$ by $B_m$ in (87), given that we are in the zero-range regime and all distances between the points $\mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$ are $\gg b$. Also using the fact that each triplet gives the same contribution by fermionic antisymmetry, we get

$$ \Delta_{2,1} \approx \frac{N_\uparrow (N_\downarrow - 1) N_\downarrow}{2} \frac{3\sqrt{3}}{8} \sum_{m,m'=-1}^{1} W_{m,m'} \int d^3\mathbf{C} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N (B_m^*(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) B_{m'}^{(1)}(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) \quad (88) $$

where

$$ W_{m,m'} := \int_{R < d_3} d^6\mathbf{R} (\Psi_m^{(0)*} \mathcal{V} \Psi_{m'})(\mathbf{R}) $$

and $\mathcal{V}(\mathbf{R})$ was defined in (76). To evaluate $W_{m,m'}$, we first use the Schrödinger equation (77) to replace $\mathcal{V}$ by $\frac{\hbar^2}{m} \Delta_R$. Since we also have $\Delta_R \Psi^{(0)} = 0$, we can write

$$ W_{m,m'} = \frac{\hbar^2}{m} \int_{R < d_3} d^5\mathbf{R} [\Psi_m^{(0)*} \Delta_R \Psi_{m'} - \Psi_{m'} \Delta_R \Psi_m^{(0)*}](\mathbf{R}). $$

Rewriting the integrand as $\nabla_R \cdot [\Psi_m^{(0)*} \nabla_R \Psi_{m'} - \Psi_{m'} \nabla_R \Psi_m^{(0)*}]$ and applying the divergence theorem yields

$$ W_{m,m'} = \frac{\hbar^2}{m} d^5 \int d^5 \Omega [\Psi_m^{(0)*} \partial_R \Psi_{m'} - \Psi_{m'} \partial_R \Psi_m^{(0)*}]_{R=d_3}. $$

We can then use the asymptotic behavior (78) of $\Psi_m$ (indeed, for $R = d_3 \gg b$, all three interparticle distance are $\gg b$ except in a small region of hyperangles). This gives $W_{m,m'} = \frac{\hbar^2}{m} 2s \delta_{m,m'} a_{2,1}^{(m)}$. Substituting this into (88) and using (17) yields $\Delta_{2,1} \approx \frac{\hbar^2}{m} 2s(s+1) \sum_{m=-1}^{1} C_{2,1}^{(m)} a_{2,1}^{(m)}$. The same reasoning gives an analogous expression for $\Delta_{1,2}$, which yields the final expression (83) for $\delta E_3 \approx \Delta_{2,1} + \Delta_{1,2}$.

## 4. 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 center-of-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
* there is a triplet of particles of spins $\uparrow\downarrow\downarrow$ which are nearby (their hyperradius is $\lesssim b$) while all interparticle distances other than the ones within that triplet are $\gg b$.

Denoting the contributions of these regions by $\Delta_2$, $\Delta_{2,1}$ and $\Delta_{1,2}$ respectively, we thus have $\Delta \approx \Delta_2 + \Delta_{2,1} + \Delta_{1,2}$. The two-nearby-particle contribution $\Delta_2$ is given by the r.h.s. of (82), as shown in App. C, in agreement with [56]. It remains to evaluate the contribution $\Delta_{2,1}$ coming from three nearby particles of spins $\uparrow\downarrow\downarrow$. We introduce a length $d_3$ satisfying (53), (54), (70)). $\Delta_{2,1}$ is then given by the contribution to the integral (86) coming from the region of $(\mathbf{r}_1, \dots, \mathbf{r}_N)$ such that there is a triplet of particles $(i, j, k)$ of spins $(\uparrow, \uparrow, \downarrow)$ of hyperradius $R_{ijk} < d_3$ while all interparticle distances other than the ones within the triplet $(i, j, k)$ are $\gg d_3$ (the result does not depend on the value of $d_3$ within the range (53), (54), as we will see). In the region where these conditions are met for the triplet (1,2,3), i.e. when (50), (52) hold, we expect a factorization

$$ \psi^{(1)}(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \sum_{m=-1}^{+1} \Psi_m(\mathbf{R}) B_m^{(1)}(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N), \quad (87) $$

as well as [cf. (15)]

$$ \psi^{(0)}(\mathbf{r}_1, \dots, \mathbf{r}_N) \approx \sum_{m=-1}^{+1} \Psi_m^{(0)}(\mathbf{R}) B_m(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) $$

with $\Psi_m^{(0)}(\mathbf{R}) := R^{s-2} \phi_m(\Omega)$. Furthermore, we can approximate $B_m^{(1)}$ by $B_m$ in (87), given that we are in the zero-range regime and all distances between the points $\mathbf{C}, \mathbf{r}_4, \dots, \mathbf{r}_N$ are $\gg b$. Also using the fact that each triplet gives the same contribution by fermionic antisymmetry, we get

$$ \Delta_{2,1} \approx \frac{N_\perp (N_\perp - 1) N_\parallel}{2} \frac{3\sqrt{3}}{8} \sum_{m,m'=-1}^{1} W_{m,m'} \int d^3\mathbf{C} d^3\mathbf{r}_4 \dots d^3\mathbf{r}_N (B_m^*(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) B_{m'}^{(1)}(\mathbf{C}; \mathbf{r}_4, \dots, \mathbf{r}_N) \quad (88) $$

where

$$ W_{m,m'} := \int_{R < d_3} d^6 \mathbf{R} (\Psi_m^{(0)*} \mathcal{V} \Psi_{m'})(\mathbf{R}) $$

and $\mathcal{V}(\mathbf{R})$ was defined in (76). To evaluate $W_{m,m'}$, we first use the Schrödinger equation (77) to replace $\mathcal{V}$ by $\frac{\hbar^2}{m} \Delta_R$. Since we also have $\Delta_R \Psi^{(0)} = 0$, we can write

$$ W_{m,m'} = \frac{\hbar^2}{m} \int_{R < d_3} d^5 \Omega [\Psi_m^{(0)*} \Delta_R \Psi_{m'} - \Psi_{m'} \Delta_R \Psi_m^{(0)*}](\mathbf{R}). $$

Rewriting the integrand as $\nabla_R \cdot [\Psi_m^{(0)*} \nabla_R \Psi_{m'} - \Psi_{m'} \nabla_R \Psi_m^{(0)*}]$ and applying the divergence theorem yields

$$ W_{m,m'} = \frac{\hbar^2}{m} d^5 \int d^5 \Omega [\Psi_m^{(0)*} \partial_R \Psi_{m'} - \Psi_{m'} \partial_R \Psi_m^{(0)*}]_{R=d_3}. $$

We can then use the asymptotic behavior (78) of $\Psi_m$ (indeed, for $R = d_3 \gg b$, all three interparticle distance are $\gg b$ except in a small region of hyperangles). This gives $W_{m,m'} = \frac{\hbar^2}{m} 2s \delta_{m,m'} a_{2,1}^{(m)}$. Substituting this into (88) and using (17) yields $\Delta_{2,1} \approx \frac{\hbar^2}{m} 2s(s+1) \sum_{m=-1}^{1} C_{2,1}^{(m)} a_{2,1}^{(m)}$. The same reasoning gives an analogous expression for $\Delta_{1,2}$, which yields the final expression (83) for $\delta E_3 \approx \Delta_{2,1} + \Delta_{1,2}$.

## 4. 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 center-of-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 $\Gamma_3$ of deeply bound dimers by three-body recombination equals $C_3 \text{Im } a_3$ times an explicit prefactor; we also showed that the finite-range correction to the energy has a contribution equal to $-\frac{\hbar}{2} C_3 a_3$ times the same prefactor. With respect to the relation between $\Gamma_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 $\text{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 $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ for the associated three-body problem, which can be further broken up into the contributions $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ for each value $m \in \{-1, 0, 1\}$ of the angular-momentum quantum number of the three-body problem. Most relations only involve $C_3$, $C_{2,1}$ or $C_{1,2}$. The only relations involving $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are the ones for finite-range interactions, in the general case where interactions are not invariant under rotation and under exchange between $\uparrow$ and $\downarrow$. In this case, the three-body parameter also depends on the spin and angular momentum indices.

Nevertheless, for a homogeneous unpolarized gas, $\Gamma_3$ simply equals $C_3 \text{Im } \bar{a}_3$ times an explicit prefactor, with $\bar{a}_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 $\text{Im } \bar{a}_3$ by measuring $\Gamma_3$. Measuring $\Gamma_3$ in the low-temperature regime would then allow one to experimentally test the power-law $\mathcal{C}_3 = \zeta_3 n^{(2s+5)/3}$ and determine the associated many-body parameter $\zeta_3$ whose computation is an open theoretical challenge.

## Acknowledgments

We are grateful to S. Tan, C. De Daniloff and F. Chevy for particularly important input. We also thank Y. Castin, N. Navon, D. Petrov, T. Yefsah and the ultracold Fermi gases team at LKB for stimulating discussions.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Appendix A. The unitary three-body problem

In this Appendix we review, for self-containedness of the article, the known solution of the three-body problem at the unitary limit [96, 176] (see also Refs. [172, 177]). Consider an eigenstate $\psi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ of the three-body problem in free space. With the convention that particles (1,2,3) have spins $(\uparrow, \uparrow, \downarrow)$, the fermionic antisymmetry reads $\psi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = -\psi(\mathbf{r}_2, \mathbf{r}_1, \mathbf{r}_3)$. Restricting to zero center-of-mass momentum, $\psi$ only depends on the relative positions between the three particles, which can be parameterized by the Jacobi coordinates

$$ \begin{gathered} \mathbf{r} = \mathbf{r}_3 - \mathbf{r}_1 \\ \frac{\sqrt{3}}{2} \boldsymbol{\rho} = \mathbf{r}_2 - \frac{\mathbf{r}_1 + \mathbf{r}_3}{2}. \end{gathered} \qquad (89) $$

It is convenient to introduce the six-dimensional vector

$$ \mathbf{R} = (\mathbf{r}, \boldsymbol{\rho}) \qquad (90) $$
range, we introduced a small three-body parameter $a_3$, and we showed that the formation rate $\Gamma_3$ of deeply bound dimers by three-body recombination equals $C_3 \text{Im } a_3$ times an explicit prefactor; we also showed that the finite-range correction to the energy has a contribution equal to $-\frac{\hbar}{2} C_3 a_3$ times the same prefactor. With respect to the relation between $\Gamma_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 $\text{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 $\uparrow\uparrow$ and $\uparrow\downarrow\downarrow$ for the associated three-body problem, which can be further broken up into the contributions $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ for each value $m \in \{-1, 0, 1\}$ of the angular-momentum quantum number of the three-body problem. Most relations only involve $C_3$, $C_{2,1}$ or $C_{1,2}$. The only relations involving $C_{2,1}^{(m)}$ and $C_{1,2}^{(m)}$ are the ones for finite-range interactions, in the general case where interactions are not invariant under rotation and under exchange between $\uparrow$ and $\downarrow$. In this case, the three-body parameter also depends on the spin and angular momentum indices.

Nevertheless, for a homogeneous unpolarized gas, $\Gamma_3$ simply equals $C_3 \text{Im } \bar{a}_3$ times an explicit prefactor, with $\bar{a}_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 $\text{Im } \bar{a}_3$ by measuring $\Gamma_3$. Measuring $\Gamma_3$ in the low-temperature regime would then allow one to experimentally test the power-law $\mathcal{C}_3 = \zeta_3 n^{(2s+5)/3}$ and determine the associated many-body parameter $\zeta_3$ whose computation is an open theoretical challenge.

## Acknowledgments

We are grateful to S. Tan, C. De Daniloff and F. Chevy for particularly important input. We also thank Y. Castin, N. Navon, D. Petrov, T. Yefsah and the *ultracold Fermi gases* team at LKB for stimulating discussions.

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## Appendix A. The unitary three-body problem

In this Appendix we review, for self-containedness of the article, the known solution of the three-body problem at the unitary limit [96, 176] (see also Refs. [172, 177]). Consider an eigenstate $\psi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3)$ of the three-body problem in free space. With the convention that particles (1,2,3) have spins $(\uparrow, \uparrow, \downarrow)$, the fermionic antisymmetry reads $\psi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = -\psi(\mathbf{r}_2, \mathbf{r}_1, \mathbf{r}_3)$. Restricting to zero center-of-mass momentum, $\psi$ only depends on the relative positions between the three particles, which can be parameterized by the Jacobi coordinates

$$ \begin{gathered} \mathbf{r} = \mathbf{r}_3 - \mathbf{r}_1 \\ \frac{\sqrt{3}}{2} \rho = \mathbf{r}_2 - \frac{\mathbf{r}_1 + \mathbf{r}_3}{2}. \end{gathered} \qquad (89) $$

It is convenient to introduce the six-dimensional vector

$$ \mathbf{R} = (\mathbf{r}, \rho) \qquad (90) $$
whose norm $R = \|\mathbf{R}\| = \sqrt{r^2 + \rho^2}$ is the hyperradius, while its direction $\mathbf{R}/R$ can be parameterized by five hyperangles denoted collectively by $\mathbf{\Omega}$,

$$ \mathbf{\Omega} \iff \mathbf{R}/R. \qquad (91) $$

The three-body Schrödinger equation then writes

$$ -\frac{\hbar^2}{m} \Delta_{\mathbf{R}} \psi(\mathbf{R}) = E \psi(\mathbf{R}) \qquad (92) $$

with the contact condition: there exists $\mathcal{A}$ such that

$$ \psi(\mathbf{R}) \underset{r \to 0}{=} \left( \frac{1}{r} - \frac{1}{a_2} \right) \mathcal{A}(\boldsymbol{\rho}) + O(r) \qquad (93) $$

(the contact condition between particles 2 and 3 then automatically holds by antisymmetry).

At the unitary limit $a_2 = \infty$, this contact condition becomes scale invariant. Accordingly, it only acts on the hyperangles; explicitly, it can be expressed as $\frac{\partial}{\partial \alpha} (\sin \alpha \psi)_{\alpha=0} = 0$ where $\alpha = \arctan(r/\rho)$. As a result, the unitary three-body problem is separable in hyperspherical coordinates: One can look for eigenstates of the factorized form

$$ \psi(\mathbf{R}) = \frac{F(R)}{R^2} \phi(\mathbf{\Omega}) $$

(where the factor $1/R^2$ is introduced for later convenience), and the three-body problem (92), (93)) separates into

* a hyperradial problem

$$ \left( -\frac{d^2}{dR^2} - \frac{1}{R} \frac{d}{dR} + \frac{s^2}{R^2} \right) F(R) = \frac{m}{\hbar^2} E F(R) \qquad (94) $$

* a hyperangular problem, defined by

$$ T_{\Omega} \phi(\Omega) = -s^2 \phi(\Omega) \qquad (95) $$

together with the contact condition

$$ \frac{\partial}{\partial \alpha} \left[ \sin \alpha \phi(\Omega) \right]_{\alpha=0} = 0 \qquad (96) $$

and the antisymmetry constraint.

Here $T_\Omega$ is a differential operator acting on the hyperangles, defined by

$$ \Delta_{\mathbf{R}} = \frac{1}{R^2} \left( \frac{\partial^2}{\partial R^2} + \frac{1}{R} \frac{\partial}{\partial R} + \frac{1}{R^2} T_{\Omega} \right) R^2. \qquad (97) $$

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_v\}$ where $v$ is a discrete index. The corresponding hyperangular eigenfunctions $\{\phi_v(\Omega)\}$ form an orthonormal basis for the hyperangular scalar product

$$ (f|g) \equiv \int f(\Omega)^* g(\Omega) d^5\Omega \qquad (98) $$

where $d^5\Omega$ denotes the differential solid angle in six-dimensional space, $d^6R = d^5\Omega R^5 dR$; 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$^{36}$ (for a mathematical proof of self-adjointness in free space, see Refs. [178–180]).

The hyperangular problem (95), (96) is analytically solvable, with two types of solutions. The first type are common solutions of the unitary and non-interacting hyperangular problems, whose wavefunction $\phi(\Omega)$ vanishes when two particles approach each other (i.e. for $\alpha \to 0$); the corresponding $s$ take integer values. The second type are the following truly interacting solutions:

<footnote>$^{36}$See [177, Sections 3.2.b and 3.2.c].</footnote>
whose norm $R = \|\mathbf{R}\| = \sqrt{r^2 + \rho^2}$ is the hyperradius, while its direction $\mathbf{R}/R$ can be parameterized by five hyperangles denoted collectively by $\mathbf{\Omega}$,

$$ \mathbf{\Omega} \iff \mathbf{R}/R. \qquad (91) $$

The three-body Schrödinger equation then writes

$$ -\frac{\hbar^2}{m} \Delta_{\mathbf{R}} \psi(\mathbf{R}) = E \psi(\mathbf{R}) \qquad (92) $$

with the contact condition: there exists $\mathcal{A}$ such that

$$ \psi(\mathbf{R}) \underset{r \to 0}{=} \left(\frac{1}{r} - \frac{1}{a_2}\right) \mathcal{A}(\boldsymbol{\rho}) + O(r) \qquad (93) $$

(the contact condition between particles 2 and 3 then automatically holds by antisymmetry).

At the unitary limit $a_2 = \infty$, this contact condition becomes scale invariant. Accordingly, it only acts on the hyperangles; explicitly, it can be expressed as $\frac{\partial}{\partial \alpha} (\sin \alpha \psi)_{\alpha=0} = 0$ where $\alpha = \arctan(r/\rho)$. As a result, the unitary three-body problem is separable in hyperspherical coordinates: One can look for eigenstates of the factorized form

$$ \psi(\mathbf{R}) = \frac{F(R)}{R^2} \phi(\mathbf{\Omega}) $$

(where the factor $1/R^2$ is introduced for later convenience), and the three-body problem (92), (93)) separates into

* a hyperradial problem

$$ \left(-\frac{d^2}{dR^2} - \frac{1}{R}\frac{d}{dR} + \frac{s^2}{R^2}\right) F(R) = \frac{m}{\hbar^2} EF(R) \qquad (94) $$

* a hyperangular problem, defined by

$$ T_{\Omega} \phi(\Omega) = -s^2 \phi(\Omega) \qquad (95) $$

together with the contact condition

$$ \frac{\partial}{\partial \alpha} \left[ \sin \alpha \phi(\Omega) \right]_{\alpha=0} = 0 \qquad (96) $$

and the antisymmetry constraint.

Here $T_\Omega$ is a differential operator acting on the hyperangles, defined by

$$ \Delta_{\mathbf{R}} = \frac{1}{R^2} \left( \frac{\partial^2}{\partial R^2} + \frac{1}{R} \frac{\partial}{\partial R} + \frac{1}{R^2} T_{\Omega} \right) R^2. \qquad (97) $$

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_v\}$ where $v$ is a discrete index. The corresponding hyperangular eigenfunctions $\{\phi_v(\Omega)\}$ form an orthonormal basis for the hyperangular scalar product

$$ (f|g) \equiv \int f(\Omega)^* g(\Omega) d^5\Omega \qquad (98) $$

where $d^5\Omega$ denotes the differential solid angle in six-dimensional space, $d^6R = d^5\Omega R^5 dR$; 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$^{36}$ (for a mathematical proof of self-adjointness in free space, see Refs. [178–180]).

The hyperangular problem (95), (96) is analytically solvable, with two types of solutions. The first type are common solutions of the unitary and non-interacting hyperangular problems, whose wavefunction $\phi(\Omega)$ vanishes when two particles approach each other (i.e. for $\alpha \to 0$); the corresponding $s$ take integer values. The second type are the following truly interacting solutions:

<footnote>$^{36}$See [177, Sections 3.2.b and 3.2.c].</footnote>