固体物理学

『1週間で学べる！Julia数値計算プログラミング』の勉強メモです。

物性物理学

物質の様々な巨視的性質を微視的な観点から研究する物理学の分野を物性物理学という。巨視的性質とは、磁性や伝熱性、伝導性などのこと。この物性物理学には、ソフトマター物理学・固体物理学・表面物理学などがある。

固体物理学

固体物理学は固体と固体内に関する関する物理現象を扱う。固体の性質の多くは電子のふるまいが決める。電子や原子を古典粒子として扱う限り固体の物性は説明できず、これらのふるまいを量子力学によって記述する。

量子化学

また、原子と電子のふるまいから物性を理論的に説明する学問分野を量子化学といい、量子化学計算では主に分子軌道法と密度汎関数理論に分かれる。

シュレディンガー方程式

１電子と１原子核の系について考える。孤立原子における１電子の運動を表すシュレディンガー方程式は

\begin{multline}
\begin{split}
\mathcal{H}\phi^a(r)=E\phi^a(r)\\
\mathcal{H}=-\frac{\hbar^2}{2m}\nabla^2+V(r)
\end{split}
\end{multline}

となる。このとき、$\phi^a$は１電子の波動関数、$E$はエネルギー、$V$はポテンシャルエネルギー、$\hbar$はプランク定数、$m$は電子の質量。
イメージはこんな感じ

ボルン-オッペンハイマー近似

電子と原子核の質量や速度の差から、電子と原子核の運動を分離して、それぞれの運動を表す近似法をボルン-オッペンハイマー近似という。
電子$N$個の運動はシュレディンガー方程式によって、

\begin{multline}
\begin{split}
&\mathcal{H}\psi(r_1,\cdots,r_N)=E\psi(r_1,\cdots,r_N)\\
&\mathcal{H}=\sum_i^N\left(-\frac{1}{2m}\nabla_i^2+V_{ion}(r_i)\right)+\frac{1}{2}\sum_{i\neq k}\frac{e^2}{|r_i-r_k|}
\end{split}
\end{multline}

とかける。ここで、原子核からのポテンシャル$V_{ion}(r)=-\sum_I Z_I\frac{e^2}{|r-R_I|}$。
電子系のエネルギーなどの物性は密度関数から計算することが可能であるという密度汎関数理論は、ホーヘンベルク-コーンの定理により示される。その計算手法はコーン-シャム方程式によって与えられる。これは、

\begin{multline}
\begin{split}
&\mathcal{H}^{KS}\psi_n(r)=E_n\psi_n(r)\\
&\mathcal{H}^{KS}=-\frac{1}{2m}\nabla^2+v_{eff}(r)\\
&v_{eff}(r)=V_{ion}(r)+\int dr'\frac{\rho(r')}{|r-r'|}+v_{XC}(\rho(r))
\end{split}
\end{multline}

と書ける。ここで、固有値$E_n$は$E_1<E_2<\cdots $となり、$E\simeq\sum_n E_n$となる。また、$\psi_n(r)$はそれぞれの固有ベクトルとなる。
コーン-シャム方程式はそのままではポテンシャルなどが考えらず、計算できない。そのため、何らかの数理モデルで近似してやる必要がある。

強束縛模型

結晶中の１電子の波動関数$\psi(r)$を孤立原子の波動関数$\phi^a(r)$の重ね合わせによって近似することを強束縛模型(強束縛近似)という。結晶中の各原子からのポテンシャルが周期的になるため、ハミルトニアン$H$は、

H(r)=\sum_{R_n}H^a(r-R_n)+\Delta U(r)

と近似できる。ここで$H^a$は孤立原子のハミルトニアン、$R_n$はポテンシャルの周期、$\Delta U$は微小量。この波動方程式は、

\psi(r)=\sum_{m,R_n}c_{m}(R_n)\phi_m^a(r-R_n)

と表す。ここで、$m$はエネルギー準位。
エネルギー準位が$N$種類あるとき、$M(<N)$本のエネルギー準位で再現する。
イメージはこんな感じ

ブロッホの定理

周期的なポテンシャルのハミルトニアンにおいて、その固有関数となるような１電子波動関数は、空間をポテンシャルの周期$R_n$だけ併進させる演算子$\hat{T_{R_n}}$に対しても固有関数になる。$\hat{T_{R_n}}$の固有値は$e^{ik{R_n}}$の形になるため、以下が成り立つ。

\psi_n(r+R_{n})=e^{ikR_n}\psi_n(r)

これをブロッホの定理という。
周期的に並んでいる結晶において、ポテンシャルは$v_{eff}(r)=v_{eff}(r+R_n)$が成り立つ。このとき、シュレディンガー方程式の解は

\psi_n^k(r)=e^{ikr}u_n^k(r)

となる。ここで、$u_n^k(r)$は並進対称性$u_n^k(r+R_n)=u_n^k(r)$をもつ周期関数、$n$はエネルギー準位、$k$は波数。この$\psi_n^k(r)$をBloch関数と呼ぶ。

ワニエ関数

Bloch関数は波数$k$に関する関数であり、強束縛模型では格子点に局在した状態をもとに近似を考えたい。格子点に局在した関数として、Bloch関数をフーリエ変換したものをWannier関数という。

\omega_{n,R}(r)=\frac{1}{V_{BZ}}\int_{BZ} e^{-ikR}\psi_n^k(r)d^3k

ここで、BZはブリアルゾーン、$R$は格子空間となる。
また、格子点の数$N$が十分大きいので、積分を和で表現することができる。

\omega_{n,R}(r)=\frac{1}{N}\sum_{BZ} e^{-ikR}\psi_n^k(r)

このWannier関数を用いて、Bloch関数は以下のようになる。

\psi_n^k(r)=\sum_{R}e^{ikR}\omega_{n,R}(r)

Wannier関数には、並進性・正規直交性・局在性の性質を持つ。また、孤立した原子の極限を考えると、Wannier関数は原子軌道に一致するはずであるため、Wannier関数の近似として原子軌道が有効となる。

LCAO近似

Wannier関数を原子軌道の線形結合(Linear Combination of Atomic Orbitals)で近似すると、

\omega_{n,R}(r)=\sum_m c_m\phi_m(r-R)

とできる。ここで、$\phi_m$は既知の関数系で一般的には直交しない。また、$c_m$は規格化条件$\sum_m|c_m|^2=1$を満たす。結晶内の電子の固有状態、Bloch関数は

\psi_n^k(r)=\sum_Re^{ikR}\sum_m c_m \phi_m^k(r-R)

となる。ここで、和の順番を入れ替えると、

\psi_n^k(r)=\sum_m c_m\sum_Re^{ikR}\phi_m(r-R)\equiv\sum_m c_m\Phi_{m}^k(r)

となる。$\sum_Re^{ikR}\phi_m(r-R)$をブロッホ和といい、ここで強束縛模型となっていることがわかる。
これをコーン-シャム方程式に代入し、両辺に左から$\Phi_l^k(r)$をかけて、全空間で積分すると、

\int dr \sum_m c_m \Phi_l^k(r) \mathcal{H}^{KS}\Phi_{m}^k(r)=E_n^k\int dr\sum_m c_m \Phi_l^k(r) \Phi_{m}^k(r)

となる。ここで下記のように行列を定義すると、

\begin{multline}
\begin{split}
(H_n^k)_{lm}&\equiv \int \Phi_l^k(r) \mathcal{H}^{KS}\Phi_{m}^k(r) dr\\
(S_n^k)_{lm}&\equiv \int \Phi_l^k(r) \Phi_{m}^k(r) dr
\end{split}
\end{multline}

係数$c_m$を並べたベクトル$\boldsymbol{c}$を使って、

H_n^k\boldsymbol{c}=E_n^k S_n^k\boldsymbol{c}

とできる。行列$S^k$を単位行列をできれば、固有値問題とでき、Bloch関数を得ることができる。

強束縛模型のハミルトニアンを解く

一番シンプルな強束縛模型のハミルトニアンは$M=1$とした1バンド模型で、

(\epsilon^k-\mu)\boldsymbol{c}^k=E_n^k\boldsymbol{c}^k

となる。ここで、$\epsilon^k$はバンド分散。固有値は$E_n^k=\epsilon^k-\mu$となる。系全体のエネルギーは

E=\int dk \sum_n E_n^k

となる。一番小さなエネルギーが基底状態となる。もし系の粒子数を考えたいなら、化学ポテンシャルの値を調整すればよい。

欲しい粒子数を与える化学ポテンシャルを計算する

一番簡単なバンド分散として、１次元のcosバンド模型$\epsilon^k=-2cos(k_xa)$を考える。$a$は原子間距離。
エネルギーは

E=\frac{1}{2\pi}\int_{-\pi/a}^{\pi/a}dk_x(\epsilon_k-\mu)\Theta(-\epsilon_k+\mu)

となる。ここで$\Theta(x)$はヘビサイド階段関数で、$x<0$のとき$0$、それ以外の時に$1$を返す。
波数空間の積分を$M$点の和に置き換える。

\frac{1}{2\pi}\int_{-\pi/a}^{\pi/a}dk_x f(k_x)=\frac{1}{M}\sum_{i=1}^M f(k_{x_i})

ここで、$k_{x_i}=\frac{2\pi}{a}\frac{i-1}{M-1}-\frac{\pi}{a}$となる。
また、ありうるすべての電子のエネルギー順位がどれくらい占有されているか、電子数の期待値をフィリングといい、

filling\equiv\frac{1}{2\pi}\int_{-\pi/a}^{\pi/a}dk_x\Theta(-\epsilon_k+\mu)

と定義される。

エネルギーとフィリングを計算

エネルギーとフィリングは以下のソースコードで求めることができる。

function calc_energy(M,myu,epsilon)
    E = 0
    ks = range(-pi,pi,length=M)
    filling = 0
    for kx in ks
        epsilonkx = epsilon(kx)-myu
        if epsilonkx < 0
            E += epsilonkx
            filling += 1
        end
    end
    return E/M,filling/M
end

ある化学ポテンシャルを与えた時、分割数$M$を増やした時のエネルギーの値をプロットすると下図のようになる。

フィリングと化学ポテンシャルの関係ををプロットすると下図のようになる。

化学ポテンシャルを求める

指定したフィリングになる化学ポテンシャルの値を求めるにはフィリング一定の直線とグラフの好転を求めればよい。二分法を用いて求める。ソースコードは以下の通り。

function bisection_method(xmin,xmax,f,eps;itamax=20)
    fmin = f(xmin)
    fmax = f(xmax)
    for i = 1:itamax
        xmid = (xmin+xmax)/2
        fmid = f(xmid)
        if abs(fmid) < eps
            return xmid,fmid
        end
        if fmid < 0
            xmin = xmid
        else
            xmax = xmid
        end
    end
end

二分法とは

解を含む区間の中間地点を求める操作を繰り返すことで方程式を解く求根アルゴリズム。
初期条件の区間$[x_{min},x_{max}]$において、$n$回繰り返した時の区間幅は、

w_n = \frac{x_{max}-x_{min}}{2^n}

となる。つまり、この値が解の精度となる。精度$10^{-m}$ほしい時、反復回数は、$n=m\log_2 10$となる。
ソート済み配列データに他する探索アルゴリズムとして二分探索があるが、同じ発想となる。

二分法の様子は下図のようになる。

強束縛模型のエネルギーを数値計算する

1バンド模型だと$\epsilon^k$が決まった瞬間に固有値と固有ベクトルが決まってしまう。しかし、非一様な系、不純物が入った系は数値計算で求める必要がある。
規則正しく並んでいる原子のうち１つを別の原子に置き換えたものだとする。このとき、不純物無しで作った強束縛模型に不純物としてポテンシャルを導入し、計算を行う。強束縛模型に不純物ポテンシャルを入れて得られる計算結果は走査型トンネル顕微鏡測定で実際に測ることができる。この測定で測ることができる量を実際に数値計算してみる。
不純物ポテンシャルを導入するために、波数空間で定義されたハミルトニアンを逆フーリエ変換して実空間表示にすればよい。１次元cosバンド模型$\epsilon^k=-2\cos(k_xa)$の逆フーリエ変換を考える。ここで、$i$番目の原資の座標を$R_i$、原子の総数を$N$とすると、

\begin{multline}
\begin{split}
\sum_{i=1}^{L_x}e^{ik_xR_i}\left(-2\cos(k_xa)-\mu\right)\boldsymbol{c}_{R_i}&=E_ne^{ikR_i}\boldsymbol{c}_{R_i}\\
\int dk_x\sum_{i=1}^{N}\left(-2\frac{e^{ik_x(R_i+a-R_j)}+e^{ik_x(R_i-a-R_j)}}{2}-\mu e^{ik_x(R_i-R_j)} \right)\boldsymbol{c}_{R_i}&=\int dk_x \sum_{i=1}^{N} E_ne^{ik_x(R_i-R_j)}\boldsymbol{c}_{R_i}\\
-\sum_{l=1}^2\boldsymbol{c}_{R_j+a_l}-\mu\boldsymbol{c}_{R_j}&=E_n \boldsymbol{c}_{R_j}
\end{split}
\end{multline}

となり、$N\times N$の行列の固有値問題となる。ここで、$a_1=a、a_2=-a$、$i=1$番目と$i=N$番目のサイトはつながっている周期境界条件とする。
ここに、座標$R_0$にある不純物ポテンシャル$V(R)=V_0\delta_{R,R_0}$を足すと、

-\sum_{l=1}^2\boldsymbol{c}_{R_j+a_l}+(V_0\delta_{R_j,R_0}-\mu)\boldsymbol{c}_{R_j}=E_n \boldsymbol{c}_{R_j}

となる。また、２次元の場合は、

-\sum_{l=1}^4\boldsymbol{c}_{R_j+a_l}+(V(R)-\mu)\boldsymbol{c}_{R_j}=E_n \boldsymbol{c}_{R_j}

となる。ソースコードは以下のようになる。

function make_H_2d(Lx,Ly,myu,V)
    N = Lx*Ly
    H = zeros(Float64,N,N)
    ds = [(1,0),(-1,0),(0,1),(0,-1)]
    for ix = 1:Lx
        for iy = 1:Ly
            i = (iy-1)*Lx+ix
            j = 0
            for d in ds
                jx = ix+d[1]
                jx += ifelse(jx > Lx,-Lx,0)
                jx += ifelse(jx < 1,Lx,0)
                jy = iy+d[2]
                jy += ifelse(jy > Ly,-Ly,0)
                jy += ifelse(jy < 1,Ly,0)
                j = (jy-1)*Lx+jx
                H[i,j] += -1
            end
            H[i,i] += -myu+V(ix,iy)
        end
    end
    return H
end

走査型トンネル顕微鏡が測定する物理量を計算する

走査型トンネル顕微鏡で測定されるトンネル電流の微分$dI/dV$を考える。これは、飼料の電子の局所状態密度に比例していることが知られており、絶対零度における定義は、

\rho(E,R)=\sum_i^N |\psi_i(R)|^2 \delta(E-E_i)

である。ここで、$i$はハミルトニアンの固有値のインデックス、$\psi_i(R)$は固有値$E_i$に関する固有ベクトルの$R$での値。デルタ関数は数値計算的には難しいため、ローレンツ関数$\frac{1}{\pi}\frac{\eta}{x^2+\eta^2}$で近似すると、

\rho(E,\eta,R)=\sum_i^N |\psi_i(R)|^2\frac{\eta}{(E-E_i)^2+\eta^2}\frac{1}{\pi}

となる。ここで、$\eta$はスメアリングパラメータと呼ばれる。
ソースコードは以下のようになる。

function calc_ldos(E,i,ene,psi,eta)
    ldos = 0.0
    for n = 1:length(ene)
        ldos += abs(psi[i,n])^2*eta/((E-ene[n])^2+eta^2)
    end
    return ldos
end

実行結果は下図になる。

デルタ関数の数値計算

ディラックのデルタ関数$\delta(x)$とは

\int_{-\infty}^{\infty}f(x)\delta(x)dx=f(0)

を満たすシュワルツ超関数のこと。

\delta(x)=
\begin{cases}
0 & x\neq 0\\
\infty & x = 0
\end{cases}

のようにあらわせる。
デルタ関数のような特異性を持つ関数は、数値計算では精度が悪くなりやすく扱いづらい。そのため、別の関数で置き換えることがよくある。

ローレンツ関数

ローレンツ関数はローレンツ分布(コーシー分布)の確率密度関数で、

f(x,x_0,\gamma)=\frac{1}{\pi}\frac{\gamma}{(x-x_0)^2+\gamma^2}

となる。ここで、$x_0$は分布の最頻値を与える位置母数、$\gamma$は半値半幅を与える尺度母数。
ローレンツ分布は期待値・分散が定義されない。
$x_0=0$で$\gamma$を変化させていった時のローレンツ関数の形は下図のようになる。

ローレンツ関数は物理の振動に関して、よくカーブフィッティングに用いられる。
ローレンツモデルにおける複素誘電率の実部(共鳴成分)がローレンツ関数の形をしており、ガウス関数よりもよく近似できるためである。

ガウス関数

ガウス関数は正規分布の確率密度関数で、

f(x,\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{-\frac{(x-\mu)^2}{2\sigma^2}}

となる。ここで、$\mu$は期待値、$\sigma^2$は分散となる。
$\mu=0$で$\sigma$を変化させていった時のガウス関数の形は下図のようになる。

自然界の事象の中には正規分布に従う数量の分布をとるものがあることが知られており、統計分野などでガウス関数がフィッティングによく使われる。

佐藤超関数論

佐藤超関数論は佐藤幹夫氏により考案された一般化関数論である。なお、佐藤幹夫氏は2023/1/9に亡くなっている。
佐藤超関数論では、一般化関数$f(x)$は、複素解析関数$F(z)$の実軸上における境界値の差

f(x)=F(x+i0)-F(x-i0)

で表す。デルタ関数は

\delta(x)=-\frac{1}{2\pi i}\left(\frac{1}{x+i0}-\frac{1}{x-i0} \right)

とできる。
これを用いた数値計算は超函数法とDE公式により行う。

離散フーリエ変換

実空間の波紋の形状をさらに解析するため、フーリエ変換を行う。周期境界条件があるため、離散フーリエ変換をすれば波数空間の情報が得られる。ある周期$N$の関数$f(x)$が$x=0,1,\cdots,N-1$の点上で定義されるとき、離散フーリエ変換は

F(k)=\sum_{x=1}^{N-1}f(x)e^{-i\frac{2\pi}{N}kx}

と定義される。離散フーリエ変換を拘束に行うアルゴリズムとして高速フーリエ変換(FFT)が知られている。
JuliaではFFTW.jlというパッケージで使用できる。使用方法は、

    ldosfft = fft(ldos)
    ldosfft[1,1] = 0
    ldosfftshift = fftshift(ldosfft)
    freq = fftshift(fftfreq(Lx,2*pi))

となる。このパッケージの大元はAbstractFFTs.jlとなっている。今回使用した関数はすべてこちらのもの。

上記計算結果にFFTを実行すると下図が得られる。

フェルミ面を計算する

フェルミ面とは

プランク定数の半整数倍($1/2,3/2,\cdots$)の強度のスピン角運動量を伴う粒子のことをフェルミ粒子と呼び、代表例として電子が挙げられる。フェルミ粒子は一つの系内において、複数の粒子が同一の量子状態を占有することがないというパウリの排他原理に従う。
フェルミ粒子の平均粒子数はフェルミ分布に従い、$T=0$で$\epsilon = \mu$を境とする階段状のステップ関数になる。このとき、フェルミ粒子は$\mu$よりも低い固有エネルギーを有する量子状態がすべて粒子により占有されている。この状態を基底状態と呼ぶ。
$T=0$において、量子状態の占有状態と非占有状態を分けるエネルギーをフェルミ準位と呼ぶ。エネルギーがフェルミ準位に等しい状態の波数を集めると曲面になり、これをフェルミ面と呼ぶ。

フェルミ面を計算する

基底状態での１つの電子のエネルギー$\epsilon=-2\cos(k_xa)-2\cos(k_ya)=\mu$となるので、等高線をプロットするcontourを使用する。

function TightBinding_Fermisurface()
    epsilonk(kx,ky) = -2*(cos(kx)+cos(ky))
    Nkx = 100
    Nky = 100
    kxs = range(-pi,pi,length=Nkx)
    kys = range(-pi,pi,length=Nky)
    myu = 0.2
    energy = zeros(Float64,Nkx,Nky)
    for (ikx,kx) in enumerate(kxs)
        for (iky,ky) in enumerate(kys)
            energy[ikx,iky] = epsilonk(kx,ky)
        end
    end
    contour(kxs,kys,energy,levels=[1,0.5,myu,-0.2,-0.5,-1],aspect_ratio=:equal,xlabel="kx",ylabel="ky",xlims=(-pi,pi),ylims=(-pi,pi),colorbar=false)
    savefig("julia_6-TightiBinding_Fermi.png")
end

結果は下図の通りとなる。

最適化パッケージで計算する

フェルミ面は$\epsilon -\mu =0$となる点の集合です。これを方程式と考え、これを満たす解$\boldsymbol{k}$を探す。最適化計算を行うには、

f(k_x,k_y)=(\epsilon-\mu)^2=(-2\cos(k_xa)-2cos(k_ya)-\mu)^2

という関数を考えればよい。($\epsilon$に関する２次方程式となり、$\epsilon=\mu$のとき最小値をとる。)
この関数の最小値を求めるために、Optimi.jlというパッケージを用いる。

function findzero(myu)
    epsilonk(k) = (-2*(cos(k[1])+cos(k[2]))-myu)^2
    k0 = 2*pi*rand(2).-pi
    res = optimize(epsilonk,k0)
    kans = Optim.minimizer(res)
    for i = 1:2
        while kans[i] > pi
            kans[i] -= 2*pi
        end
        while kans[i] < -pi
            kans[i] += 2*pi
        end
    end
    return kans
end

1000点計算した結果は下図の通り。理論値と似たような形になったことがわかる。

超伝導平均場理論を自己無撞着計算する

超伝導平均場理論

BCS理論とは超電導現象を初めて微視的に説明した理論。この理論より導出されるギャップ方程式を解き、超伝導秩序変数を求める。ギャップ方程式は、

\Delta=NV\int_0^{\hbar\omega_c}d\xi\frac{\Delta\tanh(\frac{\sqrt{\xi^2+\Delta^2}}{2k_bT})}{\sqrt{\xi^2+\Delta^2}}

となる。ここで、$N$はフェルミエネルギーにおける電子の状態密度、$V$は引力相互作用の強さ、$\hbar$はプランク定数、$\omega_c$はデバイ振動数、$k_B$はボルツマン定数。また、温度$T$は変数で、この方程式は

\Delta=G(\Delta,T)

という$\Delta$に関して非線形な方程式となっている。ここで、$\Delta$は超伝導秩序変数と呼ばれる。
左辺と右辺が等しくなるような$\Delta$を求めることを自己無撞着(左辺と右辺に矛盾が無いよう)に解く、という。

ギャップ方程式には、$xi$に関する積分があるので、１次元数値積分ができるパッケージQuadGk.jlを使用する。

function gapequation(delta,T,V)
    f(xi) = N*V*delta*tanh(sqrt(delta^2+xi^2)/(2*kB*T))/(sqrt(delta^2+xi^2))
    result = QuadGK.quadgk(f,0,hw)
    return result[1]
end

最適化パッケージを使用する

$\Delta=G(\Delta,T)$となる$\Delta$を求めるには、関数$f(\Delta)=(G(\Delta,T)-\Delta)^2$の最小値を求める問題に帰着できる。Optim.jlパッケージを使用すると以下のようなコードとなる。

function solve_optim(gapequation,T,V,eps=1e-4,maxdelta=1.5)
    gfunc(delta) = (gapequation(delta,T,V)-delta)^2
    res = optimize(gfunc,1e-10,maxdelta,rel_tol=eps)
    delta = abs(Optim.minimizer(res)[1])
    return delta
end

結果は下図のようになる。

逐次代入法で計算する

方程式は$x=G(x)$という形をしているので、ある適当な値$x_0$を用意し、$x_1=G(x_0),x_2=G(x_1),\cdots$のように次々と代入することで$x_i$がある値に収束する。この方法を逐次代入法という。以下のようなコードとなる。

function solve_direct(gapequation,T,V,eps=1e-4,initialdelta=1.5)
    delta = initialdelta
    deltaold = 100.0
    res = abs((delta-deltaold)/delta)
    while res > eps && abs(delta) > 1e-10
        delta = abs(gapequation(delta,T,V))
        res = abs((delta-deltaold)/delta)
        deltaold = delta
    end
    return delta
end

結果は下図のようになる。

線形化ギャップ方程式を解く

ギャップ方程式で超伝導を考える際、低温域で考えることとなり、$\Delta$は非常に小さい値となる。そのため、$\Delta^2$の項をすべてゼロにすることができる。そのため、ギャップ方程式は以下のように変形できる。

1=NV\int_0^{\hbar\omega_c}d\xi\frac{\tanh(\frac{xi}{2k_BT})}{xi}

これを解くことで理論値とよく一致することが確認できる。

参考