超伝導BCS理論の続きです。第二量子化入門の内容はわかっているとします。
今回は、BCS波動関数のエネルギー期待値を計算し、本当に超伝導状態のエネルギーが低いかどうかを調べます。簡単のため、一様系を考え、電子ペアはスピンシングレットの等方的なペア（s波ペア）を考えます。

前回のあらすじ

BCS波動関数は

| \psi^{BCS} \rangle =  \prod_{\boldsymbol{k}} (u_{\boldsymbol{k}}  + v_{\boldsymbol{k}} 
 c_{\boldsymbol{k} \uparrow}^{\dagger} c_{-\boldsymbol{k} \downarrow}^{\dagger})|0 \rangle,  \: \: |u_{\boldsymbol{k}} |^2 + |v_{\boldsymbol{k}}|^2 = 1

でした。ここで、

B_{\boldsymbol{k}}^{\dagger} \equiv c_{\boldsymbol{k} \uparrow}^{\dagger} c_{-\boldsymbol{k} \downarrow}^{\dagger}

を定義しておくと

| \psi^{BCS} \rangle =  \prod_{\boldsymbol{k}} (u_{\boldsymbol{k}}  + v_{\boldsymbol{k}} B_{\boldsymbol{k}}^{\dagger}) |0 \rangle,  \: \: |u_{\boldsymbol{k}} |^2 + |v_{\boldsymbol{k}}|^2 = 1

となります。

今回やること

電子がペアを組んだ時の多体の波動関数としてBCS波動関数を用意しました。この波動関数には$u_{\boldsymbol{k}}$と$v_{\boldsymbol{k}}$というパラメータがあります。実際に実現する基底状態はハミルトニアンのエネルギー期待値が最小となる状態です。ですので、ハミルトニアンをこの波動関数で挟みエネルギー期待値を計算し、その値が最小となるように$u_{\boldsymbol{k}}$と$v_{\boldsymbol{k}}$を選ぶ必要があります。この方法は変分法と呼ばれます。

BCSハミルトニアン

電子同士の相互作用が存在する時のハミルトニアンとして、以下のものを考えます。

\begin{align}
{\cal H} &= \hat{H}_1 + \hat{H}_{\rm e-e} \\
\hat{H}_1 &\equiv \sum_{\sigma=\uparrow,\downarrow} \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k} \sigma} c_{\boldsymbol{k}\sigma}^{\dagger}c_{\boldsymbol{k}\sigma} \\
\hat{H}_{\rm e-e} &\equiv \sum_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k} \boldsymbol{k}'} c_{\boldsymbol{k}' \uparrow}^{\dagger} c_{-\boldsymbol{k}' \downarrow}^{\dagger} c_{-\boldsymbol{k} \downarrow} c_{\boldsymbol{k} \uparrow}  \\
&= \sum_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k} \boldsymbol{k}'} B_{\boldsymbol{k}'}^{\dagger}  B_{\boldsymbol{k}}
\end{align}

ここで、$\xi_{\boldsymbol{k} \sigma} \equiv \epsilon_{\boldsymbol{k} \sigma} -\mu $で、$\epsilon_{\boldsymbol{k} \sigma}$は一粒子のエネルギーです。第二量子化表示の相互作用項は右から左に読み、

$c_{\boldsymbol{k} \uparrow}$: 波数kスピンアップの電子を消滅させ
$c_{\boldsymbol{k} \uparrow}$: 波数kスピンダウンの電子を消滅させ
$c_{-\boldsymbol{k}' \downarrow}^{\dagger}$: 波数kスピンダウンの電子を生成させ
$c_{\boldsymbol{k}' \uparrow}^{\dagger} $: 波数kスピンアップの電子を生成させる

という相互作用ですので、波数$k$スピンアップの電子と波数$-k$スピンダウンの電子が相互作用し、波数$k'$スピンアップと波数$-k'$スピンダウンに変化する、という意味になります。
このハミルトニアンを使って、エネルギー期待値：

E(u_{\boldsymbol{k}},v_{\boldsymbol{k}}) = \langle  \psi^{BCS}  | {\cal H} | \psi^{BCS} \rangle

を計算します。そして$E(u_{\boldsymbol{k}},v_{\boldsymbol{k}})$が最小となるような$u_{\boldsymbol{k}}$と$v_{\boldsymbol{k}}$を探します。

便利な関係式

エネルギー期待値を計算する前に、わかっていると便利な関係式を用意します。ここにあげた関係式は全てフェルミオンの生成消滅演算子に関する反交換関係:

\begin{align}
[c_{\boldsymbol{k} \sigma}^{\dagger},c_{\boldsymbol{k}' \sigma'}]_+ &= c_{\boldsymbol{k} \sigma}^{\dagger} c_{\boldsymbol{k}' \sigma'} + c_{\boldsymbol{k}' \sigma'} c_{\boldsymbol{k} \sigma}^{\dagger} = \delta_{\sigma \sigma'} \delta_{\boldsymbol{k} \boldsymbol{k}'} \\
[c_{\boldsymbol{k} \sigma},c_{\boldsymbol{k}' \sigma'}]_+ &=0 \\
[c_{\boldsymbol{k} \sigma}^{\dagger},c_{\boldsymbol{k}' \sigma'}^{\dagger}]_+ &=0
\end{align}

から導出することができます（導出は例えば丹羽　雅昭著　超伝導の基礎にあります）。

[A,B]_- = AB - BA

として、以下の関係式が成り立ちます。

\begin{align}
[c_{\boldsymbol{k} \sigma},B_{\boldsymbol{k}'}^{\dagger}]_- &= c_{-\boldsymbol{k}',\downarrow}^{\dagger} \delta_{\boldsymbol{k},\boldsymbol{k}'} \delta_{\sigma \uparrow} - c_{\boldsymbol{k}' \uparrow}^{\dagger} \delta_{\boldsymbol{k},-\boldsymbol{k}} \delta_{\sigma \downarrow} \\
[c_{\boldsymbol{k} \sigma}^{\dagger},B_{\boldsymbol{k}'}]_- &= -c_{-\boldsymbol{k}',\downarrow}^{\dagger} \delta_{\boldsymbol{k},\boldsymbol{k}'} \delta_{\sigma \uparrow} + c_{\boldsymbol{k}' \uparrow}^{\dagger} \delta_{\boldsymbol{k},-\boldsymbol{k}} \delta_{\sigma \downarrow} \\
[c_{\boldsymbol{k} \sigma},B_{\boldsymbol{k}'}]_- &= 0 \\
[c_{\boldsymbol{k} \sigma}^{\dagger}, B_{\boldsymbol{k}'}^{\dagger}]_- &= 0
\end{align}

また、以下の関係式も有用です。

\begin{align}
[c_{\boldsymbol{k}\sigma}^{\dagger} c_{\boldsymbol{k} \sigma},B_{\boldsymbol{k}'}^{\dagger}]_- &= B_{\boldsymbol{k}}^{\dagger} \delta_{\boldsymbol{k},\boldsymbol{k}'} \delta_{\sigma \uparrow} + B_{\boldsymbol{k}'}^{\dagger} \delta_{\boldsymbol{k},-\boldsymbol{k}'} \delta_{\sigma \downarrow} \\
[c_{\boldsymbol{k} \sigma}^{\dagger} c_{\boldsymbol{k} \sigma},B_{\boldsymbol{k}'}]_- &= -B_{\boldsymbol{k}'} \delta_{\boldsymbol{k},\boldsymbol{k}'} \delta_{\sigma \uparrow} - B_{\boldsymbol{k}'} \delta_{\boldsymbol{k},-\boldsymbol{k}'} \delta_{\sigma \downarrow} \\
[B_{\boldsymbol{k}},B_{\boldsymbol{k}'}^{\dagger}]_- &= (1 - (n_{\boldsymbol{k} \uparrow} + n_{\boldsymbol{k} \downarrow})) \delta_{\boldsymbol{k},\boldsymbol{k}'} \\
[B_{\boldsymbol{k}},B_{\boldsymbol{k}'}]_- &= 0 \\
[B_{\boldsymbol{k}}^{\dagger},B_{\boldsymbol{k}'}^{\dagger}]_- &= 0
\end{align}

ここで、

n_{\boldsymbol{k} \sigma} \equiv c_{\boldsymbol{k} \sigma}^{\dagger} c_{\boldsymbol{k} \sigma}

と定義しました。

エネルギー期待値

第一項の評価

それでは、エネルギー期待値を計算していきましょう。まず、ハミルトニアンの一項目$H_1$に関する期待値を計算します。これは

\begin{align}
\langle  \psi^{BCS}  | \hat{H}_1 | \psi^{BCS} \rangle &= \sum_{\sigma=\uparrow,\downarrow} \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k} \sigma} \langle  \psi^{BCS} |c_{\boldsymbol{k}\sigma}^{\dagger}c_{\boldsymbol{k}\sigma} | \psi^{BCS} \rangle 
\end{align}

となりますから、$\langle \psi^{BCS} |c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} | \psi^{BCS} \rangle $を計算します。
これは

\alpha_{\boldsymbol{k}} \equiv  u_{\boldsymbol{k}}  + v_{\boldsymbol{k}} B_{\boldsymbol{k}}^{\dagger}

を定義すると、

\begin{align}
\langle  \psi^{BCS} |c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} | \psi^{BCS} \rangle &= \langle 0| \left( \prod_{\boldsymbol{k}_1} \alpha_{\boldsymbol{k}_1}^{\dagger} \right)c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} \left( \prod_{\boldsymbol{k}_2} \alpha_{\boldsymbol{k}_2} \right) |0 \rangle 
\end{align}

となりますが、ここで、波数を$\boldsymbol{k}^1, \boldsymbol{k}^2,\cdots$とすると、

\begin{align}
&=  \langle 0| \left( \alpha_{\boldsymbol{k}^1}^{\dagger}  \alpha_{\boldsymbol{k}^2}^{\dagger} \cdots  \right)c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} \left( \alpha_{\boldsymbol{k}^1}\alpha_{\boldsymbol{k}^2} \cdots \right)  |0 \rangle
\end{align}

となります。次に、$\boldsymbol{k} \neq \boldsymbol{k}'$の時は上で導入した関係式を用いれば

\begin{align}
[c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} ,\alpha_{\boldsymbol{k}'}]_- &= 
 v_{\boldsymbol{k}} [c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} ,B_{\boldsymbol{k}'}^{\dagger}]_- =0\\
\alpha_{\boldsymbol{k}}^{\dagger}\alpha_{\boldsymbol{k}'} &= \alpha_{\boldsymbol{k}'} \alpha_{\boldsymbol{k}}^{\dagger}, \: \boldsymbol{k} \neq \boldsymbol{k}'
\end{align}

となりますので、右側の$\alpha_{\boldsymbol{k}^i}$は$\boldsymbol{k}^i = \boldsymbol{k}$となる$\boldsymbol{k}^i$以外を左に寄せることができます。その際、$\alpha_{\boldsymbol{k}^i}^{\dagger}$を入れ替えないように$\alpha_{\boldsymbol{k}^i}$の右に置くことにします。つまり、

\begin{align}
\langle  \psi^{BCS} |c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} | \psi^{BCS} \rangle &=  \langle 0| \left( \prod_{\boldsymbol{k}_1 \neq \boldsymbol{k}} \alpha_{\boldsymbol{k}_1}^{\dagger}  \alpha_{\boldsymbol{k}_1} \right)\alpha_{\boldsymbol{k}}^{\dagger} c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow}  \alpha_{\boldsymbol{k}}  |0 \rangle
\end{align}

のようにできます。次に、

\begin{align}
\alpha_{\boldsymbol{k}_1}^{\dagger}  \alpha_{\boldsymbol{k}_1} &= |u_{\boldsymbol{k}_1}|^2 + u_{\boldsymbol{k}_1}^{\ast} v_{\boldsymbol{k}_1} B_{\boldsymbol{k}_1}^{\dagger} + v_{\boldsymbol{k}_1}^{\ast} u_{\boldsymbol{k}_1} B_{\boldsymbol{k}_1} + |v_{\boldsymbol{k}_1}|^2 B_{\boldsymbol{k}_1} B_{\boldsymbol{k}_1}^{\dagger}
\end{align}

となりますが、第二項は$\langle 0 | B_{\boldsymbol{k}}^{\dagger} = 0$によって消え、第三項は右端まで動かすことができ、$B_{\boldsymbol{k}} |0 \rangle = 0$によって消えます。第一項と第四項は$|u|^2+|v|^2=1$から1となります。
以上から、

\begin{align}
\langle  \psi^{BCS} |c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} | \psi^{BCS} \rangle &=  \langle 0| \alpha_{\boldsymbol{k}}^{\dagger} c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow}  \alpha_{\boldsymbol{k}}  |0 \rangle \\
&= \langle 0| (u_{\boldsymbol{k}}^{\ast} + v_{\boldsymbol{k}}^{\ast} B_{\boldsymbol{k}}) c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow}  (u_{\boldsymbol{k}} + v_{\boldsymbol{k}} B_{\boldsymbol{k}}^{\dagger})|0 \rangle \\
&= \langle 0|\left( |u_{\boldsymbol{k}}|^2 c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} +u_{\boldsymbol{k}}^{\ast} v_{\boldsymbol{k}} c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow}  
 B_{\boldsymbol{k}}^{\dagger}+v_{\boldsymbol{k}}^{\ast} u_{\boldsymbol{k}}B_{\boldsymbol{k}}  c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow} +|v_{\boldsymbol{k}}|^2 B_{\boldsymbol{k}} c_{\boldsymbol{k}\uparrow}^{\dagger}c_{\boldsymbol{k}\uparrow}  
 B_{\boldsymbol{k}}^{\dagger} \right)|0 \rangle \\
 &=|v_{\boldsymbol{k}}|^2
\end{align}

となります。ダウンスピンに関しても同様に計算することができ、結局

\begin{align}
\langle  \psi^{BCS}  | \hat{H}_1 | \psi^{BCS} \rangle &= \sum_{\sigma=\uparrow,\downarrow} \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k} \sigma} |v_{\boldsymbol{k}}|^2
\end{align}

となります。$\xi_{\boldsymbol{k} \sigma} = \xi_{\boldsymbol{k}}$の場合には、

\begin{align}
\langle  \psi^{BCS}  | \hat{H}_1 | \psi^{BCS} \rangle &= 2 \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k}} |v_{\boldsymbol{k}}|^2
\end{align}

が得られます。

第二項の評価

次に、相互作用項の計算です。

\begin{align}
\langle  \psi^{BCS}  | \hat{H}_{\rm e-e} | \psi^{BCS} \rangle &=  \sum_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k} \boldsymbol{k}'} \langle  \psi^{BCS}  |B_{\boldsymbol{k}'}^{\dagger}  B_{\boldsymbol{k}} | \psi^{BCS} \rangle
\end{align}

となりますから、$\langle \psi^{BCS} |B_{\boldsymbol{k}'}^{\dagger} B_{\boldsymbol{k}} | \psi^{BCS} \rangle$がどうなるかをみていきます。これは

\begin{align}
\langle  \psi^{BCS}  |B_{\boldsymbol{k}'}^{\dagger}  B_{\boldsymbol{k}} | \psi^{BCS} \rangle &= \langle 0| \left( \prod_{\boldsymbol{k}_1} \alpha_{\boldsymbol{k}_1}^{\dagger} \right)B_{\boldsymbol{k}'}^{\dagger}  B_{\boldsymbol{k}}  \left( \prod_{\boldsymbol{k}_2} \alpha_{\boldsymbol{k}_2} \right) |0 \rangle 
\end{align}

となりますが、第一項の時と同様に$\boldsymbol{k}$と$\boldsymbol{k}'$以外の波数に関しては1となり、

\begin{align}
&=  \langle 0| (u_{\boldsymbol{k}}^{\ast}+v_{\boldsymbol{k}}^{\ast} B_{\boldsymbol{k}})(u_{\boldsymbol{k}'}^{\ast} + v_{\boldsymbol{k}'}^{\ast} B_{\boldsymbol{k}'}) B_{\boldsymbol{k}'}^{\dagger} B_{\boldsymbol{k}} (u_{\boldsymbol{k}'}+v_{\boldsymbol{k}'} B_{\boldsymbol{k}'}^{\dagger}) (u_{\boldsymbol{k}} + v_{\boldsymbol{k}} B_{\boldsymbol{k}}^{\dagger}) |0 \rangle
\end{align}

となります。この計算は全部で16項出てきますが、値を持つのは$c_{\boldsymbol{k}}$と$c_{\boldsymbol{k}}^{\dagger}$の数が等しいものだけです。つまり、$B_{\boldsymbol{k}}$と$B_{\boldsymbol{k}}^{\dagger}$のセットおよび$B_{\boldsymbol{k}'}$と$B_{\boldsymbol{k}'}^{\dagger}$があるものだけ残ります。つまり、

\begin{align}
&=  \langle 0| (u_{\boldsymbol{k}}^{\ast})( v_{\boldsymbol{k}'} ^{\ast}B_{\boldsymbol{k}'}) B_{\boldsymbol{k}'}^{\dagger} B_{\boldsymbol{k}} (u_{\boldsymbol{k}'}) ( v_{\boldsymbol{k}} B_{\boldsymbol{k}}^{\dagger}) |0 \rangle = \langle 0| u_{\boldsymbol{k}}^{\ast} v_{\boldsymbol{k}'}^{\ast} u_{\boldsymbol{k}'} v_{\boldsymbol{k}} B_{\boldsymbol{k}'}B_{\boldsymbol{k}'}^{\dagger} B_{\boldsymbol{k}}   B_{\boldsymbol{k}}^{\dagger} |0 \rangle \\
&= u_{\boldsymbol{k}}^{\ast} v_{\boldsymbol{k}'}^{\ast} u_{\boldsymbol{k}'} v_{\boldsymbol{k}} 
\end{align}

となりますから、

\begin{align}
\langle  \psi^{BCS}  | \hat{H}_{\rm e-e} | \psi^{BCS} \rangle &=  \sum_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k} \boldsymbol{k}'} u_{\boldsymbol{k}}^{\ast} v_{\boldsymbol{k}'}^{\ast} u_{\boldsymbol{k}'} v_{\boldsymbol{k}} 
\end{align}

が得られます。

エネルギー期待値と最小化

結局、BCSハミルトニアンによるエネルギー期待値は

\begin{align}
\langle  \psi^{BCS}  | {\cal H} | \psi^{BCS} \rangle &=  2 \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k}} |v_{\boldsymbol{k}}|^2 + \sum_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k} \boldsymbol{k}'} u_{\boldsymbol{k}}^{\ast} v_{\boldsymbol{k}'}^{\ast} u_{\boldsymbol{k}'} v_{\boldsymbol{k}} \equiv W_0 
\end{align}

です。この$W_0$を最小化するような$u_{\boldsymbol{k}}$と$v_{\boldsymbol{k}}$を探します。
ここで、$u_{\boldsymbol{k}}$は複素数ですので、

u_{\boldsymbol{k}} = |u_{\boldsymbol{k}}| e^{i \phi_{\boldsymbol{k}}}

のように書けます。しかし、$W_0$はよく見ると位相に依存していません。そこで、$u_{\boldsymbol{k}}$、$v_{\boldsymbol{k}}$を実数と仮定して最小値を探すことにします。

|u_{\boldsymbol{k}}|^2 + |v_{\boldsymbol{k}}|^2 = 1

より、

\begin{align}
u_{\boldsymbol{k}} &\equiv \cos \theta_{\boldsymbol{k}} \\
v_{\boldsymbol{k}} &\equiv \sin \theta_{\boldsymbol{k}}
\end{align}

とし、$W_0$を最小化する$\theta_{\boldsymbol{k}}$を探します。$W_0$は

\begin{align}
W_0 &= \sum_{\boldsymbol{k}} 2 \xi_{\boldsymbol{k}} \sin^2 \theta_{\boldsymbol{k}} + \sum_{\boldsymbol{k},\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \cos \theta_{\boldsymbol{k}} \sin \theta_{\boldsymbol{k}} \cos \theta_{\boldsymbol{k}'} \sin \theta_{\boldsymbol{k}'} \\
&= 2 \sum_{\boldsymbol{k}} \xi_{\boldsymbol{k}} \sin^2 \theta_{\boldsymbol{k}} + \frac{1}{4}\sum_{\boldsymbol{k},\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \sin 2\theta_{\boldsymbol{k}}   \sin 2\theta_{\boldsymbol{k}'}
\end{align}

となりますから、$\theta_{\boldsymbol{k}}$の関数です。$W_0$の最小値を見つけるには、

\frac{\partial W_0}{\partial \theta_{\boldsymbol{k}}} =0

となるような$\theta_{\boldsymbol{k}}$を見つければよいですね。よって、

\begin{align}
\frac{\partial W_0}{\partial \theta_{\boldsymbol{k}}} &= 2 \xi_{\boldsymbol{k}} \sin 2 \theta_{\boldsymbol{k}} + \frac{1}{2} \cos 2 \theta_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \sin 2 \theta_{\boldsymbol{k}'} + \frac{1}{2} \cos 2 \theta_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k}',\boldsymbol{k}} \sin 2 \theta_{\boldsymbol{k}'} \\
&= 2 \xi_{\boldsymbol{k}} \sin 2 \theta_{\boldsymbol{k}} +  \cos 2 \theta_{\boldsymbol{k}} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \sin 2 \theta_{\boldsymbol{k}'} =0
\end{align}

となります。ここで、$U_{\boldsymbol{k},\boldsymbol{k}'} = U_{\boldsymbol{k}',\boldsymbol{k}}$を使いました。この式の両辺を$2 \cos 2 \theta_{\boldsymbol{k}}$で割りますと、

\begin{align}
\xi_{\boldsymbol{k}} \tan 2 \theta_{\boldsymbol{k}} = -\frac{1}{2} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \sin 2 \theta_{\boldsymbol{k}'}
\end{align}

という式が得られます。さらに、

\Delta_{\boldsymbol{k}} \equiv -\frac{1}{2} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} \sin 2 \theta_{\boldsymbol{k}'}

という$\Delta_{\boldsymbol{k}}$を定義すると、

\begin{align}
 \tan 2 \theta_{\boldsymbol{k}} = \frac{\Delta_{\boldsymbol{k}}}{\xi_{\boldsymbol{k}}}　
\end{align}

となりますから、三角関数の公式:

\begin{align}
\cos^2 2 \theta_{\boldsymbol{k}} &= \frac{1}{1 + \tan^2 \theta_{\boldsymbol{k}}} = \frac{1}{1 + (\Delta_{\boldsymbol{k}}/\xi_{\boldsymbol{k}})^2} = \frac{\xi_{\boldsymbol{k}}^2}{ \xi_{\boldsymbol{k}}^2 + \Delta^2_{\boldsymbol{k}}}
\end{align}

を使うことで、

\begin{align}
u_{\boldsymbol{k}}^2 &= \cos^2 \theta_{\boldsymbol{k}} = \frac{1}{2} \left( 1 + \frac{\xi_{\boldsymbol{k}}}{\sqrt{\xi_{\boldsymbol{k}}^2+\Delta_{\boldsymbol{k}}^2}}  \right) \\
v_{\boldsymbol{k}}^2 &= \sin^2 \theta_{\boldsymbol{k}} = \frac{1}{2} \left( 1 - \frac{\xi_{\boldsymbol{k}}}{\sqrt{\xi_{\boldsymbol{k}}^2+\Delta_{\boldsymbol{k}}^2}}  \right)
\end{align}

が得られます。
また、

\Delta_{\boldsymbol{k}} = -\sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'}  \cos  \theta_{\boldsymbol{k}'} \sin  \theta_{\boldsymbol{k}'} = -\sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'} u_{\boldsymbol{k}'} v_{\boldsymbol{k}'}

より、

\Delta_{\boldsymbol{k}} = -\frac{1}{2} \sum_{\boldsymbol{k}'} U_{\boldsymbol{k},\boldsymbol{k}'}  \frac{\Delta_{\boldsymbol{k}'}}{\sqrt{\xi_{\boldsymbol{k}'}^2+\Delta_{\boldsymbol{k}'}^2}}

という式が得られます。
この方程式をギャップ方程式と呼びます。この式を満たすような$\Delta_{\boldsymbol{k}}$がエネルギーが最低となる$\Delta_{\boldsymbol{k}}$であり、$\Delta_{\boldsymbol{k}}$によって表現されている$u_{\boldsymbol{k}}$と$v_{\boldsymbol{k}}$が得られます。以上から、BCS波動関数の係数を決定することができました。

ギャップ方程式の解

次に、ギャップ方程式を解いてみます。ただし、そのままだと数値計算を使わないと解けませんので、簡単な相互作用を仮定して解くことにします。相互作用を

\begin{align}
U_{\boldsymbol{k},\boldsymbol{k}'} &= \left\{
\begin{array}{ll}
=U_0 & |\xi_{\boldsymbol{k}}| \le \hbar \omega_D かつ |\xi_{\boldsymbol{k}'}| \le \hbar \omega_D \\
0 & それ以外
\end{array}
\right.
\end{align}

とします。$\omega_D$とはデバイ振動数と呼ばれるもので、$\xi_{\boldsymbol{k}} = \epsilon_{\boldsymbol{k}} - \mu$ですから、フェルミエネルギー近くの$\pm \hbar \omega_D$でだけ引力相互作用が働く、という模型です。
これを使うと、ギャップ方程式は

\Delta =U_0 \sum_{\boldsymbol{k}'}' \frac{\Delta}{2 \sqrt{\xi_{\boldsymbol{k}'}^2+\Delta^2}}

　となります。ここで、$\sum_{\boldsymbol{k}'}'$は$- \hbar \omega_D \le \xi_{\boldsymbol{k}'} \le \hbar \omega_D$となる$\boldsymbol{k}'$での和をとることを意味しています。

さて、

N(\epsilon) \equiv \sum_{\boldsymbol{k}} \delta(\epsilon - \epsilon_{\boldsymbol{k}})

を導入すると、方程式は

\begin{align}
1 &=U_0 \int_{-\hbar \omega_D}^{\hbar \omega_D} d\xi \frac{N(\xi + \mu)}{2 \sqrt{\xi^2+\Delta^2}}
\end{align}

となります。さらに、$N(\epsilon)$は状態密度であることを思い返すと、三次元自由電子系であれば$N(\epsilon) = A \epsilon^{1/2}$となります。そして、$|\xi| \ll \mu$の時、$N(\xi+\mu) = N(\mu) \equiv N_0$とフェルミエネルギーでの状態密度に近似することができます。
よって、

\begin{align}
1 &=N_0 U_0 \int_{0}^{\hbar \omega_D} d\xi \frac{1}{\sqrt{\xi^2+\Delta^2}} \\
&= N_0 U_0 \ln \Big{|} \frac{\hbar \omega_D + \sqrt{(\hbar \omega_D)^2 + \Delta^2}}{\Delta} \Big{|} 
\end{align}

となります。ここで、

\int dx \frac{dx}{\sqrt{x^2+a^2}} = \ln |x + \sqrt{x^2 + a^2}|

を使いました。さらに、実際の物質では$\Delta \ll \hbar \omega_D$なので、

\sqrt{(\hbar \omega_D)^2 + \Delta^2} \sim \hbar \omega_D

とすると、$\Delta$に関して解くことができまして、

\Delta = 2 \hbar \omega_D \exp \left[ - \frac{1}{N_0 U_0}\right]

が得られます。
この結果は、$U_0$が0に近くても引力であれば解があることを意味していますから、電子間に引力があれば金属状態ではなく超伝導状態が安定であることを示しています。また、$U_0$が分母に入っていますから、$U_0$が小さいとする通常の摂動論では超伝導が理解できないこともわかります。

これで、絶対零度の時に、BCS波動関数が基底状態であることがわかりました。