贡献者： jiangnan; addis
1. 理想气体玻色爱因斯坦凝聚
首先，我们回顾一下理想玻色气体的 BEC。假设这种玻色子的色散是 $\epsilon_{ \boldsymbol{\mathbf{k}} } = \hbar^2{ \boldsymbol{\mathbf{k}} }^2/2m$（实际上长什么样在我们这里一般的分析中并没有什么关系）。我们可以看到，对于玻色-爱因斯坦分布，$n_{ \boldsymbol{\mathbf{k}} }=1/( \mathrm{e} ^{\beta(\epsilon_{ \boldsymbol{\mathbf{k}} }+\mu)}-1)$。为保证粒子占据数为正，则要求对任意动量 $ \boldsymbol{\mathbf{k}} $ 都有 $e^{\beta (\epsilon_{ \boldsymbol{\mathbf{k}} }+\mu)}-1>0$。因此当 $ \boldsymbol{\mathbf{k}} = 0$ 时，就要求 $\mu \geq 0$。读者可以看到，对于任意动量态 $ \boldsymbol{\mathbf{k}} $，在 $\mu$ 越大的情况下，粒子数占据 $n_{ \boldsymbol{\mathbf{k}} }$ 就越小，因此玻色分布有一个总粒子数的上限，即令 $\mu = 0$ 时的粒子数，
\begin{equation}
N_{max}(\beta) = \sum_{ \boldsymbol{\mathbf{k}} } \frac{1}{e^{\beta \epsilon_k}-1}~.
\end{equation}
对于量子体系，能够取到的 $k$ 的值是离散的，以一维体系为例，设体系长度为 $L$,动量 $k$ 的所有取值为 $2 n\pi/L$，即相邻动量值相差 $2\pi/L$。如果我们把这个看作是在 $k$ 空间中态的密度，则在三维情形中，单位体积里态的数目就是 $V/(2 \pi)^3$。（这里假定的波函数是有周期边界条件的，如果我们假定波函数为驻波，此处的态数目相差一个 $2^3$）。因此当体系足够大，粒子数足够多时（热力学极限下），玻色分布所能容纳的粒子数极限写为积分形式后为
\begin{equation}
N_{max}(\beta) = \frac{V}{(2\pi)^3}\int d^3 \boldsymbol{\mathbf{k}} \frac{1}{e^{\beta \epsilon_{ \boldsymbol{\mathbf{k}} }}-1}~.
\end{equation}
一旦粒子数超过这个极限，其余的粒子就会全部占据到基态，此时我们称出现了玻色-爱因斯凝聚。在发生相变临界时，如果我们将 $\beta$ 看作粒子数密度的函数，则有
\begin{equation}
\frac{1}{\beta} = k_B T_c = \frac{2\pi \hbar^2}{m}(\frac{n}{g_{3/2}(1)})^{2/3}~,
\end{equation}
其中 $g_{3/2}(1) = 2.612$ 为一玻色函数(Bose function)。该温度对应的德布罗意热波长为
\begin{equation}
\lambda_D \equiv \sqrt{\frac{2\pi \hbar^2}{m k_B T}} = (\frac{g_{3/2}(1)}{n})^{1/3}~.
\end{equation}
也即在体系的德布罗意热波长与平均粒子间距可比拟时，发生玻色爱因斯坦凝聚。在 $T< T_c$ 时，体系粒子占据分为两部分，处于热态与处于基态的粒子，
\begin{equation}
N=N_0+\frac{V}{(2\pi)^3}\int \text{d}^3{ \boldsymbol{\mathbf{k}} }\frac{1}{ \mathrm{e} ^{\beta\epsilon_{ \boldsymbol{\mathbf{k}} }}-1},\quad N_0/N\neq0~
\end{equation}
此时处于基态的粒子数呈现宏观占据。
2. 非对角长程序
看得出来我们虽然很熟悉什么是理想玻色气体的 BEC，但是对于一般的情况，如何规定它是否是处于凝聚的态呢？
首先，我们可以写出某个多体（$N$ 体）波函数 $\Psi({ \boldsymbol{\mathbf{r}} }_1,{ \boldsymbol{\mathbf{r}} }_2,\cdots,{ \boldsymbol{\mathbf{r}} }_N)$ 密度矩阵 $\rho$。它的矩阵元有 $2N$ 个指标，$\langle x_1,\cdots,x_N|\rho|y_1,\cdots,y_N\rangle$。把它 trace 掉 $N-1$ 个指标，得到
\begin{equation}
\rho({ \boldsymbol{\mathbf{r}} },{ \boldsymbol{\mathbf{r}} '}) = N\int\sum_s p_s\psi^*({ \boldsymbol{\mathbf{r}} },{ \boldsymbol{\mathbf{r}} }_2,\cdots,{ \boldsymbol{\mathbf{r}} }_N)p_s\psi({ \boldsymbol{\mathbf{r}} '},{ \boldsymbol{\mathbf{r}} }_2,\cdots,{ \boldsymbol{\mathbf{r}} }_N)\text{d}^3{ \boldsymbol{\mathbf{r}} }_2\cdots\text{d}^3{ \boldsymbol{\mathbf{r}} }_N~.
\end{equation}
二参数函数可以按照某一个基展开，
\begin{equation}
\rho({ \boldsymbol{\mathbf{r}} },{ \boldsymbol{\mathbf{r}} '}) = \sum_i N_i \varphi_i^*({ \boldsymbol{\mathbf{r}} })\varphi_i({ \boldsymbol{\mathbf{r}} '})~.
\end{equation}
对于 $N_i$，如果所有 $N_i/N=0$，那么就称这个态是 normal phase。如果有一个不等于 $0$ 的，那就是一般的 BEC，如果有多个 $N_i\neq0$ 叫 fragment BEC。后两者称之存在非对角长程序（off-diagonal long-range order，ODLRO）。
3. 超流的相互作用理论
对于正常的 BEC 的情况，$\exists N_0\gg N_{i\neq0}$，如果我们忽略两体相互作用的细节，我们就可以写出
\begin{equation}
\rho({ \boldsymbol{\mathbf{r}} , \boldsymbol{\mathbf{r}} '}) \approx N_0\varphi^*({ \boldsymbol{\mathbf{r}} })\varphi({ \boldsymbol{\mathbf{r}} '})~.
\end{equation}
其中，满足要求的系统的态是
\begin{equation}
\Psi({ \boldsymbol{\mathbf{r}} }_1,{ \boldsymbol{\mathbf{r}} }_2,\cdots,{ \boldsymbol{\mathbf{r}} }_N) = \prod_{\otimes}\varphi({ \boldsymbol{\mathbf{r}} }_i)~,
\end{equation}
我们考虑系统的 Hamiltonian。利用赝势，我们写出
\begin{equation}
\hat{H} = \sum_i^N \left(\frac{\nabla_i^2}{2m}+V({ \boldsymbol{\mathbf{r}} }_i)\right)\sum_{\langle i,j\rangle}\frac{4\pi\hbar^2 a_s}{m}\delta^3({ \boldsymbol{\mathbf{r}} }_{ij})\frac{\partial}{\partial r_{ij}}r_{ij}~.
\end{equation}
显然，我们的波函数式 9 在这里面并不存在我们上一章讨论的 $r_i\to r_j$ 的那种发散行为，因为我们实际上忽略了两体相互作用的细节。当然，这种情况下我们很容易发现，没有奇异性的时候，赝势仅仅就是一个 $\delta$ 而已，因为其形如 $r_{ij}\partial_{r_{ij}}$ 的成分在 $\delta^3({ \boldsymbol{\mathbf{r}} }_{ij})$ 下不作任何贡献：
\begin{equation}
\delta^3({ \boldsymbol{\mathbf{r}} }_{ij})\frac{\partial}{\partial r_{ij}}r_{ij} = \delta^3({ \boldsymbol{\mathbf{r}} }_{ij})~,
\end{equation}
我们可以通过对 $\langle \varphi|H+\mu N|\varphi\rangle$ 取极值得到本征值方程。
\begin{equation}
\varepsilon = N\int \text{d}^3{ \boldsymbol{\mathbf{r}} }\varphi^*({ \boldsymbol{\mathbf{r}} })\left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\varphi({ \boldsymbol{\mathbf{r}} })+\frac{N(N-1)}{2}\int \text{d}^3{ \boldsymbol{\mathbf{r}} }\frac{4\pi\hbar^2a_s}{m}|\varphi({ \boldsymbol{\mathbf{r}} })|^4~.
\end{equation}
以 $\varphi^*$ 作为变量，$N\to\infty$ 而且 $\phi = \sqrt{N}\varphi, U\equiv \dfrac{4\pi\hbar^2a_s}{m}$，我们有
\begin{equation}
\Rightarrow \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi + U|\phi|^2\phi = \mu\phi~,
\end{equation}
这就是著名的
不含时 Gross-Pitaevskii 方程。
如果我们考虑的不是一个定态问题，而是 $\phi = \phi({ \boldsymbol{\mathbf{r}} },t)$，那么我们的这个方程实际上要稍作改动。符合直觉的我们可以得到 含时 Gross-Pitaevskii 方程
\begin{equation}
\mathrm{i} \hbar\frac{\partial\phi}{\partial t} = \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi + U|\phi|^2\phi~,
\end{equation}
接下来我们做两个魔法操作。首先将 $\phi^*\times$ 式 14 减去它的共轭，显然其左手边为
\begin{equation}
\mathrm{i} \hbar\phi^*\frac{\partial\phi}{\partial t} + \mathrm{i} \hbar\phi\frac{\partial\phi^*}{\partial t} = \mathrm{i} \hbar\frac{\partial\phi^*\phi}{\partial t} = \mathrm{i} \hbar\frac{\partial|\phi|^2}{\partial t}~,
\end{equation}
而右手边则为
\begin{equation}
\begin{split}
&\left(\phi^* \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi + U|\phi|^4\right) -\left(\phi \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi^* + U|\phi|^4\right) \\
=& \frac{\hbar^2}{2m}\nabla\cdot\left(\phi\nabla\phi^* - \phi^*\nabla\phi\right)~.
\end{split}
\end{equation}
于是，就得到
\[ \mathrm{i} \hbar\frac{\partial|\phi|^2}{\partial t} = \frac{\hbar^2}{2m}\nabla\cdot\left(\phi\nabla\phi^* - \phi^*\nabla\phi\right) ~.\]
定义 $\phi({ \boldsymbol{\mathbf{r}} }) = \sqrt{\rho} \mathrm{e} ^{ \mathrm{i} \theta}$，其中 $\rho = \rho({ \boldsymbol{\mathbf{r}} }),\theta = \theta({ \boldsymbol{\mathbf{r}} })$。
我们从而有
\begin{equation}\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho{ \boldsymbol{\mathbf{v}} }),\quad { \boldsymbol{\mathbf{v}} }\overset{\text{def}}{\equiv}\frac{\hbar}{m}\frac{\partial\theta}{\partial t}~, \end{equation}
这就得到了连续性方程。当然，还得补充一个无旋方程
\begin{equation}
\nabla\times{ \boldsymbol{\mathbf{v}} }=0~,
\end{equation}
另一个魔法是将 $\phi^*\times$
式 14 它的共轭。左手边为
\begin{equation}
\mathrm{i} \hbar\phi^*\frac{\partial\phi}{\partial t} - \mathrm{i} \hbar\phi\frac{\partial\phi^*}{\partial t} = \mathrm{i} \hbar\phi^{*2}\frac{\partial \phi/\phi^*}{\partial t} = -2\hbar\rho\frac{\partial\theta}{\partial t}~,
\end{equation}
而右边就有点复杂了，下面的最后一项会很复杂：
\begin{equation}
\begin{split}
&\left(\phi^* \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi + U|\phi|^4\right) + \left(\phi \left(-\frac{\hbar^2\nabla^2}{2m}+V({ \boldsymbol{\mathbf{r}} })\right)\phi^* + U|\phi|^4\right) \\
=& 2V\rho + 2U\rho^2 - \frac{\hbar^2}{2m}\left(\phi^*\nabla^2\phi + \phi\nabla^2\phi^*\right)~.
\end{split}
\end{equation}
我们显然希望能够提取出来一个 $\rho$，所以我们给出详细的计算过程：
\begin{equation}
\begin{split}
\nabla^2\phi &= \nabla^2\left(\sqrt\rho \mathrm{e} ^{ \mathrm{i} \theta}\right)\\
&=\nabla\cdot\left(\frac{ \mathrm{e} ^{ \mathrm{i} \theta}}{2\sqrt\rho}\nabla\rho+i\sqrt{\rho} \mathrm{e} ^{ \mathrm{i} \theta}\nabla\theta\right)\\
&=\frac{ \mathrm{e} ^{ \mathrm{i} \theta}}{2\sqrt\rho}\nabla^2\rho+\frac{ \mathrm{i} \mathrm{e} ^{ \mathrm{i} \theta}}{\sqrt\rho}\nabla\theta\nabla\rho - \frac{ \mathrm{e} ^{ \mathrm{i} \theta}}{4\rho^{3/2}}(\nabla\rho)^2-\sqrt\rho \mathrm{e} ^{ \mathrm{i} \theta}(\nabla\theta)^2+i\sqrt\rho \mathrm{e} ^{ \mathrm{i} \theta}\nabla^2\theta~.
\end{split}
\end{equation}
类似有
\begin{equation}
\begin{split}
\nabla^2\phi^* &= \nabla^2\left(\sqrt\rho \mathrm{e} ^{- \mathrm{i} \theta}\right)\\
&=\nabla\cdot\left(\frac{ \mathrm{e} ^{- \mathrm{i} \theta}}{2\sqrt\rho}\nabla\rho-i\sqrt{\rho} \mathrm{e} ^{- \mathrm{i} \theta}\nabla\theta\right)\\
&=\frac{ \mathrm{e} ^{- \mathrm{i} \theta}}{2\sqrt\rho}\nabla^2\rho-\frac{ \mathrm{i} \mathrm{e} ^{- \mathrm{i} \theta}}{\sqrt\rho}\nabla\theta\nabla\rho - \frac{ \mathrm{e} ^{- \mathrm{i} \theta}}{4\rho^{3/2}}(\nabla\rho)^2-\sqrt\rho \mathrm{e} ^{- \mathrm{i} \theta}(\nabla\theta)^2-i\sqrt\rho \mathrm{e} ^{- \mathrm{i} \theta}\nabla^2\theta~,
\end{split}
\end{equation}
从而有
\begin{equation}
\phi^*\nabla^2\phi+\phi\nabla^2\phi^* = \nabla^2\rho - \frac{1}{2\rho}(\nabla\rho)^2-2\rho(\nabla\theta)^2 = 2\sqrt\rho \nabla^2\sqrt\rho - 2\rho(\nabla\theta)^2~.
\end{equation}
从而我们得到原本的式子右边等于
\begin{equation}
-\frac{\hbar^2}{2m}\frac{1}{\sqrt\rho}+\frac{\hbar^2}{2m}(\nabla\theta)^2+V+U\rho~.
\end{equation}
将 $\nabla$ 作用到等式两边，得到
\begin{equation}
m\frac{\partial \boldsymbol{\mathbf{v}} }{\partial t} = -\nabla\cdot\left[-\frac{\hbar^2}{2m}\frac{1}{\sqrt\rho}\nabla^2\sqrt\rho + \frac{1}{2}m{ \boldsymbol{\mathbf{v}} }^2+V+U\rho\right]~.
\end{equation}
以上的推导都是零温情况，只包含超流成分；有限温的问题可以加入正常流体的成分，从而变成二流体模型。考虑一个均匀系统 $V( \boldsymbol{\mathbf{r}} ) = 0$，基态有均匀的密度 $\rho_0$，且 ${ \boldsymbol{\mathbf{v}} } = 0$。我们可以对非基态的态的密度函数作展开：$\rho = \rho_0+\delta\rho$。只取 $\delta\rho,{ \boldsymbol{\mathbf{v}} }$ 和 $k$ 的领头阶，我们可以简化刚刚得到的方程：
\begin{equation}
m\frac{\partial \boldsymbol{\mathbf{v}} }{\partial t} = -U\nabla\cdot\delta\rho,\quad \frac{\partial\rho}{\partial t} = -\nabla(\rho{ \boldsymbol{\mathbf{v}} }) \approx-\rho_0\nabla \boldsymbol{\mathbf{v}} ~.
\end{equation}
联合起来就得到
\begin{equation}
\frac{\partial^2\delta\rho}{\partial t^2} = \frac{U\rho_0}{m}\nabla^2\delta\rho~,
\end{equation}
即得到了声子色散
\begin{equation}
\omega = \sqrt{\frac{U\rho_0}{m}}k~.
\end{equation}
声速得到为 $c=\sqrt{\dfrac{U\rho_0}{m}}$，是一个低能的线性 gapless 激发（从而无色散）。这与 $\text{U}(1)$ 对称性是联系起来的。
4. 超流体的朗道判据
最后我们来说一下超流与临界速度。这个临界速度就和我们 BCS 里面的临界磁场有点像：超过这个速度，超流就不再是超流了。
考虑一个质量为 $m_0$，速度 ${ \boldsymbol{\mathbf{v}} }_i$ 的杂质，在与某处于 BEC 态的物体相互作用。如果有摩擦，也就是说，这个杂质被散射到某一速度 ${ \boldsymbol{\mathbf{v}} }_f$ 了，我们可以在初始 BEC 质心系写能动量守恒
\begin{equation}
m_0{ \boldsymbol{\mathbf{v}} }_i + { \boldsymbol{\mathbf{q}} } = m_0{ \boldsymbol{\mathbf{v}} }_f,\quad \frac{m_0{ \boldsymbol{\mathbf{v}} }_i^2}{2} = \frac{m_0{ \boldsymbol{\mathbf{v}} }_i^2}{2} + c|{ \boldsymbol{\mathbf{q}} }|~.
\end{equation}
得到
\begin{equation}
{ \boldsymbol{\mathbf{v}} }_i\cdot{ \boldsymbol{\mathbf{q}} }-c|{ \boldsymbol{\mathbf{q}} }| = \frac{{ \boldsymbol{\mathbf{q}} }^2}{2m}~.
\end{equation}
可见右边等式大于 $0$，如果 ${ \boldsymbol{\mathbf{v}} }_i< c$，则左边无法满足也为正，即无法发生散射，也就是真正意义的 “超流”（不会产生摩擦作用），这个临界速度也就是声速 $c$。
5. 朗道判据的理解(comment)
朗道在判定超流的时候，假定在静止状态下，超流体的能量为 $E_g$，如果超流体中产生一个声子（元激发），动量为 p，那么超流体的能量就会增加为 $E_g + \epsilon(p)$。一般情况下，激发一个声子的能量总是要比没有声子的能量要高的，但这个情形在流动起来的超流体中并不总是成立。当超流体以速度 $V$ 运动起来后，没有激发时超流体的能量就变为 $E_g + MV^2/2$,但这时再激发一个动量为 $p$ 的声子对体系增加的能量就不再是 $\epsilon(p)$ 了。这个声子动量将会与超流体的速度耦合起来，我们首先要假定声子是一种波，这种波的动量与能量满足 $E = \hbar \omega, P = \hbar k$,这种波在传播时的波函数可以写为 $\psi(\vec{x},t) = e^{i (\omega t - \vec{k}\vec{x})}$。在进行伽利略变换时，可以认为声子波与超流体之间是没有相互作用的两个实体，它们分别进行伽利略变换，$t' = t, x' = x+vt$。变换后的波函数写为
\begin{equation}
\psi(x',t') = e^{i(\omega t'- \vec{k}\cdot (\vec{x'}-\vec{v}t))} = e^{i(\omega + \vec{k}\cdot \vec{v})t' -i \vec{k}\cdot\vec{x'}}~,
\end{equation}
这个变换等价于令 $\epsilon(p)$ 变为 $\epsilon(p)+\vec{v}\cdot \vec{p}$,也即在流动的超流体中，激发出一个声子后体系总能量变为 $E_g + MV^2/2+\epsilon(p) + \vec{v}\cdot\vec{p}$。可以看到，如果 $\epsilon(p)+\vec{v}\cdot\vec{p}<0$，那么体系激发出一个声子的能量就会小于没有声子激发的能量。朗道假定在超流体中存在一种线性激发，它的能谱为 $\epsilon(p)= c|p|$,此时激发声子能量为 $c|p|+\vec{v}\cdot\vec{p}$,如果 $v>c$ 时，体系就会倾向激发与 $v$ 相反的声子，由此每激发一个声子超流体的动量都会减小，在传播时就会产生与容器壁的摩擦，进而产生耗散。在处理元激发的时候我们有时会把元激发看做一种粒子，但是这种粒子必须要与凝聚态体系发生相互作用，这种相互作用一般而言是难以刻画的，而 $\vec{p}\cdot\vec{v}$ 就是粒子与凝聚态发生相互作用的形式，在超流体中由于可以激发声子，如果不考虑超流体与器壁的相互作用，超流体的动量 $MV$ 与声子动量 $p$ 都不再是守恒量，但是它们的总动量确实一个守恒量，且能量也是一个守恒量，所以在构造能量的时候就需要以 $(M\vec{V}+\vec{p})^2/2M$ 这样的方式构造，自然就会包含这种相互作用。当然，这是一个 argue 而非证明。