最近做东大过去问复习四大力学，偶遇磁场中玻色爱因斯坦凝聚（BEC），拼尽全力无法战胜，看了俩篇文献总算是明白哪里想错了，特此记录。这也是本人第一次真正看凝聚态相关的研究论文。
BEC回顾
考虑理想玻色气体，满足Bose-Einstein分布：
\[n(\epsilon)=\frac{1}{\mathrm{e}^\frac{\epsilon-\mu}{k_B T}-1}\]比较关键的就是此时的态密度：
\[D(\epsilon) = \frac{2\pi V}{h^3}(2m)^{3/2}\epsilon^{1/2}\]这里我们假设自旋为0，如果是$s$，只需要乘上一个$2s+1$就好了。[^ 1]然后我们naive写下下面这个式子计算总粒子数：
\[N=\int_{0}^{\infty}d\epsilon \frac{D(\epsilon)}{\mathrm{e}^\frac{\epsilon-\mu}{k_B T}-1}\]但是$D(0)=0$，而显然$n(0)\neq 0$，也就是说基态的贡献直接被我们忽略了！，而真正的粒子数应当用下面的式子计算：
\[N = \int_{0}^{\infty}d\epsilon \frac{D(\epsilon)}{\mathrm{e}^\frac{\epsilon-\mu}{k_B T}-1} + \frac{1}{\mathrm{e}^\frac{-\mu}{k_B T}-1}\]如果温度比较高，也就是弱简并条件下，这个时候$N_0/N1\ll 1 $，所以可以忽略，但是在温度降低时越来越多的粒子跃迁到基态，直到某个临界温度$T_c$，化学势升高为$0$1，粒子全部聚集在基态，这时就称发生了玻色爱因斯坦凝聚。为了计算所有的热力学量，不妨计算体系巨配分函数：
\[\log \Xi = -\frac{V}{\lambda_T^3}\int_{0}^{\infty}\log (1-e^{-\alpha-x})x^{1/2}dx\]其中已进行积分换元$x\equiv \beta \epsilon$，其中：
\[\alpha\equiv - \frac{\mu}{k_B T},\quad \lambda_T \equiv \frac{h}{(2\pi m k_B T)^{1/2}}\]将积分中的对数展开：
\[\log(1-e^{-\alpha-x})=-\sum_{j=1}^{\infty}\frac{e^{-j(\alpha + x)}}{j}\]然后逐项积分得到：
\[\log\Xi = \frac{V}{\lambda_T^3}F_{5/2}(-\frac{\mu}{k_B T})\]这里$F_s$用对数积分函数$\mathrm{Li}_s$定义：
\[F_s(x) = \mathrm{Li}_s(e^{-x}) = \sum_{j=1}^{\infty}\frac{(e^{-x})^j}{j^s}\]这个“让数学家疯狂的操作”实际上是不疯狂的，这个级数在$\mu\leq 0$的范围内实际上一致收敛，所以上面逐项积分的办法实际上是完全ok的。当然别忘了上面的计算依旧由于态密度在基态为0，而没考虑基态的贡献，所以巨配分函数应该修正为：
\[\boxed{ \log\Xi = \frac{V}{\lambda_T^3}F_{5/2}\left(-\frac{\mu}{k_B T}\right)-\log(1-e^{-\alpha}) }\]由此便得到：
\[N=k_B T\frac{\partial}{\partial \mu}\log\Xi = \frac{V}{\lambda_T^3}F_{3/2}\left(-\frac{\mu}{k_B T}\right) + \frac{1}{e^\alpha -1}\]上式第一项是$\epsilon>0$的贡献，后面一项是基态的贡献，那么我们就可以通过直接把第二项忽视求出BEC时的临界温度，而且这个时候$\mu = 0$：
\[T_c = \frac{2\pi}{\zeta(3/2)^{3/2}}\frac{\hbar^2}{mk_B}\left(\frac{N}{V}\right)^{2/3}\]注意这里利用了$F_s(0)=\zeta(s)$。在$T_c$附近，由于$F_{3/2}(x)$的图像长下面的样子：2
所以在临界温度附近，也就是$x\approx 1$的区域，$F_s$其实对宗量的取值并不感冒，所以在临界温度附近考虑时我们完全可以继续认为$\mu = 0$，只是温度在变化。这一点后面也会用到，所以得到凝聚到基态的粒子数的比例为：
\[\boxed{ \frac{N_0}{N} = 1-\left(\frac{T}{T_c}\right)^{3/2} }\]磁场中的BEC
现在考虑磁场中的BEC，考虑粒子的自旋为$1$，完全可以考虑更高自旋的粒子，分析方法完全一样，只需要在后面的公式中多添加几项就好了。塞曼效应会导致自旋简并度消除，而且量子力学中学习过，这还和磁场大小有关，因为精细结构和磁场的修正如果是同阶，塞曼效应就会变得比较复杂。这里我们考虑磁场比较小的最简单情况，这个时候塞曼能量修正只和自旋磁量子数有关：
\[\Delta E= c B m_z,\quad m_z = 0,\pm 1\]这只是一阶修正，二阶修正正比于$(B m_z)^2$，这里我们就不考虑了3。
原先的自旋带来的三重简并直接劈裂成三个不同的能级。显然现在要关注态密度的变化，有如下分析：4
考虑能量$\epsilon$，他有三个可能来源，一个是$\epsilon$未劈裂，另两个来源于$\epsilon \pm cB$劈裂后的能级，态密度也就应当对应原先的这三个能级的态密度之和：
\[\tilde{D}(\epsilon)= D(\epsilon) + D(\epsilon + cB)+D(\epsilon-c B)\]这里：
\[D(\epsilon)\equiv\begin{cases}\frac{2\pi V}{h^3}(2m)^{3/2}\epsilon^{1/2},&\epsilon>0\\0,&\epsilon<0\end{cases}\]
然后就可以开始计算粒子数了，注意这里基态能量不是0，而是$-cB$：
\[\begin{aligned} N&=\int_{-cB}^{+\infty}\frac{\tilde{D}(\epsilon)}{\exp(\frac{\epsilon-\mu}{k_BT})-1}d\epsilon+N_0\\ &=\int_{0}^{+\infty}\frac{D(\epsilon)}{\exp(\frac{\epsilon-\mu}{k_BT})-1}d\epsilon+\int_{-cB}^{+\infty}\frac{D(\epsilon+cB)}{\exp(\frac{\epsilon-\mu}{k_BT})-1}d\epsilon+\int_{+cB}^{+\infty}\frac{D(\epsilon-cB)}{\exp(\frac{\epsilon-\mu}{k_BT})-1}d\epsilon+N_0\\ &=\int_{0}^{+\infty}\frac{D(\epsilon)}{\exp(\frac{\epsilon-\mu}{k_BT})-1}d\epsilon+\int_{0}^{+\infty}\frac{D(\epsilon+cB)}{\exp(\frac{\epsilon+cB-\mu-cB}{k_BT})-1}d(\epsilon+cB)+\int_{0}^{+\infty}\frac{D(\epsilon-cB)}{\exp(\frac{\epsilon-cB-\mu+cB}{k_BT})-1}d(\epsilon-cB)+N_0\\ &=\frac{V}{\lambda_T^3}\left[F_{3/2}\left(-\frac{\mu}{k_B T}\right)+F_{3/2}\left(-\frac{\mu+cB}{k_B T}\right)+F_{3/2}\left(-\frac{\mu-cB}{k_B T}\right)\right]+N_0 \end{aligned}\]注意，这个时候基态不为$0$，而BEC被忽略的是基态，不难看到确实有：
\[\boxed{ \tilde{D}(-cB)=D(0)=0 }\]在发生BEC时，基态粒子数目应当非常大，也就是说：
\[N_0=\frac{1}{\exp(\frac{-cB-\mu(T_c)}{k_BT})-1}\to+\infty\Rightarrow \boxed{\mu(T_c)=-cB}\]分析出这一点是非常重要的，由此便可以写下$T_c$满足的方程：
\[N=\frac{V}{\lambda_T^3}\left[F_{3/2}\left(\frac{cB}{k_B T_c}\right)+F_{3/2}\left(0\right)+F_{3/2}\left(\frac{2cB}{k_B T_c}\right)\right]\]前面没有磁场时的结果我们记为$T_c^0$，并引入记号$t\equiv T/T_c^c$，从而得到：
\[\begin{aligned} 1&=\frac{t_c^{3/2}}{3\zeta(\frac{3}{2})}\left[\zeta\left(\frac32\right)+F_{3/2}\left(\frac{cB}{k_BT_c}\right)+F_{3/2}\left(\frac{2cB}{k_BT_c}\right)\right]\\ &=\frac{t_c^{3/2}}{3\zeta(\frac{3}{2})}\left[\zeta\left(\frac32\right)+F_{3/2}\left(\frac{b}{t_c}\right)+F_{3/2}\left(2\frac{b}{t_c}\right)\right] \end{aligned}\]这里
\[b\equiv \frac{cB}{k_BT_c^0}\]$t_c$关于$b$的渐进行为如下：5
\[t_c\sim \begin{cases} 1+\frac{4\sqrt{\pi}(1+\sqrt{2})}{9\zeta(3/2)}\sqrt{b},&b\ll 1\\ 3^{2/3}-\frac{2}{3^{1/3}\zeta(3/2)}\exp(-b/3^{2/3}),&b\gg 1 \end{cases}\]图像也可以大致画出来：
这里红色和紫色是渐进关系图像，蓝色是直接ContourPlot出来的图像，可见BEC临界温度随着磁场的加大而加大，所以如果磁场缓慢绝热的变化，那么根据量子力学中绝热近似，如果磁场增大前体系发生BEC，那么今后粒子也将都处于基态继续发生BEC，而磁场加大导致BEC临界温度升高，所以体系实际上会因为磁场的增大而升高温度。
同样，BEC临界温度附近$N_0/N$可以根据下式计算：6
\[\frac{N_0}{N} = 1-\frac{t^{3/2}}{3\zeta(3/2)}\left[\zeta(3/2)+F_{3/2}\left(b/t\right)+\left(2b/t\right)\right]\]图像大致长这样：7
图中从左到右，$b=0,0.1,1,10$。文献53中计算了更多的热力学量，感兴趣的读者可自行查阅。
参考
-
对于玻色体系，由于$n(\epsilon)>0$，所以$\mu$恒小于0，温度降低，不断升高，直到临界温度发生BEC，$\mu(T_c) = 0$ ↩
-
MMA里面$\mathrm{Li}$叫
PolyLog↩ -
Tao C, Jiao Z, Gu Q. Reentrance of Bose-Einstein condensation in spinor atomic gases in a magnetic field. Phys Rev A 2009;79:043614. https://doi.org/10.1103/PhysRevA.79.043614. ↩ ↩2
-
注意，一直到现在我们都假设能级是准连续的 ↩
-
Simkin MV, Cohen EGD. Magnetic properties of a Bose-Einstein condensate. Phys Rev A 1999;59:1528–32. https://doi.org/10.1103/PhysRevA.59.1528. ↩ ↩2
-
同样的方法，$\mu\approx -cB$，F内宗量只随$T$改变 ↩
-
本来打算搞个动图，但是懒得折腾MMA输出GIF的方式了 ↩