# Comptes Rendus
Physique

Shimpei Endo et Yvan Castin

## Quatrièmes coefficients d'amas et du viriel d'un gaz unitaire de fermions pour un rapport de masse quelconque

Volume 23 (2022), p. 41-110

Published online: 12 August 2022

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

<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/

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

Les Comptes Rendus. Physique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1878-1535
coefficients of order $n + 1$ for an infinity-mass impurity fermion, matching the conjecture for $n = 3$. Finally, in a harmonic potential, we predict a non-monotonic behavior of the $3 + 1$ cluster coefficient with trapping frequency, near mass ratios where this coefficient vanishes in the homogeneous case. A multilingual version is available in separate files on the open archive HAL at https://hal.archives-ouvertes.fr/hal-03592961.

**Mots-clés.** gaz de fermions, limite unitaire, invariance d'échelle, développement du viriel, développement en amas.

**Keywords.** Fermi gases, unitary limit, scale invariance, virial expansion, cluster expansion.

Manuscrit reçu le 28 février 2022, révisé le 11 avril 2022, accepté le 17 mai 2022.

**The English version is available after the French version**

# 1. Présentation du problème et principaux résultats

## 1.1. Le système

Notre objet d'étude est un gaz tridimensionnel de fermions neutres non relativistes à deux états de spin $\uparrow$ et $\downarrow$ dans le régime d'interaction maximale permise en phase gazeuse c'est-à-dire dans la limite dite unitaire [1, 2] : il n'y a pas d'interaction entre fermions dans le même état de spin mais une interaction binaire de portée nulle, exclusivement dans l'onde $s$ et de longueur de diffusion infinie entre fermions d'états de spin différents $\uparrow$ et $\downarrow$.<fnref n="1" /> Nous considérons le cas général où les masses $m_\uparrow$ et $m_\downarrow$ d'un fermion peuvent être différentes dans les deux états de spin, ce qui fait l'originalité de notre travail. Ce système, loin d'être une vue de l'esprit, est réalisable en laboratoire avec des gaz d'atomes froids piégés. Ses principales propriétés à l'équilibre ont été déterminées expérimentalement pour des masses égales $m_\uparrow = m_\downarrow$, en particulier sa superfluidité à basse température dans le cas non polarisé (à nombres de fermions égaux dans $\uparrow$ et $\downarrow$) [3, 4] et son équation d'état à température et polarisation quelconques [5-7]. Le cas $m_\uparrow \neq m_\downarrow$ reste à explorer; une façon naturelle de l'obtenir est d'utiliser un mélange de deux espèces atomiques fermioniques complètement polarisées, comme il a été fait dans les références [8-10].

## 1.2. Développements en amas et du viriel

Compte tenu de la force des interactions, qui ne fournit aucun petit paramètre évident, il y a peu d'outils théoriques permettant des prédictions quantitativement fiables sur le gaz unitaire qu'on puisse comparer aux mesures. L'un d'entre eux est la simulation de Monte-Carlo diagrammatique du problème à N corps sur ordinateur [11, 12]. Un autre est le développement en amas ou du viriel [13], c'est-à-dire le développement de la pression P du gaz infini spatialement homogène en puissances des fugacités $z_σ = \exp(\beta\mu_σ)$ ou des densités dans l'espace des phases $ρ_σ\lambda_σ^3$ dans le régime fortement non dégénéré où elles tendent vers zéro, avec $ρ_σ$ la densité de la composante

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_{n_\uparrow+n_\downarrow}) = A_{ij}(a^{-1} - r_{ij}^{-1}) + O(r_{ij}) $$

où $a$ est la longueur de diffusion, et ceci pour tout $i$ et pour tout $j$. Ici, nous sommes dans la limite unitaire $1/a = 0$.

<footnote n="1">¹ Cette situation correspond au remplacement du potentiel d'interaction par les conditions de contact de Wigner-Bethe-Peierls sur la fonction d'onde $\psi$ du système dans chaque configuration de spin $|\uparrow\dots\uparrow\downarrow\dots\downarrow\rangle = |\uparrow\rangle^{\otimes n_\uparrow} |\downarrow\rangle^{\otimes n_\downarrow}$ : lorsque la distance $r_{ij}$ entre la particule $i$ dans l'état de spin $\uparrow$ et la particule $j$ dans l'état de spin $\downarrow$ tend vers zéro à position fixée $\mathbf{R}_{ij}$ de leur centre de masse, les positions $\mathbf{r}_k$ des autres particules étant fixées à des valeurs différentes de $\mathbf{R}_{ij}$, il existe une constante $A_{ij}$ (fonction de $\mathbf{R}_{ij}$ et des $\mathbf{r}_k$) telle que</footnote>
$\sigma = \uparrow, \downarrow$ dans l'espace réel, $\mu_\sigma$ son potentiel chimique, $\lambda_\sigma = (2\pi\hbar^2/m_\sigma k_B T)^{1/2}$ sa longueur d'onde thermique de de Broglie à la température $T$ et $\beta = 1/k_B T$ [14] :

$$ \frac{P\lambda^3}{k_B T} = \sum_{(n_\uparrow, n_\downarrow) \in \mathbb{N}^{2*}} b_{n_\uparrow, n_\downarrow} z_\uparrow^{n_\uparrow} z_\downarrow^{n_\downarrow} = \sum_{(n_\uparrow, n_\downarrow) \in \mathbb{N}^{2*}} c_{n_\uparrow, n_\downarrow} (\rho_\uparrow \lambda_\uparrow^3)^{n_\uparrow} (\rho_\downarrow \lambda_\downarrow^3)^{n_\downarrow} \quad (1) $$

Pour adimensionner la pression, il a fallu introduire une longueur d'onde thermique de de Broglie de référence $\lambda = (2\pi\hbar^2/\tilde{m}k_B T)^{1/2}$ dépendant d'une masse moyenne $\tilde{m}$ à préciser (voir l'équation (12)). On peut alors tenter d'extrapoler au régime non trivial $z_\sigma \approx 1$ avec des recettes heuristiques comme l'approximant de Padé [15] ou des méthodes de resommation optimisées [16] qui tiennent compte du comportement des coefficients aux grands ordres. Le développement en amas ou du viriel a notre préférence ici sur la simulation Monte-Carlo car il est plus proche du calcul analytique, les coefficients d'ordre $n$ se déduisant de la solution d'un problème à au plus $n$ fermions en interaction, à petit nombre de corps donc : pour obtenir $b_{n_\uparrow, n_\downarrow}$, il suffit de déterminer les fonctions de partition canoniques de tous les systèmes ayant un nombre de fermions inférieur ou égal à $n_\sigma$ dans chaque état de spin $\sigma = \uparrow, \downarrow$.

### 1.3. Simplicité de la limite unitaire

L'invariance d'échelle du gaz unitaire simplifie considérablement le calcul des coefficients d'amas pour $n > 2$ (l'ordre deux, donné en toute généralité par la formule de Beth-Uhlenbeck [17–19], ne fait pas débat). Ainsi, les coefficients d'ordre trois sont connus analytiquement, même s'il y a brisure de l'invariance d'échelle au niveau du trois corps par l'effet Efimov [20, 21], en net contraste avec le modèle d'une interaction de sphère dure de rayon $a$ où les coefficients ne sont connus analytiquement (pour des bosons) que dans les cas limites $\lambda/a \ll 1$ [22–25] ou $\gg 1$ [26–29]. La méthode du régulateur harmonique [30–33], consistant à piéger chaque composante du gaz dans un potentiel harmonique isotrope fictif, $U_\sigma(\mathbf{r}) = m_\sigma\omega^2 r^2/2$, dont on fait tendre vers zéro la pulsation de piégeage $\omega$ (commune aux deux états de spin) à la fin des calculs, permet de tirer parti au mieux de l'invariance d'échelle puisque le spectre à $n$ corps dans le piège se déduit de l'ensemble discret des exposants d'échelle $s_i$ des états propres d'énergie nulle dans l'espace libre [1, 34]. Plus précisément, on généralise le développement en amas au cas piégé, en remplaçant la pression par le grand potentiel $\Omega$ et en la ramenant à la fonction de partition $Z_1 = 1/[2\text{sh}(\tilde{\omega}/2)]^3$ d'un fermion seul,<fnref n="2" /> si bien que

$$ \frac{-\Omega}{k_B T Z_1} = \sum_{(n_\uparrow, n_\downarrow) \in \mathbb{N}^{2*}} B_{n_\uparrow, n_\downarrow}(\tilde{\omega}) z_\uparrow^{n_\uparrow} z_\downarrow^{n_\downarrow} \quad (2) $$

Les coefficients du cas piégé, repérés par une lettre majuscule pour éviter toute confusion, dépendent seulement du rapport sans dimension $\tilde{\omega} = \hbar\omega/k_B T$ par invariance d'échelle du gaz unitaire. On les relie à ceux du cas homogène au moyen de l'approximation d'homogénéité locale, exacte dans la limite $\omega \to 0$ [32, 33], pour obtenir [21]

$$ b_{n_\uparrow, n_\downarrow} = \frac{(n_\uparrow m_\uparrow + n_\downarrow m_\downarrow)^{3/2}}{\tilde{m}^{3/2}} B_{n_\uparrow, n_\downarrow} (0^+) \quad (3) $$

Cette méthode conduit à une expression analytique intégrale des coefficients d'ordre trois car la fonction transcendante d'Efimov $\Lambda(s)$, dont les exposants d'échelle $s_i$ sont les racines, est connue explicitement [35–38] ce qui permet d'exprimer les coefficients sous forme d'une intégrale de contour autour de $\mathbb{R}^+$ au moyen du théorème des résidus puis de rabattre le contour sur l'axe

<footnote n="2">² Le numérateur au premier membre de l'équation (1) peut être vu, dans une boîte de quantification de volume $V$ arbitrairement grand, comme le rapport entre l'opposé $PV$ du grand potentiel du gaz et la fonction de partition $V/\lambda^3$ d'une particule seule fictive de masse $\tilde{m}$.</footnote>
imaginaire pur par analyticité sur $C \backslash \mathbb{R}$ [20]. En revanche, pour les coefficients d'ordre quatre non triviaux,<fnref n="3" /> la fonction d'Efimov est le déterminant d'opérateurs $M_{3,1}(s)$ [39] ou $M_{2,2}(s)$ [40] dépendant paramétriquement de s, que l'on sait calculer numériquement seulement sur l'axe imaginaire pur $s = iS$ et dont les propriétés analytiques dans le plan complexe, imparfaitement connues, ne permettent pas une utilisation sûre du théorème des résidus. Aussi la référence [41] n'a-t-elle pu produire qu'une conjecture, que nous rappelons brièvement.

## 1.4. Notre conjecture de 2016 sur les quatrièmes coefficients d'amas

Comme dans la référence [40], on introduit l'expression intégrale calquée sur celle des coefficients d'ordre trois (on y remplace $\Lambda(s)$ par un déterminant d'opérateur) :

$$I_{n_1, n_1}(\bar{\omega}) = \sum_{\ell=0}^{+\infty} \sum_{\epsilon} (\ell + 1/2) \int_{-\infty}^{+\infty} \frac{dS}{2\pi} \frac{\sin(\bar{\omega}S)}{\sin\bar{\omega}} \frac{d}{dS} \ln \det \left[ M_{n_1, n_1}^{(\ell, \epsilon)}(iS) \right] \quad (4)$$

où la somme porte sur le moment cinétique $\ell$ interne (c'est-à-dire après séparation du centre de masse) des états propres à quatre corps et sur leur parité interne $\epsilon = \pm 1$ (limitée à $\epsilon = 1$ pour $\ell = 0$), et où l'opérateur $M_{n_1, n_1}(iS)$ est restreint au sous-espace $(\ell, \epsilon)$ correspondant. Alors $I_{n_1, n_1}(\bar{\omega})$ donne exactement $B_{n_1, n_1}(\bar{\omega})$ (comme c'est le cas à l'ordre trois) lorsque les objets asymptotiques découplés (des paquets indépendants non monoatomiques de fermions fortement corrélés par les interactions, en termes desquels s'expriment les états propres très excités) sont discernables; sinon, il manque l'effet de statistique quantique lié à l'indiscernabilité de ces objets, qui contribue au coefficient d'amas même s'ils n'interagissent pas entre eux, pour la même raison que les coefficients d'amas du gaz parfait quantique diffèrent de ceux du gaz parfait classique. Dans la configuration $(n_1, n_1) = (3, 1)$, dite désormais 3 + 1 pour abréger, les objets asymptotiques possibles sont un triplon $\uparrow\uparrow\downarrow$ ou un pairon $\uparrow\downarrow$ de fermions fortement corrélés; il ne peut y en avoir qu'un à la fois, ce qui exclut tout effet de statistique quantique. En revanche, dans la configuration 2 + 2, les fermions peuvent également se découpler en deux pairons $\uparrow\downarrow$ de fermions corrélés, qui sont des bosons indiscernables; comme ces bosons n'interagissent pas, on calcule facilement leur fonction de partition [41]. La conjecture de la référence [41] s'écrit en définitive :

$$B_{3,1}(\bar{\omega}) = I_{3,1}(\bar{\omega}) \quad ; \quad B_{2,2}(\bar{\omega}) = I_{2,2}(\bar{\omega}) + \frac{1}{32} \frac{1}{\mathrm{ch}\bar{\omega}\mathrm{ch}^3(\bar{\omega}/2)} \quad (5)$$

le cas 1 + 3 se déduisant du cas 3 + 1 par échange des deux états de spin donc changement du rapport de masse $\alpha = m_1/m_1$ en son inverse $1/\alpha$. Pour $\alpha = 1$, la conjecture (5) est en accord avec le calcul de Monte-Carlo quantique de la référence [42] jusqu'à la valeur minimale de $\bar{\omega}$ accessible, $\bar{\omega} \approx 1$; elle est en accord avec un calcul numérique récent plus performant jusqu'à des valeurs de $\bar{\omega} \ll 1$ [43, 44]. Après utilisation de la relation (3), elle est aussi en accord avec le même calcul numérique performant effectué directement dans le cas spatialement homogène d'une boîte de quantification [45]. La conjecture (5) est donc confirmée pour des masses égales.

## 1.5. Contenu de l'étude

L'idée du présent travail est de croire à la validité de la conjecture (5) pour un rapport de masse quelconque $\alpha$ et d'accéder aux quatrièmes coefficients d'amas et du viriel correspondants, en

<footnote n="3">³ Si l'un des $n_\sigma$ est nul, les fermions du problème à $n$ corps sont sans interaction et le coefficient $B$ correspondant se réduit à celui du gaz parfait. On trouve ainsi $B_{n,0}(0^+) = B_{0,n}(0^+) = (-1)^{n+1}/n^4$ et, au moyen de la relation (3),</footnote>

<footnote>$$b_{n,0} = (m_1/\bar{m})^{3/2} (-1)^{n+1} / n^{5/2} \quad \text{et} \quad b_{0,n} = (m_1/\bar{m})^{3/2} (-1)^{n+1} / n^{5/2}.$$</footnote>
calculant numériquement le déterminant des opérateurs $M_{n_1, n_1}^{(\ell, \epsilon)}(iS)$, puis en intégrant sur $S$ et en sommant sur $\ell$ et $\epsilon$ dans l'expression (4). En réalité, il nous faut quand même mettre des bornes sur le rapport de masse, car la méthode de résolution du problème à quatre corps dans un piège harmonique, à la base de l'expression (5), suppose la séparabilité de l'équation de Schrödinger interne en coordonnées hypersphériques; ceci n'est vrai que s'il n'y a pas de brisure d'invariance d'échelle par effet Efimov au niveau du trois corps, ce qui nous contraint aux intervalles

$$ 
\begin{array}{l@{\qquad}l}
\text{pour } B_{3,1}: \alpha < \alpha_c^{3\text{corps}} \approx 13,60697 & \text{pour } B_{1,3}: \alpha > \frac{1}{\alpha_c^{3\text{corps}}} \approx 0,0734917 \\
& \text{pour } B_{2,2}: \frac{1}{\alpha_c^{3\text{corps}}} < \alpha < \alpha_c^{3\text{corps}}
\end{array}
\qquad (6) $$

où $\alpha_c^{3\text{corps}}$, encore noté $\alpha_c^{2,1}$, est le seuil d'apparition de l'effet Efimov à trois corps dans le système $\uparrow\uparrow\downarrow$, qui se produit dans la voie de moment cinétique interne à trois corps $L = 1$ [35, 46]. Sous la forme (5), la conjecture suppose aussi qu'il n'y a pas d'effet Efimov à quatre corps<fnref n="4" />; un tel effet ne se produit que dans la configuration $3+1$ ou $1+3$ [39, 40], au rapport de masse $\alpha_c^{4\text{corps}} = \alpha_c^{3,1}$ ou à son inverse, d'où les conditions supplémentaires

$$ \text{pour } B_{3,1}: \alpha < \alpha_c^{4\text{corps}} \approx 13,3842 \quad ; \quad \text{pour } B_{1,3}: \alpha > \frac{1}{\alpha_c^{4\text{corps}}} \approx 0,074715 \qquad (7) $$

Une fois le développement en amas connu, les coefficients du viriel $c_{n_1, n_1}$ s'en déduisent aisément, sous la forme de fonctions rationnelles des coefficients d'amas $b_{n'_1, n'_1}$ d'ordre total $n'_1 + n'_1$ inférieur ou égal à $n_1 + n'_1$ :<fnref n="5" />

$$ \rho_{\sigma} \lambda^3 = \sum_{(n_1, n_1) \in \mathbb{N}^{2*}} n_{\sigma} b_{n_1, n_1} z_1^{n_1} z_1^{n_1} $$

tiré de la relation thermodynamique $\rho_{\sigma} = \partial P / \partial \mu_{\sigma}$, et l'on ajuste les $c_{n_1, n_1}$ ordre par ordre pour retrouver le développement au second membre de l'équation (1). Il est plus simple, dans un premier temps, d'écrire le développement du viriel sous la forme

$$ P \frac{\lambda^3}{k_B T} = \sum_{(n_1, n_1) \in \mathbb{N}^{2*}} a_{n_1, n_1} (\rho_1 \lambda^3)^{n_1} (\rho_1 \lambda^3)^{n_1} $$

puis, une fois les coefficients $a_{n_1, n_1}$ calculés, d'utiliser la relation

$$ c_{n_1, n_1} = a_{n_1, n_1} / (t_1^{n_1} t_1^{n_1}) $$

où l'on a posé $t_{\sigma} = (\bar{m}/m_{\sigma})^{3/2}$. On notera que $b_{1,0} t_1 = b_{0,1} t_1 = 1$, et que $a_{1,0} = a_{0,1} = 1$ quelle que soit la masse de référence $\bar{m}$ en vertu de la loi du gaz parfait $P = (\rho_1 + \rho_1) k_B T$. En remplaçant les $b_{n,0}$ et $b_{0,n}$ par leurs expressions données dans la note 3, il vient finalement

$$ c_{3,1} = \left(-\frac{3}{4} + \frac{1}{\sqrt{3}}\right) b_{1,1} - \left(\frac{3}{\sqrt{2}}\right) b_{2,1} - 3b_{3,1} + 3b_{1,1} \left(b_{2,1} + b_{1,1}/\sqrt{8}\right) t_1 - b_{1,1}^3 t_1^2 $$

$$ c_{2,2} = -3 \left[ b_{2,2} + \left(b_{1,2} + b_{2,1}\right)/\sqrt{8} + b_{1,1}/8 \right] + \left(9\sqrt{2}/8\right) b_{1,1}^2 (t_1 + t_1) + 6b_{1,1} \left(b_{2,1} t_1 + b_{1,2} t_1\right) - 3b_{1,1}^3 t_1 t_1 $$

(les coefficients du gaz parfait $c_{4,0}$ et $c_{0,4}$ sont donnés sous forme explicite plus loin dans le texte principal). On tire facilement $b_{1,1}$ de l'équation (3) sachant que $B_{1,1}(0^+) = 1/2$. En revanche, $b_{2,1}$ et $b_{1,2}$ sont calculés numériquement à partir des expressions intégrales explicites de la référence [21].

<footnote n="4">⁴ De plus, notre hypothèse d'une interaction strictement invariante d'échelle devient difficile à satisfaire expérimentalement si $\alpha$, bien que du côté non efimovien, est trop proche du seuil d'apparition des tétramères. Pour s'affranchir de ces limitations, on pourrait compléter la condition de contact à 3 + 1 corps de l'habituel modèle d'interaction de portée nulle de la note 1 au moyen d'une longueur appelée « paramètre à 3+1 corps » dans la voie $(\ell, \epsilon) = (1, +1)$ [39] où se produit l'effet Efimov à 3+1 corps. Dans ce cas, $B_{3,1}(0^+)$ et $b_{3,1}$ seraient des fonctions lisses de $\alpha$ en $\alpha = \alpha_c^{4\text{corps}}$. La référence [21] met en œuvre ce programme pour $B_{2,1}(0^+)$ c'est-à-dire dans le problème à 2 + 1 corps.</footnote>

<footnote n="5">⁵ Pour les obtenir, on remplace dans le développement de la pression au troisième membre de l'équation (1) les densités par leur développement du viriel</footnote>
$$c_{1,0} = b_{1,0} \ ; \ c_{2,0} = -b_{2,0} \ ; \ c_{3,0} = \frac{4b_{2,0}^2}{b_{1,0}} - 2b_{3,0} \ ; \qquad (8)$$

$$c_{4,0} = -3b_{4,0} - \frac{20b_{2,0}^3}{b_{1,0}^2} + \frac{18b_{3,0}b_{2,0}}{b_{1,0}} \ ; \ c_{1,1} = -b_{1,1}$$

$$c_{2,1} = -2b_{2,1} + \frac{4b_{1,1}b_{2,0}}{b_{1,0}} + \frac{b_{1,1}^2}{b_{0,1}} \ ; \qquad (9)$$

$$c_{3,1} = -3b_{3,1} - \frac{b_{1,1}^3}{b_{0,1}^2} - \frac{6b_{1,1}^2 b_{2,0}}{b_{1,0} b_{0,1}} - \frac{24b_{1,1} b_{2,0}^2}{b_{1,0}^2} + \frac{3b_{1,1} b_{2,1}}{b_{0,1}} + \frac{12b_{2,0} b_{2,1}}{b_{1,0}} + \frac{9b_{1,1} b_{3,0}}{b_{1,0}}$$

$$c_{2,2} = -3b_{2,2} - \frac{3b_{1,1}^3}{b_{0,1} b_{1,0}} - 9b_{1,1}^2 \left( \frac{b_{0,2}}{b_{0,1}^2} + \frac{b_{2,0}}{b_{1,0}^2} \right) + 6b_{1,1} \left( \frac{b_{1,2}}{b_{0,1}} - \frac{2b_{0,2}b_{2,0}}{b_{0,1}b_{1,0}} + \frac{b_{2,1}}{b_{1,0}} \right) + \frac{6b_{0,2}b_{2,1}}{b_{0,1}} + \frac{6b_{1,2}b_{2,0}}{b_{1,0}} \quad (10)$$

plus les équations obtenues en échangeant les deux états de spin ↑ et ↓. Nos expressions des $c_{n,0}$ sont en accord avec l'équation (10.33) de la référence [13], si ce n'est qu'elles l'étendent à un rapport $m_↑/\tilde{m}$ différent de un au travers du coefficient $b_{1,0} = \lambda^3/\lambda_↑^3$.

Avant de présenter nos résultats, nous devons choisir la masse de référence $\tilde{m}$. Pour ce faire, nous nous appuyons sur la forme particulière du développement en amas communément utilisée dans la littérature pour des fugacités égales, qui met en facteur de la série le nombre d'états de spin des fermions,

$$\frac{P\lambda^3}{k_B T} \stackrel{z_\sigma=z}{=} 2 \sum_{n=1}^{+\infty} b_n z^n \quad \text{avec} \quad b_n = \frac{1}{2} \sum_{n_1=0}^{n} b_{n_1, n-n_1} \qquad (11)$$

voir en particulier les références [5, 7] où $b_4$ est mesuré pour des masses égales<fnref n="6" />. Avoir sorti un tel facteur n'a d'intérêt que si le coefficient $b_1$ prend la valeur très simple égale à un. De la note 3 en bas de page, nous tirons donc le choix naturel<fnref n="7" />

$$b_1 = 1 \implies \tilde{m}^{3/2} = \frac{1}{2} (m_\uparrow^{3/2} + m_\downarrow^{3/2}) \quad \text{d'où} \quad \frac{1}{\lambda^3} = \frac{1}{2} \left( \frac{1}{\lambda_\uparrow^3} + \frac{1}{\lambda_\downarrow^3} \right) \qquad (12)$$

Ce choix rend également très simple le développement du viriel dans le cas de densités dans l'espace des phases égales dans les deux états de spin : on dispose alors de la belle expression $\rho_\sigma \lambda_\sigma^3 = \rho \lambda^3 / 2$ pour la valeur commune, où $\rho = \rho_\uparrow + \rho_\downarrow$ est la densité totale, et l'on pose comme dans l'équation (11) :

$$\frac{P\lambda^3}{k_B T} \stackrel{\rho_\sigma \lambda_\sigma^3 = \rho \lambda^3 / 2}{=} 2 \sum_{n=1}^{+\infty} c_n (\rho \lambda^3 / 2)^n \quad \text{avec} \quad c_n = \frac{1}{2} \sum_{n_1=0}^{n} c_{n_1, n-n_1} \qquad (13)$$

<footnote n="6">⁶Les valeurs de $b_4$ mesurées à l'ENS et au MIT sont en accord entre elles mais à un facteur ≈ 2 et respectivement à ≈ 2σ et à ≈ 3,5σ de la valeur désormais acceptée [45], où σ est l'incertitude expérimentale. Dans le cas de l'ENS, l'écart est dû au fait que toutes les données utilisables sur la pression sont de fugacité z > 1 et que leur extrapolation naïve à z = 0 échoue. En effet, la série (11) converge lentement et ne peut être tronquée avec une bonne approximation à n = 4 même pour une valeur aussi petite que z = 0,22, voir la version augmentée [47] de la publication [41] qui tire cette conclusion de l'équation d'état à haute précision obtenue par Monte-Carlo diagrammatique [11, 16].</footnote>

<footnote n="7">⁷Dans le contexte de la note 2, ce choix revient à prendre comme fonction de partition à un corps de référence $\tilde{Z}_1$ la moyenne arithmétique de $Z_{11}$ et $Z_{1↓}$. Il se réduit bien à $\tilde{m} = m$ dans le cas de masses égales $m_↓ = m_↑ = m$. Il assure que le facteur faisant passer du cas piégé au cas homogène dans l'équation (3) est une fonction uniformément bornée des masses (majorée par $2(n_↑^3 + n_↓^3)^{1/2}$ d'après l'inégalité de Hölder de paramètres p = 3/2 et q = 3). Ces conditions ne sont pas satisfaites si l'on prend pour $\tilde{m}$ la masse réduite de deux fermions ↑ et ↓ comme dans la référence [21].</footnote>
## 1.6. *Les résultats*

Nous représentons les quatrièmes coefficients d'amas du gaz unitaire non triviaux (différents de ceux du gaz parfait) en fonction du rapport de masse $\alpha$ sur la figure 1a (les valeurs $b_{4,0} = -(m_1/\bar{m})^{3/2}/32$ et $b_{0,4} = -(m_1/\bar{m})^{3/2}/32$ tirées de la note 3 ne sont donc pas représentées mais contribuent bien entendu à $b_4$). L'échelle logarithmique utilisée sur l'axe des abscisses met en évidence la symétrie $\alpha \leftrightarrow 1/\alpha$. Pour un rapport de masse $\alpha = 1$, nous retrouvons la valeur $b_4 = 0,030(1)$ conjecturée par la référence [41] et confirmée numériquement par la référence [45]. Lorsqu'on s'éloigne de ce point, le comportement de $b_4$ est d'abord dicté par $b_{2,2}$ qui le fait changer de signe (la courbe noire et la courbe verte sont presque parallèles), avant que $b_{3,1}$ ou $b_{1,3}$ ne l'emporte aux rapports de masse élevés $\alpha$ ou $1/\alpha$ et ne fasse remonter $b_4$ à des valeurs grandes et positives. Aux seuils $\alpha_c^{\text{4corps}}$ ou $1/\alpha_c^{\text{4corps}}$ de l'effet Efimov à 3+1 ou 1+3 corps, repérés par des lignes verticales pointillées, $b_4$ admet une limite finie mais une dérivée infinie comme $b_{3,1}$ ou $b_{1,3}$, voir la section 2 ci-dessous et la note 4 en bas de page. En revanche, $b_{2,2}$ reste une fonction lisse, et ne présenterait une dérivée infinie qu'aux seuils de l'effet Efimov à trois corps, voir la section 3. Sur la figure 1b, on constate que les quatrièmes coefficients du viriel non triviaux (autres que $c_{4,0} = A_4(m_1/\bar{m})^{3/2}$ et $c_{0,4} = A_4(m_1/\bar{m})^{3/2}$, avec $A_4 = (18+15\sqrt{2}-16\sqrt{6})/192 \approx 1,11 \times 10^{-4}$) ont une structure similaire même s'ils sont de signe constant. On remarque aussi que le coefficient $c_4$ dépend peu ou très peu des corrélations à quatre corps induites par les interactions : une approximation de $c_4$ négligeant les contributions des coefficients d'amas $b_{3,1}$, $b_{2,2}$ et $b_{1,3}$, représentée en tireté sur la figure 1b, est partout proche ou très proche de $c_4$. Même si la possibilité assez récente de préparer des gaz d'atomes froids spatialement homogènes dans une boîte de potentiel à fond plat [48,49] rend la mesure des coefficients du viriel $c_n$ très naturelle, puisque les densités $\rho_\sigma$ y sont directement accessibles, cette dernière prédiction la rend moins motivante que la mesure des coefficients d'amas $b_n$ dans un gaz inhomogène selon la technique spécifique au cas harmoniquement piégé de la référence [5], dans laquelle ce sont les potentiels chimiques $\mu_\sigma$ qui s'imposent comme variables pertinentes<fnref n="8" />. À toutes fins utiles, nous donnons aussi les quatrièmes coefficients d'amas et du viriel sous forme numérique dans la table 1.

## 1.7. *Plan du reste de l'article*

Dans la suite, nous expliquons comment nous avons pu obtenir des résultats précis en un temps de calcul raisonnable, en particulier au moyen d'une accélération de convergence de la somme sur $\ell$ dans l'expression (4) s'appuyant sur un développement asymptotique du sommande et jouant un rôle essentiel près des seuils efimoviens. Il nous a fallu pour cela généraliser la méthode analytique de la référence [40], mise en œuvre pour 3 + 1 corps, au cas beaucoup plus difficile de 2 + 2 corps. Dans le cas à 3 + 1 corps, nous mettons aussi en évidence une dépendance non monotone inattendue de $B_{3,1}(\bar{\omega})$ en $\bar{\omega}$, avec changement de signe, pour des rapports de masse $\alpha$ proches de 3,5 ou 6,6. La section 2 traite du cas à 3 + 1 corps et la section 3 du cas à 2 + 2 corps dans un potentiel harmonique.

<footnote n="8">⁸Dans un potentiel de piégeage $U(r)$ commun aux deux états de spin et très allongé selon l'axe propre $Oz$, les gaz d'atomes froids sont bien décrits par l'approximation d'homogénéité locale, et ont donc des potentiels chimiques locaux sur l'axe bien défini $\mu_\sigma^{\text{loc}}(z) = \mu_\sigma - U(0, 0, z)$, ce qui rend très directe la mesure de l'équation d'état grand-canonique, donc des coefficients d'amas, en suivant la proposition astucieuse de la référence [50], qui est utilisable cependant seulement dans le cas où $U(r)$ est harmonique.</footnote>
<figure><img src="image_3.png" /></figure>

**FIGURE 1** – Pour un gaz unitaire spatialement homogène de fermions à deux états de spin $\uparrow$ et $\downarrow$, quatrièmes coefficients d'amas (a) et du viriel (b) définis par les équations ((1), (11), (13)) et prédits par la conjecture (5) de la référence [41], en fonction du rapport de masse $\alpha = m_1/m_1$. L'interaction $\uparrow\downarrow$ de portée nulle est supposée être invariante d'échelle dans le problème à quatre corps, ce qui oblige à se limiter aux rapports de masse compris entre les valeurs critiques de l'effet Efimov à 3 + 1 et à 1 + 3 corps, indiquées par les pointillés verticaux. En (a), les traits verticaux fins repèrent les points où $b_4$ s'annule (toujours avec changement de signe), à savoir $\alpha \approx 1,535$ et $\alpha \approx 10,355$ du côté $\alpha > 1$, et la vignette est un agrandissement (les points calculés y sont représentés par des disques sur des lignes d'interpolation). En (b), la courbe tiretée néglige dans $c_4$ la composante véritablement à quatre corps issue des interactions, celle $c_4^{\text{corr}} = -3(b_{1,3} + b_{2,2} + b_{3,1})/2$ qui dépend de $b_{3,1}$, $b_{2,2}$ ou $b_{1,3}$, pour représenter $c_4^{\text{sans corr}} = c_4 - c_4^{\text{corr}}$.

## 2. Coefficient d'amas à 3 + 1 fermions dans un piège

Nous expliquons dans cette section comment calculer numériquement, mais de manière efficace, la quantité $I_{3,1}(\bar{\omega})$ définie par l'équation (4), pour une valeur quelconque (nulle ou strictement positive) de la pulsation de piégeage réduite $\bar{\omega} = \hbar\omega/k_B T$. Le coefficient d'amas $B_{3,1}(\bar{\omega})$ cherché en découle directement d'après la conjecture (5).

### 2.1. Formulation du problème

Rappelons d'abord l'expression de l'opérateur $M_{3,1}^{(\ell, \epsilon)}(iS)$ intervenant dans l'équation (4), telle qu'elle apparaît dans la référence [39]. Nous avons affaire à la somme d'une partie diagonale $\mathcal{D}_{1,3}$ et d'un opérateur à noyau $\mathcal{K}_{3,1}$ agissant sur les fonctions $f_{m_z}(x, u)$ de deux variables continues, le logarithme du rapport des normes de deux vecteurs d'onde $x \in \mathbb{R}^+$ et le cosinus de l'angle entre eux $u = \cos\theta \in [-1, 1]$, et d'une variable discrète, le nombre quantique magnétique $m_z$ selon l'axe de quantification $Oz$, variant par pas de deux entre $-\ell$ et $\ell$ pour la parité $\epsilon = (-1)^{\ell}$, et entre $-\ell + 1$ et $\ell - 1$ pour la parité $\epsilon = (-1)^{\ell-1}$ et $\ell \neq 0$. Il vient dans une notation mixte de Dirac et de Schrödinger :

$$ \begin{aligned} & \langle x, u | \langle \ell, m_z | M_{3,1}^{(\ell, \epsilon)} (iS) | f \rangle \\ & = \mathcal{D}_{3,1}(x, u) f_{m_z}(x, u) + \int_0^{+\infty} dx' \int_{-1}^{1} du' \sum_{\substack{m'_z | (-1)^{m'_z} = \epsilon \\ m'_z \neq 0}} \mathcal{K}_{3,1}^{(\ell)}(x, u, m_z; x', u', m'_z) f_{m'_z}(x', u') \end{aligned} \quad (14) $$
<table><caption>**TABLE 1** – Valeurs numériques des quatrièmes coefficients d'amas $b_{n,4-n}$ et du viriel $c_{n,4-n}$ du gaz unitaire spatialement homogène de fermions à deux états de spin ↑ et ↓, tabulées en fonction du rapport de masse $\alpha = m_↓/m_↑$, pour le choix (12) de la masse de référence $\tilde{m}$. Comme on passe de $b_{n,4-n}$ à $b_{4-n,n}$ et de $c_{n,4-n}$ à $c_{4-n,n}$ en changeant $\alpha$ en $1/\alpha$, nous nous sommes limités à $\alpha \ge 1$. Nous ne donnons pas les coefficients associés aux entiers $n = 0$ et $n = 4$ car ils sont identiques à ceux du gaz parfait, mais nous donnons les coefficients moyens $b_4$ et $c_4$ utiles dans le cas de fugacités ou de densités dans l'espace des phases égales dans les deux états de spin, voir les équations (11) et (13). Les incertitudes sur $b_{3,1}$, $b_{1,3}$ et $b_{2,2}$, non précisées, sont inférieures au pour cent. Les valeurs données sont représentées graphiquement sur la figure 1.</caption><thead><tr><th>α</th><th>1</th><th>1,5</th><th>2</th><th>2,5</th><th>3</th><th>3,5</th><th>4</th><th>4,5</th><th>5</th><th>5,5</th></tr></thead><tbody><tr><td>b<sub>3,1</sub></td><td>0,1837</td><td>0,1604</td><td>0,1247</td><td>0,08238</td><td>0,03726</td><td>-0,007203</td><td>-0,04762</td><td>-0,07982</td><td>-0,09941</td><td>-0,1016</td></tr><tr><td>b<sub>1,3</sub></td><td>0,18374</td><td>0,18672</td><td>0,18199</td><td>0,17641</td><td>0,17133</td><td>0,16696</td><td>0,16325</td><td>0,16010</td><td>0,15740</td><td>0,15508</td></tr><tr><td>b<sub>2,2</sub></td><td>-0,2445</td><td>-0,2781</td><td>-0,3469</td><td>-0,4321</td><td>-0,5283</td><td>-0,6328</td><td>-0,7448</td><td>-0,8637</td><td>-0,9890</td><td>-1,1209</td></tr><tr><td>b<sub>4</sub></td><td>0,03026</td><td>0,003256</td><td>-0,05135</td><td>-0,1179</td><td>-0,19110</td><td>-0,2678</td><td>-0,3458</td><td>-0,4229</td><td>-0,4968</td><td>-0,5649</td></tr><tr><td>c<sub>3,1</sub></td><td>-2,2558</td><td>-4,4056</td><td>-7,3318</td><td>-11,100</td><td>-15,773</td><td>-21,407</td><td>-28,050</td><td>-35,737</td><td>-44,495</td><td>-54,332</td></tr><tr><td>c<sub>1,3</sub></td><td>-2,2558</td><td>-1,2265</td><td>-0,8296</td><td>-0,6292</td><td>-0,5114</td><td>-0,4350</td><td>-0,3820</td><td>-0,3433</td><td>-0,3139</td><td>-0,2910</td></tr><tr><td>c<sub>2,2</sub></td><td>-7,1913</td><td>-7,3897</td><td>-7,7991</td><td>-8,3081</td><td>-8,8748</td><td>-9,4800</td><td>-10,111</td><td>-10,758</td><td>-11,415</td><td>-12,078</td></tr><tr><td>c<sub>4</sub></td><td>-5,8513</td><td>-6,5108</td><td>-7,9801</td><td>-10,018</td><td>-12,580</td><td>-15,661</td><td>-19,271</td><td>-23,419</td><td>-28,112</td><td>-33,350</td></tr><tr><td>α</td><td>6</td><td>6,5</td><td>7</td><td>7,5</td><td>8</td><td>8,5</td><td>9</td><td>9,5</td><td>10</td><td>10,5</td></tr><tr><td>b<sub>3,1</sub></td><td>-0,07925</td><td>-0,02715</td><td>0,06412</td><td>0,2022</td><td>0,4010</td><td>0,6724</td><td>1,0339</td><td>1,5059</td><td>2,1195</td><td>2,9095</td></tr><tr><td>b<sub>1,3</sub></td><td>0,15305</td><td>0,15128</td><td>0,14972</td><td>0,14833</td><td>0,14708</td><td>0,14597</td><td>0,14496</td><td>0,14404</td><td>0,14321</td><td>0,14244</td></tr><tr><td>b<sub>2,2</sub></td><td>-1,2590</td><td>-1,4034</td><td>-1,5542</td><td>-1,7113</td><td>-1,8749</td><td>-2,0452</td><td>-2,2223</td><td>-2,4077</td><td>-2,6003</td><td>-2,8015</td></tr><tr><td>b<sub>4</sub></td><td>-0,6239</td><td>-0,6709</td><td>-0,7014</td><td>-0,7116</td><td>-0,6947</td><td>-0,6447</td><td>-0,5530</td><td>-0,4101</td><td>-0,2001</td><td>0,09397</td></tr><tr><td>c<sub>3,1</sub></td><td>-65,248</td><td>-77,214</td><td>-90,190</td><td>-104,10</td><td>-118,85</td><td>-134,28</td><td>-150,20</td><td>-166,33</td><td>-182,32</td><td>-197,64</td></tr><tr><td>c<sub>1,3</sub></td><td>-0,2725</td><td>-0,2574</td><td>-0,2449</td><td>-0,2343</td><td>-0,2252</td><td>-0,2174</td><td>-0,2105</td><td>-0,2045</td><td>-0,1991</td><td>-0,1944</td></tr><tr><td>c<sub>2,2</sub></td><td>-12,740</td><td>-13,397</td><td>-14,044</td><td>-14,676</td><td>-15,288</td><td>-15,873</td><td>-16,423</td><td>-16,926</td><td>-17,374</td><td>-17,747</td></tr><tr><td>c<sub>4</sub></td><td>-39,130</td><td>-45,434</td><td>-52,239</td><td>-59,505</td><td>-67,180</td><td>-75,184</td><td>-83,417</td><td>-91,733</td><td>-99,948</td><td>-107,79</td></tr><tr><td>α</td><td>11</td><td>11,5</td><td>12</td><td>12,5</td><td>12,75</td><td>13</td><td>13,1</td><td>13,2</td><td>13,3</td><td>13,3842</td></tr><tr><td>b<sub>3,1</sub></td><td>3,9362</td><td>5,2839</td><td>7,0993</td><td>9,7004</td><td>11,522</td><td>14,042</td><td>15,402</td><td>17,152</td><td>19,728</td><td>26,101</td></tr><tr><td>b<sub>1,3</sub></td><td>0,14174</td><td>0,14109</td><td>0,14049</td><td>0,13993</td><td>0,13967</td><td>0,13941</td><td>0,13932</td><td>0,13922</td><td>0,13912</td><td>0,13904</td></tr><tr><td>b<sub>2,2</sub></td><td>-3,0120</td><td>-3,2360</td><td>-3,4725</td><td>-3,7278</td><td>-3,8654</td><td>-4,0133</td><td>-4,0767</td><td>-4,1430</td><td>-4,2140</td><td>-4,2784</td></tr><tr><td>b<sub>4</sub></td><td>0,5017</td><td>1,0632</td><td>1,8524</td><td>3,0250</td><td>3,8667</td><td>5,0528</td><td>5,7008</td><td>6,5428</td><td>7,7951</td><td>10,950</td></tr><tr><td>c<sub>3,1</sub></td><td>-211,58</td><td>-223,04</td><td>-230,30</td><td>-230,37</td><td>-225,99</td><td>-217,03</td><td>-211,71</td><td>-205,24</td><td>-198,00</td><td>-200,77</td></tr><tr><td>c<sub>1,3</sub></td><td>-0,1901</td><td>-0,1863</td><td>-0,1828</td><td>-0,1796</td><td>-0,1781</td><td>-0,1766</td><td>-0,1761</td><td>-0,1756</td><td>-0,1750</td><td>-0,1746</td></tr><tr><td>c<sub>2,2</sub></td><td>-18,025</td><td>-18,163</td><td>-18,114</td><td>-17,767</td><td>-17,414</td><td>-16,864</td><td>-16,561</td><td>-16,189</td><td>-15,720</td><td>-15,213</td></tr><tr><td>c<sub>4</sub></td><td>-114,90</td><td>-120,69</td><td>-124,30</td><td>-124,16</td><td>-121,79</td><td>-117,04</td><td>-114,22</td><td>-110,80</td><td>-106,95</td><td>-108,08</td></tr></tbody></table>

La partie diagonale est indépendante du moment cinétique et de l'exposant d'échelle iS, et le noyau $\mathcal{K}_{3,1}$ s'obtient par symétrisation d'un noyau primitif $K_{3,1}$ car on a tiré parti de la symétrie d'échange fermionique du vecteur d'état pour se restreindre à $x > 0$ :<fnref n="9" />

$$ \mathcal{D}_{3,1}(x, u) = \frac{(1 + 2\alpha + \alpha u / \text{ch } x)^{1/2}}{1 + \alpha} ; \qquad (15) \\ \mathcal{K}_{3,1}^{(\ell)}(x, u, m_z; x', u', m_z') = \sum_{\eta, \eta' = \pm 1} (\eta\eta')^{\ell+1} K_{3,1}^{(\ell)}(\eta x, u, \eta m_z; \eta' x', u', \eta' m_z') $$

<footnote n="9">⁹ En notation de Dirac, $\mathcal{K}_{3,1}^{(\ell)} = (1+U)K_{3,1}^{(\ell)}(1+U)$ avec l'involution hermitienne $U = -P_x e^{i\pi L_x/\hbar}$, l'opérateur de parité $P_x$ changeant $|x\rangle$ en $|-x\rangle$ et l'opérateur de rotation d'axe $Ox$ d'angle $\pi$ tel que $\exp(i\pi L_x/\hbar)|\ell, m_z\rangle = (-1)^{\ell}|\ell, -m_z\rangle$ [40].</footnote>
Les éléments de matrice du noyau primitif dépendent de iS et contiennent une intégrale sur un angle de rotation $\phi$ autour de l'axe Ox, $L_x$ étant l'opérateur moment cinétique selon cette direction :

$$ K_{3,1}^{(\ell)}(x, u, m_z; x', u', m_z') = \frac{(\lambda \lambda')^{3/2}}{[(1 + \lambda^2)(1 + \lambda'^2)]^{1/4}} \int_0^{2\pi} \frac{d\phi}{2\pi^2} \frac{(1 + \lambda^2)^{iS/2} e^{-im_z\theta/\hbar} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m_z' \rangle e^{im_z'\theta'/2} (1 + \lambda'^2)^{-iS/2}}{1 + \lambda^2 + \lambda'^2 + \frac{2\alpha}{1+\alpha} [\lambda u + \lambda' u' + \lambda\lambda'(uu' + vv' \cos\phi)]} \quad (16) $$

où, comme dans la référence [39], nous avons posé pour abréger :

$$ \lambda = e^{x}, \quad \lambda' = e^{x'}, \quad \theta = \arccos u \in [0, \pi], \quad \theta' = \arccos u' \in [0, \pi], \quad v = \sin \theta, \quad v' = \sin \theta' \quad (17) $$

Pour évaluer $I_{3,1}(\bar{\omega})$ sur ordinateur, il faut d'abord remplacer l'opérateur par une matrice de taille finie, en tronquant la variable $x$ à $x_{max}$ et en la discrétisant suivant la méthode d'intégration du point milieu, puis en discrétisant la variable $\theta$ (que nous préférons à la variable $u$ car elle conduit à un intégrande lisse) suivant la méthode d'intégration de Gauss-Legendre. On calcule ensuite le déterminant de la matrice en la mettant sous la forme de Cholesky, pour tirer parti du fait que l'opérateur $M_{3,1}^{(\ell, \epsilon)}(iS)$ est défini positif en l'absence d'effet Efimov à 3 + 1 corps. Enfin, on calcule l'intégrale sur S dans l'intervalle $[0, S_{max}]$ par la méthode du point milieu (il faut connaître le logarithme du déterminant aux multiples entiers du pas d'intégration dS pour obtenir sa dérivée aux multiples demi-entiers) et on tient compte de la contribution de l'intervalle omis $[S_{max}, +\infty[$ au moyen d'une approximation exponentielle $A\exp(-BS)$ de la dérivée logarithmique du déterminant, justifiée par la référence [40] et dont les coefficients A et B sont estimés par ajustement sur un voisinage de $S_{max}$, ici l'intervalle $[S_{max} - 5/2, S_{max}]$.

## 2.2. Approximant asymptotique et applications

Il reste à tenir compte de la troncature sur le moment cinétique $\ell$ à un certain $\ell_{max}$, inévitable dans un calcul numérique. En pratique, il n'est pas raisonnable d'aller au-delà de $\ell_{max} = 15$, les valeurs élevées de $\ell$ conduisant à une complexité $O(\ell^3)$ et étant très coûteuses en temps. Malheureusement, cette coupure ne suffit pas pour de grands rapports de masse $\alpha \gg 1$ si l'on vise une précision sur $B_{3,1}(0^+)$ meilleure que le pour cent. L'idée salvatrice est de déterminer un approximant asymptotique $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ de la contribution $I_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ de moment cinétique $\ell$ et de parité $\epsilon$ à la quantité cherchée $I_{3,1}(\bar{\omega})$. Alors, au lieu de négliger complètement les termes $I_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ pour $\ell > \ell_{max}$, on les remplace par $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ comme suit,

$$ I_{3,1}(\bar{\omega}) \approx \sum_{\ell=0}^{\ell_{max}} \sum_{\epsilon} I_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) + \sum_{\ell=\ell_{max}+1}^{+\infty} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) = \sum_{\ell=0}^{\ell_{max}} \sum_{\epsilon} \left[ I_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) - J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) \right] + \sum_{\ell=0}^{+\infty} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) \quad (18) $$

ce qui constitue une méthode d'accélération de convergence de la série, l'erreur commise tendant vers zéro plus rapidement avec $\ell_{max}$. On obtient un approximant asymptotique exact à l'ordre sous-dominant en $\ell$ en généralisant à $\bar{\omega} \neq 0$ l'application de la méthode de la référence [40], qui prend le noyau $\mathcal{K}_{3,1}^{(\ell, \epsilon)}$ comme petit paramètre formel et développe le logarithme du déterminant de $M_{3,1}^{(\ell, \epsilon)}$ au second ordre :

$$ \begin{aligned} \ln \det M_{3,1}^{(\ell, \epsilon)} &= \ln \det^{\epsilon} [\mathcal{D}_{3,1} + \mathcal{K}_{3,1}^{(\ell)}] \\ &= \ln \det^{\epsilon} [\mathcal{D}_{3,1}] + \ln \det^{\epsilon} [\mathbb{1} + \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)}] \\ &= \text{cte} + \text{Tr}^{\epsilon} \left[ \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} - \frac{1}{2} \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} + \dots \right] \end{aligned} \quad (19) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Contraste entre la contribution $I_{3,1}^{(\ell, \epsilon)}$ de la voie de moment cinétique $\ell$ et de parité $\epsilon$ à $I_{3,1}$ et son approximant asymptotique $J_{3,1}^{(\ell, \epsilon)}$ donné par l'équation (47), en fonction du rapport de masse $\alpha = m_1/m_\ell$ et de $\ell$, dans la limite $\tilde{\omega} = 0^+$. On entend ici par contraste entre deux quantités réelles $a$ et $b$ le rapport $|a-b|/(|a|+|b|)$. (a) Cas $\epsilon = (-1)^\ell$. (b) Cas $\epsilon = (-1)^{(\ell+1)}$ (ce qui impose $\ell > 0$). Les paramètres de troncature et de discrétisation numériques sont ceux de la figure 3. Le contraste sature à un (comme ceci se produit pour $\ell = 0$ et $\ell = 5$ sur (a), pour $\ell = 3$ sur (b)) lorsque les deux quantités comparées sont de signes opposés.</figcaption></figure>

où le symbole $\epsilon$ en exposant de la trace ou du déterminant signifie qu'on se limite au sous-espace des valeurs de $m_\alpha$ compatibles avec la parité $\epsilon$. Le calcul, un peu long, est effectué dans l'annexe A, et l'expression correspondante de $J_{3,1}^{(\ell, \epsilon)}(\tilde{\omega})$ est donnée dans l'équation (47) sous forme d'intégrales multiples. Pour ce choix déjà très élaboré, les premières valeurs de l'approximant $(0 \le \ell \le \ell_{max})$ restent faciles à calculer numériquement avec la même troncature et discrétisation que pour le déterminant complet, mais il serait pénible d'aller au-delà, ce qu'exige pourtant la seconde somme au deuxième membre de l'équation (18); heureusement, la série infinie au troisième membre de l'équation (18) admet une expression intégrale simple, voir l'équation (59), que l'on peut même exprimer analytiquement pour $\tilde{\omega} = 0^+$ en termes de fonctions connues comme la fonction dilogarithme, voir l'équation (60). Le lecteur pourra juger de la qualité de notre approximant asymptotique sur la figure 2, où il est comparé au résultat numérique. La figure montre, en fonction du rapport de masse $\alpha$, dans quelles voies de moment cinétique $\ell$ l'approximant dévie de plus d'un pour cent de la valeur exacte (c'est la précision cible sur le coefficient d'amas); ces voies doivent absolument être incluses dans la somme numérique de $0 \text{ à } \ell_{max}$.

Une retombée intéressante de notre méthode d'accélération de convergence est l'obtention d'un équivalent asymptotique de la contribution de moment cinétique $\ell$ et de parité $\epsilon$ au coefficient d'amas $B_{3,1}(0^+)$; il suffit de garder la contribution linéaire en $\mathcal{K}_{3,1}^{(\ell)}$ dans l'équation (19) et de déterminer son comportement dominant à grand $\ell$ avec le théorème intégral de Cauchy, comme il est fait dans l'annexe A. On trouve :

$$I_{3,1}^{(\ell, \epsilon)}(0^{+}) \sim \frac{1}{2} \left\{ \frac{(1+\alpha)^2}{\sqrt{\alpha(1+3\alpha)}} \operatorname{Im} \left[ \left( \frac{\ell}{6\pi C_0} \right)^{1/2} \frac{z_0}{1+z_0} z_0^{\ell} \right] + \epsilon \frac{(1+\alpha)^2}{\alpha^{3/2}} \operatorname{Re} \left[ \left( \frac{\ell}{2\pi C_1} \right)^{1/2} \frac{(1-z_1)}{\sqrt{1-\cos\xi_1}} z_1^{\ell} \right] \right\} \quad (20)$$
où

$$
\begin{equation}
\begin{aligned}
C_0 &= \frac{(1 - \cos \xi_0) (5/4 + \cos \xi_0)}{(1/2 + \cos \xi_0) (z_0 - 1/z_0)} \\
\text{et } C_1 &= 4 (z_1 - 1/z_1) (\cos \xi_1 - 1/2) \left[ (1 + \alpha^{-1}) \cos \xi_1 + \frac{1 + 3\alpha^{-1}}{4} \right]
\end{aligned}
\tag{21}
\end{equation}
$$

avec le raccourci de notation cosξₙ ≡ (zₙ + 1/zₙ)/2, n = 0 ou 1, le nombre complexe z₀ donné par l'équation (66) et le nombre complexe z₁ = -z₀*. C'est la généralisation à 3 + 1 fermions d'un résultat obtenu pour trois bosons dans la référence [20], voir l'équation (42) de cette référence. Ainsi, à un facteur en loi de puissance près, I<sub>3,1</sub><sup>(ℓ,ε)</sup> (0<sup>+</sup>) tend exponentiellement vers zéro avec ℓ, avec des oscillations irrégulières dues au fait que l'argument de z₀ n'est en général pas commensurable à π.

2.3. Régime de Born–Oppenheimer

Dans la limite $\alpha \to 0$, on remarque que $z_0$ tend vers zéro dans l'équation (66), $z_0 \sim i\sqrt{\alpha/3}$, donc on s'attend à ce que la suppression exponentielle du sommande dans (4) devienne très rapide pour 3+1 fermions. Cette attente est confirmée numériquement et, pour $\alpha = 0$, seule la voie $\ell = 0$ contribue. Dans ce cas, l'unique fermion $\downarrow$ est infiniment massif et se comporte pour les fermions $\uparrow$ comme un diffuseur ponctuel fixe, de longueur de diffusion infinie dans l'onde $s$ et placé au centre du piège. L'approximation de Born-Oppenheimer devient exacte et donne l'équation de Schrödinger indépendante du temps sur la fonction d'onde $\Psi(\mathbf{r}_\downarrow)$ de la particule lourde<fnref n="10" /> :

$$
E_{\text{BO}} \Psi(\mathbf{r}_{\downarrow}) = \left[ -\frac{\hbar^2}{2m_{\downarrow}} \Delta_{\mathbf{r}_{\downarrow}} + \frac{1}{2} m_{\downarrow} \omega^2 r_{\downarrow}^2 + W(\mathbf{0}) \right] \Psi(\mathbf{r}_{\downarrow}) \quad (m_{\downarrow} \to +\infty) \quad (22)
$$

Dans le problème piégé à $n_\uparrow + 1$ fermions, $W(\mathbf{0})$ est l'énergie d'un état propre de $n_\uparrow$ fermions $\uparrow$
sans interaction en présence du centre diffuseur. Un fermion $\uparrow$ de nombres quantiques orbitaux
$(n, \ell, m_z)$ ne voit le centre diffuseur que si $\ell = 0$, auquel cas son spectre est abaissé de $\hbar\omega$, donc
admet les niveaux d'énergie

$$
\epsilon_{n,\ell,m_z} = \begin{cases} (2n + \ell + 3/2)\hbar\omega & \text{si } \ell > 0 \\ (2n + 1/2)\hbar\omega & \text{si } \ell = 0 \end{cases} \quad (n \in \mathbb{N}, -\ell \le m_z \le \ell) \quad (23)
$$

Comme l'énergie propre de Born-Oppenheimer E<sub>BO</sub> est somme de W(0) et d'un niveau d'énergie vibrationnel d'une particule ↓ dans le piège, nous concluons que

$$
\lim_{\alpha \to 0^+} Z_{n_1, 1} = Z_1 Z_{n_1}^{\text{diff}} \qquad (24)
$$

où $Z_{n_1}^{\text{diff}}$ est la fonction de partition canonique d'un gaz parfait de $n_1$ fermions de spectre (23), en présence donc du diffuseur fixe, $Z_1$ est comme dans l'équation (2) la fonction de partition d'un fermion seul et $Z_{n_1,1}$ est celle du gaz unitaire de $n_1+1$ fermions piégés. Il est commode de calculer son écart $\Delta Z_{n_1,1}$ à celle du gaz parfait de $n_1+1$ fermions piégés, puisque (23) diffère du spectre ordinaire seulement dans la voie de moment cinétique nul. En tenant compte de la statistique de

<footnote n="10">¹⁰Les particules sont à température fixée *T* donc, dans la limite *m*<sub>↓</sub> → +∞, le fermion ↓ occupe une région autour du centre du piège de rayon *R* = *O*((*k*<sub>B</sub>*T*)/*m*<sub>↓</sub>*ω*<sup>2</sup>)<sup>1/2</sup>). En revanche, le potentiel de Born–Oppenheimer *W*(*r*<sub>↓</sub>) a une échelle d’énergie *k*<sub>B</sub>*T* et varie avec une échelle de longueur *λ*<sub>↑</sub>, la longueur d’onde thermique de Broglie des fermions légers, si *k*<sub>B</sub>*T* > *ħ*ω, et est de l’ordre de *ħ*ω et varie avec une échelle de longueur (*ħ*/*m*<sub>↑</sub>*ω*)<sup>1/2</sup>, la taille de l’état vibrationnel fondamental d’un fermion léger, sinon. La variation de *W*(*r*<sub>↓</sub>) devient négligeable dans tous les cas et on peut le remplacer par *W*(*O*). Pour la même raison, le potentiel dit scalaire ou topologique, qui s’ajoute au potentiel *W* dans l’approximation adiabatique complète [51], peut être omis.</footnote>
Fermi via le principe d'exclusion de Pauli, et [41, équations (80)–(83)] qui relient les coefficients d'amas du système piégé aux fonctions de partition à petit nombre de corps, nous trouvons que

$$B_{1,1}(\bar{\omega}) = Z_1^{-1} \Delta Z_{1,1} \underset{\alpha \to 0^+}{\to} \frac{1}{2 \operatorname{ch}(\bar{\omega}/2)} ; \quad (25)$$

$$B_{2,1}(\bar{\omega}) = Z_1^{-1} \Delta Z_{2,1} - Z_1 B_{1,1} \underset{\alpha \to 0^+}{\to} - \frac{1}{8 \operatorname{ch}^2(\bar{\omega}/2) \operatorname{ch} \bar{\omega}}$$

$$B_{3,1}(\bar{\omega}) = Z_1^{-1} \Delta Z_{3,1} - Z_{2,0} B_{1,1} - Z_1 B_{2,1} \underset{\alpha \to 0^+}{\to} \frac{1}{16 \operatorname{ch}^3(\bar{\omega}/2) (\operatorname{ch} \bar{\omega})(2 \operatorname{ch} \bar{\omega} - 1)} \quad (26)$$

La valeur (26) de la limite est en accord parfait avec notre calcul numérique de l'expression (4), non seulement pour $\bar{\omega} = 0^+$, où $B_{3,1}(0^+) \to 1/16 = 0,0625$, comme on le voit sur la figure 3a, mais, comme nous l'avons vérifié, pour tout $\bar{\omega}$, ce qui constitue un test supplémentaire de la conjecture (5) dans le cas à 3 + 1 corps.<fnref n="11" /> Signalons en passant qu'il est possible d'aller aux ordres suivants en utilisant la version grand-canonique de l'équation (24),

$$\Xi_{z_1 \to 0 \text{ puis } \alpha \to 0^+} = \Xi_\uparrow + z_1 Z_1 \Xi_\downarrow^{\text{diff}} + O(z_1^2) \quad (27)$$

où $\Xi_\downarrow (\Xi_\downarrow^{\text{diff}})$ est la grande fonction de partition du gaz parfait de fermions $\uparrow$ de fugacité $z_\downarrow$ en absence (en présence) du centre diffuseur et $\Xi$ est celle du gaz unitaire à deux composantes. En développant le grand potentiel $\Omega = -k_B T \ln \Xi$ comme dans l'équation (2), nous obtenons<fnref n="12" />

$$\sum_{n_1=0}^{+\infty} z_1^{n_1} \lim_{\alpha \to 0^+} B_{n_1,1}(\bar{\omega}) = \frac{\Xi_\downarrow^{\text{diff}}}{\Xi_\uparrow} = \prod_{s=0}^{+\infty} \frac{1+z_\downarrow e^{-(2n+1/2)\bar{\omega}}}{1+z_\downarrow e^{-(2n+3/2)\bar{\omega}}} = \exp \left[ \sum_{s=1}^{+\infty} \frac{(-1)^{s+1} z_\downarrow^s}{2s \operatorname{ch}(s\bar{\omega}/2)} \right] \quad (28)$$

Dans le cas particulier $\bar{\omega} = 0^+$, l'exponentielle au quatrième membre de (28) se réduit à $\sqrt{1+z_\downarrow}$, facile à développer en puissances de $z_\downarrow$, d'où le résultat à tous les ordres en termes de la fonction $\Gamma$ d'Euler :

$$\lim_{\alpha \to 0^+} B_{n_1,1}(0^+) = \frac{(-1)^{n_1} \Gamma(n_1 - 1/2)}{n_1! \Gamma(-1/2)} \quad \forall n_1 \in \mathbb{N} \quad (29)$$

## 2.4. Les résultats

Nos résultats numériques pour $I_{3,1}(0^+)$, donc pour le coefficient d'amas $B_{3,1}(0^+)$ d'après l'équation (5), sont donnés explicitement dans la table 2 et sont représentés en fonction du rapport de masse sur la figure 3a, sans barres d'erreur mais avec une incertitude inférieure au pour cent. Nous trouvons que $B_{3,1}(0^+)$ est positif, sauf sur un intervalle $[\alpha_1, \alpha_2]$, comme on le voit mieux sur l'agrandissement inclus dans la figure. Aux grandes valeurs de $\alpha$, $B_{3,1}(0^+)$ présente une croissance marquée que nous attribuons à l'effet Efimov à trois corps, mais admet une limite finie au seuil de l'effet Efimov à 3 + 1 corps où doit s'arrêter notre calcul, atteinte avec une pente infinie et que nous déterminons par extrapolation, voir les courbes en tireté sur la figure 3a et les explications qui suivent,

$$B_{3,1}(0^+) \underset{\alpha \to (\alpha_c^{\text{corps}})^{-}}{\to} 2,47 \pm 0,03 \quad (30)$$

<footnote n="11">¹¹Dans le cas à 2 + 1 corps, en partant de l'expression analytique intégrale de $B_{2,1}(\bar{\omega})$ de la référence [21] en termes de la fonction transcendante d'Efimov $\Lambda_{2,1}^{(\ell)}$ sur l'axe imaginaire pur, et en montrant dans la limite $\alpha \to 0^+$ que $\Lambda_{2,1}^{(\ell)}(iS)$ tend vers 1 pour $\ell > 0$ et tend vers $1+1/\operatorname{ch}(\pi S/2)$ pour $\ell = 0$ (par exemple au moyen de [38, équation (46)]), on retrouve le second résultat dans l'équation (25).</footnote>

<footnote n="12">¹²On a fait apparaître une somme sur un entier $s$ dans le quatrième membre en prenant le logarithme du produit infini dans le troisième membre, en développant autour de 1 la fonction $\ln$ en série entière d'indice $s$ puis en sommant sur $n$ les séries géométriques qui en résultent.</footnote>
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** – (a) Quatrième coefficient d'amas $B_{3,1}(0^+)$ du gaz unitaire de fermions dans un piège harmonique de raideur infinitésimale, en fonction du rapport de masse $\alpha = m_\uparrow/m_\downarrow$. Disques (reliés par un trait fin) : notre calcul numérique de l'expression (4) pour $(n_\uparrow, n_\downarrow) = (3,1)$ et application de la conjecture (5), avec les paramètres de troncature $x_{max} = 5$, $S_{max} = 25$, $\ell_{max} = 15$ et de discrétisation $dx = 1/10$, $dS = 1/8$, $n_\theta = 37$ (c'est le nombre de valeurs de $\theta$ dans la méthode d'intégration de Gauss-Legendre); pour $\alpha > 10$ et $0 \le \ell \le 5$, nous prenons plutôt $x_{max} = 20$ et $n_\theta = 25$ pour maintenir l'erreur inférieure au pour cent. On utilise la méthode d'accélération de convergence (18) et on extrapole à $S_{max} = +\infty$ comme expliqué dans le texte. Tireté bleu : ajustement cubique en la variable $X = (\alpha_c^{4\text{corps}} - \alpha)^{1/2}$ sur l'intervalle $10 \le \alpha \le 13,3$, à savoir $B_{3,1}(0^+) = 2,4401 - 2,1627X + 0,64137X^2 - 0,069596X^3$, où $\alpha_c^{4\text{corps}} \approx 13,3842$ est le seuil efimovien à $3+1$ corps (trait plein vertical). Tireté rouge : idem sur l'intervalle $12,5 \le \alpha \le 13,3$, à savoir $B_{3,1}(0^+) = 2,4637 - 2,2834X + 0,82318X^2 - 0,15269X^3$. La vignette est un agrandissement montrant mieux les annulations avec changement de signe de $B_{3,1}(0^+)$ en $\alpha_1 \approx 3,412$ et $\alpha_2 \approx 6,678$ (pointillés noirs verticaux). (b) Idem pour la demi-dérivée seconde $(1/2)B''_{3,1}(0^+)$ du coefficient d'amas par rapport à la pulsation de piégeage réduite $\tilde{\omega} = \hbar\omega/k_B T$, aux différences près que (i) nous ne donnons pas d'ajustement en la variable $X$, (ii) nous mettons des barres d'erreur car l'incertitude relative peut maintenant dépasser le pour cent, (iii) nous incluons un agrandissement dans l'agrandissement pour bien faire voir que le second point d'annulation (avec changement de signe) de $(1/2)B''_{3,1}(0^+)$ est très proche mais distinct de celui de $B_{3,1}(0^+)$ (toujours indiqué par un pointillé vertical).</figcaption></figure>

<table><caption>**TABLE 2** – Valeurs numériques des quatrièmes coefficients d'amas $B_{3,1}(0^+)$ et $B_{1,3}(0^+)$ du gaz unitaire de fermions dans un piège de raideur infinitésimale, tabulées en fonction du rapport de masse $\alpha = m_\uparrow/m_\downarrow$. Comme on passe d'un coefficient à l'autre en changeant $\alpha$ en $1/\alpha$, nous nous sommes limités à $\alpha \ge 1$. Les incertitudes, non précisées, sont inférieures au pour cent. Ces valeurs sont représentées graphiquement sur la figure 3a.</caption><thead><tr><th>α</th><th>1</th><th>1,5</th><th>2</th><th>2,5</th><th>3</th><th>3,5</th><th>4</th><th>4,5</th><th>5</th><th>5,5</th></tr></thead><tbody><tr><td>B<sub>3,1</sub></td><td>0,02297</td><td>0,01764</td><td>0,01289</td><td>0,008233</td><td>0,003651</td><td>-0,000697</td><td>-0,004571</td><td>-0,007623</td><td>-0,00946</td><td>-0,00964</td></tr><tr><td>B<sub>1,3</sub></td><td>0,02297</td><td>0,02775</td><td>0,03116</td><td>0,03387</td><td>0,03611</td><td>0,03802</td><td>0,03967</td><td>0,04110</td><td>0,04236</td><td>0,04349</td></tr><tr><td>α</td><td>6</td><td>6,5</td><td>7</td><td>7,5</td><td>8</td><td>8,5</td><td>9</td><td>9,5</td><td>10</td><td>10,5</td></tr><tr><td>B<sub>3,1</sub></td><td>-0,00751</td><td>-0,00257</td><td>0,006065</td><td>0,01912</td><td>0,0379</td><td>0,06354</td><td>0,09769</td><td>0,1423</td><td>0,2003</td><td>0,275</td></tr><tr><td>B<sub>1,3</sub></td><td>0,04449</td><td>0,04539</td><td>0,04621</td><td>0,04695</td><td>0,04763</td><td>0,04825</td><td>0,04882</td><td>0,04935</td><td>0,04983</td><td>0,05029</td></tr><tr><td>α</td><td>11</td><td>11,5</td><td>12</td><td>12,5</td><td>12,75</td><td>13</td><td>13,1</td><td>13,2</td><td>13,3</td><td>20</td></tr><tr><td>B<sub>3,1</sub></td><td>0,3721</td><td>0,4996</td><td>0,6714</td><td>0,9176</td><td>1,090</td><td>1,329</td><td>1,457</td><td>1,623</td><td>1,867</td><td></td></tr><tr><td>B<sub>1,3</sub></td><td>0,05071</td><td>0,05110</td><td>0,05147</td><td>0,05182</td><td>0,05198</td><td>0,05214</td><td>0,05220</td><td>0,05227</td><td>0,05233</td><td>0,05523</td></tr></tbody></table>
## 2.5. Près du seuil efimovien à 3 + 1 corps

La limite $\alpha \to (\alpha_c^{4\text{corps}})^-$ dans l'équation (30) est difficile à atteindre numériquement avec pré-
cision pour plusieurs raisons. Premièrement, si $\alpha \to +\infty$, le taux de décroissance exponentielle
$\ln(1/|z_0|)$ dans l'équation (20) tend vers zéro comme $\alpha^{-1/2}$ donc la convergence de la série $I_{3,1}(0^+)$
devient de plus en plus lente avec $\ell$; heureusement, comme le rapport de masse $\alpha$ reste infé-
rieur au seuil efimovien à quatre corps, $\ln(1/|z_0|)$ ne devient pas si petit que ça (il reste supérieur
à 0,267) et le problème est résolu par l'accélération de convergence (18) si l'on se contente d'une
erreur sur $B_{3,1}(0^+)$ inférieure au pour cent (voir la figure 2). Deuxièmement, l'effet Efimov à 3 + 1
corps sur le point d'apparaître dans la voie $\ell = 1, \varepsilon = +1$ oblige à augmenter $x_{\text{max}}$ aux valeurs
considérées dans la référence [39], du moins dans les voies de moment cinétique $\ell \le 5$, ce qui
augmente beaucoup le temps de calcul et conduit à une singularité de type

$$ (\alpha_c^{4\text{corps}} - \alpha)^{1/2} $$

dont il faut tenir compte dans l'extrapolation de $B_{3,1}(0^+)$ à $\alpha = \alpha_c^{4\text{corps}}$, voir les tiretés sur la
figure 3a et la légende de la figure<fnref n="13" />. Troisièmement, le seuil de l'effet Efimov à trois corps
$\alpha_c^{3,1} \approx 13,60697$ est proche et affecte toutes les voies de moment cinétique $\ell$ du problème à
quatre corps [40], même s'il se produit dans la voie de moment cinétique $L = 1$ du problème à
$2+1$ fermions. En effet, le spectre continu de l'opérateur $M_{3,1}^{(\ell,\varepsilon)}(iS)$ est la réunion de continuums
correspondant à $2+1$ fermions fortement corrélés par les interactions et un fermion ↑ découplé,
donc paramétrés par un moment cinétique $L$ à trois corps [39]. La contribution du continuum de
moment cinétique $L$ au coefficient d'amas $3+1$ dans la voie $(\ell, \varepsilon)$ s'écrit à un facteur près [40] :

$$ I_{3,1}^{(\ell, \varepsilon)}(0^+) |_{C^{0, L}} \propto \int_{\mathbb{R}} dS \int_{\mathbb{R}^+} dk S \partial_S \theta_L(k, S) \frac{d}{dk} \ln \Lambda_{2,1}^{(L)}(ik) \quad (31) $$

Les modes propres du continuum ont asymptotiquement une structure d'onde plane dans l'es-
pace de la variable $x$, c'est-à-dire qu'ils sont superposition, lorsque $x \to +\infty$, d'une onde inci-
dente $\exp(-ikx)$ et d'une onde réfléchie $-\exp[i\theta_L(k, S)]\exp(ikx)$ où le déphasage $\theta_L(k, S)$ est fonc-
tion du nombre d'onde $k > 0$ (fictif car $x$ n'est pas une vraie position), de l'exposant d'échelle ima-
ginaire pur iS et du moment cinétique à trois corps $L$. On voit apparaître dans l'expression (31) la
fonction transcendante d'Efimov $\Lambda_{2,1}^{(L)}$ dont les exposants d'échelle du problème unitaire à $2+1$
corps de moment cinétique $L$ sont les racines. Ceci est dû au fait que les modes du continuum
de nombre d'onde $k$ sont associés à la valeur propre $\Lambda_{2,1}^{(L)}(ik)$ de $M_{3,1}^{(\ell, \varepsilon)}(iS)$ [39]<fnref n="14" />. Le point clé est
maintenant que le bord inférieur du continuum de $L=1$, à savoir le minimum de $\Lambda_{2,1}^{L=1}(ik)$ par
rapport à la variable $k$, tend vers zéro lorsque $\alpha \to \alpha_c^{2,1}$ à cause de l'effet Efimov à $2+1$ corps. Ceci

$$ \frac{dB_{3,1}(0^+)}{dX} = -(\ell_0 + 1/2) c_0^{1/2} \approx -2,2 \quad \text{en} \quad X \equiv (\alpha_c^{3,1} - \alpha)^{1/2} = 0. $$

Cette contrainte est assez bien satisfaite par les ajustements de la figure 3a.

<footnote n="13">¹³Cette singularité est présente dans l'un des exposants d'échelle $s_i$ du problème unitaire à 3 + 1 corps, celui $s_0$ dont le carré s'annule en changeant de signe au seuil efimovien, dans la voie $\ell_0 = 1$ [39]. Or $B_{3,1}(0^+)$ dépend linéairement des exposants d'échelle. On le voit bien sur la composante universelle du troisième coefficient d'amas du gaz unitaire de bosons dans l'équation (36) de la référence [20] en faisant tendre $\bar{\omega}$ (noté $x$ dans cette référence) vers zéro; on le voit aussi pour $2+1$ fermions dans la référence [21]. Plus précisément, $B_{3,1}(0^+) = -(\ell_0 + 1/2)s_0 + \dots$ où l'ellipse est une fonction lisse de $\alpha$ au voisinage de $\alpha_c^{3,1}$. Comme $s_0^2 \approx c_0(\alpha_c^{3,1} - \alpha)$ avec $c_0 \approx 2,2$ près du seuil [39], on trouve que</footnote>

<footnote n="14">¹⁴ En l'absence de coupure dans l'espace des $x$ ($x_{\text{max}} = +\infty$), on tombe sur un paradoxe : $k$ décrit l'ensemble continu $\mathbb{R}^+$ indépendamment de la variable $S$ et la dérivée de la valeur propre $\Lambda_{2,1}^{(L)}(ik)$ par rapport à $S$ est nulle, donc le continuum ne devrait pas contribuer à $I_{3,1}$ dans l'équation (4). La bonne façon de raisonner est de mettre une coupure $x_{\text{max}}$ que l'on fait tendre vers l'infini à la fin, avec la condition que les modes propres s'annulent en $x = x_{\text{max}}$. L'équation $\theta_L(k, S) = -2kx_{\text{max}}$ [$2\pi$] qui en résulte quantifie $k$, c'est-à-dire le restreint à un ensemble discret, et le fait dépendre de $S$ si le déphasage $\theta_L(k, S)$ en dépend [40].</footnote>
a deux conséquences. D'abord, une conséquence pratique dans le calcul de $I_{3,1}$ sur ordinateur : comme le continuum est rendu discret par la troncature de $x$ à $x_{max}$, son bord inférieur s'écarte du vrai bord d'environ $1/x_{max}^{2}$<fnref n="15" />; cet artefact numérique oblige à augmenter significativement $x_{max}$, heureusement seulement pour les voies de moment cinétique à quatre corps $\ell \le 5$ tant que $\alpha \le 13,3$. Ensuite, une conséquence physique : une seconde singularité apparaît dans $I_{3,1}(0^+)$, de la forme $(\alpha_c^{2,1} - \alpha)^{1/2}$, donc elle aussi en racine carrée mais centrée sur le seuil efimovien à trois corps<fnref n="16" />. Cette singularité à trois corps se produit en un point très proche de la singularité à quatre corps, ce qui jette un doute sur la précision de l'extrapolation faite dans la figure 3a qui n'en teneait pas compte. Nous y remédions par des fonctions d'ajustement incluant les deux singularités, polynomiales en les quantités $X = (\alpha_c^{3,1} - \alpha)^{1/2}$ et $Y = (\alpha_c^{2,1} - \alpha)^{1/2}$ :

$$ 
\begin{gathered}
I_{3,1}(0^+) = A_0 + A_1 X + A_2 (Y - Y_0) + A_3 X^2 \\
\text{et } I_{3,1}(0^+) = A_0 + A_1 X + A_2 (Y - Y_0) + A_3 X^2 + A_4 X^3 + A_5 (Y^3 - Y_0^3)
\end{gathered}
\qquad (32) $$

où $Y_0 = (Y^2 - X^2)^{1/2} = (\alpha_c^{2,1} - \alpha_c^{3,1})^{1/2}$ est aussi la valeur de $Y$ en $X = 0$. Sur l'intervalle $10 \le \alpha \le 13,3$, ceci conduit à $A_0 = 2,495$ et $A_0 = 2,447$ d'où le résultat final (30), qui conforte celui plus naïf de la figure 3a.

## 2.6. $\operatorname{Cas} \bar{\omega} \neq 0$

Pour terminer cette section, étudions brièvement la dépendance du coefficient d'amas $B_{3,1}(\bar{\omega})$ en la raideur du piège. Une façon d'en rendre compte qui soit utile aux expériences sur les gaz d'atomes froids est de calculer le premier écart à la limite de raideur nulle, de degré pair en $\bar{\omega}$ puisque l'intégrande de l'équation (4) est une fonction paire de $\bar{\omega}$ :

$$ B_{3,1}(\bar{\omega}) = \lim_{\bar{\omega} \to 0^+} B_{3,1}(0^+) + \frac{1}{2} B''_{3,1}(0^+) \bar{\omega}^2 + O(\bar{\omega}^4) \qquad (33) $$

Ceci permet en effet de quantifier l'erreur due à l'approximation d'homogénéité locale, systématiquement utilisée dans l'expérience, et qui revient à garder seulement le premier terme du développement de Taylor (33). Il est d'ailleurs facile de voir que le petit paramètre contrôlant cette approximation ne peut être que $\bar{\omega}$ à la limite unitaire : l'approximation n'a de sens que si la longueur de corrélation à l'équilibre du gaz homogène dans chaque état de spin $\sigma$, ici nécessairement la longueur thermique de de Broglie $\lambda_{\sigma}$ par invariance d'échelle dans le cas non dégénéré, est beaucoup plus petite que le rayon spatial du gaz piégé $R_{\sigma} = (k_B T / m_{\sigma} \omega^2)^{1/2}$; on a bien $\lambda_{\sigma} / R_{\sigma} \propto \hbar \omega / k_B T = \bar{\omega} \ll 1$. Pour obtenir une expression intégrale de la dérivée seconde $B''_{3,1}(0^+)$ et mettre en œuvre une accélération de convergence, il suffit de dériver deux fois par rapport à $\bar{\omega}$ l'équation (4) sous le signe intégral, l'équation (18) sous le signe somme et les expressions (47), (59)) puis de faire tendre $\bar{\omega}$ vers zéro.

$$ \propto \int_{0}^{+\infty} dk \left[ \frac{k^2}{(k^2 + \sigma_0^2)} - \frac{k^2}{(k^2 + 1)} \right] \propto 1 - \sigma_0 \quad \text{où} \quad \sigma_0 \sim \gamma_0 (\alpha_c^{2,1} - \alpha)^{1/2} \quad \text{et} \quad \gamma_0 \approx 0,438. $$

<footnote n="15">¹⁵En effet, le nombre d'onde $k$ varie par pas de l'ordre de $1/x_{max}$, comme le montre la note 14, et $\Lambda_{2,1}^{L=1}(ik)$ varie quadratiquement près de son minimum.</footnote>

<footnote n="16">¹⁶Å gauche de ce seuil, on écrit à la Weierstrass $\ln \Lambda_{2,1}^{(L=1)}(ik) = \ln[(k^2 + \sigma_0^2)/(k^2 + 1)] + \dots$ où l'ellipse est une fonction régulière de $\alpha$ même au seuil et $\sigma_0 > 0$ est l'exposant d'échelle du problème à 2 + 1 corps dont le carré s'annule en changeant de signe au seuil. Tout près du seuil, $\sigma_0 \ll 1$ et, en admettant comme dans la référence [40] que $\theta_{L=1}(k, S) \sim b(S)k$ lorsque $k \to 0$, où le rayon de diffusion $b(S)$ est une fonction inconnue de $S$, on trouve que $I_{3,1}^{(\ell, \epsilon) \neq (1,+)}(0^+) \text{ contient une singularité}$</footnote>

<footnote>Nous avons confirmé cette prédiction par un calcul numérique spécifique de $I_{3,1}(0^+)$ dans la voie $\ell = 0$ tout près du seuil, allant jusqu'à $\alpha = 13,59$; une belle loi linéaire en $Y$ est observée et un ajustement donne $I_{3,1}^{(\ell=0)}(0^+) \approx 0,026 + 0,040 Y$ où la variable $Y$ est celle de l'équation (32).</footnote>
Le résultat est représenté en fonction du rapport de masse sur la figure 3b<fnref n="17" />. Cherchons à l'interpréter dans un scénario très simple : comme dans le cas des masses égales $m_↑ = m_↓$ étudié dans la référence [41], $B_{3,1}(\bar{\omega})$ serait une fonction monotone de $\bar{\omega}$, bien entendu de limite nulle à l'infini. La dérivée seconde en l'origine serait alors toujours de signe opposé à la valeur en $\bar{\omega} = 0$. Ceci explique la figure 3b si on la regarde à grande échelle : $B''_{3,1}(0^+)$ semble effectivement s'annuler en changeant de signe aux mêmes rapports de masse $\alpha = \alpha_1$ et $\alpha = \alpha_2$ que $B_{3,1}(0^+)$. On voit cependant sur un premier agrandissement, en vignette dans la figure 3b, que ce scénario est pris en défaut en $\alpha_1$ (il semble toujours tenir la route en $\alpha_2$). En effet, comme le montre la figure 4a, lorsque $\alpha$ se rapproche de $\alpha_1 \approx 3,412$ par valeurs inférieures (courbes de haut en bas), la fonction $\bar{\omega} \mapsto B_{3,1}(\bar{\omega})$ cesse d'être monotone, devient convexe près de l'origine et atteint un maximum absolu (positif) en un point $\bar{\omega}_0 > 0$ parti de zéro avant que $B_{3,1}(0^+)$ ne devienne < 0. Ce maximum absolu se maintient quand $\alpha$ continue de croître au-delà de $\alpha_1$, mais sa position se déplace vers $+\infty$ et sort de la figure 4a. Un second agrandissement, dans la vignette de la figure 3b, montre que le scénario simple échoue également (mais de peu) au second point d'annulation $\alpha_2 \approx 6,678$ : lorsque $\alpha$ se rapproche de $\alpha_2$ par valeurs supérieures (courbes de haut en bas sur la figure 4b<fnref n="18" />), un minimum absolu (négatif) venu du côté $\bar{\omega} = +\infty$ se rapproche de l'origine et devient plus prononcé, avant que $B_{3,1}(0^+)$ ne devienne à son tour < 0. Par commodité, nous donnons les valeurs de $B_{3,1}(\bar{\omega})$ et de sa demi-dérivée seconde en $\bar{\omega} = 0^+$ sous forme numérique dans la table 3, pour les rapports de masse de la figure 4. Soulignons pour terminer que le scénario simple de départ était en réalité hautement improbable car il impliquait que le coefficient d'amas à $3+1$ corps dans le piège était nul pour toute raideur, $B_{3,1}(\bar{\omega}) \equiv 0$, aux rapports de masse où $B_{3,1}(0^+) = 0$.

<table><caption>TABLE 3 – Valeurs numériques du quatrième coefficient d'amas $B_{3,1}(0^+)$ et de sa demi-dérivée seconde $(1/2)B''_{3,1}(0^+)$ par rapport à $\bar{\omega}$, données avec leurs chiffres significatifs, pour les rapports de masse de la figure 4. Les demi-dérivées sont représentées graphiquement (pour certaines avec une barre d'erreur) sur la figure 3b.</caption><thead><tr><th>α</th><th>2,5</th><th>3</th><th>3,25</th><th>3,375</th><th>3,5</th><th>4</th><th>4,5</th><th>5</th><th>5,5</th></tr></thead><tbody><tr><td>B<sub>3,1</sub></td><td>0,00823</td><td>0,00365</td><td>0,00143</td><td>0,00035</td><td>-0,00070</td><td>-0,00457</td><td>-0,00763</td><td>-0,00947</td><td>-0,00963</td></tr><tr><td>1/2 B''<sub>3,1</sub></td><td>-0,00660</td><td>0,00055</td><td>0,00420</td><td>0,00596</td><td>0,00766</td><td>0,0140</td><td>0,0190</td><td>0,0217</td><td>0,0212</td></tr><tr><td>α</td><td>6,5</td><td>6,5625</td><td>6,625</td><td>6,65625</td><td>6,6875</td><td>6,71875</td><td>6,75</td><td>6,8125</td><td>6,875</td></tr><tr><td>B<sub>3,1</sub></td><td>-0,00260</td><td>-0,00174</td><td>-0,00083</td><td>-0,00035</td><td>0,00015</td><td>0,00066</td><td>0,00119</td><td>0,00229</td><td>0,00346</td></tr><tr><td>1/2 B''<sub>3,1</sub></td><td>0,0047</td><td>0,0027</td><td>0,0007</td><td>-0,0003</td><td>-0,0013</td><td>-0,0026</td><td>-0,0037</td><td>-0,0063</td><td>-0,0088</td></tr></tbody></table>

### 3. Coefficient d'amas à $2+2$ fermions dans un piège

Nous expliquons dans cette section comment effectuer un calcul numérique efficace de la quantité $I_{2,2}(\bar{\omega})$ définie par l'équation (4), en nous limitant pour simplifier au cas $\bar{\omega} = 0^+$, c'est-à-dire à un piège de raideur infinitésimale. Le coefficient d'amas $B_{2,2}(0^+)$ du gaz unitaire s'en

$$[I_{3,1}^{(\ell, \epsilon)}]''(0^+) - [J_{3,1}^{(\ell, \epsilon)}]''(0^+) \text{ pour } \ell > \ell_{\text{max}} \text{ par } A\ell^{\gamma} \text{Re}(e^{i\psi} z_0^{\ell})$$

<footnote n="17">¹⁷Pour les valeurs de $\alpha$ proches de $\alpha_2$, nous améliorons la méthode d'accélération de convergence en approximant</footnote>

<footnote>plutôt que par 0 comme dans l'équation (18); les paramètres réels $A, \psi$ et $\gamma$ sont obtenus par ajustement sur l'intervalle $6 \le \ell \le \ell_{\text{max}}$, et le nombre complexe $z_0$ est celui de la loi asymptotique (20).</footnote>

<footnote n="18">¹⁸En faisant varier tous les paramètres de troncature et de discrétisation dans le calcul numérique, nous estimons à $2 \times 10^{-5}$ l'incertitude sur les courbes de la figure 4b les plus proches de l'axe horizontal (par exemple pour $\alpha = 6,6875$); les dépendances en $\bar{\omega}$ représentées sont donc significatives.</footnote>
<figure><img src="image_6.png" /><figcaption>FIGURE 4 – Quatrième coefficient d'amas $B_{3,1}(\tilde{\omega})$ du gaz unitaire de fermions piégé en fonction de la raideur réduite $\tilde{\omega}^2 = (\hbar\omega/k_B T)^2$ des potentiels harmoniques, pour différents rapports de masse $\alpha = m_1/m_\downarrow$ proches des points d'annulation $\alpha_1 \approx 3,412$ (a) et $\alpha_2 \approx 6,678$ (b) de $B_{3,1}(0^+)$, ce qui est utile à la compréhension des agrandissements sur la figure 3b, c'est-à-dire du rapport entre les signes de la valeur et de la courbure de $B_{3,1}(\tilde{\omega})$ en l'origine. Dans (b), les pentes des tangentes à l'origine (en pointillé) sont tirées de la figure 3b.</figcaption></figure>

déduit par la conjecture (5). Comme les deux états de spin $\uparrow$ et $\downarrow$ jouent ici des rôles parfaitement symétriques, $B_{2,2}(0^+)$ est invariant par changement du rapport de masse $\alpha$ en son inverse $1/\alpha$ et l'on peut limiter les calculs numériques au cas $\alpha \ge 1$. Au contraire de la section 2 précédente, on peut aller jusqu'au seuil efimovien à trois corps, $\alpha < \alpha_c^{2,1}$, puisqu'il n'existe pas d'effet Efimov à 2+2 corps [40].

### 3.1. Formulation du problème

Rappelons l'expression de l'opérateur hermitien $M_{2,2}^{(\ell, \epsilon)}(iS)$ donnée dans la référence [40] pour le moment cinétique $\ell$, la parité $\epsilon = \pm$ (limitée à $+\text{si } \ell = 0$) et l'exposant d'échelle imaginaire pur iS, avec la même notation de Schrödinger-Dirac que dans l'équation (14) mais la différence que la variable $x$ varie sur la droite réelle tout entière :

$$ \langle x, u | \langle \ell, m_z | M_{2,2}^{(\ell, \epsilon)} (iS) | f \rangle = \mathcal{D}_{2,2}(x, u) f_{m_z}(x, u) + \int_{-\infty}^{+\infty} dx' \int_{-1}^{1} du' \sum_{\substack{m'_z | (-1)^{m'_z} = \epsilon \\ m'_z \neq 0}} K_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) f_{m'_z}(x', u') \quad (34) $$

avec une partie diagonale indépendante de iS, du moment cinétique et du nombre quantique magnétique $m_z$ :

$$ \mathcal{D}_{2,2}(x, u) = \left[ \frac{\alpha}{(1+\alpha)^2} \left(1 + \frac{u}{\text{ch}x}\right) + \frac{e^{-x} + \alpha e^x}{2(\alpha+1)\text{ch}x} \right]^{1/2} \quad (35) $$

et un noyau intégral matriciel décomposé en trois contributions $K_{2,2} = K_1 + K_2 + K_3$ écrites ligne par ligne dans cet ordre :
$$
\begin{align}
K_{2,2}^{(\ell)}(x, u, m_z; x', u', m_z') & \\
&= \left(\frac{e^x \operatorname{ch} x'}{e^{x'} \operatorname{ch} x}\right)^{iS/2} \left(\frac{e^{x+x'}}{4 \operatorname{ch} x \operatorname{ch} x'}\right)^{1/4} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{e^{-i m_z \theta/2} \langle \ell, m_z | e^{i \phi L_x / \hbar} | \ell, m_z' \rangle e^{i m_z' \theta'/2}}{\operatorname{ch}(x-x') + \frac{1}{1+\alpha} [(u+e^{-x})(u'+e^{-x'}) + v v' \cos \phi]} \nonumber \\
&\quad + \left(\frac{e^{-x} \operatorname{ch} x'}{e^{-x'} \operatorname{ch} x}\right)^{iS/2} \left(\frac{e^{-x-x'}}{4 \operatorname{ch} x \operatorname{ch} x'}\right)^{1/4} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{e^{i m_z \theta/2} \langle \ell, m_z | e^{i \phi L_x / \hbar} | \ell, m_z' \rangle e^{-i m_z' \theta'/2}}{\operatorname{ch}(x-x') + \frac{\alpha}{1+\alpha} [(u+e^x)(u'+e^{x'}) + v v' \cos \phi]} \nonumber \\
&\quad + \frac{(-1)^{\ell+1}}{4\pi [(u+\operatorname{ch} x)(u'+\operatorname{ch} x') \operatorname{ch} x \operatorname{ch} x']^{1/4}} \left(\frac{\operatorname{ch} x' \operatorname{ch} x}{\operatorname{ch} x \operatorname{ch} x'}\right)^{iS/2} \nonumber \\
&\quad - \frac{e^{i m_z \gamma(x, u)} \langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m_z' \rangle e^{-i m_z' \gamma(x', u')}}{\left(\frac{e^{-x'} + \alpha e^{x'}}{1+\alpha}\right) (u+\operatorname{ch} x) + \left(\frac{e^{-x} + \alpha e^x}{1+\alpha}\right) (u'+\operatorname{ch} x')} \tag{36}
\end{align}
$$

avec l'angle $\gamma(x, u) = \arctan\{\operatorname{th}(x/2)[(1-u)/(1+u)]^{1/2}\}$, l'angle $\theta \in [0, \pi]$ tel que $u = \cos\theta$ et la notation $v = \sin\theta$. La troisième contribution $K_3$ est non nulle seulement dans le secteur de parité $\epsilon = (-1)^\ell$ (dans l'autre secteur, on a $\langle \ell, m_x = 0 | \ell, m_z \rangle = 0$). Contrairement aux deux premières, c'est une fonction singulière des coordonnées $(x, u)$ et $(x', u')$, divergente au point $(x=0, u=-1)$, ou si l'on préfère $(x=0, \theta=\pi)$. Dans le calcul numérique, limité au rectangle $(x, \theta) \in [-x_{\max}, x_{\max}] \times [0, \pi]$, on isole la singularité par un demi-disque de rayon $R$ et de centre $(0, \pi)$, dans lequel on utilise une grille logarithmico-polaire comme dans la référence [40]<fnref n="19" />, et au-dehors duquel on utilise le même type de grille que dans la section 2.

## 3.2. Approximant asymptotique

Il reste à mettre en œuvre la même technique d'accélération de convergence que dans l'équation (18) dans le cas particulier $\bar{\omega} = 0^+$,

$$ I_{2,2}(0^+) \approx \sum_{\ell=0}^{\ell_{\max}} \sum_{\epsilon} \left[ I_{2,2}^{(\ell, \epsilon)}(0^+) - J_{2,2}^{(\ell, \epsilon)}(0^+) \right] + \sum_{\ell=0}^{+\infty} \sum_{\epsilon} J_{2,2}^{(\ell, \epsilon)}(0^+) \quad (37) $$

où, comme dans la section 2, $I_{2,2}^{(\ell, \epsilon)}(0^+)$ est la contribution de la voie de moment cinétique $\ell$ et de parité $\epsilon$ à la quantité cherchée $I_{2,2}(0^+)$, et $J_{2,2}^{(\ell, \epsilon)}(0^+)$ en est un approximant asymptotique, valable à grand moment cinétique. La référence [40] donne un tel approximant seulement dans le secteur de parité gentil $\epsilon = (-1)^{l-1}$ (l'opposé de la parité naturelle) où $K_3 \equiv 0$, en traitant perturbativement $K_1$ et $K_2$ au second ordre dans le logarithme du déterminant de $M_{2,2}^{(\ell, \epsilon)}(\mathrm{iS})$. Le calcul ressemble beaucoup à l'équation (19), avec la simplification que les termes linéaires et les termes carrés sont indépendants de $S$ et peuvent être absorbés dans la constante, si bien qu'il ne reste que les termes rectangles [40] :

$$ \ln \det \left[ M_{2,2}^{(\ell, \epsilon)} (\mathrm{iS}) \right] \stackrel{\epsilon=(-1)^{\ell-1}}{=} \mathrm{cte} - \mathrm{Tr}^{\epsilon} \left[ \mathcal{D}_{2,2}^{-1} K_1 \mathcal{D}_{2,2}^{-1} K_2 \right] + \dots \quad (38) $$

Dans le secteur de parité méchant $\epsilon = (-1)^\ell$ (celui de la parité naturelle), une difficulté surgit : à cause de son caractère divergent, on ne peut pas traiter $K_3$ perturbativement. On effectue

<footnote n="19">¹⁹ On pose $(x, \pi - \theta) = \mathrm{Re}^\ell(\cos \psi, \sin \psi)$, où $t_{\min} < t < 0$ et $0 < \psi < \pi$. On a choisi $R = 2/5$. L'intégration sur $\psi$ est faite par la méthode de Gauss-Legendre à 15 points (série $0 \le \ell \le 4$) ou 25 points (série $0 \le \ell \le 6$), celle sur $t$ avec la méthode du point milieu et un pas $dt = 0,1$ ou $dt = 0,2$. On extrapole quadratiquement à $t_{\min} = -\infty$ à partir des trois choix $t_{\min} = -10$, $t_{\min} = -20$ et $t_{\min} = -40$.</footnote>
d'abord un changement de jauge sur les $K_j$ éliminant la dépendance en $S$ de $K_3$ sans changer
le déterminant, ce que l'on signale par un tilde,

$$
\tilde{K}_j (x, u, m_z; x', u', m'_z) \equiv \frac{e^{-im_z \gamma(x, u)}}{[(u + \operatorname{ch} x) \operatorname{ch} x']^{-iS/2}} K_j (x, u, m_z; x', u', m'_z) \frac{e^{im'_z \gamma(x', u')}}{[(u' + \operatorname{ch} x') \operatorname{ch} x']^{iS/2}} \quad (39)
$$

puis on isole une partie purement « externe » de $\tilde{K}_3$, agissant sur l'espace « orbital » des $(x, u)$
mais pas sur celui des $|\ell, m_z\rangle$, donc telle que $\tilde{K}_3 = \tilde{K}_3^{\text{ext}} \otimes |\ell, m_x = 0\rangle\langle\ell, m_x = 0|$. Enfin, on effectue
un développement du logarithme du déterminant en puissances de $K_1$ et $K_2$ sans hypothèse
sur $\tilde{K}_3^{\text{ext}}$ :

$$
\begin{align}
\ln \det \left[ M_{2,2}^{(\ell, \epsilon)} (\mathrm{iS}) \right]^{\epsilon = (-1)^\ell} \mathrm{cte} + \ln \det \left[ 1 + \frac{1}{\mathcal{D}_{2,2}^{-1} + \tilde{K}_3^{\mathrm{ext}} \otimes |\ell, m_x = 0\rangle \langle \ell, m_x = 0|} (\tilde{K}_1 + \tilde{K}_2) \right] \nonumber \\
&= \mathrm{cte} - \mathrm{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\mathrm{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle] - \mathrm{Tr}^\epsilon [\mathcal{D}_{2,2}^{-1} \tilde{K}_1 \mathcal{D}_{2,2}^{-1} \tilde{K}_2] \nonumber \\
&\quad + \mathrm{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\mathrm{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) \mathcal{D}_{2,2}^{-1} (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle] \nonumber \\
&\quad - \frac{1}{2} \mathrm{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\mathrm{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle] \nonumber \\
&\quad + \dots \tag{40}
\end{align}
$$

Dans cette expression apparaît une forme resommée de la partie externe :

$$
\tilde{K}_3^{\text{ext}} \equiv \mathcal{D}_{22} [\mathcal{D}_{22}^{-1} - (\mathcal{D}_{22} + \tilde{K}_3^{\text{ext}})^{-1}] \mathcal{D}_{22} = \tilde{K}_3^{\text{ext}} - \tilde{K}_3^{\text{ext}} (\mathcal{D}_{22} + \tilde{K}_3^{\text{ext}})^{-1} \tilde{K}_3^{\text{ext}} \quad (41)
$$

Une écriture optimale de $\tilde{K}_3^{\text{ext}}$ est obtenue par reparamétrage de l'ansatz de Faddeev du problème
à 2 + 2 corps : on considère que les fonctions sur lesquelles agit l'opérateur $M_{2,2}$ dépendent
désormais des vecteurs d'onde relatif et du centre de masse des particules 2 et 4 comme dans
la référence [52] au lieu des vecteurs d'onde individuels $\mathbf{k}_2$ et $\mathbf{k}_4$ comme dans les références [39,
40], ce qui permet d'éviter le demi-disque autour de la singularité et nous conduit à inverser
numériquement un opérateur agissant sur une seule variable réelle au lieu des deux variables
$(x, u)$. Ceci amène un gain de temps de calcul et une simplification considérables. Le lecteur est
renvoyé à l'annexe C pour plus de détails. On trouve numériquement que le second terme au
troisième membre de l'équation (40), formellement du premier ordre, est en réalité du même
ordre de grandeur que le troisième terme (leurs contributions à $J_{2,2}^{(\ell,\epsilon)}(0^+)$ tendent vers zéro
exponentiellement en $\ell$ avec le même taux, voir la figure 5a). Le noyau resommé $\tilde{K}_3^{\text{ext}}$ est donc,
comme $K_1$ et $K_2$, un infiniment petit du premier ordre; par conséquent, nous négligeons les
quatrième et cinquième termes de l'équation (40) pour garder

$$
\ln \det \left[ M^{(\ell, \epsilon)} (\mathrm{iS}) \right]^{\epsilon = (-1)^\ell} \mathrm{cte} - \mathrm{Tr}_{x,u} \left[ \mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\mathrm{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle \right] - \mathrm{Tr}^\epsilon \left[ \mathcal{D}_{2,2}^{-1} \tilde{K}_1 \mathcal{D}_{2,2}^{-1} \tilde{K}_2 \right] + \dots \quad (42)
$$

Si l'on veut, on peut défaire le changement de jauge (39) dans le troisième terme sans changer
la trace, pour retrouver la contribution (38). Dans l'annexe B, nous donnons une expression plus
explicite de l'approximant asymptotique $J_{2,2}^{(\ell,\epsilon)}(0^+)$ déduit des développements (38) et (42), voir
les équations ((71), (72), (78), (80)), ainsi que de sa somme sur $\ell$ et sur $\epsilon$, voir les équations ((81),
(86)).

3.3. *Application et résultats*

La somme sur $\ell$ dans $I_{2,2}(0^+)$ présente le même phénomène de convergence lente que dans $I_{3,1}(0^+)$ aux grandes valeurs du rapport de masse $\alpha$. Nous trouvons à nouveau que l'approximant asymptotique, ici $J_{2,2}^{(\ell,\epsilon)}(0^+)$, tend vers zéro exponentiellement en $\ell$, avec un taux $c$ que nous
<figure><img src="image_7.png" /><figcaption>FIGURE 5 – Quelques résultats numériques sur l’approximant asymptotique $J_{2,2}^{(\ell,\epsilon)}(0^+)$ de la contribution $I_{2,2}^{(\ell,\epsilon)}(0^+)$ de moment cinétique $\ell$ et de parité $\epsilon$ à la quantité $I_{2,2}(0^+)$ de l’équation (4). (a) Pour un rapport de masse $\alpha = 5$ et le secteur de parité méchant $\epsilon = (-1)^\ell$ (celui de la parité naturelle), valeur absolue des trois termes de $J_{2,2}^{(\ell,\epsilon)}(0^+)$ dans l’expression (71) (symboles ◦, + et × dans cet ordre) en fonction du moment cinétique $\ell$. (b) Taux $c$ de décroissance exponentielle en $\ell$ de $J_{2,2}^{(\ell,\epsilon)}(0^+)$ comme dans l’équation (43), en fonction de la racine carrée du rapport de masse inverse. Les barres d’erreur donnent l’écart entre les valeurs de $c$ tirées de diverses fonctions d’ajustement (pour $\epsilon = (-1)^{\ell-1}$, $\ln|J_{2,2}^{(\ell,\epsilon)}(0^+)|$ est pris comme fonction affine de $\ell$ plus un terme en $1/\ell$, ou plus un terme en $\ln\ell$, ou plus les deux). La droite tiretée est une extrapolation affine à $\alpha^{-1/2} = 0$. (c) Pour $\alpha = 13,3$ et $\epsilon = (-1)^\ell$, valeur exacte $I_{2,2}^{(\ell,\epsilon)}(0^+)$ obtenue numériquement (carrés) et approximant asymptotique $J_{2,2}^{(\ell,\epsilon)}(0^+)$ (étoiles) en fonction de $\ell$.</figcaption></figure>

calculons numériquement par commodité dans le secteur de parité gentil $\epsilon = (-1)^{\ell-1}$ et pour un rapport de masse $\alpha \gtrsim 3$:<fnref n="20" />

$$J_{2,2}^{(\ell,\epsilon)}(0^+) \stackrel{\epsilon=(-1)^{\ell-1}}{\underset{\ell \to +\infty}{\equiv}} \epsilon \exp[-c\ell + O(\ln \ell)] \quad (43)$$

Le taux $c$ est représenté en fonction du rapport de masse sur la figure 5b. On y observe, comme dans la section 2, une loi $c \propto \alpha^{-1/2}$ aux grands rapports de masse. Dans ce régime, la méthode d'accélération de convergence est un allié indispensable du calcul numérique, comme le montre la figure 5c. Elle nous permet d'obtenir le quatrième coefficient d'amas $B_{2,2}(0^+)$ du système piégé, donné sous forme numérique dans la table 4 et représenté en fonction du rapport de masse sur la figure 6, avec une incertitude inférieure au pour cent en un temps de calcul raisonnable. À cause de l'effet Efimov à 2 + 1 corps, ce coefficient présente une singularité en $(\alpha_c^{2,1} - \alpha)^{1/2}$ près du seuil, pour la même raison que celle exposée autour de l'équation (31) dans la section 2. Nous en tenons compte dans les ajustements en tireté sur la figure 6, pour obtenir l'extrapolation

$$B_{2,2}(0^+) \underset{\alpha \to (\alpha_c^{3\text{corps}})^{-}}{\stackrel{\alpha \to -\infty}{\longrightarrow}} -0,737 \pm 0,007 \quad (44)$$

<footnote n="20">²⁰Les comportements dominants dans les deux secteurs de parité diffèrent seulement par leurs lois de puissance en $\ell$. Quand le rapport de masse est trop proche de 1, $J_{2,2}^{(\ell,\epsilon=(-1)^{\ell-1})}(0^+)$ tend vers zéro en oscillant ce qui rend l'extraktion numérique du taux plus difficile.</footnote>
<figure><img src="image_8.png" /></figure>

**FIGURE 6** – Quatrième coefficient d'amas $B_{2,2}(0^+)$ du gaz unitaire de fermions dans un piège harmonique de raideur infinitésimale, en fonction du rapport de masse $\alpha = m_1/m_1$ restreint à $\alpha \ge 1$ par symétrie $\alpha \leftrightarrow 1/\alpha$. Disques (reliés par un trait fin) : notre calcul numérique de l'expression (4) pour $(n_1, n_1) = (2,2)$ et application de la conjecture (5), avec les paramètres $x_{max} = -x_{min} = 20$ et $n_\theta = 15$ dans la série $0 \le l \le 6$, $x_{max} = -x_{min} = 5$ et $n_\theta = 25$ dans la série $0 \le l \le 6$ (le recouvrement des deux séries sert de vérification), et dans tous les cas $dx = 1/15$, $S_{max} = 12$, $dS = 12/100$ (voir la note 19 pour le traitement de la singularité du noyau intégral en $(x, \theta) = (0, \pi)$). On utilise la méthode d'accélération de convergence (37) et la même extrapolation à $S_{max} = +\infty$ que dans la section 2. Trait vertical fin : seuil efimovien à trois corps $\alpha = \alpha_c^{2,1}$. Tireté bleu : ajustement cubique en la variable $Y = (\alpha_c^{2,1} - \alpha)^{1/2}$ sur l'intervalle $10 \le \alpha \le 13,5$, à savoir $B_{2,2}(0^+) = -0,73742 + 0,061673Y + 0,074061Y^2 - 0,0084386Y^3$. Tireté rouge : idem sur l'intervalle $12,5 \le \alpha \le 13,5$, à savoir $B_{2,2}(0^+) = -0,73761 + 0,062206Y + 0,073901Y^2 - 0,0086514Y^3$.

**TABLE 4** – Valeurs numériques du quatrième coefficient d'amas $B_{2,2}(0^+)$ du gaz unitaire de fermions dans un piège de raideur infinitésimale, tabulées en fonction du rapport de masse $\alpha = m_1/m_1$. Comme le coefficient est invariant par changement de $\alpha$ en $1/\alpha$, nous nous sommes limités à $\alpha \ge 1$. Les incertitudes, non précisées, sont inférieures au pour cent. Ces valeurs sont représentées graphiquement sur la figure 6.

<table><thead><tr><th>α</th><th>1</th><th>1,5</th><th>2</th><th>2,5</th><th>3</th><th>3,5</th><th>4</th><th>4,5</th><th>5</th><th>5,5</th></tr></thead><tbody><tr><th>B<sub>2,2</sub></th><td>-0,03056</td><td>-0,03529</td><td>-0,04518</td><td>-0,05777</td><td>-0,07233</td><td>-0,08845</td><td>-0,1060</td><td>-0,1248</td><td>-0,1449</td><td>-0,1662</td></tr><tr><th>α</th><td>6</td><td>6,5</td><td>7</td><td>7,5</td><td>8</td><td>8,5</td><td>9</td><td>9,5</td><td>10</td><td>10,5</td></tr><tr><th>B<sub>2,2</sub></th><td>-0,1886</td><td>-0,2122</td><td>-0,2370</td><td>-0,2629</td><td>-0,2900</td><td>-0,3183</td><td>-0,3478</td><td>-0,3788</td><td>-0,4110</td><td>-0,4448</td></tr><tr><th>α</th><td>11</td><td>11,5</td><td>12</td><td>12,5</td><td>12,75</td><td>13</td><td>13,1</td><td>13,2</td><td>13,3</td><td>13,4</td></tr><tr><th>B<sub>2,2</sub></th><td>-0,4801</td><td>-0,5177</td><td>-0,5575</td><td>-0,6004</td><td>-0,6235</td><td>-0,6484</td><td>-0,6590</td><td>-0,6701</td><td>-0,6819</td><td>-0,6948</td></tr></tbody></table>

## Remerciements

Ce travail a bénéficié du soutien financier japonais « Grants-in-Aid for Scientific Research » KAKENHI sous les numéros 21H00116 et 22K03492.
# Annexe A. Approximant asymptotique de $I_{3,1}^{(\ell,\epsilon)}(\bar{\omega})$, sa somme sur $\ell$ et $\epsilon$, son composition dominant

## A.1. L'approximant asymptotique

Pour obtenir une approximation $J_{3,1}^{(\ell,\epsilon)}(\bar{\omega})$ à grand $\ell$ de la contribution $I_{3,1}^{(\ell,\epsilon)}(\bar{\omega})$ de moment cinétique $\ell$ et de parité $\epsilon$ à la quantité $I_{3,1}(\bar{\omega})$ de l'équation (4), nous utilisons les notations de la section 2 et partons du développement (19). Un calcul astucieux de la trace sur $x > 0$ permet de regrouper des contributions faisant apparaître $x$ et $-x$ et de se ramener à des intégrales sur la droite réelle tout entière, dans lesquelles il ne reste à prendre que les traces sur les variables $u$ et $m_z$ :

$$ \ln \det M_{3,1}^{(\ell,\epsilon)} = \text{cte} - \int_{\mathbb{R}} dx \, e^{iSx} \text{Tr}_{u,m_z}^{\epsilon} \left[ \langle x | \bar{K}_{3,1}^{(\ell)} | -x \rangle e^{i\pi L_x / \hbar} \right] \\ + \int_{\mathbb{R}^2} dx \, dx' \, e^{iSx} \text{Tr}_{u,m_z}^{\epsilon} \left[ \langle x | \bar{K}_{3,1}^{(\ell)} | x' \rangle \langle x' | \bar{K}_{3,1}^{(\ell)} | -x \rangle e^{i\pi L_x / \hbar} \right] \\ - \frac{1}{2} \int_{\mathbb{R}^2} dx \, dx' \, e^{iS(x-x')} \text{Tr}_{u,m_z}^{\epsilon} \left[ \langle x | \bar{K}_{3,1}^{(\ell)} | x' \rangle e^{i\pi L_x / \hbar} \langle -x' | \bar{K}_{3,1}^{(\ell)} | -x \rangle e^{i\pi L_x / \hbar} \right] + \dots \quad (45) $$

où l'opérateur de rotation d'axe Ox d'angle $\pi$ provient de la note 9 en bas de page et nous avons introduit le noyau primitif (16) divisé à gauche par la partie diagonale (15) et pris à exposant d'échelle nul, c'est-à-dire $\bar{K}_{3,1}^{(\ell)} \equiv \mathcal{D}_{3,1}^{-1} K_{3,1}^{(\ell)}(iS = 0)$. Les contributions gardées au second membre de l'équation (45) constituent une fonction $-\phi^{(\ell,\epsilon)}(S)$ de S paire, régulière, à valeurs réelles et à décroissance rapide. Aussi peut-on intégrer par parties sur S dans l'équation (4) et faire apparaître la composante de Fourier de $\phi^{(\ell,\epsilon)}(S)$ à la pulsation $\bar{\omega}$ :

$$ J_{3,1}^{(\ell,\epsilon)}(\bar{\omega}) = \frac{(\ell + 1/2)\bar{\omega}}{\text{sh}\bar{\omega}} \int_{-\infty}^{+\infty} \frac{dS}{2\pi} e^{-i\bar{\omega}S} \phi^{(\ell,\epsilon)}(S) \quad (46) $$

facile à obtenir à partir de (45) au moyen de l'identité au sens des distributions $\int_{\mathbb{R}} \exp(ikS) dS = 2\pi\delta(k)$. Un calcul long mais sans difficulté, traitant les contributions à $\phi^{(\ell,\epsilon)}(S)$ dans l'ordre où elles apparaissent, donne finalement :

$$ J_{3,1}^{(\ell,\epsilon)}(\bar{\omega}) = \frac{(\ell + 1/2)\bar{\omega}}{\pi \text{sh}\bar{\omega} \sqrt{2\text{ch}\bar{\omega}}} \int_{-1}^{1} \frac{du}{\mathcal{D}_{3,1}(\bar{\omega}, u)} \int_{0}^{2\pi} \frac{d\phi}{2\pi} \frac{\mathcal{F}^{(\ell,\epsilon)}(\theta; \phi + \pi)}{1 + 2\text{ch}(2\bar{\omega}) + \frac{2\alpha}{1+\alpha}(2u\text{ch}\bar{\omega} + u^2 + v^2 \cos\phi)} \\ - \frac{(\ell + 1/2)\bar{\omega}}{\pi^2 \text{sh}\bar{\omega} \sqrt{2\text{ch}\bar{\omega}}} \int_{\mathbb{R}} dx' \int_{-1}^{1} du \, du' \int_{0}^{2\pi} \frac{d\phi \, d\phi'}{(2\pi)^2} \\ \frac{(\lambda'^{2+1/2}/\sqrt{2\text{ch}x'}) [\mathcal{D}_{3,1}(\bar{\omega}, u) \mathcal{D}_{3,1}(x', u')]^{-1} \mathcal{F}^{(\ell,\epsilon)}(\theta; \phi + \pi)}{[1 + \lambda_{\bar{\omega}}^2 + \lambda'^2 + \frac{2\alpha}{1+\alpha}(u\lambda_{\bar{\omega}} + u'\lambda' + \lambda_{\bar{\omega}}\lambda'(uu' + vv' \cos\phi))] [(\lambda_{\bar{\omega}}\phi, \phi) \rightarrow (\lambda_{\bar{\omega}}^{-1}\phi, \phi')] \\ + \frac{(\ell + 1/2)\bar{\omega}}{2\pi^2 \text{sh}\bar{\omega}} \int_{\mathbb{R}} dX \int_{-1}^{1} du \, du' \int_{0}^{2\pi} \frac{d\phi \, d\phi'}{(2\pi)^2} \\ \frac{[2\sqrt{\text{ch}X_+ \text{ch}X_-} \mathcal{D}_{3,1}(X_+, u) \mathcal{D}_{3,1}(X_-, u')]^{-1} \mathcal{F}^{(\ell,\epsilon)}(\theta, \theta'; \phi, \phi')}{[1 + \lambda_+^2 + \lambda_-^2 + \frac{2\alpha}{1+\alpha}(\lambda_+ u + \lambda_- u' + \lambda_+ \lambda_-(uu' + vv' \cos\phi))] [(\lambda_\pm, \phi) \rightarrow (\lambda_\pm^{-1}, \phi')]} \quad (47) $$
avec les notations $X_{\pm} = X \pm \bar{\omega}/2$, $\lambda_{\bar{\omega}} = \exp(\bar{\omega})$, $\lambda_{\pm} = \exp(X_{\pm})$ (ces deux dernières calquées sur $\lambda = \exp(x)$) complétant celles de l'équation (17) et les fonctions angulaires<fnref n="21" />

$$ \mathcal{F}^{(\ell, \varepsilon)}(\theta; \phi) \equiv \mathrm{Tr}_{m_z}^{\varepsilon} \left( e^{-i\theta L_z / \hbar} e^{i\phi L_x / \hbar} \right) \\ \text{et } \mathcal{F}^{(\ell, \varepsilon)}(\theta, \theta'; \phi, \phi') \equiv \mathrm{Tr}_{m_z}^{\varepsilon} \left( e^{-i\theta L_z / \hbar} e^{i\phi L_x / \hbar} e^{i\theta' L_z / \hbar} e^{i\phi' L_x / \hbar} \right) \qquad (48) $$

Les sommes sur la parité $\varepsilon$ des quantités définies dans (48) ont des expressions simples en termes d'angles $\xi \in [0, \pi]$, voir la référence [40]<fnref n="22" />,

$$ \mathcal{F}^{(\ell)}(\theta; \phi) = \frac{\sin[(\ell + 1/2)\xi]}{\sin(\xi/2)} \quad \text{avec} \quad 1 + 2\cos\xi = u(1 + \cos\phi) + \cos\phi \qquad (49) $$

$$ \mathcal{F}^{(\ell)}(\theta; \phi) = \frac{\sin[(\ell + 1/2)\xi]}{\sin(\xi/2)} \quad \text{avec} \quad (50) $$
$$ 1 + 2\cos\xi = uu'(1 + \cos\phi\cos\phi') - (u+u')\sin\phi\sin\phi' + vv'(\cos\phi + \cos\phi') + \cos\phi\cos\phi' $$

qui montre aussi comment revenir à parité fixée, par exemple

$$ \mathcal{F}^{(\ell, \varepsilon)}(\theta; \phi) = \frac{1}{2} \sum_{n=0}^{1} \varepsilon^n \mathcal{F}^{(\ell)}(\theta + n\pi; \phi) \qquad (51) $$

Nous préférons cependant réserver les équations ((49), (50), (51)) à des études analytiques; pour le calcul numérique, nous évaluons les traces de (48) dans la base propre de $L_z$ (sous la contrainte $(-1)^{m_z} = \varepsilon$) après insertion de relations de fermeture dans la base propre de $L_x$ (obtenue par diagonalisation numérique de sa matrice tridiagonale dans la base des $|\ell, m_z\rangle$) à l'endroit des opérateurs de rotation d'angle $\phi$ ou $\phi'$, puis utilisons la valeur de l'intégrale<fnref n="23" />

$$ \int_0^{2\pi} \frac{d\phi}{2\pi} \frac{e^{in\phi}}{b_0 + b_1 \cos\phi} = \frac{1}{\sqrt{b_0^2 - b_1^2}} \left( \frac{-b_1}{b_0 + \sqrt{b_0^2 - b_1^2}} \right)^{|n|} \quad \forall n \in \mathbb{Z}, \forall b_0 > 0, \forall b_1 \in ]-b_0, b_0] \quad (52) $$

## A.2. Sa somme sur $\ell$ et $\varepsilon$

La somme de $\int_{3,1}^{(\ell,\varepsilon)} (\bar{\omega})$ sur la parité $\varepsilon$ est tout de suite faite, voir les équations ((49), (50)). La somme sur $\ell \in \mathbb{N}$ fait apparaître, au sens des distributions, la série de Fourier d'un peigne de Dirac<fnref n="24" /> donc se calcule grâce à l'identité :

$$ \sum_{\ell=0}^{+\infty} (2\ell + 1) \frac{\sin[(\ell + 1/2)\xi]}{\sin(\xi/2)} = \frac{\pi\delta(\xi)}{\sin^2(\xi/2)} \qquad (53) $$

$$ \sum_{\ell \in \mathbb{N}} (\ell + 1/2) \{\cos(\ell\xi) - \cos[(\ell + 1)\xi]\} = (-1/2) + \sum_{\ell \in \mathbb{N}} \cos(\ell\xi). $$

<footnote n="21">21 Un lecteur attentif objectera que la définition de $\mathcal{F}^{(\ell, \varepsilon)}(\theta, \theta'; \phi, \phi')$ devrait en principe comporter un projecteur orthogonal sur le sous-espace de parité $(-1)^{L_z/\hbar} = \varepsilon$ à côté de l'opérateur $e^{i\theta' L_z/\hbar}$. La parité en $\phi$ et en $\phi'$ du dénominateur dans la troisième contribution à (47) permet cependant d'en faire l'économie. Par exemple, seule la partie paire de $e^{i\phi L_x/\hbar}$ contribue après intégration sur $\phi$; or, elle préserve la parité de $L_z/\hbar$.</footnote>

<footnote n="22">22 En passant aux arcs moitiés, on peut remarquer qu'on a plus simplement $\cos(\xi/2) = \cos(\theta/2)|\cos(\phi/2)|$ dans l'équation (49).</footnote>

<footnote n="23">23 Pour optimiser, on peut (i) remplacer $\exp[i(-m_z\theta + m'_z\theta')]$ par sa partie réelle $\cos(m_z\theta)\cos(m'_z\theta') + \sin(m_z\theta)\sin(m'_z\theta')$ et tabuler les sinus et cosinus correspondants, (ii) se ramener à une boucle sur un seul indice $m_x = m'_x$ et utiliser la symétrie par changement de $m_x$ en $-m_x$, (iii) utiliser la symétrie par changement de $(m_z, m'_z)$ en $(-m_z, -m'_z)$ (on se limite à $(-1)^{m'_z} = \varepsilon$, voir la note 21), (iv) tabuler les puissances de degré $|m_x|$ apparaissant dans (52).</footnote>

<footnote n="24">24 On pourra simplement penser à écrire $\sin[(\ell + 1/2)\xi] = \{\cos(\ell\xi) - \cos[(\ell + 1)\xi]\}/[2\sin(\xi/2)]$ pour faire apparaître la somme partiellement télescopique</footnote>
Il faut donc trouver les cas d'annulation de l'angle $\xi$, sachant que les angles polaires $\theta, \theta'$ (désormais pris comme variables d'intégration de préférence à $u, u'$) décrivent $[0, \pi]$ et que les angles azimutaux $\phi, \phi'$ peuvent être pris commodément dans $[-\pi, \pi]$. Dans l'équation (49), on trouve comme seul point d'annulation $(\theta, \phi) = (0, 0)$; au voisinage de ce point, $\xi^2 \approx \theta^2 + \phi^2$ si bien que l'action de $\delta(\xi)$ s'évalue facilement en coordonnées polaires $(\rho, \gamma)$ dans le plan $(\theta, \phi)$, pour une fonction régulière f quelconque :

$$ \int_{0}^{\pi} d\theta \sin\theta \int_{-\pi}^{\pi} \frac{d\phi}{2\pi} f(\theta, \phi) \frac{\pi\delta(\xi)}{\sin^2(\xi/2)} = f(0,0) \int_{0}^{+\infty} d\rho \rho \int_{-\pi/2}^{\pi/2} d\gamma \rho \cos\gamma \frac{\delta(\rho)}{2(\rho/2)^2} = 2f(0,0) \quad (54) $$

Dans l'équation (50), on trouve trois lignes d'annulation<fnref n="25" />:

$$ l_1 : \theta = \theta' = 0, \phi + \phi' = 0 \quad ; \quad l_2 : \theta = \theta' = \pi, \phi - \phi' = 0; \quad l_3 : \theta = \theta', \phi = \phi' = 0 \quad (55) $$

Dans l'intégrale à laquelle nous sommes confrontés, qui fait intervenir une fonction régulière g des quatre angles,

$$ \mathcal{J} = \int_{0}^{\pi} d\theta \sin\theta \int_{0}^{\pi} d\theta' \sin\theta' \int_{-\pi}^{\pi} \frac{d\phi}{2\pi} \int_{-\pi}^{\pi} \frac{d\phi'}{2\pi} g(\theta, \theta', \phi, \phi') \frac{\pi\delta(\xi)}{\sin^2(\xi/2)} \quad (56) $$

les lignes $l_1$ et $l_2$ ont une contribution nulle à cause des jacobiens $\sin\theta$ et $\sin\theta'$. Au voisinage de la ligne $l_3$, on quadratise la quantité $\xi$ en les variables $\phi, \phi'$ et $\theta' - \theta$; l'intégrale à $\theta$ fixé

$$ \int_{-\infty}^{+\infty} d\phi \int_{-\infty}^{+\infty} d\phi' \int_{-\infty}^{+\infty} d(\theta' - \theta) \frac{\delta\left(\sqrt{\phi^2 + \phi'^2 + 2\phi\phi' \cos\theta + (\theta' - \theta)^2}\right)}{\phi^2 + \phi'^2 + 2\phi\phi' \cos\theta + (\theta' - \theta)^2} = \frac{2\pi}{\sin\theta} \quad (57) $$

se calcule bien dans la base propre de la forme quadratique apparaissant au dénominateur et sous la racine carrée (elle est de valeurs propres 1 et $1 \pm \cos\theta$). Il reste simplement

$$ \mathcal{J} = 2 \int_{0}^{\pi} d\theta \sin\theta g(\theta, \theta, 0, 0) \quad (58) $$

d'où l'expression de la somme désirée (en respectant l'ordre des contributions dans l'équation (47))<fnref n="26" />:

$$ \begin{aligned} \sum_{\ell \in \mathbb{N}} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) = & \frac{\bar{\omega}}{\pi \sinh \bar{\omega} \sqrt{2 \cosh \bar{\omega}}} \frac{[\mathcal{D}_{3,1}(\bar{\omega}, u=1)]^{-1}}{1 + 2 \cosh(2\bar{\omega}) + \frac{2\alpha}{1+\alpha}(1 + 2 \cosh \bar{\omega})} \\ & - \frac{\bar{\omega}}{2\pi^2 \sinh \bar{\omega} \sqrt{\cosh \bar{\omega}}} \int_{\mathbb{R}} dx' \int_{-1}^{1} du' \frac{\lambda'^{2+1/2}}{\sqrt{\cosh x'}} \frac{[\mathcal{D}_{3,1}(\bar{\omega}, u=1) \mathcal{D}_{3,1}(x', u')]^{-1}}{[1 + \lambda_{\bar{\omega}}^2 + \lambda'^2 + \frac{2\alpha}{1+\alpha}(\lambda_{\bar{\omega}} + u'\lambda' + \lambda_{\bar{\omega}}\lambda' u')] [\lambda_{\bar{\omega}} \to \lambda_{\bar{\omega}}^{-1}]} \\ & + \frac{\bar{\omega}}{4\pi^2 \sinh \bar{\omega}} \int_{\mathbb{R}} dX \int_{-1}^{1} du \frac{[\sqrt{\cosh X_+ \cosh X_-} \mathcal{D}_{3,1}(X_+, u) \mathcal{D}_{3,1}(X_-, u)]^{-1}}{[1 + \lambda_+^2 + \lambda_-^2 + \frac{2\alpha}{1+\alpha}((\lambda_+ + \lambda_-)u + \lambda_+\lambda_-)] [\lambda_\pm \to \lambda_\pm^{-1}]} \end{aligned} \quad (59) $$

<footnote n="25">²⁵Si $\xi = 0$, $\mathcal{F}^{(1)}(\theta, \theta'; \phi, \phi') = 3$; or, c'est la trace d'une matrice de rotation dans $\mathbb{R}^3$ qui doit donc se réduire à l'identité : on a $\mathcal{R}_z(\theta)\mathcal{R}_x(-\phi)\mathcal{R}_z(-\theta')\mathcal{R}_x(-\phi') = \mathbf{1}$ où $\mathcal{R}_x(\theta)$ est la rotation d'angle $\theta$ d'axe $Ox$, etc. En prenant l'élément de matrice $zx$ de cette relation, on aboutit à la condition $\sin\theta' \sin\phi = 0$. De même, en changeant l'ordre des opérateurs par permutation circulaire sous la trace, on aboutit à $\sin\theta \sin\phi' = 0$. Les lignes d'annulation $l_1$ et $l_2$ dans (55) correspondent à $\sin\theta = \sin\theta' = 0$ c'est-à-dire $\theta = \theta' = 0|\pi$, les cas $\theta = 0, \theta' = \pi$ et $\theta = \pi, \theta' = 0$ ne convenant trivalement pas. La ligne $l_3$ correspond à $\sin\phi = \sin\phi' = 0$ c'est-à-dire $\phi = \phi' = 0|\pi$, qui se réduit en pratique à $\phi = \phi' = 0$; en effet, les cas $\phi = 0, \phi' = \pi$ et $\phi = \pi, \phi' = 0$ ne conviennent visiblement pas, et le cas $\phi = \phi' = \pi$ impose $\cos(\theta + \theta') = 1$ donc $\theta = \theta' = 0$ ou $\pi$, ce qui correspond à des points d'annulation, de contribution nulle à l'intégrale sur $\theta, \theta', \phi, \phi'$. Les cas restants $\sin\theta = \sin\phi = 0$ ou $\sin\theta' = \sin\phi' = 0$ ne conduisent qu'à des points d'annulation et ne contribuent pas non plus.</footnote>

<footnote n="26">²⁶Là où apparaît $\phi + \pi$ en argument de $\mathcal{F}^{(\ell)}$, on effectue le changement de variable $\phi \to \phi - \pi$.</footnote>
En pratique, le cas le plus important correspond à $\tilde{\omega} = 0^+$, car c'est celui qui donne accès au co-efficient d'amas $b_{3,1}$ du gaz homogène. Nous donnons l'expression explicite de (59) correspondante, toujours en respectant l'ordre des contributions, en termes de fonctions trigonométriques et du dilogarithme ou fonction de Bose $g_2(z) = \sum_{n>0} z^n/n^2$ :

$$ 
\begin{aligned}
& \sum_{\ell \in \mathbb{N}} \sum_{\varepsilon} J_{3,1}^{(\ell, \varepsilon)}(\tilde{\omega} = 0^+) \\
&= \frac{1}{3\pi\sqrt{2}} \frac{(1+\alpha)^2}{(1+3\alpha)^{3/2}} - \frac{(1+\alpha)^5}{\sqrt{2}\pi^2\alpha(1+2\alpha)(1+3\alpha)^{3/2}} \left[ \arctan \frac{\alpha}{\sqrt{1+2\alpha}} - \frac{\alpha \operatorname{argh} \frac{\sqrt{1+3\alpha}}{\sqrt{2(1+2\alpha)}}}{(\alpha+1)\sqrt{2(1+3\alpha)}} \right] \\
&\quad + \frac{(1+\alpha)^4}{2\pi^2\alpha(1+3\alpha)^2} \left\{ \pi \arcsin \frac{\alpha}{1+2\alpha} + 2 \operatorname{Re} \left[ g_2\left(\frac{-1}{\zeta}\right) - g_2\left(\frac{1}{\zeta}\right) \right] \right. \\
&\qquad \left. - \operatorname{Re} \left[ g_2\left(\frac{2\zeta}{1+\zeta}\right) - g_2\left(\frac{-2\zeta}{1-\zeta}\right) + g_2\left(\frac{2}{1-\zeta}\right) - g_2\left(\frac{2}{1+\zeta}\right) \right] \right\}
\end{aligned} 
\quad (60) \\
\text{où } \zeta = (2\alpha + i\sqrt{2(1+3\alpha)})/(1+\alpha).
$$

### A.3. *Son comportement dominant*

Pour obtenir un équivalent de $J_{3,1}^{(\ell, \varepsilon)}(\tilde{\omega})$ ou, ce qui revient au même, de $I_{3,1}^{(\ell, \varepsilon)}(\tilde{\omega})$ à grand moment cinétique, il suffit de se limiter à la première contribution au second membre de l'équation (47), qui provient du terme d'ordre un en le noyau $\mathcal{K}_{3,1}^{(\ell)}$ dans le développement (19). On l'écrit alors sous forme d'une intégrale sur l'angle $\xi$ en insérant un delta de Dirac $\delta(\xi - \xi(u, \phi))$ reliant $\xi$ à $u$ et $\phi$ comme dans l'équation (49), puis en calculant explicitement l'intégrale sur $\phi$ puis sur $u$.<fnref n="27" />
Pour simplifier, nous nous limitons à $\tilde{\omega} = 0^+$ et trouvons :

$$ 
\begin{gathered}
I_{3,1}^{(\ell, \varepsilon)}(0^+) \underset{\ell \to +\infty}{\sim} \int_0^\pi d\xi \frac{\sin[(\ell + 1/2)\xi]}{2 \sin(\xi/2)} [\rho_0(\xi) + \varepsilon \rho_1(\xi)] \\
\text{où } \rho_n(\xi) = (2\ell + 1) A_n \frac{\sin^2(\xi/2)}{w_n(\cos\xi)} F_n(f_n(\cos\xi)) \quad (61)
\end{gathered}
$$

avec la notation $\beta = 1/\alpha$, les préfacteurs et les fonctions auxiliaires<fnref n="28" />

$$ 
\begin{gathered}
A_0 = \frac{\sqrt{2}(1+\alpha)^2}{\pi^2\alpha\sqrt{1+3\alpha}} \quad ; \quad w_0(Z) = 3\beta + (1+2Z)^2 \quad ; \\
F_0(Z) = \frac{\arctan\sqrt{Z}}{\sqrt{Z}} \quad ; \quad f_0(Z) = \frac{(1-Z)(5+4Z)}{w_0(Z)}
\end{gathered} 
\quad (62) 
$$

$$ 
\begin{gathered}
A_1 = \frac{\sqrt{2}(1+\alpha)^{5/2}}{\pi^2\alpha^2} \quad ; \\
w_1(Z) = 8(1+\beta)Z^2 - (5\beta+7)Z + \beta(3\beta+1) \quad ;
\end{gathered} 
\quad (63) 
$$

$$ F_1(Z) = Z \left( \frac{\pi}{2} - \arctan Z \right) $$

<footnote n="27">²⁷Après s'être ramené à une intégrale sur $\phi \in [0, \pi]$ par périodicité et parité de l'intégrande, on utilise la relation $\delta(2\cos\xi + 1 - u(1 + \cos\phi) - \cos\phi) = \delta(\phi - \phi_0)/[(1+u)\sin\phi_0]$ avec $\sin\phi_0 = 2(u - \cos\xi)^{1/2}(1 + \cos\xi)^{1/2}/(1+u)$ où la racine $\phi_0$ est dans $[0, \pi]$ si et seulement si $0 \le \theta \le \xi$ compte tenu de la note 22. On utilise aussi la relation (51).</footnote>

<footnote n="28">²⁸ Quand nous passerons dans le plan complexe, il faudra comprendre que $\arctan Z = \ln[(1+iZ)/(1-iZ)]/(2i)$. Du coup, la ligne de coupure de la fonction $F_0(Z)$ est $]-\infty, -1]$ et celle de $F_1(Z)$ est $]-\infty, -1] \cup i[1, +\infty]$. On notera que, malgré les apparences, $F_0(Z)$ est analytique dans le disque unité, comme le montre son développement en série en $Z=0$. Remarquons aussi qu'avec cette définition de la fonction arctan, on a $\pi/2 - \arctan(ix) = -i\argth(1/x) \forall x \in \mathbb{R} \setminus [-1, 1]$, ce qui sert déjà dans l'intégrale sur l'axe réel (61) car la quantité sous la dernière racine carrée au dénominateur de (64) peut devenir négative; on a fait ici les choix $\ln(-1) = i\pi$ et $\sqrt{-1} = i$.</footnote>
$$f_1(Z) = \frac{w_1(Z)}{4\sqrt{1+\beta}\sqrt{1-Z}\left\{\left[\frac{3\beta+1}{4} + (1+\beta)Z\right]\left[3\beta + (1-2Z)^2\right]\right\}^{1/2}} \quad (64)$$

Pour se ramener à des intégrales de contour sur le cercle unité C, de la forme

$$\mathcal{I}_n = \int_C \frac{dz}{iz} (z-1) z^\ell \frac{F_n\left(f_n\left(\frac{z+1/z}{2}\right)\right)}{w_n\left(\frac{z+1/z}{2}\right)} \quad (65)$$

on procède en trois étapes : (i) on élimine les facteurs $\sin(\xi/2)$ dans l'intégrande de (61), d'abord par simplification du sinus au dénominateur devant les crochets avec un sinus au numérateur de $\rho_n(\xi)$, puis en faisant apparaître $\cos[(\ell + 1)\xi] - \cos(\ell\xi)$ avec le sinus restant, comme dans la note 24; (ii) on étend l'intégrale sur $\xi$ à $[-\pi, \pi]$ par parité de l'intégrande; (iii) on prend comme nouvelle variable d'intégration $z = \exp(i\xi)$, qui décrit le cercle unité, si bien que $\cos\xi = (z+1/z)/2$, et l'on peut écrire $\cos[(\ell + 1)\xi] - \cos(\ell\xi) = \text{Re}(z^{\ell+1} - z^\ell)$ et sortir la fonction Re de l'intégrale puisque le reste de l'intégrande est réel sur C. Pour les déterminations habituelles du logarithme et de la racine carrée dans le plan complexe (ligne de coupure sur $\mathbb{R}^-$), on trouve que l'intégrande de $\mathcal{I}_0$ admet comme singularité dans le disque unité une ligne de coupure en forme de trident, réunion d'un segment $OO'$ et d'un arc $AA'$, voir la figure 7a, et celui de $\mathcal{I}_1$ une ligne de coupure en forme de fleur, réunion de trois arcs $OB, OB'$ et $BB'$, voir la figure 7b$^{29}$. Le point O est ici l'origine des coordonnées, A et $A'$ sont les points d'affixes $z_0$ et $z_0^*$, B et $B'$ les points d'affixes $z_1 = -z_0^*$ et $z_1^*$, et $z_0$ est l'unique solution de l'équation $z+z^{-1}+1+i\sqrt{3\beta}=0$ dans le disque unité (si bien que $w_0(\frac{z_0+1/z_0}{2})=0$) :

$$z_0 = \frac{1}{2} \left[ i\sqrt{1 - i\sqrt{3\beta}}\sqrt{3 + i\sqrt{3\beta}} - (1 + i\sqrt{3\beta}) \right] \quad (66)$$

À l'aide du théorème intégral de Cauchy, nous resserrons le contour d'intégration jusqu'à ce qu'il épouse les lignes de coupure, sans que cela change la valeur de $\mathcal{I}_n$. Dans la limite $\ell \to +\infty$, à cause du facteur $z^\ell$, l'intégrale est alors dominée par le voisinage des points de singularité les plus éloignés de O, à savoir A et $A'$ pour $\mathcal{I}_0, B$ et $B'$ pour $\mathcal{I}_1$, ceci valant pour tout rapport de masse $\alpha$. $\mapsto$ Pour le calcul de $\mathcal{I}_0$ au voisinage de $z_0$, on pose $z = z_0(1 + \eta u)$, où $u$ est un nombre complexe et $\eta > 0$ est un infiniment petit, et on approxime la corne du trident par sa demi-tangente à droite au sommet A; on dispose alors des équivalents $f_0((z+1/z)/2) \sim C_0/(\eta u)$ et $z^\ell \sim z_0^\ell \exp(\ell\eta u)$, et la demi-tangente devient la demi-droite d'origine O de direction $-C_0$ dans l'espace des $u$ (compte tenu de la note 28) donc la demi-droite $\mathbb{R}^-$ dans l'espace des $v = \eta u/C_0$, si bien que

$$\mathcal{I}_0^2 \sim \frac{(z_0 - 1) z_0^\ell}{2\sqrt{3\beta}(z_0 - 1/z_0)} \int_{\mathcal{C}_u} \frac{du}{u} e^{\ell \eta u} F_0\left(\frac{C_0}{\eta u}\right) = \frac{z_0^{\ell+1}}{2\sqrt{3\beta}(1 + z_0)} \int_{\mathcal{C}_v} \frac{dv}{v} e^{\ell C_0 v} F_0(1/v) \quad (67)$$

où les lacets $\mathcal{C}_u$ et $\mathcal{C}_v$ entourent les demi-tangentes en question dans le sens trigonométrique dans l'espace des $u$ ou des $v$. L'expression de $C_0$ est donnée dans l'équation (21). L'intégrale au dernier membre de (67) se calcule aisément :

$$\begin{align}
\int_{\mathcal{C}_v} \frac{dv}{v} e^{\ell C_0 v} F_0(1/v) &= \int_{-\infty}^{0} \frac{dx}{x} e^{\ell x C_0} \left[ F_0\left(\frac{1}{x - i0^+}\right) - F_0\left(\frac{1}{x + i0^+}\right) \right] \nonumber \\
&= \int_{-1}^{0} \frac{dx}{x} e^{\ell x C_0} (-i\pi|x|^{1/2}) \stackrel{\ell \to +\infty}{\sim} \frac{i\pi^{3/2}}{\sqrt{\ell C_0}} \tag{68}
\end{align}$$

<footnote>$^{29}$Dans ce second cas, on s'attend à avoir trois lignes de coupure, celle $L_1$ venant de la fonction $F_1(Z)$, celles $L_3$ et $L_2$ venant de $\sqrt{1-Z}$ et de l'autre racine carrée au dénominateur de la fonction $f_1(Z)$. En réalité, $L_3$ est incluse dans $L_1$. De plus, les intersections $L_1 \cap L_2$ et $L_1 \cap L_3$ ne sont pas des lignes de coupure de l'intégrande, pour la raison que $z \mapsto z^{1/2}(\pi/2 - \arctan z^{1/2})$ admet seulement $[-1,0]$ comme ligne de coupure (pour $z \in [-\infty, -1]$, le changement de signe de $\pi/2$-arctan compense celui de la racine carrée). Les arcs $OB$ et $OB'$ sont constitués des points se trouvant seulement sur $L_2$, l'arc $BB'$ regroupe les points appartenant seulement à $L_1$.</footnote>
La corne $A'$ d'affiche $z_0^*$ donne une contribution complexe conjuguée de la corne A. Nous retrouvons la première partie de l'équation (20). $\rightarrow$ Pour le calcul de $J_1$ au voisinage de $z_1$, on procède de même, en posant $z = z_1(1 + \eta u)$, $\eta \to 0^+$. Cette fois le comportement local $f_1((z+1/z)/2) \sim A/(B\eta u)^{1/2}$ admet lui-même une ligne de coupure à cause de la racine carrée et est caractérisé par deux amplitudes complexes $A$ et $B$, d'où la structure à deux lignes de coupure, celle provenant de $f_1$ et celle provenant de arctan $f_1$ dans $F_1(f_1)^{30}$ approximée par deux demi-tangentes représentées dans l'espace des $v = \eta B u / A^2$ simplement par $A^{-2} \mathbb{R}^-$ et, compte tenu de la note 28, par $\mathbb{R}^-$. Cette fois

$$ J_1^{(\ell, \epsilon)} \xrightarrow{\ell \to +\infty} \frac{(z_1-1) z_1^{\ell} A^2}{i \mathcal{B} w_1 \left(\frac{z_1+1/z_1}{2}\right)} \int_{\mathcal{C}_v} dv e^{\ell A^2 v / \mathcal{B}} \frac{A}{\sqrt{v A^2}} \left( \frac{\pi}{2} - \arctan \frac{A}{\sqrt{v A^2}} \right) \quad (69) $$

où le lacet $\mathcal{C}_v$ entoure (dans le troisième quadrant) la demi-droite $A^{-2}\mathbb{R}^-$ dans le sens trigonométrique puis rejoint verticalement la demi-droite $\mathbb{R}^-$ par en dessous pour l'entourer elle aussi dans le sens trigonométrique. Le calcul conduit à$^{31}$

$$ \begin{aligned} & \int_{\mathcal{C}_v} dv e^{\ell A^2 v / \mathcal{B}} F_1 \left( \frac{A}{\sqrt{v A^2}} \right) \\ &= \int_{-\infty}^{0} \frac{dx}{A^2} e^{\ell x / \mathcal{B}} \left[ F_1 \left( \frac{A}{\sqrt{x - i0^+}} \right) - F_1 \left( \frac{A}{\sqrt{x + i0^+}} \right) \right] \\ &+ \int_{-\infty}^{0} dx e^{\ell A^2 x / \mathcal{B}} \left[ F_1 \left( \frac{A}{\sqrt{(x - i0^+)^2 A^2}} \right) - F_1 \left( \frac{A}{\sqrt{(x + i0^+)^2 A^2}} \right) \right] \xrightarrow{\ell \to +\infty} \\ & \qquad \int_{-\infty}^{0} \frac{dx}{A} e^{\ell x / \mathcal{B}} \frac{i\pi}{\sqrt{-x}} + \int_{-\infty}^{0} dx e^{\ell A^2 x / \mathcal{B}} \frac{(-i\pi)}{\sqrt{-x}} = \frac{i\pi^{3/2}}{\sqrt{\ell}} \left( \frac{\sqrt{\mathcal{B}}}{A} - \frac{1}{\sqrt{A^2 / \mathcal{B}}} \right) = \frac{2i\pi^{3/2} \mathcal{B}}{A \sqrt{\ell \mathcal{B}}} \end{aligned} \quad (70) $$

ce qui redonne la seconde partie de l'équation (20), le coefficient $\mathcal{B}$ étant appelé $C_1$ et donné par l'équation (21).

## Annexe B. Approximant asymptotique de $I_{2,2}^{(\ell, \epsilon)}(0^+)$, sa somme sur $\ell$ et $\epsilon$

### B.1. *L'approximant asymptotique*

Il s'agit ici d'écrire le plus explicitement possible une approximation $J_{2,2}^{(\ell, \epsilon)}(0^+)$ à grand $\ell$ de la contribution $I_{2,2}^{(\ell, \epsilon)}(0^+)$ de moment cinétique $\ell$ et de parité $\epsilon$ à la quantité $I_{2,2}(0^+)$ de l'équation (4). Nous prenons comme point de départ les développements ((38), (42)) selon la valeur de $\epsilon$. Leur report dans l'intégrande de (4) conduit, après intégration sur $S$, à la structure

$$ J_{2,2}^{(\ell, \epsilon)}(0^+) = \mathcal{J}_{K_1, K_2}^{(\ell, \epsilon)} + \delta_{\epsilon, (-1)^\ell} \left[ \mathcal{J}_{K_1, K_3}^{(\ell)} + \mathcal{J}_{K_2, K_3}^{(\ell)} \right] \quad (71) $$

où l'indice indique qu'il s'agit d'une contribution croisée des noyaux intégraux $K_i$ et $K_j$. Le cas $(i, j) = (1, 2)$ a été traité dans la référence [40] pour $\epsilon = (-1)^{\ell-1}$, en utilisant les formes (36) des

<footnote>$^{30}$Dans le cas précédent $n = 0$, ceci ne se produit pas : dans la fonction $F_0(z)$, de part et d'autre de la ligne de coupure $\mathbb{R}^-$ de la racine carrée, on hésite certes entre deux valeurs opposées pour $\sqrt{z}$ mais arctan est une fonction impaire. Aussi la ligne de coupure de $F_0(z)$ est-elle dictée par la fonction arctan et vaut $]-\infty, -1]$, en accord avec la note 28.</footnote>

<footnote>$^{31}$On vérifie en effet que $\sqrt{A^2/\mathcal{B}} = -A/\sqrt{\mathcal{B}}$ pour tout rapport de masse $\alpha$.</footnote>
<figure><img src="image_9.png" /><figcaption>FIGURE 7 – Lignes de coupure de l’intégrande de l’équation (65) dans le disque unité (trait rouge épais) dans le cas (a) $n = 0$ et (b) $n = 1$, obtenues numériquement pour un rapport de masse $\alpha = 1$ pris comme exemple.</figcaption></figure>

noyaux; nous nous contentons ici de l'étendre au cas $\varepsilon = (-1)^{\ell}$, ce qui n'est qu'un jeu d'écriture, et de recopier le résultat :

$$ \mathcal{J}_{K_1, K_2}^{(\ell, \varepsilon)} = \int_{\mathbb{R}} dx \int_{-1}^{1} du \, du' \int_{0}^{2\pi} \frac{d\phi d\phi'}{(2\pi)^2} \frac{[(2\ell+1)/(4\pi)^2][\operatorname{ch} x \mathcal{D}_{2,2}(x, u)\mathcal{D}_{2,2}(x, u')]^{-1} \mathcal{F}^{(\ell, \varepsilon)}(\theta, \theta'; \phi, \phi')}{\{1 + \frac{1}{1+\alpha}[(u+e^{-x})(u'+e^{-x})+vv' \cos \phi]\}\{1 + \frac{\alpha}{1+\alpha}[(u+e^x)(u'+e^x)+vv' \cos \phi']\}} \quad (72) $$

avec les notations (17) et (48). Le cas $(i, j) = (1, 3)$ est beaucoup plus simple à traiter (en particulier numériquement) en utilisant les formes (96) des noyaux, repérées par un accent tchèque. Un seul facteur dépend de l'exposant d'échelle iS, sur lequel il faut intégrer dans (4). Dans les notations de (96) il vient après une intégration par parties :

$$ \begin{aligned} \int_{\mathbb{R}} \frac{dS}{2\pi} (|z'_1|/|z_1|)^{iS} &= 2\delta(\ln|z'_1|^2 - \ln|z_1|^2) = 2|z_1| |z'_1| \delta(|z_1|^2 - |z'_1|^2) \\ &= 2|z_1| |z'_1| \delta(e^{2x} - e^{2x'} - \sqrt{2\alpha}e^x u + \sqrt{2\alpha}e^{x'} u') \end{aligned} \quad (73) $$

Nous décidons d'intégrer la distribution de Dirac sur $u'$ à $u$, $x$, $x'$ fixés, ce qui lie la valeur de $u'$ à ces autres variables :

$$ u' = e^{x-x'} u + \frac{e^{2x'} - e^{2x}}{\sqrt{2\alpha}e^{x'}} \quad (74) $$

Comme $u'$ doit être compris entre $-1$ et $1$, l'intervalle d'intégration sur $u$ est contraint et se réduit à $[u_{\min}, u_{\max}]$ avec

$$ \begin{aligned} u_{\min} &= \max \left( -1, \frac{e^{2x} - e^{2x'}}{\sqrt{2\alpha}e^x} - e^{x'-x} \right) ; \\ u_{\max} &= \min \left( 1, \frac{e^{2x} - e^{2x'}}{\sqrt{2\alpha}e^x} + e^{x'-x} \right) ; \quad u_{\min} < u_{\max} \Leftrightarrow |e^{x'} - e^x| < \sqrt{2\alpha} \end{aligned} \quad (75) $$

Pour calculer le noyau resommé (41), nous reconnaissions dans la partie externe de $\tilde{K}_3$ (dernière contribution de (96)) une forme factorisée en un opérateur sur l'espace des $x$ et un projecteur
orthogonal de rang un sur l'espace des $u$, ce qui réduit le problème à l'inversion d'un opérateur sur l'espace des $x$ seulement$^{32}$ :

$$ \check{K}_3^{\text{ext}} = \hat{k} \otimes |\chi\rangle\langle\chi| \quad \text{et} \quad \langle\chi|\chi\rangle = 1 \implies \check{K}_3^{\text{ext}} = [\mathbf{1} - (\mathbf{1} + \hat{k})^{-1}] \otimes |\chi\rangle\langle\chi| \quad (76) $$

avec

$$ \langle u | \chi \rangle = \frac{1}{\sqrt{2}} \quad \text{et} \quad \langle x | \hat{k} | x' \rangle = \frac{(-1)^{l+1}}{\pi} \frac{2e^{x+x'}}{e^{2x} + e^{2x'} + 1} \frac{1}{(1 + e^{-2x})^{1/4} (1 + e^{-2x'})^{1/4}} \quad (77) $$

Nous obtenons finalement

$$ \mathcal{J}_{K_1, K_3}^{(\ell)} = \frac{(1 + \alpha)^3}{8\pi\alpha^2} \int_{-\infty}^{+\infty} dx \int_{-\infty}^{+\infty} dx' \frac{\mathcal{R}(x, x') \langle x' | \mathbf{1} - (\mathbf{1} + \hat{k})^{-1} | x \rangle}{(1 + e^{-2x})^{1/4} (1 + e^{-2x'})^{1/4}} \quad (78) $$

où nous avons introduit la fonction symétrique de $x$ et $x'$ tenant compte de (75) par une fonction de Heaviside $Y$ :

$$ \mathcal{R}(x, x') \equiv Y(\sqrt{2\alpha} - |e^{x'} - e^x|) \int_{u_{\min}}^{u_{\max}} du \frac{e^x}{|z_1|} \int_0^{2\pi} \frac{d\phi}{2\pi} \frac{(2\ell + 1) \langle l, m_x = 0 | e^{i\tau_1 L_z/\hbar} e^{i\phi L_x/\hbar} e^{-i\tau'_1 L_z/\hbar} | l, m_x = 0 \rangle}{|z_2|^2 + |z'_2|^2 + (1 + \alpha^{-1}) (\cos \tau_1 \cos \tau'_1 + \sin \tau_1 \sin \tau'_1 \cos \phi)} \quad (79) $$

Les nombres complexes $z_i$ et $z'_i$ et leurs arguments $\tau_i$ et $\tau'_i$ sont définis après l'équation (94) en fonction des angles $\theta$ et $\theta'$ dans l'intervalle $[0, \pi]$ tels que $u = \cos\theta$ et $u' = \cos\theta'$, l'angle $\theta'$ étant relié à $\theta$ par l'équation (74); l'écriture du dénominateur de l'intégrande dans (79) tient compte de l'égalité des modules $|z'_i| = |z_1|$.33 Enfin, la dernière contribution dans (71) se déduit de la seconde contribution en changeant partout $\alpha$ en son inverse $1/\alpha$ :

$$ \mathcal{J}_{K_2, K_3}^{(\ell)}(\alpha) = \mathcal{J}_{K_1, K_3}^{(\ell)}(1/\alpha) \quad (80) $$

<footnote>$^{32}$Il n'y a pas ici de changement de jauge de type (39) à effectuer car $\check{K}_3$ est directement indépendant de l'exposant d'échelle. On remarque aussi sur l'équation (95) que $\mathcal{G}_{22} \equiv 1$.</footnote>

<footnote>$^{33}$En pratique, l'intégrale angulaire dans (79) est faite numériquement sur l'angle $\theta$ (plutôt que sur $u$) avec la méthode de Gauss-Legendre à 41 points; si $x$ est assez proche mais distinct du point de singularité $x_{\text{sing}} = \ln \sqrt{\alpha/2}$ où le premier dénominateur $|z_1|$ peut s'annuler, c'est-à-dire $0 < |x - x_{\text{sing}}| < \theta_c$ avec $\theta_c = 0, 15$ par exemple, une structure étroite en $\theta$ apparaît et l'intégration sur $[\theta_{\min}, \theta_c]$ est effectuée avec le changement de variable $\theta = |x - x_{\text{sing}}|$ sh $t$, avec la méthode du point milieu sur la variable $t$ discrétisée en 100 points équidistants. En revanche, le second dénominateur dans (79) reste toujours supérieur à $1 + \min(\alpha, \alpha^{-1})$ sur le domaine d'intégration et ne peut s'annuler. L'intégration sur $x$ et $x'$ est faite avec la méthode du point milieu sur un intervalle de troncature $[x_{\min}, x_{\max}]$, en s'arrangeant pour que $x_{\text{sing}}$ soit le bord de l'un des sous-intervalles de largeur $dx$; on fait le choix fixe $x_{\min} = -5$ car $\mathcal{R}(x, x') = O(\min(\exp(x), \exp(x')))$ lorsque $x, x' \to -\infty$, mais on extrapole à $x_{\max} = +\infty$ cubiquement en $1/x_{\max}$ à partir des quatre valeurs $x_{\max}^{(0)} = 7$, $x_{\max}^{(1)} = 30$, $x_{\max}^{(2)} = 52$ et $x_{\max}^{(3)} = 75$. Pour tirer parti du fait que $\mathcal{R}(x, x') = O(\exp(-2x))$ sur un support de plus en plus étroit $x - x' = O(\exp(-x))$ lorsque $x \to +\infty$, on sépare la matrice $A$, version discrétisée de l'opérateur $\mathbf{1} + \hat{k}$, en blocs $A_{ij}$ entre les intervalles numéro un $[x_{\min}, x_{\max}^{(0)}]$ et numéro deux $[x_{\max}^{(0)}, x_{\max}^{(n)}]$, on néglige les blocs 12, 21 et 22 de son inverse $A^{-1}$ et on calcule son bloc 11 par l'expression $(A^{-1})_{11} = [A_{11} - A_{12}(A_{22})^{-1}A_{21}]^{-1}$ qui inclut l'effet du couplage non diagonal sous forme d'un hamiltonien effectif du second ordre (il serait incorrect de le négliger complètement car $\langle x | \hat{k} | x' \rangle$ ne tend pas vers zéro lorsque $x, x' \to +\infty$ à $x - x'$ fixé). Comme la fonction $\mathcal{R}(x, x')$ présente près de $x = x' = x_{\text{sing}}$ à gauche des oscillations de plus en plus rapides lorsque $\ell$ augmente, on ne peut pas prendre une valeur fixe du pas $dx$ mais on doit utiliser une méthode itérative de type Romberg : à partir du choix initial $dx = 1/20$, on réduit le pas $dx$ d'un facteur 2 et on extrapole à pas nul linéairement en $(dx)^2$ jusqu'à ce que la valeur extrapolée soit stable au pour cent près (ou que $dx$ descend en dessous de la valeur très faible $2 \times 10^{-3}$). L'élément de matrice dans $|\ell, m_x = 0\rangle$ est calculé comme avant l'équation (52) avec les mêmes astuces de symétrie que dans la note 23 et l'intégrale sur $\phi$ est déduite de l'équation (52).</footnote>
## B.2. Sa somme sur $\ell$ et $\varepsilon$

Il reste à calculer la somme sur toutes les voies ($\ell, \varepsilon$) de l'approximant asymptotique (71). La contribution de type $K_1 K_2$ se traite exactement comme dans l'annexe A :

$$ 
\begin{aligned}
& \sum_{\ell \in \mathbb{N}} \sum_{\varepsilon} \mathcal{J}_{K_1, K_2}^{(\ell, \varepsilon)} \\
&= \int_{\mathbb{R}} dx \int_{-1}^{1} du \frac{[8\pi^2 \operatorname{ch} x \mathcal{D}_{2,2}^2(x, u)]^{-1}}{[1 + \frac{1}{1+\alpha}(1 + e^{-2x} + 2ue^{-x})][1 + \frac{\alpha}{1+\alpha}(1 + e^{2x} + 2ue^{x})]} \\
&= \frac{(1+\alpha)^2}{8\pi^2 \alpha} \left\{ 2\pi \operatorname{arcsin}\left(\frac{\alpha}{(1+2\alpha)(\alpha+2)}\right)^{1/2} \right. \\
& \quad - \operatorname{Re}\left[g_2\left(\frac{2z}{1+z}\right) - g_2\left(\frac{-2z}{1-z}\right) + g_2\left(\frac{2}{1-z}\right) - g_2\left(\frac{2}{1+z}\right) + 2g_2\left(\frac{1}{z}\right) - 2g_2\left(\frac{-1}{z}\right)\right] \\
& \quad \left. - \operatorname{Re}\left[g_2\left(\frac{2z'}{1+z'}\right) - g_2\left(\frac{-2z'}{1-z'}\right) + g_2\left(\frac{2}{1-z'}\right) - g_2\left(\frac{2}{1+z'}\right) + 2g_2\left(\frac{1}{z'}\right) - 2g_2\left(\frac{-1}{z'}\right)\right] \right\}
\end{aligned}
\quad (81) $$

avec $z = \sqrt{\alpha} + i\sqrt{1+\alpha}$, $z' = \sqrt{\alpha^{-1}} + i\sqrt{1+\alpha^{-1}}$ et $g_2$ la fonction dilogarithme. Dans la contribution de type $K_1 K_3$, sommons dans un premier temps sur les moments cinétiques $\ell$ de parité donc de $(-1)^\ell = \eta$ fixés, afin que l'opérateur $\hat{k}$ dans (78) prenne la valeur fixe $\hat{k}_\eta$. Transformons la moyenne quantique dans l'état $|\ell, m_x = 0\rangle$ au numérateur de l'intégrande de (79) en injectant une relation de fermeture sur les états propres de $L_x$ et en utilisant l'expression des éléments de matrice correspondants déduite de la référence [53, équations (7.2–9) page 101, (7.3–15) page 105 et (7.4–7) page 109] :

$$ \langle \ell, m_x = 0 | e^{i\tau L_z/\hbar} | \ell, m_x \rangle e^{i\psi m_x} = \left(\frac{4\pi}{2\ell+1}\right)^{1/2} Y_{\ell}^{m_x}(\tau, \psi) \quad \forall \tau \in [0, \pi], \forall \psi \in \mathbb{R} \quad (82) $$

où $Y_\ell^m$ est l'habituelle harmonique sphérique, ce qui conduit à<fnref n="34" /> <fnref n="35" />

$$ (2\ell + 1) \langle l, m_x = 0 | e^{i\tau_1 L_z/\hbar} e^{i\phi L_x/\hbar} e^{-i\tau'_1 L_z/\hbar} | l, m_x = 0 \rangle = 4\pi \sum_{m_x = -\ell}^{\ell} Y_{\ell}^{m_x}(\tau_1, \phi) [Y_{\ell}^{m_x}(\tau'_1, 0)]^* \quad (83) $$

Il reste à invoquer la relation de fermeture (8.6–10) page 146 de la référence [53] sur les harmoniques sphériques et la propriété de parité spatiale $Y_\ell^m(\pi - \theta, \phi + \pi) = (-1)^\ell Y_\ell^m(\theta, \phi)$ pour obtenir la relation de fermeture à parité de $\ell$ fixée<fnref n="36" /> :

<footnote n="34">³⁴On utilise la relation (82) deux fois, sous sa forme directe avec $(\tau, \psi) = (\tau_1, \phi)$ et sous sa forme conjuguée avec $(\tau, \psi) = (\tau'_1, 0)$. Si $\tau_1$ se trouve dans $[-\pi, 0]$, la relation ne s'applique pas, mais il suffit alors de changer $\tau_1$ en $-\tau_1$, ce qui revient à changer la variable d'intégration $\phi$ en $\phi + \pi$ dans (79) (en effet, $\exp(-i\tau_1 L_z/\hbar) = \exp(i\pi L_x/\hbar)\exp(i\tau_1 L_z/\hbar)\exp(i\pi L_x/\hbar)$) mais ne modifie pas la valeur de l'intégrale. On procède de même si $\tau'_1 \in [-\pi, 0]$. Dans la suite, nous pouvons donc supposer que $\tau_1$ et $\tau'_1$ sont dans l'intervalle $[0, \pi]$.</footnote>

<footnote n="35">³⁵En utilisant le théorème d'addition des harmoniques sphériques, voir [53, équation (8.6–3) page 145], on montre que le second membre de (83) s'écrit aussi $(2\ell + 1)P_\ell(\cos\delta)$ où $P_\ell$ est un polynôme de Legendre et $\delta$ est l'angle entre les vecteurs unitaires de coordonnées polaires $(\tau_1, \phi)$ et $(\tau'_1, 0)$; son cosinus $\cos\delta = \cos\tau_1 \cos\tau'_1 + \sin\tau_1 \sin\tau'_1 \cos\phi$ apparaît justement au dénominateur de l'intégrande de (79).</footnote>

<footnote n="36">³⁶Dans les distributions de Dirac $\delta(\phi - \phi_0)$, $\phi$ a un sens modulo $2\pi$; on peut donc remplacer l'intervalle d'intégration arbitraire $[0, 2\pi]$ de l'équation (79) par l'intervalle de longueur $2\pi$ centré sur $\phi_0$.</footnote>
$$
\begin{align}
& \sum_{\ell \ge 0} \sum_{(-1)^\ell = \eta} \sum_{m_x = -\ell}^{\ell} Y_\ell^{m_x}(\tau_1, \phi) [Y_\ell^{m_x}(\tau'_1, 0)]^* \nonumber \\
&= \frac{1}{2} [\delta(\cos \tau_1 - \cos \tau'_1) \delta(\phi) + \eta \delta(\cos \tau_1 + \cos \tau'_1) \delta(\phi - \pi)] \nonumber \\
&= \frac{|z_1|}{4} \left[ \sqrt{2\alpha} e^{-2x} \delta(x - x') \delta(\phi) + \eta e^{-x} \delta(u - u_0) \delta(\phi - \pi) \right] \tag{84} \\
& \text{où } u_0 \equiv \frac{2\alpha + e^{2x} - e^{2x'}}{2\sqrt{2\alpha}e^x} \nonumber
\end{align}
$$

où nous avons remplacé dans le troisième membre cos τ₁ ± cos τ'₁ par sa valeur en nous souvenant
que |z₁| = |z'₁| et en utilisant (74). L'intégration sur φ est immédiate compte tenu de la note 36.
Pour intégrer sur u, il faut savoir si la racine u₀ se trouve dans l'intervalle [u_min, u_max]. À cette fin,
on partage le support de la fonction de Heaviside dans (79) en quatre zones distinctes :

$$
\exp(x) < \exp(x') < \exp(x) + \sqrt{2\alpha} \quad \text{et} \quad \exp(x) + \exp(x') > \sqrt{2\alpha}, \qquad (i)
$$

$$
\exp(x') < \exp(x) < \exp(x') + \sqrt{2\alpha} \quad \text{et} \quad \exp(x) + \exp(x') > \sqrt{2\alpha}, \qquad (ii)
$$

$$
\exp(x) < \exp(x') < \sqrt{2\alpha} - \exp(x), \qquad (iii)
$$

$$
\exp(x') < \exp(x) < \sqrt{2\alpha} - \exp(x'). \quad (iv)
$$

En appelant $u_{\min}^{\exp}$ et $u_{\max}^{\exp}$ les expressions dépendant de $x$ dans les définitions (75), on trouve qu'on a systématiquement $u_{\min} = -1 < u_0 < u_{\max} = u_{\max}^{\exp}$ dans la zone (i), $u_{\min} = u_{\min}^{\exp} < u_0 < u_{\max} = 1$ dans la zone (ii), $u_{\min} = -1 < u_{\max} = 1 < u_0$ dans la zone (iii), $u_{\min} = u_{\min}^{\exp} < u_{\max} = u_{\max}^{\exp} < u_0$ dans la zone (iv). Autrement dit, l'intégrale de $\delta(u-u_0)$ sur $u$ vaut toujours un dans les deux premières zones et zéro dans les deux dernières. Nous en déduisons la somme de la quantité $\mathcal{R}(x, x')$ sur tous les $\ell$ de parité fixée :

$$
S_{\eta}(x, x') = \frac{1}{2} \left[ \begin{aligned}
& \frac{1}{2} \alpha e^{-2x} \ln \frac{e^{2x} + \sqrt{2/\alpha} e^x + \alpha^{-1} + 1/2}{e^{2x} - \sqrt{2/\alpha} e^x + \alpha^{-1} + 1/2} \delta(x - x') \\
& + \eta \frac{Y(\sqrt{2\alpha} - |e^x - e^{x'}|) Y(e^x + e^{x'} - \sqrt{2\alpha})}{e^{2x} + e^{2x'} + 1}
\end{aligned} \right] \quad (85)
$$

Il reste à sommer sur $\eta = \pm 1$ pour arriver au résultat cherché :

$$
\sum_{\ell=0}^{+\infty} \mathcal{J}_{K_1, K_3}^{(\ell)} = \frac{(1+\alpha)^3}{8\pi\alpha^2} \sum_{\eta=\pm} \int_{\mathbb{R}^2} dx dx' \frac{S_{\eta}(x, x') \langle x' | \mathbf{1} - (\mathbf{1} + \hat{k}_{\eta})^{-1} | x \rangle}{(1+e^{-2x})^{1/4} (1+e^{-2x'})^{1/4}} \quad (86)
$$

Répétons-le, l'opérateur $\hat{k}_{\eta}$ se déduit de l'équation (77) en remplaçant au second membre $(-1)^{\ell+1}$ par $(-\eta)$; l'inversion numérique des opérateurs $\mathbf{1}+\hat{k}_{\eta}$ et l'intégration sur $x$ et $x'$ sont faites avec les mêmes techniques et astuces que dans la note 33 (en revanche, il n'y a plus d'intégration à faire sur $u$). Enfin, comme le montre l'équation (80), $\sum_{\ell \in \mathbb{N}} \mathcal{J}_{K_2, K_3}^{(\ell)}$ se déduit de l'expression (86) en changeant partout $\alpha$ en $1/\alpha$ (y compris dans $S_{\eta}(x, x')$).

**Annexe C.** *L'opérateur* $M_{2,2}^{(\ell,\varepsilon)}$*(iS) *dans la formulation de la référence* [52]

Pour obtenir l'opérateur $M_{2,2}^{(\ell,\varepsilon)}$ (iS) à la base de la conjecture (4), (5)) sur le coefficient d'amas $B_{2,2}(\tilde{\omega})$ du système piégé, il faut d'abord écrire un ansatz de Faddeev pour un état propre du problème unitaire à 2 + 2 corps $\uparrow\uparrow\downarrow$ d'énergie nulle et de quantité de mouvement nulle dans l'espace libre. Cet ansatz s'exprime en termes d'une fonction inconnue D de deux vecteurs
d'onde. La prise en compte des conditions de contact de Wigner-Bethe-Peierls à deux corps
↑↓ (comme dans la note 1) conduit à une équation intégrale sur cette fonction. Ensuite, on
utilise l'invariance par rotation pour projeter l'équation sur le sous-espace de moment cinétique
ℓ et de parité ε comme la référence [54] l'explique en détail. Enfin, on utilise l'invariance
d'échelle du problème unitaire (en l'absence d'effet Efimov à trois corps) pour choisir une
fonction D d'exposant d'échelle s bien défini. L'équation intégrale se résume alors à la condition
detM<sub>2,2</sub><sup>(ℓ,ε)</sup> (s) = 0 où M<sub>2,2</sub><sup>(ℓ,ε)</sup> (s) est un opérateur à noyau, qu'il reste à prolonger sur l'axe imaginaire
pur s = iS pour évaluer l'expression (4).

L'expression ((34), (35), (36)) de $M_{2,2}^{(\ell,\epsilon)}$ (iS) correspond au choix de fonction inconnue $D(\mathbf{k}_2, \mathbf{k}_4)$ où $\mathbf{k}_2$ et $\mathbf{k}_4$ sont les vecteurs d'onde de deux fermions de spins opposés, comme dans la référence [40] (les vecteurs d'onde $\mathbf{k}_1$ et $\mathbf{k}_3$ disparaissent dans la prise de limite traduisant la condition de contact). L'équation intégrale de départ est donnée par l'équation (13) de la référence [40]. Un autre choix est fait dans la référence [52], correspondant au changement de fonction

$$
D(\mathbf{k}_2, \mathbf{k}_4) = F \left( \mathbf{u} \equiv -(\mathbf{k}_2 + \mathbf{k}_4), \mathbf{v} \equiv \frac{m_\uparrow \mathbf{k}_2 - m_\downarrow \mathbf{k}_4}{m_\uparrow + m_\downarrow} \right) \quad (87)
$$

et ayant l'avantage de fournir une expression beaucoup plus simple du noyau intégral singulier $K_3$, c'est-à-dire de la troisième contribution dans (36) (mais pas des noyaux $K_1$ et $K_2$, ce qui justifie en définitive de garder (36) dans le calcul numérique de (4)). Au signe près, les nouvelles variables sont simplement le vecteur d'onde du centre de masse et du mouvement relatif des particules 2 et 4. Dans la suite, nous utilisons la notation $\alpha = m_\uparrow / m_\downarrow$ pour abréger.

Déroulons les étapes précédemment énoncées. L'équation intégrale s'écrit dans le paramé-
trage (87) [52]

$$
\frac{\kappa}{4\pi} F(\mathbf{u}, \mathbf{v}) + \int \frac{d^3 u'}{(2\pi)^3} \left[ \frac{F(\mathbf{u}', \mathbf{v}_{14})}{\kappa^2 + \sigma_{14}^2} + \frac{F(\mathbf{u}', \mathbf{v}_{23})}{\kappa^2 + \sigma_{23}^2} \right] - \int \frac{d^3 v'}{(2\pi)^3} \frac{F(-\mathbf{u}, -\mathbf{v}')}{\kappa^2 + v'^2} = 0 \quad (88)
$$

avec les notations tirées de la référence [52]<sup>37</sup>,

$$
\begin{equation}
\begin{aligned}
\kappa &= \left( v^2 + \frac{2\alpha u^2}{(1+\alpha)^2} \right)^{1/2} ; \quad \sigma_{14} = \mathbf{v} + \frac{1-\alpha}{1+\alpha} \mathbf{u} + \mathbf{u}' ; \quad \mathbf{v}_{14} = \mathbf{v} - \frac{\alpha}{1+\alpha} \mathbf{u} + \frac{\alpha}{1+\alpha} \mathbf{u}' ; \\
\sigma_{23} &= \mathbf{v} + \frac{1-\alpha}{1+\alpha} \mathbf{u} - \mathbf{u}' ; \quad \mathbf{v}_{23} = \mathbf{v} + \frac{\mathbf{u}}{1+\alpha} - \frac{\mathbf{u}'}{1+\alpha}
\end{aligned}
\tag{89}
\end{equation}
$$

en corrigeant ce qui nous semble être une erreur de signe dans le coefficient du vecteur **k** (ici appelé **u'**) dans la définition de **v**<sub>23</sub> (troisième ligne et première colonne de la table III de cette référence) et dans l'expression de *E*<sub>coll</sub> (notation non introduite ici) juste avant l'équation (141) de cette référence. Comme dans la référence [40], pour rendre plus simples les transformations à venir, nous adoptons une formulation variationnelle de l'équation (88), δ*E* / δ*F*<sup>*</sup>(**u**, **v**) = 0 avec la fonctionnelle de *F* et *F*<sup>*</sup> qui suit :

$$
\begin{equation}
\begin{aligned}
\mathcal{E} \equiv \int & d^3 u d^3 v \frac{\kappa}{4\pi} F^*(\mathbf{u}, \mathbf{v}) F(\mathbf{u}, \mathbf{v}) \\
& + \int \frac{d^3 u d^3 v d^3 u' d^3 v'}{(2\pi)^3} F^*(\mathbf{u}, \mathbf{v}) F(\mathbf{u}', \mathbf{v}') \left[ \frac{\delta(\mathbf{v}' - \mathbf{v}_{14})}{\kappa^2 + \sigma_{14}^2} + \frac{\delta(\mathbf{v}' - \mathbf{v}_{23})}{\kappa^2 + \sigma_{23}^2} - \varepsilon \frac{\delta(\mathbf{u} - \mathbf{u}')}{\kappa^2 + v'^2} \right]
\end{aligned}
\tag{90}
\end{equation}
$$

où la parité ε = ±1 de la solution a été introduite pour faire disparaître le signe – devant **u** et **v**'
dans F(-**u**, -**v**'). Tenons maintenant compte de l'invariance par rotation, en nous restreignant
au sous-espace de moment cinétique total ℓ ∈ N, avec une composante de moment cinétique
nulle selon l'axe de quantification Oz. La solution F(**u**, **v**) s'exprime alors en termes de 2ℓ + 1
fonctions ∫<sub>m<sub>z</sub></sub><sup>(ℓ)</sup>(**u**, **v**, **w**) (-ℓ ≤ m<sub>z</sub> ≤ ℓ) de trois variables réelles seulement, les modules **u** et **v** des

<footnote n="37"><sup>37</sup>La référence [52] range les fermions dans l'ordre ↑↓↑↓. Nous renumérotons les particules en conséquence.</footnote>
deux vecteurs et l'angle non orienté $\theta = \langle \mathbf{u}, \mathbf{v} \rangle \in [0, \pi]$ entre eux ou, ce qui revient au même, son cosinus $w = \cos\theta$; fixer la parité à $\varepsilon$ revient à imposer $(-1)^{m_z} = \varepsilon$ donc à découpler le problème en $\ell$ et en $\ell + 1$ fonctions inconnues, ce que nous indiquons par un exposant $\varepsilon$ sur le signe somme ci-dessous. Notre ansatz pour $F$ est celui de [40, équation (14)]. Nous l'insérons dans la fonctionnelle (90) et intégrons sur les variables autres que les arguments des fonctions $f_{m_z}^{(\ell)}$ de la même façon que dans cette référence. Par exemple, appliquons [40, équation (45)] à la dernière contribution de la fonctionnelle $\mathcal{E}$, celle comportant $\varepsilon$ en facteur, qui donne naissance au noyau singulier $K_3$. On choisit d'abord un trièdre direct de référence $\mathcal{T}$ d'axe polaire $\mathbf{u}$ et un autre $\mathcal{T}'$ d'axe polaire $\mathbf{u}'$. L'intégration sur $\mathbf{v}$ (ou sur $\mathbf{v}'$) est effectuée dans les coordonnées sphériques associées à $\mathcal{T}$ (ou à $\mathcal{T}'$), le cosinus de l'angle polaire correspondant étant justement $w$ (ou $w'$). Intégrer sur les autres variables revient alors à prendre la moyenne sur les orientations de $\mathcal{T}$ et $\mathcal{T}'$ ce qui peut se faire en fixant la direction du vecteur $\mathbf{u}$ selon la direction de $\mathbf{u}'$ par une intégration dans le groupe SO(3) sur la rotation $\mathcal{R}$ faisant passer de $\mathcal{T}$ à $\mathcal{T}'$ avec une mesure invariante, explicite $d\mathcal{R} = da(\sin b)dbdc/8\pi^2$ dans le paramétrage d'Euler $\mathcal{R} = \mathcal{R}_Z(a)\mathcal{R}_Y(b)\mathcal{R}_Z(c)$, où les angles $a$ et $c$ décrivent un intervalle de longueur $2\pi$, l'angle $b$ décrit $[0, \pi]$ et le repère cartésien direct OXYZ est d'orientation quelconque par rapport à celui de référence Oxyz (voir [53, section 8.2]); il en résulte que

$$ 
\begin{aligned}
& \int \frac{d^3 u d^3 v d^3 u' d^3 v'}{(2\pi)^3} F^*(\mathbf{u}, \mathbf{v}) F(\mathbf{u}', \mathbf{v}') \frac{\delta(\mathbf{u}-\mathbf{u}')}{\kappa^2 + v'^2} \\
&= \sum_{m_z, m'_z}^{\varepsilon} 2 \int_0^{+\infty} du u^2 dv v^2 du' u'^2 dv' v'^2 \int_{-1}^1 dw dw' f_{m_z}^{(l)*}(u, v, w) f_{m'_z}^{(l)}(u', v', w') \\
&\qquad \times \int_{\text{SO}(3)} d\mathcal{R} \left( (\ell, m_z | R | \ell, m'_z) \right)^* \frac{\delta(\mathbf{u}\mathbf{e}_x - \mathbf{u}'\mathcal{R}\mathbf{e}_x)}{v^2 + v'^2 + \frac{2\alpha u^2}{(1+\alpha)^2}}
\end{aligned} 
\quad (91) $$

où l'opérateur $\mathcal{R}$ représente la rotation $\mathcal{R}$ dans l'espace des états d'une particule quantique. Le choix des axes $OZ = Oz$ et $OY = Ox$ conduit à $(\sin b)\delta(\mathbf{u}\mathbf{e}_x - \mathbf{u}'\mathcal{R}\mathbf{e}_x) = (\sin b)\delta(\mathbf{u}\mathcal{R}_z(-a)\mathbf{e}_x - \mathbf{u}'\mathcal{R}_x(b)\mathcal{R}_z(c)\mathbf{e}_x) = \delta(\mathbf{u}\cos a - \mathbf{u}'\cos c)\delta(-\mathbf{u}\sin a - \mathbf{u}'\sin c\cos b)\delta(-\mathbf{u}'\sin c) = [\delta(\mathbf{u} - \mathbf{u}')/(\mathbf{u}\mathbf{u}')]\delta(a)\delta(c) + \delta(a - \pi)\delta(c - \pi)]$ où nous avons successivement utilisé l'invariance par rotation de la distribution de Dirac tridimensionelle et décomposé son action en distributions de Dirac selon $Ox$, $Oy$ et $Oz$. L'intégration dans SO(3) se réduit simplement à la ligne $\mathcal{R} = \mathcal{R}_x(b)$ si $a = c = 0$ et à la ligne $\mathcal{R} = \mathcal{R}_x(-b)$ si $a = c = \pi$, ce dont on tient compte par une intégration sur $b$ étendue à $[-\pi, \pi]$, ce qui fait apparaître un projecteur sur l'état de moment cinétique nul selon $Ox$, $\int_{-\pi}^{\pi} \frac{db}{2\pi} \langle \ell, m_z | e^{-ibL_x/\hbar} | \ell, m'_z \rangle = \langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m'_z \rangle$. Finalement

$$ 
\begin{aligned}
\mathcal{E} = & \sum_{m_z}^{\varepsilon} \int_0^{+\infty} du u^2 dv v^2 \int_{-1}^{1} dw f_{m_z}^{(l)*}(u, v, w) f_{m_z}^{(l)}(u, v, w) \frac{1}{2} \left( v^2 + \frac{2\alpha u^2}{(1+\alpha)^2} \right)^{1/2} \\
& + \sum_{m_z, m'_z}^{\varepsilon} \int_0^{+\infty} du u^2 dv v^2 du' u'^2 dv' v'^2 \int_{-1}^{1} dw dw' f_{m_z}^{(l)*}(u, v, w) f_{m'_z}^{(l)}(u', v', w') \\
& \qquad \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \left[ \frac{\delta(|Z_1| - |Z'_1|)(|Z_1| - |Z'_1|)^{-1} e^{iT_1 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-iT'_1 m'_z}}{|Z_2|^2 + |Z'_2|^2 + \frac{2\alpha u u'}{1+\alpha} (\cos T_1 \cos T'_1 + \sin T_1 \sin T'_1 \cos \phi)} \right. \\
& \qquad \left. + \frac{\delta(|Z_2| - |Z'_2|)(|Z_2| - |Z'_2|)^{-1} e^{iT_2 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-iT'_2 m'_z}}{|Z_1|^2 + |Z'_1|^2 + \frac{2\alpha u u'}{1+\alpha} (\cos T_2 \cos T'_2 + \sin T_2 \sin T'_2 \cos \phi)} \right] \\
& + \frac{(-\varepsilon)\delta(\mathbf{u}-\mathbf{u}') \langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m'_z \rangle}{uu'} \frac{v^2 + v'^2 + \frac{2\alpha u^2}{(1+\alpha)^2}}
\end{aligned} 
\quad (92) $$
où nous avons introduit les nombres complexes $Z_1 \equiv \alpha u/(1+\alpha) - v \exp(i\theta) = |Z_1|\exp(iT_1)$,
$Z_2 = u/(1+\alpha) + v \exp(i\theta) = |Z_2|\exp(iT_2)$ et leurs homologues pour les variables primées, $Z'_1 \equiv$
$\alpha u'(1+\alpha) - v' \exp(i\theta') = |Z'_1|\exp(iT'_1)$, $Z'_2 = u'(1+\alpha) + v' \exp(i\theta') = |Z'_2|\exp(iT'_2)$. Enfin, tirons
parti de l'invariance d'échelle au moyen de l'ansatz

$$
f_{mz}^{(l)}(u, v, w) = \frac{(\operatorname{ch} x)^{s+3/2} (1 + e^{-2x})^{-s/2} \Phi_{mz}^{(l)}(x, w)}{(u^2 + v^2/\zeta^2)^{(s+7/2)/2}} \quad \text{où } v = \zeta e^x u \text{ et } \zeta = \frac{\sqrt{2\alpha}}{1+\alpha} \quad (93)
$$

choisi astucieusement pour que la partie diagonale de la fonctionnelle (première contribu-
tion dans (92)) devienne scalaire (indépendante de toute variable) et pour que la partie sin-
gulière (dernière contribution) ne dépende pas du rapport de masse ni même de l'exposant
d'échelle s. L'intégration sur u fait sortir comme dans la référence [40] un facteur constant in-
fini, ici (ζ^4/16) ∫_0^+∞ du/u, pour donner la fonctionnelle régularisée

$$
\begin{equation}
\begin{split}
\bar{\mathcal{E}} = \sum_{m_z}^{\mathcal{E}} & \int_{-\infty}^{+\infty} dx \int_{-1}^{1} dw \Phi_{m_z}^{(l)*}(x, w) \Phi_{m_z}^{(l)}(x, w) \\
& + \frac{2}{\zeta^3} \sum_{m_z, m'_z}^{\mathcal{E}} \int_{-\infty}^{+\infty} dx dx' \int_{-1}^{1} dw dw' \frac{e^{x+x'} \Phi_{m_z}^{(l)*}(x, w) \Phi_{m'_z}^{(l)}(x', w')}{(1+e^{-2x})^{1/4} (1+e^{-2x'})^{1/4}} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \\
& \qquad \left[ \frac{((|z'_1|/|z_1|)^s (|z_1| |z'_1|)^{-1/2} e^{i\tau_1 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-i\tau'_1 m'_z}}{|z'_1|^2 |z_2|^2 + |z_1|^2 |z'_2|^2 + (1+\beta) (\operatorname{Re} z_1 \operatorname{Re} z'_1 + \operatorname{Im} z_1 \operatorname{Im} z'_1 \cos \phi)} \right. \\
& \qquad \left. + \frac{((|z'_2|/|z_2|)^s (|z_2| |z'_2|)^{-1/2} e^{i\tau_2 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-i\tau'_2 m'_z}}{|z'_2|^2 |z_1|^2 + |z_2|^2 |z'_1|^2 + (1+\alpha) (\operatorname{Re} z_2 \operatorname{Re} z'_2 + \operatorname{Im} z_2 \operatorname{Im} z'_2 \cos \phi)} \right. \\
& \qquad \left. + (-\epsilon) \zeta^3 \frac{\langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m'_z \rangle}{e^{2x} + e^{2x'} + 1} \right]
\end{split}
\tag{94}
\end{equation}
$$

où l'on a posé $\beta = 1/\alpha$, $z_1 \equiv (\alpha/2)^{1/2} - \exp(x+i\theta) = |z_1|\exp(i\tau_1)$, $z_2 \equiv (\beta/2)^{1/2} + \exp(x+i\theta) = |z_2|\exp(i\tau_2)$ et les relations équivalentes pour les variables primées, $z'_1 \equiv (\alpha/2)^{1/2} - \exp(x'+i\theta') = |z'_1|\exp(i\tau'_1)$, $z'_2 \equiv (\beta/2)^{1/2} + \exp(x'+i\theta') = |z'_2|\exp(i\tau'_2)$ et où l'on rappelle que $\theta = \arccos w$ et $\theta' = \arccos w'$. Une simple dérivation fonctionnelle de (94) par rapport à $\Phi_{mz}^{(\ell)*}(x, w)$, la prise d'un exposant d'échelle imaginaire pur et le retour à la notation $u = \cos\theta$ du corps de l'article donnent la nouvelle forme de l'opérateur $M_{2,2}^{(\ell,\epsilon)}(iS)$, repérée par un accent tchèque pour éviter toute confusion avec ((34), (35), (36)) :

$$
\begin{equation}
\begin{split}
& \langle x, u | \langle \ell, m_z | \check{M}_{2,2}^{(\ell,\epsilon)} (iS) | f \rangle \\
& \qquad = f_{m_z}(x, u) + \int_{-\infty}^{+\infty} dx' \int_{-1}^{1} du' \sum_{m'_z \mid (-1)^{m'_z} = \epsilon} \check{K}_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) f_{m'_z}(x', u')
\end{split}
\tag{95}
\end{equation}
$$

avec un noyau intégral séparé en trois contributions $\check{K} = \check{K}_1 + \check{K}_2 + \check{K}_3$ écrites ligne par ligne :

$$
\begin{align}
& \check{K}_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) \nonumber \\
& = \frac{2\zeta^{-3}e^{x+x'}(|z_1||z'_1|)^{-1/2}}{(1+e^{-2x})^{1/4}(1+e^{-2x'})^{1/4}} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{|z_1|^{-iS} e^{i\tau_1 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-i\tau'_1 m'_z} |z'_1|^{iS}}{|z'_1|^2|z_2|^2 + |z_1|^2|z'_2|^2 + (1+\beta)(\operatorname{Re}z_1 \operatorname{Re}z'_1 + \operatorname{Im}z_1 \operatorname{Im}z'_1 \cos\phi)} \nonumber \\
& \quad + \frac{2\zeta^{-3}e^{x+x'}(|z_2||z'_2|)^{-1/2}}{(1+e^{-2x})^{1/4}(1+e^{-2x'})^{1/4}} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{|z_2|^{-iS} e^{i\tau_2 m_z} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m'_z \rangle e^{-i\tau'_2 m'_z} |z'_2|^{iS}}{|z'_2|^2|z_1|^2 + |z_2|^2|z'_1|^2 + (1+\alpha)(\operatorname{Re}z_2 \operatorname{Re}z'_2 + \operatorname{Im}z_2 \operatorname{Im}z'_2 \cos\phi)} \nonumber \\
& \quad + \frac{(-1)^{\ell+1}}{\pi} \frac{e^{x+x'}}{e^{2x} + e^{2x'} + 1} \frac{\langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m'_z \rangle}{(1+e^{-2x})^{1/4}(1+e^{-2x'})^{1/4}} \tag{96}
\end{align}
$$
où l'on s'est souvenu que $\langle l, m_z | l, m_x = 0 \rangle = 0$ si $\varepsilon \neq (-1)^l$. Pour vérifier que l'équation (96) est écrite dans le même ordre que l'équation (36), c'est-à-dire que le noyau $\tilde{K}_i$ est précisément le noyau $K_i$ écrit dans la formulation de la référence [52], nous montrons numériquement que $\int_{\mathbb{R}} dS \text{Tr}[\mathcal{D}_{2,2}^{-1}K_1(iS)\mathcal{D}_{2,2}^{-1}K_3(iS)] = \int_{\mathbb{R}} dS \text{Tr}[\tilde{K}_1(iS)\tilde{K}_3]$.

## English version

# 1. Introduction to the problem and main results

## 1.1. The system

Our object of study is a three-dimensional gas of non-relativistic neutral fermions with two spin states $\uparrow$ and $\downarrow$ in the regime of maximum interaction allowed in the gas phase, i.e. in the so-called unitary limit [1, 2]: there is no interaction between fermions in the same spin state and a binary zero-range interaction, exclusively in the s-wave and of infinite scattering length between fermions of different spin states $\uparrow$ and $\downarrow$.<fnref n="1" /> We consider the general case where the fermion masses $m_\uparrow$ and $m_\downarrow$ can be different for the two spin states, which is the originality of our work. This system is realizable in laboratory with trapped cold atomic gases. Indeed, its main properties at equilibrium have been determined experimentally for equal masses $m_\uparrow = m_\downarrow$: its superfluidity at low temperature has been observed in the unpolarized case (with equal numbers of $\uparrow$ and $\downarrow$ fermions) [3, 4] and its equation of state has been measured at any temperature and polarization [5–7]. The case $m_\uparrow \neq m_\downarrow$ remains to be explored; a natural way to obtain it is to use a mixture of two fully polarised fermionic atomic species, as has been done in references [8–10].

## 1.2. Cluster and virial expansions

Given the strength of the interactions, which provides no obvious small parameter, there are few theoretical tools for quantitatively reliable predictions on the unitary gas that can be compared to measurements. One of them is the diagrammatic Monte Carlo simulation of the many-body problem on a computer [11, 12]. Another is the cluster or virial expansion [13] of the pressure $P$ of the spatially homogeneous infinite gas into powers of the fugacities $z_\sigma = \exp(\beta\mu_\sigma)$ or of the phase space densities $\rho_\sigma\lambda_\sigma^3$ respectively, in the strongly non-degenerate regime where they tend to zero, with $\rho_\sigma$ the density of the $\sigma = \uparrow, \downarrow$ component in real space, $\mu_\sigma$ its chemical potential, $\lambda_\sigma = (2\pi\hbar^2/m_\sigma k_B T)^{1/2}$ its de Broglie thermal wavelength at temperature $T$, and $\beta = 1/k_B T$ [14]:

$$ \frac{P \lambda^3}{k_B T} = \sum_{(n_\uparrow, n_\downarrow) \in \mathbb{N}^{2*}} b_{n_\uparrow, n_\downarrow} z_\uparrow^{n_\uparrow} z_\downarrow^{n_\downarrow} = \sum_{(n_\uparrow, n_\downarrow) \in \mathbb{N}^{2*}} c_{n_\uparrow, n_\downarrow} (\rho_\uparrow \lambda_\uparrow^3)^{n_\uparrow} (\rho_\downarrow \lambda_\downarrow^3)^{n_\downarrow} \quad (1) $$

$$ \psi(\mathbf{r}_1, \dots, \mathbf{r}_{n_\uparrow+n_\downarrow}) = A_{ij} (a^{-1} - r_{ij}^{-1}) + O(r_{ij}) $$

<footnote n="1">¹ This situation corresponds to the replacement of the interaction potential by the Wigner-Bethe-Peierls contact conditions on the wave function $\psi$ of the system in each spin configuration $|\uparrow\dots\uparrow\downarrow\dots\downarrow\rangle = |\uparrow\rangle^{\otimes n_\uparrow} |\downarrow\rangle^{\otimes n_\downarrow}$: when the distance $r_{ij}$ between particle $i$ in spin state $\uparrow$ and particle $j$ in spin state $\downarrow$ tends to zero at a fixed position $\mathbf{R}_{ij}$ of their center of mass, the positions $\mathbf{r}_k$ of the other particles being fixed at values different from $\mathbf{R}_{ij}$, there exists a constant $A_{ij}$ (function of $\mathbf{R}_{ij}$ and of the $\mathbf{r}_k$'s) such that</footnote>

<footnote>where $a$ is the scattering length, and this for all $i$ and for all $j$. Here, we are in the unitary limit $1/a = 0$.</footnote>
To scale the pressure, it was necessary to introduce a reference de Broglie thermal wavelength $\lambda = (2\pi\hbar^2/\bar{m}k_B T)^{1/2}$ depending on a mean mass $\bar{m}$ to be specified (see equation (12)). One can then try to extrapolate to the non-trivial regime $z_σ \approx 1$ with heuristic recipes such as the Padé approximant [15] or optimized resummation methods [16] which take into account the behavior of the coefficients at large orders. The cluster or virial expansion has the advantage over the Monte Carlo simulation that it is closer to the analytical calculation, since the coefficients of order $n$ can be obtained from the solution of a problem with at most $n$ interacting fermions, i.e. with few bodies: to obtain $b_{n_1,n_1}$, it suffices to determine the canonical partition functions of all systems with $n_σ$ or less fermions in each spin state $σ = \uparrow, \downarrow$.

### 1.3. The unitary limit

The scale invariance of the unitary gas simplifies considerably the calculation of cluster coefficients for $n > 2$ (the order two, given by the Beth–Uhlenbeck formula [17–19], is not debated). Thus, the third-order coefficients are known analytically, even if the scale invariance at the three-body level is broken by the Efimov effect [20,21], in sharp contrast to the model of hard sphere interaction of radius $a$ where the coefficients are known analytically (for bosons) only in the limiting cases $\lambda/a \ll 1$ [22–25] or $\gg 1$ [26–29]. The harmonic regulator method [30–33], consisting in trapping each component of the gas in a fictitious isotropic harmonic potential, $U_σ(\mathbf{r}) = m_σ\omega^2r^2/2$, whose trapping frequency $\omega$ (common to both spin states) is made to tend to zero at the end of the calculations, allows us to take full advantage of the scale invariance since the $n$-body spectrum in the trap is obtained from the discrete set of scale exponents $s_i$ of the zero energy $E = 0$ eigenstates in free space [1,34]. More precisely, we generalize the cluster expansion to the trapped case, replacing the pressure by the grand potential $\Omega$ and taking its ratio to the partition function $Z_1 = 1/[2\text{sh}(\tilde{\omega}/2)]^3$ of a single fermion,<fnref n="2" /> so that

$$ \frac{-\Omega}{k_B T Z_1} = \sum_{(n_1, n_\downarrow) \in \mathbb{N}^{2*}} B_{n_1, n_\downarrow}(\tilde{\omega}) z_1^{n_\downarrow} z_\downarrow^{n_\downarrow} \quad (2) $$

The coefficients of the trapped case, marked by a capital letter to avoid confusion, depend only on the dimensionless ratio $\tilde{\omega} = \hbar\omega/k_B T$ due to the scale invariance of the unitary gas. They are related to those of the homogeneous case by means of the local density approximation, exact in the limit $\omega \to 0$ [32, 33] and giving [21]:

$$ b_{n_1, n_\downarrow} = \frac{(n_1 m_\downarrow + n_\downarrow m_1)^{3/2}}{\bar{m}^{3/2}} B_{n_1, n_\downarrow}(0^+) \quad (3) $$

This method leads to an analytic integral expression of the third order coefficients because the transcendental Efimov function $\Lambda(s)$, whose roots are the scale exponents $s_i$, is known explicitly [35–38]. This allows to express the coefficients as a contour integral around $\mathbb{R}^+$ by means of the residue theorem and then to unfold the contour on the pure imaginary axis by analyticity on $C \setminus \mathbb{R}$ [20]. On the other hand, for non-trivial fourth order coefficients,<fnref n="3" /> the Efimov function is the determinant of operators $M_{3,1}(s)$ [39] or $M_{2,2}(s)$ [40] parametrically dependent on $s$; we know

$$ b_{n,0} = (m_\downarrow / \bar{m})^{3/2} (-1)^{n+1} / n^{5/2} \quad \text{and} \quad b_{0,n} = (m_\uparrow / \bar{m})^{3/2} (-1)^{n+1} / n^{5/2}. $$

<footnote n="2">² Indeed, the numerator in the first side of equation (1) can be seen, in a quantization box of arbitrarily large volume $V$, as the ratio between $PV$, i.e. minus the gas grand potential, and the partition function $V/\lambda^3$ of a fictitious single particle of mass $\bar{m}$.</footnote>

<footnote n="3">³ If one of the $n_σ$ is zero, the fermions of the $n$-body problem are non-interacting and the corresponding $B$ coefficient reduces to that of the ideal gas. We thus find $B_{n,0}(0^+) = B_{0,n}(0^+) = (-1)^{n+1}/n^4$ and, by means of relation (3),</footnote>
how to compute it numerically only on the pure imaginary axis $s = iS$ and its imperfectly known analytical properties in the complex plane do not guarantee a safe use of the residue theorem. Therefore reference [41] could only produce a conjecture, which we briefly recall.

## 1.4. Our 2016 conjecture on the fourth cluster coefficients

As in reference [40], we introduce the integral expression modeled on that of the third order coefficients ($\Lambda(s)$ is replaced by an operator determinant):

$$I_{n_1, n_1}(\bar{\omega}) = \sum_{\ell=0}^{+\infty} \sum_{\epsilon} (\ell + 1/2) \int_{-\infty}^{+\infty} \frac{dS}{2\pi} \frac{\sin(\bar{\omega}S)}{\sinh\bar{\omega}} \frac{d}{dS} \ln \det [M_{n_1, n_1}^{(\ell, \epsilon)}(iS)] \quad (4)$$

where the sum is taken on the internal angular momentum $\ell$ (i.e. after separation of the center of mass) of the four-body eigenstates and their internal parity $\epsilon = \pm 1$ (limited to $\epsilon = 1$ for $\ell = 0$), and where the operator $M_{n_1, n_1}(iS)$ is restricted to the corresponding $(\ell, \epsilon)$ subspace. Then $I_{n_1, n_1}(\bar{\omega})$ gives exactly $B_{n_1, n_1}(\bar{\omega})$ (as is the case at order three) when the asymptotic decoupled objects (independent non-monoatomic packets of fermions strongly correlated by interactions, in terms of which the highly excited eigenstates are expressed) are distinguishable; otherwise, the quantum statistical effect of the indistinguishability of these objects, which contributes to the cluster coefficient even if they do not interact with each other, is missing, for the same reason that the cluster coefficients of the quantum ideal gas differ from those of the classical ideal gas. In the $(n_1, n_1) = (3, 1)$ configuration, called $3+1$ for short from here on, the possible asymptotic objects are a $\uparrow\uparrow\downarrow$ triplon or a $\uparrow\downarrow$ pairon of strongly correlated fermions; there can only be one at a time, which rules out any quantum statistical effect. On the other hand, in the $2+2$ configuration, the fermions can also decouple into two $\uparrow\downarrow$ pairons of correlated fermions, which are indistinguishable bosons; since these bosons do not interact, their partition function is easily calculated [41]. The conjecture of reference [41] is finally written

$$B_{3,1}(\bar{\omega}) = I_{3,1}(\bar{\omega}) \quad ; \quad B_{2,2}(\bar{\omega}) = I_{2,2}(\bar{\omega}) + \frac{1}{32} \frac{1}{\cosh\bar{\omega} \cosh^3(\bar{\omega}/2)} \quad (5)$$

the $1+3$ case being deduced from the $3+1$ case by exchanging the two spin states thus changing the mass ratio $\alpha = m_1/m_1$ to its inverse $1/\alpha$. For $\alpha = 1$, conjecture (5) is in agreement with the quantum Monte Carlo calculation of reference [42] down to the minimal accessible value of $\bar{\omega}$, $\bar{\omega} \approx 1$; it is also in agreement with a recent, more powerful numerical calculation down to values of $\bar{\omega} \ll 1$ [43, 44]. After using relation (3), it is also in agreement with the same powerful numerical calculation performed directly in the spatially homogeneous case of a quantization box [45]. Conjecture (5) is thus confirmed for equal masses.

## 1.5. Content of the study

The idea of the present work is to believe in the validity of conjecture (5) for any mass ratio $\alpha$ and to access the corresponding fourth cluster and virial coefficients, by numerically calculating the determinant of the operators $M_{n_1, n_1}^{(\ell, \epsilon)}(iS)$, and then integrating over $S$ and summing over $\ell$ and $\epsilon$ in expression (4). In reality, we still have to put bounds on the mass ratio, because the method of solving the four-body problem in a harmonic trap, at the basis of expression (5), assumes separability of the internal Schrödinger equation in hyperspherical coordinates; this is true only
if there is no Efimov scale invariance breaking at the three-body level, which constrains us to the intervals

$$
\begin{align}
& \text{for } B_{3,1}: \alpha < \alpha_c^{3\text{body}} \approx 13.60697 \ ; \quad \text{for } B_{1,3}: \alpha > \frac{1}{\alpha_c^{3\text{body}}} \approx 0.0734917 \ ; \nonumber \\
& \text{for } B_{2,2}: \frac{1}{\alpha_c^{3\text{body}}} < \alpha < \alpha_c^{3\text{body}} \tag{6}
\end{align}
$$

where $\alpha_c^{3\text{body}}$, also noted $\alpha_c^{2,1}$, is the threshold of the three-body Efimov effect in the $\uparrow\uparrow\downarrow$ system, which occurs in the three-body internal momentum channel $L=1$ [35, 46]. In using form (5) of the conjecture, it is also important to ensure that there is no four-body Efimov effect;<fnref n="4" /> indeed, it was shown that such an effect occurs only in the $3+1$ or $1+3$ configuration [39, 40], at mass ratio $\alpha_c^{4\text{body}} = \alpha_c^{3,1}$ or its inverse, hence the additional conditions

$$
\text{for } B_{3,1}: \alpha < \alpha_c^{4\text{body}} \approx 13.3842 \ ; \quad \text{for } B_{1,3}: \alpha > \frac{1}{\alpha_c^{4\text{body}}} \approx 0.074715 \quad (7)
$$

Once the cluster expansion is known, the virial coefficients $c_{n_1, n_1}$ are easily deduced, as rational functions of the cluster coefficients $b_{n'_1, n'_1}$ of total order $n'_1 + n'_1$ less than or equal to $n_1 + n_1$:<fnref n="5" />

$$
\begin{align}
c_{1,0} &= b_{1,0} \quad ; & c_{2,0} &= -b_{2,0} \quad ; \nonumber \\
c_{3,0} &= \frac{4b_{2,0}^2}{b_{1,0}} - 2b_{3,0} \quad ; & c_{4,0} &= -3b_{4,0} - \frac{20b_{2,0}^3}{b_{1,0}^2} + \frac{18b_{3,0}b_{2,0}}{b_{1,0}} \quad ; & c_{1,1} &= -b_{1,1} \tag{8}
\end{align}
$$

$$
\rho_{\sigma} \Lambda^3 = \sum_{(n_1, n_1) \in \mathbb{N}^{2*}} n_{\sigma} b_{n_1, n_1} z_1^{n_1} z_1^{n_1}
$$

from the thermodynamic relation $\rho_\sigma = \partial P / \partial \mu_\sigma$, and we adjust the $c_{n_1, n_1}$ order by order to find the expansion in the second side of equation (1). It is simpler, at first, to write the virial expansion in the form of

$$
P \lambda^3 / k_B T = \sum_{(n_1, n_1) \in \mathbb{N}^{2*}} a_{n_1, n_1} (\rho_1 \lambda^3)^{n_1} (\rho_1 \lambda^3)^{n_1}
$$

and then, once the coefficients $a_{n_1, n_1}$ have been calculated, to use the relation

$$
c_{n_1, n_1} = a_{n_1, n_1} / (t_1^{n_1} t_1^{n_1})
$$

where we have put $t_\sigma = (\bar{m}/m_\sigma)^{3/2}$. Note that $b_{1,0}t_1 = b_{0,1}t_1 = 1$, and that $a_{1,0} = a_{0,1} = 1$ whatever the reference mass $\bar{m}$ according to the ideal gas law $P = (\rho_1 + \rho_2)k_B T$. Replacing $b_{n,0}$ and $b_{0,n}$ by their expressions given in Footnote 3, we finally get

$$
\begin{align*}
c_{3,1} &= \left(-\frac{3}{4} + \frac{1}{\sqrt{3}}\right) b_{1,1} - \left(\frac{3}{\sqrt{2}}\right) b_{2,1} - 3b_{3,1} + 3b_{1,1} \left(b_{2,1} + \frac{b_{1,1}}{\sqrt{8}}\right) t_1 - b_{1,1}^3 t_1^2, \\
c_{2,2} &= -3 \left[b_{2,2} + \left(b_{1,2} + b_{2,1}\right)/\sqrt{8} + b_{1,1}/8\right] + \left(9\sqrt{2}/8\right) b_{1,1}^2 (t_1+t_2) + 6b_{1,1} \left(b_{2,1}t_1 + b_{1,2}t_2\right) - 3b_{1,1}^3 t_1 t_2
\end{align*}
$$

(the ideal gas coefficients $c_{4,0}$ and $c_{0,4}$ are given in explicit form later in the main text). We easily derive $b_{1,1}$ from equation (3) knowing that $B_{1,1}(0^+) = 1/2$. On the other hand, $b_{2,1}$ and $b_{1,2}$ are computed numerically from explicit integral expressions of reference [21].

<footnote n="4">⁴ Moreover, our assumption of a strictly scale-invariant interaction becomes difficult to satisfy experimentally if $\alpha$, although on the non-Efimovian side, is too close to the tetramer appearance threshold. All this can be improved. One could complement the $3+1$-body contact condition of the usual zero range interaction model of Footnote 1 by means of a length called “3+1-body parameter” in the $(\ell, \epsilon) = (1, +1)$ [39] channel where the $3+1$-body Efimov effect occurs. In this case, $B_{3,1}(0^+)$ and $b_{3,1}$ would be smooth functions of $\alpha$ at $\alpha = \alpha_c^{4\text{body}}$. Reference [21] proves and implements this for $B_{2,1}(0^+)$ that is in the $2+1$ body problem.</footnote>

<footnote n="5">⁵ To obtain them, we replace in the expansion of the pressure in the third side of equation (1) the densities by their virial expansion</footnote>
$$c_{2,1} = -2b_{2,1} + \frac{4b_{1,1}b_{2,0}}{b_{1,0}} + \frac{b_{1,1}^2}{b_{0,1}}; \qquad (9)$$

$$c_{3,1} = -3b_{3,1} - \frac{b_{1,1}^3}{b_{0,1}^2} - \frac{6b_{1,1}^2 b_{2,0}}{b_{1,0} b_{0,1}} - \frac{24b_{1,1}b_{2,0}^2}{b_{1,0}^2} + \frac{3b_{1,1}b_{2,1}}{b_{0,1}} + \frac{12b_{2,0}b_{2,1}}{b_{1,0}} + \frac{9b_{1,1}b_{3,0}}{b_{1,0}}$$

$$c_{2,2} = -3b_{2,2} - \frac{3b_{1,1}^3}{b_{0,1} b_{1,0}} - 9b_{1,1}^2 \left( \frac{b_{0,2}}{b_{0,1}^2} + \frac{b_{2,0}}{b_{1,0}^2} \right) + 6b_{1,1} \left( \frac{b_{1,2}}{b_{0,1}} - \frac{2b_{0,2}b_{2,0}}{b_{0,1}b_{1,0}} + \frac{b_{2,1}}{b_{1,0}} \right) + \frac{6b_{0,2}b_{2,1}}{b_{0,1}} + \frac{6b_{1,2}b_{2,0}}{b_{1,0}} \quad (10)$$

plus the equations obtained by exchanging the two spin states ↑ and ↓. Our expressions of $c_{n,0}$ are in agreement with [13, equation (10.33)], except that they extend it to a ratio $m_↑/\bar{m}$ different from unity through the coefficient $b_{1,0} = \lambda^3/\lambda_↑^3$.

Before presenting our results, we need to choose the reference mass $\bar{m}$. To do so, we rely on the particular form of the cluster expansion commonly used in the literature for equal fugacities, which puts the number of spin states of the fermions as a factor of the series,

$$\frac{P\lambda^3}{k_B T} \stackrel{z_\sigma=z}{\equiv} 2 \sum_{n=1}^{+\infty} b_n z^n \quad \text{with} \quad b_n = \frac{1}{2} \sum_{n_\uparrow=0}^{n} b_{n_\uparrow, n-n_\uparrow} \qquad (11)$$

see in particular references [5,7] where $b_4$ is measured for equal masses.<fnref n="6" /> Pulling out such a factor makes sense if the coefficient $b_1$ takes the very simple value equal to one. From Footnote 3, we thus derive the natural choice<fnref n="7" />

$$b_1 = 1 \quad \Rightarrow \quad \bar{m}^{3/2} = \frac{1}{2} (m_{\uparrow}^{3/2} + m_{\downarrow}^{3/2}) \quad \text{hence} \quad \frac{1}{\lambda^3} = \frac{1}{2} \left( \frac{1}{\lambda_{\uparrow}^3} + \frac{1}{\lambda_{\downarrow}^3} \right) \qquad (12)$$

This choice also makes the virial expansion very simple in the case of equal phase space densities in the two spin states: one then has the nice expression $\rho_\sigma \lambda_\sigma^3 = \rho \lambda^3/2$ for the joint value, where $\rho = \rho_↑ + \rho_↓$ is the total density, and one sets as in equation (11):

$$\frac{P\lambda^3}{k_B T} \stackrel{\rho_\sigma \lambda_\sigma^3 = \rho \lambda^{3/2}}{\equiv} 2 \sum_{n=1}^{+\infty} c_n (\rho \lambda^{3/2})^n \quad \text{with} \quad c_n = \frac{1}{2} \sum_{n_\uparrow=0}^{n} c_{n_\uparrow, n-n_\uparrow} \qquad (13)$$

## 1.6. The results

We plot the unitary gas non-trivial fourth cluster coefficients (we do not plot here those of the ideal gas) as functions of the mass ratio $\alpha$ in Figure 1a (the values $b_{4,0} = -(m_↑/\bar{m})^{3/2}/32$ and $b_{0,4} = -(m_↓/\bar{m})^{3/2}/32$ taken from Footnote 3 are thus not plotted but of course contribute to $b_4$). The logarithmic scale used on the x-axis highlights the $\alpha \leftrightarrow 1/\alpha$ symmetry. For a mass ratio

<footnote n="6">⁶ The values of $b_4$ measured at ENS and MIT are in agreement with each other but within a factor ≈ 2 and respectively within ≈ 2σ and ≈ 3.5σ of the now accepted value [45], where σ is the experimental uncertainty. In the case of ENS, the discrepancy is due to the fact that all usable pressure data are of fugacity $z > 1$ and their naive extrapolation to $z = 0$ fails. Indeed, series (11) converges slowly and cannot be truncated to a good approximation at $n = 4$ even for a value as small as $z = 0.22$, see the augmented version [47] of publication [41] which draws this conclusion from the high precision equation of state obtained by diagrammatic Monte Carlo [11, 16].</footnote>

<footnote n="7">⁷ In the context of Footnote 2, this choice amounts to taking as a reference one-body partition function $\bar{Z}_1$ the arithmetic mean of $Z_{1↑}$ and $Z_{1↓}$. It reduces indeed to $\bar{m} = m$ in the case of equal masses $m_↑ = m_↓ = m$. It ensures that the factor relating the trapped case to the homogeneous case in equation (3) is a uniformly bounded function of the masses (bounded from above by $2(n_↑^3 + n_↓^3)^{1/2}$ according to Hölder's inequality of parameters $p = 3/2$ and $q = 3$). These conditions are not satisfied if we take for $\bar{m}$ the reduced mass of two ↑ and ↓ fermions as in reference [21].</footnote>
$\alpha = 1$, we find the value $b_4 = 0.030(1)$ conjectured by Reference [41] and confirmed numerically by Reference [45]. As we move away from this point, the behavior of $b_4$ is first imposed by $b_{2,2}$ which causes it to change sign (the black curve and the green curve are almost parallel), before $b_{3,1}$ or $b_{1,3}$ prevails at high mass ratios $\alpha$ or $1/\alpha$ and causes $b_4$ to rise to large and positive values. At the $\alpha_c^{4\text{body}}$ or $1/\alpha_c^{4\text{body}}$ thresholds of the 3+1- or 1+3-body Efimov effect, marked by dotted vertical lines, $b_4$ has a finite limit but an infinite derivative like $b_{3,1}$ or $b_{1,3}$, see Section 2 below and Footnote 4. On the other hand, $b_{2,2}$ remains a smooth function, and would only exhibit an infinite derivative at the thresholds of the three-body Efimov effect, see Section 3. In Figure 1b, we see that the non-trivial fourth virial coefficients (we do not plot here the trivial ones $c_{4,0} = A_4(m_↓/\bar{m})^{3/2}$ and $c_{0,4} = A_4(m_↓/\bar{m})^{3/2}$, with $A_4 = (18+15\sqrt{2}-16\sqrt{6})/192 \approx 1.11 \times 10^{-4}$) have a similar structure even though they are of constant sign. We also notice that coefficient $c_4$ depends little on the four-body correlations induced by the interactions: an approximation of $c_4$ neglecting the contributions of the cluster coefficients $b_{3,1}$, $b_{2,2}$ and $b_{1,3}$, plotted in dashed line in Figure 1b, is everywhere close or very close to $c_4$. The fairly recent possibility to prepare spatially homogeneous cold atomic gases in a flat-bottom potential box [48, 49] makes the measurement of the virial coefficients $c_n$ very natural, since the densities $\rho_\sigma$ are directly accessible there. This last prediction makes it however less motivating than the measurement of the cluster coefficients $b_n$ in an inhomogeneous gas according to the specific technique for the harmonically trapped case [5], in which the chemical potentials $\mu_\sigma$ are the relevant variables to be used.<fnref n="8" /> For all practical purposes, we also give the fourth cluster and virial coefficients in numerical form in Table 1.

## 1.7. *Outline of the rest of the article*

In the following, we explain how we were able to obtain accurate results in a reasonable computation time, in particular by means of a convergence acceleration of the sum on $\ell$ in expression (4) relying on an asymptotic expansion of the summand and playing an essential role near the Efimovian thresholds. To do so, we had to generalize the analytical method of reference [40], implemented for 3+1 bodies, to the more difficult case of 2+2 bodies. In the 3+1-body case, we also highlight an unexpected non-monotonic dependence of $B_{3,1}(\bar{\omega})$ on $\bar{\omega}$, with change of sign, for mass ratios $\alpha$ close to 3.5 or 6.6. Section 2 deals with the 3+1-body case and Section 3 with the 2+2-body case in a harmonic potential.

# 2. Cluster coefficient for 3 + 1 fermions in a trap

We explain in this section how to perform an efficient numerical calculation of the quantity $I_{3,1}(\bar{\omega})$ defined by equation (4), for any value (zero or positive) of the reduced trapping frequency $\bar{\omega} = \hbar\omega/k_B T$. The desired cluster coefficient $B_{3,1}(\bar{\omega})$ follows directly from conjecture (5).

## 2.1. *Formulation of the problem*

Let us first recall the expression of the operator $M_{3,1}^{(\ell,\epsilon)}(\text{iS})$ involved in equation (4), as was obtained in reference [39]. We are dealing with the sum of a diagonal part $\mathcal{D}_{1,3}$ and a kernel

<footnote n="8">⁸In a trapping potential $U(\mathbf{r})$ common to both spin states and very elongated along the eigenaxis Oz, the cold-atom gases are well described by the local density approximation, and thus have well-defined local chemical potentials on the axis $\mu_\sigma^{\text{loc}}(z) = \mu_\sigma - U(0,0,z)$, which makes the measurement of the grand-canonical equation of state, and thus of the cluster coefficients, very straightforward, following the clever proposal of reference [50], which is usable however only if $U(\mathbf{r})$ is harmonic.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 1** – For a spatially homogeneous unitary Fermi gas with two spin states ↑ and ↓, fourth cluster (a) and virial (b) coefficients defined by equations ((1), (11), (13)) and predicted by conjecture (5) of reference [41], as functions of the mass ratio $\alpha = m_↑/m_↓$. The ↑↓ zero-range interaction is assumed to be scale-invariant in the four-body problem, which forces one to restrict oneself to the mass ratios between the critical values of the 3+1- and 1+3-body Efimov effect, indicated by the vertical dotted lines. In (a), the thin vertical lines mark the points where $b_4$ vanishes (always with sign change), namely $\alpha \approx 1.535$ and $\alpha \approx 10.355$ on the $\alpha > 1$ side, and the inset is an enlargement (the calculated points are represented by disks on interpolation lines). In (b), the dashed curve neglects in $c_4$ the true interaction-induced four-body component, the one $c_4^{\text{corr}} = -3(b_{1,3} + b_{2,2} + b_{3,1})/2$ that depends on $b_{3,1}$, $b_{2,2}$ or $b_{1,3}$, to represent $c_4^{\text{no corr}} = c_4 - c_4^{\text{corr}}$.</figcaption></figure>

operator $\mathcal{K}_{3,1}$ acting on functions $f_{m_z}(x, u)$ of two continuous variables, the logarithm $x \in \mathbb{R}^+$ of the norm ratio of two wavevectors and the cosine $u = \cos\theta \in [-1, 1]$ of the angle between them, and a discrete variable, the magnetic quantum number $m_z$ along the quantization axis $Oz$, varying in steps of two between $-\ell$ and $\ell$ for parity $\varepsilon = (-1)^{\ell}$, and between $-\ell+1$ and $\ell-1$ for parity $\varepsilon = (-1)^{\ell-1}$ and $\ell \neq 0$. It reads in a mixed Dirac and Schrödinger notation:

$$ \langle x, u | \langle \ell, m_z | M_{3,1}^{(\ell, \varepsilon)} (\text{iS}) | f \rangle = \mathcal{D}_{3,1}(x, u) f_{m_z}(x, u) + \int_0^{+\infty} dx' \int_{-1}^{1} du' \sum_{m_z' \mid (-1)^{m_z'} = \varepsilon} \mathcal{K}_{3,1}^{(\ell)}(x, u, m_z; x', u', m_z') f_{m_z'}(x', u') \quad (14) $$

The diagonal part is independent of angular momentum and scaling exponent iS. As we took advantage of the fermionic exchange symmetry of the state vector to restrict to $x > 0$, the kernel $\mathcal{K}_{3,1}$ is obtained by symmetrization of a primitive kernel $K_{3,1}$:<fnref n="9" />

$$ \mathcal{D}_{3,1}(x, u) = \frac{(1 + 2\alpha + \alpha u / \text{ch } x)^{1/2}}{1 + \alpha} ; \quad (15) \\ \mathcal{K}_{3,1}^{(\ell)}(x, u, m_z; x', u', m_z') = \sum_{\eta, \eta' = \pm 1} (\eta\eta')^{\ell+1} K_{3,1}^{(\ell)}(\eta x, u, \eta m_z; \eta' x', u', \eta' m_z') $$

<footnote n="9">⁹ In Dirac notation, $\mathcal{K}_{3,1}^{(\ell)} = (1+U)K_{3,1}^{(\ell)}(1+U)$ with the Hermitian involution $U = -P_x e^{i\pi L_x/\hbar}$, the parity operator $P_x$ changing $|x\rangle$ into $|-x\rangle$ and the rotation operator of axis $Ox$ of angle $\pi$ such that $\exp(i\pi L_x/\hbar)|\ell, m_z\rangle = (-1)^{\ell}|\ell, -m_z\rangle$ [40].</footnote>
<table><caption>**Table 1** – Numerical values of the fourth cluster coefficients $b_{n,4-n}$ and virial coefficients $c_{n,4-n}$ of the spatially homogeneous unitary gas of two-spin-state ↑ and ↓ fermions, tabulated as functions of the mass ratio $\alpha = m_↑/m_↓$, for choice (12) of the reference mass $\bar{m}$. As we go from $b_{n,4-n}$ to $b_{4-n,n}$ and from $c_{n,4-n}$ to $c_{4-n,n}$ by changing $\alpha$ to $1/\alpha$, we restrict to $\alpha \ge 1$. We do not give the coefficients associated with the integers $n=0$ and $n=4$ because they are identical to those of the ideal gas, but we give the mean coefficients $b_4$ and $c_4$ useful in the case of equal fugacities or phase-space densities in the two spin states, see equations (11) and (13). The uncertainties on $b_{3,1}$, $b_{1,3}$ and $b_{2,2}$, not specified, are less than one percent. The given values are shown graphically in Figure 1.</caption><thead><tr><th>α</th><th>1</th><th>1.5</th><th>2</th><th>2.5</th><th>3</th><th>3.5</th><th>4</th><th>4.5</th><th>5</th><th>5.5</th></tr></thead><tbody><tr><td>b<sub>3,1</sub></td><td>0.1837</td><td>0.1604</td><td>0.1247</td><td>0.08238</td><td>0.03726</td><td>-0.007203</td><td>-0.04762</td><td>-0.07982</td><td>-0.09941</td><td>-0.1016</td></tr><tr><td>b<sub>1,3</sub></td><td>0.18374</td><td>0.18672</td><td>0.18199</td><td>0.17641</td><td>0.17133</td><td>0.16696</td><td>0.16325</td><td>0.16010</td><td>0.15740</td><td>0.15508</td></tr><tr><td>b<sub>2,2</sub></td><td>-0.2445</td><td>-0.2781</td><td>-0.3469</td><td>-0.4321</td><td>-0.5283</td><td>-0.6328</td><td>-0.7448</td><td>-0.8637</td><td>-0.9890</td><td>-1.1209</td></tr><tr><td>b<sub>4</sub></td><td>0.03026</td><td>0.003256</td><td>-0.05135</td><td>-0.1179</td><td>-0.19110</td><td>-0.2678</td><td>-0.3458</td><td>-0.4229</td><td>-0.4968</td><td>-0.5649</td></tr><tr><td>c<sub>3,1</sub></td><td>-2.2558</td><td>-4.4056</td><td>-7.3318</td><td>-11.100</td><td>-15.773</td><td>-21.407</td><td>-28.050</td><td>-35.737</td><td>-44.495</td><td>-54.332</td></tr><tr><td>c<sub>1,3</sub></td><td>-2.2558</td><td>-1.2265</td><td>-0.8296</td><td>-0.6292</td><td>-0.5114</td><td>-0.4350</td><td>-0.3820</td><td>-0.3433</td><td>-0.3139</td><td>-0.2910</td></tr><tr><td>c<sub>2,2</sub></td><td>-7.1913</td><td>-7.3897</td><td>-7.7991</td><td>-8.3081</td><td>-8.8748</td><td>-9.4800</td><td>-10.111</td><td>-10.758</td><td>-11.415</td><td>-12.078</td></tr><tr><td>c<sub>4</sub></td><td>-5.8513</td><td>-6.5108</td><td>-7.9801</td><td>-10.018</td><td>-12.580</td><td>-15.661</td><td>-19.271</td><td>-23.419</td><td>-28.112</td><td>-33.350</td></tr><tr><td>α</td><td>6</td><td>6.5</td><td>7</td><td>7.5</td><td>8</td><td>8.5</td><td>9</td><td>9.5</td><td>10</td><td>10.5</td></tr><tr><td>b<sub>3,1</sub></td><td>-0.07925</td><td>-0.02715</td><td>0.06412</td><td>0.2022</td><td>0.4010</td><td>0.6724</td><td>1.0339</td><td>1.5059</td><td>2.1195</td><td>2.9095</td></tr><tr><td>b<sub>1,3</sub></td><td>0.15305</td><td>0.15128</td><td>0.14972</td><td>0.14833</td><td>0.14708</td><td>0.14597</td><td>0.14496</td><td>0.14404</td><td>0.14321</td><td>0.14244</td></tr><tr><td>b<sub>2,2</sub></td><td>-1.2590</td><td>-1.4034</td><td>-1.5542</td><td>-1.7113</td><td>-1.8749</td><td>-2.0452</td><td>-2.2223</td><td>-2.4077</td><td>-2.6003</td><td>-2.8015</td></tr><tr><td>b<sub>4</sub></td><td>-0.6239</td><td>-0.6709</td><td>-0.7014</td><td>-0.7116</td><td>-0.6947</td><td>-0.6447</td><td>-0.5530</td><td>-0.4101</td><td>-0.2001</td><td>0.09397</td></tr><tr><td>c<sub>3,1</sub></td><td>-65.248</td><td>-77.214</td><td>-90.190</td><td>-104.10</td><td>-118.85</td><td>-134.28</td><td>-150.20</td><td>-166.33</td><td>-182.32</td><td>-197.64</td></tr><tr><td>c<sub>1,3</sub></td><td>-0.2725</td><td>-0.2574</td><td>-0.2449</td><td>-0.2343</td><td>-0.2252</td><td>-0.2174</td><td>-0.2105</td><td>-0.2045</td><td>-0.1991</td><td>-0.1944</td></tr><tr><td>c<sub>2,2</sub></td><td>-12.740</td><td>-13.397</td><td>-14.044</td><td>-14.676</td><td>-15.288</td><td>-15.873</td><td>-16.423</td><td>-16.926</td><td>-17.374</td><td>-17.747</td></tr><tr><td>c<sub>4</sub></td><td>-39.130</td><td>-45.434</td><td>-52.239</td><td>-59.505</td><td>-67.180</td><td>-75.184</td><td>-83.417</td><td>-91.733</td><td>-99.948</td><td>-107.79</td></tr><tr><td>α</td><td>11</td><td>11.5</td><td>12</td><td>12.5</td><td>12.75</td><td>13</td><td>13.1</td><td>13.2</td><td>13.3</td><td>13.3842</td></tr><tr><td>b<sub>3,1</sub></td><td>3.9362</td><td>5.2839</td><td>7.0993</td><td>9.7004</td><td>11.522</td><td>14.042</td><td>15.402</td><td>17.152</td><td>19.728</td><td>26.101</td></tr><tr><td>b<sub>1,3</sub></td><td>0.14174</td><td>0.14109</td><td>0.14049</td><td>0.13993</td><td>0.13967</td><td>0.13941</td><td>0.13932</td><td>0.13922</td><td>0.13912</td><td>0.13904</td></tr><tr><td>b<sub>2,2</sub></td><td>-3.0120</td><td>-3.2360</td><td>-3.4725</td><td>-3.7278</td><td>-3.8654</td><td>-4.0133</td><td>-4.0767</td><td>-4.1430</td><td>-4.2140</td><td>-4.2784</td></tr><tr><td>b<sub>4</sub></td><td>0.5017</td><td>1.0632</td><td>1.8524</td><td>3.0250</td><td>3.8667</td><td>5.0528</td><td>5.7008</td><td>6.5428</td><td>7.7951</td><td>10.950</td></tr><tr><td>c<sub>3,1</sub></td><td>-211.58</td><td>-223.04</td><td>-230.30</td><td>-230.37</td><td>-225.99</td><td>-217.03</td><td>-211.71</td><td>-205.24</td><td>-198.00</td><td>-200.77</td></tr><tr><td>c<sub>1,3</sub></td><td>-0.1901</td><td>-0.1863</td><td>-0.1828</td><td>-0.1796</td><td>-0.1781</td><td>-0.1766</td><td>-0.1761</td><td>-0.1756</td><td>-0.1750</td><td>-0.1746</td></tr><tr><td>c<sub>2,2</sub></td><td>-18.025</td><td>-18.163</td><td>-18.114</td><td>-17.767</td><td>-17.414</td><td>-16.864</td><td>-16.561</td><td>-16.189</td><td>-15.720</td><td>-15.213</td></tr><tr><td>c<sub>4</sub></td><td>-114.90</td><td>-120.69</td><td>-124.30</td><td>-124.16</td><td>-121.79</td><td>-117.04</td><td>-114.22</td><td>-110.80</td><td>-106.95</td><td>-108.08</td></tr></tbody></table>

The matrix elements of the primitive kernel depend on iS and contain an integral over a rotation angle $\phi$ around axis Ox as follows,

$$ K_{3,1}^{(\ell)}(x, u, m_z; x', u', m_z') = \frac{(\lambda\lambda')^{3/2}}{[(1+\lambda^2)(1+\lambda'^2)]^{1/4}} \int_0^{2\pi} \frac{d\phi}{2\pi^2} \frac{(1+\lambda^2)^{iS/2} e^{-im_z\theta/\hbar} \langle l, m_z | e^{i\phi L_x/\hbar} | l, m_z' \rangle e^{im_z'\theta'/2} (1+\lambda'^2)^{-iS/2}}{1+\lambda^2+\lambda'^2+\frac{2\alpha}{1+\alpha}[\lambda u + \lambda' u' + \lambda\lambda'(uu' + v\nu'\cos\phi)]} \quad (16) $$

where $L_x$ is the angular momentum operator along Ox. Here, as in reference [39], we have put for abbreviation:

$$ \lambda = e^x, \quad \lambda' = e^{x'}, \quad \theta = \arccos u \in [0, \pi], \quad \theta' = \arccos u' \in [0, \pi], \quad v = \sin \theta, \quad v' = \sin \theta' \quad (17) $$

To evaluate $I_{3,1}(\bar{\omega})$ numerically, one must first replace the operator by a finite size matrix, by truncating the variable $x$ to $x_{max}$ and discretizing it according to the midpoint integration method, then by discretizing the variable $\theta$ (which we prefer to the variable $u$ because it leads to a smooth integrand) according to the Gauss-Legendre integration method. Then we compute
the determinant of the matrix by putting it in the Cholesky form, to take advantage of the fact that the operator $M_{3,1}^{(\ell, \epsilon)}$ (iS) is positive in the absence of 3 + 1-body Efimov effect. Finally, we compute the integral over S in the interval $[0, S_{max}]$ by the midpoint method (we need to know the logarithm of the determinant at integer multiples of the integration step dS to obtain its derivative at half-integer multiples). We take into account the contribution of the omitted interval $[S_{max}, +\infty[$ by means of an exponential approximation $A\exp(-BS)$ of the logarithmic derivative of the determinant justified by reference [40], where we have estimated the coefficients A and B by fitting on a neighborhood of $S_{max}$, in practice the interval $[S_{max} - 5/2, S_{max}]$.

## 2.2. Asymptotic approximant and applications

It remains to take into account the truncation on the angular momentum $\ell$ at some $\ell_{max}$, which is unavoidable in a numerical calculation. In practice, it is not reasonable to go beyond $\ell_{max} = 15$, because high values of $\ell$ lead to a complexity $O(\ell^3)$ and their numerical calculation would be very expensive in time. Unfortunately, this cut-off is not yet enough for large mass ratios $\alpha \gg 1$ if one aims at an accuracy on $B_{3,1}(0^+)$ better than one percent. The solution to this difficulty is to determine an asymptotic approximant $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ of the contribution $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ of angular momentum $\ell$ and parity $\epsilon$ to the desired quantity $I_{3,1}(\bar{\omega})$. Then, instead of neglecting completely the terms $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ for $\ell > \ell_{max}$, we replace them by $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ as follows,

$$I_{3,1}(\bar{\omega}) \approx \sum_{\ell=0}^{\ell_{\max}} \sum_{\epsilon} I_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) + \sum_{\ell=\ell_{\max}+1}^{+\infty} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) = \sum_{\ell=0}^{\ell_{\max}} \sum_{\epsilon} \left[ I_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) - J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) \right] + \sum_{\ell=0}^{+\infty} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)}(\bar{\omega}) \quad (18)$$

By using this method, we accelerate the convergence of the series: the error tends to zero more quickly with $\ell_{max}$. An exact asymptotic approximant to subleading order in $\ell$ is obtained by generalizing to $\bar{\omega} \neq 0$ the method of reference [40], which takes the kernel $\mathcal{K}_{3,1}^{(\ell, \epsilon)}$ as a small formal parameter and expands the logarithm of the determinant of $M_{3,1}^{(\ell, \epsilon)}$ to second order:

$$ \begin{aligned} \ln \det M_{3,1}^{(\ell, \epsilon)} &= \ln \det^{\epsilon} [\mathcal{D}_{3,1} + \mathcal{K}_{3,1}^{(\ell)}] \\ &= \ln \det^{\epsilon} [\mathcal{D}_{3,1}] + \ln \det^{\epsilon} [\mathbf{1} + \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)}] \\ &= \text{const} + \text{Tr}^{\epsilon} \left[ \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} - \frac{1}{2} \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} \mathcal{D}_{3,1}^{-1} \mathcal{K}_{3,1}^{(\ell)} + \dots \right] \end{aligned} \quad (19) $$

where the symbol $\epsilon$ in superscript of the trace and determinant means that we restrict ourselves to the subspace of values of $m_z$ compatible with parity $\epsilon$. The computation is done in Appendix A, and the corresponding expression of $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ is given in equation (47) in the form of multiple integrals. For this already very elaborate choice, the first values of the approximant ($0 \le \ell \le \ell_{max}$) are still easy to compute numerically with the same truncation and discretization as for the full determinant, but it would be tedious to go to larger values of $\ell$, which is however required by the second sum in the second side of equation (18); fortunately, the infinite series in the third side of equation (18) has a simple integral expression, see equation (59), which can even be expressed analytically for $\bar{\omega} = 0^+$ in terms of known functions such as the dilogarithm function, see equation (60). The efficiency of our asymptotic approximant can be seen in Figure 2, where it is compared to the numerical result. The figure shows, as a function of the mass ratio $\alpha$, in which angular momentum channels $\ell$ the approximant deviates by more than one percent from the exact value (this is the accuracy on the cluster coefficient we aim for in this work); these channels must be included in the numerical sum from 0 to $\ell_{max}$.

An interesting by-product of our convergence acceleration method is to obtain an asymptotic equivalent of the angular momentum $\ell$ and parity $\epsilon$ contribution to the cluster coefficient
<figure><img src="image_11.png" /><figcaption>**Figure 2** – Contrast between the contribution $I_{3,1}^{(\ell, \epsilon)}$ of the angular momentum $\ell$ and parity $\epsilon$ channel to $I_{3,1}$ and its asymptotic approximant $J_{3,1}^{(\ell, \epsilon)}$ given by equation (47), as a function of the mass ratio $\alpha = m_↑/m_↓$ and $\ell$, in the limit $\tilde{\omega} = 0^+$. The contrast between two real quantities $a$ and $b$ is here the ratio $|a-b|/(|a|+|b|)$. (a) Case $\epsilon = (-1)^\ell$. (b) Case $\epsilon = (-1)^{(\ell+1)}$ (which imposes $\ell > 0$). The numerical truncation and discretization parameters are those of Figure 3. The contrast saturates to one (as it happens for $\ell = 0$ and $\ell = 5$ in (a), for $\ell = 3$ in (b)) when the two compared quantities are of opposite signs.</figcaption></figure>

$B_{3,1}(0^+)$; it suffices to keep the contribution linear in $\mathcal{K}_{3,1}^{(\ell)}$ in equation (19) and to determine its dominant behavior at large $\ell$ with Cauchy's integral theorem, as done in Appendix A. We find:

$$ \frac{1}{2} \left\{ \frac{(1+\alpha)^2}{\sqrt{\alpha(1+3\alpha)}} \operatorname{Im} \left[ \left( \frac{\ell}{6\pi C_0} \right)^{1/2} \frac{z_0}{1+z_0} z_0^\ell \right] + \varepsilon \frac{(1+\alpha)^2}{\alpha^{3/2}} \operatorname{Re} \left[ \left( \frac{\ell}{2\pi C_1} \right)^{1/2} \frac{(1-z_1)}{\sqrt{1-\cos\xi_1}} z_1^\ell \right] \right\} \quad (20) $$

where

$$ C_0 = \frac{(1 - \cos \xi_0) (5/4 + \cos \xi_0)}{(1/2 + \cos \xi_0) (z_0 - 1/z_0)} \quad (21) $$

$$ \text{and} \quad C_1 = 4(z_1 - 1/z_1)(\cos\xi_1 - 1/2) \left[ (1 + \alpha^{-1})\cos\xi_1 + \frac{1 + 3\alpha^{-1}}{4} \right] $$

with the shorthand notation $\cos\xi_n \equiv (z_n + 1/z_n)/2$, $n = 0$ or $1$, the complex number $z_0$ given by equation (66) and the complex number $z_1 = -z_0^*$. This is the generalization to $3+1$ fermions of a result obtained for three bosons in reference [20], see equation (42) of this reference. Thus, to within a power law factor, $I_{3,1}^{(\ell, \epsilon)}(0^+)$ tends exponentially to zero with $\ell$, with irregular oscillations due to the fact that the argument of $z_0$ is in general not commensurable to $\pi$.

## 2.3. Born-Oppenheimer regime

In the limit $\alpha \to 0$, we notice that $z_0$ tends to zero in equation (66), $z_0 \sim i\sqrt{\alpha/3}$, so we expect the exponential suppression of the summand in (4) to become very fast for $3+1$ fermions. This expectation is confirmed numerically and, for $\alpha = 0$, only the $\ell = 0$ channel contributes. In this case, the single spin-↓ fermion is infinitely massive and behaves for the spin-↑ fermions as a fixed pointlike scatterer, of infinite s-wave scattering length and placed at the center of the trap. The
Born–Oppenheimer approximation becomes exact and gives the time-independent Schrödinger equation on the wave function $\Psi(\mathbf{r}_\downarrow)$ of the heavy particle:<fnref n="10" />

$$E_{\text{BO}} \Psi(\mathbf{r}_{\downarrow}) = \left[ -\frac{\hbar^2}{2m_{\downarrow}} \Delta_{\mathbf{r}_{\downarrow}} + \frac{1}{2} m_{\downarrow} \omega^2 r_{\downarrow}^2 + W(0) \right] \Psi(\mathbf{r}_{\downarrow}) \quad (m_{\downarrow} \to +\infty) \quad (22)$$

In the $n_\downarrow + 1$ fermion trapped problem, $W(0)$ is the energy of an eigenstate of $n_\downarrow$ non-interacting fermions in the presence of the scattering center. A spin-↑ fermion of orbital quantum numbers $(n, \ell, m_z)$ sees the scattering center only if $\ell = 0$, in which case its spectrum is lowered by $\hbar\omega$, thus having energy levels

$$\epsilon_{n,\ell,m_z} = \begin{cases} (2n + \ell + 3/2)\hbar\omega & \text{if } \ell > 0 \\ (2n + 1/2)\hbar\omega & \text{if } \ell = 0 \end{cases} \quad (n \in \mathbb{N}, -\ell \le m_z \le \ell) \quad (23)$$

Since the Born–Oppenheimer energy $E_{\text{BO}}$ is the sum of $W(0)$ and a vibrational energy level of a spin-↓ particle in the trap, we conclude that

$$\lim_{\alpha \to 0^+} Z_{n_\downarrow, 1} = Z_1 Z_{n_\downarrow}^{\text{scat}} \quad (24)$$

where $Z_{n_\downarrow}^{\text{scat}}$ is the canonical partition function of a ideal gas of $n_\downarrow$ fermions with spectrum (23), i.e. in the presence of the fixed scatterer, $Z_1$ is as in equation (2) the partition function of a single fermion and $Z_{n_\downarrow, 1}$ is that of the unitary gas of $n_\downarrow + 1$ trapped fermions. It is convenient to calculate its deviation $\Delta Z_{n_\downarrow, 1}$ from that of the ideal gas of $n_\downarrow + 1$ trapped fermions, since (23) differs from the ordinary spectrum only in the zero angular momentum channel. Taking into account Fermi statistics via the Pauli exclusion principle, and [41, equations (80)–(83)] which relate the cluster coefficients of the trapped system to the few-body partition functions, we find that

$$\begin{align} B_{1,1}(\bar{\omega}) &= Z_1^{-1} \Delta Z_{1,1} \xrightarrow{\alpha \to 0^+} \frac{1}{2 \operatorname{ch}(\bar{\omega}/2)} ; \notag \\ B_{2,1}(\bar{\omega}) &= Z_1^{-1} \Delta Z_{2,1} - Z_1 B_{1,1} \xrightarrow{\alpha \to 0^+} -\frac{1}{8 \operatorname{ch}^2(\bar{\omega}/2) \operatorname{ch}\bar{\omega}} \end{align} \quad (25)$$

$$B_{3,1}(\bar{\omega}) = Z_1^{-1} \Delta Z_{3,1} - Z_{2,0} B_{1,1} - Z_1 B_{2,1} \xrightarrow{\alpha \to 0^+} \frac{1}{16 \operatorname{ch}^3(\bar{\omega}/2) (\operatorname{ch}\bar{\omega})(2 \operatorname{ch}\bar{\omega} - 1)} \quad (26)$$

The value (26) of the limit is in perfect agreement with our numerical calculation of expression (4) not only for $\bar{\omega} = 0^+$, where $B_{3,1}(0^+) \to 1/16 = 0.0625$, as seen in Figure 3a, but, as we have verified, for all $\bar{\omega}$, which is an additional test of conjecture (5) in the 3 + 1-body case.<fnref n="11" /> As a side remark, let us point out that it is possible to go to the next orders by using the grand-canonical version of equation (24),

$$\Xi_{z_\downarrow \to 0 \text{ then } \alpha \to 0^+} = \Xi_\uparrow + z_\downarrow Z_1 \Xi_\downarrow^{\text{scat}} + O(z_\downarrow^2) \quad (27)$$

<footnote n="10">¹⁰The particles are at fixed temperature $T$ so, in the $m_\downarrow \to +\infty$ limit, the spin-↓ fermion occupies a region around the center of the trap of radius $R = O((k_B T/m_\downarrow \omega^2)^{1/2})$. On the other hand, the Born–Oppenheimer potential $W(\mathbf{r}_\downarrow)$ has an energy scale $k_B T$ and varies with a length scale $\lambda_\downarrow$, the de Broglie thermal wavelength of the light fermions, if $k_B T > \hbar\omega$, and it is of the order of $\hbar\omega$ and varies with a length scale $(\hbar/m_\downarrow\omega)^{1/2}$, the size of the ground vibrational state of a light fermion, otherwise. The variation of $W(\mathbf{r}_\downarrow)$ becomes negligible in all cases and we can replace it by $W(0)$. For the same reason, the so-called scalar or topological potential, which is added to the potential $W$ in the full adiabatic approximation [51], can be omitted.</footnote>

<footnote n="11">¹¹In the 2 + 1-body case, starting from the analytic integral expression of $B_{2,1}(\bar{\omega})$ of reference [21] in terms of the transcendental Efmov function $\Lambda_{2,1}^{(\ell)}$ on the pure imaginary axis, and by showing in the limit $\alpha \to 0^+$ that $\Lambda_{2,1}^{(\ell)}(iS)$ tends to 1 for $\ell > 0$ and tends to $1 + 1/\operatorname{ch}(\pi S/2)$ for $\ell = 0$ (for example by means of [38, equation (46)]), we find the second result in equation (25).</footnote>