<table><caption>**Table 2** – Numerical values of the fourth cluster coefficients $B_{3,1}(0^+)$ and $B_{1,3}(0^+)$ of the unitary Fermi gas in an infinitesimal-stiffness trap, tabulated as functions of the mass ratio $\alpha = m_\uparrow/m_\downarrow$. Since we go from one coefficient to the other by changing $\alpha$ to $1/\alpha$, we limited ourselves to $\alpha \ge 1$. The uncertainties, not specified, are less than one percent. These values are shown graphically in 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>

where $\Xi_\uparrow (\Xi_\downarrow^{\text{scat}})$ is the grand partition function of the ideal gas of spin-↑ fermions of fugacity $z_\uparrow$ in the absence (presence) of the scattering center and $\Xi$ is that of the two-component unitary gas. Expanding the grand potential $\Omega = -k_B T \ln \Xi$ as in equation (2), we obtain<fnref n="12" />

$$ \sum_{n_\uparrow=0}^{+\infty} z_\uparrow^{n_\uparrow} \lim_{\alpha \to 0^+} B_{n_\uparrow, 1}(\bar{\omega}) = \frac{\Xi_\uparrow^{\text{scat}}}{\Xi_\downarrow} = \prod_{n=0}^{+\infty} \frac{1+z_\uparrow 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_\uparrow^s}{2s \operatorname{ch}(s\bar{\omega}/2)} \right] \quad (28) $$

In the special case $\bar{\omega} = 0^+$, the exponential in the fourth side of (28) reduces to $\sqrt{1+z_\uparrow}$, easy to expand into powers of $z_\uparrow$, hence the result at all orders in terms of Euler's $\Gamma$ function:

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

## 2.4. The results

Our numerical results for $I_{3,1}(0^+)$, hence for the cluster coefficient $B_{3,1}(0^+)$ from equation (5), are given explicitly in Table 2 and are plotted as functions of the mass ratio in Figure 3a, with no error bars but with an uncertainty of less than one percent. We find that $B_{3,1}(0^+)$ is positive, except over an interval $[\alpha_1, \alpha_2]$, as best seen in the enlargement included in the figure. At large values of $\alpha$, $B_{3,1}(0^+)$ shows a marked growth that we attribute to the three-body Efimov effect; it has a finite limit at the threshold of the 3+1-body Efimov effect where our computation stops, reached with an infinite slope and that we determine by extrapolation, see the dashed curves in Figure 3a and the explanations which follow,

$$ B_{3,1}(0^+) \underset{\alpha \to (\alpha_c^{\text{body}})^-}{\longrightarrow} 2.47 \pm 0.03 \quad (30) $$

## 2.5. Near the 3+1-body Efimovian threshold

The limit $\alpha \to (\alpha_c^{\text{body}})^-$ in equation (30) is difficult to achieve numerically with precision for several reasons. First, if $\alpha \to +\infty$, the exponential decay rate $\ln(1/|z_0|)$ in equation (20) tends to zero as $\alpha^{-1/2}$ so the convergence of the series $I_{3,1}(0^+)$ becomes slower and slower with $\ell$;

<footnote n="12">¹²A sum over an integer $s$ in the fourth side was introduced by taking the logarithm of the infinite product in the third side, expanding around 1 the function $\ln$ into an integer series of index $s$ and then summing over $n$ the resulting geometric series.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 3** – (a) Fourth cluster coefficient $B_{3,1}(0^+)$ of the unitary Fermi gas in an infinitesimally stiff harmonic trap, as a function of the mass ratio $\alpha = m_1/m_1$. Disks (connected by a thin line): Our numerical calculation of expression (4) for $(n_1, n_1) = (3, 1)$ and application of conjecture (5), with truncation parameters $x_{max} = 5$, $S_{max} = 25$, $\ell_{max} = 15$ and discretization $dx = 1/10$, $dS = 1/8$, $n_\theta = 37$ (this is the number of values of $\theta$ in the Gauss–Legendre integration method); for $\alpha > 10$ and $0 \le \ell \le 5$, we take instead $x_{max} = 20$ and $n_\theta = 25$ to keep the error below one percent. We use the convergence acceleration method (18) and extrapolate to $S_{max} = +\infty$ as explained in the text. Blue dashed line: cubic fit in the variable $X = (\alpha_c^{4\text{body}} - \alpha)^{1/2}$ on the interval $10 \le \alpha \le 13.3$, i.e. $B_{3,1}(0^+) = 2.4401 - 2.1627X + 0.64137X^2 - 0.069596X^3$, where $\alpha_c^{4\text{body}} \approx 13.3842$ is the 3 + 1-body Efimovian threshold (vertical solid line). Red dashed line: same on the interval $12.5 \le \alpha \le 13.3$, i.e. $B_{3,1}(0^+) = 2.4637 - 2.2834X + 0.82318X^2 - 0.15269X^3$. The inset is an enlargement showing better $B_{3,1}(0^+)$ vanishing with sign change at $\alpha_1 \approx 3.412$ and $\alpha_2 \approx 6.678$ (vertical black dotted lines). (b) Same for the half-second derivative $(1/2)B_{3,1}''(0^+)$ of the cluster coefficient with respect to the reduced trapping frequency $\bar{\omega} = \hbar\omega/k_B T$, with the differences that (i) we do not give a fit in the variable $X$, (ii) we put error bars because the relative uncertainty may now exceed one percent, (iii) we include a magnification in the enlargement to make it clear that the second zero (with sign change) of $(1/2)B_{3,1}''(0^+)$ is very close to but distinct from that of $B_{3,1}(0^+)$ (always indicated by a vertical dotted line).</figcaption></figure>

fortunately, as the mass ratio $\alpha$ remains below the four-body Efimovian threshold, $\ln(1/|z_0|)$ does not become so small (it remains above 0.267) and the problem is solved by the convergence acceleration (18) if one is satisfied with an error on $B_{3,1}(0^+)$ less than one percent (see Figure 2). Secondly, the 3 + 1-body Efimov effect about to appear in the $\ell = 1, \varepsilon = +1$ channel forces to increase $x_{max}$ to the values considered in Reference [39], at least in the $\ell \le 5$ angular momentum channels; this increases the computation time a lot and leads to a singularity of type

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

which has to be taken into account in the extrapolation of $B_{3,1}(0^+)$ to $\alpha = \alpha_c^{4\text{body}}$, see dashed lines in Figure 3a and the legend of the figure.<fnref n="13" /> Third, the threshold for the three-body Efimov

<footnote n="13">¹³This singularity is present in one of the scaling exponents $s_i$ of the 3 + 1-body unitary problem, the one $s_0$ whose square vanishes by changing sign at the Efimovian threshold, in the $\ell_0 = 1$ channel [39]. Now $B_{3,1}(0^+)$ depends linearly on the scaling exponents. We see it well on the universal component of the third cluster coefficient of the unitary Bose gas in equation (36) of reference [20] by making $\bar{\omega}$ (denoted $x$ in this reference) tend to zero; we also see it for $2 + 1$</footnote>
effect $\alpha_c^{2,1} \approx 13.60697$ is close by and affects all angular momentum channels $\ell$ of the four-body problem [40], even though it occurs in the $L=1$ angular momentum channel of the $2+1$ fermion problem. Indeed, the continuum spectrum of the $M_{3,1}^{(\ell,\epsilon)}$ (iS) operator is the union of continua corresponding to $2+1$ fermions strongly correlated by the interactions and a decoupled spin-↑ fermion, thus parametrized by a three-body angular momentum $L$ [39]. The contribution of the continuum of angular momentum $L$ to the $3+1$ cluster coefficient in the $(\ell, \epsilon)$ channel is written up to a factor [40]:

$$I_{3,1}^{(\ell,\epsilon)}(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)$$

The eigenmodes of the continuum have asymptotically a plane wave structure in the space of the variable $x$, i.e. they are superposition, when $x \to +\infty$, of an incident wave $\exp(-ikx)$ and a reflected wave $-\exp[i\theta_L(k, S)]\exp(ikx)$. Here the phase shift $\theta_L(k, S)$ is a function of the wavenumber $k > 0$ (fictitious because $x$ is not a real position), of the pure imaginary scaling exponent $iS$ and of the three-body angular momentum $L$. We see in expression (31) the transcendental Efimov function $\Lambda_{2,1}^{(L)}$ whose roots are the scaling exponents of the unitary $2+1$-body problem with angular momentum $L$. This is because the continuum modes of wavenumber $k$ correspond to the eigenvalue $\Lambda_{2,1}^{(L)}(ik)$ of $M_{3,1}^{(\ell,\epsilon)}$ (iS) [39].$^{14}$ The key point now is that the lower edge of the $L=1$ continuum, namely the minimum of $\Lambda_{2,1}^{L=1}(ik)$ with respect to the variable $k$, tends to zero when $\alpha \to \alpha_c^{2,1}$ because of the $2+1$-body Efimov effect. This has two consequences. First, a practical consequence in the calculation of $I_{3,1}$ on a computer: as the continuum is made discrete by truncating $x$ at $x_{max}$, its lower edge deviates from the true edge by about $1/x_{max}^{2}$. $^{15}$ this numerical artifact forces to significantly increase $x_{max}$, fortunately only for four-body angular momentum channels $\ell \le 5$ as long as $\alpha \le 13.3$. Then, a physical consequence: a second singularity appears in $I_{3,1}(0^{+})$, of the form $(\alpha_c^{2,1} - \alpha)^{1/2}$, again a square root, now centered on the three-body Efimovian threshold.$^{16}$ This three-body singularity occurs at a point very close to the four-body singularity, which casts doubt on the accuracy of the extrapolation made in Figure 3a

$$d B_{3,1}(0^{+}) / dX = -(\ell_0 + 1/2) c_0^{1/2} \approx -2.2 \quad \text{at} \quad X \equiv (\alpha_c^{3,1} - \alpha)^{1/2} = 0.$$

This constraint is fairly well satisfied by the fits in Figure 3a.

$$\propto \int_{0}^{+\infty} dk \left[ k^2 / (k^2 + \sigma_0^2) - k^2 / (k^2 + 1) \right] \propto 1 - \sigma_0 \quad \text{where} \quad \sigma_0 \sim \gamma_0 (\alpha_c^{2,1} - \alpha)^{1/2} \quad \text{and} \quad \gamma_0 = 0.438.$$

<footnote>fermions in reference [21]. More precisely, $B_{3,1}(0^+) = -(\ell_0 + 1/2)s_0 + \dots$ where the ellipse is a smooth function of $\alpha$ in the neighborhood of $\alpha_c^{3,1}$. As $s_0^2 = c_0(\alpha_c^{3,1} - \alpha)$ with $c_0 \approx 2.2$ near the threshold [39], we find that</footnote>

<footnote>$^{14}$ In the absence of a cutoff in $x$ space ($x_{max} = +\infty$), we fall into a paradox: $k$ spans the continuous set $\mathbb{R}^+$ independently of the variable $S$ and the derivative of the eigenvalue $\Lambda_{2,1}^{(L)}(ik)$ with respect to $S$ is zero, so the continuum should not contribute to $I_{3,1}$ in equation (4). The right way to reason is to put a cutoff $x_{max}$ that we make tend to infinity at the end, with the condition that the eigenmodes vanish at $x = x_{max}$. The resulting equation $\theta_L(k, S) = -2kx_{max}[2\pi]$ quantizes $k$, i.e. restricts it to a discrete set, and makes it $S$-dependent as the phase shift $\theta_L(k, S)$.</footnote>

<footnote>$^{15}$ In fact, the wavenumber $k$ varies in steps of the order of $1/x_{max}$, as shown in Footnote 14, and $\Lambda_{2,1}^{L=1}(ik)$ varies quadratically near its minimum.</footnote>

<footnote>$^{16}$ To the left of this threshold, we write à la Weierstrass $\ln \Lambda_{2,1}^{(L=1)}(ik) = \ln[(k^2 + \sigma_0^2)/(k^2 + 1)] + \dots$ where the ellipse is a regular function of $\alpha$ even at the threshold and $\sigma_0 > 0$ is the scaling exponent of the $2+1$-body problem whose square vanishes by changing sign at the threshold. Very close to the threshold, $\sigma_0 \ll 1$ and, assuming as in reference [40] that $\theta_{L=1}(k, S) \sim b(S)k$ when $k \to 0$, where the scattering radius $b(S)$ is an unknown function of $S$, we find that $I_{3,1}^{(\ell,\epsilon)\neq(1,+)}(0^+)$ contains a singularity</footnote>

<footnote>We have confirmed this prediction by a specific numerical calculation of $I_{3,1}(0^+)$ in the $\ell=0$ channel very close to the threshold, going up to $\alpha = 13.59$; a nice linear law in $Y$ is observed and a fit gives $I_{3,1}^{(\ell=0)}(0^+) \approx 0.026 + 0.040Y$ where the variable $Y$ is the one in equation (32).</footnote>
which did not take it into account. We remedy this by fitting functions including the two singularities, polynomial in the quantities $X = (\alpha_c^{3,1} - \alpha)^{1/2}$ and $Y = (\alpha_c^{2,1} - \alpha)^{1/2}$:

$$ 
\begin{aligned}
I_{3,1}(0^+) &= A_0 + A_1 X + A_2 (Y - Y_0) + A_3 X^2 \\
\text{and} \quad 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{aligned} 
\quad (32) $$

where $Y_0 = (Y^2 - X^2)^{1/2} = (\alpha_c^{2,1} - \alpha_c^{3,1})^{1/2}$ is also the value of $Y$ at $X = 0$. On the interval $10 \le \alpha \le 13.3$, this leads to $A_0 = 2.495$ and $A_0 = 2.447$ hence the final result (30), which supports the more naive one in Figure 3a.

## 2.6. Case $\bar{\omega} \neq 0$

To conclude this section, let us briefly study the dependence on trap stiffness of the cluster coefficient $B_{3,1}(\bar{\omega})$. One experimentally useful way to account for this is to calculate the first deviation from the zero stiffness limit, of even degree in $\bar{\omega}$ since the integrand of equation (4) is an even function of $\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) \quad (33) $$

This allows to quantify the error due to the local density approximation, systematically used in the experiment and which amounts to keeping only the first term of Taylor expansion (33). It is easy to see that the small parameter controlling this approximation can only be $\bar{\omega}$ in the unitary limit: the approximation only makes sense if the equilibrium correlation length of the homogeneous gas in each spin state $\sigma$, i.e. the de Broglie thermal length $\lambda_\sigma$ in the nondegenerate case due to scale invariance, is much smaller than the spatial radius of the trapped gas $R_\sigma = (k_B T / m_\sigma \omega^2)^{1/2}$; one has indeed $\lambda_\sigma / R_\sigma \propto \hbar\omega / k_B T = \bar{\omega} \ll 1$. To obtain an integral expression for the second derivative $B_{3,1}''(0^+)$ and implement convergence acceleration, we take twice the derivative with respect to $\bar{\omega}$ of equation (4) under the integral sign, equation (18) under the sum sign and expressions (47), (59)) and then make $\bar{\omega}$ tend to zero.

The result is plotted as a function of the mass ratio in Figure 3b.<fnref n="17" /> Let us try to interpret it in a naive scenario: as in the case of equal masses $m_\uparrow = m_\downarrow$ studied in reference [41], $B_{3,1}(\bar{\omega})$ would simply be a monotonic function of $\bar{\omega}$, of course of zero limit at infinite $\bar{\omega}$. The second derivative at the origin would then always have the opposite sign to the value at $\bar{\omega} = 0$. This explains Figure 3b if we look on a large scale: $B_{3,1}''(0^+)$ seems indeed to vanish by changing sign at the same mass ratios $\alpha = \alpha_1$ and $\alpha = \alpha_2$ as $B_{3,1}(0^+)$. However, we can see on a first enlargement, in inset in Figure 3b, that this scenario fails at $\alpha_1$ (it still seems to hold at $\alpha_2$). Indeed, as shown in Figure 4a, when $\alpha$ approaches $\alpha_1 \approx 3.412$ from lower values (top to bottom curves), the function $\bar{\omega} \mapsto B_{3,1}(0^+)$ ceases to be monotonic, becomes convex near the origin and reaches an absolute (positive) maximum at a point $\bar{\omega}_0 > 0$ that departed from zero before $B_{3,1}(0^+)$ becomes $< 0$. This absolute maximum persists when $\alpha$ continues to grow beyond $\alpha_1$, but its position moves toward $+\infty$ and out of Figure 4a. A second enlargement, in the inset of Figure 3b, shows that the simple scenario also fails (but just barely) at the second nodal point $\alpha_2 \approx 6.678$: when $\alpha$ approaches $\alpha_2$ by higher values (top-down curves in Figure 4b<fnref n="18" />), an absolute (negative) minimum coming

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

<footnote n="17">¹⁷For values of $\alpha$ close to $\alpha_2$, we improve the convergence acceleration method by approximating</footnote>

<footnote>rather than by 0 as in equation (18); the real parameters $A$, $\psi$ and $\gamma$ are obtained by fitting on the interval $6 \le \ell \le \ell_{\max}$, and the complex number $z_0$ is that of the asymptotic law (20).</footnote>

<footnote n="18">¹⁸By varying all truncation and discretization parameters in the numerical calculation, we estimate the uncertainty on the curves in Figure 4b closest to the horizontal axis to be $2 \times 10^{-5}$ (e.g., for $\alpha = 6.6875$); the plotted dependencies on $\bar{\omega}$ are thus significant.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 4** – Fourth cluster coefficient $B_{3,1}(\bar{\omega})$ of the trapped unitary Fermi gas as a function of the reduced stiffness $\bar{\omega}^2 = (\hbar\omega/k_B T)^2$ of the harmonic potentials, for different mass ratios $\alpha = m_\uparrow/m_\downarrow$ close to the nodal points $\alpha_1 \approx 3.412$ (a) and $\alpha_2 \approx 6.678$ (b) of $B_{3,1}(0^+)$, which Figure 3b, i.e., the relationship between the signs of the value and curvature of $B_{3,1}(\bar{\omega})$ at the origin. In (b), the slopes of the tangents at the origin (dotted) are taken from Figure 3b.</figcaption></figure>

from the $\bar{\omega} = +\infty$ side approaches the origin and becomes more pronounced, before $B_{3,1}(0^+)$ in turn becomes $< 0$. For convenience, we give the values of $B_{3,1}(\bar{\omega})$ and its half-second derivative at $\bar{\omega} = 0^+$ in numerical form in Table 3, for the mass ratios of Figure 4. To conclude, let us finally point out that the simple scenario assuming monotonicity of $B_{3,1}(\bar{\omega})$ was actually highly improbable because it implied that the 3 + 1-body cluster coefficient in the trap was zero for any stiffness, $B_{3,1}(\bar{\omega}) \equiv 0$, at mass ratios where $B_{3,1}(0^+) = 0$.

<table><caption>**Table 3** – Numerical values of the fourth cluster coefficient $B_{3,1}(0^+)$ and its half-second derivative $(1/2)B''_{3,1}(0^+)$ with respect to $\bar{\omega}$, given with their significant digits, for the mass ratios of Figure 4. The half-derivatives are plotted (some with an error bar) in 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><sup>1</sup>&frasl;<sub>2</sub>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><td>6.9375</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><td>0.00469</td></tr><tr><td><sup>1</sup>&frasl;<sub>2</sub>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><td>-0.0115</td></tr></tbody></table>

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

We explain in this section how to perform an efficient numerical calculation of the quantity $I_{2,2}(\bar{\omega})$ defined by equation (4), limiting ourselves for simplicity to the case $\bar{\omega} = 0^+$, i.e. to a trap of infinitesimal stiffness. The unitary gas cluster coefficient $B_{2,2}(0^+)$ can be deduced by conjecture (5). As the two spin states $\uparrow$ and $\downarrow$ play perfectly symmetric roles here, $B_{2,2}(0^+)$ is invariant by changing the mass ratio $\alpha$ into its inverse $1/\alpha$ and one can limit the numerical calculations to the case $\alpha \ge 1$. In contrast to the previous Section 2, one can go up to the three-body Efimovian threshold, $\alpha < \alpha_c^{2,1}$, since there is no 2 + 2-body Efimov effect [40].
### 3.1. Formulation of the problem

Recall the expression of the Hermitian operator $M_{2,2}^{(\ell, \epsilon)}$ (iS) given in reference [40] for angular momentum $\ell$, parity $\epsilon = \pm$ (restricted to + if $\ell = 0$) and the pure imaginary scale exponent iS, with the same Schrödinger–Dirac notation as in equation (14) but with the difference that the variable $x$ varies on the whole real axis:

$$ \begin{aligned} & \left\langle x, u \middle| \left\langle \ell, m_z \middle| M_{2,2}^{(\ell, \epsilon)} \middle| f \right\rangle \right. \\ & = \mathcal{D}_{2,2}(x, u) f_{m_z}(x, u) + \int_{-\infty}^{+\infty} dx' \int_{-1}^{1} du' \sum_{\substack{m'_z \\ m'_z | (-1)^{m'_z} = \epsilon}} K_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) f_{m'_z}(x', u') \quad (34) \end{aligned} $$

with a diagonal part independent of iS, angular momentum and magnetic quantum number $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) $$

and a matrix integral kernel decomposed into three contributions $K_{2,2} = K_1 + K_2 + K_3$ written line by line in this order:

$$ \begin{aligned} K_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) & = \left(\frac{e^x \text{ch } x'}{e^{x'} \text{ch } x}\right)^{iS/2} \left(\frac{e^{x+x'}}{4 \text{ch } x \text{ch } x'}\right)^{1/4} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{e^{-im_z\theta/2} \langle \ell, m_z | e^{i\phi L_x/\hbar} | \ell, m'_z \rangle e^{im'_z\theta'/2}}{\text{ch}(x-x') + \frac{1}{1+\alpha} [(u+e^{-x})(u'+e^{-x'}) + v\nu'\cos\phi]} \\ & \quad + \left(\frac{e^{-x} \text{ch } x'}{e^{-x'} \text{ch } x}\right)^{iS/2} \left(\frac{e^{-x-x'}}{4 \text{ch } x \text{ch } x'}\right)^{1/4} \int_0^{2\pi} \frac{d\phi}{(2\pi)^2} \frac{e^{im_z\theta/2} \langle \ell, m_z | e^{i\phi L_x/\hbar} | \ell, m'_z \rangle e^{-im'_z\theta'/2}}{\text{ch}(x-x') + \frac{\alpha}{1+\alpha} [(u+e^x)(u'+e^{x'}) + v\nu'\cos\phi]} \\ & \quad + \frac{(-1)^{\ell+1}}{4\pi[(u+\text{ch }x)(u'+\text{ch }x')\text{ch }x\text{ch }x']^{1/4}} \frac{\left(\frac{u'+\text{ch }x'}{u+\text{ch }x}\right)\text{ch }x'}{\left(\frac{u'+\text{ch }x'}{u+\text{ch }x}\right)\text{ch }x} \left(\frac{e^{im_z\gamma(x,u)} \langle \ell, m_z | \ell, m_x = 0 \rangle \langle \ell, m_x = 0 | \ell, m'_z \rangle e^{-im'_z\gamma(x',u')}}{\left(\frac{e^{-x'}+\alpha e^{x'}}{1+\alpha}\right)(u+\text{ch }x) + \left(\frac{e^{-x}+\alpha e^x}{1+\alpha}\right)(u'+\text{ch }x')} \right)^{iS/2} \end{aligned} \quad (36) $$

with the angle $\gamma(x, u) = \arctan\{\text{th}(x/2)[(1-u)/(1+u)]^{1/2}\}$, the angle $\theta \in [0, \pi]$ such that $u = \cos\theta$ and the notation $v = \sin\theta$. The third contribution $K_3$ is nonzero only in the parity sector $\epsilon = (-1)^\ell$ (in the other sector, we have $\langle \ell, m_x = 0 | \ell, m_z \rangle = 0$). Contrary to the first two contributions, it is a singular function of the coordinates $(x, u)$ and $(x', u')$, diverging at the point $(x = 0, u = -1)$, or equivalently $(x = 0, \theta = \pi)$. In the numerical calculation, limited to the rectangle $(x, \theta) \in [-x_{\text{max}}, x_{\text{max}}] \times [0, \pi]$, we isolate the singularity by a half disk of radius $R$ and center $(0, \pi)$, in which we use a logarithmic-polar grid as in reference [40]$^{19}$, and outside of which we use the same type of grid as in Section 2.

### 3.2. Asymptotic approximant

It remains to implement the same convergence acceleration technique as in equation (18) in the particular case $\tilde{\omega} = 0^+$,

$$ I_{2,2}(0^+) \approx \sum_{\ell=0}^{\ell_{\text{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) $$

<footnote>$^{19}$ We set $(x, \pi - \theta) = R e^\ell (\cos \psi, \sin \psi)$, where $t_{\text{min}} < t < 0$ and $0 < \psi < \pi$. We have chosen $R = 2/5$. The integration on $\psi$ is done by the Gauss–Legendre method with 15 points (series $0 \le \ell \le 4$) or 25 points (series $0 \le \ell \le 6$), the integration on $t$ with the midpoint rule and a step $dt = 0.1$ or $dt = 0.2$. We extrapolate quadratically to $t_{\text{min}} = -\infty$ from the three choices $t_{\text{min}} = -10$, $t_{\text{min}} = -20$ and $t_{\text{min}} = -40$.</footnote>
where, as in Section 2, $I_{2,2}^{(\ell, \epsilon)}(0^+)$ is the contribution of the angular momentum $\ell$ and parity $\epsilon$ channel to the desired quantity $I_{2,2}(0^+)$, and $J_{2,2}^{(\ell, \epsilon)}(0^+)$ is an asymptotic approximant, valid at large angular momentum. Reference [40] gives such an approximant only in the tractable unnatural-parity sector $\epsilon = (-1)^{l-1}$ where $K_3 \equiv 0$, by perturbatively treating $K_1$ and $K_2$ to second order in the logarithm of the determinant of $M_{2,2}^{(\ell, \epsilon)}(\text{iS})$. The calculation is very similar to equation (19), with the simplification that the linear terms and the square terms are independent of $S$ and can be absorbed into the constant, so that only the crossed terms remain [40]:

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

In the more challenging natural-parity sector $\epsilon = (-1)^{\ell}$, a difficulty arises: because of its divergent character, we cannot treat $K_3$ perturbatively. We first perform a gauge transform on $K_j$ eliminating the $S$-dependence of $K_3$ without changing the determinant, which is indicated by a tilde,

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

then we isolate a purely “external” part of $\tilde{K}_3$, acting on the “orbital” space of $(x, u)$ but not on that of the $|\ell, m_z\rangle$'s, i.e. such that $\tilde{K}_3 = \tilde{K}_3^{\text{ext}} \otimes |\ell, m_x = 0\rangle\langle\ell, m_x = 0|$. Finally, we carry out an expansion of the logarithm of the determinant in powers of $K_1$ and $K_2$ without any hypothesis on $\tilde{K}_3^{\text{ext}}$:

$$ \begin{aligned} \ln \det [M_{2,2}^{(\ell, \epsilon)} (\text{iS})] &\stackrel{\epsilon = (-1)^{\ell}}{=} \text{const} + \ln \det \left[ \mathbf{1} + \frac{1}{\mathcal{D}_{2,2} + \tilde{K}_3^{\text{ext}} \otimes |\ell, m_x = 0\rangle\langle\ell, m_x = 0|} (\tilde{K}_1 + \tilde{K}_2) \right] \\ &= \text{const} - \text{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\text{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle] - \text{Tr}^{\epsilon} [\mathcal{D}_{2,2}^{-1} \tilde{K}_1 \mathcal{D}_{2,2}^{-1} \tilde{K}_2] \\ &\quad + \text{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\text{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] \\ &\quad - \frac{1}{2} \text{Tr}_{x,u} [\mathcal{D}_{2,2}^{-1} \tilde{K}_3^{\text{ext}} \mathcal{D}_{2,2}^{-1} \langle \ell, m_x = 0 | (\tilde{K}_1 + \tilde{K}_2) | \ell, m_x = 0 \rangle] \\ &\quad + \dots \end{aligned} \quad (40) $$

In this expression, a resummed form of the external part appears:

$$ \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) $$

An optimal writing of $\tilde{K}_3^{\text{ext}}$ is obtained by reparameterization of the Faddeev ansatz of the 2+2-body problem: we consider that the functions on which operator $M_{2,2}$ acts now depend on the relative wave vectors and the center of mass of the particles 2 and 4 as in reference [52] instead of the single-particle wave vectors $\mathbf{k}_2$ and $\mathbf{k}_4$ as in references [39,40]; this avoids the half-disk around the singularity and leads us to numerically invert an operator acting on a single real variable instead of the two variables $(x, u)$. This leads to a considerable saving of computation time and simplification. The reader is referred to Appendix C for more details. We find numerically that the second term in the third side of equation (40), formally of the first order, is actually of the same order of magnitude as the third term (their contributions to $J_{2,2}^{(\ell, \epsilon)}(0^+)$ tend to zero exponentially with $\ell$ with the same rate, see Figure 5a). The resummed kernel $\tilde{K}_3^{\text{ext}}$ is thus, like $K_1$ and $K_2$, a first-order infinitesimal; consequently, we neglect the fourth and fifth terms of equation (40) to keep

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

If one wants, one can undo the gauge transform (39) in the third term without changing the trace, to recover contribution (38). In Appendix B, we give a more explicit expression of the asymptotic
approximant $J_{2,2}^{(\ell, \epsilon)}(0^+)$ deduced from expansions (38) and (42), see equations ((71), (72), (78), (80)), as well as its sum over $\ell$ and over $\epsilon$, see equations ((81), (86)).

### 3.3. Application and results

The sum over $\ell$ in $I_{2,2}(0^+)$ shows the same phenomenon of slow convergence as in $I_{3,1}(0^+)$ at large values of the mass ratio $\alpha$. We again find that the asymptotic approximant, here $J_{2,2}^{(\ell, \epsilon)}(0^+)$, tends to zero exponentially in $\ell$, with a rate $c$ that we compute numerically for convenience in the tractable parity sector $\epsilon = (-1)^{\ell-1}$ and for a mass ratio $\alpha \gtrsim 3$: $^{20}$

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

The rate $c$ is plotted as a function of the mass ratio in Figure 5b. Here, as in Section 2, a $c \propto \alpha^{-1/2}$ law is observed at large mass ratios. In this regime, the convergence acceleration method is an indispensable aid to numerical calculation, as shown in Figure 5c. It allows us to obtain the fourth cluster coefficient $B_{2,2}(0^+)$ of the trapped system, given in numerical form in Table 4 and plotted as a function of the mass ratio in Figure 6, with an uncertainty of less than one percent in a reasonable computation time. Because of the 2+1-body Efimov effect, this coefficient has a $(\alpha_c^{2,1} - \alpha)^{1/2}$ singularity near the threshold, for the same reason as discussed around equation (31) in Section 2. We take this into account in the dashed fits in Figure 6, to obtain the extrapolation

$$B_{2,2}(0^+) \xrightarrow{\alpha \to (\alpha_c^{\text{body}})^{-1}} -0.737 \pm 0.007 \quad (44)$$

<figure><img src="image_14.png" /><figcaption>**Figure 5** – Some numerical results on the asymptotic approximant $J_{2,2}^{(\ell, \epsilon)}(0^+)$ of the contribution $I_{2,2}^{(\ell, \epsilon)}(0^+)$ of angular momentum $\ell$ and parity $\epsilon$ to the quantity $I_{2,2}(0^+)$ in equation (4). (a) For a mass ratio $\alpha = 5$ and the challenging natural-parity sector $\epsilon = (-1)^{\ell}$, absolute value of the three terms of $J_{2,2}^{(\ell, \epsilon)}(0^+)$ in expression (71) (symbols ◦, +, and × in that order) as a function of angular momentum $\ell$. (b) Rate $c$ of exponential decay with $\ell$ of $J_{2,2}^{(\ell, \epsilon)}(0^+)$ as in equation (43), as a function of the square root of the inverse mass ratio. The error bars give the deviation between values of $c$ from various fitting functions (for $\epsilon = (-1)^{\ell-1}$, $\ln|J_{2,2}^{(\ell, \epsilon)}(0^+)|$ is taken as an affine function of $\ell$ plus one term $\propto 1/\ell$, or plus one term $\propto \ln \ell$, or plus both). The dashed line is an affine extrapolation to $\alpha^{-1/2} = 0$. (c) For $\alpha = 13.3$ and $\epsilon = (-1)^{\ell}$, exact value $I_{2,2}^{(\ell, \epsilon)}(0^+)$ obtained numerically (squares) and asymptotic approximant $J_{2,2}^{(\ell, \epsilon)}(0^+)$ (stars) as functions of $\ell$.</figcaption></figure>

<footnote>$^{20}$The dominant behaviors in the two parity sectors differ only in their power laws in $\ell$. When the mass ratio is too close to 1, $J_{2,2}^{(\ell, \epsilon=(-1)^{\ell-1})}(0^+)$ tends to zero while oscillating which makes numerical rate extraction more difficult.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 6** – Fourth cluster coefficient $B_{2,2}(0^+)$ of the unitary Fermi gas in a harmonic trap of infinitesimal stiffness, as a function of the mass ratio $\alpha = m_1/m_\downarrow$ restricted to $\alpha \ge 1$ by $\alpha \leftarrow 1/\alpha$ symmetry. Disks (connected by a thin line): Our numerical calculation of expression (4) for $(n_1, n_\downarrow) = (2,2)$ and application of conjecture (5), with parameters $x_{max} = -x_{min} = 20$ and $n_\theta = 15$ in the series $0 \le \ell \le 4$, $x_{max} = -x_{min} = 5$ and $n_\theta = 25$ in the series $0 \le \ell \le 6$ (the overlap of the two series serves as a check), and in all cases $dx = 1/15$, $S_{max} = 12$, $dS = 12/100$ (see footnote 19 for the treatment of the singularity of the integral kernel at $(x, \theta) = (0, \pi)$). We use the convergence acceleration method (37) and the same extrapolation to $S_{max} = +\infty$ as in Section 2. Thin vertical line: three-body Efimovian threshold $\alpha = \alpha_c^{2,1}$. Blue dashed line: cubic fit in the variable $Y = (\alpha_c^{2,1} - \alpha)^{1/2}$ on the interval $10 \le \alpha \le 13.5$, i.e. $B_{2,2}(0^+) = -0.73742 + 0.061673Y + 0.074061Y^2 - 0.0084386Y^3$. Red dashed line: same on the interval $12.5 \le \alpha \le 13.5$, i.e. $B_{2,2}(0^+) = -0.73761 + 0.062206Y + 0.073901Y^2 - 0.0086514Y^3$.</figcaption></figure>

<table><caption>**Table 4** – Numerical values of the fourth cluster coefficient $B_{2,2}(0^+)$ of the unitary Fermi gas in an infinitesimal-stiffness trap, tabulated as functions of the mass ratio $\alpha = m_1/m_\downarrow$. Since the coefficient is invariant by change of $\alpha$ to $1/\alpha$, we restricted ourselves to $\alpha \ge 1$. The uncertainties, not specified, are less than one percent. These values are shown graphically in Figure 6.</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>2,2</sub></td><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><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>2,2</sub></td><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><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.4</td></tr><tr><td>B<sub>2,2</sub></td><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>

## Acknowledgements

This work was supported by Japanese "Grants-in-Aid for Scientific Research" KAKENHI Grant Numbers 21H00116 and 22K03492.
# Appendix A. Asymptotic approximant of $I_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$, its sum over $\ell$ and $\epsilon$, its dominant behavior

## A.1. The asymptotic approximant

To obtain an approximation $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ at large $\ell$ of the contribution $I_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ of angular momentum $\ell$ and parity $\epsilon$ to the quantity $I_{3,1}(\bar{\omega})$ of equation (4), we use the notations of Section 2 and start from expansion (19). A clever calculation of the trace on $x > 0$ allows us to collect contributions making $x$ and $-x$ appear and to reduce to integrals on the whole real axis, in which we only have to take the traces on the variables $u$ and $m_z$:

$$ \begin{aligned} \ln \det M_{3,1}^{(\ell, \epsilon)} ={}& \text{const} - \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 \end{aligned} \quad (45) $$

where the rotation operator of axis Ox of angle $\pi$ comes from Footnote 9 and we have introduced the primitive kernel (16) divided on the left by the diagonal part (15) and taken with zero scale exponent, i.e. $\bar{K}_{3,1}^{(\ell)} \equiv \mathcal{D}_{3,1}^{-1} K_{3,1}^{(\ell)} (iS = 0)$. The contributions kept to the second side of equation (45) constitute an even, regular, real-valued and rapidly decreasing function of S, that is $-\phi^{(\ell, \epsilon)}(S)$. So we can integrate by parts over S in equation (4) and we recognize the Fourier component of $-\phi^{(\ell, \epsilon)}(S)$ at frequency $\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) $$

which is easy to obtain from (45) by means of the identity in the sense of distributions $\int_{\mathbb{R}} \exp(ikS) dS = 2\pi\delta(k)$. A long but not difficult calculation, treating the contributions to $\phi^{(\ell, \epsilon)}(S)$ in the order they appear, finally gives:

$$ \begin{aligned} 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} \\ & \qquad \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) \to (\lambda_{\bar{\omega}}^{-1}, \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} \\ & \qquad \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) \to (\lambda_{\pm}^{-1}, \phi')] } \end{aligned} \quad (47) $$
with the notations $X_{\pm} = X \pm \bar{\omega}/2$, $\lambda_{\bar{\omega}} = \exp(\bar{\omega})$, $\lambda_{\pm} = \exp(X_{\pm})$ (these last two modeled on $\lambda = \exp(x)$) completing those of equation (17) and the angular functions$^{21}$

$$ \begin{gathered} \mathcal{F}^{(\ell, \epsilon)}(\theta; \phi) \equiv \mathrm{Tr}_{m_z}^{\epsilon} \left( e^{-i\theta L_z/\hbar} e^{i\phi L_x/\hbar} \right) \\ \text{and} \quad \mathcal{F}^{(\ell, \epsilon)}(\theta, \theta'; \phi, \phi') \equiv \mathrm{Tr}_{m_z}^{\epsilon} \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) \end{gathered} \quad (48) $$

The sums over parity $\epsilon$ of the quantities defined in (48) have simple expressions in terms of angles $\xi \in [0, \pi]$, see reference [40],$^{22}$

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

$$ \mathcal{F}^{(\ell)}(\theta; \phi) = \frac{\sin[(\ell + 1/2)\xi]}{\sin(\xi/2)} \quad \text{with} \quad (50) $$

$$ 1 + 2 \cos \xi = uu'(1 + \cos \phi \cos \phi') - (u+u') \sin \phi \sin \phi' + \nu \nu' (\cos \phi + \cos \phi') + \cos \phi \cos \phi' $$

which also shows how to return to fixed parity, for example

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

However, we prefer to reserve equations ((49), (50), (51)) to analytical studies; for numerical computation, we evaluate the traces of (48) in the eigenbasis of $L_z$ (under the constraint $(-1)^{m_z} = \epsilon$) after insertion of closure relations in the eigenbasis of $L_x$ (obtained by numerical diagonalization of its tridiagonal matrix in the $|\ell, m_z\rangle$ basis) at the location of the $\phi$- or $\phi'$-angle rotation operators, then we use the value of the integral$^{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. Its sum on $\ell$ and $\epsilon$

The sum of $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ on the parity $\epsilon$ is straightforwardly done, see equations ((49), (50)). The sum on $\ell \in \mathbb{N}$ gives rise, in the sense of distributions, to the Fourier series of a Dirac comb$^{24}$ and is calculated thanks to the identity:

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

It thus remains to find the cases of cancellation of the angle $\xi$, knowing that the polar angles $\theta, \theta'$ (from now on taken as integration variables in preference to $u, u'$) span $[0, \pi]$ and that the

$$ \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>$^{21}$ One may object that the definition of $\mathcal{F}^{(\ell, \epsilon)}(\theta, \theta'; \phi, \phi')$ should in principle include an orthogonal projector on the subspace of parity $(-1)^{L_z/\hbar} = \epsilon$ next to the operator $e^{i\theta' L_z/\hbar}$. The parity in $\phi$ and $\phi'$ of the denominator in the third contribution to (47), however, allows us to do without it. For example, only the even part of $e^{i\phi L_x/\hbar}$ contributes after integration over $\phi$; it indeed preserves the parity of $L_z/\hbar$.</footnote>

<footnote>$^{22}$ Turning to the half-angles, we note that we have more simply $\cos(\xi/2) = \cos(\theta/2)\cos(\phi/2)$ in equation (49).</footnote>

<footnote>$^{23}$ For optimization, we can (i) replace $\exp[i(-m_z\theta + m'_z\theta')]$ by its real part $\cos(m_z\theta)\cos(m'_z\theta') + \sin(m_z\theta)\sin(m'_z\theta')$ and tabulate the corresponding sines and cosines, (ii) reduce to a single index loop $m_x = m'_x$ and use the symmetry $m_x \leftrightarrow -m_x$, (iii) use the symmetry $(m_z, m'_z) \leftrightarrow (-m_z, -m'_z)$ (we restrict ourselves to $(-1)^{m'_z} = \epsilon$, see Footnote 21), (iv) tabulate the powers of degree $|m_x|$ appearing in (52).</footnote>

<footnote>$^{24}$ One simply writes $\sin[(\ell + 1/2)\xi] = \{\cos(\ell\xi) - \cos[(\ell + 1)\xi]\}/[2\sin(\xi/2)]$ to introduce the partially telescopic sum</footnote>
azimuthal angles $\phi, \phi'$ can be conveniently taken in $[-\pi, \pi]$. In equation (49), we find as the only nodal point $(\theta, \phi) = (0, 0)$; in the neighborhood of this point, $\xi^2 \approx \theta^2 + \phi^2$ so that the action of $\delta(\xi)$ is easily evaluated in polar coordinates $(\rho, \gamma)$ in the plane $(\theta, \phi)$, for any regular function $f$:

$$ \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 \int_{-\pi/2}^{\pi/2} d\gamma \rho \cos\gamma \frac{\delta(\rho)}{2(\rho/2)^2} = 2f(0,0) \quad (54) $$

In equation (50), we find three nodal lines<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) $$

In the integral we are dealing with, which involves a regular function $g$ of the four 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) $$

the lines $l_1$ and $l_2$ have a zero contribution because of the Jacobians $\sin\theta$ and $\sin\theta'$. In the vicinity of line $l_3$, we quadratize the quantity $\xi$ into the variables $\phi, \phi'$ and $\theta' - \theta$; the integral at fixed $\theta$

$$ \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) $$

is well computed in the eigenbasis of the quadratic form appearing in the denominator and under the square root (it has eigenvalues 1 and $1 \pm \cos\theta$). It remains simply

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

from which the expression of the desired sum (respecting the order of the contributions in equation (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 \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 \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}} \rightarrow \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 \rightarrow \lambda_\pm^{-1}]} \end{aligned} \quad (59) $$

In practice, the most important case corresponds to $\bar{\omega} = 0^+$, because it gives access to the cluster coefficient $b_{3,1}$ of the homogeneous gas. We give the corresponding explicit expression of (59),

<footnote n="25">²⁵If $\xi = 0$, $\mathcal{T}^{(1)}(\theta, \theta'; \phi, \phi') = 3$; however, this is the trace of a rotation matrix in $\mathbb{R}^3$ which must therefore reduce to the identity: we have $\mathcal{R}_z(\theta)\mathcal{R}_x(-\phi)\mathcal{R}_z(-\theta')\mathcal{R}_x(-\phi') = \mathbf{1}$ where $\mathcal{R}_x(\theta)$ is the rotation of angle $\theta$ of axis $Ox$, etc. By taking the $zx$ matrix element of this relation, we get the condition $\sin\theta' \sin\phi = 0$. Similarly, by changing the order of the operators by circular permutation under the trace, we end up with $\sin\theta \sin\phi' = 0$. The nodal lines $l_1$ and $l_2$ in (55) correspond to $\sin\theta = \sin\theta' = 0$ that is $\theta = \theta' = 0[\pi]$, the cases $\theta = 0, \theta' = \pi$ and $\theta = \pi, \theta' = 0$ being trivially not suitable. The line $l_3$ corresponds to $\sin\phi = \sin\phi' = 0$ i.e. $\phi = \phi' = 0[\pi]$, which in practice reduces to $\phi = \phi' = 0$; indeed, the cases $\phi = 0, \phi' = \pi$ and $\phi = \pi, \phi' = 0$ are obviously not suitable, and the case $\phi = \phi' = \pi$ imposes $\cos(\theta + \theta') = 1$ so $\theta = \theta' = 0$ or $\pi$, which corresponds to nodal points, of zero contribution to the integral on $\theta, \theta', \phi, \phi'$. The remaining cases $\sin\theta = \sin\phi = 0$ or $\sin\theta' = \sin\phi' = 0$ lead only to nodal points and do not contribute either.</footnote>

<footnote n="26">²⁶ Where $\phi + \pi$ appears in the argument of $\mathcal{T}^{(\ell)}$, we perform the change of variable $\phi \rightarrow \phi - \pi$.</footnote>
always respecting the order of the contributions, in terms of trigonometric functions and the dilogarithm or Bose function $g_2(z) = \sum_{n>0} z^n / n^2$:

$$ 
\begin{aligned}
& \sum_{\ell \in \mathbb{N}} \sum_{\epsilon} J_{3,1}^{(\ell, \epsilon)} (\bar{\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) $$

where $\zeta = (2\alpha + i\sqrt{2(1+3\alpha)})/(1+\alpha)$.

### A.3. Its dominant behavior

To obtain an asymptotic equivalent of $J_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ or, what amounts to the same thing, of $I_{3,1}^{(\ell, \epsilon)}(\bar{\omega})$ for large angular momentum, it suffices to restrict ourselves to the first contribution to the second side of equation (47), which comes from the term of order one in the kernel $\mathcal{K}_{3,1}^{(\ell)}$ in expansion (19). We then write it as an integral over the angle $\xi$ by inserting a Dirac delta $\delta(\xi - \xi(u, \phi))$ linking $\xi$ to $u$ and $\phi$ as in equation (49), then explicitly calculating the integral over $\phi$ and then over $u$.<fnref n="27" /> To simplify, we restrict ourselves to $\bar{\omega} = 0^+$ and find:

$$ I_{3,1}^{(\ell, \epsilon)}(0^+) \underset{\ell \to +\infty}{\sim} \int_0^\pi \mathrm{d}\xi \frac{\sin[(\ell + 1/2)\xi]}{2\sin(\xi/2)} [\rho_0(\xi) + \epsilon\rho_1(\xi)] \quad \text{where} \quad \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) $$

with the notation $\beta = 1/\alpha$, the prefactors and the auxiliary functions<fnref n="28" />

$$ 
\begin{aligned}
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{aligned}
\quad (62) $$

$$ 
\begin{aligned}
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{aligned}
\quad (63) $$

$$ 
\begin{aligned}
F_1(Z) &= Z \left( \frac{\pi}{2} - \arctan Z \right) \\
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] [3\beta + (1-2Z)^2] \right\}^{1/2}}
\end{aligned}
\quad (64) $$

<footnote n="27">²⁷After having reduced to an integral on $\phi \in [0, \pi]$ by periodicity and parity of the integrand, we use the relation $\delta(2\cos\xi + 1 - u(1 + \cos\phi) - \cos\phi) = \delta(\phi - \phi_0)/[(1+u)\sin\phi_0]$ with $\sin\phi_0 = 2(u - \cos\xi)^{1/2}(1 + \cos\xi)^{1/2}/(1+u)$ where the root $\phi_0$ is in $[0, \pi]$ if and only if $0 \le \theta \le \xi$ given Footnote 22. We also use relation (51).</footnote>

<footnote n="28">²⁸When we move to the complex plane, we mean that $\arctan Z = \ln[(1+iZ)/(1-iZ)]/(2i)$. As a result, the branch cut of the function $F_0(Z)$ is $]-\infty, -1]$ and that of $F_1(Z)$ is $]-\infty, -1] \cup i[1, +\infty]$. Note that, despite appearances, $F_0(Z)$ is analytic in the unit disk, as shown by its series expansion at $Z = 0$. Note also that with this definition of the arctan function, we have $\pi/2 - \arctan(ix) = -i\argh(1/x) \forall x \in \mathbb{R} \setminus [-1, 1]$, which is already used in the integral (61) on the real axis because the quantity under the last square root in the denominator of (64) can become negative; we have made here the choices $\ln(-1) = i\pi$ and $\sqrt{-1} = i$.</footnote>
To reduce to contour integrals on the unit circle C, of the form

$$ \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) $$

we proceed in three steps: (i) we eliminate the factors $\sin(\xi/2)$ in the integrand of (61), first by simplifying the sine in the denominator in front of the brackets with a sine in the numerator of $\rho_n(\xi)$, then by using the remaining sine to make $\cos[(\ell+1)\xi] - \cos(\ell\xi)$ appear as in Footnote 24; (ii) we extend the integral on $\xi$ to $[-\pi, \pi]$ by parity of the integrand; (iii) we take as new integration variable $z = \exp(i\xi)$, which spans the unit circle, so that $\cos\xi = (z+1/z)/2$, and we can write $\cos[(\ell+1)\xi] - \cos(\ell\xi) = \text{Re}(z^{\ell+1} - z^\ell)$ and take out the function $\text{Re}$ of the integral since the remainder of the integrand is real-valued on $C$. For the usual determinations of the logarithm and the square root in the complex plane (branch cut on $\mathbb{R}^{-}$), we find that the integrand of $\mathcal{I}_0$ has as singularity in the unit disk a trident-shaped branch cut, union of a segment $OO'$ and an arc $AA'$, see Figure 7a, and that of $\mathcal{I}_1$ a flower-shaped branch cut, union of three arcs $OB$, $OB'$ and $BB'$, see Figure 7b.<fnref n="29" /> The point $O$ is here the origin of the coordinates, $A$ and $A'$ are the points of affixes $z_0$ and $z_0^*$, $B$ and $B'$ are the points of affixes $z_1 = -z_0^*$ and $z_1^*$, and $z_0$ is the unique solution of the equation $z + z^{-1} + 1 + i\sqrt{3\beta} = 0$ in the unit disk (so that $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) $$

Using Cauchy's integral theorem, we shrink the integration contour until it fits the branch cuts, without changing the value of $\mathcal{I}_n$. In the limit $\ell \to +\infty$, because of the factor $z^\ell$, the integral is then dominated by the neighborhood of the singularity points farthest from $O$, namely $A$ and $A'$ for $\mathcal{I}_0$, $B$ and $B'$ for $\mathcal{I}_1$, this being true for any mass ratio $\alpha$. $\mapsto$ For the computation of $\mathcal{I}_0$ in the neighborhood of $z_0$, let $z = z_0(1 + \eta u)$, where $u$ is a complex number and $\eta > 0$ is an infinitesimal, and approximate the horn of the trident by its right semitangent at the vertex $A$; we then have the equivalents $f_0((z+1/z)/2) \sim C_0/(\eta u)$ and $z^\ell \sim z_0^\ell \exp(\ell\eta u)$, and the semitangent becomes the half line of origin $O$ of direction $-C_0$ in the $u$ space (considering Footnote 28), that is the half line $\mathbb{R}^-$ in the space of $v = \eta u/C_0$, so that

$$ \mathcal{I}_0^\mathbb{R} \xrightarrow{\ell \to +\infty} \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) $$

where the paths $\mathcal{C}_u$ and $\mathcal{C}_v$ surround the half-tangents counterclockwise in $u$ or $v$ space. The expression for $C_0$ is given in equation (21). The integral in the last side of (67) is easily computed:

$$ \begin{aligned} \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] \\ &= \int_{-1}^{0} \frac{dx}{x} e^{\ell x C_0} (-i\pi|x|^{1/2}) \xrightarrow{\ell \to +\infty} \frac{i\pi^{3/2}}{\sqrt{\ell C_0}} \end{aligned} \quad (68) $$

The horn $A'$ of affix $z_0^*$ gives a complex conjugate contribution of the horn $A$. We get the first part of equation (20). $\hookrightarrow$ For the calculation of $\mathcal{I}_1$ in the neighborhood of $z_1$, we proceed in the same way, by setting $z = z_1(1 + \eta u)$, $\eta \to 0^+$. Now the local behavior $f_1((z+1/z)/2) \sim A/(\beta\eta u)^{1/2}$ has itself a branch cut because of the square root and is characterized by two complex amplitudes $A$

<footnote n="29">²⁹In this second case, we expect to have three branch cuts, the one $L_1$ coming from the function $F_1(Z)$, the ones $L_3$ and $L_2$ coming from $\sqrt{1-Z}$ and the other square root in the denominator of the function $f_1(Z)$. In reality, $L_3$ is included in $L_1$. Moreover, the intersections $L_1 \cap L_2$ and $L_1 \cap L_3$ are not branch cuts of the integrand, for the reason that $z \to z^{1/2}(\pi/2 - \arctan z^{1/2})$ has only $[-1, 0]$ as a branch cut (for $z \in [-1, 0]$), the change of sign of $\pi/2$-arctan compensates for that of the square root). The arcs $OB$ and $OB'$ are made up of the points lying only on $L_2$, and the arc $BB'$ collects the points belonging only to $L_1$.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 7** – Branch cuts of the integrand of equation (65) in the unit disk (thick red line) in the case (a) $n = 0$ and (b) $n = 1$, obtained numerically for a mass ratio $\alpha = 1$ taken as example.</figcaption></figure>

and $\mathcal{B}$, hence the structure with two branch cuts, the one coming from $f_1$ and the one coming from arctan $f_1$ in $F_1(f_1)$.<fnref n="30" /> They are approximated by two semitangents corresponding in the space of $v = \eta \mathcal{B} u / \mathcal{A}^2$ simply to $\mathcal{A}^{-2} \mathbb{R}^-$ and, in light of Footnote 28, to $\mathbb{R}^-$. One has

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

where the path $\mathcal{C}_v$ surrounds (in the third quadrant) the half line $\mathcal{A}^{-2} \mathbb{R}^-$ counterclockwise and then vertically joins the half line $\mathbb{R}^-$ from below to surround it also counterclockwise. The calculation leads to<fnref n="31" />

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

which gives the second part of equation (20), where the coefficient $\mathcal{B}$ is called $C_1$ and given by equation (21).

## Appendix B. Asymptotic approximant of $I_{2,2}^{(\ell, \varepsilon)}(0^+)$, its sum over $\ell$ and $\varepsilon$

<footnote n="30">³⁰In the previous case $n = 0$, this does not occur: In the function $F_0(z)$, on either side of the branch cut $\mathbb{R}^-$ of the square root, $\sqrt{z}$ can of course take two opposite values but arctan is an odd function. As a consequence, the branch cut of $F_0(z)$ is imposed by the function arctan and is $| -\infty, -1]$, in agreement with Footnote 28.</footnote>

<footnote n="31">³¹One indeed has $\sqrt{\mathcal{A}^2 / \mathcal{B}} = -\mathcal{A} / \sqrt{\mathcal{B}}$ for any mass ratio $\alpha$.</footnote>
## B.1. The asymptotic approximant

The aim here is to write as explicitly as possible an approximation $J_{2,2}^{(\ell, \epsilon)}(0^+)$ at large $\ell$ of the contribution $I_{2,2}^{(\ell, \epsilon)}(0^+)$ of angular momentum $\ell$ and parity $\epsilon$ to the quantity $I_{2,2}(0^+)$ of equation (4). We take as a starting point expansions ((38),(42)) according to the value of $\epsilon$. Their insertion into the integrand of (4) leads, after integration on $S$, to the structure

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

where the subscript indicates that it is a cross contribution of the integral kernels $K_i$ and $K_j$. The case $(i, j) = (1, 2)$ was treated in reference [40] for $\epsilon = (-1)^{\ell-1}$, using forms (36) of the kernels; here we simply extend it to the case $\epsilon = (-1)^{\ell}$, which is just a rewriting, and we copy the result:

$$\mathcal{J}_{K_1, K_2}^{(\ell, \epsilon)} = \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] [\text{ch} x \mathcal{D}_{2,2}(x, u) \mathcal{D}_{2,2}(x, u')]^{-1} \mathcal{F}^{(\ell, \epsilon)}(\theta, \theta'; \phi, \phi')}{\{1 + \frac{1}{1+\alpha} [(u+e^{-x})(u'+e^{-x}) + \nu\nu' \cos\phi]\} \{1 + \frac{\alpha}{1+\alpha} [(u+e^x)(u'+e^x) + \nu\nu' \cos\phi']\}} \quad (72)$$

with notations (17) and (48). The case $(i, j) = (1, 3)$ is much simpler to handle (especially numerically) by using forms (96) of the kernels, marked by a Czech accent. Only one factor depends on the scaling exponent $iS$, on which we have to integrate in (4). In the notations of (96) it comes after an integration by parts:

$$\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)$$

We decide to integrate the Dirac distribution on $u'$ at fixed $u, x, x'$, which links the value of $u'$ to these other variables:

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

As $u'$ must be between -1 and 1, the integration interval on $u$ is constrained and reduces to $[u_{\min}, u_{\max}]$ with

$$\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)$$

To compute the resummed kernel (41), we recognize in the external part of $\tilde{K}_3$ (last contribution of (96)) a form factorized into an operator on $x$-space and an orthogonal projector of rank one on $u$-space, which reduces the problem to the inversion of an operator on $x$-space only: 32

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

with

$$\langle u | \chi \rangle = \frac{1}{\sqrt{2}} \quad \text{and} \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)$$

<footnote>32There is no gauge transform of type (39) to be made here because $\tilde{K}_3$ is directly independent of the scale exponent. We also notice in equation (95) that $\hat{\mathcal{D}}_{22} \equiv 1$.</footnote>
We finally obtain

$$ \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) $$

where we have introduced the symmetric function of $x$ and $x'$ taking into account (75) by a Heaviside function $Y$:

$$ \mathcal{R}(x, x') = Y \left( \sqrt{2\alpha} - |e^{x'} - e^x| \right) \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) $$

The complex numbers $z_i$ and $z'_i$ and their arguments $\tau_i$ and $\tau'_i$ are defined below equation (94) as functions of the angles $\theta$ and $\theta'$ in the interval $[0, \pi]$ such that $u = \cos\theta$ and $u' = \cos\theta'$, where the angle $\theta'$ is related to $\theta$ by equation (74); the writing of the denominator of the integrand in (79) takes into account the equality of the moduli $|z'_i| = |z_1|$.33 Finally, the last contribution in (71) is deduced from the second contribution by changing everywhere $\alpha$ into its inverse $1/\alpha$:

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

## B.2. Its sum on $\ell$ and $\epsilon$

It remains to compute the sum on all channels $(\ell, \epsilon)$ of the asymptotic approximant (71). The contribution of type $K_1 K_2$ is treated exactly as in Appendix A:

<footnote n="33">33 In practice, the angular integral in (79) is done numerically on the angle $\theta$ (rather than on $u$) with the 41-point Gauss–Legendre method; if $x$ is close enough to but different from the singularity point $x_{\text{sing}} = \ln\sqrt{\alpha/2}$ where the first denominator $|z_1|$ can vanish, i.e. $0 < |x - x_{\text{sing}}| < \theta_c$ with $\theta_c = 0.15$ for example, a narrow structure in $\theta$ appears and the integration on $[\theta_{\min}, \theta_c]$ is performed with the change of variable $\theta = |x - x_{\text{sing}}| \sin t$, with the midpoint rule on the variable $t$ discretized in 100 equidistant points. On the other hand, the second denominator in (79) always remains greater than $1 + \min(\alpha, \alpha^{-1})$ on the integration domain and cannot vanish. The integration on $x$ and $x'$ is done with the midpoint method on a truncation interval $[x_{\min}, x_{\max}]$, by arranging that $x_{\text{sing}}$ is the edge of one of the subintervals of width $dx$; we make the fixed choice $x_{\min} = -5$ since $\mathcal{R}(x, x') = O(\min(\exp(x), \exp(x'))$ when $x, x' \to -\infty$, but we extrapolate to $x_{\max} = +\infty$ cubically in $1/x_{\max}$ from the four values $x_{\max}^{(0)} = 7$, $x_{\max}^{(1)} = 30$, $x_{\max}^{(2)} = 52$ and $x_{\max}^{(3)} = 75$. To take advantage of the fact that $\mathcal{R}(x, x') = O(\exp(-2x))$ on an increasingly narrow support $x - x' = O(\exp(-x))$ when $x \to +\infty$, we separate the matrix $A$, discretized version of the operator $\mathbf{1} + \hat{k}$, in blocks $A_{ij}$ between the intervals number one $[x_{\min}, x_{\max}^{(0)}]$ and number two $[x_{\max}^{(0)}, x_{\max}^{(n)}]$, neglecting the blocks 12, 21 and 22 of its inverse $A^{-1}$ and we then compute its block 11 by the expression $(A^{-1})_{11} = [A_{11} - A_{12}(A_{22})^{-1}A_{21}]^{-1}$ which includes the effect of the non-diagonal coupling in the form of a second-order effective Hamiltonian (it would be incorrect to neglect it completely because $\langle x|\hat{k}|x' \rangle$ does not tend to zero when $x, x' \to +\infty$ at $x - x'$ fixed). As the function $\mathcal{R}(x, x')$ shows near $x = x' = x_{\text{sing}}$ on the left more and more rapid oscillations when $\ell$ increases, we cannot take a fixed value of the step $dx$ but we have to use an iterative method of Romberg type: From the initial choice $dx = 1/20$, we reduce the step $dx$ by a factor of 2 and extrapolate at zero step linearly in $(dx)^2$ until the extrapolated value is stable at the percent level (or $dx$ falls below the very small value $2 \times 10^{-3}$). The matrix element in $|\ell, m_x = 0\rangle$ is computed as above equation (52) with the same symmetry tricks as in Footnote 23 and the integral over $\phi$ is deduced from equation (52).</footnote>
$$
\begin{align}
& \sum_{\ell \in \mathbb{N}} \sum_{\epsilon} \mathcal{J}_{K_1, K_2}^{(\ell, \epsilon)} \nonumber \\
&= \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})]} \nonumber \\
&= \frac{(1+\alpha)^2}{8\pi^2\alpha} \left\{ \begin{aligned}
& 2\pi \operatorname{arcsin} \left( \frac{\alpha}{(1+2\alpha)(\alpha+2)} \right)^{1/2} \\
& - \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] \\
& - \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]
\end{aligned} \right\} \tag{81}
\end{align}
$$

with $z = \sqrt{\alpha} + i\sqrt{1+\alpha}$, $z' = \sqrt{\alpha^{-1}} + i\sqrt{1+\alpha^{-1}}$ and $g_2$ the dilogarithm function. In the $K_1 K_3$ type contribution, let us first sum over the angular momentum $\ell$ of fixed parity $(-1)^\ell = \eta$, so that the operator $\hat{k}$ in (78) takes the fixed value $\hat{k}_\eta$. Let us transform the quantum average in state $|\ell, m_x = 0\rangle$ in the numerator of the integrand of (79) by inserting a closure relation in the eigenbasis of $L_x$ and using the expression of the corresponding matrix elements deduced from [53, equations (7.2–9) on page 101, (7.3–15) on page 105 and (7.4–7) on 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)
$$

where $Y_\ell^m$ are the usual spherical harmonics. It then leads to<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)
$$

It remains to invoke the closure relation (8.6–10) on page 146 of reference [53] on spherical
harmonics and the spatial parity property $Y_{\ell}^{m}(\pi - \theta, \phi + \pi) = (-1)^{\ell} Y_{\ell}^{m}(\theta, \phi)$ to obtain the closure
relation with fixed $\ell$ parity:<fnref n="36" />

<footnote n="34">³⁴The relation (82) is used twice, in its direct form with $(\tau, \psi) = (\tau_1, \phi)$ and in its conjugated form with $(\tau, \psi) = (\tau'_1, 0)$. If $\tau_1$ is in $[-\pi, 0]$, the relation does not apply, but it is then sufficient to change $\tau_1$ into $-\tau_1$, which amounts to changing the integration variable $\phi$ into $\phi+\pi$ in (79) (indeed, $\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)$) and does not modify the value of the integral. We proceed in the same way if $\tau'_1 \in [-\pi, 0]$. In the following, we can therefore assume that $\tau_1$ and $\tau'_1$ are in the interval $[0, \pi]$.</footnote>

<footnote n="35">³⁵Using the addition theorem for spherical harmonics, see [53, equation (8.6–3) on page 145], we show that the second side of (83) is also written $(2\ell + 1)P_\ell(\cos\delta)$ where $P_\ell$ is a Legendre polynomial and $\delta$ is the angle between the unit vectors of polar coordinates $(\tau_1, \phi)$ and $(\tau'_1, 0)$; its cosine $\cos\delta = \cos\tau_1 \cos\tau'_1 + \sin\tau_1 \sin\tau'_1 \cos\phi$ appears in the denominator of the integrand of (79).</footnote>

<footnote n="36">³⁶ In the Dirac distributions $\delta(\phi - \phi_0)$, $\phi$ has a meaning modulo $2\pi$; thus, one can replace the arbitrary integration interval $[0, 2\pi]$ of equation (79) by the interval of length $2\pi$ centered on $\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{where } u_0 \equiv \frac{2\alpha + e^{2x} - e^{2x'}}{2\sqrt{2\alpha}e^x} \nonumber
\end{align}
$$

where we have replaced in the third side $\cos \tau_1 \pm \cos \tau'_1$ by its value, remembering that $|z_1| = |z'_1|$ and using (74). The integration on $\phi$ is straightforward in view of Footnote 36. To integrate over $u$, we need to know if the root $u_0$ is in the interval $[u_{\min}, u_{\max}]$. For this purpose, we divide the support of the Heaviside function in (79) into four distinct areas:

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

$$ \exp(x') < \exp(x) < \exp(x') + \sqrt{2\alpha} \quad \text{and} \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'). \qquad (iv) $$

If the expressions depending on $x$ in the definitions (75) are denoted by $u_{\min}^{\exp}$ and $u_{\max}^{\exp}$, we find that we systematically have $u_{\min} = -1 < u_0 < u_{\max} = u_{\max}^{\exp}$ in zone (i), $u_{\min} = u_{\min}^{\exp} < u_0 < u_{\max} = 1$ in zone (ii), $u_{\min} = -1 < u_{\max} = 1 < u_0$ in zone (iii), $u_{\min} = u_{\min}^{\exp} < u_{\max} = u_{\max}^{\exp} < u_0$ in zone (iv). In other words, the integral of $\delta(u - u_0)$ over $u$ is always equal to one in the first two areas and to zero in the last two. We deduce the sum of the quantity $\mathcal{R}(x, x')$ on all $\ell$ of fixed parity:

$$
S_{\eta}(x, x') = \frac{1}{2} \left[ \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') \right. \\
\left. \qquad + \eta \frac{Y(\sqrt{2\alpha} - |e^x - e^{x'}|) Y(e^x + e^{x'} - \sqrt{2\alpha})}{e^{2x} + e^{2x'} + 1} \right] \quad (85)
$$

It remains to sum on $\eta = \pm 1$ to arrive at the desired result:

$$ \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) $$

Again, the operator $\hat{k}_{\eta}$ is deduced from equation (77) by replacing in the second side $(-1)^{\ell+1}$ by $(-\eta)$; the numerical inversion of the operators $\mathbf{1} + \hat{k}_{\eta}$ and the integration over $x$ and $x'$ are done with the same techniques and tricks as in Footnote 33 (on the other hand, there is no more integration to do on $u$). Finally, as shown in equation (80), $\sum_{\ell \in \mathbb{N}} \mathcal{J}_{K_2, K_3}^{(\ell)}$ is deduced from expression (86) by changing everywhere $\alpha$ to $1/\alpha$ (also in $S_{\eta}(x, x')$).

## Appendix C. The operator $M_{2,2}^{(\ell, \varepsilon)}(\mathrm{i}S)$ in the formulation of reference [52]

To obtain the operator $M_{2,2}^{(\ell, \varepsilon)}(\mathrm{i}S)$ at the basis of conjecture (4),(5)) on the cluster coefficient $B_{2,2}(\hat{\omega})$ of the trapped system, we first write a Faddeev ansatz for an eigenstate of the unitary $2+2$-body $\uparrow\uparrow\uparrow\uparrow$ problem of zero energy and zero momentum in free space. This ansatz is expressed in terms of an unknown function $D$ of two wave vectors. Taking into account the Wigner-Bethe- Peierls two-body $\uparrow\downarrow$ contact conditions (as in Footnote 1) leads to an integral equation for this
function. Then, we use the rotational invariance to project the equation on the subspace of
angular momentum $\ell$ and parity $\epsilon$ as reference [54] explains it in detail. Finally, we use the
scale invariance of the unitary problem (in the absence of three-body Efimov effect) to choose
a function $D$ with a well-defined scale exponent $s$. The integral equation is then reduced to the
condition $\det M_{2,2}^{(\ell,\epsilon)}(s) = 0$ where $M_{2,2}^{(\ell,\epsilon)}(s)$ is a kernel operator, which we have to extend on the
pure imaginary axis $s = iS$ to evaluate expression (4).

Expression ((34), (35), (36)) of $M_{2,2}^{(\ell,\epsilon)}$ (iS) corresponds to the choice of unknown function $D(\mathbf{k}_2, \mathbf{k}_4)$ where $\mathbf{k}_2$ and $\mathbf{k}_4$ are the wave vectors of two opposite-spin fermions, as in reference [40] (the wave vectors $\mathbf{k}_1$ and $\mathbf{k}_3$ disappear in the limit expressing the contact condition). The starting integral equation is given by [40, equation (13)]. Another choice is made in reference [52], corresponding to the change of function

$$
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)
$$

It has the advantage of providing a much simpler expression of the singular integral kernel $K_3$, i.e. of the third contribution in (36) (but not of the kernels $K_1$ and $K_2$, which justifies in the end keeping (36) in the numerical calculation of (4)). Up to a sign, the new variables are simply the wave vector of the center of mass and the relative motion of particles 2 and 4. In the following, we use the notation $\alpha = m_\uparrow / m_\downarrow$.

Let's follow the previously stated steps. The starting integral equation is written in parameter-
ization (87) [52]

$$
\frac{\kappa}{4\pi} F(\mathbf{u}, \mathbf{v}) + \int \frac{d^3 \mathbf{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 \mathbf{v}'}{(2\pi)^3} \frac{F(-\mathbf{u}, -\mathbf{v}')}{\kappa^2 + \nu'^2} = 0 \quad (88)
$$

with the notations taken from reference [52],<fnref n="37" />

$$
\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}
$$

by correcting what seems to us to be a sign error in the coefficient of the vector **k** (here called **u**')
in the definition of **v**<sub>23</sub> (third row and first column of Table III of this reference) and in the
expression of *E*<sub>coll</sub> (notation not introduced here) just above equation (141) of this reference. As
in reference [40], to make future transformations simpler, we adopt a variational formulation of
equation (88), δ*E* / δ*F*<sup>*</sup>(**u**, **v**) = 0 with the functional of *F* and *F*<sup>*</sup> that follows:

$$
\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 + \nu'^2} \right] \quad (90)
$$

where the parity $\epsilon = \pm 1$ of the solution has been introduced to make the sign – disappear in
front of **u** and **v**' in $F(-\mathbf{u}, -\mathbf{v}')$. Let us now consider rotational invariance, restricting ourselves
to the subspace of total angular momentum $\ell \in \mathbb{N}$, with a zero angular momentum component
along the quantization axis Oz. The solution $F(\mathbf{u}, \mathbf{v})$ is then expressed in terms of $2\ell + 1$ functions
$f_{m_z}^{(\ell)}(u, v, w)$ ($-\ell \le m_z \le \ell$) of only three real variables, the moduli u and v of the two vectors and
the non-oriented angle $\theta = \langle \mathbf{u}, \mathbf{v} \rangle \in [0, \pi]$ between them or, what amounts to the same thing, its
cosine $w = \cos\theta$; to fix the parity at $\epsilon$ is to impose $(-1)^{m_z} = \epsilon$ thus to decouple the problem into
$\ell$ and $\ell + 1$ unknown functions, which we indicate by an exponent $\epsilon$ on the sum sign below. Our

<footnote n="37"><sup>37</sup>Reference [52] arranges the fermions in the order ↑↓↑↓. We renumber the particles accordingly.</footnote>
ansatz for $F$ is that of [40, equation (14)]. We insert it into functional (90) and integrate over the variables other than the arguments of the functions $f_{m_z}^{(\ell)}$ in the same way as in that reference. For example, let us apply [40, equation (45)] to the last contribution of functional $\mathcal{E}$, the one with $\epsilon$ in factor, that gives rise to the singular kernel $K_3$. First, a direct reference trihedron $\mathcal{T}$ with polar axis $\mathbf{u}$ and another $\mathcal{T}'$ with polar axis $\mathbf{u}'$ are chosen. The integration on $\mathbf{v}$ (or on $\mathbf{v}'$) is performed in the spherical coordinates associated to $\mathcal{T}$ (or to $\mathcal{T}'$), the cosine of the corresponding polar angle being $w$ (or $w'$). Integrating over the other variables then amounts to taking the average over the orientations of $\mathcal{T}$ and $\mathcal{T}'$ which can be done by fixing the direction of the vector $\mathbf{u}$ along the convenient direction $\mathcal{O}\mathbf{x}$ and replacing the integration over the direction of $\mathbf{u}'$ by an integration over rotation $\mathcal{R}$ mapping $\mathcal{T}$ to $\mathcal{T}'$ in the SO(3) group with an invariant measure, explicit $d\mathcal{R} = d\mathcal{a}(\sin b)dbdc/8\pi^2$ in the Euler parameterization $\mathcal{R} = \mathcal{R}_Z(a)\mathcal{R}_Y(b)\mathcal{R}_Z(c)$, where angles $a$ and $c$ span an interval of length $2\pi$, angle $b$ spans $[0, \pi]$ and the direct Cartesian reference frame $\mathcal{O}\mathbf{X}\mathbf{Z}$ is of any orientation with respect to the reference frame $\mathcal{O}\mathbf{x}\mathbf{z}$ (see [53, Section 8.2]); it follows that

$$ 
\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} \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}(\langle l, m_z | R | l', m'_z \rangle)^* \frac{\delta(u\mathbf{e}_x - u'\mathcal{R}\mathbf{e}_x)}{v^2 + v'^2 + \frac{2\alpha u^2}{(1+\alpha)^2}}
\end{aligned} 
\quad (91) $$

where operator $\mathcal{R}$ represents rotation $\mathcal{R}$ in the Hilbert space of a quantum particle. The choice of axes $\mathcal{O}\mathbf{Z} = \mathcal{O}\mathbf{z}$ and $\mathcal{O}\mathbf{Y} = \mathcal{O}\mathbf{x}$ leads to $(\sin b)\delta(u\mathbf{e}_x - u'\mathcal{R}\mathbf{e}_x) = (\sin b)\delta(u\mathcal{R}_z(-a)\mathbf{e}_x - u'\mathcal{R}_x(b)\mathcal{R}_z(c)\mathbf{e}_x) = \delta(u\cos a - u'\cos c)\delta(-u\sin a - u'\sin c\cos b)\delta(-u'\sin c) = [\delta(u - u')/(u u')][\delta(a)\delta(c) + \delta(a - \pi)\delta(c - \pi)]$ where we successively used the rotational invariance of the three dimensional Dirac distribution and decomposed its action into Dirac distributions along $\mathcal{O}\mathbf{x}$, $\mathcal{O}\mathbf{y}$ and $\mathcal{O}\mathbf{z}$. Integration in SO(3) simply reduces to the line $\mathcal{R} = \mathcal{R}_x(b)$ if $a=c=0$ and to the line $\mathcal{R} = \mathcal{R}_x(-b)$ if $a=c=\pi$, which is taken into account by an integration on $b$ extended to $[-\pi, \pi]$, giving rise to a projector on the state of zero angular momentum along $\mathcal{O}\mathbf{x}$,

$$ \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. $$

We get

$$ 
\begin{aligned}
\mathcal{E} = & \sum_{m_z}^{\epsilon} \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}^{\epsilon} \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 \frac{\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{2uu'}{1+\alpha} (\cos T_1 \cos T'_1 + \sin T_1 \sin T'_1 \cos \phi)} \right]}{\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}} \\
& + \frac{\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_1|^2 + |Z'_1|^2 + \frac{2\alpha uu'}{1+\alpha} (\cos T_2 \cos T'_2 + \sin T_2 \sin T'_2 \cos \phi)} \right]}{\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}} \\
& + \frac{(-\epsilon)\delta(u-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) $$

where we have introduced the complex numbers $Z_1 = \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)$ and their counterparts for the primed variables $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)$. Finally, let us take advantage of scale invariance by means of the ansatz

$$f_{m_z}^{(l)}(u, v, w) = \frac{(\mathrm{ch} x)^{s+3/2} (1 + e^{-2x})^{-s/2} \Phi_{m_z}^{(l)}(x, w)}{(u^2 + v^2/\zeta^2)^{(s+7/2)/2}} \quad \text{where } v = \zeta e^x u \quad \text{and} \quad \zeta = \frac{\sqrt{2\alpha}}{1+\alpha} \quad (93)$$

chosen cleverly so that the diagonal part of the functional (first contribution in (92)) becomes scalar (independent of any variable) and so that the singular part (last contribution) does not depend on the mass ratio or even on the scaling exponent $s$. The integration on $u$ brings out as in reference [40] an infinite constant factor, here $(\zeta^4/16) \int_0^{+\infty} du/u$, to give the regularized functional

$$\tilde{\mathcal{E}} = \sum_{m_z}^{\epsilon} \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 \to -\infty}^{\epsilon} dx dx' \int_{-1}^{+\infty} 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} \\ \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) (\mathrm{Re} \, z_1 \mathrm{Re} \, z'_1 + \mathrm{Im} \, z_1 \mathrm{Im} \, z'_1 \cos \phi)} \\ + \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) (\mathrm{Re} \, z_2 \mathrm{Re} \, z'_2 + \mathrm{Im} \, z_2 \mathrm{Im} \, z'_2 \cos \phi)} \\ + (-\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] \quad (94)$$

where we have set $\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)$ and the equivalent relations for the primed variables, $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)$ and where we recall that $\theta = \arccos w$ and $\theta' = \arccos w'$. A simple functional derivation of (94) with respect to $\Phi_{m_z}^{(\ell, \epsilon)}(x, w)$, taking a pure imaginary scale exponent and returning to the notation $u = \cos\theta$ of the main article gives the new form of the operator $M_{2,2}^{(\ell, \epsilon)}(iS)$, marked with a Czech accent to avoid confusion with ((34), (35), (36)):

$$\langle x, u | \langle \ell, m_z | \check{M}_{2,2}^{(\ell, \epsilon)}(iS) | f \rangle = 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') \quad (95)$$

with an integral kernel separated into three contributions $\check{K} = \check{K}_1 + \check{K}_2 + \check{K}_3$ written line by line:

$$\check{K}_{2,2}^{(\ell)}(x, u, m_z; x', u', m'_z) = \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)(\mathrm{Re}\,z_1\mathrm{Re}\,z'_1 + \mathrm{Im}\,z_1\mathrm{Im}\,z'_1 \cos\phi)} \\ + \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)(\mathrm{Re}\,z_2\mathrm{Re}\,z'_2 + \mathrm{Im}\,z_2\mathrm{Im}\,z'_2 \cos\phi)} \\ + \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}} \quad (96)$$
where we remembered that $\langle \ell, m_z | \ell, m_x = 0 \rangle = 0$ if $\epsilon \neq (-1)^{\ell}$. To verify that equation (96) is written in the same order as equation (36), i.e. that the kernel $\tilde{K}_i$ is precisely the kernel $K_i$ written in the formulation of reference [52], we show numerically that

$$ \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]. $$

Références

[1] Y. Castin, F. Werner, « The Unitary Gas and its Symmetry Properties », in BCS-BEC Crossover and the Unitary Fermi gas (W. Zwerger, éd.), Lecture Notes in Physics, vol. 836, Springer, 2011, p. 127-191.

[2] M. Randeria, E. Taylor, « Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the Unitary Fermi Gas », Ann. Rev. Cond. Matter Phys. 5 (2014), n° 1, article no. 209.

[3] M. W. Zwierlein, J. R. Abo-Shaeer, A. Schirotzek, C. H. Schunck, W. Ketterle, « Vortices and superfluidity in a strongly interacting Fermi gas », Nature 435 (2005), p. 1047-1051.

[4] L. A. Sidorenkov, M. K. Tey, R. Grimm, Y.-H. Hou, L. Pitaevskii, S. Stringari, « Second sound and the superfluid fraction in a Fermi gas with resonant interactions », Nature 498 (2013), p. 78-81.

[5] S. Nascimbène, N. Navon, K. J. Jiang, F. Chevy, C. Salomon, « Exploring the thermodynamics of a universal Fermi gas », Nature 463 (2010), p. 1057-1060.

[6] M. Horikoshi, S. Nakajima, M. Ueda, T. Mukaiyama, « Measurement of universal thermodynamic functions for a unitary Fermi gas », Science 327 (2010), n° 5964, p. 442-445.

[7] M. J. Ku, A. T. Sommer, L. W. Cheuk, M. W. Zwierlein, « Revealing the superfluid lambda transition in the universal thermodynamics of a unitary Fermi gas », Science 335 (2012), n° 6068, p. 563-567.

[8] C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. Bruun, F. Schreck, R. Grimm, « Metastability and coherence of repulsive polarons in a strongly interacting Fermi mixture », Nature 485 (2012), p. 615-618.

[9] C. Ravensbergen, E. Soave, V. Corre, M. Kreyer, B. Huang, E. Kirilov, R. Grimm, « Resonantly Interacting Fermi-Fermi Mixture of $^{161}$Dy and $^{40}$K », Phys. Rev. Lett. 124 (2020), n° 20, article no. 203402.

[10] E. Neri, A. Ciamei, C. Simonelli, I. Goti, M. Inguscio, A. Trenkwalder, M. Zaccanti, « Realization of a cold mixture of fermionic chromium and lithium atoms », Phys. Rev. A 101 (2020), n° 6, article no. 063602.

[11] K. Van Houcke, F. Werner, E. Kozik, N. Prokof'ev, B. Svistunov, M. J. H. Ku, A. T. Sommer, L. W. Cheuk, A. Schirotzek, M. W. Zwierlein, « Feynman diagrams versus Fermi-gas Feynman emulator », Nat. Phys. 8 (2012), p. 366-370.

[12] R. Rossi, « Contributions to unbiased diagrammatic methods for interacting fermions », Thèse, Université Paris sciences et lettres, Paris, France, 2017, thèse en ligne tel-01704724v2, https://tel.archives-ouvertes.fr/ tel-01704724v2.

[13] K. Huang, Statistical Mechanics, 2nd éd., John Wiley & Sons, New York, 1987.

[14] X.-J. Liu, « Virial expansion for a strongly correlated Fermi system and its application to ultracold atomic Fermi gases », Phys. Rep. 524 (2013), n° 2, p. 37-83.

[15] X.-J. Liu, H. Hu, « Virial expansion for a strongly correlated Fermi gas with imbalanced spin populations », Phys. Rev. A 82 (2010), n° 4, article no. 043626.

[16] R. Rossi, T. Ohgoe, K. Van Houcke, F. Werner, « Resummation of Diagrammatic Series with Zero Convergence Radius for Strongly Correlated Fermions », Phys. Rev. Lett. 121 (2018), n° 13, article no. 130405.

[17] E. Beth, G. E. Uhlenbeck, « The quantum theory of the non-ideal gas I. Deviations from the classical theory », Physica 3 (1936), n° 8, p. 729-745.

[18] E. Beth, G. E. Uhlenbeck, « The quantum theory of the non-ideal gas. II. Behaviour at low temperatures », Physica 4 (1937), n° 10, p. 915-924.

[19] L. Landau, E. Lifchitz, Physique statistique - 1ère partie, troisième éd., Éditions Mir, Moscou, 1984.

[20] Y. Castin, F. Werner, « Le troisième coefficient du viriel du gaz de Bose unitaire », Canadian Journal of Physics 91 (2013), n° 5, p. 382-389.

[21] C. Gao, S. Endo, Y. Castin, « The third virial coefficient of a two-component unitary Fermi gas across an Efimov-effect threshold », Eur. Phys. Lett. 109 (2015), n° 1, article no. 16003.

[22] P. C. Hemmer, « The hard core quantum gas at high temperatures », Phys. Lett., A 27 (1968), n° 6, p. 377-378.

[23] B. Jancovici, « Quantum-Mechanical Equation of State of a Hard-Sphere Gas at High Temperature », Phys. Rev. 178 (1969), n° 1, p. 295-297.

[24] B. Jancovici, « Quantum-Mechanical Equation of State of a Hard-Sphere Gas at High Temperature. II », Phys. Rev. 184 (1969), n° 1, p. 119-123.
[25] B. Jancovici, S. P. Merkuriev, «Quantum-mechanical third virial coefficient of a hard-sphere gas at high temperature», *Phys. Rev. A* **12** (1975), n° 6, p. 2610-2621.

[26] T. D. Lee, C. N. Yang, «Many-Body Problem in Quantum Statistical Mechanics. II. Virial Expansion for Hard-Sphere Gas», *Phys. Rev. 116* (1959), n° 1, p. 25-31.

[27] A. Pais, G. E. Uhlenbeck, «On the Quantum Theory of the Third Virial Coefficient», *Phys. Rev. 116* (1959), n° 2, p. 250-269.

[28] S. K. Adhikari, R. D. Amado, «Low-Temperature Behavior of the Quantum Cluster Coefficients», *Phys. Rev. Lett.* **27** (1971), n° 8, p. 485-487.

[29] W. G. Gibson, «Low-Temperature Expansion of the Third-Cluster Coefficient of a Quantum Gas», *Phys. Rev. A* **6** (1972), n° 6, p. 2469-2477.

[30] A. Comtet, Y. Georgelin, S. Ouvry, «Statistical aspects of the anyon model», *J. Phys. A, Math. Gen.* **22** (1989), n° 18, p. 3917-3925.

[31] J. McCabe, S. Ouvry, «Perturbative three-body spectrum and the third virial coefficient in the anyon model», *Phys. Lett. B* **260** (1991), p. 113-119.

[32] X.-J. Liu, H. Hu, P. D. Drummond, «Virial Expansion for a Strongly Correlated Fermi Gas», *Phys. Rev. Lett.* **102** (2009), n° 16, article no. 160401.

[33] X.-J. Liu, H. Hu, P. D. Drummond, «Three attractively interacting fermions in a harmonic trap: Exact solution, ferromagnetism, and high-temperature thermodynamics», *Phys. Rev. A* **82** (2010), n° 2, article no. 023619.

[34] F. Werner, Y. Castin, «Unitary gas in an isotropic harmonic trap: Symmetry properties and applications», *Phys. Rev. A* **74** (2006), n° 5, article no. 053604.

[35] V. Efimov, «Energy levels of three resonantly interacting particles», *Nucl. Phys. A* **210** (1973), n° 1, p. 157-188.

[36] F. Werner, Y. Castin, «Unitary Quantum Three-Body Problem in a Harmonic Trap», *Phys. Rev. Lett.* **97** (2006), n° 15, article no. 150401.

[37] S. T. Rittenhouse, N. P. Mehta, C. H. Greene, «Green's functions and the adiabatic hyperspherical method», *Phys. Rev. A* **82** (2010), n° 2, article no. 022706.

[38] Y. Castin, E. Tignone, «Trimers in the resonant (2+1)-fermion problem on a narrow Feshbach resonance: Crossover from Efimovian to hydrogenoid spectrum», *Phys. Rev. A* **84** (2011), n° 6, article no. 062704.

[39] Y. Castin, C. Mora, L. Pricoupenko, «Four-Body Efimov Effect for Three Fermions and a Lighter Particle», *Phys. Rev. Lett.* **105** (2010), n° 22, article no. 223201.

[40] S. Endo, Y. Castin, «Absence of a four-body Efimov effect in the 2+2 fermionic problem», *Phys. Rev. A* **92** (2015), n° 5, article no. 053624.

[41] S. Endo, Y. Castin, «The interaction-sensitive states of a trapped two-component ideal Fermi gas and application to the virial expansion of the unitary Fermi gas», *J. Phys. A, Math. Gen.* **49** (2016), n° 26, article no. 265301.

[42] Y. Yan, D. Blume, «Path integral Monte Carlo Determination of the Fourth-Order Virial Coefficient for Unitary Two-Component Fermi Gas with Zero-Range Interactions», *Phys. Rev. Lett.* **116** (2016), n° 23, article no. 230401.

[43] Y. Hou, K. J. Morrell, A. J. Czejdo, J. E. Drut, «Fourth- and fifth-order virial expansion of harmonically trapped fermions at unitarity», *Phys. Rev. Research* **3** (2021), n° 3, article no. 033099.

[44] A. J. Czejdo, J. E. Drut, Y. Hou, K. J. Morrell, «Toward an automated-algebra framework for high orders in the virial expansion of quantum matter», *Condens. Matter* **7** (2022), n° 1, article no. 13.

[45] Y. Hou, J. E. Drut, «Fourth- and Fifth-Order Virial Coefficients from Weak Coupling to Unitarity», *Phys. Rev. Lett.* **125** (2020), n° 5, article no. 050403.

[46] D. S. Petrov, «Three-body problem in Fermi gases with short-range interparticle interaction», *Phys. Rev. A* **67** (2003), n° 1, article no. 010703.

[47] S. Endo, Y. Castin, «The interaction-sensitive states of a trapped two-component ideal Fermi gas and application to the virial expansion of the unitary Fermi gas», post-publication hal-01246611v5, https://hal.archives-ouvertes.fr/ hal-01246611, 2021.

[48] A. L. Gaunt, T. F. Schmidutz, I. Gotlibovych, R. P. Smith, Z. Hadzibabic, «Bose-Einstein condensation of atoms in a uniform potential», *Phys. Rev. Lett.* **110** (2013), n° 20, article no. 200406.

[49] B. Mukherjee, Z. Yan, P. B. Patel, Z. Hadzibabic, T. Yefsah, J. Struck, M. W. Zwierlein, «Homogeneous Atomic Fermi Gases», *Phys. Rev. Lett.* **118** (2017), n° 2, article no. 123401.

[50] T.-L. Ho, Q. Zhou, «Obtaining the phase diagram and thermodynamic quantities of bulk systems from the densities of trapped gases», *Nat. Phys.* **6** (2010), n° 2, p. 131-134.

[51] R. Dum, M. Olshanii, «Gauge Structures in Atom-Laser Interaction: Bloch Oscillations in a Dark Lattice», *Phys. Rev. Lett.* **76** (1996), n° 11, p. 1788-1797.

[52] L. Pricoupenko, «Isotropic contact forces in arbitrary representation: Heterogeneous few-body problems and low dimensions», *Phys. Rev. A* **83** (2011), n° 6, article no. 062711.

[53] W.-K. Tung, *Group Theory in Physics*, World Scientific, Philadelphia, 1985.

[54] C. Mora, Y. Castin, L. Pricoupenko, «Integral equations for the four-body problem», *C. R. Physique* **12** (2011), n° 1, p. 71-85.