# Fonctions de distribution de paires d'un gaz superfluide de fermions de spin 1/2 en interaction de contact dans la théorie BCS dépendant du temps linéarisée

*Pair distribution functions of a superfluid spin 1/2*

Fermi gas with contact interactions in the linearized time-dependent BCS theory

Yvan Castin<sup>a</sup>

<sup>a</sup> Laboratoire Kastler Brossel, ENS-Université PSL, CNRS, Université Sorbonne et Collège de France, 24 rue Lhomond, 75231 Paris, France

Courriel : yvan.castin@lkb.ens.fr (Y. Castin)

**Résumé.** Nous montrons que la théorie de champ moyen minimale à utiliser pour le calcul des fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ d'un gaz superfluide de fermions de spin 1/2 spatialement homogène non polarisé n'est pas la théorie BCS statique ordinaire, mais la théorie BCS dépendant du temps linéarisée, mise en œuvre par le truchement du théorème de fluctuation-dissipation. En effet, la première ignore totalement la branche d'excitation acoustique – les phonons – du superfluide, alors que la seconde en tient compte explicitement, ainsi que des fluctuations quantiques induites par le continuum de paires brisées. Contrairement à la première, la seconde théorie (i) répercute l'effet de ces excitations collectives sur l'équation d'état du système, y compris à température nulle, (ii) permet à la fonction $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ de descendre à distance assez grande strictement en dessous de sa valeur asymptotique ($\rho/2)^2$ où $\rho$ est la densité du gaz, comme il se doit d'après l'hydrodynamique quantique de Landau et Khalatnikov à basse température, et (iii) prédit dans la fonction $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ à courte distance des contributions sous-dominantes en $|\mathbf{r}-\mathbf{r}'|^2 \ln|\mathbf{r}-\mathbf{r}'|$ à 3D et en $|\mathbf{r}-\mathbf{r}'|^2 \ln(-\ln|\mathbf{r}-\mathbf{r}'|)$ à 2D, à côté des contributions dominantes en $|\mathbf{r}-\mathbf{r}'|$ à 3D et en $|\mathbf{r}-\mathbf{r}'|^2 \ln|\mathbf{r}-\mathbf{r}'|$ à 2D déjà présentes dans la théorie BCS statique mais avec un coefficient plus faible. Cette discussion est pertinente pour les travaux théoriques récents d'Obeso-Jureidini et de Romero-Rochín, et pour les expériences en cours sur les gaz d'atomes froids à l'ENS et au MIT.

**Abstract.** We show that the minimal mean-field theory to use for calculating the pair distribution functions $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ of a spatially homogeneous, unpolarized spin-1/2 superfluid Fermi gas is not the ordinary static BCS theory, but the linearized time-dependent BCS theory implemented via the fluctuation-dissipation theorem. Indeed, the former completely ignores the acoustic excitation branch – the phonons – of the superfluid, while the latter explicitly takes it into account, as well as the quantum fluctuations induced by the broken-pair continuum. Unlike the first, the second theory (i) reflects the effect of these collective excitations on the system's equation of state, including at zero temperature, (ii) allows the function $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ to go at sufficiently large distances strictly below its asymptotic value $(\rho/2)^2$ where $\rho$ is the gas density, as
expected according to the quantum hydrodynamics of Landau and Khalatnikov at low temperatures, and (iii) predicts in the function $g_{||}(\mathbf{r}, \mathbf{r}')$ at short distances subdominant contributions $|\mathbf{r} - \mathbf{r}'|^2 \ln|\mathbf{r} - \mathbf{r}'|$ in 3D and $|\mathbf{r} - \mathbf{r}'|^2 \ln(-\ln|\mathbf{r} - \mathbf{r}'|)$ in 2D, alongside the dominant contributions $|\mathbf{r} - \mathbf{r}'|$ in 3D and $|\mathbf{r} - \mathbf{r}'|^2 \ln|\mathbf{r} - \mathbf{r}'|$ in 2D already present in static BCS theory but with a lower coefficient. This discussion is relevant to the recent theoretical work of Obeso-Jureidini and Romero-Rochín, and to the ongoing experiments on cold atomic gases at ENS and MIT.

**Mots-clés.** Gaz de fermions, densité de paires, facteur de structure, théorème de fluctuation-dissipation, superfluidité, limite unitaire, interactions de contact.

**Keywords.** Fermi gases, pair density, structure factor, fluctuation-dissipation theorem, superfluidity, unitary limit, contact interactions.

Table des matières

<table><tr><td>1. Introduction, motivations, insuffisance de la théorie BCS statique</td><td>3</td></tr><tr><td>2. Obtention des g<sub>σσ'</sub>(r, r')</td><td>6</td></tr><tr><td>2.1. Interaction, état du système et brefs rappels sur la théorie BCS statique</td><td>7</td></tr><tr><td>2.2. L'obtention des δg<sub>σσ'</sub><sup>dtl</sup>(r, r')</td><td>11</td></tr><tr><td>2.3. Isoler la partie collective dans les δg<sub>σσ'</sub><sup>dtl</sup>(r, r')</td><td>13</td></tr><tr><td>2.4. Différentes formulations en termes d'intégrales curvilignes et applications</td><td>14</td></tr><tr><td>2.5. Quelle équation d'état dans notre théorie BCS dépendant du temps linéarisée?</td><td>18</td></tr><tr><td>3. Résultats explicites sur les g<sub>σσ'</sub><sup>dtl</sup>(r, r')</td><td>20</td></tr><tr><td>3.1. Comportement des δg<sub>σσ'</sub><sup>dtl</sup>(r, r') à grande distance et basse température</td><td>20</td></tr><tr><td>3.1.1. Ce qu'on déduit de l'intégration par parties</td><td>20</td></tr><tr><td>3.1.2. Ce qu'on déduit d'une formulation par intégrale de contour</td><td>22</td></tr><tr><td>3.1.3. Ce qu'en dit l'hydrodynamique quantique</td><td>24</td></tr><tr><td>3.2. Résultats numériques</td><td>25</td></tr><tr><td>3.2.1. Cas tridimensionnel; comparaison au Monte-Carlo quantique de Giorgini</td><td>25</td></tr><tr><td>3.2.2. Cas bidimensionnel; comparaison aux expériences de l'ENS</td><td>29</td></tr><tr><td>3.3. Comportement des g<sub>σσ'</sub><sup>dtl</sup>(r, r') à courte distance</td><td>30</td></tr><tr><td>3.3.1. Étude de δg̃<sub>coll</sub>(q) pour q → +∞</td><td>30</td></tr><tr><td>3.3.2. Étude de g<sub>||</sub><sup>dtl</sup>(r, r') pour |r - r'| → 0+</td><td>32</td></tr><tr><td>3.3.3. Étude de g<sub>||</sub><sup>dtl</sup>(r, r') pour |r - r'| → 0+</td><td>34</td></tr><tr><td>4. Conclusion générale</td><td>36</td></tr><tr><td>Annexe A. Développement à faible q de δg̃<sub>coll</sub>(q) à l'ordre zéro en exp(-βE<sub>gap</sub>)</td><td>38</td></tr><tr><td>A.1. Ce que donne le continuum de paire brisée</td><td>38</td></tr><tr><td>A.2. Ce que donne le phonon de Bogoliubov-Anderson</td><td>39</td></tr><tr><td>A.3. Applications</td><td>40</td></tr><tr><td>Annexe B. Développement à grand q de δg̃<sub>coll</sub>(q) à l'ordre zéro en exp(-βE<sub>gap</sub>)</td><td>41</td></tr><tr><td>B.1. Deux astuces simplificatrices et une table récapitulative</td><td>41</td></tr><tr><td>B.2. Cas bidimensionnel</td><td>44</td></tr><tr><td>B.3. Cas tridimensionnel</td><td>45</td></tr><tr><td>Annexe C. Comportement à courte distance de certaines intégrales de Fourier</td><td>50</td></tr><tr><td>C.1. Cas bidimensionnel</td><td>50</td></tr><tr><td>C.2. Cas tridimensionnel</td><td>51</td></tr><tr><td>Références</td><td>52</td></tr></table>
# 1. Introduction, motivations, insuffisance de la théorie BCS statique

Les fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ d'un gaz de fermions à deux états internes $\uparrow, \downarrow$ en interaction de contact, à l'équilibre thermique dans le cas spatialement homogène et non polarisé (même densité $\rho_\sigma = \rho/2$ dans chaque composante), définies comme

$$g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}') \equiv \rho_{\sigma}\rho_{\sigma'} + \delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}') \equiv \langle \hat{\psi}_{\sigma}^{\dagger}(\mathbf{r}) \hat{\psi}_{\sigma'}^{\dagger}(\mathbf{r}') \hat{\psi}_{\sigma'}(\mathbf{r}') \hat{\psi}_{\sigma}(\mathbf{r}) \rangle \quad (1)$$

où $\hat{\psi}_\sigma(\mathbf{r})$ est l'opérateur de champ fermionique pour le spin $\sigma \in \{\uparrow, \downarrow\}$ et $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ sont les écarts aux valeurs asymptotiques, ont fait l'objet récemment d'études approfondies [1, 2] en dimension $d=2$ ou $d=3$ à température $T=0$ au moyen de la théorie BCS ordinaire [3], dite *statique*, une théorie variationnelle réduisant l'état fondamental du gaz à un état cohérent de paires liées $\uparrow\downarrow$ (les équivalents pour notre gaz neutre des paires de Cooper des supraconducteurs, ici dans l'état singulet de spin). Nous nous proposons cependant, dans cet article, de souligner les insuffisances de cette théorie et d'en utiliser une meilleure, dite *dépendant du temps linéarisée*, pour le calcul des fonctions $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$.

En effet, la théorie BCS statique a le génie de tenir compte du phénomène de condensation par paires du gaz, de donner une borne variationnelle sur l'énergie et de prédire un spectre d'excitation $\epsilon_k$ de quasi-particules fermioniques avec une bande interdite de largeur $E_{gap} > 0$, le fameux spectre BCS. Mais elle n'est pas adaptée à toutes les observables, en particulier à des observables aussi multimodes que $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$. On le voit bien dans la limite CBE où l'on fait tendre vers zéro la densité $\rho$ dans le cas où deux fermions $\uparrow$ et $\downarrow$ admettent un état lié dans le vide : la théorie BCS statique se réduit alors à un condensat pur de dimères bosoniques; or, dans un gaz de bosons condensé, un ansatz de condensat pur conduit à une fonction de distribution de paires $g_2^B(\mathbf{r}, \mathbf{r}')$ triviale, partout égale – après prise de la limite thermodynamique – à la densité $\rho_B$ des bosons au carré, donc à un écart $\delta g_2^B(\mathbf{r}, \mathbf{r}')$ à la valeur asymptotique $\rho_B^2$ identiquement nul. La théorie multimode de Bogolioubov [4], qui tient compte correctement de la fraction non condensée et des modes de quasi-particules bosoniques – dites justement quasi-particules de Bogolioubov – montre toutefois qu'il n'en est rien; en particulier, à température nulle, elle prédit que $\delta g_2^B(\mathbf{r}, \mathbf{r}')$ admet une queue en $-1/|\mathbf{r} - \mathbf{r}'|^{d+1}$, et que l'intégrale de $\delta g_2^B(\mathbf{r}, \mathbf{r}' = 0)$ sur tout l'espace $\mathbf{r} \in \mathbb{R}^d$ prend la valeur strictement négative $-\rho_B$, voir la section 8.7.2 de la référence [5]. Rappelons l'adage : zéro est une mauvaise approximation pour toute quantité non nulle...

Mais étendons cette objection aux fermions dans le cas général. L'insuffisance de la théorie BCS statique, à température nulle ou pas, se voit bien sur la distribution de paires de fermions de spins opposés : comme l'ansatz variationnel est un état gaussien pour le champ fermionique, le théorème de Wick s'applique et l'on a

$$\delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = |\langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}}|^2 \quad (2)$$

où « sta » mis en indice ou en exposant signifie « BCS statique ». Or, en appliquant le critère de Penrose et Onsager [6] à la matrice densité à deux corps $\langle \hat{\psi}_{\sigma_1}^{\dagger}(\mathbf{r}_1) \hat{\psi}_{\sigma_2}^{\dagger}(\mathbf{r}_2) \hat{\psi}_{\sigma_2'}(\mathbf{r}_2') \hat{\psi}_{\sigma_1'}(\mathbf{r}_1') \rangle_{\text{sta}}$ dans le volume de quantification $[0, L]^d$ comme dans les références [7, 8], c'est-à-dire en identifiant son vecteur propre macroscopiquement peuplé – correctement normalisé dans $[0, L]^{2d}$ – à la fonction
d'onde $\phi_0(\mathbf{r}_1, \mathbf{r}_2)$ de paire liée $\uparrow\downarrow$ du condensat BCS, il vient<fnref n="1" />

$$ \phi_0(\mathbf{r}_1, \mathbf{r}_2) = \langle \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle_{\text{sta}} / (\bar{N}_0^{\text{paires}})^{1/2} \quad \text{où} \quad \bar{N}_0^{\text{paires}} = \int_{[0,L]^d} d^d \mathbf{r}_1 \int_{[0,L]^d} d^d \mathbf{r}_2 |\langle \hat{\psi}_\downarrow(\mathbf{r}_2) \hat{\psi}_\uparrow(\mathbf{r}_1) \rangle_{\text{sta}}|^2 \quad (3) $$

est le nombre moyen de paires condensées, si bien que

$$ \delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \bar{N}_0^{\text{paires}} |\phi_0(\mathbf{r}, \mathbf{r}')|^2 \quad (4) $$

La prédiction (2) décrit donc seulement la structure interne d'une paire liée $\uparrow\downarrow$ mais correspond à une distribution de paires liées plate, donc à une décorrélation complète entre les positions desdites paires, voir la figure 1 purement illustrative. Il en résulte qu'aux distances relatives $|\mathbf{r}-\mathbf{r}'| \gg \ell_{\text{paire}}$, où $\ell_{\text{paire}}$ est la taille d'une paire liée, la fonction $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ se réduit tristement au produit des densités $\rho_\uparrow^{\text{sta}} \rho_\downarrow^{\text{sta}}$ à un exponentiellement petit près, la moyenne anormale $\langle \hat{\psi}_\downarrow(\mathbf{r}') \hat{\psi}_\uparrow(\mathbf{r}) \rangle_{\text{sta}}$ étant à décroissance rapide pour la variable $|\mathbf{r}-\mathbf{r}'|$ [1, 2].<fnref n="2" /> En particulier, il saute

$$ \langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} \stackrel{d=3}{\sim} \frac{m\Delta}{2\pi^2\hbar^2|\mathbf{r}-\mathbf{r}'|} \mathrm{Im} \int_0^{+\infty} d\epsilon \frac{\exp[-i\mathbf{k}\epsilon](\mathbf{r}-\mathbf{r}')}{(\Delta^2 + \epsilon^2)^{1/2}} \mathrm{Re} \left[ i\kappa_0^{1/2} \exp(-i\kappa_0|\mathbf{r}-\mathbf{r}'|) \right] \quad (5) $$

avec $k_0 = [2m(\mu - i\Delta)/\hbar^2]^{1/2}$, $\mu$ le potentiel chimique du gaz, $m$ la masse d'un fermion, $\Delta > 0$ le paramètre d'ordre écrit à l'ordre zéro en exp($-\beta E_{\text{gap}}$) (voir la section 2.1). Pour cela, nous effectuons l'intégrale angulaire dans l'expression (25) et utilisons la parité formelle en $k$ de l'intégrande pour écrire le résultat sous la forme $\langle \hat{\psi}_\downarrow(\mathbf{r}') \hat{\psi}_\uparrow(\mathbf{r}) \rangle_{\text{sta}} = [\Delta/(8i\pi^2|\mathbf{r}-\mathbf{r}'|)] \int_{-\infty}^{+\infty} dk \, k \exp(-ik|\mathbf{r}-\mathbf{r}'|)/\epsilon_k$. La fonction $k \to \epsilon_k$ prolongée analytiquement aux valeurs de $k$ complexes admet dans chaque quadrant une ligne de coupure partant de $k_0^*$, $-k_0$, $-k_0^*$ ou $k_0$. Déformons alors comme dans [2] le chemin d'intégration en le rabattant dans le demi-plan inférieur sur les lignes de coupure issues de $k_0$ et $-k_0^*$, la contribution de la seconde étant le complexe conjugué de la première. Il reste à paramétrer la première comme suit, $\epsilon \in [0, +\infty] \mapsto k(\epsilon) = [2m(\mu - i\sqrt{\Delta^2 + \epsilon^2})/\hbar^2]^{1/2}$, et à remarquer que $\epsilon_{k(\epsilon)}$ saute de $-i\epsilon$ à $+i\epsilon$ si l'on franchit la ligne de coupure de haut en bas, pour aboutir à

$$ \langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} \stackrel{d=3}{=} \frac{m\Delta}{2\pi^2\hbar^2|\mathbf{r}-\mathbf{r}'|} \mathrm{Im} \int_0^{+\infty} d\epsilon \frac{\exp[-i\mathbf{k}\epsilon](\mathbf{r}-\mathbf{r}')}{(\Delta^2 + \epsilon^2)^{1/2}} \quad (6) $$

Si $|\mathbf{r}-\mathbf{r}'| \to +\infty$, cette intégrale est dominée par l'intervalle $\epsilon = O(1/|\mathbf{r}-\mathbf{r}'|^{1/2})$ donc par un voisinage du point de branchement $k_0$; il reste à développer $k(\epsilon)$ à l'ordre $\epsilon^2$ sous l'exponentielle et à négliger $\epsilon$ sous la racine carrée au dénominateur de l'intégrande, pour obtenir – après intégration gaussienne – l'équivalent annoncé. La même procédure donne pour la moyenne normale (26) :

$$ \langle \hat{\psi}_{\downarrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\uparrow}(\mathbf{r}') \rangle_{\text{sta}} \stackrel{d=3}{=} \frac{-m}{2\pi^2 \hbar^2 |\mathbf{r}-\mathbf{r}'|} \mathrm{Re} \int_0^{+\infty} d\epsilon \frac{\exp[-i\mathbf{k}\epsilon](\mathbf{r}-\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'| + +\infty} \sim \frac{(m\Delta/\hbar^2)^{1/2}}{(2\pi|\mathbf{r}-\mathbf{r}'|)^{3/2}} \mathrm{Im} \left[ i\kappa_0^{1/2} \exp(-i\kappa_0|\mathbf{r}-\mathbf{r}'|) \right] \quad (7) $$

En dimension $d=2$, toujours à $T=0$, nous utilisons les expressions analytiques (30,31) des moyennes normale et anormale établies dans la référence [1], le comportement asymptotique connu des fonctions de Bessel et la propriété spécifique à 2D que $k_0 = [2m(\mu - i\Delta)/\hbar^2]^{1/2}$ s'écrit aussi $k_0 = k_F^{\text{sta}} - iq_{\text{dim}}$ avec les notations et les relations (29) de la section 2.1. Alors, après regroupement astucieux des deux moyennes en parties réelle et imaginaire d'une même quantité complexe (comme on pourrait le faire aussi à 3D), on trouve que

$$ \langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} + i \langle \hat{\psi}_{\downarrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\uparrow}(\mathbf{r}') \rangle_{\text{sta}} \stackrel{d=2}{\sim} \frac{(m\Delta/\hbar^2)^{1/2}}{2\pi|\mathbf{r}-\mathbf{r}'|} \exp[i(\pi/4 - k_0|\mathbf{r}-\mathbf{r}'|)] \quad (8) $$

ce qui est compatible avec les équations (40,41) moins précises de la référence [1]. À température non nulle $T \neq 0$, cas non considéré dans les références [1, 2], la fonction $k \to 1/\epsilon_k$ dans les intégrandes est remplacée par $k \to (1-2f_k)/\epsilon_k = \text{th}(\beta c_k/2)/\epsilon_k$ qui n'admet plus ni (i) ligne de coupure pour $k \in \mathbb{C}$ ni (ii) pôle pour $k$ racine de l'équation $\epsilon_k = 0$. En effet, (i) (th$z$)/$z$ est une fonction paire de $z$, qui garde la même valeur à la traversée d'une ligne de coupure de $k \to \epsilon_k$ – on peut dire aussi que son développement en série entière ne comporte que des puissances paires de $z$, et (ii) (th$z$)/$z$ tend vers 1 lorsque $z \to 0$, c'est une fonction analytique de $z$ sur un voisinage de l'origine. Ne restent donc comme singularités

<footnote n="1">¹L'état fondamental du gaz dans la théorie BCS statique est un état cohérent de paires $|\Psi_{BCS}\rangle \propto \exp\left[\int d^d\mathbf{r}_1 d^d\mathbf{r}_2 \alpha\phi_{BCS}(\mathbf{r}_1, \mathbf{r}_2) \hat{\psi}_\uparrow^\dagger(\mathbf{r}_1) \hat{\psi}_\uparrow^\dagger(\mathbf{r}_2)\right]|0\rangle$ où $\alpha$ est une amplitude complexe, $\phi_{BCS}(\mathbf{r}_1, \mathbf{r}_2)$ est une fonction d'onde de paire normalisée – nous reprenons la notation des références [1, 2] – et $|0\rangle$ est le vide de fermions. Il serait cependant incorrect d'assimiler $\phi_{BCS}$ à $\phi_0$. Ces fonctions ne sont pas proportionnelles; elles sont même bien distinctes, sauf dans la limite CBE où elles se réduisent à la fonction d'onde d'un dimère d'impulsion nulle. Ce piège à éviter illustre les effets subtils de l'antisymérisation fermionique.</footnote>

<footnote n="2">²En dimension $d=3$, nous trouvons à $T=0$ que</footnote>
<figure><img src="image_1.png" /></figure>

FIGURE 1 – Représentation purement schématique de la distribution spatiale des fermions ↑ et ↓ dans un gaz 2D condensé par paires, dans un régime d'interaction intermédiaire – entre limite CBE et limite BCS – où la taille $\ell_{\text{paire}}$ d'une paire liée est plus faible mais du même ordre de grandeur que la distance moyenne $\rho^{-1/d}$ entre particules. Les fermions ↑ (étoiles noires) et ↓ (étoiles rouges) dans un même cercle appartiennent à la même paire liée (les cercles sont centrés sur les barycentres des paires liées). Comme le montre l'équation (4), la fonction de distribution $g_{↑↓}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ de la théorie BCS statique tient compte des corrélations spatiales entre deux fermions de spins opposés au sein d'une même paire liée, mais ignore celles entre deux paires liées différentes, car elle néglige les modes acoustiques du superfluide et les fluctuations de densité qu'ils induisent.

aux yeux que $g_{↑↓}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ est partout supérieure à sa valeur asymptotique.

En réalité, dans les gaz de fermions en interaction à courte portée, il existe à basse température (ici à $T=0$ pour simplifier) des excitations collectives bosoniques totalement ignorées par la théorie BCS statique. Elles correspondent à une branche d'excitation acoustique – par ondes sonores – dont le départ est linéaire en nombre d'onde,

$$ \omega_{\mathbf{q}} = c q + O(q^3) \quad (10) $$

où $c$ est la vitesse du son dans l'état fondamental du gaz. Les quanta associés sont donc des phons de vecteur d'onde $\mathbf{q}$ et d'énergie $\hbar\omega_{\mathbf{q}}$ ($\omega_{\mathbf{q}}$ est la pulsation propre de l'onde sonore correspondante). Ils jouent le rôle, dans le cas des superfluides de fermions, des quasi-particules de Bogolioubov des gaz de bosons. Bogolioubov – encore lui – et Anderson, en contemporains avisés de Bardeen, Cooper et Schrieffer, avaient bien vu l'énorme lacune de la référence [3] et y avaient

$$ \langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} \left| 1 + i\pi \frac{k_B T}{\Delta} \right| + i \langle \hat{\psi}_{\downarrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\uparrow}(\mathbf{r}') \rangle_{\text{sta}} \frac{T^{\mp 0}}{|\mathbf{r}-\mathbf{r}'|} - \frac{mk_B T/\hbar^2}{(\pi|\mathbf{r}-\mathbf{r}'|)^{(d-1)/2}} \times \begin{cases} i e^{-ik_0(T)|\mathbf{r}-\mathbf{r}'|} & \text{en } d=3 \\ \left[ \frac{2}{k_0(T)} \right]^{1/2} e^{i[\pi/4 - k_0(T)|\mathbf{r}-\mathbf{r}'|]} & \text{en } d=2 \end{cases} \quad (9) $$

<footnote>dans le plan complexe que les pôles thermiques fermioniques, tels que $\exp(\beta\epsilon_k) = -1$; ceux du quadrant inférieur droit s'écrivent $k_n(T) = [2m(\mu - i)\Delta + i\pi(2n + 1)k_B T]]^{1/2}$, $n \in \mathbb{N}$, les autres s'en déduisent par changement de signe global et conjugaison complexe. La formulation astucieuse par intégrale de contour s'applique toujours (le facteur thermique $f_k$ ne change pas la parité formelle en $k$ de l'intégrande et n'empêche pas de refermer le chemin d'intégration par un demi-cercle à l'infini dans le demi-plan inférieur). À grande distance, d'après le théorème des résidus, les moyennes anormale et normale sont dominées par les pôles les plus proches de l'axe réel, à savoir $k_0(T)$ et $-k_0^*(T)$; en suivant la section 8.7.3.2 de la référence [5], en particulier pour le cas plus délicat $d=2$, on obtient les comportements asymptotiques</footnote>
rapidement remédié [9, 10] au moyen d'un formalisme assez lourd, celui de l'Approximation de la Phase Aléatoire (RPA) pour Anderson. D'où le nom de « Bogoliubov-Anderson » souvent accolé à cette branche acoustique et à ces phonons [11].

Il a été compris depuis que l'on peut obtenir la même relation de dispersion $q \mapsto \omega_q$ plus simplement, en généralisant trivialement la théorie BCS au cas dépendant du temps (à $T=0$, l'état cohérent est dans un mode spatial de paire $\updownarrow$ dépendant du temps) puis en la linéarisant autour de son état stationnaire, c'est-à-dire autour de la théorie BCS statique [12]. Une vérification et une comparaison soigneuse à la RPA ont notamment été faites dans la référence [13]. L'analogie avec les gaz de bosons devient évidente : la théorie BCS dépendant du temps est l'analogue de l'équation de Gross-Pitayevski dépendant du temps ; la linéarisation de la première reproduit le spectre d'excitation phononique de Bogoliubov-Anderson, tout comme la linéarisation de la seconde reproduit celui des quasi-particules de Bogoliubov. Dans les deux cas, une théorie de champ classique reproduit les pulsations propres d'une théorie de champ quantique.

Le point central est maintenant, le lecteur l'aura compris, que les amplitudes des modes de phonons – des ondes sonores dans le gaz – sont le siège de fluctuations thermiques et quantiques, ces dernières subsistant à $T=0$ même s'il n'y a pas de phonons présents (une affirmation rendue tangible et quantitative par l'hydrodynamique quantique de Landau et Khalatnikov dans la section 3.1). Ceci entraîne inévitablement des fluctuations de la densité des paires liées donc l'apparition d'une structure non triviale dans $g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$ aux distances $|\mathbf{r}-\mathbf{r}'| \gg \ell_{\text{paire}}$. C'est cet effet qui conduit à l'existence des parties négatives de $\delta g_2^B(\mathbf{r}, \mathbf{r}')$, mentionnées plus haut, dans un gaz de bosons assez froid. Dans le cas des fermions, à suffisamment basse température, on doit donc s'attendre à perdre la propriété $\delta g_{\updownarrow}(\mathbf{r}, \mathbf{r}') \ge 0$ imprudemment déduite du théorème de Wick, puisque les paires liées $\updownarrow$ ressemblent à des bosons pour des phonons de grandes longueurs d'onde $\lambda_\phi \gg \ell_{\text{paire}}$.

La morale de cette histoire est que la théorie minimale acceptable pour le calcul des fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ dans un gaz de fermions de spin $1/2$ en interaction de contact doit tenir compte de leur branche d'excitation acoustique $\omega_q$. Ce ne peut être la théorie BCS statique, qui ne renferme que les excitations élémentaires de quasi-particules fermioniques $\epsilon_k$. Nous proposons ici comme candidat naturel son amélioration directe, la théorie BCS dépendant du temps linéarisée autour de sa solution stationnaire, qui garde la simplicité d'une théorie de champ moyen. Dans la suite, nous montrons comment obtenir les fonctions $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ dans cette théorie, en reliant fonctions de réponse et fonctions de corrélation au travers du théorème de fluctuation-dissipation dans la section 2, puis nous analysons en détail les prédictions obtenues dans la section 3, en comparant aussi à des résultats numériques au-delà du champ moyen et à des mesures expérimentales dans des gaz d'atomes froids. Nous concluons dans la section 4.

## 2. Obtention des $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ dans la théorie BCS dépendant du temps linéarisée

Dans cette section, nous mettons en place dans l'ensemble grand-canonique le formalisme BCS dépendant du temps linéarisé et ses fonctions de réponse densité-densité. Après la définition du problème et des rappels ciblés sur la théorie BCS statique dans la section 2.1, nous obtenons au moyen du théorème de fluctuation-dissipation dans la section 2.2 les écarts $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ des fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ à leurs valeurs asymptotiques sous forme d'intégrales des susceptibilités densité-densité $\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, \epsilon + i0^+)$ sur le vecteur d'onde $\mathbf{q}$ et l'énergie $\epsilon$. Nous isolons dans la section 2.3 le noyau dur de la théorie, c'est-à-dire une correction commune $\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+)$ de ces susceptibilités à leurs approximations BCS statique, donc une correction commune sur les écarts $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$. L'intégrale de $\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+)$ sur $\epsilon$ admet dans la section 2.4 d'utiles réécritures en termes d'intégrales de contour dans le plan complexe, qui en facilitent
l'évaluation numérique ou permettent de séparer les contributions du phonon de Bogolioubov-Anderson et du continuum de paire brisée. Enfin, nous montrons dans la section 2.5 que, lorsqu'on change les fonctions de corrélation de la densité par rapport à BCS statique comme nous le faisons ici, il faut aussi paradoxalement changer sa valeur moyenne, c'est-à-dire l'équation d'état grand-canonique du gaz, si l'on veut obtenir des valeurs sensées de $g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}')$ à courte distance; ceci aboutit à une belle réécriture de la correction commune des fonctions de distribution de paires à leurs approximations BCS statique.

## 2.1. Interaction, état du système et brefs rappels sur la théorie BCS statique

Nos fermions de spin 1/2 et de masse m vivent en dimension $d \in \{2;3\}$ avec des conditions aux limites périodiques dans la boîte de quantification $[0, L]^d$, dont on fait tendre le volume $L^d$ vers l'infini à la fin des calculs.<fnref n="3" /> Ils n'interagissent pas lorsqu'ils sont dans le même état interne $\sigma \in \{\uparrow, \downarrow\}$ (pas de résonance de diffusion dans l'onde $p$ entre $\uparrow$ et $\downarrow$, ni entre $\downarrow$ et $\uparrow$); dans le cas contraire, ils interagissent par un potentiel de portée négligeable et de longueur de diffusion $a_{dD}$ (résonance de diffusion dans l'onde $s$ entre $\uparrow$ et $\downarrow$). Nous suivons la convention que l'état stationnaire de diffusion à deux corps $\uparrow\downarrow$ à énergie nulle s'écrit (en dehors du potentiel)

$$ \psi_0(\mathbf{r}, \mathbf{r}') \stackrel{d=3}{=} 1 - a_{3D} / |\mathbf{r} - \mathbf{r}'| \quad \text{et} \quad \psi_0(\mathbf{r}, \mathbf{r}') \stackrel{d=2}{=} \ln(|\mathbf{r} - \mathbf{r}'| / a_{2D}) \quad (11) $$

À 2D, le problème à deux corps admet toujours un état lié, d'énergie de liaison

$$ \epsilon_{\text{dim}} = \hbar^2 q_{\text{dim}}^2 / m \quad \text{avec} \quad q_{\text{dim}} = 2 \exp(-\gamma) / a_{2D} \quad (12) $$

où $\gamma = 0,577215...$ est la constante d'Euler-Mascheroni. À 3D, il n'en admet un, d'énergie de liaison $\hbar^2/ma_{3D}^2$, que pour $a_{3D} > 0$.

Le gaz est pris à l'équilibre thermique dans l'ensemble grand-canonique à une température $T$ non nulle mais inférieure à la température de transition superfluide, $0 < T < T_c$. Comme il est non polarisé, il a même potentiel chimique dans les deux états internes, $\mu_\uparrow = \mu_\downarrow = \mu$. La théorie BCS statique approxime son état par un opérateur densité à N corps gaussien le plus général, $\hat{\sigma}_{\text{sta}} \propto \exp(-\beta \hat{H}_{\text{quad}})$, où $\beta = 1/k_B T$ et où $\hat{H}_{\text{quad}}$, forme quadratique en les opérateurs de champ fermioniques $\hat{\psi}_\sigma(\mathbf{r})$ et leurs hermitiens conjugués, contient non seulement les termes normaux $\hat{\psi}_\sigma^\dagger \hat{\psi}_\sigma$ mais aussi les termes anormaux $\hat{\psi}_\uparrow^\dagger \hat{\psi}_\downarrow$ et $\hat{\psi}_\downarrow^\dagger \hat{\psi}_\uparrow$ décrivant la condensation par paires $\uparrow\downarrow$.<fnref n="4" />,<fnref n="5" /> Après optimisation, on trouve que $\hat{H}_{\text{quad}} = \hat{H}_{\text{sta}}$ où l'hamiltonien BCS $\hat{H}_{\text{sta}}$ se déduit de l'hamiltonien grand-canonique complet $\hat{H}_{\text{GC}}$ en remplaçant le terme quartique d'interaction par des termes quadratiques au moyen de contractions de Wick incomplètes (voir la section Hartree-Fock-Bogolioubov dans la référence [16]; le cas limite $T = 0$ est traité en détail dans la section 3.2 de la référence [17]).

Cet hamiltonien quadratique est diagonalisé par la transformation de Bogolioubov, ici un développement des opérateurs de champ ou de leurs hermitiens conjugués sur les modes propres des quasi-particules BCS :

$$ \left( \begin{pmatrix} \hat{\psi}_{\uparrow} \\ \hat{\psi}_{\downarrow} \end{pmatrix} \begin{pmatrix} \mathbf{r} \\ \mathbf{r} \end{pmatrix} \right) = \frac{1}{L^{d/2}} \sum_{\mathbf{k}} \hat{\gamma}_{\mathbf{k}\uparrow} \begin{pmatrix} U_{\mathbf{k}} \\ V_{\mathbf{k}} \end{pmatrix} \exp(i\mathbf{k}\cdot\mathbf{r}) + \gamma_{\mathbf{k}\downarrow}^{\dagger} \begin{pmatrix} -V_{\mathbf{k}} \\ U_{\mathbf{k}} \end{pmatrix} \exp(-i\mathbf{k}\cdot\mathbf{r}) \quad \Rightarrow \quad \hat{H}_{\text{sta}} = E_0 + \sum_{\mathbf{k},\sigma} \epsilon_{\mathbf{k}} \hat{\gamma}_{\mathbf{k}\sigma}^{\dagger} \hat{\gamma}_{\mathbf{k}\sigma} \quad (13) $$

<footnote n="3">³Le cas $d=1$ est traité en champ moyen et au-delà par exemple dans la référence [14], et relève en tout cas de la théorie des systèmes intégrables par ansatz de Bethe [15].</footnote>

<footnote n="4">⁴En l'absence (i) de couplage cohérent de Rabi entre $\uparrow$ et $\downarrow$ et (ii) d'appariement dans l'onde $p$ entre mêmes états de spin, il est inutile d'inclure des termes (i) en $\hat{\psi}_\sigma^\dagger \hat{\psi}_{-\sigma}$ et (ii) en $\hat{\psi}_\sigma \hat{\psi}_\sigma$.</footnote>

<footnote n="5">⁵À 2D, il n'y a en réalité pas de condensation des paires liées à la limite thermodynamique (sauf si $T=0$), mais cette limitation de la théorie BCS statique importe peu ici car les fonctions de distribution de paires $g_{\sigma\sigma'}$ sont des observables locales – comprendre, insensibles à la portée infinie ou pas de l'ordre non diagonal. Voir aussi la référence [14], le même problème se posant à 1D (cette fois y compris à $T=0$).</footnote>
Les opérateurs de création $\hat{\gamma}_{k\sigma}^{\dagger}$ et d'annihilation $\hat{\gamma}_{k\sigma}$ d'une quasi-particule BCS de vecteur d'onde $\mathbf{k}$, de spin $\sigma$ et d'énergie propre $\epsilon_{\mathbf{k}}$, qui sont les coefficients de ce développement, obéissent aux relations d'anticommutation fermioniques discrètes habituelles, par exemple $\{\hat{\gamma}_{k\sigma}, \hat{\gamma}_{k'\sigma'}^{\dagger}\} = \delta_{\sigma\sigma'}\delta_{kk'}$ avec $\delta$ le delta de Kronecker. L'expression de l'énergie $E_0$ du vide de quasi-particules n'est pas requise ici. On trouve pour les amplitudes modales et les énergies propres

$$ U_{\mathbf{k}} = \frac{1}{\sqrt{2}} \left( 1 + \frac{\xi_{\mathbf{k}}}{\epsilon_{\mathbf{k}}} \right)^{1/2} ; \quad V_{\mathbf{k}} = \frac{1}{\sqrt{2}} \left( 1 - \frac{\xi_{\mathbf{k}}}{\epsilon_{\mathbf{k}}} \right)^{1/2} ; \quad \epsilon_{\mathbf{k}} = (\xi_{\mathbf{k}}^2 + \Delta^2)^{1/2} \quad (14) $$

avec les raccourcis de notation habituels

$$ \xi_{\mathbf{k}} = E_{\mathbf{k}} - \mu \quad \text{et} \quad E_{\mathbf{k}} = \frac{\hbar^2 k^2}{2m} \quad (15) $$

Comme annoncé, le spectre d'excitation $\epsilon_{\mathbf{k}}$ présente une bande d'énergie interdite $[0, E_{\text{gap}}]$ avec $E_{\text{gap}} = \Delta$ si $\mu \ge 0$, $E_{\text{gap}} = (\mu^2 + \Delta^2)^{1/2}$ si $\mu < 0$. Le paramètre d'ordre $\Delta$, pris ici réel positif, satisfait selon la dimensionnalité de l'espace à

$$ -\frac{m}{4\pi\hbar^2 a_{3D}} \stackrel{d=3}{=} \int_{\mathbb{R}^3} \frac{d^3 k}{(2\pi)^3} \left( \frac{1-2f_{\mathbf{k}}}{2\epsilon_{\mathbf{k}}} - \frac{1}{2E_{\mathbf{k}}} \right) ; \quad 0 \stackrel{d=2}{=} \int_{\mathbb{R}^2} \frac{d^2 k}{(2\pi)^2} \left( \frac{1-2f_{\mathbf{k}}}{2\epsilon_{\mathbf{k}}} - \frac{1}{2E_{\mathbf{k}} + \epsilon_{\text{dim}}} \right) \quad (16) $$

où la loi de Fermi-Dirac $f_{\mathbf{k}} = 1/[\exp(\beta\epsilon_{\mathbf{k}}) + 1]$ donne le nombre thermique moyen $\langle \hat{\gamma}_{k\sigma}^{\dagger} \hat{\gamma}_{k\sigma} \rangle_{\text{sta}}$ de quasi-particules fermioniques $k\sigma$ et l'on a fait usage de l'identité $U_{\mathbf{k}}V_{\mathbf{k}} = \Delta/2\epsilon_{\mathbf{k}}$.<fnref n="6" />

Puisque l'opérateur densité variationnel $\hat{\sigma}_{\text{sta}}$ est gaussien, le théorème de Wick s'applique et les moyennes d'observables dans la théorie BCS statique – repérées par l'indice « sta » – se déduisent des moyennes bilinéaires suivantes, dites anormales ou anormales selon qu'elles brisent la symétrie $U(1)$ ou pas :

$$ \langle \hat{\psi}_{\downarrow}(\mathbf{r}', t') \psi_{\downarrow}(\mathbf{r}, t) \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} U_{\mathbf{k}} V_{\mathbf{k}} [f_{\mathbf{k}} \exp(-i\phi_{\mathbf{k}}) - (1-f_{\mathbf{k}}) \exp(i\phi_{\mathbf{k}})] \exp[i\mathbf{k} \cdot (\mathbf{r}-\mathbf{r}')] \quad (17) $$

$$ \langle \psi_{\downarrow}(\mathbf{r}, t) \hat{\psi}_{\downarrow}(\mathbf{r}', t') \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} U_{\mathbf{k}} V_{\mathbf{k}} [(1-f_{\mathbf{k}}) \exp(-i\phi_{\mathbf{k}}) - f_{\mathbf{k}} \exp(i\phi_{\mathbf{k}})] \exp[i\mathbf{k} \cdot (\mathbf{r}-\mathbf{r}')] \quad (18) $$

$$ \langle \hat{\psi}_{\uparrow}^{\dagger}(\mathbf{r}, t) \hat{\psi}_{\uparrow}(\mathbf{r}', t') \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} [U_{\mathbf{k}}^2 f_{\mathbf{k}} \exp(i\phi_{\mathbf{k}}) + V_{\mathbf{k}}^2 (1-f_{\mathbf{k}}) \exp(-i\phi_{\mathbf{k}})] \exp[i\mathbf{k} \cdot (\mathbf{r}-\mathbf{r}')] \quad (19) $$

$$ \langle \hat{\psi}_{\uparrow}(\mathbf{r}, t) \hat{\psi}_{\uparrow}^{\dagger}(\mathbf{r}', t') \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} [U_{\mathbf{k}}^2 (1-f_{\mathbf{k}}) \exp(-i\phi_{\mathbf{k}}) + V_{\mathbf{k}}^2 f_{\mathbf{k}} \exp(i\phi_{\mathbf{k}})] \exp[i\mathbf{k} \cdot (\mathbf{r}-\mathbf{r}')] \quad (20) $$

où la dépendance en temps des opérateurs de champ est celle du point de vue de Heisenberg pour l'hamiltonien $\hat{H}_{\text{sta}}$, et l'on a posé $\phi_{\mathbf{k}} = \epsilon_{\mathbf{k}}(t-t')/\hbar$ pour abréger. Faut-il le préciser, pour notre gaz non polarisé sans couplage de Rabi, les moyennes normales (19,20) ne dépendent pas de l'état interne $\sigma = \uparrow, \downarrow$ choisi et les moyennes normales croisées (entre deux états de spin différents) sont nulles.

En spécialisant la relation (19) à des positions et à des temps égaux, on obtient l'équation d'état de la theorie BCS statique, reliant la densité moyenne totale du gaz $\rho^{\text{sta}} = \rho_{\uparrow}^{\text{sta}} + \rho_{\downarrow}^{\text{sta}}$ au potentiel chimique $\mu$ et à la température $T$,

$$ \rho^{\text{sta}} = 2\rho_{\sigma}^{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left[ 1 - (1 - 2f_{\mathbf{k}}) \frac{\xi_{\mathbf{k}}}{\epsilon_{\mathbf{k}}} \right] \quad (21) $$

<footnote n="6">⁶Nous ne spécifions pas le modèle microscopique de l'interaction, en pratique d'une portée $b$ non nulle qu'on fait tendre vers zéro à la fin des calculs. Une possibilité est le modèle sur réseau $\mathbf{r} \in b\mathbb{Z}^d$ de la référence [18], avec interaction sur site $V(\mathbf{r}_i, \mathbf{r}_j) = g_0\delta_{\mathbf{r}_i, \mathbf{r}_j}/b^d$ et relation de dispersion libre $\mathbf{k} \in \mathcal{D} \mapsto E_{\mathbf{k}}(\mathcal{D} = [-\pi/b, \pi/b]^d$ est la première zone de Brillouin). Dans ce cas, on a simplement $\Delta = g_0\langle\hat{\psi}_{\downarrow}(\mathbf{r})\hat{\psi}_{\uparrow}(\mathbf{r}')\rangle_{\text{sta}}$; en ajustant la constante de couplage nue $g_0$ pour reproduire la longueur de diffusion $a_{\text{D}}$ (ou plus commodément à 2D l'énergie de liaison du dimère $\epsilon_{\text{dim}}$), on arrive à l'équation (16).</footnote>
On introduit par la même occasion le nombre d'onde de Fermi $k_F^{\text{sta}}$ du gaz parfait non polarisé de température nulle mais de même densité que le gaz en interaction :

$$k_F^{d=3} = (3\pi^2\rho)^{1/3} \quad ; \quad k_F^{d=2} = (2\pi\rho)^{1/2} \qquad (22)$$

ce qui permet de définir l'énergie et la température de Fermi, $E_F = k_B T_F = \hbar^2 k_F^2 / 2m$. Enfin, donnons les fonctions de distribution de paires et leurs écarts à leurs valeurs asymptotiques dans la théorie BCS statique, dont nous avons souligné les insuffisances dans la section 1 :

$$\delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') \equiv g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') - \rho_{\uparrow}^{\text{sta}} \rho_{\downarrow}^{\text{sta}} = +|\langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}}|^2 \qquad (23)$$

$$\delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') \equiv g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') - \rho_{\uparrow}^{\text{sta}} \rho_{\downarrow}^{\text{sta}} = -|\langle \hat{\psi}_{\uparrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\downarrow}(\mathbf{r}') \rangle_{\text{sta}}|^2 \qquad (24)$$

avec les moyennes anormale et normale tirées des équations (17,19) prises à temps égaux :

$$\langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{(-\Delta)}{2\epsilon_k} (1 - 2f_k) \exp[i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')] \qquad (25)$$

$$\langle \hat{\psi}_{\uparrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\downarrow}(\mathbf{r}') \rangle_{\text{sta}} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{1}{2} \left[ 1 - (1 - 2f_k) \frac{\xi_k}{\epsilon_k} \right] \exp[i\mathbf{k} \cdot (\mathbf{r} - \mathbf{r}')] \qquad (26)$$

Signalons pour terminer qu'à une température $T \ll E_{\text{gap}}/k_B$, on peut négliger la dissociation thermique des paires liées $\uparrow\downarrow$, donc se placer à l'ordre zéro en le petit paramètre $\exp(-\beta E_{\text{gap}})$ sans que cela soit très contraignant physiquement vu sa décroissance exponentielle, et omettre les nombres d'occupation thermiques $f_k$. Nous verrons dans la section 3.2 que cette approximation est justifiée dans les gaz d'atomes froids aux températures minimales accessibles. Des formes analytiques peuvent alors être obtenues pour l'équation sur le paramètre d'ordre (16) et l'équation d'état (21). En dimension $d=3$, elles font intervenir des intégrales elliptiques [19]; inspiré par la référence [20], nous les simplifions quelque peu au moyen d'un paramétrage hyperbolique $\tau = \text{argsh}(\mu/\Delta)$ et des identités $K(ix) = K(x/\sqrt{x^2+1})/\sqrt{x^2+1}$ et $E(ix) = E(x/\sqrt{1+x^2})\sqrt{1+x^2}$ $(\forall x \in \mathbb{R})$ sur les intégrales elliptiques complètes de première espèce $K(k) = F(\pi/2, k)$ et de seconde espèce $E(k) = E(\pi/2, k)$ [21] :

$$-\frac{\pi}{4a_{3D}} \left( \frac{\hbar^2}{m\Delta} \right)^{1/2} \stackrel{d=3}{=} e^{+\tau/2} [\text{ch}\tau K(ie^\tau) - e^{-\tau} E(ie^\tau)] \qquad (27)$$

$$\frac{3}{2} \pi^2 \rho^{\text{sta}} \left( \frac{\hbar^2}{m\Delta} \right)^{3/2} \stackrel{d=3}{=} e^{-\tau/2} [\text{ch}\tau K(ie^\tau) + \text{sh}\tau E(ie^\tau)] \qquad (28)$$

En dimension $d=2$, en revanche, elles sont explicites [22] (il suffit d'utiliser $\xi_k$ comme variable d'intégration radiale) :

$$\begin{gathered} 0 \stackrel{d=2}{=} \text{argsh} \frac{\mu}{\Delta} + \ln \frac{\epsilon_{\text{dim}}}{\Delta} \quad \text{et} \quad \rho^{\text{sta}} \stackrel{d=2}{=} \frac{m}{2\pi\hbar^2} (\mu + \sqrt{\Delta^2 + \mu^2}) \\ \Rightarrow \quad \Delta \stackrel{d=2}{=} [\epsilon_{\text{dim}}(2\mu + \epsilon_{\text{dim}})]^{1/2} = \frac{\hbar^2}{m} k_F^{\text{sta}} q_{\text{dim}} \quad \text{et} \quad \rho^{\text{sta}} \stackrel{d=2}{=} \frac{m}{2\pi\hbar^2} (2\mu + \epsilon_{\text{dim}}) \end{gathered} \qquad (29)$$

Dans ce cas, la référence [1] parvient à exprimer les moyennes (25,26) en termes de fonctions de Bessel, ce qui est remarquable :

$$\langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} \stackrel{d=2}{=} -\frac{m\Delta}{2\pi\hbar^2} J_0(k_F^{\text{sta}}|\mathbf{r}-\mathbf{r}'|) K_0(q_{\text{dim}}|\mathbf{r}-\mathbf{r}'|) \qquad (30)$$

$$\langle \hat{\psi}_{\uparrow}^{\dagger}(\mathbf{r}) \hat{\psi}_{\downarrow}(\mathbf{r}') \rangle_{\text{sta}} \stackrel{d=2}{=} \frac{m\Delta}{2\pi\hbar^2} J_1(k_F^{\text{sta}}|\mathbf{r}-\mathbf{r}'|) K_1(q_{\text{dim}}|\mathbf{r}-\mathbf{r}'|) \qquad (31)$$

Mentionnons un intérêt particulier de l'équation d'état à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ pour notre problème : dans la théorie BCS dépendant du temps linéarisée, la pente $c$ (vitesse du son à
température nulle) de la relation de dispersion (10) des phonons, dont la section 1 a souligné l'importance, souscrit à la relation exacte de l'hydrodynamique des superfluides

$$mc^2(T=0) = \rho \frac{d\mu}{d\rho} \qquad (32)$$

pour l'équation d'état *approchée* $\rho^{\text{sta}}(\mu)$ de la théorie BCS statique, à cause de la nature variationnelle non perturbative de cette théorie [23]. D'où les expressions de la référence [19] à 3D (là aussi simplifiée par nos soins) et des références [19, 24] à 2D :

$$mc^2 \stackrel{d=3}{=} \frac{2}{3} \Delta \frac{[E(ie^T) - K(ie^T)] [\text{ch} \tau K(ie^T) + \text{sh} \tau E(ie^T)]}{|E(ie^T) + (ie^T - 1) K(ie^T)|^2} = \frac{2}{3} \mu \frac{1 + \frac{2\Delta^2}{3\mu\rho^{\text{sta}}g_{3D}}}{1 + \left(\frac{2\Delta}{3\rho^{\text{sta}}g_{3D}}\right)^2} \quad (33)$$

$$mc^2 \stackrel{d=2}{=} \frac{1}{2} \left( \mu + \sqrt{\Delta^2 + \mu^2} \right) = \mu + \frac{1}{2} \epsilon_{\text{dim}} = \frac{\pi \hbar^2}{m} \rho^{\text{sta}} = E_F^{\text{sta}} \quad (34)$$

L'expression au troisième membre de (33), qui fait intervenir la constante de couplage effective $g_{3D} = 4\pi\hbar^2 a_{3D}/m$, est absente de la référence [19]; nous l'avons obtenue en éliminant les intégrales elliptiques par résolution des équations (27,28) considérées formellement comme un système linéaire sur les inconnues $E(ie^T)$ et $K(ie^T)$. Dans la limite CBE où, à 3D, $a_{3D} > 0$ est fixé, $\mu \to -\hbar^2/2ma_{3D}^2$, $\rho^{\text{sta}} \to 0^+$ et $\Delta^2/\rho^{\text{sta}}g_{3D} \to \hbar^2/ma_{3D}^2$ [25], elle permet de vérifier facilement que $mc^2 \sim \rho^{\text{sta}}g_{3D}/4$, ce qui est de la forme attendue $m_d c^2 = \rho_d g_{3D}^{dd}$ pour un gaz bosonique de dimères $\uparrow\downarrow$ en interaction faible de masse $m_d$, de densité $\rho_d$ et de constante de couplage effective $g_{3D}^{dd}$; la valeur de $g_{3D}^{dd}$ prédite est cependant assez mauvaise car la théorie BCS statique ignore la branche d'excitation acoustique du superfluide.<fnref n="7" /> À 2D, en revanche, l'expression (34) donc l'équation d'état (29) sont catastrophiquement fausses dans cette limite, car elles ne reproduisent pas la bonne loi d'échelle $m_d c^2 \sim 2\pi\hbar^2\rho_d/[m_d \ln(\rho_d^{1/2} a_{2D}^{dd})]$ [29, 30], avec $a_{2D}^{dd}$ la longueur de diffusion dimère-dimère. Cette insuffisance est due là encore à l'omission des phonons par la théorie BCS statique; si l'on en tient compte, on retrouve le logarithme de la densité attendu [31].<fnref n="8" />

Signalons pour terminer que la température de transition $T_c$ prédite à 3D par la théorie BCS statique, solution de l'équation (16) pour $\Delta = 0$ donc pour $\epsilon_k = |\xi_k|$, est à peu près correcte $\propto \Delta(T=0)/k_B$ dans la limite BCS, du bon ordre de grandeur $\propto T_F/k_B$ à la limite unitaire, mais catastrophiquement fausse dans la limite CBE où elle reproduit – à un facteur près – la température de dissociation $\propto \hbar^2/[ma_{3D}^2 k_B \ln(k_F a_{3D})]$ des dimères dans le gaz (le logarithme est une contribution entropique favorisant la brisure de paires) plutôt que la température de transition $\propto T_F$ du gaz parfait de dimères bosoniques [25].<fnref n="9" /> Ce sera sans grande pertinence ici car nous nous limiterons en pratique à des températures $T \ll T_F$, $\Delta/k_B$ et nous nous mettrons à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ dans nos applications et résultats explicites, voir la section 3. La référence [32] montre comment aller au-delà de la théorie BCS statique dans le calcul de $T_c$;

<footnote n="7">⁷Par définition, on a $g_{3D}^{dd} = 4\pi\hbar^2 a_{3D}^{dd}/m_d$ avec $a_{3D}^{dd}$ la longueur de diffusion dimère-dimère tridimensionnelle. La théorie BCS statique prédit donc $a_{3D}^{dd} = 2a_{3D}$, ce qui est assez éloigné de la valeur exacte $a_{3D}^{dd} \approx 0,60a_{3D}$ obtenue par résolution directe de l'équation de Schrödinger pour quatre fermions $\uparrow\uparrow\uparrow\uparrow$ en interaction de contact [26] ou par une méthode diagrammatique de la théorie quantique des champs [27]. L'amélioration de la théorie BCS statique au moyen de l'approximation des fluctuations gaussiennes, qui tient compte des modes de phonons, voir notre note 20, donne le résultat bien meilleur $a_{3D}^{dd} \approx 0,55a_{3D}$ [28].</footnote>

<footnote n="8">⁸L'approximation des fluctuations gaussiennes, voir notre note 20, reproduit le logarithmique avec une longueur de diffusion dimère-dimère bidimensionnelle $a_{2D}^{dd} \approx 0,56a_{2D}$.</footnote>

<footnote n="9">⁹Introduisons l'énergie de liaison du dimère tridimensionnel $\epsilon_d = \hbar^2/ma_{3D}^2$ pour abréger. Alors, comme le dit la référence [25], l'équation (16) montre que $x_c \equiv \epsilon_d/k_B T_c^{\text{sta}} \to +\infty$ dans la limite CBE et le relief à $\delta\mu \equiv \mu + \epsilon_d/2 \to 0^+$; l'équation d'état (21) écrite pour $\Delta = 0$ relie ensuite $x_c$ à la densité. On trouve que $\int_{\mathbb{R}^3} [d^3k/(2\pi)^3] [2\epsilon_d\delta\mu/(2E_k + \epsilon_d)]^2 \sim \int_{\mathbb{R}^3} [d^3k/(2\pi)^3] 2f_k|\Delta=0 = \rho^{\text{sta}}|_{\Delta=0}$ donc en définitive que $(\pi/4)\delta\mu/\epsilon_d \sim (\pi/2)^{1/2} \exp(-x_c/2)/x_c^{3/2} \sim \pi^2 \rho^{\text{sta}} a_{3D}^3$, d'où la loi d'échelle annoncée sur $T_c^{\text{sta}}$.</footnote>
en particulier, elle retrouve exactement la valeur attendue dans la limite CBE. Le point clé est encore et toujours la prise en compte de la branche d'excitation acoustique $\hbar\omega_q$ ignorée par BCS statique mais qui, dans la limite CBE à nombre d'onde $q > 0$ fixé, se réduit triviallement à la mise en mouvement du centre de masse d'une paire liée avec un vecteur d'onde $\mathbf{q}$ et une énergie purement cinétique $\hbar^2 q^2 / 4m$, donc reproduit le spectre d'excitation d'un condensat de bosons sans interaction de masse $2m$ et la fonction de partition associée.

## 2.2. L'obtention des $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ en elle-même

Cette obtention résulte d'un enchaînement élémentaire de résultats connus, ne reposant au départ sur aucune approximation.

Premièrement, le facteur de structure dynamique $S_{\sigma\sigma'}(\mathbf{q}, \omega)$ du gaz à l'équilibre thermique, ici un tenseur à deux indices puisque nos fermions sont de spin non nul, est par définition la transformée de Fourier (au vecteur d'onde $\mathbf{q}$ et à la pulsation $\omega$) de la fonction de corrélation spatio-temporelle de la densité, si bien que

$$ \langle \hat{\rho}_{\sigma}(\mathbf{r}, t) \hat{\rho}_{\sigma'}(\mathbf{r}', t') \rangle - \langle \hat{\rho}_{\sigma}(\mathbf{r}, t) \rangle \langle \hat{\rho}_{\sigma'}(\mathbf{r}', t') \rangle = \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \int_{\mathbb{R}} \frac{d\omega}{2\pi} \exp[i\{\mathbf{q} \cdot (\mathbf{r} - \mathbf{r}') - \omega(t - t')\}] S_{\sigma\sigma'}(\mathbf{q}, \omega) \quad (35) $$

où la moyenne est prise dans l'état thermique exact, l'évolution de Heisenberg correspond à l'hamiltonien exact et l'on a bien entendu $\hat{\rho}_{\sigma}(\mathbf{r}, t) = \hat{\psi}_{\sigma}^{\dagger}(\mathbf{r}, t)\hat{\psi}_{\sigma}(\mathbf{r}, t)$ en seconde quantification. En spécialisant l'équation (35) aux temps égaux et en utilisant les relations d'anticommutation fermioniques sur les opérateurs de champ pour les mettre dans l'ordre normal comme dans la définition (1) des fonctions de distribution de paires, on obtient

$$ \delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}') = -\rho_{\sigma}\delta(\mathbf{r}-\mathbf{r}') + \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \int_{\mathbb{R}} \frac{d\omega}{2\pi} \exp[i\mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')] S_{\sigma\sigma'}(\mathbf{q}, \omega) \quad (36) $$

Deuxièmement, le théorème de fluctuation-dissipation (voir la section §124 de la référence [33]), dit aussi théorème de Callen-Welton [34], relie le facteur de structure dynamique à la susceptibilité $\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, \hbar\omega)$, elle-même transformée de Fourier de la fonction de réponse linéaire densité-densité,<fnref n="10" /> comme suit :

$$ \operatorname{Im}[\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, \hbar\omega + i0^+)] = -\frac{1}{2\hbar} S_{\sigma\sigma'}(\mathbf{q}, \omega) [1 - \exp(-\beta\hbar\omega)] \quad (39) $$

$$ \chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t') = \langle [\hat{\rho}_{\sigma}(\mathbf{r}, t), \hat{\rho}_{\sigma'}(\mathbf{r}', t')] \rangle \Theta(t-t') \quad (37) $$

Dans cette formule, la moyenne et l'évolution de Heisenberg sont celles du système non perturbé et $\Theta$ est la fonction de Heaviside. Comme le gaz est initialement spatialement homogène et à l'équilibre, $\chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t')$ est une fonction seulement des variables relatives $\mathbf{r}-\mathbf{r}'$ et $t-t'$, dont la susceptibilité est la transformée de Fourier régularisée comme suit :

$$ \tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, z) = \int_{\mathbb{R}^d} d^d\mathbf{r} \int_{\mathbb{R}} dt \chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t') e^{i[\mathbf{z}(t-t')/\hbar - \mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')]} \quad (38) $$

avec $z = \hbar\omega + i\eta$, $\eta \to 0^+$. Le facteur $\exp[-\eta(t-t')]$ ainsi introduit assure la convergence de l'intégrale sur le temps. L'expression (38) permet d'étendre $\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, z)$ au demi-plan complexe supérieur $\operatorname{Im} z > 0$, et la propriété $[\chi_{\sigma\sigma'}(\mathbf{q}, z)]^* = \chi_{\sigma\sigma'}(\mathbf{q}, z^*)$ au demi-plan inférieur. La fonction obtenue est analytique sur $C \setminus \mathbb{R}$ mais admet des singularités (pôles, lignes de coupure) sur l'axe réel.

<footnote n="10">¹⁰L'application d'un potentiel extérieur infinitésimal $U$ dépendant du temps dans l'état interne $\sigma'$ induit dans l'état interne $\sigma$, à l'instant $t$ et au point $\mathbf{r}$, une modification de densité moyenne $\delta\rho_{\sigma}(\mathbf{r}, t) = \int d^d\mathbf{r}' dt' \chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t') U(\mathbf{r}', t')$, où la fonction de réponse est donnée par la formule de Kubo,</footnote>
Il en résulte que<fnref n="11" />, après le changement de variable $\omega = \epsilon/\hbar$:

$$ \delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}') = -\rho_{\sigma}\delta(\mathbf{r}-\mathbf{r}') + \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \exp[i\mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')] \int_{\mathbb{R}} \frac{d\epsilon}{(-\pi)} \frac{\operatorname{Im}[\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, \epsilon + i0^+)]}{1 - \exp(-\beta\epsilon)} \quad (40) $$

Troisièmement, les susceptibilités densité-densité ont déjà été calculées de façon approchée au moyen de la théorie BCS dépendant du temps linéarisée autour de son état stationnaire.<fnref n="12" /> En généralisant à température non nulle (comme dans la référence [36]<fnref n="13" />) les résultats des références [20, 37, 38], nous trouvons au vecteur d'onde $\mathbf{q}$ et à l'énergie complexe $z$ que

$$ \tilde{\chi}_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{q}, z) = \frac{1}{2}(\Sigma_{33} - \Sigma_{44}) + \frac{2\Sigma_{13}\Sigma_{23}\Sigma_{12} - \Sigma_{13}^2\Sigma_{22} - \Sigma_{23}^2\Sigma_{11}}{2(\Sigma_{11}\Sigma_{22} - \Sigma_{12}^2)} \quad (41) $$

$$ \tilde{\chi}_{\uparrow\uparrow}^{\text{dtl}}(\mathbf{q}, z) = \frac{1}{2}(\Sigma_{33} + \Sigma_{44}) + \frac{2\Sigma_{13}\Sigma_{23}\Sigma_{12} - \Sigma_{13}^2\Sigma_{22} - \Sigma_{23}^2\Sigma_{11}}{2(\Sigma_{11}\Sigma_{22} - \Sigma_{12}^2)} \quad (42) $$

où l'exposant « dtl » est mis pour « théorie BCS dépendant du temps linéarisée ». Nous avons introduit au second membre les fonctions suivantes, en sous-entendant leur dépendance en les variables $(\mathbf{q}, z) \in \mathbb{R}^d \times (C \setminus \mathbb{R})$ pour alléger :

$$ \Sigma_{11} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{(\epsilon_+ + \epsilon_-)(\epsilon_+ \epsilon_- + \xi_+ \xi_- + \Delta^2)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{(\epsilon_+ - \epsilon_-)(-\epsilon_+ \epsilon_- + \xi_+ \xi_- + \Delta^2)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} + \frac{1}{2\epsilon_k} \right\} \quad (43) $$

$$ \Sigma_{22} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{(\epsilon_+ + \epsilon_-)(\epsilon_+ \epsilon_- + \xi_+ \xi_- - \Delta^2)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{(\epsilon_+ - \epsilon_-)(-\epsilon_+ \epsilon_- + \xi_+ \xi_- - \Delta^2)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} + \frac{1}{2\epsilon_k} \right\} \quad (44) $$

$$ \Sigma_{33} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{(\epsilon_+ + \epsilon_-)(\epsilon_+ \epsilon_- - \xi_+ \xi_- + \Delta^2)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{(\epsilon_+ - \epsilon_-)(-\epsilon_+ \epsilon_- - \xi_+ \xi_- + \Delta^2)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} \right\} \quad (45) $$

$$ \Sigma_{44} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{(\epsilon_+ + \epsilon_-)(\epsilon_+ \epsilon_- - \xi_+ \xi_- - \Delta^2)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{(\epsilon_+ - \epsilon_-)(-\epsilon_+ \epsilon_- - \xi_+ \xi_- - \Delta^2)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} \right\} \quad (46) $$

$$ \Sigma_{12} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{z(\epsilon_+ \xi_- + \epsilon_- \xi_+)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{z(\epsilon_+ \xi_- - \epsilon_- \xi_+)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} \right\} \quad (47) $$

$$ \Sigma_{13} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{z\Delta(\epsilon_+ + \epsilon_-)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{z\Delta(\epsilon_+ - \epsilon_-)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} \right\} \quad (48) $$

$$ \Sigma_{23} = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{\Delta(\epsilon_+ + \epsilon_-)(\xi_+ + \xi_-)(1 - f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ + \epsilon_-)^2]} + \frac{\Delta(\epsilon_+ - \epsilon_-)(\xi_+ + \xi_-)(f_+ - f_-)}{2\epsilon_+ \epsilon_- [z^2 - (\epsilon_+ - \epsilon_-)^2]} \right\} \quad (49) $$

<footnote n="11">¹¹Pour un système de taille finie, il faudrait en général tenir compte, dans le facteur de structure dynamique, d'une contribution en $\delta(\omega)$ non contrainte par (39), due en particulier à un écart entre la susceptibilité thermodynamique (qui traduit le changement de densité à l'équilibre thermique en présence d'un réservoir à la température $T$) et la limite de pulsation nulle de la susceptibilité dynamique $\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, \hbar\omega + i0^+)$ (qui traduit la réponse d'un système isolé purement hamiltonien), voir le chapitre 8 de la référence [35], l'équation (66) de la référence [14] et son annexe C.</footnote>

<footnote n="12">¹²La littérature donne plutôt les susceptibilités sur la densité totale et sur la densité différentielle. Celles qui nous intéressent s'en déduisent par bilinéarité et invariance par échange de $\uparrow$ et $\downarrow$: $\chi_{\hat{\rho}_1+\epsilon\hat{\rho}_1, \hat{\rho}_1+\epsilon\hat{\rho}_1} = \chi_{\uparrow\uparrow} + \epsilon\chi_{\uparrow\downarrow} + \epsilon\chi_{\downarrow\uparrow} = 2(\chi_{\uparrow\uparrow} + \epsilon\chi_{\uparrow\downarrow})$.</footnote>

<footnote n="13">¹³Les susceptibilités de la référence [36], voir ses équations (43,46), extensives plutôt qu'intensives, différent des nôtres par un facteur $L^d/2$ pour la réponse densité totale-densité totale (notée $\tilde{\chi}_{33}$ dans [36]), et par un facteur $-L^d/2$ pour la réponse densité différentielle-densité différentielle (notée $\tilde{\chi}_{44}$ dans [36]). Signalons que les références [20,37] ont fait l'impasse sur la densité différentielle.</footnote>
avec les notations très compactes<fnref n="14" />

$$ \mathbf{k}_{\pm} = \mathbf{k}_{\pm} \mathbf{q}/2 \quad ; \quad \xi_{\pm} = \xi_{\mathbf{k}_{\pm}} \quad ; \quad \epsilon_{\pm} = \epsilon_{\mathbf{k}_{\pm}} \quad ; \quad f_{\pm} = f_{\mathbf{k}_{\pm}} \qquad (50) $$

On notera que les fonctions $\Sigma_{ij}$ sont de parité bien définie pour la variable $z$ (elles sont paires, sauf $\Sigma_{12}$ et $\Sigma_{13}$ qui sont impaires), et qu'elles ne sont pas indépendantes, au sens où

$$ \Sigma_{11}(\mathbf{q}, z) - \Sigma_{22}(\mathbf{q}, z) = \Sigma_{33}(\mathbf{q}, z) - \Sigma_{44}(\mathbf{q}, z) \quad \text{et} \quad \Sigma_{13}(\mathbf{q}, z) = \frac{z}{2\Delta} [\Sigma_{11}(\mathbf{q}, z) - \Sigma_{22}(\mathbf{q}, z)] \quad \forall (\mathbf{q}, z) \quad (51) $$

On a bien entendu $\tilde{\chi}_{\uparrow\downarrow}^{\text{dtl}} \equiv \tilde{\chi}_{\uparrow\downarrow}^{\text{dtl}}$ et $\tilde{\chi}_{\downarrow\uparrow}^{\text{dtl}} \equiv \tilde{\chi}_{\downarrow\uparrow}^{\text{dtl}}$ par symétrie entre les deux états de spin. Toutes ces susceptibilités sont des fonctions paires de $z$.

Le report des expressions (41,42) dans l'équation (40) conduit aux fonctions $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ de la théorie BCS dépendant du temps linéarisée. Nous les étudions en détail dans la suite.

## 2.3. Isoler la partie collective dans les $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$

Il se trouve que les premières contributions au second membre de (41,42), avant les grosses fractions, ne sont autres que les susceptibilités de la théorie BCS statique. Pour le voir, il suffit d'appliquer la formule de Kubo (37) à cette théorie puis de faire directement usage du théorème de Wick et des fonctions de corrélation à deux temps (17,18,19,20), sans chercher à calculer le commutateur. Des manipulations élémentaires<fnref n="15" /> conduisent à

$$ \chi_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, t; \mathbf{r}', t') = \frac{1}{i\hbar} \int_{\mathbb{R}^d} \frac{d^d k_1}{(2\pi)^d} \int_{\mathbb{R}^d} \frac{d^d k_2}{(2\pi)^d} U_{\mathbf{k}_1} V_{\mathbf{k}_1} U_{\mathbf{k}_2} V_{\mathbf{k}_2} \left\{ \left[ e^{i(\phi_{\mathbf{k}_1} + \phi_{\mathbf{k}_2})} - e^{-i(\phi_{\mathbf{k}_1} + \phi_{\mathbf{k}_2})} \right] (f_{\mathbf{k}_1} + f_{\mathbf{k}_2} - 1) \right. \\ \left. + \left[ e^{i(\phi_{\mathbf{k}_1} - \phi_{\mathbf{k}_2})} - e^{-i(\phi_{\mathbf{k}_1} - \phi_{\mathbf{k}_2})} \right] (f_{\mathbf{k}_2} - f_{\mathbf{k}_1}) \right\} \exp[i(\mathbf{k}_1 - \mathbf{k}_2) \cdot (\mathbf{r} - \mathbf{r}')] \quad (52) $$

$$ \chi_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, t; \mathbf{r}', t') = \frac{1}{i\hbar} \int_{\mathbb{R}^d} \frac{d^d k_1}{(2\pi)^d} \int_{\mathbb{R}^d} \frac{d^d k_2}{(2\pi)^d} \left\{ \left[ e^{i(\phi_{\mathbf{k}_1} + \phi_{\mathbf{k}_2})} - e^{-i(\phi_{\mathbf{k}_1} + \phi_{\mathbf{k}_2})} \right] \frac{1}{2} (U_{\mathbf{k}_1}^2 V_{\mathbf{k}_1}^2 + U_{\mathbf{k}_2}^2 V_{\mathbf{k}_2}^2) (f_{\mathbf{k}_1} + f_{\mathbf{k}_2} - 1) \right. \\ \left. + \left[ e^{i(\phi_{\mathbf{k}_1} - \phi_{\mathbf{k}_2})} - e^{-i(\phi_{\mathbf{k}_1} - \phi_{\mathbf{k}_2})} \right] (f_{\mathbf{k}_1} - f_{\mathbf{k}_2}) \frac{1}{2} (U_{\mathbf{k}_1}^2 U_{\mathbf{k}_2}^2 + V_{\mathbf{k}_1}^2 V_{\mathbf{k}_2}^2) \right\} \exp[i(\mathbf{k}_1 - \mathbf{k}_2) \cdot (\mathbf{r} - \mathbf{r}')] \quad (53) $$

Il reste à prendre la transformée de Fourier spatio-temporelle comme dans l'équation (38) après le changement de variables de jacobien unité $\mathbf{k}_{1,2} = \mathbf{k} \pm \mathbf{q}/2$ pour trouver le résultat annoncé :

$$ \tilde{\chi}_{\uparrow\downarrow}^{\text{sta}}(\mathbf{q}, z) = \frac{1}{2} [\Sigma_{33}(\mathbf{q}, z) - \Sigma_{44}(\mathbf{q}, z)] \quad ; \quad \tilde{\chi}_{\uparrow\downarrow}^{\text{sta}}(\mathbf{q}, z) = \frac{1}{2} [\Sigma_{33}(\mathbf{q}, z) + \Sigma_{44}(\mathbf{q}, z)] \quad (54) $$

Ainsi, la grande fraction dans (41,42) – c’est exactement la même dans les deux cas – tient compte des excitations collectives (les phonons de Bogolioubov-Anderson) ignorées par la théorie BCS statique et par les références [1,2]. En vue des études qui vont suivre, isolons sa contribution aux fonctions de distribution de paires en posant (un exposant « dtl » serait redondant sur la contribution collective) :

$$ \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) = \frac{2\Sigma_{13}\Sigma_{23}\Sigma_{12} - \Sigma_{12}^2\Sigma_{22} - \Sigma_{23}^2\Sigma_{11}}{2(\Sigma_{11}\Sigma_{22} - \Sigma_{12}^2)} \quad \text{et} \quad \delta\tilde{g}_{\text{coll}}(\mathbf{q}) = \int_{\mathbb{R}^d} \frac{d\epsilon}{(-\pi)} \frac{\text{Im}[\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+)]}{1 - \exp(-\beta\epsilon)} \quad (55) $$

<footnote n="14">¹⁴En partant des définitions (14), nous avons éliminé les amplitudes modales $U_{\pm} = U_{\mathbf{k}_{\pm}}$ et $V_{\pm} = V_{\mathbf{k}_{\pm}}$ au moyen des relations $(U_{+}V_{-} + sU_{-}V_{+})^2 = (\epsilon_{+}\epsilon_{-} - \xi_{+}\xi_{-} + s\Delta^2)/2\epsilon_{+}\epsilon_{-}$, $(U_{+}U_{-} + sV_{+}V_{-})^2 = (\epsilon_{+}\epsilon_{-} + \xi_{+}\xi_{-} + s\Delta^2)/2\epsilon_{+}\epsilon_{-}$, $(U_{+}V_{-} + sU_{-}V_{+})(U_{+}U_{-} + sV_{+}V_{-}) = \Delta(\epsilon_{+} + s\epsilon_{-})/2\epsilon_{+}\epsilon_{-}$ et $(U_{+}V_{-} + sU_{-}V_{+})(U_{+}U_{-} - sV_{+}V_{-}) = \Delta(\xi_{+} + s\xi_{-})/2\epsilon_{+}\epsilon_{-}$ ($s \in \{-1; +1\}$). On remarque que, dans les fonctions $\Sigma_{ij}$, le second terme se déduit du premier si l'on fait le changement formel $\epsilon_{-} \to -\epsilon_{-}$ (qui a pour effet de changer $f_{-}$ en $1-f_{-}$) sans toucher ni à $\epsilon_{+}$ ni aux $\xi_{\pm}$.</footnote>

<footnote n="15">¹⁵L'invariance par parité permet de changer à volonté le signe des vecteurs d'onde sous les exponentielles de Fourier, donc en particulier – comme dans (53) – de symétrer l'intégrande dans les intégrales sur deux vecteurs d'onde au moyen des identités $\int d^d k_1 \int d^d k_2 f(\mathbf{k}_1, \mathbf{k}_2) \exp[i(\mathbf{k}_1 - \mathbf{k}_2) \cdot (\mathbf{r} - \mathbf{r}')] = \int d^d k_1 \int d^d k_2 f(\mathbf{k}_2, \mathbf{k}_1) \exp[i(\mathbf{k}_2 - \mathbf{k}_1) \cdot (\mathbf{r} - \mathbf{r}')] = \int d^d k_1 \int d^d k_2 (1/2)[f(\mathbf{k}_1, \mathbf{k}_2) + f(-\mathbf{k}_2, -\mathbf{k}_1)] \exp[i(\mathbf{k}_1 - \mathbf{k}_2) \cdot (\mathbf{r} - \mathbf{r}')]$. Dans notre cas, les amplitudes modales et les énergies propres, donc les fonctions $f$, non contentes d'être paires, dépendent seulement des modules des vecteurs d'onde.</footnote>
Il en résulte que

$$ \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') = \delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}') + \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') \quad (56) $$

avec la contribution collective commune

$$ \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') = \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) \exp[i\mathbf{q} \cdot (\mathbf{r} - \mathbf{r}')] \quad (57) $$

Nous en déduisons en particulier les belles relations naturellement équivalentes

$$ \delta g_{\uparrow\uparrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - \delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') = \delta g_{\uparrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') - \delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') \quad (58) $$

$$ \delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - \delta g_{\uparrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - \delta g_{\uparrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') \quad (59) $$

sur lesquelles nous reviendrons dans l'équation (75). Pour être complet, calculons la transformée de Fourier des $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ au sens de l'équation (57) en appliquant le théorème de fluctuation-dissipation – plus précisément l'équation (40) qui en résulte – aux susceptibilités (54) de la théorie BCS statique :<fnref n="16" />

$$ \delta \tilde{g}_{\uparrow\downarrow}^{\text{sta}}(\mathbf{q}) = + \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{4\epsilon_+\epsilon_-} (1 - 2f_+) (1 - 2f_-) \quad (60) $$

$$ \delta \tilde{g}_{\uparrow\downarrow}^{\text{sta}}(\mathbf{q}) = - \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{1}{4} \left[ 1 - \frac{\xi_+}{\epsilon_+} (1 - 2f_+) \right] \left[ 1 - \frac{\xi_-}{\epsilon_-} (1 - 2f_-) \right] \quad (61) $$

Les formes factorisées de l'intégrande – en une fonction de $\epsilon_+$ fois la même fonction de $\epsilon_-$ – annoncent les expressions de Wick (23,24) des $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$, qui sont des carrés parfaits, et que l'on retrouve par le changement de variables de jacobien unité $(\mathbf{k}, \mathbf{q}) \mapsto (\mathbf{k}_+, \mathbf{k}_-)$ dans la prise de la transformée de Fourier inverse des $\delta \tilde{g}_{\sigma\sigma'}(\mathbf{q})$.

## 2.4. Différentes formulations en termes d'intégrales curvilignes et applications

L'écriture de la contribution collective $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ dans (55), reflétant directement celle originelle du théorème de fluctuation-dissipation (39), est malcommode à utiliser, analytiquement ou numériquement, car son intégrande est très singulier.

En effet, comme on le voit sur leurs définitions (43–49), les fonctions $z \mapsto \Sigma_{ij}(\mathbf{q}, z)$ admettent des lignes de coupure sur l'axe réel, là où les dénominateurs s'annulent pour certaines valeurs du vecteur d'intégration $\mathbf{k}$, en pratique de son module $k$ et de son angle $\theta$ avec le vecteur d'onde $\mathbf{q}$ (tel que $\mathbf{k} \cdot \mathbf{q} = kq \cos\theta$). Le dénominateur des secondes contributions entre accolades s'annule si $z = \pm(\epsilon_+ - \epsilon_-)$; c'est la condition de résonance du processus de transfert d'une quasi-particule fermionique préexistante – thermique – du vecteur d'onde $\mathbf{k}_- = \mathbf{k} - \mathbf{q}/2$ au vecteur d'onde $\mathbf{k}_+ = \mathbf{k} + \mathbf{q}/2$ ou le contraire par un potentiel extérieur d'excitation de vecteur d'onde $\mathbf{q}$ et de pulsation $z = \hbar\omega$ (une excitation de Bragg dans les expériences d'atomes froids [39]). Or, lorsque $\mathbf{k}$ décrit $\mathbb{R}^d$, la différence d'énergie $\epsilon_+ - \epsilon_-$ décrit $\mathbb{R}^d$;<fnref n="17" /> les fonctions $\Sigma_{ij}$ admettent donc l'axe réel tout entier comme ligne de coupure.

Heureusement, la situation se simplifie un peu à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ où nous nous placerons – sauf exception – dans la section 3 : les nombres d'occupation fermioniques $f_\pm$ sont alors mis à zéro, et les fonctions $\Sigma_{ij}$ se réduisent aux premières contributions entre accolades.

<footnote n="16">¹⁶En gros, il suffit d'utiliser la propriété $\text{Im}\{\epsilon_0 / [(\epsilon + i0^+) - \epsilon_0^2]\} = (\pi/2)[\delta(\epsilon + \epsilon_0) - \delta(\epsilon - \epsilon_0)]$ et de ne pas passer à côté de simplifications remarquables comme $(1 - f_+ - f_-)[\exp[\beta(\epsilon_+ + \epsilon_-)] + 1]/[\exp[\beta(\epsilon_+ + \epsilon_-)] - 1] - (f_+ - f_-)[\exp[\beta(\epsilon_+ - \epsilon_-)] + 1]/[\exp[\beta(\epsilon_+ - \epsilon_-)] - 1] = 1$. Pour faire disparaître le Dirac dans (40), on se souvient de l'expression (21) de $\rho^{\text{sta}}$ et on réalise que $\int d^d q \exp[i\mathbf{q} \cdot (\mathbf{r} - \mathbf{r}')] \int d^d k f(\mathbf{k}_+) = (2\pi)^d \delta(\mathbf{r} - \mathbf{r}') \int d^d k f(\mathbf{k})$ pour toute fonction $f$ d'intégrale convergente sur $\mathbb{R}^d$.</footnote>

<footnote n="17">¹⁷Pour le voir, il suffit, à $\mathbf{q}$ non nul fixé, de prendre les limites $k \to +\infty$ et $\cos\theta \to 0$ en maintenant constant le produit $k \cos\theta$; alors $\epsilon_+ - \epsilon_- \sim \xi_+ - \xi_- = \hbar qk(\cos\theta)/m$ tend vers toute valeur souhaitée.</footnote>
Dans ce cas, les fonctions $\Sigma_{ij}$ présentent une ligne de coupure seulement sur les demi-droites réelles $]-\infty, -E_{\text{bord}}(\mathbf{q})]$ et $[E_{\text{bord}}(\mathbf{q}), +\infty[$, où $E_{\text{bord}}(\mathbf{q}) = \text{inf}_k \epsilon_+ + \epsilon_-$ est le bord inférieur du continuum de paire brisée à vecteur d'onde total $\mathbf{q}$ fixé (sous l'effet d'un potentiel d'excitation transmettant une impulsion $\hbar\mathbf{q}$ bien définie, une paire liée $\uparrow\downarrow$ du condensat peut se briser en deux quasi-particules fermioniques de quantité de mouvement totale fixée $\hbar\mathbf{q}$ mais de vecteur d'onde relatif $\mathbf{k}$ libre dans $\mathbb{R}^d$; aussi la variation d'énergie correspondante décrit-elle un intervalle de $\mathbb{R}$). On a explicitement [23, 40] :<fnref n="18" />

$$E_{\text{bord}}(\mathbf{q}) = \begin{cases} 2\Delta & \text{si } \mu \ge 0 \text{ et } q < 2(2m\mu/\hbar^2)^{1/2} \\ 2\epsilon_{\mathbf{q}/2} & \text{sinon} \end{cases} \quad (62)$$

où, dans le cas $\mu > 0$, $k_{\min} = (2m\mu/\hbar^2)^{1/2}$ est la position du minimum $E_{\text{gap}} = \Delta$ de la relation de dispersion $k \mapsto \epsilon_k$ des quasi-particules fermioniques (pour $\mu < 0$, le minimum $E_{\text{gap}} = (\mu^2 + \Delta^2)^{1/2}$ est atteint en $k = 0$). De plus, si le nombre d'onde $q$ appartient au domaine d'existence $q \in \mathcal{E}_{\text{son}}$ de la branche d'excitation acoustique dans la théorie BCS dépendant du temps linéarisée,<fnref n="19" /> la susceptibilité collective (55) admet un pôle en l'énergie propre $\hbar\omega_q \in ]0, E_{\text{bord}}(\mathbf{q})[$ du phonon de Bogoliubov-Anderson correspondant, et en son opposé $-\hbar\omega_q$, puisque l'on a par définition<fnref n="20" />,<fnref n="21" />

$$[\Sigma_{11}(\mathbf{q}, z)\Sigma_{22}(\mathbf{q}, z) - \Sigma_{12}^2(\mathbf{q}, z)] \Big|_{z=\pm\hbar\omega_q} = 0 \quad (63)$$

Tirons les fruits de cette discussion, et remplaçons dans l'expression (55) de $\delta\tilde{\chi}_{\text{coll}}(\mathbf{q})$ l'axe réel par d'utiles chemins d'intégration restant à distance non nulle des singularités dans le plan complexe.

Le premier choix, applicable à tous les ordres en $\exp(-\beta E_{\text{gap}})$ et utilisé dans nos calculs numériques de la section 3.2, tire parti de l'analyticité de $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ sur $\mathbb{C} \setminus \mathbb{R}$ et de la propriété $[\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)]^* = \tilde{\chi}_{\text{coll}}(\mathbf{q}, z^*)$ (valables en réalité au-delà de l'approximation BCS dépendant du temps linéarisée) pour écrire

$$\mathrm{Im}[\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+)] = [\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+) - \tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon - i0^+)] / (2i) \quad (64)$$

et pour séparer les deux domaines d'intégration apparus $\mathbb{R} + i0^+$ et $\mathbb{R} - i0^+$ en deux chemins $-\infty + i\eta \to +\infty + i\eta$ et $+\infty - i\eta \to +\infty - i\eta$, parcourus en sens inverse et à distance non infinitésimale

<footnote n="18">¹⁸Il serait trop naïf de croire que les lignes de coupure admettent comme seuls points de branchement leurs extrémités $\pm E_{\text{bord}}(\mathbf{q})$. Suivant les valeurs de $q$ et de la force des interactions, il peut y avoir un ou deux autres points de branchement dans $]-E_{\text{bord}}(\mathbf{q}), +\infty[$ et dans $]-\infty, -E_{\text{bord}}(\mathbf{q})[$, que l'on identifie comme points anguleux dans des densités spectrales $\rho_{ij}(\mathbf{q}, \epsilon)$ (telles que $\Sigma_{ij}(\mathbf{q}, z) \propto \int d\epsilon \rho_{ij}(\mathbf{q}, \epsilon)/(z-\epsilon)$ ou, si l'on préfère, $\text{Im}\,\Sigma_{ij}(\mathbf{q}, \epsilon + i0^+) \propto \rho_{ij}(\mathbf{q}, \epsilon)$) [20]. Le même phénomène se produit au-delà de l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, avec la possibilité de points de branchement supplémentaires, purement thermiques [41]. Les références [20, 41] sont à 3D mais leurs considérations s'appliquent à 2D, sauf bien sûr l'expression finale des densités spectrales en termes d'intégrales elliptiques dans [20], qui résulte de la forme $\int d(\cos\theta)$ de l'intégrale polaire en dimension trois.</footnote>

<footnote n="19">¹⁹À 2D, un mode de phonon existe à tout nombre d'onde $q > 0$, $\mathcal{E}_{\text{son}}^{d=2} = ]0, +\infty[$ [24]. À 3D, le domaine d'existence est borné connexe, $\mathcal{E}_{\text{son}}^{d=3} = ]0, q_{\text{sup}}[$, pour $a_{3D} < 0$ c'est-à-dire pour $\Delta/\mu \in ]0, 1, 162[$; il est à deux composantes connexes, $\mathcal{E}_{\text{son}}^{d=3} = ]0, q_{\text{sup}}[\cup \{q_{\text{inf}}, +\infty[ \text{ si } \Delta/\mu \in ]1, 162; 1, 729[$; il vaut enfin le demi-axe réel tout entier, $\mathcal{E}_{\text{son}}^{d=3} = ]0, +\infty[$, si $\Delta/\mu > 1, 729$ ou $\mu < 0$ [23, 42].</footnote>

<footnote n="20">²⁰On obtient la même équation aux valeurs propres dans l'approximation des fluctuations gaussiennes, c'est-à-dire dans une formulation par intégrale de chemin où, après introduction d'un champ auxiliaire dans la transformation de Hubbard-Stratonovich puis sommation sur les degrés de liberté fermioniques, on quadratise l'action obtenue autour de son point stationnaire BCS statique [19, 25, 28, 32, 43] (on tombe ainsi sur une intégrale de chemin gaussienne qu'on sait calculer mais c'est une autre histoire). Pour convaincre le lecteur, comparons explicitement à [28] : dans les notations de cette référence, l'équation aux valeurs propres s'écrit $M_{11}(\mathbf{q}, z) M_{22}(\mathbf{q}, z) = M_{12}(\mathbf{q}, z) M_{21}(\mathbf{q}, z)$; or, nous trouvons que $M_{11} = (\Sigma_{11} + \Sigma_{22})/2 + \Sigma_{12}$, $M_{22} = (\Sigma_{11} + \Sigma_{22})/2 - \Sigma_{12}$ et $M_{12} = M_{21} = (\Sigma_{22} - \Sigma_{11})/2$ en omettant à nouveau les variables $(\mathbf{q}, z)$. Il suffit de développer et simplifier pour retrouver (63).</footnote>

<footnote n="21">²¹Si l'on quitte l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, $z \mapsto \tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ n'admet plus de pôle; en revanche, son prolongement analytique à travers sa ligne de coupure du demi-plan supérieur au demi-plan inférieur en admet, et les énergies propres du son correspondantes sont complexes – les phonons deviennent des résonances [41].</footnote>
$\eta > 0$ de l'axe réel, appelés respectivement $C_h$ (chemin haut) et $C_b$ (chemin bas) sur la figure 2a. On aboutit alors à

$$ \begin{aligned} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) &= -k_B T \tilde{\chi}_{\text{coll}}(\mathbf{q}, 0) + \int_{C_h \cup C_b} \frac{dz}{(-2i\pi)} \frac{\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)}{1 - \exp(-\beta z)} \\ &= -k_B T \tilde{\chi}_{\text{coll}}(\mathbf{q}, 0) + \text{Im} \int_{-\infty}^{+\infty} \frac{d\epsilon}{(-\pi)} \frac{\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i\eta)}{1 - \exp[-\beta(\epsilon + i\eta)]} \end{aligned} \quad (65) $$

Deux remarques s'imposent : (i) la séparation en deux chemins a fait apparaître un pôle non physique en $z = 0$, dont le contre-terme $\propto T$ au second membre de (65) compense la contribution déduite du théorème des résidus;<fnref n="22" /> ce pôle ne posait aucun problème dans la formulation originelle (55) car $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z = 0)$ est de valeur réelle, donc de partie imaginaire nulle; (ii) le dénominateur $1 - \exp(-\beta z)$ fait apparaître les pôles thermiques $z = 2i\pi nk_B T$, $n \in \mathbb{Z}^*$; ceux d'indices $n = \pm 1$ limitent l'excursion verticale de $C_h$ et $C_b$; dans le numérique, pour éviter de trop s'en approcher, nous prenons $\eta = \pi k_B T$ selon la doctrine du « juste milieu », et calculons en pratique seulement la contribution de $C_h$, celle de $C_b$ s'en déduisant par changement de signe et conjugaison complexe, comme au troisième membre de (65).<fnref n="23" />

À l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, un deuxième choix de chemin d'intégration s'impose. Puisqu'il n'y a plus de ligne de coupure entre $-E_{\text{bord}}(\mathbf{q})$ et $E_{\text{bord}}(\mathbf{q})$, on peut enlever dans $C_h \cup C_b$ une portion centrale entourant le pôle non physique $z = 0$ comme seule singularité, au prix de l'ajout de segments verticaux AD et BC reliant $C_h$ à $C_b$, voir la figure 2a; on écrit donc $C_h \cup C_b = C_g \cup ABCD \cup C_d$ où le rectangle ABCD est parcouru dans le sens des aiguilles d'une montre, le contour $C_g$ est à gauche et le contour $C_d$ est à droite, voir la figure 2b. Les segments verticaux ajoutés ont une contribution nulle dans l'ensemble, car ils sont parcourus en des sens opposés;

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) = -k_B T \tilde{\chi}_{\text{coll}}(\mathbf{q}, 0) - 2k_B T \text{Re} \sum_{n=1}^{+\infty} \tilde{\chi}_{\text{coll}}(\mathbf{q}, 2i\pi nk_B T) = -k_B T \sum_{n=-\infty}^{+\infty} \tilde{\chi}_{\text{coll}}(\mathbf{q}, 2i\pi nk_B T) \quad (66) $$

<footnote n="22">²²On retrouve ce résultat avec une astuce moins-plus et la théorie des distributions : pour que l'intégrande dans (55) devienne une fonction seulement de $\epsilon + i0^+$ (et pas de $\epsilon$), on écrit $\text{Im}[\tilde{\chi}(\mathbf{q}, \epsilon + i0^+)]/[1 - \exp(-\beta\epsilon)] = \text{Im}[\tilde{\chi}(\mathbf{q}, \epsilon + i0^+)]/[1 - \exp(-\beta\epsilon)] = \text{Im}[\tilde{\chi}(\mathbf{q}, \epsilon + i0^+)] + \text{Im}[\tilde{\chi}(\mathbf{q}, \epsilon + i0^+)]\phi(\epsilon)]$ en omettant l'indice « coll » pour alléger. La distribution introduite, $\phi(\epsilon) = 1/[1 - \exp(-\beta\epsilon)] - 1/[1 - \exp(-\beta(\epsilon + i0^+))]$, n'a une action non nulle que sur un voisinage de l'origine; on peut donc linéariser l'exponentielle et écrire $\phi(\epsilon) = k_B[1/\epsilon - 1/(\epsilon + i0^+)] = k_B T[1/\epsilon - \text{v.p.}(1/\epsilon) + i\pi\delta(\epsilon)]$ avec v.p. la distribution valeur principale et $\delta$ la distribution de Dirac. Comme $\text{Im}[\tilde{\chi}(\mathbf{q}, \epsilon + i0^+)]$ tend vers zéro en $\epsilon = 0$, la distribution $1/\epsilon - \text{v.p.}(1/\epsilon)$ est sans effet, seule la contribution du Dirac est non nulle, et l'on reproduit le contre-terme de (65).</footnote>

<footnote n="23">²³ Dans ce troisième membre, refermons le chemin d'intégration par un demi-cercle à l'infini dans le demi-plan complexe supérieur (en admettant que l'intégrale soit nulle sur le demi-cercle) puis appliquons le théorème des résidus : comme la susceptibilité est analytique pour $\text{Im } z > 0$, les seules singularités entourées sont les pôles thermiques et il vient</footnote>

<footnote>après usage de la propriété $[\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)]^* = \tilde{\chi}_{\text{coll}}(\mathbf{q}, z^*)$ dans la seconde égalité. Nous aboutissons ainsi à une série de Matsubara bosonique [44] (les $2i\pi nk_B T$ sont précisément les énergies complexes de Matsubara pour des bosons) sans jamais avoir introduit de fonction de Green à $N$ corps ni de conditions aux limites périodiques en temps imaginaire! Pour notre calcul numérique de la section 3.2, qui met une coupure $k < k_{\text{max}}$ sur le vecteur d'onde interne des paires dans les fonctions $\Sigma_{ij}$, la représentation (66) aurait, par rapport à l'intégration sur $\epsilon$ dans (65), l'avantage d'un pas $\delta\epsilon$ plus grand mais l'inconvénient d'une convergence asymptotique beaucoup plus lente : l'erreur commise en tronquant la série de Matsubara à $n = n_{\text{max}}$ tend vers zéro comme $1/n_{\text{max}}$, au lieu de s'annuler exactement au-delà d'un certain $\epsilon_{\text{max}}$ dans (65), voir la note 38; il y aurait bien une parade (on pourrait estimer analytiquement le bout négligé $n > n_{\text{max}}$ à l'aide d'un développement asymptotique du sommande) mais elle serait laborieuse; aussi n'avons nous pas utilisé (66). Pour terminer, remarquons que, dans la limite d'une température nulle $T \to 0^+$, la série de Matsubara – vue comme une somme de Riemann – est remplacée par une intégrale, celle de l'expression (70) avec $\eta = +\infty$.</footnote>
<figure><img src="image_2.png" /></figure>

FIGURE 2 – Différents chemins d’intégration utilisés à température $T > 0$ dans la réécriture de la quantité $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ de l’équation (55) sous forme d’intégrales curvilignes dans le plan complexe, ayant l’avantage – aussi bien numérique qu’analytique – de rester à distance non infinitésimale des singularités de l’intégrande. (a) Dans le cas général, à tous les ordres en $\exp(-\beta E_{\text{gap}})$, voir l’équation (65). (b) À l’ordre zéro en $\exp(-\beta E_{\text{gap}})$ dans la susceptibilité $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ mais pas dans le dénominateur thermique de l’intégrande de $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$, voir l’équation (67). (c) À l’ordre zéro en $\exp(-\beta E_{\text{gap}})$ partout, voir l’équation (68). Lignes hachurées : lignes de coupure de la susceptibilité. Symboles $\times$: pôles $\pm\hbar\omega_q$ de la susceptibilité à l’ordre zéro en $\exp(-\beta E_{\text{gap}})$ ($\omega_q > 0$ est la pulsation propre du phonon de Bogolioubov-Anderson, solution de l’équation (63), si elle existe au nombre d’onde $q$). Étoiles : pôles de Bose $z = 2i\pi nk_B T$, $n$ entier relatif (racines du dénominateur thermique) dans (55); le pôle non physique $n = 0$ est en rouge. Le bord $E_{\text{bord}}(\mathbf{q})$ du continuum de paire brisée à vecteur d’onde $\mathbf{q}$ fixé est donné par l’équation (62). Dans les noms des chemins, les indices h et b signifient haut et bas, les indices g et d signifient gauche et droite, C signifie chemin, c signifie contour, L signifie lacet. Les flèches marquent l’orientation des chemins. Le rectangle ABCD intervient dans le raisonnement faisant passer de la forme (65) à la forme (67). En (a), les chemins $C_h$ et $C_b$ doivent passer entre l’axe réel et les pôles de Bose $n = \pm 1$ qui en sont les plus proches. Cette contrainte est bien sûr inopérante en (b) et en (c).

l’intégrale curviligne sur ABCD compense le contre-terme d’après le théorème des résidus; il reste

$$ \boxed{\delta \tilde{g}_{\text{coll}}(\mathbf{q}) = \int_{C_g \cup C_d} \frac{\mathrm{d}z}{(-2i\pi)} \frac{\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)}{1 - \exp(-\beta z)} = \int_{C_d} \frac{\mathrm{d}z}{(-i\pi)} \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \left( \frac{1}{2} + \frac{1}{\exp(\beta z) - 1} \right)} \quad (67) $$

On a simplifié encore l’expression au troisième membre, en utilisant le fait que $C_g$ se déduit de $C_d$ par changement de $z$ en $-z$; il en resulte que $\int_{C_g} dz f(z) = \int_{C_d} dz [-f(-z)]$ pour toute fonction $f$ de la variable complexe, et donc que seule la partie impaire de l’intégrande subsiste, ce qui fait apparaître une distribution de Bose-Einstein.<fnref n="24" /> Cependant, l’ordre zéro en $\exp(-\beta E_{\text{gap}})$ n’est

<footnote n="24">²⁴La susceptibilité $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ est une fonction paire de $z$ – voir les remarques après (50,51) – et l’on a $1/[1 - \exp(-\beta z)] - 1/[1 - \exp(\beta z)] = 1 + 2/[\exp(\beta z) - 1]$.</footnote>
pas implémenté de façon entièrement cohérente, puisqu'on a omis les nombres d'occupation fermioniques $f_{\pm}$ dans $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ tout en gardant la loi de Bose thermique pour des énergies complexes $z$ de partie réelle arbitrairement grande dans $C_d$.

Nous effectuons donc un dernier changement de chemin d'intégration, en tirant parti de l'absence de singularité entre le pôle de phonon $z = \hbar\omega_q$ et le point de branchement $z = E_{\text{bord}}(\mathbf{q})$ du continuum, pour remplacer $C_d$ par $c_d \cup L_d$ où le contour circulaire $c_d$ entoure le pôle de phonon et le lacet $L_d$ entoure la ligne de coupure – le continuum de paire brisée à $\mathbf{q}$ fixé – dans le sens des aiguilles d'une montre comme sur la figure 2c. Comme les points du lacet sont de partie réelle $\text{Re} z \gtrsim E_{\text{bord}}(\mathbf{q}) \ge E_{\text{gap}}$, il faut y négliger la loi de Bose thermique à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ où nous sommes. La même approximation ne peut être faite de manière générale sur $c_d$ car l'énergie propre du phonon $\hbar\omega_q$ devient arbitrairement inférieure à $k_B T > 0$ à faible nombre d'onde $q$. Nous séparons en définitive la transformée de Fourier de $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ en une contribution du phonon $\phi$ de Bogolioubov-Anderson (s'il existe au nombre d'onde $\mathbf{q}$ considéré, voir la note 19 sur le domaine d'existence $\mathcal{E}_{\text{son}}$ de la branche acoustique, sinon $\delta \tilde{g}_{\text{coll}}^{\phi}(\mathbf{q}) = 0$) et une contribution collective du continuum de paire brisée,

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) = \delta \tilde{g}_{\text{coll}}^{\phi}(\mathbf{q}) + \delta \tilde{g}_{\text{coll}}^{C^0}(\mathbf{q}) \quad \text{où} \quad \left\{ \begin{array}{l} \delta \tilde{g}_{\text{coll}}^{\phi}(\mathbf{q}) = \int_{c_d} \frac{dz}{(-i\pi)} \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \left( \frac{1}{2} + \frac{1}{\exp(\beta z) - 1} \right) = \left( \tilde{n}_{\mathbf{q}} + \frac{1}{2} \right) \Phi(\mathbf{q}) \\ \\ \delta \tilde{g}_{\text{coll}}^{C^0}(\mathbf{q}) = \int_{L_d} \frac{dz}{(-2i\pi)} \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \end{array} \right. \qquad (68) $$

Ici, en vertu du théorème des résidus appliqué à l'intégrale sur le contour $c_d$, $\tilde{n}_q$ est le nombre d'occupation thermique du mode de phonon $\mathbf{q}$ et $\Phi(\mathbf{q})$ est deux fois le résidu de la susceptibilité collective en $z = \hbar\omega_q$ :

$$ \tilde{n}_{\mathbf{q}} = \frac{1}{\exp(\beta\hbar\omega_{\mathbf{q}}) - 1} \quad \text{et} \quad \Phi(\mathbf{q}) = \frac{-s_{11} (|\Sigma_{11}|^{1/2} \Sigma_{23} - s_{11} s_{12} |\Sigma_{22}|^{1/2} \Sigma_{13})^2}{\frac{d}{dz} (\Sigma_{11} \Sigma_{22} - \Sigma_{12}^2)} \Big|_{(\mathbf{q}, z=\hbar\omega_{\mathbf{q}})} \quad (69) $$

avec $s_{ij}$ le signe de $\Sigma_{ij}(\mathbf{q}, \hbar\omega_q)$. L'expression de $\Phi(\mathbf{q})$ se vérifie en développant le carré au numérateur et en reconnaissant la définition (55) au moyen de l'équation (63), selon laquelle $s_{11} = s_{22}$ et $\Sigma_{12} = s_{12}|\Sigma_{11}\Sigma_{22}|^{1/2}$ au point considéré.

Dans les formulations (67,68), les distances à l'axe réel sur les chemins $C_d$, $c_d$ et $L_d$ ne sont plus contraintes par la température (elles n'ont pas à être inférieures à $2\pi k_B T$), mais le rayon de $c_d$ doit rester inférieur à $\hbar\omega_q$ et à $E_{\text{bord}}(\mathbf{q}) - \hbar\omega_q$. Si l'on veut que le chemin d'intégration reste à distance non infinitésimale des singularités de l'intégrande, (i) le chemin $c_d \cup L_d$ ne convient pas près d'un point $q_*$ où la branche acoustique rejoint le continuum de paire brisée $(E_{\text{bord}}(\mathbf{q}) - \hbar\omega_q \to 0)$ si $q \to q_*$ et il vaut mieux revenir au chemin $C_d$ et à l'expression (67), (ii) même $C_d$ ne convient plus dans la limite $q \to 0$ où $\hbar\omega_q \to 0$ et il faut revenir à $C_h$ et à l'expression (65). Pour conclure, expliquons la particularité du cas $T = 0$ qui nous a conduit à l'éviter dans nos calculs numériques : la définition (55) de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ se réduit certes à

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) \stackrel{T=0}{=} \int_{0}^{+\infty} \frac{d\epsilon}{(-\pi)} \text{Im}[\tilde{\chi}_{\text{coll}}(\mathbf{q}, \epsilon + i0^+)] = \text{Im} \left[ \left( \int_{0}^{i\eta} + \int_{i\eta}^{i\eta+\infty} \right) \frac{dz}{(-\pi)} \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \right], \quad (70) $$

il n'y a plus trace de pôle non physique en $z = 0$ et l'écart vertical $\eta > 0$ à la ligne de coupure peut être choisi librement au troisième membre, mais le chemin d'intégration, partant forcément de l'origine, ne peut rester à distance non infinitésimale du pôle de phonon $z = \hbar\omega_q$ si $q \to 0$.

## 2.5. Quelle équation d'état dans notre théorie BCS dépendant du temps linéarisée?

Par essence, le théorème de fluctuation-dissipation donne accès aux écarts $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ des fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ à leurs valeurs asymptotiques $\rho_{\sigma}\rho_{\sigma'}$. Or, dans les gaz
d'atomes froids, ce sont précisément ces fonctions $g_{\sigma\sigma'}$ qui intéressent les numériciens et les expérimentateurs, voir notre section 3.2. Il nous faut donc connaître la densité du gaz, ici non polarisé ($\rho_\uparrow = \rho_\downarrow = \rho/2$), pour compléter les $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ en les $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$.

Grâce à la théorie BCS dépendant du temps linéarisée, nous avons amélioré le calcul des $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ au moyen des fonctions de réponse densité-densité $\chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t')$ mais nous n'avons pas touché au reste de la théorie BCS statique, en particulier pas aux observables à un corps. Il semblerait donc entendu que l'équation d'état grand-canonique $\rho(\mu, T)$ reste la même, et que l'on ait

$$g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{?}{=} \rho_{\sigma}^{\text{sta}} \rho_{\sigma'}^{\text{sta}} + \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \quad (71)$$

où $\rho_{\sigma}^{\text{sta}}$ est donnée par (21). En réalité, il en va autrement, car la fonction $\delta g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}')$ contient de façon un peu cachée une information sur la densité moyenne. En effet, les particules dans l'état interne $\uparrow$ sont des fermions identiques sans interaction directe entre eux (pas de résonance dans l'onde $p$), si bien que leur densité de paires doit être de limite nulle à distance nulle :

$$\lim_{|\mathbf{r}-\mathbf{r}'|\to 0^+} g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}') = 0 \quad \text{donc} \quad \lim_{|\mathbf{r}-\mathbf{r}'|\to 0^+} \delta g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}') = -\rho_{\uparrow} \rho_{\downarrow} \quad (72)$$

Nous en déduisons la bonne valeur de la densité du gaz non polarisé dans notre théorie BCS dépendant du temps linéarisée et les bonnes expressions des distributions de paires :<fnref n="25" />

$$\frac{1}{4}[\rho^{\text{dtl}}(\mu, T)]^2 = \frac{1}{4}[\rho^{\text{sta}}(\mu, T)]^2 - \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) \quad \text{et} \quad g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') = \rho_{\sigma}^{\text{dtl}} \rho_{\sigma'}^{\text{dtl}} + \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \quad (73)$$

Dans les exemples numériques de la section 3.2, à très basse température $k_B T \ll mc^2$, $E_{\text{gap}}$ donc pratiquement dans l'état fondamental, nous trouvons que la fonction $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$, définie dans l'équation (55), est partout négative, voir l'équation (97); on a alors $\rho^{\text{dtl}} > \rho^{\text{sta}}$ à potentiel chimique $\mu$ fixé ou, si l'on préfère, $\mu^{\text{dtl}} < \mu^{\text{sta}}$ à densité $\rho$ fixée (la stabilité thermodynamique du gaz en phase spatialement homogène exige que $\partial_\mu \rho > 0$, voir l'équation (32)). Voilà qui va dans le bon sens, celui d'une amélioration de la théorie BCS variationnelle.

Grâce à l'équation d'état (73), nous pouvons maintenant préciser la valeur commune des deux différences dans la belle relation (59) :

$$g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) [e^{i\mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')} - 1] \quad (74)$$

Rappelons-le, les deux théories (statique ou pas) doivent être prises dans (74) au même potentiel chimique, pas à la même densité; c'est évident dans la présente section 2 puisque nous sommes dans un point de vue grand-canonique, mais ça ne le sera pas toujours dans la suite, voir la section 3.2. Nous sommes parvenu au résultat (74) – qui servira en particulier dans la section 3.3 – au moyen de la chaîne d'égalités évidentes au vu des équations (1,56,57,73) :

$$g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \rho_{\sigma}^{\text{dtl}} \rho_{\sigma'}^{\text{dtl}} - \rho_{\sigma}^{\text{sta}} \rho_{\sigma'}^{\text{sta}} + \left[ \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - \delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}') \right] = \frac{1}{4} [\rho^{\text{dtl}}(\mu, T)]^2 \\ - \frac{1}{4} [\rho^{\text{sta}}(\mu, T)]^2 + \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') = - \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) + \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) e^{i\mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')} \quad (75)$$

<footnote n="25">²⁵Nous avons utilisé la belle relation (58), la propriété (72) pour BCS statique et BCS dépendant du temps, puis la séparation (56) de $\delta g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}')$ en partie statique et collective; l'intégrale sur $\mathbf{q}$ dans (73) n'est autre que la limite de $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ lorsque $|\mathbf{r}-\mathbf{r}'| \to 0^+$, limite dont l'existence et l'expression annoncée sont établies dans la section 3.3 pour la théorie à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$.</footnote>
### 3. Résultats explicites sur les $g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$, numériques et analytiques

Dans cette section, nous présentons l'ensemble de nos résultats sur les fonctions de distribution de paires dans la théorie BCS dépendant du temps linéarisé. D'abord à grande distance dans la section 3.1, où nous obtenons analytiquement à basse température une belle loi raccordant le régime asymptotique quantique $|\mathbf{r}-\mathbf{r}'| \to +\infty$ mais $|\mathbf{r}-\mathbf{r}'| \ll \lambda_\phi$ au régime asymptotique thermique $|\mathbf{r}-\mathbf{r}'| \gg \lambda_\phi \to +\infty$ où $\lambda_\phi = \hbar c / 2\pi k_B T$ est la longueur d'onde thermique des phonons du superfluide. Ensuite numériquement aux distances intermédiaires dans la section 3.2, d'une fraction de $k_F^{-1}$ à quelques $k_F^{-1}$ dans le régime d'interaction forte, domaine de prédilection des simulations de Monte-Carlo quantique et des mesures expérimentales dans les gaz d'atomes froids. Enfin à courte distance donc aux grands nombre d'onde $q$ dans la section 3.3, où des lois simples décrivent la divergence de $g_{11}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ et le départ de $g_{11}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ de la valeur nulle. Tous ces résultats sont obtenus à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, c'est-à-dire à suffisamment basse température pour qu'on puisse mettre omettre les nombres d'occupation $f_k$ des quasi-particules fermioniques dans la théorie BCS statique et dans les fonctions $\Sigma_{ij}$ donnant les susceptibilités densité-densité $\tilde{\chi}_{\sigma\sigma'}(\mathbf{q}, z)$.

#### 3.1. Comportement des $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ à grande distance et basse température

##### 3.1.1. Ce qu'on déduit de l'intégration par parties

Une première technique pour déterminer le comportement à grande distance d'une intégrale de Fourier isotrope comme celle (57) donnant $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ consiste, après moyenne angulaire $\langle\exp[i\mathbf{q}\cdot(\mathbf{r}-\mathbf{r}')]\rangle_{\text{ang}} = \sin(q|\mathbf{r}-\mathbf{r}'|)/(q|\mathbf{r}-\mathbf{r}'|)$ en dimension $d=3$ ou $\langle\exp[i\mathbf{q}\cdot(\mathbf{r}-\mathbf{r}')]\rangle_{\text{ang}} = J_0(q|\mathbf{r}-\mathbf{r}'|)$ en dimension $d=2$, à intégrer par parties de façon répétée [45]. Nous reprenons ici les équations (8.300,8.307) de la référence [5] :

$$ \int_0^{+\infty} dq \phi(q) \sin(q|\mathbf{r}-\mathbf{r}'|) \Big|_{|\mathbf{r}-\mathbf{r}'| \to +\infty} = \sum_{s=0}^{S} (-1)^s \frac{\phi^{(2s)}(0)}{|\mathbf{r}-\mathbf{r}'|^{2s+1}} + O\left(\frac{1}{|\mathbf{r}-\mathbf{r}'|^{2S+3}}\right) \quad (76) $$

$$ \int_0^{+\infty} dq J_0(q|\mathbf{r}-\mathbf{r}'|)\phi(q) \Big|_{|\mathbf{r}-\mathbf{r}'| \to +\infty} = \sum_{s=0}^{S} (-1)^s \frac{\phi^{(2s)}(0)\Gamma(s+1/2)}{\Gamma(1/2)s!|\mathbf{r}-\mathbf{r}'|^{2s+1}} + O\left(\frac{1}{|\mathbf{r}-\mathbf{r}'|^{2S+3}}\right) \quad (77) $$

où $S$ est un entier naturel quelconque, $\phi(q) = q\delta\tilde{g}_{\text{coll}}(\mathbf{q})/2\pi^2|\mathbf{r}-\mathbf{r}'|$ en dimension $d=3$ et $\phi(q) = q\delta\tilde{g}_{\text{coll}}(\mathbf{q})/2\pi$ en dimension $d=2$.

Or, pour $\text{Im}z \neq 0$, les fonctions $\Sigma_{ij}(\mathbf{q}, z)$ de la section 2.2 sont formellement des fonctions paires du nombre d'onde $q$, au sens où elles sont des intégrales sur un vecteur d'onde interne $\mathbf{k}$ de fonctions lisses de $\mathbf{k}$ invariantes par changement de $\mathbf{q}$ en $-\mathbf{q}$; cette propriété s'étend à la susceptibilité $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ donnant $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ dans (55) et, après intégration sur un chemin restant toujours à distance non infinitésimale des singularités de l'intégrande, à la fonction $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ elle-même. À température non nulle, cette contrainte est satisfaite par la représentation (65), si bien que $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ n'admet que des puissances paires dans son développement limité en $q=0$ :

$$ \delta\tilde{g}_{\text{coll}}(\mathbf{q}) \stackrel{T>0}{=} \sum_{q \to 0}^{S} c_{2S} q^{2S} + O(q^{2S+2}) \quad (78) $$

avec des coefficients qui importent peu pour l'instant. À l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, le même raisonnement sur le lacet $L_d$ dans (68) montre que la contribution du continuum – indépendante de la température – jouit de la même propriété,

$$ \delta\tilde{g}_{\text{coll}}^{C_0}(\mathbf{q}) = \sum_{q \to 0}^{S} c_{2S}^{C_0} q^{2S} + O(q^{2S+2}) \quad (79) $$
avec les deux premiers coefficients calculés analytiquement dans l'annexe A. Bien que l'équation aux valeurs propres (63) soit formellement une fonction paire de $q$ et de $\omega_q$, la pulsation propre $\omega_q$ du phonon $q$ dans la théorie BCS dépendant du temps linéarisée est de départ linéaire comme dans (10) donc n'admet que des puissances impaires,<fnref n="26" />

$$ \hbar\omega_q = \hbar c q \left[ 1 + \sum_{s=1}^{S} \zeta_s \left( \frac{\hbar q}{mc} \right)^{2s} + O(q^{2S+2}) \right], \quad (80) $$

dont les premiers coefficients $\zeta_s$ sont connus [24, 42].<fnref n="27" />,<fnref n="28" /> Écrivons tout de même le coefficient $\zeta_2$ à 3D, car l'astuce d'écriture (33) nous permet de le simplifier grandement :<fnref n="29" />

$$ \zeta_2^{\text{dtdl}} = \frac{32 + 64z^2 + 52z^4 + 4xz(24 + 53z^2 + 39z^4) + x^2(-35 - 56z^2 + 13z^4 + 54z^6)}{-8 \times 135x^2(1 + z^2)^3} \bigg|_{x=\frac{\Delta}{\mu}, z=\frac{2\Delta}{3\rho^{\text{sta}}g_{3D}}} \quad (81) $$

De façon similaire, comme son expression (69) est formellement paire en $q$ mais impaire en $\omega_q$ (à cause de la dérivée $d/dz$ au dénominateur), la fonction $\Phi(q)$ dans la contribution à $\delta\tilde{g}_{\text{coll}}(q)$ du phonon de Bogolioubov-Anderson n'admet que des puissances impaires :<fnref n="30" />

$$ \Phi(\mathbf{q}) = \sum_{q \to 0}^{S} \Phi^{(2s+1)}(0) q^{2s+1} + O(q^{2S+3}) \quad (82) $$

La dérivée première en zéro est calculée dans l'annexe A et s'écrit en dimension quelconque :

$$ \Phi'(0) = \frac{\hbar \rho^{\text{sta}}}{4mc} \quad (83) $$

À température non nulle, les fonctions $\phi(q)$ dans (76,77) n'admettent donc dans leur développement que des puissances impaires de $q$. On en déduit que toutes leurs dérivées d'ordre pair $\phi^{(2s)}(q)$ sont nulles en $q=0$, et que l'écart $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ de $g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ à sa valeur asymptotique est à décroissance rapide (plus rapide que toute loi de puissance $1/|\mathbf{r}-\mathbf{r}'|^{2S+3}$). Comme c'était déjà le cas pour les prédictions $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ de la théorie BCS statique, voir les références [1, 2] et notre note 2,<fnref n="31" /> cette conclusion s'étend aux écarts $\delta g_{\sigma\sigma'}^{\text{dtd}}(\mathbf{r}, \mathbf{r}')$ en vertu des relations (56).

En revanche, à température nulle, la loi de Bose $\bar{n}_q$ disparaît dans (68) et $\delta\tilde{g}_{\text{coll}}(q)$ acquiert des puissances impaires de $q$, celles de la fonction $\Phi(q)$; il n'y a pas de contradiction avec le raisonnement conduisant à (78) car la condition de distance non infinitésimale aux singularités

<footnote n="26">²⁶La pulsation au carré n'a certes que des puissances paire, $\omega_q^2 = \omega_0^2 + c^2 q^2 + O(q^4)$, mais la prise de la racine carrée change la donne puisque $\omega_0 = 0$.</footnote>

<footnote n="27">²⁷À 3D, les fluctuations quantiques du champ bosonique des phonons, ignorées dans la théorie de champ classique qu'est BCS dépendant du temps linéarisée, font apparaître dans $\omega_q$ un terme en $q^5 \ln q$, et même un terme imaginaire pur en $q^5$ si la branche acoustique est de départ convexe ($\zeta_2 > 0$) [46, 47].</footnote>

<footnote n="28">²⁸On sait aussi en dimension quelconque que, dans la limite CBE où le spectre d'excitation acoustique se réduit à celui de Bogolioubov pour un condensat ou un quasi-condensat de dimères bosoniques, $\hbar\omega_q^{\text{Bog}} = [(\hbar^2 q^2/2m_d)(m_d c^2 + \hbar^2 q^2/2m_d)]^{1/2}$ avec $m_d = 2m$ la masse d'un dimère, on doit avoir $\zeta_s \to (-1)^s \Gamma(s-1/2)/[\Gamma(-1/2)s!4^{2s}]$ en vertu de l'équation (8.303) de la référence [5] et de la relation $\Gamma(1/2) = (-1/2)\Gamma(-1/2)$.</footnote>

<footnote n="29">²⁹Il y a un piège à éviter : en comparant l'équation (12) de la référence [42], $mc^2 = (2\mu/3)(1+xy)/(1+y^2)$ où $y = \delta\mu\Delta|a_{3D}|$, à notre équation (33), on pourrait croire que $y = z$. En réalité – on le voit bien à la limite unitaire – il faut prendre pour y la racine non évidente de l'équation $(1+xy)/(1+y^2) = (1+xz)/(1+z^2)$, à savoir $y = (x-z)/(1+xz)$; son report dans l'équation (16) de [42] donne notre expression (81).</footnote>

<footnote n="30">³⁰Voilà qui est cohérent avec (78) car, par la vertu du terme $1/2$, $\bar{n}_q + 1/2 = 1/2 \text{th}(\beta\hbar\omega_q/2)$ en facteur de $\Phi(q)$ dans $\delta\tilde{g}_{\text{coll}}^{\phi}(q)$ n'admet à température $T > 0$ que des puissances impaires dans son développement.</footnote>

<footnote n="31">³¹Les références [1, 2] sont à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ mais la propriété reste vraie au-delà; il suffit d'appliquer (76,77) aux fonctions $\delta\tilde{g}_{\sigma\sigma'}^{\phi}(q)$ des équations (60,61). D'ailleurs, notre note 2, qui va à tous les ordres en $\exp(-\beta E_{\text{gap}})$ dans sa seconde partie, le confirme explicitement.</footnote>
ne peut être satisfaite dans l'état fondamental du système lorsque $q \to 0$, voir l'équation (70). L'ordre dominant suffit ici ($S = 1$ dans (76,77)) :

$$ \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{T=0}{\sim}_{|\mathbf{r}-\mathbf{r}'|\to+\infty} \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') \stackrel{T=0}{\sim}_{|\mathbf{r}-\mathbf{r}'|\to+\infty} \begin{cases} -\frac{\Phi'(0)}{2\pi^2|\mathbf{r}-\mathbf{r}'|^4} & \text{en } d=3 \\ -\frac{\Phi'(0)}{4\pi|\mathbf{r}-\mathbf{r}'|^3} & \text{en } d=2 \end{cases} \quad (84) $$

Puisque $\Phi'(0) > 0$, la fonction de distribution de paires de fermions de spins opposés tend vers sa valeur asymptotique par en dessous, en contradiction avec la théorie BCS statique mais en accord avec la discussion physique de la section 1.

### 3.1.2. Ce qu'on déduit d'une formulation par intégrale de contour

Passons maintenant à une seconde technique d'étude à grande distance des intégrales de Fourier isotropes, utile lorsque le résultat est à décroissance rapide, donc à température non nulle, voir par exemple la référence [45] et la section 8.7.3.2 de la référence [5].

Commençons par le cas $d = 3$, plus facile. La fonction $\phi(q)$ dans (76) n'admet que des puissances impaires de $q$; on admet qu'elle peut être prolongée formellement en une fonction lisse et impaire de $q$ sur $\mathbb{R}$ tout entier. La fonction sinus étant impaire, l'intégrande est pair et l'on peut écrire

$$ \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') = \int_{-\infty}^{+\infty} \frac{dq}{2i} \phi(q) \exp(iq|\mathbf{r}-\mathbf{r}'|) \quad (85) $$

En admettant ensuite que $\phi(q)$ admet un prolongement analytique à $q \in \mathbb{C}$ en $O(1/q^\alpha)$ à l'infini, avec $\alpha > 1$,32 on transforme (85) en intégrale de contour en refermant le chemin d'intégration par un grand demi-cercle dans le demi-plan supérieur, avec contournement des éventuelles lignes de coupure. Le théorème de Cauchy permet enfin de décaler vers le haut la partie horizontale du chemin d'intégration, jusqu'à atteindre la première singularité $k_0$ de l'intégrande, que ce soit un pôle ou un point de branchement. On s'attend alors à ce que $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ décroisse à l'infini comme $\exp(-\kappa_0|\mathbf{r}-\mathbf{r}'|)$ avec $\kappa_0 = \text{Im } k_0 > 0$, à un facteur en loi de puissance près.33

Malheureusement, au contraire de la note 2 où nous avons pu appliquer jusqu'au bout cette technique à la théorie BCS statique, la fonction $\phi(q)$ mise en jeu ici n'a pas d'expression simple et ses singularités dans le plan complexe ne sont pas connues analytiquement. Pour avancer, nous nous plaçons donc dans la limite où la température tend vers zéro, en pratique dans le régime

$$ k_B T \ll mc^2, \Delta, E_F \quad (86) $$

Les singularités les plus proches de l'axe réel, qui conduisent donc aux exponentielles de $|\mathbf{r}-\mathbf{r}'|$ décroissant le plus lentement, sont alors les pôles de Bose $z_n = 2i\pi nk_B T (n \in \mathbb{N}^*)$ du nombre d'occupation thermique $\tilde{n}_q$ du mode de phonon $q$ dans l'expression (68), et le théorème des résidus donne finalement à 3D :34

$$ \delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{T=0^+}{\sim}_{|\mathbf{r}-\mathbf{r}'|\to+\infty} \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') \stackrel{T=0^+}{\sim}_{|\mathbf{r}-\mathbf{r}'|\to+\infty} -\frac{2\pi}{|\mathbf{r}-\mathbf{r}'|} \left(\frac{k_B T}{\hbar c}\right)^3 \Phi'(0) \sum_{n=1}^{+\infty} n^2 \exp\left(-\frac{2\pi n k_B T}{\hbar c} |\mathbf{r}-\mathbf{r}'|\right) \quad (87) $$

Comme à température nulle, le résultat est à valeurs négatives. Physiquement, pour entrer dans le régime des grandes distances, il faut que $|\mathbf{r}-\mathbf{r}'| \gg 1/k_F^{\text{sta}}$, mais aussi que $|\mathbf{r}-\mathbf{r}'| \gg \hbar/mc$ dans la limite CBE (c'est la longueur de relaxation du condensat de dimères $\xi$) et que $|\mathbf{r}-\mathbf{r}'| \gg \hbar k_F/m\Delta$

<footnote n="32">32 L'annexe B montre que $\Phi(q) = O(1/q^4)$ donc que $\phi(q) = O(1/q^3)$ sur l'axe réel.</footnote>

<footnote n="33">33 En particulier, on a sur le contour d'intégration la majoration $|\exp(iq|\mathbf{r}-\mathbf{r}'|)| \le \exp(-\kappa_0|\mathbf{r}-\mathbf{r}'|)$.</footnote>

<footnote n="34">34 La sommation peut être effectuée explicitement, $\sum_{n \ge 1} n^2 z^n = g_{-2}(z) = z(z+1)/(1-z)^3$ où $g_\alpha$ est une fonction de Bose ou polylogarithme.</footnote>
dans la limite BCS (c'est la taille d'une paire liée $\ell_{paire}$). En revanche, on n'a pas à mettre de contrainte sur le rapport entre $|\mathbf{r}-\mathbf{r}'|$ et la longueur d'onde thermique des phonons $\lambda_{\phi} = \hbar c/k_B T$. En particulier, les deux situations extrêmes sont bien décrites et continûment raccordées par (87), ce qui en fait toute l'originalité : (i) le cas où $|\mathbf{r}-\mathbf{r}'|$ diverge plus vite que $\lambda_{\phi}$ et où la somme est dominée par le terme $n=1$, qui donne à un facteur en loi de puissance près la décroissance exponentielle attendue avec $\kappa_0 = 2\pi k_B T/\hbar c$ (régime asymptotique thermique), et (ii) le cas où $|\mathbf{r}-\mathbf{r}'|$ diverge plus lentement que $\lambda_{\phi}$ et où la somme sur $n$ peut être remplacée par une intégrale (régime asymptotique quantique). Dans le second cas, les fonctions de distribution de paires ne voient pas que le système est à température non nulle, et l'on retrouve exactement l'équivalent algébrique (84) en $d=3$.<fnref n="35" />

Terminons par le cas $d=2$, plus technique. La difficulté est que $J_0(q|\mathbf{r}-\mathbf{r}'|)$ est une fonction paire de $q$, ce qui rend l'intégrande de (77) impair et empêche d'étendre l'intégration à $\mathbb{R}$ tout entier. Nous nous sortons de ce mauvais pas en reprenant l'idée de la section 8.7.3.2 de la référence [5] et en l'étendant à tous les ordres en $1/|\mathbf{r}-\mathbf{r}'|$. Pour cela, nous remplaçons la fonction de Bessel par son développement asymptotique

$$J_0(q|\mathbf{r}-\mathbf{r}'|) \underset{|\mathbf{r}-\mathbf{r}'| \to +\infty}{\approx} \left( \frac{2}{\pi q|\mathbf{r}-\mathbf{r}'|} \right)^{1/2} \left[ \cos\left(q|\mathbf{r}-\mathbf{r}'| - \frac{\pi}{4}\right) u(q|\mathbf{r}-\mathbf{r}'|) - \sin\left(q|\mathbf{r}-\mathbf{r}'| - \frac{\pi}{4}\right) \frac{v(q|\mathbf{r}-\mathbf{r}'|)}{q|\mathbf{r}-\mathbf{r}'|} \right] \quad (88)$$

où les fonctions $u(x)$ et $v(x)$ sont données dans l'alinéa §8.451(1) de la référence [21] sous forme de séries ne présentant que des puissances paires de $1/x$. On montre alors comme dans la référence [5] que, pour toute fonction $f_-(q)$ impaire sur $\mathbb{R}$ et pour toute fonction $f_+(q)$ paire sur $\mathbb{R}$, on a grâce à la racine carrée $q^{1/2}$:<fnref n="36" />

$$\int_0^{+\infty} \frac{dq}{q^{1/2}} f_-(q) \cos(q|\mathbf{r}-\mathbf{r}'| - \pi/4) = e^{-i\pi/4} \int_{-\infty+i0^+}^{+\infty+i0^+} \frac{dz}{2z^{1/2}} f_-(z) e^{iz|\mathbf{r}-\mathbf{r}'|} \quad (89)$$

$$\int_0^{+\infty} \frac{dq}{q^{1/2}} f_+(q) \sin(q|\mathbf{r}-\mathbf{r}'| - \pi/4) = -e^{i\pi/4} \int_{-\infty+i0^+}^{+\infty+i0^+} \frac{dz}{2z^{1/2}} f_+(z) e^{iz|\mathbf{r}-\mathbf{r}'|} \quad (90)$$

ce qu'il reste à appliquer à $f_-(q) = \phi(q)u(q|\mathbf{r}-\mathbf{r}'|)$ et à $f_+(q) = \phi(q)v(q|\mathbf{r}-\mathbf{r}'|)/q$. La suite ressemble au 3D. Il faut seulement penser à remarquer que $u(ix) + v(ix)/x = (2x/\pi)^{1/2} \exp(x)K_0(x)$, en vertu de la relation §8.451(6) de la référence [21]. On trouve en définitive à 2D que

$$\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \underset{|\mathbf{r}-\mathbf{r}'| \to +\infty}{\stackrel{T \to 0^+}{\approx}} \delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}') \underset{|\mathbf{r}-\mathbf{r}'| \to +\infty}{\stackrel{T \to 0^+}{\approx}} -4\pi \left(\frac{k_B T}{\hbar c}\right)^3 \Phi'(0) \sum_{n=1}^{+\infty} n^2 K_0\left(\frac{2\pi n k_B T}{\hbar c} |\mathbf{r}-\mathbf{r}'|\right) \quad (91)$$

Comme à 3D, cet approximant, à notre connaissance original, est à valeurs négatives (on a $\Phi'(0) > 0$ et $K_0(x) > 0 \forall x > 0$) et raccorde continûment le régime asymptotique thermique au régime asymptotique quantique. En particulier, (i) en limitant la somme au terme $n=1$ et en remplaçant la fonction $K_0$ par un équivalent asymptotique $K_0(x) \sim (\pi/2x)^{1/2} \exp(-x)$, on retrouve la décroissance exponentielle attendue (à un facteur en loi de puissance près)

$$\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \underset{\text{puis } T \to 0^+}{\stackrel{|\mathbf{r}-\mathbf{r}'| \to +\infty}{\sim}} - \frac{2\pi}{|\mathbf{r}-\mathbf{r}'|^{1/2}} \left(\frac{k_B T}{\hbar c}\right)^{5/2} \Phi'(0) \exp(-\kappa_0 |\mathbf{r}-\mathbf{r}'|) \quad \text{avec} \quad \kappa_0 = 2\pi k_B T / \hbar c; \quad (92)$$

(ii) en remplaçant dans (91) la somme sur $n$ par une intégrale, on retrouve exactement l'équivalent algébrique à température nulle (84) en $d=2$.<fnref n="37" />

<footnote n="35">³⁵On donne $\int_0^{+\infty} dx x^2 \exp(-x) = 2$.</footnote>

<footnote n="36">³⁶La première relation est l'équation (8.349) de [5]; la seconde s'en déduit par dérivation par rapport à $|\mathbf{r}-\mathbf{r}'|$. L'idée intuitive est que la racine carrée d'un nombre réel négatif contient un facteur $\exp(i\pi/2)$, c'est-à-dire précisément la phase relative entre les composantes de Fourier de $\cos(|\mathbf{r}-\mathbf{r}'| - \pi/4)$ ou de $\sin(q|\mathbf{r}-\mathbf{r}'| - \pi/4)$ (au signe près pour le sinus).</footnote>

<footnote n="37">³⁷On donne $\int_0^{+\infty} dx x^2 K_0(x) = \pi/2$.</footnote>
### 3.1.3. *Ce qu'en dit l'hydrodynamique quantique*

Prenons maintenant un peu de recul, et demandons-nous dans quelle mesure les beaux approximants (87,91) sont spécifiques de la théorie BCS dépendant du temps linéarisée. Il apparaît que, dans la double limite des basses températures $T \to 0^+$ et des grandes distances $|\mathbf{r}-\mathbf{r}'| \to +\infty$ considérée, on peut faire usage d'une théorie effective de basse énergie, l'hydrodynamique quantique de Landau et Khalatnikov [48]. On s'attend à ce que cette théorie soit exacte à l'ordre dominant en température pour tout système – même fortement corrélé comme l'hélium liquide – (i) siège d'interactions à courte portée et (ii) entièrement superfluide dans son état fondamental. Notre gaz de fermions de spin 1/2 non polarisé et en interaction de contact satisfait à ces deux conditions, voir plus de détails dans l'annexe B de la référence [12] et une explication pédagogique dans la référence [49]. Alors, (i) le fluide admet dans son état fondamental une branche d'excitation acoustique de départ linéaire $\omega_q \sim cq$, avec une vitesse $c$ reliée à son équation d'état exacte par la relation hydrodynamique (32); (ii) à suffisamment basse température, on peut assimiler l'état d'équilibre des phonons, qui sont les quanta de ces excitations acoustiques, à celui d'un gaz parfait de bosons de potentiel chimique nul $\mu_\phi = 0$ donc d'opérateur densité à N corps $\hat{\sigma}_\phi \propto \exp(-\beta \sum_q \hbar\omega_q \hat{b}_q^\dagger \hat{b}_q)$; (iii) dans une description à gros grains, les fluctuations de densité du fluide à valeur opérateur sont dominées – toujours à suffisamment basse température – par les phonons et admettent le développement modal

$$ \delta\hat{\rho}(\mathbf{r}) = \frac{1}{L^{d/2}} \sum_{\mathbf{q} \neq 0} \rho_{\mathbf{q}} (\hat{b}_{\mathbf{q}} + \hat{b}_{-\mathbf{q}}^{\dagger}) \exp(i\mathbf{q} \cdot \mathbf{r}) \quad \text{avec} \quad \rho_{\mathbf{q}} = \left(\frac{\hbar\rho q}{2mc}\right)^{1/2} \qquad (93) $$

Ici, $\rho$ est la densité moyenne exacte du fluide dans son état fondamental, $m$ est la masse d'une particule du fluide et les opérateurs de création $\hat{b}_q^\dagger$ et d'annihilation $\hat{b}_q$ d'un phonon de vecteur d'onde $\mathbf{q}$ obéissent aux relations de commutation bosoniques, par exemple $[b_q, \hat{b}_q^\dagger] = \delta_{qq'}$. D'où la fonction de corrélation de la densité totale à la limite thermodynamique dans cette théorie,

$$ \langle \delta\hat{\rho}(\mathbf{r})\delta\hat{\rho}(\mathbf{r}') \rangle_{\text{hydro}} = \int_{q<\Lambda} \frac{d^d q}{(2\pi)^d} 2\rho_{\mathbf{q}}^2 \left( \tilde{n}_{\mathbf{q}} + \frac{1}{2} \right) \exp[i\mathbf{q} \cdot (\mathbf{r}-\mathbf{r}')] \qquad (94) $$

où la coupure ultraviolette $\Lambda$ est $\gg k_B T/\hbar c$ mais $\ll mc/\hbar$. Or, les deux états internes $\sigma' = \uparrow, \downarrow$ jouant des rôles symétriques pour un spin $\sigma$ très distant, et la dissociation thermique des paires $\uparrow\downarrow$ étant exponentiellement négligeable dans notre gaz de fermions à basse température, on doit avoir à grande distance (à distance en tout cas beaucoup plus grande que la taille d'une paire liée, $|\mathbf{r}-\mathbf{r}'| \gg \ell_{\text{paire}}$):

$$ \delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}') \approx \frac{1}{4} \langle \delta\hat{\rho}(\mathbf{r})\delta\hat{\rho}(\mathbf{r}') \rangle_{\text{hydro}} \qquad (95) $$

La comparaison à (68) impose

$$ \Phi(\mathbf{q}) \sim \frac{\hbar\rho q}{4mc} \qquad (96) $$

L'expression (83) de $\Phi'(0)$ et les approximants asymptotiques (84,87,91) qui en dérivent sont donc exacts et très généraux, pourvu que l'on y remplace la densité et la vitesse du son BCS par leurs vraies valeurs dans le gaz.

En conclusion, dans notre gaz de fermions non polarisé à suffisamment basse température, les fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ sont en certains points, en particulier à suffisamment grande distance $|\mathbf{r}-\mathbf{r}'|$, strictement en dessous de leurs valeurs asymptotiques. Cette propriété, en contradiction directe avec la théorie BCS statique pour $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$, voir l'équation (2), résulte des fluctuations quantiques ou thermiques des positions des centres de masse des paires liées induites par les phonons, grands absents de cette théorie. Elle vaut aussi bien pour la théorie BCS dépendant du temps linéarisée que pour la réalité physique.
## 3.2. Résultats numériques

Nous calculons numériquement les fonctions de distribution de paires $g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ de la section 2 au moyen de leur décomposition (56), puis de la représentation (65) de la partie collective $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ – en réalité de sa transformée de Fourier $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ – par intégrale curviligne dans le plan complexe. Pour simplifier, nous nous plaçons à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ c'est-à-dire que nous omettons les facteurs thermiques $f_{\pm}$ dans les fonctions $\Sigma_{ij}(\mathbf{q}, z)$ donc dans la susceptibilité $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$.<fnref n="38" />

Pour tester le calcul numérique de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$, nous comparons aux comportements aux limites (157) et (109) obtenus analytiquement; dans les exemples considérés ici, qui sont tous à très basse température, la limite à faible $q$ est négative comme celle à grand $q$, et nous trouvons que $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ raccorde de façon monotone l'une à l'autre, en restant donc partout négative :

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q})^{\text{num}} < 0 \quad \forall \mathbf{q} \qquad (97) $$

Les motivations de notre étude numérique sont doubles : (i) il s'agit de voir si $\delta g_{\updownarrow}(\mathbf{r}, \mathbf{r}')/\rho_{\uparrow}\rho_{\downarrow}$ peut prendre des valeurs négatives significatives (par rapport à l'unité) dans la théorie BCS dépendant du temps linéarisée, question à laquelle l'étude asymptotique de la section 3.1 ne peut répondre; (ii) il s'agit de comparer cette théorie à d'autres approches ou à des mesures dans un gaz d'atomes froids fermioniques.

### 3.2.1. Cas tridimensionnel; comparaison au Monte-Carlo quantique de Giorgini

En dimension $d=3$, nous nous plaçons à la limite unitaire $a_{3D}^{-1} = 0$; c'est le régime d'interaction maximale accessible en phase gazeuse, donc a priori le plus difficile à décrire, mais c'est aussi le plus étudié, théoriquement et expérimentalement dans les gaz d'atomes froids, et le plus

<footnote n="38">³⁸ Dans le numérique, nous avons introduit une coupure isotrope suffisamment grande $k < k_{\text{max}}$ sur le vecteur d'onde interne $\mathbf{k}$ dans les fonctions $\Sigma_{ij}(\mathbf{q}, z)$, avec $k_{\text{max}} = 273(m\Delta/\hbar^2)^{1/2}$ sur la figure 3a et $k_{\text{max}} = 157q_{\text{dim}}$ sur la figure 3b. L'intérêt est aussi que la ligne de coupure commençant à $E_{\text{bord}}(\mathbf{q})$ ne va pas au-delà de $\epsilon_{\text{max}} \approx \hbar^2(k_{\text{max}}^2 + q^4/4) - 2\mu$, puisque les dénominateurs d'énergie dans les $\Sigma_{ij}(\mathbf{q}, z)$ ne peuvent plus s'annuler si $z = \epsilon > 0$ est trop grand; le chemin d'intégration $C_h$ peut alors rejoindre l'axe réel par un segment vertical entre $z = \epsilon_{\text{max}} + i\eta$ (avec $\eta = \pi k_B T$ comme il est dit après l'équation (65)) et le point terminal $z = \epsilon_{\text{max}}$. Quant à la ligne de coupure aboutissant à $-E_{\text{bord}}(\mathbf{q})$, il est inutile de l'explorer en deçade de $\text{Re } z = -2\Delta$ à cause de la suppression exponentielle par le dénominateur thermique dans (65). Dans l'intégration sur $\epsilon$, nous avons pris un pas $\delta\epsilon = 0,2 \min(\max(|\mu|, \Delta), \Delta, \hbar c q, k_B T, \eta)$ (voir à ce propos la note 23). Dans l'intégration sur le module de $\mathbf{k}$, nous avons pris un pas $\delta\mathbf{k} = 0,2 \min(\max((m\mu)^{1/2}/\hbar, (m\Delta)^{1/2}/\hbar), m\Delta/\hbar^2/(m\mu/\hbar^2)^{1/2})$ pour $\mu > 0$ et – cas non exploré – $\delta\mathbf{k} = 0,2 \max((-m\mu)^{1/2}/\hbar, (m\Delta)^{1/2}/\hbar)$ pour $\mu \le 0$. Dans l'intégrale sur l'angle $\theta$ entre $\mathbf{k}$ et $\mathbf{q}$, réduit à l'intervalle $[0, \pi/2]$ par symétrie, nous avons pris 30 points. Nous avons vérifié que les résultats changent peu si l'on divise par deux l'un des pas d'intégration $\delta\epsilon, \delta k$ ou $\delta\theta$. Dans le calcul numérique de la transformée de Fourier inverse donnant $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$ à partir de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$, nous mettons une coupure $q_{\text{max}} = k_{\text{max}}^{1/2}$ sur le nombre d'onde $q$, et nous prenons un pas $\delta q = 0,1(m\Delta/\hbar^2)^{1/2}$ à 3D et $\delta q = 0,13q_{\text{dim}}$ à 2D. La coupure $q_{\text{max}}$ choisie serait très insuffisante pour $g_{\updownarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ si l'on calculait directement la transformée de Fourier inverse de $\delta \tilde{g}_{\updownarrow}^{\text{dtl}}(\mathbf{q})$ sans utiliser la décomposition (56); en effet, la contribution $\delta g_{\updownarrow}^{\text{sta}}(\mathbf{q})$ de la théorie BCS statique décroît lentement avec $q$, voir l'équation (122), alors que la contribution collective $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ décroît comme $1/q^4$, voir l'équation (109). La partie $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ s'obtient facilement au travers des moyennes anormale et anormale comme dans (23,24); à 2D, ces moyennes sont connues analytiquement, voir les équations (30,31); à 3D, elles sont données par les intégrales sur $\mathbf{k}$ (25,26) dont l'évaluation numérique poserait aussi un problème de troncature numérique sur $\mathbf{k}$ [2] si l'intégrande n'était pas connu analytiquement : on peut évaluer le bout omis numériquement au moyen d'un développement asymptotique de l'intégrande, ou carrément utiliser une représentation astucieuse par intégrale de contour, voir notre note 2 et la référence [2]. Enfin, à 3D, la précision numérique sur $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ à grand $q$ ne suffit pas à bien décrire le départ de $g_{\updownarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ à faible $|\mathbf{r} - \mathbf{r}'|$; nous préférons calculer le bout $q > 15(m\Delta/\hbar^2)^{1/2}$ de la transformée de Fourier inverse après la substitution $\delta \tilde{g}_{\text{coll}}(\mathbf{q}) \approx C/q^4 + D/q^5 + E/q^6$; les coefficients $C$ et $D$ sont donnés dans l'équation (110), le coefficient $E$ s'obtient par ajustement.</footnote>
<figure><img src="image_3.png" /></figure>

FIGURE 3 – Dans un gaz non polarisé de fermions de spin 1/2, fonctions de distribution de paires $g_{\uparrow\downarrow}$ (rangée du haut) et $g_{\downarrow\uparrow}$ (rangée du bas) en dimension $d = 3$ à la limite unitaire (colonne de gauche) et en dimension $d = 2$ pour une interaction de paramètre $\ln(k_F/q_{\text{dim}}) = 1$ (colonne de droite). Tireté : calcul des écarts $\delta g_{\sigma\sigma'}$ des $g_{\sigma\sigma'}$ à leurs valeurs asymptotiques par la théorie BCS dépendant du temps linéarisé, puis complétion naïve (71) des $\delta g_{\sigma\sigma'}$ en les $g_{\sigma\sigma'}$ au moyen de la densité de la théorie BCS statique. Trait plein noir : même calcul des $\delta g_{\sigma\sigma'}$ dans la théorie BCS dépendant du temps linéarisée, mais complétion correcte (73) au moyen de la densité modifiée $\rho^{\text{dtl}}$ dans cette théorie. Trait plein rouge : prédictions analytiques (126,129) pour le départ de $g_{\uparrow\downarrow}$, toujours dans cette théorie. La densité totale $\rho = \rho_{\uparrow} + \rho_{\downarrow} = 2\rho_{\sigma}$ et le nombre d'onde de Fermi $k_F$ (22) qui en découle sont considérés comme fixés sur la figure; le potentiel chimique $\mu$ dans notre théorie grand-canonique n'a donc pas même valeur pour les courbes en tireté (pour lesquelles $\rho = \rho^{\text{sta}}(\mu, T)$) et en trait plein (pour lesquelles $\rho = \rho^{\text{dtl}}(\mu, T)$). Disques noirs : résultats d'un calcul de Monte-Carlo quantique 3D par diffusion à surface nodale fixée (température nulle) [50, 51] lissés par ajustement polynomial quintique pour la colonne de gauche; mesures expérimentales au moyen d'un microscope dans un gaz bidimensionnel d'atomes froids [52] pour la colonne de droite. La théorie BCS dépendant du temps linéarisée est mise en œuvre à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, à 3D à une température $T = 0,144\Delta/k_B \approx 0,251mc^2/k_B$ (assez basse pour permettre une comparaison à la température nulle) et à 2D à une température $T = 0,09T_F$ (comme dans l'expérience).

intéressant car il est invariant d'échelle. Prenons dans le numérique une température suffisamment basse pour pouvoir nous placer à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$,

$$ \frac{k_B T}{\Delta} = 0,144 \qquad (98) $$

où $\Delta$ est le paramètre d'ordre BCS. Le potentiel chimique étant positif, on a alors $E_{\text{gap}} = \Delta$ et $\exp(-\beta E_{\text{gap}}) \approx 10^{-3} \ll 1$. Des équations (27,28) de la théorie BCS statique et de la définition (22)
du nombre d'onde de Fermi, nous tirons $\mu/\Delta = 0,860\,436...$ et $\mu/E_F^{\text{sta}} = 0,590\,605...$. Le rapport $\mu/E_{\text{gap}}$ prédit est proche de la valeur mesurée aux températures les plus basses accessibles mais le rapport $\mu/E_F$ en est assez loin, puisqu'on trouve dans les expériences $E_{\text{gap}}/E_F \approx 0,44$ [53] et $\mu/E_F \approx 0,376(4)$ [54], ce qui conduit à $\mu/\Delta \approx 0,86$ (en admettant qu'on ait aussi $E_{\text{gap}} = \Delta$ dans la réalité physique). Heureusement, l'équation d'état modifiée (73) de la théorie BCS dépendant du temps linéarisée fait nettement mieux : elle ne change certes pas le rapport $\mu/\Delta$ mais donne

$$ \frac{\mu}{E_F^{\text{dtl}}} = 0,401 \qquad (99) $$

ce qui, par ailleurs, est en parfait accord avec l'approximation des fluctuations gaussiennes, voir les références [28,43] et notre note 20. Notre choix de température correspond alors à $k_B T/mc^2 \approx 0,251 \ll 1$, si bien que l'effet des phonons thermiques devrait être assez faible pour que l'équation d'état et les fonctions de distribution de paires $g_{\sigma\sigma'}$ soient proches de leurs valeurs dans l'état fondamental (l'invariance d'échelle impose $mc^2 = 2\mu/3$ dans l'expression hydrodynamique (32) de la vitesse du son).<fnref n="39" /> Il correspond aussi à $T/T_F^{\text{dtl}} \approx 0,067$, donc à une température nettement plus basse que celle $T \approx 0,1 T_F$ de la référence [54] mais encore supérieure à celle $T \approx 0,05 T_F$ de la référence [55], le record dans un gaz d'atomes froids fermioniques.<fnref n="40" />

Sur la colonne de gauche de la figure 3, nous représentons les fonctions de distribution de paires réduites $g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')/\rho_\sigma\rho_{\sigma'} = 1+\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')/(\rho/2)^2$ en fonction de la distance adimensionnée $k_F|\mathbf{r}-\mathbf{r}'|$, avec $k_F = (3\pi^2\rho)^{1/3}$ comme dans (22). La figure se place dans le point de vue d'une densité $\rho$ fixée, comme il est dit dans la légende. Cependant, sa colonne de gauche ne dépend pas de la valeur précise de $\rho$ choisie parce qu'elle se trouve à la limite unitaire, invariante d'échelle; de plus, comme on y a fixé le rapport $k_B T/\Delta$ ou, ce qui revient au même ici, le rapport $k_B T/\mu$, plutôt que $T/T_F$, nous jugeons plus commode d'en effectuer la discussion à potentiel chimique $\mu$ fixé (il en ira différemment pour la colonne de droite). Après détermination numérique des écarts $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ dans la théorie BCS dépendant du temps linéarisée grand-canonique (le gros morceau

<footnote n="39">³⁹ Les phonons changent le grand potentiel du système d'une quantité $\delta\Omega_\phi = k_B T \sum_q \ln[1-\exp(-\beta\hbar\omega_q)]$ donc la densité du gaz de fermions grand-canonique d'une quantité $\delta\rho_\phi = -\partial_\mu\delta\Omega_\phi/L^d \sim -\int_{\mathbb{R}^d} [d^d q/(2\pi)^d]\hbar(\partial_\mu c)q/[\exp(\beta\hbar c q) - 1] = -[\Omega_d/(2\pi)^d]\hbar(\partial_\mu c)(k_B T/\hbar c)^{d+1}\Gamma(d+1)\zeta(d+1)$ où la relation de dispersion $\omega_q$ est prise linéaire (10) à suffisamment basse température $T \to 0^+$, $\Omega_d$ est l'angle total de l'espace $(2\pi \text{ si } d=2, 4\pi \text{ si } d=3)$, $\Gamma$ est la fonction Gamma d'Euler et $\zeta$ est la fonction zêta de Riemann. À 3D à la limite unitaire, nous en déduisons que $\delta\rho_\phi/\rho \sim -(\pi^4/30)(k_B T/mc^2)^4(\xi_B/3)^{3/2}$ avec $\xi_B = \mu/E_F$ le paramètre de Bertsch. Pour nos paramètres, ceci donne dans la théorie BCS dépendant du temps linéarisée $\delta\rho_\phi/\rho \approx -6 \cdot 10^{-4} \ll 1$. Pour voir si l'équation d'état modifiée (73) donne la bonne valeur de $\delta\rho_\phi$ à basse température, supposons que l'intégrale sur $q$ dans (73) soit une petite correction à $\rho^{\text{sta}}$; dans ce cas, on peut écrire $\rho^{\text{dtl}} = \rho^{\text{sta}} - (2/\rho^{\text{sta}}) \int_{\mathbb{R}^d} [d^d q/(2\pi)^d] \delta\tilde{g}_{\text{coll}}(q)$. On trouve alors que $\delta\rho_\phi^{\text{dtl}} = -\int_{\mathbb{R}^d} [d^d q/(2\pi)^d](2\Phi'(0)/\rho^{\text{sta}})q/[\exp(\beta\hbar c q) - 1]$; vu l'expression (83) de $\Phi'(0)$, il y a accord si $\hbar\partial_\mu c = \hbar/2mc$ donc si $\partial_\mu(mc^2) = 1$. Cette condition n'est en général pas parfaitement satisfaite, sauf à 2D, voir l'équation (34). Dans le cas où $\rho^{\text{dtl}}$ diffère significativement de $\rho^{\text{sta}}$, l'expression précédente de $\delta\rho_\phi^{\text{dtl}}$ doit être multipliée par $\rho^{\text{sta}}/\rho^{\text{dtl}}$ et la condition devient $(\rho^{\text{dtl}}/\rho^{\text{sta}})\partial_\mu(mc^2) = 1$, ce qui donne à 3D à la limite unitaire $(2/3)(\xi_B^{\text{sta}}/\xi_B^{\text{dtl}})^{3/2} = 1$ soit 1,19#1.</footnote>

<footnote n="40">⁴⁰ Pour $T/T_F \approx 0,1$ et $\xi_B \approx 0,376$ comme dans [54], on trouve un faible effet des phonons $\delta\rho_\phi/\rho \approx -0,0036$, voir notre note 39, mais un facteur de dissociation des paires $\exp(-\beta E_{\text{gap}}) \approx 0,013$ a priori pas complètement négligeable. Pour en avoir le cœur net, nous avons déterminé le changement de densité $\delta\rho_\gamma$ du gaz de fermions grand-canonique dû aux quasi-particules fermioniques d'hamiltonien (13). En procédant comme dans la note 39 et en approximant $\epsilon_k$ au voisinage de son minimum par $\epsilon_k \approx E_{\text{gap}} + \hbar^2(k-k_0)^2/2m_*$, ce qui est légitime à basse température $T \to 0^+$, nous trouvons pour $k_0 > 0$ que $\delta\rho_\gamma = (-2) \int_{\mathbb{R}^d} [d^d k/(2\pi)^d] \partial_\mu \epsilon_k / [1 + \exp(\beta\epsilon_k)] \frac{d^3}{T-0^+} (-2(\partial_\mu E_{\text{gap}}) \exp(-\beta E_{\text{gap}})) k_0^{2/3} (2\pi m_* k_B T/\hbar^2)^{1/2} / (2\pi^2)$. À la limite unitaire, où $\partial_\mu E_{\text{gap}} = E_{\text{gap}}/\mu$ par invariance d'échelle, les valeurs mesurées $E_{\text{gap}} = 0,44 E_F$, $k_0 \approx 0,92 k_F$ [53] et la prédiction théorique $m_*/m \approx 0,56$ [56] donnent $\delta\rho_\gamma/\rho \approx -1,4 \cdot 10^{-3}$ à la température de notre calcul numérique $T = 0,067 T_F$ mais $\delta\rho_\gamma/\rho \approx -0,015$ à la température expérimentale $T = 0,1 T_F$. Dans l'expérience, ceci entraîne une surestimation du paramètre de Bertsch $\xi_B = (2m\mu/\hbar^2)/(3\pi^2\rho)^{2/3}$ de $\delta\xi_B/\xi_B = (-2/3)\delta\rho_\gamma/\rho \approx 0,01$ en valeur relative, encore dans les barres d'erreur de la référence [54].</footnote>
est bien entendu $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$), il faut tout de même calculer la densité $\rho$ (et donc le nombre d'onde de Fermi $k_F$) pour former les quantités adimensionnées sur l'axe horizontal et l'axe vertical. La courbe en tireté correspond au choix tentant mais incorrect de l'équation d'état BCS statique $\rho = \rho^{\text{sta}}(\mu, T)$ (21), la courbe en trait plein au choix correct d'une équation d'état corrigée par les excitations collectives du superfluide, $\rho = \rho^{\text{dtl}}(\mu, T)$ (73), comme l'explique la section 2.5 (on se met dans les deux cas à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$, ce qui revient formellement à prendre $\rho^{\text{sta}}$ et la susceptibilité $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ à température nulle). Commentons brièvement ces résultats.

Dans le cas $\sigma = \sigma' = \uparrow$ (rangée du bas), la courbe en tireté est clairement catastrophiquement fausse, car elle conduit à des valeurs négatives de $g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}')$ sur un voisinage de la distance nulle; or, on voit bien sur sa définition (1) que $g_{\uparrow\uparrow}(\mathbf{r}, \mathbf{r}')$ est positive – elle est de la forme $\langle \hat{A}^\dagger \hat{A} \rangle$ où $A$ est un opérateur. La courbe en trait plein est bien plus convaincante : elle tend vers zéro à distance nulle, comme il se doit pour des fermions identiques $\uparrow\uparrow$ sans interaction directe, et elle reproduit l'allure générale des résultats Monte-Carlo (disques noirs). Elle en diffère cependant de manière spectaculaire par son départ linéaire en $|\mathbf{r} - \mathbf{r}'|$, plutôt que quadratique dans le Monte-Carlo [50]. Comme nous le verrons dans notre section 3.3, ce départ linéaire n'est pas l'apanage de la théorie BCS dépendant du temps linéarisée; il existe déjà dans la théorie BCS statique (avec une pente plus faible il est vrai), un fait passé sous silence dans les références [50, 51].

Dans le cas $\sigma = \uparrow, \sigma' = \downarrow$, la courbe en tireté, en accord acceptable avec les résultats Monte-Carlo, nous fait une fausse joie, tandis que la courbe en trait plein, avec son creux bien moins prononcé dans $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')/\rho_{\uparrow}\rho_{\downarrow}$, nous procure une certaine déception. Pour comprendre pourquoi un changement somme toute assez modeste de l'équation d'état ($\mu/E_F$ passe de 0,59 dans BCS statique à 0,40 dans BCS dépendant du temps linéarisé) affecte autant le creux, écrivons comme suit sa profondeur (comptée négativement) sur la figure 3 :

$$ \min_{|\mathbf{r}-\mathbf{r}'|} \frac{g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')}{\rho_{\uparrow}\rho_{\downarrow}} - 1 = \frac{(3\pi^2)^2}{2} \left(\frac{\mu}{E_F}\right)^3 \min_{|\mathbf{r}-\mathbf{r}'|} \frac{\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')}{(m\mu/\hbar^2)^3} \quad (100) $$

où la densité a été exprimée en termes du nombre d'onde de Fermi. La valeur du minimum au second membre, fonction des seules variables thermodynamiques $\mu$ et $T$, plus précisément à la limite unitaire de leur rapport $k_B T/\mu$ ici fixé, est indépendante de l'équation d'état; en revanche, son préfacteur, proportionnel au cube du rapport $\mu/E_F$, prend mécaniquement un facteur $(0,40/0,59)^3 \approx 0,31$ lorsqu'on passe du choix trop naïf $\rho = \rho^{\text{sta}}$ au bon choix $\rho = \rho^{\text{dtl}}$. La profondeur relative passe ainsi de $-0,16$ à $-0,05$, qui est la vraie prédiction de notre théorie, à comparer à $-0,2$ pour le Monte-Carlo. La même analyse montre pourquoi la position du minimum $|\mathbf{r}-\mathbf{r}'|_0$ change au contraire assez peu (il est légèrement décalé vers la droite) : dans l'écriture

$$ k_F |\mathbf{r} - \mathbf{r}'|_0 = \left(\frac{E_F}{\mu}\right)^{1/2} \left[ (2m\mu/\hbar^2)^{1/2} |\mathbf{r} - \mathbf{r}'|_0 \right] \quad (101) $$

la quantité entre crochets est indépendante de l'équation d'état, et le préfacteur est multiplié seulement par $(0,59/0,40)^{1/2} \approx 1,21$ entre la courbe en tireté et la courbe en trait plein. La position réduite $k_F|\mathbf{r} - \mathbf{r}'|_0$ du minimum passe ainsi de 1,82 à 2,21, qui est la vraie prédiction de notre théorie, à comparer à 2,05 pour le Monte-Carlo.

Une autre façon de raisonner est de faire le chemin inverse, c'est-à-dire de grand-canoniser les résultats des références [50,51] en les exprimant dans l'unité de nombre d'onde grand-canonique naturelle $(m\mu/\hbar^2)^{1/2}$. Comme on a $\mu/E_F \approx 0,42$ dans ces références, on trouve par inversion de la relation (100) que

$$ \min_{|\mathbf{r}-\mathbf{r}'|} \frac{\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')}{(m\mu/\hbar^2)^3} \approx \begin{cases} -6,0 \cdot 10^{-3} & \text{dans le calcul Monte-Carlo de [50,51]} \\ -1,8 \cdot 10^{-3} & \text{dans notre théorie BCS dépendant du temps linéarisée} \end{cases} \quad (102) $$
Le désaccord est flagrant, sans qu'il y ait à débattre de l'équation d'état, et l'accord acceptable entre la courbe en tireté et les points Monte-Carlo sur le panneau supérieur gauche de la figure 3 n'est qu'une illusion.

### 3.2.2. Cas bidimensionnel; comparaison aux expériences de l'ENS

En dimension $d = 2$, nous disposons de mesures expérimentales récentes des fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ dans un gaz spatialement homogène d'atomes froids fermioniques en interaction [52, 57], aussi bien pour des spins égaux que des spins opposés, grâce à une technique de microscopie initialement mise au point pour des atomes dans un réseau optique [58–60] puis étendue avec succès à un espace continu [61] (voir la référence [62] pour les bosons). Choisissons là aussi une interaction assez forte, sinon maximale : à 2D, il n'y a pas d'équivalent strict de la limite unitaire mais on peut avoir une longueur de diffusion $a_{2D}$ ou, si l'on préfère, une taille de dimère dans l'espace libre $q_{\dim}^{-1}$ de l'ordre de la distance moyenne entre particules; parmi les forces d'interaction disponibles dans [52], prenons donc$^{41}$

$$ \eta \equiv \ln(k_F/q_{\dim}) = 1 \quad \text{c'est-à-dire} \quad \left\{ \begin{array}{l} \rho q_{\dim}^{-2} = \frac{1}{2\pi} \exp(2\eta) \approx 1,1760 \\ \\ \frac{\epsilon_{\dim}}{E_F} = 2 \exp(-2\eta) \approx 0,271 \end{array} \right. \qquad (103) $$

où le nombre d'onde $q_{\dim}$ est relié à $a_{2D}$ et à l'énergie de liaison du dimère $\epsilon_{\dim}$ comme dans (12), et où $k_F$ est le nombre d'onde de Fermi (22). Nous prenons la même température que dans l'expérience,

$$ \frac{T}{T_F} = 0,09 \quad \text{soit} \quad \frac{k_B T}{\epsilon_{\dim}} = \frac{1}{2} \exp(2\eta) \frac{T}{T_F} \approx 0,3325 \qquad (104) $$

Elle a le bon goût de conduire à un taux de dissociation thermique des paires liées $\uparrow\downarrow$ négligeable,

$$ \exp(-\beta E_{\text{gap}}) \approx 10^{-3} \qquad (105) $$

si l'on se souvient de l'expression (29) du paramètre d'ordre et si l'on anticipe sur la valeur (108) du potentiel chimique, ce qui nous permet de restreindre la théorie à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$. Contrairement à la limite unitaire tridimensionnelle, la densité $\rho$ et le nombre d'onde de Fermi $k_F$ sont maintenant fixés par le choix (103); après adimensionnement des nombres d'onde par $q_{\dim}$, des longueurs par $q_{\dim}^{-1}$ et des énergies par $\epsilon_{\dim}$, il nous faut donc résoudre numériquement l'équation

$$ \rho q_{\dim}^{-2} = \check{\rho}(\check{\mu} = \mu / \epsilon_{\dim}, \check{T} = k_B T / \epsilon_{\dim}) \qquad (106) $$

où le signe diacritique (un accent tchèque) marque l'adimensionnement, $\check{\rho}(\check{\mu}, \check{T})$ est la fonction équation d'état réduite et $\check{\mu}$ est la seule inconnue (les autres quantités sont imposées par (103,104)). L'équation d'état $\rho^{\text{sta}}(\mu, T)$ de la théorie BCS statique est connue analytiquement (29), ce qui donne immédiatement

$$ \check{\mu}^{\text{sta}} = \frac{1}{2} [\exp(2\eta) - 1] \approx 3,1945 \qquad (107) $$

Pour celle $\rho^{\text{dtl}}(\mu, T)$ (73) de la théorie BCS linéarisée, il faut calculer numériquement $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ et son intégrale sur $\mathbf{q}$ pour plusieurs valeurs du potentiel chimique; extrapolations linéaires et interpolations quadratiques donnent finalement

$$ \check{\mu}^{\text{dtl}} \approx 2,0545 \qquad (108) $$

<footnote>$^{41}$La notation $\eta$ reprend celle de la référence [52]; en revanche, la quantité notée $a_{2D}$ dans cette référence n'est pas notre longueur de diffusion bidimensionnelle (11) mais plutôt la taille de dimère $q_{\dim}^{-1}$.</footnote>
ce qui est fort différent et conduit dans l'expression (34), soit dit en passant, à $mc^2/\epsilon_{dim} \approx 2,5545$ donc à $k_B T/mc^2 \approx 0,130 \ll 1$ (l'effet des phonons thermiques est très faible, $q_{dim}^{-2} \delta\rho_\phi \sim [-\zeta(3)/2\pi](k_B T/mc^2)^3 mc^2/\epsilon_{dim} \approx -10^{-3}$ d'après la note 39, et les points expérimentaux sont pratiquement à température nulle). Faire le choix naïf et incorrect (107) dans le calcul grand-canonique BCS dépendant du temps linéarisé des fonctions de distribution de paires conduit aux courbes en tireté sur la moitié droite de la figure 3, le choix éclairé (108) aux courbes en trait plein. Pour la fonction $g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$, la courbe en tireté est disqualifiée comme à 3D (elle prend des valeurs négatives), mais la courbe en trait plein – notre véritable prédiction – est en assez bon accord avec l'expérience. Pour la fonction $g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$, le creux du tireté a une profondeur (comptée négativement) $\min_{|\mathbf{r}-\mathbf{r}'|} g_{\updownarrow}(\mathbf{r}, \mathbf{r}')/\rho_{\up} \rho_{\downarrow} - 1 \approx -0,0964$ plus proche de la valeur expérimentale $\approx -0,13$ que celle du trait plein $\approx -0,0592$, mais une position du minimum $k_F|\mathbf{r}-\mathbf{r}'|_0 \approx 1,53$ plus éloignée de la valeur expérimentale $\approx 2,2$ que celle du trait plein $\approx 1,738$.42,43

### 3.3. Comportement des $g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ à courte distance

#### 3.3.1. Étude de $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ pour $q \to +\infty$

Pour déterminer le comportement à courte distance d'une fonction définie par sa transformée de Fourier, il faut connaître le comportement de la transformée de Fourier aux grands nombres d'onde $q$.

Cette remarque liminaire étant faite, repartons, dans notre théorie grand-canonique BCS dépendant du temps linéarisée, de la décomposition (56) des fonctions de distribution de paires $g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$, plus précisément des écarts $\delta g_{\sigma\sigma'}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ aux valeurs asymptotiques $\rho_\sigma^{\text{dtl}} \rho_{\sigma'}^{\text{dtl}}$, en la somme des écarts $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ dans la théorie BCS statique de même potentiel chimique $\mu$ [mais pas de même densité $\rho = 2\rho_\sigma$, voir l'équation (73)] et d'une unique contribution collective $\delta g_{\text{coll}}(\mathbf{r}, \mathbf{r}')$. Dans la suite, nous travaillerons à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ et négligerons tous les nombres d'occupation de Fermi-Dirac $f_k$ dans la théorie BCS statique et dans les fonctions $\Sigma_{ij}$ (43-49) de la théorie BCS dépendant du temps.

Les écarts $\delta g_{\sigma\sigma'}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ s'expriment, grâce au théorème de Wick, en termes des moyennes anormale et normale comme dans (23,24), moyennes dont les transformées de Fourier sont connues explicitement, voir les équations (25,26), et dont le comportement à grand nombre d'onde $k$ s'obtient trivalement. Leur étude à courte distance sera simple à faire, surtout à 2D où on sait les écrire en termes de fonctions de Bessel (30,31), ce qui évite même le détour par l'espace de Fourier [1].

En revanche, la transformée de Fourier $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ de la partie collective n'a rien de simple, puisque son expression (67) utile ici (en pratique le troisième membre) est elle-même une in-

<footnote>$^{42}$Comme les données expérimentales sont assez bruitées sur $g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$, nous les remplaçons par un ajustement cubique sur l'intervalle $1,4 \le k_F|\mathbf{r}-\mathbf{r}'| \le 3,5$ afin de déterminer la position et la valeur du minimum.</footnote>

<footnote>$^{43}$Nous ne reproduisons pas ici le jeu de la grand-canonisation (102) car nous ne connaissons pas bien le potentiel chimique du gaz de fermions : il n'a pas été mesuré dans la référence [52] et les valeurs qu'on pourrait tirer de la littérature théorique ou expérimentale pour une force des interactions $\eta = 1$ c'est-à-dire $\ln(k_F a_{2D}) \approx 1,116$ sont fort dispersées : $\mu/E_F \approx 0,43$ d'après la figure 3a de [24] (approximation des fluctuations gaussiennes), $\mu/E_F \approx 0,35$ par interpolation des résultats de la table I dans [63] (Monte-Carlo quantique à surface nodale fixée), $\mu/E_F \approx 0,29$ d'après les polynômes d'interpolation de [64] (Monte-Carlo quantique à champ auxiliaire sur réseau), $\mu/E_F = 0,14(2)$ dans les données de la figure 1 de [65] (expérience d'atomes froids) [à titre de comparaison, notre équation (108) correspond à $\mu/E_F \approx 0,556$].</footnote>
tégrale de contour sur l'énergie complexe z d'une susceptibilité collective $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ assez rébarbative (55). Donnons tout de suite le résultat des développements asymptotiques de l'annexe B :

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) = \underset{q \to +\infty}{\overset{\text{en } d=3}{\left\{ \begin{array}{l} \frac{C_{\text{coll}}}{q^4} + \frac{D_{\text{coll}}}{q^5} + O\left(\frac{1}{q^6}\right) \\ \\ \frac{C_{\text{coll}}}{q^4} \frac{\ln(A_{\text{coll}}q)}{\ln(B_{\text{coll}}q)} + O\left(\frac{\ln^\alpha q}{q^6}\right) \end{array} \right.}} \quad \text{en } d=2 \quad (109) $$

avec les coefficients suivants :

$$ d=3: \quad C_{\text{coll}} = -\frac{4\pi}{3\sqrt{3}} \frac{m^2 \Delta^2}{\hbar^4} \rho^{\text{sta}} ; \quad D_{\text{coll}} = -2 \left(\frac{4\pi}{9} - \frac{2\sqrt{3}}{3}\right) \frac{m^2 \Delta^2}{\hbar^4} a_{3D}^{-1} \rho^{\text{sta}} + \left(\frac{m\Delta}{\hbar^2}\right)^4 \mathcal{D} \quad (110) $$

$$ d=2: \quad C_{\text{coll}} = -\frac{2m^2 \Delta^2}{\hbar^4} \rho^{\text{sta}} ; \quad A_{\text{coll}} = \frac{2}{\sqrt{3}q_{\text{dim}}} \neq B_{\text{coll}} = \frac{\sqrt{3}}{2q_{\text{dim}}} \quad (111) $$

Ici $\rho^{\text{sta}}$ et $\Delta$ sont la densité (21) et le paramètre d'ordre (16) dans la théorie BCS statique au potentiel chimique $\mu$, $a_{3D}$ est la longueur de diffusion tridimensionnelle (11) et $q_{\text{dim}}$ le nombre d'onde de dimère bidimensionnel (12). Nous donnons la constante $\mathcal{D}$ de l'équation (110) sous forme numérique,

$$ \mathcal{D} = -0,549686\dots \quad (112) $$

et dans l'annexe B sous forme intégrale (202); les intégrales mises en jeu pourraient être calculées analytiquement, mais au prix d'un effort qui nous a semblé excessif.

Les calculs laborieux de l'annexe B ne sont en réalité pas dénués d'intérêt, car les résultats (109) sont de prime abord assez surprenants. Cherchons en effet à les retrouver qualitativement par un argument simple de comptage de puissances $q^n$. Dans la limite $q \to +\infty$, l'énergie propre $\hbar\omega_q$ du mode de phonon (s'il existe à grand $q$, voir la note 19) et le bord $E_{\text{bord}}(\mathbf{q})$ (62) du continuum de paire brisée divergent comme $q^2$; on peut donc prendre une énergie complexe $z$ d'ordre $q^2$ dans le contour d'intégration $C_d$ de l'équation (67) [ceci autorise d'ailleurs à négliger la contribution thermique $\sim \exp(-\beta z)$ exponentiellement petite, ce qui explique pourquoi les coefficients (110,111) ne dépendent pas de la température]. On s'attend dès lors à ce que le nombre d'onde $k$ typique dans les intégrales (43-49) donnant les $\Sigma_{ij}(\mathbf{q}, z)$ soit d'ordre $q$, et les énergies $\epsilon_{\pm}$ et $\xi_{\pm}$ d'ordre $q^2$. D'où les lois de puissance lorsque $q \to +\infty$ à $Z = z/E_q$ fixé [avec $E_q = \hbar^2 q^2/2m$ comme dans (15)] :

$$ \Sigma_{11}(\mathbf{q}, z), \Sigma_{12}(\mathbf{q}, z), \Sigma_{22}(\mathbf{q}, z) \approx q^{d-2} ; \quad \Sigma_{13}(\mathbf{q}, z), \Sigma_{23}(\mathbf{q}, z) \approx q^{d-4} ; \quad \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \approx q^{d-6} \quad (113) $$

et en définitive

$$ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) \approx \underset{q \to +\infty}{\overset{\text{en } d=2}{\int}} \frac{\mathrm{d}z}{(-2i\pi)} \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \approx q^{d-4} \quad (114) $$

puisque $dz = E_q dZ$. Voilà qui est en désaccord avec les résultats (109), fort heureusement d'ailleurs. Sinon, l'intégrale de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ aurait une divergence ultraviolette logarithmique à 2D, quadratique à 3D et l'équation d'état modifiée (73) n'aurait aucun sens!

Expliquons cet échec du comptage de puissances d'abord en dimension $d=2$. Le calcul complet de l'annexe B confirme, à l'ordre dominant, que $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ est bien $\approx q^{-4}$ (à des facteurs en $\ln q$ près$^{44}$), mais chaque contribution – considérée comme une fonction de $z$ – relève de l'un des deux cas suivants, (i) elle n'a pas de singularité (pôle ou ligne de coupure) sur le demi-axe réel positif, ou (ii) elle présente des singularités, mais seulement sur le demi-axe réel positif (en particulier, pas sur le demi-axe réel négatif). Dans le premier cas, le lacet $C_d$ entoure un domaine

<footnote>$^{44}$ Le logarithme $\ln q$ résulte du fait que $\int_{k<\eta q} d^2k/(2\pi)^2(1/\epsilon_k) = (m/\pi\hbar^2)\ln(\eta q/q_{\text{dim}}) + O(1/q^2)$ lorsque $q \to +\infty$ à $\eta \ll 1$ fixé; cette intégrale est en effet une sous-partie des $\Sigma_{ij}$, elle intervient dans leur restriction à un voisinage de $\mathbf{k}_- = \mathbf{0}$ de rayon $\ll |\mathbf{k}_+ - \mathbf{k}_-| = q$ (sur ce voisinage, on peut prendre $k_+ = q$).</footnote>
$o(|\mathbf{r} - \mathbf{r}'|^{\alpha})$ ($\alpha > 0$) près à 2D, la contribution dominante à ladite moyenne divergeant comme $|\mathbf{r} - \mathbf{r}'|^{-1}$ dans le premier cas, et comme $\ln|\mathbf{r} - \mathbf{r}'|$ dans le second. Nous trouvons en définitive que

$$ g_{\downarrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{|\mathbf{r}-\mathbf{r}'|\to 0^+}{=} \begin{cases} \frac{1}{4} [\rho^{\text{sta}}(\mu)]^2 + \left(\frac{m\Delta}{4\pi\hbar^2}\right)^2 \left[ \left(\frac{1}{|\mathbf{r}-\mathbf{r}'|} - \frac{1}{a_{3D}}\right)^2 - \frac{2m\mu}{\hbar^2} \right] + o(1) & \text{en } d=3 \\ \frac{1}{4} [\rho^{\text{sta}}(\mu)]^2 + \left(\frac{m\Delta}{2\pi\hbar^2}\right)^2 \ln^2\left(\frac{|\mathbf{r}-\mathbf{r}'|}{a_{2D}}\right) + o(1) & \text{en } d=2 \end{cases} \quad (119) $$

À 3D, il a suffi d'effectuer quelques manipulations sur la représentation de Fourier (25) de la moyenne anormale, écrite ici à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ comme nous l'avons dit,

$$ \begin{aligned} \langle \hat{\psi}_{\downarrow}(\mathbf{r}') \hat{\psi}_{\uparrow}(\mathbf{r}) \rangle_{\text{sta}} = & - \int_{\mathbb{R}^3} \frac{d^3 k}{(2\pi)^3} \frac{\Delta}{2E_k} e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')} - \int_{\mathbb{R}^3} \frac{d^3 k}{(2\pi)^3} \left(\frac{\Delta}{2\epsilon_k} - \frac{\Delta}{2E_k}\right) \\ & - \int_{\mathbb{R}^3} \frac{d^3 k}{(2\pi)^3} \left(\frac{\Delta}{2\epsilon_k} - \frac{\Delta}{2E_k}\right) \left(e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')} - 1\right) \Big|_{|\mathbf{r}-\mathbf{r}'|\to 0^+} \\ = & - \frac{m\Delta}{4\pi\hbar^2} \left(\frac{1}{|\mathbf{r}-\mathbf{r}'|} - \frac{1}{a_{3D}}\right) + \frac{m^2\mu\Delta}{4\pi\hbar^2} |\mathbf{r}-\mathbf{r}'|[1+o(1)] \end{aligned} \quad (120) $$

Dans la deuxième égalité, nous nous sommes souvenu de l'équation (16) sur le paramètre d'ordre $\Delta$ et de la transformée de Fourier inverse de $1/k^2$, $\int_{\mathbb{R}^3} [d^3 k/(2\pi)^3](1/k^2) \exp[i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')] = 1/4\pi|\mathbf{r}-\mathbf{r}'|$, puis nous avons pu appliquer le théorème de l'annexe C à la troisième intégrale du second membre car $(2\epsilon_k)^{-1} - (2E_k)^{-1} = C/k^4 + O(1/k^6)$ lorsque $k \to +\infty$, avec $C = 2m^2\mu/\hbar^4$. Une simple élévation au carré donne (119) en dimension $d=3$, ce qui est bien (qualitativement) de la forme générale attendue pour un gaz de fermions de spin 1/2 en interaction de contact, voir la relation (3a) de la table V dans la référence [66].

À 2D, il a suffi de reporter les développements connus des fonctions de Bessel $K_0(x) = -\ln(xe^\gamma/2) + O(x^2 \ln x)$ ($\gamma = 0,577...$ est la constante d'Euler-Mascheroni) et $J_0(x) = 1 + O(x^2)$ dans l'expression analytique (30) de la moyenne anormale, puis de relier $q_{\text{dim}}$ à la longueur de diffusion bidimensionnelle comme dans (12) pour obtenir (119) en dimension $d=2$, ce qui est malheureusement en désaccord (même qualitatif) avec la référence [66], voir la relation (3b) de sa table V : celle-ci affirme en toute généralité pour des interactions de contact qu'il existe une constante $\mathcal{C}$, appelée coefficient de contact sur le modèle des cas 1D [67] et 3D [68], telle que

$$ g_{\downarrow\downarrow}(\mathbf{r}, \mathbf{r}') \stackrel{|\mathbf{r}-\mathbf{r}'|\to 0^+}{=} \frac{\mathcal{C}}{(2\pi)^2} \ln^2\left(\frac{|\mathbf{r}-\mathbf{r}'|}{a_{2D}}\right) + O(r^2 \ln^2 r) \quad (121) $$

Par souci d'exhaustivité, demandons-nous dans quelle mesure nous aurions pu obtenir les lois (119) sans passer par la moyenne anormale mais en travaillant directement sur la représentation de Fourier (60) de $\delta g_{\downarrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$. Le développement à grand $q$ de l'intégrale sur le vecteur d'onde interne $\mathbf{k}$ donnant $\delta \tilde{g}_{\downarrow\downarrow}^{\text{sta}}(\mathbf{q})$ n'est pas complètement trivial : il ne suffit pas, après la mise à l'échelle évidente $\mathbf{k} = q\mathbf{K}$, de développer l'intégrande en puissances de $q^{-1}$ à $\mathbf{K}$ fixé, car ceci fait apparaître une singularité en $\mathbf{K} = \mp \hat{\mathbf{q}}/2$ ($\hat{\mathbf{q}} = \mathbf{q}/q$ est la direction de $\mathbf{q}$) suffisamment sévère à 2D pour n'être pas intégrable (on a en effet $\epsilon_{\pm} \sim E_{\mathbf{q}}(\mathbf{K} \pm \hat{\mathbf{q}}/2)^2$). Le voisinage des singularités $k_{\pm} < \eta q, \eta \ll 1$, doit en réalité être traité à part, en développant l'intégrande en puissances de $q$ à $\mathbf{k}$ fixé, même en
dimension $d=3$. Nous renvoyons le lecteur à la section 4.3.5.4 de la référence [5], et donnons ici directement le résultat :$^{47}$

$$
\delta \tilde{g}_{\uparrow\downarrow}^{\text{sta}}(\mathbf{q}) \underset{q \to +\infty}{=} \begin{cases} \frac{m^2 \Delta^2}{8 \hbar^4 q} - \frac{m^2 \Delta^2}{2 \pi \hbar^4 a_{3D} q^2} + O\left(\frac{1}{q^4}\right) & \text{en } d=3 \\ \frac{m^2 \Delta^2}{\pi \hbar^4 q^2} \ln\left(\frac{q}{q_{\text{dim}}}\right) + O\left(\frac{1}{q^4}\right) & \text{en } d=2 \end{cases} \quad (122)
$$

Une transformation de Fourier inverse redonne (119) mais seulement à des constantes additives inconnues près, celles provenant du $O(1/q^4)$ négligé.<fnref n="48" /> Il serait difficile d’aller au-delà.

### 3.3.3. Étude de $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ pour $|\mathbf{r}-\mathbf{r}'| \to 0^+$

Déterminons maintenant à l'ordre quadratique le comportement à courte distance de la fonction de distribution de paires de fermions de même spin $\uparrow$ dans notre théorie BCS dépendant du temps linéarisé. À cette fin, spécialisons la relation (74) au cas $\sigma = \sigma' = \uparrow$,

$$ g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') - g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta\tilde{g}_{\text{coll}}(\mathbf{q}) [e^{i\mathbf{q}\cdot(\mathbf{r}-\mathbf{r}')} - 1] \quad (123) $$

et écrivons la moyenne normale (26) intervenant dans l'expression de Wick (24) de $g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ sous une forme similaire (ici à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$) et après multiplication haut et bas par la quantité conjuguée $\epsilon_k + \xi_k$ pour faire apparaître un numérateur constant sous le signe intégral) :

$$ g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = - \left\{ \rho^{\text{sta}} + \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{2\epsilon_k(\epsilon_k + \xi_k)} [e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')} - 1] \right\} \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{2\epsilon_k(\epsilon_k + \xi_k)} [e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}')} - 1] \quad (124) $$

Pour développer ces expressions à l'ordre $|\mathbf{r}-\mathbf{r}'|^2$, on peut utiliser les mêmes théorèmes de l'annexe C – initialement pensés pour (123) – puisque l'intégrande de (124) décroît lui aussi comme l'inverse du nombre d'onde à la puissance quatre; comme le montrera l'équation (133), l'analogie avec (109) est complète si l'on pose en dimension quelconque

$$ C_{\text{sta}} \stackrel{\forall d}{\rightleftharpoons} - \frac{m^2 \Delta^2}{\hbar^4} \rho^{\text{sta}} \quad (125) $$

Aussi serons-nous bref.

En dimension $d=3$ nous obtenons le développement

$$ g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{d=3}{\rightleftharpoons} -\frac{1}{8\pi} (C_{\text{sta}} + C_{\text{coll}}) |\mathbf{r}-\mathbf{r}'| + \frac{1}{6} |\mathbf{r}-\mathbf{r}'|^2 \left[ \mathcal{B}_{\text{sta}} + \mathcal{B}_{\text{coll}} + \frac{D_{\text{coll}}}{2\pi^2} \ln(k_{\Delta}|\mathbf{r}-\mathbf{r}'|) + o(1) \right] \quad (126) $$

<footnote n="47">⁴⁷L'intégrale sur le complémentaire du voisinage, $k_± > ηq$ ou plutôt $K_± = |K_± \hat{q}/2| > η$ puisqu'on y a mis $\mathbf{k}$ à l'échelle, est souvent difficile à calculer, comme l'avait remarqué la référence [5]. À 3D, on tombe sur $I_{3D} = \int_0^{+\infty} dK_- \int_0^{\theta_m(K_-)} \sin\theta d\theta / (1 + 2K_- \cos\theta + K_-^2)$ où $\theta_m(K) = \arccos(\max(-1, -1/2K))$; le calcul formel donne un résultat explicite $I_{3D} = [2\pi^2 - 6\ln^2 2 + \ln 3 \ln(64/27) - 3g_2(1/4) - 6g_2(1/3)]/6$ en termes de la fonction de Bose ou dilogarithme $g_2(z)$; des relations connues sur $g_2(z)$ nous ont permis de montrer qu'on a plus simplement $I_{3D} = \pi^2/4$. À 2D, on tombe sur la partie difficile $I_{2D} = \int_{1/2}^{+\infty} (dK_- / K_-) \int_{-\theta_m(K_-)}^{+\theta_m(K_-)} d\theta / (1 + 2K_- \cos\theta + K_-^2)$; en calculant dans la limite $\eta \to 0^+$ l'intégrale auxiliaire $J = \int_{K \cdot q > 0} d^2 K / [((K + \hat{q}/2)^2 + \eta^2) [((K - \hat{q}/2)^2 + \eta^2)]]$ de deux façons différentes, en variable $K$ puis en variable $K_-$, on trouve que $J = -2\pi \ln\eta + o(1)$ puis que $J = -\pi \ln 3 + I_{2D} - 2\pi \ln\eta + o(1)$, d'où $I_{2D} = \pi \ln 3$.</footnote>

<footnote n="48">⁴⁸À 3D, on régularise la transformée de Fourier inverse (TFI) de $1/q$ par l'habituel facteur $\exp(-\eta q)$ où $\eta \to 0^+$, pour obtenir TFI(1/q) = 1/2π²r². À 2D, on introduit comme intermédiaire la fonction d'onde du dimère $\phi_{\text{dim}}(\mathbf{r})$, dont la transformée de Fourier (TF) est connue, $\tilde{\phi}_{\text{dim}}(\mathbf{k}) = 2\pi^{1/2} q_{\text{dim}}/(k^2 + q_{\text{dim}}^2)$. Ceci permet de calculer analytiquement la TF de la densité de probabilité du dimère, $\tilde{\rho}_{\text{dim}}(\mathbf{q}) = \text{TF}(q_{\text{dim}}^2(\mathbf{r})) = \ln[1 + Z + Y(Y + 4 + Z)/2]/Z$ avec $Z = \sqrt{Y(Y + 4)}$ et $Y = q^2/q_{\text{dim}}^2$; on vérifie alors que $\tilde{\rho}_{\text{dim}}(\mathbf{q}) = 2\ln Y/Y + O(\ln Y/Y^2) = 4\ln(q/q_{\text{dim}})/(q^2/q_{\text{dim}}^2) + O(1/q^4)$ lorsque $q \to +\infty$. Or, $\phi_{\text{dim}}(\mathbf{r}) = (q_{\text{dim}}/\pi^{1/2}) K_0(q_{\text{dim}} r) = (q_{\text{dim}}/\pi^{1/2}) [-\ln(r/a_{2D})] + O(r^2 \ln r)$ lorsque $r \to 0^+$, ce qui permet de conclure.</footnote>
avec le nombre d'onde $k_Δ = (mΔ/ħ^2)^{1/2}$ et les coefficients suivants, provenant de la théorie BCS statique (indice « sta ») ou de la partie collective de la théorie BCS dépendant du temps linéarisée (indice « coll ») :

$$ \mathcal{B}_{\text{sta}} \stackrel{d=3}{=} \frac{3k_0 C_{\text{sta}}}{8\pi} - 6 \left( \frac{C_{\text{sta}}}{8\pi \rho^{\text{sta}}} \right)^2 + \int_{\mathbb{R}^3} \frac{d^3 k}{(2\pi)^3} k^2 \left[ \frac{1}{2} \left( 1 - \frac{\xi \mathbf{k}}{\epsilon \mathbf{k}} \right) \rho^{\text{sta}} + \frac{C_{\text{sta}}}{(k^2 + k_0^2)^2} \right] \quad (127) $$

$$ \mathcal{B}_{\text{coll}} \stackrel{d=3}{=} \frac{3q_0 C_{\text{coll}}}{8\pi} + \frac{D_{\text{coll}}}{2\pi^2} \left( -\frac{13}{12} + \gamma + \ln \frac{q_0}{k_{\Delta}} \right) - \int_{\mathbb{R}^3} \frac{d^3 q}{(2\pi)^3} q^2 \left[ \delta \tilde{g}_{\text{coll}}(\mathbf{q}) - \frac{C_{\text{coll}}}{(q^2 + q_0^2)^2} - \frac{q D_{\text{coll}}}{(q^2 + q_0^2)^3} \right] \quad (128) $$

Rappelons que $\gamma = 0,577215...$ est la constante d'Euler-Mascheroni. Par construction, les quantités (127,128) ne dépendent pas des nombres d'onde arbitraires $k_0$ ou $q_0 > 0$ (on peut s'en convaincre en vérifiant que leurs dérivées par rapport à $k_0$ ou $q_0$ sont identiquement nulles). À la limite unitaire, nous trouvons numériquement que $\mathcal{B}_{\text{sta}} \approx 5,07 \cdot 10^{-3} k_{\Delta}^8$ et que $\mathcal{B}_{\text{coll}} \approx -1,29 \cdot 10^{-2} k_{\Delta}^8$, ce qui est utile pour le panneau inférieur gauche de la figure 3 - la loi (126) y est représentée en rouge. Il mérite d'être souligné que le départ linéaire de $g_{\|\|}(\mathbf{r}, \mathbf{r}')$ n'est pas propre à la théorie BCS dépendant du temps linéarisé; il préexiste dans la théorie BCS statique, il est vrai avec une pente significativement plus faible puisque $C_{\text{coll}}/C_{\text{sta}} = 4\pi/3\sqrt{3} \approx 2,42$ indépendamment de la force des interactions. Le départ linéaire est donc significativement renforcé mais pas créé par les excitations collectives. En revanche, la contribution en $\ln|\mathbf{r}-\mathbf{r}'|$ dans le coefficient du terme quadratique est une nouveauté, provenant du terme en $1/q^5$ (109) dans $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$ à grand $q$, absent de $\delta\tilde{g}_{\|\|}(\mathbf{q})$ comme nous le verrons dans l'équation (133).

En dimension $d=2$ nous obtenons le développement

$$ g_{\|\|}^{\text{dtl}}(\mathbf{r}, \mathbf{r}') \stackrel{d=2}{=} \frac{1}{4} |\mathbf{r}-\mathbf{r}'|^2 \left\{ \begin{aligned} & \mathcal{B}_{\text{sta}} + \frac{C_{\text{sta}}}{2\pi} \left[ \ln(q_{\text{dim}}|\mathbf{r}-\mathbf{r}'|/2) + \gamma - 1 \right] \\ & + \mathcal{B}_{\text{coll}} + \frac{C_{\text{coll}}}{2\pi} \left[ \ln(q_{\text{dim}}|\mathbf{r}-\mathbf{r}'|/2) - \ln\left(\frac{A_{\text{coll}}}{B_{\text{coll}}}\right) \ln(-\ln(q_{\text{dim}}|\mathbf{r}-\mathbf{r}'|)) + \gamma - 1 \right] + o(1) \end{aligned} \right\} \quad (129) $$

Comme précédemment, l'indice des coefficients indique leur provenance. Les deux premières contributions entre accolades sont donc celles de la théorie BCS statique; le reste provient de la partie collective de BCS dépendant du temps linéarisé. Il nous a fallu introduire les quantités

$$ \mathcal{B}_{\text{sta}} \stackrel{d=2}{=} \frac{C_{\text{sta}}}{2\pi} \ln\left(\frac{k_0}{q_{\text{dim}}}\right) + \int_{\mathbb{R}^2} \frac{d^2 k}{(2\pi)^2} \left[ \frac{k^2}{2} \left(1 - \frac{\xi \mathbf{k}}{\epsilon \mathbf{k}}\right) \rho^{\text{sta}} + \frac{C_{\text{sta}}}{k^2 + k_0^2} \right] = \frac{C_{\text{sta}}}{2\pi} \left( \frac{1}{2} - \frac{\pi \rho^{\text{sta}}}{2 q_{\text{dim}}^2} \right) \quad (130) $$

$$ \mathcal{B}_{\text{coll}} \stackrel{d=2}{=} \frac{C_{\text{coll}}}{2\pi} \left[ \ln\left(\frac{q_0}{q_{\text{dim}}}\right) + \ln\left(\frac{A_{\text{coll}}}{B_{\text{coll}}}\right) \ln\left(\ln(B_{\text{coll}}q_0)\right) \right. \\ \left. - \int_{\mathbb{R}^2} \frac{d^2 q}{(2\pi)^2} \left[ q^2 \delta \tilde{g}_{\text{coll}}(\mathbf{q}) - \frac{C_{\text{coll}}}{q^2 + q_0^2} \frac{\ln[A_{\text{coll}}^2(q^2 + q_0^2)]}{\ln[B_{\text{coll}}^2(q^2 + q_0^2)]} \right] \right] \quad (131) $$

Indépendantes du nombre d'onde $k_0 > 0$ ou du nombre d'onde $q_0$ pourvu qu'on ait $q_0 > 1/B$, ces quantités sont les équivalents bidimensionnels de (127,128), d'où la notation à l'identique. Nous avons donné l'expression explicite très simple de $\mathcal{B}_{\text{sta}}$ (à 3D, elle ferait apparaître des intégrales elliptiques).49 Pour une force des interactions $\ln(k_F^{\text{dtl}}/q_{\text{dim}}) = 1$, nous trouvons numériquement que $\mathcal{B}_{\text{coll}} \approx -0,291 q_{\text{dim}}^6$, ce qui nous a permis de représenter (129) en rouge sur le panneau inférieur droit de la figure 3. Il y avait déjà dans BCS statique une contribution logarithmique au

<footnote n="49">49 L'application directe des résultats de l'annexe C à la fonction $\tilde{g}_{\|\|}^{\text{sta}}(\mathbf{q})$ de l'équation (61) écrite à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ conduirait compte tenu de (133) à $\mathcal{B}_{\text{sta}} = (C_{\text{sta}}/2\pi) \ln(q_0/q_{\text{dim}}) - \int_{\mathbb{R}^2} [d^2 q/(2\pi)^2] [q^2 \tilde{g}_{\|\|}^{\text{sta}}(\mathbf{q}) - C_{\text{sta}}/(q^2 + k_0^2)]$, expression indépendante de $q_0 > 0$ et en principe équivalente, mais nettement plus difficile à calculer; la même remarque pourrait être faite à 3D.</footnote>
coefficient du terme quadratique; elle est quand même triplée dans BCS dépendant du temps linéarisée, puisque $C_{\text{coll}} = 2C_{\text{sta}}$. Au contraire, la contribution en double logarithme est une nouveauté, due aux facteurs logarithmiques dans la queue en $1/q^4$ (109) de $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$, absents de $\delta\tilde{g}_{\updownarrow\uparrow}(\mathbf{q})$ comme nous le verrons dans l'équation (133).

Par souci d'exhaustivité et de symétrie entre $g_{\updownarrow\uparrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ et $g_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$, voyons dans quelle mesure nous pouvons développer $g_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ à courte distance sans passer par la moyenne normale (26) mais en appliquant directement les théorèmes de l'annexe C à sa représentation intégrale, qui tient compte automatiquement de sa limite nulle à distance nulle :

$$g_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}') = \int_{\mathbb{R}^d} \frac{d^d q}{(2\pi)^d} \delta\tilde{g}_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{q}) [e^{i\mathbf{q} \cdot (\mathbf{r} - \mathbf{r}')} - 1] \quad (132)$$

où la transformée de Fourier $\delta\tilde{g}_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{q})$ est donnée par l'équation (61). Par la même méthode que celle ayant conduit à (122), nous trouvons indépendamment de la dimension que

$$\delta \tilde{g}_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{q}) \stackrel{\frac{\gamma d}{q-+\infty}}{=} \frac{C_{\text{sta}}}{q^4} + O\left(\frac{1}{q^6}\right) \quad (133)$$

L'introduction du coefficient $C_{\text{sta}}$ avait été anticipée par l'équation (125). Pour le reste, la théorie BCS statique est plus pauvre aux grands nombres d'onde que la théorie BCS dépendant du temps linéarisée au sens où, dans les notations des lois asymptotiques (109) sur $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$, la loi (133) correspond en dimension $d=3$ à un coefficient $D_{\text{sta}}=0$ (pas de contribution sous-dominante en $1/q^5$) et en dimension $d=2$ à des coefficients égaux $A_{\text{sta}} = B_{\text{sta}}$ (pas de facteurs logarithmiques dans le terme dominant $1/q^4$).50 L'application des théorèmes redonne en tout cas directement les contributions comprenant un facteur $C_{\text{sta}}$ dans (126,129) mais conduit à des expressions des $\mathcal{B}_{\text{sta}}$ peu utilisables en pratique car difficiles à intégrer, voir la note 49.

## 4. Conclusion générale

Nous avons étudié à basse température les fonctions de distribution de paires $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ d'un gaz bidimensionnel ou tridimensionnel de fermions de spin 1/2 en interaction de contact, non polarisé et spatialement homogène, dans ce qui nous paraît être la théorie de champ moyen minimale acceptable pour leur calcul, la théorie BCS dépendant du temps linéarisée, dont les fonctions de réponse densité-densité $\chi_{\sigma\sigma'}(\mathbf{r}, t; \mathbf{r}', t)$ déduites d'une dynamique de champ classique conduisent à des prédictions entièrement quantiques sur les $g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$ par la grâce du théorème de fluctuation-dissipation. En revanche, la théorie BCS statique, qui ignore totalement les excitations collectives, en particulier la branche acoustique – et donc les phonons – de Bogolioubov-Anderson, est disqualifiée. L'étude est menée dans l'ensemble grand-canonique de potentiel chimique $\mu$ et de température $T$, et la densité du gaz à l'équilibre s'en déduit par l'équation d'état $\rho(\mu, T)$. L'insuffisance la plus criante de BCS statique apparaît dans la fonction de distribution de paires de fermions de spins opposés $g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$ et son écart $\delta g_{\updownarrow}(\mathbf{r}, \mathbf{r}')$ à sa valeur

<footnote n="50">50 La situation est plus simple que celle de $\delta\tilde{g}_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{q})$ (voir la note 47 et le texte qui l'appelle) car la fonction $1 - \xi_k / \epsilon_k = \Delta^2 / \epsilon_k (\epsilon_k + \xi_k)$ est d'intégrale convergente sur $\mathbb{R}^d$ même pour $d=3$; l'intégrale sur $\mathbf{k}$ définissant $\delta\tilde{g}_{\updownarrow\uparrow}^{\text{sta}}(\mathbf{q})$ dans (61) est alors dominée aux grands $q$ par des voisinages de $\mathbf{k}_{\pm} = 0$ de rayon $\eta q$, $\eta \ll 1$ fixé : on peut remplacer $\mathbf{k}_{\mp}$ par $\mathbf{q}$ dans le reste de l'intégrande et on sort une loi en $q^{-4}$. Par comptage de puissances, on vérifie en effet que le domaine d'intégration $k_{\pm} \approx q$ donne une contribution sous-dominante en $O(q^{d-8})$ (on fait le changement de variable $\mathbf{k} = q\mathbf{K}$ puis on passe à la limite $q \to +\infty$ à $\mathbf{K}$ fixé). Cet argument suffit à montrer que $D_{\text{sta}} = 0$ à 2D, mais pas à 3D ($3-8=-5$). Pour $d=3$, il faut faire le calcul complet, l'annulation de $D_{\text{sta}}$ résultant d'une compensation subtile entre la zone de proximité $k_+$ ou $k_- < \eta q$ et la zone asymptotique complémentaire $k_+$ et $k_- > \eta q$. Mathématiquement, ceci provient du fait a priori non évident qu'une certaine expression intégrale $I(\eta)$ est de limite nulle lorsque le paramètre de coupure $\eta < 1/2$ séparant les deux zones tend vers $0^+$, à savoir $I(\eta) = -2\eta^{-1} + \int_{\eta}^{+\infty} (dK_- / K^2) \int_0^{\theta_m(K^-)} (\sin\theta d\theta) / (K^2 + 2K_- \cos\theta + 1)^2 = -\eta / (1 - \eta^2) - 3 \arg\eta$ où l'angle $\theta_m(K)$ est défini dans la note 47 (le terme $-2\eta^{-1}$ provient de la zone de proximité, le reste de la zone asymptotique).</footnote>
asymptotique $|\rho(\mu, T)/2|^2$ : cette théorie décrit en fait seulement la structure spatiale d'une paire liée $\uparrow\downarrow$ du condensat de paires, puisque $\delta g_{\uparrow\downarrow}^{\text{sta}}(\mathbf{r}, \mathbf{r}')$ est proportionnelle à $|\phi_0(\mathbf{r}, \mathbf{r}')|^2$ avec $\phi_0(\mathbf{r}, \mathbf{r}')$ la fonction d'onde d'une paire liée, voir la section 1; ainsi $\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ serait une fonction partout positive, tendant de façon monotone vers zéro, avec une constante de décroissance exponentielle $\kappa_0$ de limite non nulle à température nulle, voir la note 2. La théorie BCS dépendant du temps linéarisée montre que, sous l'effet des fluctuations quantiques ou thermiques dans les modes de phonons, la réalité est bien différente : à basse température,

(i) La fonction $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ admet un minimum absolu strictement inférieur à sa valeur asymptotique, un creux donc, avec une profondeur assez proche de celle du Monte-Carlo quantique de la référence [50] pour le gaz unitaire à 3D, et de celle dans l'expérience de la référence [52] à 2D pour une force des interactions $\ln(k_F/q_{\text{dim}}) = 1$, si l'on utilise l'équation d'état BCS statique $\rho^{\text{sta}}(\mu, T)$ pour passer de $\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ à $g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$; l'accord se dégrade cependant, de façon catastrophique à 3D, si l'on utilise – comme il se doit en principe – l'équation d'état modifiée sous l'effet des modes collectifs d'excitation $\rho^{\text{dtl}}(\mu, T)$, voir la section 3.2; ceci résulte d'une absence d'autocohérence de la théorie BCS dépendant du temps linéarisée.

(ii) La fonction $\delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ tend vers zéro à l'infini par valeurs strictement négatives, avec une constante $\kappa_0 \sim 2\pi k_B T / \hbar c$ où $c$ est la vitesse du son; à température nulle $\delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ tend vers zéro comme l'inverse de la distance $|\mathbf{r}-\mathbf{r}'|$ à la puissance $d+1$ où $d \in \{2;3\}$ est la dimension de l'espace, là aussi par valeurs strictement négatives.

Dans la limite $T \to 0^+$ et $|\mathbf{r}-\mathbf{r}'| \to +\infty$, la section 3.1 donne plus précisément un approximant original de $\delta g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ sous forme d'une série partout négative raccordant le régime asymptotique quantique $|\mathbf{r}-\mathbf{r}'| \ll \hbar c / (2\pi k_B T)$, où le système semble être dans son état fondamental, au régime asymptotique thermique $|\mathbf{r}-\mathbf{r}'| \gg \hbar c / (2\pi k_B T)$. À une renormalisation près de la densité et de la vitesse du son, l'approximant est en accord avec l'hydrodynamique quantique de Landau et Khatnikov, une théorie effective de basse énergie que l'on pense être exacte à l'ordre dominant en température. Il vaut aussi pour $\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$, qui est équivalent à $\delta g_{\uparrow\downarrow}(\mathbf{r}, \mathbf{r}')$ aux distances beaucoup plus grandes que la taille d'une paire liée. En revanche, à courte distance $|\mathbf{r}-\mathbf{r}'| \to 0$, le comportement – divergeant bien entendu vers $+\infty$ – de $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ reste le même que celui de la théorie BCS statique de même potentiel chimique, l'écart n'étant qu'un $o(1)$, voir la section 3.3.

En ce qui concerne la fonction $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ prise à courte distance, l'écart entre BCS statique et BCS dépendant du temps linéarisé de même potentiel chimique est seulement quantitatif à l'ordre dominant en $|\mathbf{r}-\mathbf{r}'|$: à 3D, les coefficients du départ linéaire $|\mathbf{r}-\mathbf{r}'|$ diffèrent par un facteur $1+4\pi/3\sqrt{3} \approx 3,4$; à 2D, ceux du départ logarithmico-quadratique $|\mathbf{r}-\mathbf{r}'|^2 \ln|\mathbf{r}-\mathbf{r}'|$ diffèrent par un facteur 3. Les coefficients les plus faibles sont ceux de BCS statique, et les rapports donnés valent indépendamment de la force des interactions. Les différences qualitatives apparaissent à l'ordre sous-dominant, la théorie BCS dépendant du temps linéarisée comportant des termes en $|\mathbf{r}-\mathbf{r}'|^2 \ln|\mathbf{r}-\mathbf{r}'|$ à 3D et en $|\mathbf{r}-\mathbf{r}'|^2 \ln(-\ln|\mathbf{r}-\mathbf{r}'|)$ à 2D, totalement absents de BCS statique, voir la section 3.3. Là aussi, il faut utiliser l'équation d'état modifiée dans le calcul de $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$, le théorème de fluctuation-dissipation donnant seulement les écarts $\delta g_{\sigma\sigma'}(\mathbf{r}, \mathbf{r}')$. Le changement d'équation d'état de $\rho^{\text{sta}}(\mu, T)$ en $\rho^{\text{dtl}}(\mu, T)$ est cette fois bienvenu : au lieu de dégrader l'accord avec la théorie ou l'expérience, il permet d'éviter une absurdité, celle de valeurs de $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ strictement négatives sur un voisinage de la distance nulle! Il ne règle cependant pas un net désaccord avec le Monte-Carlo quantique sur le gaz unitaire à 3D, selon lequel le départ de $g_{\uparrow\downarrow}^{\text{dtl}}(\mathbf{r}, \mathbf{r}')$ est quadratique plutôt que linéaire, voir la section 3.2.

Le présent article n'est qu'une première incursion dans le sujet des fonctions de distribution de paires des gaz de fermions. Il pose les bases de la théorie BCS minimale acceptable mais laisse en suspens deux questions importantes, (i) une comparaison à une théorie plus élaborée à
base d'intégrale de chemin sur un champ bosonique auxiliaire (après sommation sur les degrés de liberté fermioniques), prise dans l'approximation dite des fluctuations gaussiennes [25, 32], (ii) un calcul explicite à tous les ordres en le facteur de dissociation thermique des paires liées $\exp(-\beta E_{gap})$, en particulier au voisinage de la transition superfluide, dont on peut se demander si elle laisse une empreinte observable dans les fonctions de distribution de paires, par exemple au travers d'une atteinte exponentielle de leurs valeurs asymptotiques de plus en plus lente lorsque $T \to T_c$, comme c'est le cas pour la condensation de Bose-Einstein du gaz parfait.<fnref n="51" /> L'article nous a semblé assez long pour que la résolution de ces questions soit remise à plus tard.

Annexe A. Développement à faible $q$ de $\delta\tilde{g}_{coll}(q)$ à l'ordre zéro en $\exp(-\beta E_{gap})$

À l'ordre zéro en $\exp(-\beta E_{gap})$ (les nombres d'occupation de Fermi-Dirac étant mis à zéro), l'équation (68) sépare la transformée de Fourier $\delta\tilde{g}_{coll}(q)$ de la partie collective des fonctions de distribution de paires en la contribution $\Phi(q)(\tilde{n}_q + 1/2)$ du mode de phonon et la contribution $\delta\tilde{g}_{coll}^{C0}(q)$ du continuum de paire brisée, en décomposant le chemin d'intégration sur l'énergie complexe $z$ en un cercle $c_d$ entourant le pôle de phonon $z = \hbar\omega_q$ et un lacet $L_d$ entourant la ligne de coupure $z \in [E_{bord}(q), +\infty[$ du continuum. Ici, nous déterminons le comportement des fonctions $\Phi(q)$ et $\delta\tilde{g}_{coll}^{C0}(q)$ à faible nombre d'onde $q$; rappelons que, dans cette limite, l'existence du mode de phonon est garantie, voir la note 19.

A.1. *Ce que donne le continuum de paire brisée*

Commençons par la contribution du continuum de paire brisée et partons de l'expression (55) de la susceptibilité collective $\tilde{\chi}_{coll}(q, z)$. Comme le lacet $L_d$ reste à distance non infinitésimale des singularités sur l'axe réel, nous pouvons développer les fonctions $\Sigma_{ij}(q, z)$ en puissances de $q$ à $z$ fixé, ce qui ne fait intervenir que des puissances paires (les intégrandes sont des fonctions lisses et paires du vecteur d'onde $q$) :

$$ \Sigma_{ij}(q, z) \stackrel{z \text{ fixé}}{=} \sum_{ij}^{(0)}(z) + q^2 \sum_{ij}^{(2)}(z) + O(q^4) \quad (134) $$

On remarque que les ordres zéro ne font intervenir que deux fonctions indépendantes :

$$ \sigma_0(z) = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{1}{\epsilon_k(z^2 - 4\epsilon_k^2)} ; \quad \sigma_1(z) = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\xi_k}{\epsilon_k(z^2 - 4\epsilon_k^2)} \quad (135) $$

On a en effet

$$ \Sigma_{11}^{(0)}(z) = \frac{1}{2} z^2 \sigma_0(z) \quad ; \quad \Sigma_{22}^{(0)}(z) = \frac{1}{2} (z^2 - 4\Delta^2) \sigma_0(z) \quad ; \quad \Sigma_{12}^{(0)}(z) = z \sigma_1(z) \quad (136) $$

$$ \Sigma_{13}^{(0)}(z) = z\Delta\sigma_0(z) \quad ; \quad \Sigma_{23}^{(0)}(z) = 2\Delta\sigma_1(z) \quad (137) $$

Ceci conduit à une belle simplification dans le développement de $\tilde{\chi}_{coll}(q, z)$ :

$$ \tilde{\chi}_{\text{coll}}(\mathbf{q}, z) \stackrel{\text{z fixé}}{=} -\Delta^2 \sigma_0(z) - \frac{q^2}{z^2} \left[ (z^2 - 2\Delta^2) \Sigma_{11}^{(2)}(z) - z^2 \Sigma_{22}^{(2)}(z) \right] + O(q^4) \quad (138) $$

Comme l'ordre zéro n'est autre que $-(1/2)[\Sigma_{11}^{(0)}(z) - \Sigma_{22}^{(0)}(z)]$, on en déduit en particulier que $\tilde{\chi}_{11}^{\text{dtl}}(q, z)$ a un départ quadratique $\approx q^2$, si l'on se souvient des équations (41,51). Il reste à intégrer (138) sur $L_d$ comme dans (68). Pour cela, on contracte le lacet sur l'axe réel et on utilise les relations $\Sigma_{ij}(q, z = \epsilon + i0^+) - \Sigma_{ij}(q, z = \epsilon - i0^+) = 2i \text{Im}[\Sigma_{ij}(q, z = \epsilon + i0^+)]$, $\epsilon \in \mathbb{R}$ et $\text{Im}[1/(x + i0^+)] =$

<footnote n="51">⁵¹Mathématiquement, nous n'excluons pas que, dans le cas des fermions, cette relaxation exponentielle soit oscillante, c'est-à-dire de la forme $\exp(ik_0|r-r'|)$ plus complexe conjugué avec une partie réelle Re $k_0$ non infinitésimale mais une partie imaginaire $\kappa_0 = \text{Im} k_0 > 0$ tendant vers zéro à la transition</footnote>
$-\pi\delta(x)$, $x \in \mathbb{R}$ (après décomposition de l'intégrande de $\Sigma_{ij}(\mathbf{q}, z)$ en éléments simples vis-à-vis de la variable $z$) pour calculer l'intégrale sur $\epsilon \in [E_{\text{bord}}(q), +\infty[$.$^{52}$ Dans les intégrales sur le vecteur d'onde $\mathbf{k}$ restantes, on développe en puissances de $\mathbf{q}$ sous le signe somme puis on intègre terme à terme pour obtenir les intégrales curvilignes des $\Sigma_{ij}^{(n)}(z)$ utiles à l'équation (138) :

$$ \int_{L_d} \frac{dz}{(-2i\pi)} [\Sigma_{11}(\mathbf{q}, z) - \Sigma_{22}(\mathbf{q}, z)] = +\Delta^2 E_q \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{2u^2 E_k (\xi_k^2 - \Delta^2) - \xi_k^2 \xi_k}{4\epsilon_k^6} + O(q^4) \quad (139) $$

$$ \begin{aligned} \int_{L_d} \frac{dz}{(-2i\pi)} \frac{\Sigma_{11}(\mathbf{q}, z)}{z^2} &= \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\epsilon_+ \epsilon_- + \xi_+ \xi_- + \Delta^2}{4\epsilon_+ \epsilon_- (\epsilon_+ + \epsilon_-)^2} \\ &= \int_{q \to 0^+} \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{1}{8\epsilon_k^2} - E_q \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\epsilon_k^2 \xi_k + 4\Delta^2 E_k u^2}{16\epsilon_k^6} + O(q^4) \end{aligned} \quad (140) $$

où $u$ est le cosinus de l'angle entre $\mathbf{q}$ et $\mathbf{k}$ et l'on rappelle que $E_{\mathbf{k}} = \hbar^2 k^2 / 2m$. Après regroupement dans (138) il vient

$$ \delta \tilde{g}_{\text{coll}}^{C0}(\mathbf{q}) = \int_{q \to 0^+} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{4\epsilon_k^2} + E_q \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2 \xi_k (\epsilon_k^2 - 4u^2 E_k \xi_k)}{8\epsilon_k^6} + O(q^4) \quad (141) $$

ce qui répond à la question. Nous ne pouvons nous empêcher de remarquer que le résultat est beaucoup plus simple pour la distribution de paires de fermions de spins opposés, où l'on a

$$ \delta \tilde{g}_{\uparrow\uparrow}^{C0}(\mathbf{q}) = -E_q \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2 u^2 E_k}{4\epsilon_k^4} + O(q^4) \quad (142) $$

## A.2. *Ce que donne le phonon de Bogolioubov-Anderson*

Passons maintenant au développement de la fonction $\Phi(\mathbf{q})$ dans la contribution du phonon de Bogolioubov-Anderson à $\delta\tilde{g}_{\text{coll}}(\mathbf{q})$, en partant de son expression (69). On réutilise les développements (134) des $\Sigma_{ij}(\mathbf{q}, z)$ en puissances de $q$ en y ajoutant un développement en puissances de $z = \hbar\omega_{\mathbf{q}} = \hbar cq + O(q^3)$ (c est la vitesse du son), ce que l'on peut faire puisque $z \ll E_{\text{bord}}(\mathbf{q})$ est très éloigné de toute singularité des $\Sigma_{ij}$. En se limitant aux ordres utiles ici, on obtient :

$$ \Sigma_{11}(\mathbf{q}, z) \stackrel{z=O(q)}{=} \frac{1}{2} z^2 \sigma_0(0) + q^2 \Sigma_{11}^{(2)}(0) + O(q^4) \quad ; \quad \Sigma_{22}(\mathbf{q}, z) \stackrel{z=O(q)}{=} -2\Delta^2 \sigma_0(0) + O(q^2) \quad (143) $$

$$ \Sigma_{12}(\mathbf{q}, z) \stackrel{z=O(q)}{=} \frac{z\sigma_1(0) + O(q^3)}{q \to 0^+} \quad ; \quad \Sigma_{13}(\mathbf{q}, z) \stackrel{z=O(q)}{=} \frac{z\Delta\sigma_0(0) + O(q^3)}{q \to 0^+} \quad (144) $$

$$ \Sigma_{23}(\mathbf{q}, z) \stackrel{z=O(q)}{=} \frac{2\Delta\sigma_1(0) + O(q^2)}{q \to 0^+} \quad (145) $$

L'équation (63) sur l'énergie propre du mode de phonon se réduit alors à

$$ 0 = \left[ -|\Delta\sigma_0(0) + i\sigma_1(0)|^2 z^2 - 2\Delta^2\sigma_0(0)\Sigma_{11}^{(2)}(0)q^2 \right] \Big|_{z=\hbar cq} + O(q^4) \quad (146) $$

ce qui permet d'éliminer le coefficient $\Sigma_{11}^{(2)}(0)$ dans la suite; on en déduit en particulier que

$$ \Sigma_{11}(\mathbf{q}, z = \hbar\omega_{\mathbf{q}}) = \frac{(\hbar c)^2 \sigma_1^2(0)}{2\Delta^2 |\sigma_0(0)|} q^2 + O(q^4) \quad (147) $$

et que son signe $s_{11}$ est égal à 1. On a pu remplacer la quantité $\sigma_0(0)$ par $-|\sigma_0(0)|$ puisqu'elle est négative, comme on le voit bien sur sa définition (135); en revanche, le signe de la quantité $\sigma_1(0)$ donc de $\Sigma_{12}(\mathbf{q}, \hbar\omega_{\mathbf{q}})$ est indéterminé, mais c'est aussi celui de $\Sigma_{23}(\mathbf{q}, \hbar\omega_{\mathbf{q}})$ donc il se factorise et

<footnote>$^{52}$En pratique, il suffit de savoir que $\text{Im}\{((\epsilon_+ + \epsilon_-)/[(\epsilon + i0^+)^2 - (\epsilon_+ + \epsilon_-)^2])\} = (\pi/2)[\delta(\epsilon + \epsilon_+ + \epsilon_-) - \delta(\epsilon - \epsilon_+ - \epsilon_-)]$.</footnote>
disparaît dans le carré au numérateur de $\Phi(\mathbf{q})$ dans (69). Le dénominateur de $\Phi(\mathbf{q})$ dans (69) s'obtient par simple dérivation par rapport à $z$ de l'expression entre crochets dans (146). On aboutit en définitive à l'équivalent de $\Phi(\mathbf{q})$ cherché :

$$ \Phi(\mathbf{q}) = \frac{\hbar c q}{|\sigma_0(0)| |\Delta\sigma_0(0) + i\sigma_1(0)|^2 + O(q^3)} \quad (148) $$

avec, rappelons-le pour la commodité de la lecture,

$$ \sigma_0(0) = - \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{1}{4\epsilon_k^3} ; \quad \sigma_1(0) = - \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\xi_k}{4\epsilon_k^3} \quad (149) $$

Pour obtenir la forme (83) beaucoup plus belle de $\Phi'(0)$ utilisée dans la section 3.1, il suffit de prendre la dérivée par rapport au potentiel chimique de l'équation (16) sur le paramètre d'ordre et de l'équation (21) sur la densité à longueur de diffusion $a_{dD}$ fixée, puis de se souvenir de la relation hydrodynamique (32) satisfaite par la théorie BCS statique :

$$ \Delta \frac{d}{d\mu} \Delta = \frac{\sigma_1(0)}{\sigma_0(0)} ; \quad \frac{\rho^{\text{sta}}}{mc^2} = \frac{d}{d\mu} \rho^{\text{sta}} = -4 \left[ \Delta^2 \sigma_0(0) + \Delta \left( \frac{d}{d\mu} \Delta \right) \sigma_1(0) \right] = \frac{4}{|\sigma_0(0)|} |\Delta\sigma_0(0) + i\sigma_1(0)|^2 \quad (150) $$

L'équation (148) se simplifie alors comme annoncé :

$$ \Phi(\mathbf{q}) = \frac{\hbar \rho^{\text{sta}}}{4mc} q + O(q^3) \quad (151) $$

### A.3. Applications

En première application, nous vérifions que, dans notre théorie BCS dépendant du temps linéarisée, le phonon de Bogoliubov-Anderson sature la relation sommatoire bien connue sur la susceptibilité densité-densité [69]

$$ \mathcal{R}(\mathbf{q}) = \int_{\mathbb{R}} \frac{d\epsilon}{(-\pi)} \epsilon \operatorname{Im}[\tilde{\chi}_{\rho\rho}(\mathbf{q}, \epsilon + i0^+)] = \rho \frac{\hbar^2 q^2}{m} \quad (152) $$

à l'ordre dominant en $q$ (la relation ne serait catastrophiquement pas vérifiée par la théorie BCS statique à faible $q$ [11]). De la note 12 puis de l'équation (51) nous déduisons que $\tilde{\chi}_{\rho\rho}(\mathbf{q}, z) = 2\Sigma_{33}(\mathbf{q}, z) + 4\tilde{\chi}_{\text{coll}}(\mathbf{q}, z) = 2(\Sigma_{11} - \Sigma_{22} + \Sigma_{44})(\mathbf{q}, z) + 4\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$. La fonction $z \mapsto \tilde{\chi}_{\rho\rho}(\mathbf{q}, z)$ admet donc un pôle en $z = \hbar\omega_q$ de résidu $2\Phi(\mathbf{q})$ et un pôle en $z = -\hbar\omega_q$ de résidu opposé (c'est une fonction paire), si l'on revient à la définition de $\Phi(\mathbf{q})$ dans le texte précédant l'équation (69). Comme $\operatorname{Im}[1/(\epsilon \mp \hbar\omega_q + i0^+)] = -\pi\delta(\epsilon \mp \hbar\omega_q)$, la contribution du phonon à la relation sommaire s'écrit

$$ \mathcal{R}^{\phi}(\mathbf{q}) = 4\Phi(\mathbf{q})\hbar\omega_q = \rho^{\text{sta}} \frac{\hbar^2 q^2}{m} + O(q^4) \quad (153) $$

Il reste donc à montrer que la contribution du continuum est un $O(q^4)$. À l'aide du développement (138) et de l'expression (46) de $\Sigma_{44}(\mathbf{q}, z)$, nous trouvons que $\Sigma_{44}^{(0)}(z) = 0$ et que

$$ \frac{1}{4} \tilde{\chi}_{\rho\rho}(\mathbf{q}, z) \stackrel{\text{z fixé}}{=} \frac{2q^2 \Delta^2}{z^2} \sum_{11}^{(2)}(z) + \frac{q^2}{2} \left[ -\sum_{11}^{(2)}(z) + \sum_{22}^{(2)}(z) + \sum_{44}^{(2)}(z) \right] + O(q^4) \quad (154) $$

En intégrant sur $\epsilon$ comme précédemment (voir la note 52), nous aboutissons à

$$ \int_{\mathbb{R}^d} [-E_{\text{bord}}(\mathbf{q}), E_{\text{bord}}(\mathbf{q})] \frac{d\epsilon}{(-\pi)} \operatorname{Im}[\tilde{\chi}_{\rho\rho}(\mathbf{q}, \epsilon + i0^+)] = \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \left\{ \frac{4\Delta^2(\epsilon_+\epsilon_- + \xi_+\xi_- + \Delta^2)}{(\epsilon_+ + \epsilon_-)\epsilon_+\epsilon_-} + \frac{(\epsilon_+ + \epsilon_-)}{\epsilon_+\epsilon_-} \left[ -(\epsilon_+\epsilon_- + \xi_+\xi_- + \Delta^2) + (\epsilon_+\epsilon_- + \xi_+\xi_- - \Delta^2) + (\epsilon_+\epsilon_- - \xi_+\xi_- - \Delta^2) \right] \right\} + O(q^4) \quad (155) $$
Ceci est nul à l'ordre $q^2$, même avant intégration sur l'angle et sur le module de $k$, comme le montre un développement limité de l'intégrande. On a donc bien

$$ \mathcal{R}^{C_0}(\mathbf{q}) = O(q^4) \qquad (156) $$

En seconde application, nous donnons la limite de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ à vecteur d'onde nul, par simplement regroupement des résultats (142) et (151) dans l'expression (68) :

$$ \lim_{q \to 0^+} \delta \tilde{g}_{\text{coll}}(\mathbf{q}) = \frac{k_B T}{4mc^2} \rho_{\text{sta}} - \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{4c_{\mathbf{k}}^2} \qquad (157) $$

avec

$$ \int_{\mathbb{R}^d} \frac{d^d k}{(2\pi)^d} \frac{\Delta^2}{4c_{\mathbf{k}}^2} = \begin{cases} \frac{m\Delta}{8\pi\hbar^2} \left(\frac{\pi}{2} + \arctan \frac{\mu}{\Delta}\right) & \text{si } d=2 \\ \frac{m^{3/2}\Delta}{8\pi\hbar^3} \left(\mu + \sqrt{\mu^2 + \Delta^2}\right)^{1/2} & \text{si } d=3 \end{cases} \qquad (158) $$

La limite est donc négative à suffisamment basse température.

## Annexe B. Développement à grand $q$ de $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$

Comme il est expliqué dans la section 3.3, nous sommes à l'ordre zéro en $\exp(-\beta E_{\text{gap}})$ (tous les nombres d'occupation de Fermi-Dirac $f_k$, $f_\pm$ sont mis à zéro) et nous devons développer la partie collective $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ (55) des susceptibilités densité-densité dans la limite $q \to +\infty$ en puissances de $q^{-1}$ à énergie complexe réduite $Z = z/E_q$ fixée (ici $E_q = \hbar^2 q^2/2m$ comme dans l'équation (15) et l'on prend $Z \in \mathbb{C} \setminus \mathbb{R}$ pour éviter les singularités – pôles ou lignes de coupure – sur l'axe réel), en allant à un ordre assez élevé pour que la contribution correspondante à $\delta \tilde{g}_{\text{coll}}(\mathbf{q})$ – obtenue par intégration curviligne sur le chemin $C_d$ comme dans (114) – soit non nulle.

Nous utilisons pour cela la méthode de développement multiéchelle de la section 4.3.5.4 de la référence [5], déjà esquissée dans notre section 3.3, (i) en limitant l'intégration sur $\mathbf{k}$ au demi-espace $\mathbf{k} \cdot \mathbf{q} > 0$ au prix d'un facteur deux global, (ii) en développant l'intégrande en puissances de $q^{-1}$ à $\bar{\mathbf{k}}_- \equiv \mathbf{k} - \mathbf{q}/2$ fixé dans la boule $k_- < \eta q$, $\eta \ll 1$, (iii) en développant l'intégrande sur le complémentaire de la boule en puissances de $q^{-1}$ à $\mathbf{K}_- = \mathbf{k}_-/q$ fixé, (iv) en redéveloppant l'ensemble en puissances de $q^{-1}$, en particulier en écrivant $\int_0^{\eta q} dq f(q) = \int_0^{+\infty} dq f(q) - \int_{\eta q}^{+\infty} dq f(q)$ et en remplaçant $f(q)$ par son développement asymptotique dans la troisième intégrale, et enfin (iv) en prenant la limite $\eta \to 0^+$ dans les coefficients de $q^{-n}$.<fnref n="53" />

### B.1. Deux astuces simplificatrices et une table récapitulative

Les calculs à venir sont grandement simplifiés par deux astuces.

Premièrement, on peut transformer avantageusement l'expression de $\tilde{\chi}_{\text{coll}}(\mathbf{q}, z)$ en introduisant des demi-sommes et des demi-différences comme suit :

$$ \alpha(\mathbf{q}, z) \equiv \frac{1}{2} [\Sigma_{23}(\mathbf{q}, z) + \Sigma_{13}(\mathbf{q}, z)] \quad ; \quad D_{12}(\mathbf{q}, z) \equiv \frac{1}{2} [\Sigma_{11}(\mathbf{q}, z) - \Sigma_{22}(\mathbf{q}, z)] \quad ; \quad (159) $$

$$ S_{12}(\mathbf{q}, z) \equiv \frac{1}{2} [\Sigma_{11}(\mathbf{q}, z) + \Sigma_{22}(\mathbf{q}, z)] + \Sigma_{12}(\mathbf{q}, z) \qquad (160) $$

<footnote n="53">⁵³Pour mettre en œuvre cette méthode, nous remplaçons au besoin le contre-terme $1/2\epsilon_{\mathbf{k}}$ par $1/4\epsilon_+ + 1/4\epsilon_-$ dans la définition (43,44) de $\Sigma_{11}(\mathbf{q}, z)$ et $\Sigma_{22}(\mathbf{q}, z)$, afin que les seuls points singuliers dans la limite $q \to +\infty$ à $\mathbf{k}/q$ fixé soient les $\mathbf{k} = \pm \mathbf{q}/2$ (et pas $\mathbf{k} = 0$). Nous y sommes autorisé par l'identité $\int_{\mathbb{R}^d} [d^d k/(2\pi)^d](2e_{\mathbf{k}}^{-1} - e_{+}^{-1} - e_{-}^{-1}) = 0$ que nous démontrons en introduisant $\Phi(\mathbf{r}) = \int_{\mathbb{R}^d} [d^d k/(2\pi)^d](2e_{\mathbf{k}}^{-1} - e_{+}^{-1} - e_{-}^{-1}) \exp(\mathbf{i}\mathbf{k} \cdot \mathbf{r}) = 2[1 - \cos(\mathbf{q} \cdot \mathbf{r}/2)]\phi(\mathbf{r})$ avec $\phi(\mathbf{r}) = \int_{\mathbb{R}^d} [d^d k/(2\pi)^d] \exp(\mathbf{i}\mathbf{k} \cdot \mathbf{r})/\epsilon_{\mathbf{k}}$, puis en remarquant que $\Phi(\mathbf{r}) \to 0$ lorsque $r \to 0$ puisque $\phi(\mathbf{r}) = O(\ln r)$ à 2D et $\phi(\mathbf{r}) = O(1/r)$ à 3D.</footnote>