BCS Gap Equation and the Superconducting Gap

Derivation

BCS Gap Equation and the Superconducting Gap

D-262 Home PU-303 Threads energy · matter · symmetry · chance Depends on Cooper Pair Instability, Quantization of Lattice Vibrations

Statement

Starting from the reduced BCS pairing Hamiltonian with an attractive interaction \(-V\) between time-reversed electron pairs \((\mathbf{k}\uparrow,-\mathbf{k}\downarrow)\), a mean-field decoupling and a Bogoliubov canonical transformation diagonalise the Hamiltonian in terms of quasiparticles of energy \(E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+\Delta^{2}}\), yielding the self-consistent gap equation \(1=V N(0)\int_{0}^{\hbar\omega_D}\frac{\tanh(\beta E/2)}{2E}\,d\xi\), whose \(T=0\) solution is \(\Delta_0=2\hbar\omega_D e^{-1/N(0)V}\) and whose vanishing-gap limit fixes the critical temperature \(k_B T_c=1.13\,\hbar\omega_D e^{-1/N(0)V}\).

Why it matters

The gap equation is the quantitative heart of superconductivity. It converts the qualitative Cooper-pair instability into a full thermodynamic order parameter \(\Delta(T)\), predicts the energy gap seen in tunnelling and infrared absorption, and yields the celebrated universal ratio \(2\Delta_0/k_BT_c\approx 3.53\) that holds across weak-coupling superconductors independent of material detail.

It is also the prototype of every mean-field theory of a fermionic condensate: the same Bogoliubov-de Gennes structure reappears in superfluid \(^3\)He, neutron-star pairing, colour superconductivity in dense quark matter, and the BCS side of the BCS-BEC crossover in ultracold atoms. Mastering this derivation is mastering the template for spontaneously broken \(U(1)\) symmetry in a fermion system.

Assumptions

Reduced (BCS) interaction.Only zero-total-momentum, opposite-spin pairs \((\mathbf{k}\uparrow,-\mathbf{k}\downarrow)\) are retained; drop this and the mean field is not a single complex scalar \(\Delta\) but a full momentum-dependent gap function requiring the Eliashberg treatment.

Separable, constant attraction.\(V_{\mathbf{k}\mathbf{k}'}=-V\) for \(|\xi|,|\xi'|<\hbar\omega_D\) and zero otherwise; drop the cutoff and the \(\xi\)-integral diverges logarithmically with no gap scale set.

Constant density of states \(N(0)\).The DOS is taken flat over the thin Debye shell around \(E_F\); drop it and \(N(\xi)\) must stay inside the integral, distorting the \(T_c\) prefactor.

Mean-field factorisation.Fluctuations of \(b_{\mathbf{k}}=c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\) about \(\langle b_{\mathbf{k}}\rangle\) are neglected to quadratic order; drop this in low dimensions and phase fluctuations (Berezinskii-Kosterlitz-Thouless physics) destroy the mean-field \(T_c\).

Real, \(\mathbf{k}\)-independent order parameter.The overall \(U(1)\) phase is gauged to zero and \(\Delta\) taken isotropic (s-wave); drop isotropy and one must solve a matrix gap equation with nodes, as in d-wave cuprates.

Derivation

\[ \hat H=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}\,c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}-V\sum_{\mathbf{k}\mathbf{k}'}c_{\mathbf{k}\uparrow}^{\dagger}c_{-\mathbf{k}\downarrow}^{\dagger}c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow} \]

Reduced pairing Hamiltonian, energies measured from \(E_F\) so \(\xi_{\mathbf{k}}=\epsilon_{\mathbf{k}}-\mu\); scatters a pair from \(\mathbf{k}'\) to \(\mathbf{k}\). A

\[ \Delta\equiv V\sum_{\mathbf{k}'}\langle c_{-\mathbf{k}'\downarrow}c_{\mathbf{k}'\uparrow}\rangle,\qquad c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}=\langle\cdots\rangle+\delta b_{\mathbf{k}} \]

Define the gap as the interaction-weighted pair amplitude and split each pair operator into its expectation plus a fluctuation \(\delta b_{\mathbf{k}}\). B

\[ \hat H_{MF}=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}-\sum_{\mathbf{k}}\left(\Delta\,c_{\mathbf{k}\uparrow}^{\dagger}c_{-\mathbf{k}\downarrow}^{\dagger}+\Delta^{*}c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\right)+\frac{|\Delta|^{2}}{V} \]

Insert step 2 and discard the \(O(\delta b^{2})\) term; the constant \(|\Delta|^2/V\) restores the double-counted mean-field energy. C

\[ \begin{pmatrix}\gamma_{\mathbf{k}\uparrow}\\ \gamma_{-\mathbf{k}\downarrow}^{\dagger}\end{pmatrix}=\begin{pmatrix}u_{\mathbf{k}} & -v_{\mathbf{k}}\\ v_{\mathbf{k}} & u_{\mathbf{k}}\end{pmatrix}\begin{pmatrix}c_{\mathbf{k}\uparrow}\\ c_{-\mathbf{k}\downarrow}^{\dagger}\end{pmatrix},\quad u_{\mathbf{k}}^{2}+v_{\mathbf{k}}^{2}=1 \]

Bogoliubov canonical transformation with real \(u,v\); the constraint keeps \(\{\gamma,\gamma^{\dagger}\}=1\), i.e. the new operators are fermions. B

\[ 2u_{\mathbf{k}}v_{\mathbf{k}}\,\xi_{\mathbf{k}}=\Delta\,(u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2}) \]

Require the anomalous (off-diagonal) \(\gamma\gamma\) and \(\gamma^{\dagger}\gamma^{\dagger}\) terms to cancel in \(\hat H_{MF}\); this single condition fixes the mixing angle. C

\[ u_{\mathbf{k}}^{2}=\frac{1}{2}\left(1+\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\right),\quad v_{\mathbf{k}}^{2}=\frac{1}{2}\left(1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\right),\quad E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+\Delta^{2}} \]

Solve step 5 with \(u^2+v^2=1\); the quasiparticle spectrum \(E_{\mathbf{k}}\) shows the minimum gap \(\Delta\) at the Fermi surface \(\xi=0\). B

\[ \hat H_{MF}=E_{0}+\sum_{\mathbf{k}\sigma}E_{\mathbf{k}}\,\gamma_{\mathbf{k}\sigma}^{\dagger}\gamma_{\mathbf{k}\sigma},\qquad E_{0}=\sum_{\mathbf{k}}(\xi_{\mathbf{k}}-E_{\mathbf{k}})+\frac{|\Delta|^{2}}{V} \]

Substituting steps 4-6 diagonalises the Hamiltonian: independent quasiparticles above a condensate ground-state energy \(E_0\). C

\[ \langle c_{-\mathbf{k}\downarrow}c_{\mathbf{k}\uparrow}\rangle=u_{\mathbf{k}}v_{\mathbf{k}}\left(1-2f(E_{\mathbf{k}})\right)=\frac{\Delta}{2E_{\mathbf{k}}}\tanh\!\left(\frac{\beta E_{\mathbf{k}}}{2}\right) \]

Evaluate the pair amplitude with \(\langle\gamma^{\dagger}\gamma\rangle=f(E)\) the Fermi function and \(2u_{\mathbf{k}}v_{\mathbf{k}}=\Delta/E_{\mathbf{k}}\); \(1-2f=\tanh(\beta E/2)\). B

\[ \Delta=V\sum_{\mathbf{k}}\frac{\Delta}{2E_{\mathbf{k}}}\tanh\!\left(\frac{\beta E_{\mathbf{k}}}{2}\right)\;\Longrightarrow\;1=V\sum_{\mathbf{k}}\frac{\tanh(\beta E_{\mathbf{k}}/2)}{2E_{\mathbf{k}}} \]

Reinsert step 8 into the definition of \(\Delta\) (step 2) and cancel the common non-zero \(\Delta\); this is the self-consistency condition. A

\[ 1=N(0)V\int_{0}^{\hbar\omega_{D}}\frac{\tanh\!\left(\sqrt{\xi^{2}+\Delta^{2}}/2k_BT\right)}{\sqrt{\xi^{2}+\Delta^{2}}}\,d\xi \]

Convert \(\sum_{\mathbf{k}}\to N(0)\int d\xi\) over the Debye shell (constant DOS) and use the even integrand to fold \(\int_{-\hbar\omega_D}^{\hbar\omega_D}\to 2\int_0^{\hbar\omega_D}\). B

\[ 1=N(0)V\int_{0}^{\hbar\omega_{D}}\frac{\tanh\!\left(\dfrac{\sqrt{\xi^{2}+\Delta^{2}}}{2k_BT}\right)}{\sqrt{\xi^{2}+\Delta^{2}}}\,d\xi \]

Reading. The gap \(\Delta(T)\) is fixed self-consistently: an electron pair feels an attraction \(V\) but only pairs where a coherent gap \(\Delta\) already exists, weighted by how many quasiparticle states are thermally excited via \(\tanh\). At \(T=0\) the \(\tanh\to1\) and the integral gives \(\Delta_0=2\hbar\omega_D\,e^{-1/N(0)V}\); as \(T\to T_c\) the gap vanishes and \(k_BT_c=1.13\,\hbar\omega_D\,e^{-1/N(0)V}\).

Units check. \(N(0)\) has units \([\text{energy}]^{-1}\) (states per unit energy per spin), \(V\) is \([\text{energy}]\), so \(N(0)V\) is dimensionless — correct for an exponent. The integrand \(d\xi/\sqrt{\xi^2+\Delta^2}\) is dimensionless (energy over energy), so the right side is dimensionless, matching the \(1\) on the left. \(\Delta_0\sim\hbar\omega_D\) carries energy units as required.

Limiting cases

\(T\to0\): \(\tanh\to1\), giving \(\Delta_0=2\hbar\omega_D e^{-1/N(0)V}\) — the maximum gap.
\(T\to T_c^-\): set \(\Delta=0\) in the integrand, yielding \(k_BT_c=(2e^{\gamma}/\pi)\hbar\omega_D e^{-1/N(0)V}\approx1.134\,\hbar\omega_D e^{-1/N(0)V}\) with \(\gamma\) the Euler-Mascheroni constant.
Weak coupling \(N(0)V\ll1\): \(\Delta_0/E_F\ll1\); pairing is a thin Fermi-surface effect and the exponential makes \(T_c\) non-analytic in \(V\) (no perturbation theory reaches it).
Universal ratio: dividing the two limits, \(2\Delta_0/k_BT_c=2\pi/e^{\gamma}\approx3.528\), independent of \(\omega_D\), \(N(0)\), and \(V\).
Near \(T_c\): \(\Delta(T)\approx3.06\,k_BT_c\sqrt{1-T/T_c}\), the mean-field square-root onset of the order parameter.

Breaks when

Strong coupling \(\lambda=N(0)V\gtrsim0.5\): retardation and the finite phonon lifetime matter; the instantaneous \(-V\) is replaced by the frequency-dependent Eliashberg kernel and \(2\Delta_0/k_BT_c\) rises above 3.53 (e.g. Pb, Hg).
Low dimensions (1D/2D): the neglected phase fluctuations diverge; long-range order is destroyed by the Mermin-Wagner/BKT mechanism and the true transition is not the mean-field \(T_c\).
Anisotropic / nodal pairing: for d-wave (cuprates) or p-wave (\(^3\)He) the scalar \(\Delta\) becomes \(\Delta_{\mathbf{k}}\) with sign changes and nodes; the flat-DOS s-wave gap equation gives the wrong \(T_c\) and thermodynamics.
Magnetic impurities / pair breaking: spin-flip scattering breaks time-reversed pairs (Abrikosov-Gor'kov); gapless superconductivity appears and \(T_c\) is suppressed faster than the clean equation predicts.

Failure modes

Cancelling \(\Delta\) from both sides of step 9 while forgetting that \(\Delta=0\) is always a (normal-state) solution — the non-trivial branch requires \(N(0)V>0\).
Writing the quasiparticle energy as \(E_{\mathbf{k}}=\xi_{\mathbf{k}}+\Delta\) instead of \(\sqrt{\xi_{\mathbf{k}}^2+\Delta^2}\); the gap is the minimum of \(E\), not an additive shift.
Integrating \(\xi\) from \(0\) to \(\infty\); the Debye cutoff \(\hbar\omega_D\) is essential — without it the log diverges and no gap scale emerges.
Using \(\hbar\omega_D\) as the \(T=0\) prefactor for both \(\Delta_0\) and \(k_BT_c\) with the same numerical coefficient — the coefficients differ (\(2\) vs \(1.13\)), which is the whole content of the universal ratio.
Confusing the single-particle gap \(\Delta\) with the condensation energy \(\tfrac12N(0)\Delta_0^2\); they scale differently in \(V\).
Treating \(f(E_{\mathbf{k}})\) as the electron occupation rather than the quasiparticle occupation, giving the wrong sign in \(1-2f\).

Discussion

The physical content of the gap equation is spontaneous breaking of the global \(U(1)\) gauge symmetry associated with particle-number conservation. The mean field \(\Delta=|\Delta|e^{i\varphi}\) picks a phase \(\varphi\), and the ground state is a coherent superposition of states with different particle number — the BCS wavefunction \(\prod_{\mathbf{k}}(u_{\mathbf{k}}+v_{\mathbf{k}}c_{\mathbf{k}\uparrow}^{\dagger}c_{-\mathbf{k}\downarrow}^{\dagger})|0\rangle\). The rigidity of \(\varphi\) is precisely what carries dissipationless current, tying this derivation directly to the Meissner effect and flux quantisation.

The exponential \(e^{-1/N(0)V}\) is the derivation's most striking feature: it is non-analytic at \(V=0\), so no finite order of perturbation theory in the electron-phonon coupling ever produces superconductivity. This is why the phenomenon resisted explanation for nearly fifty years after its 1911 discovery — one must sum an infinite ladder of scattering (the Cooper channel) to see the instability that this mean field then stabilises. The prior Cooper-pair result supplies exactly that log-divergent pair susceptibility; the gap equation is its non-perturbative resummation.

The quasiparticles \(\gamma\) are coherent superpositions of an electron and a hole (\(u\) electron-like, \(v\) hole-like), so charge is not a good quantum number for excitations — only \(\pm1\) modulo the condensate. This particle-hole mixing is what tunnelling spectroscopy measures as the symmetric gap edges at \(\pm\Delta\), and it is the direct microscopic ancestor of Bogoliubov-de Gennes quasiparticles and, in the topological context, of Majorana zero modes where \(u=v\).

At the deepest level the gap function is an anomalous propagator \(F(\mathbf{k},\omega)=\langle T\,c_{\mathbf{k}\uparrow}c_{-\mathbf{k}\downarrow}\rangle\) that is nonzero only in the broken-symmetry phase; the mean-field gap equation is the leading (Hartree-Fock-Gor'kov) truncation of the Dyson equation for the \(2\times2\) Nambu Green's function. Promoting \(-V\) to the retarded phonon propagator \(D(\mathbf{k}-\mathbf{k}',\omega-\omega')\) and keeping the electron self-energy gives the Eliashberg equations, which reduce to BCS in the limit \(\omega_D\ll E_F\) with the coupling \(\lambda\) small — placing this entire result as one controlled limit of a fully retarded, strong-coupling field theory.

Common misconceptions. The gap is not a rigid band gap like a semiconductor's: it is a self-consistent, temperature-dependent order parameter that vanishes continuously at \(T_c\), and its very existence is the presence of the condensate, not a property of the bare band structure. Nor does \(\Delta\) equal the pairing binding energy of a single Cooper pair; \(\Delta\) is a collective mean field to which every pair on the Fermi surface contributes.

Worked examples

Estimate \(T_c\) for aluminium: \(\Theta_D=428\,\text{K}\), \(N(0)V=0.18\).

\[ \hbar\omega_D=k_B\Theta_D,\qquad k_BT_c=1.13\,k_B\Theta_D\,e^{-1/N(0)V} \]

Use the \(T_c\) limit with the Debye temperature as the phonon cutoff. A

\[ T_c=1.13\times428\,\text{K}\times e^{-1/0.18}=1.13\times428\times e^{-5.556} \]

\[ T_c=1.13\times428\times3.86\times10^{-3}\,\text{K}=1.87\,\text{K} \]

\[ T_c\approx1.9\,\text{K} \]

Reading. The measured \(T_c\) of Al is \(1.18\,\text{K}\); the weak-coupling estimate is the right order of magnitude, the overshoot reflecting the crudeness of \(N(0)V\) and Coulomb repulsion (the \(\mu^*\) correction) which the bare BCS formula omits.

Units check. \(k_B\Theta_D\) carries energy; the exponential is dimensionless; dividing through by \(k_B\) returns a temperature in kelvin.

Given \(T_c=9.3\,\text{K}\) for niobium, predict the \(T=0\) gap \(\Delta_0\) and the tunnelling edge in meV.

\[ \frac{2\Delta_0}{k_BT_c}=3.53\;\Longrightarrow\;\Delta_0=1.76\,k_BT_c \]

Apply the universal weak-coupling ratio, independent of material parameters. A

\[ \Delta_0=1.76\times(8.617\times10^{-5}\,\text{eV/K})\times9.3\,\text{K} \]

\[ \Delta_0=1.76\times8.01\times10^{-4}\,\text{eV}=1.41\times10^{-3}\,\text{eV} \]

\[ \Delta_0\approx1.41\,\text{meV},\qquad 2\Delta_0\approx2.8\,\text{meV} \]

Reading. A superconductor-insulator-superconductor tunnel junction of Nb shows a conductance gap edge at \(2\Delta_0\approx2.8\,\text{meV}\); the measured value is closer to \(3.0\,\text{meV}\) because Nb is a moderate-coupling material with \(2\Delta_0/k_BT_c\approx3.8\), slightly above the weak-coupling 3.53.

Units check. \(k_B\) in eV/K times \(T\) in K gives eV; scaling by \(10^3\) gives meV, the natural scale for superconducting gaps.

Problems

Show that the constant term \(|\Delta|^2/V\) in \(\hat H_{MF}\) equals \(N(0)\Delta^2\times(1/N(0)V)\) and comment on why it is needed.
Solution
With \(\Delta=V\sum_{\mathbf{k}}\langle b_{\mathbf{k}}\rangle\), the mean-field factorisation \(b^{\dagger}b\to b^{\dagger}\langle b\rangle+\langle b^{\dagger}\rangle b-\langle b^{\dagger}\rangle\langle b\rangle\) leaves the subtraction \(-\sum\langle b^{\dagger}\rangle\langle b\rangle=-|\Delta|^2/V\) with sign flipped by the leading minus of the interaction, giving \(+|\Delta|^2/V\). Writing \(1/V=N(0)/(N(0)V)\) gives \(|\Delta|^2/V=N(0)\Delta^2/(N(0)V)\). It is needed to avoid double-counting the interaction energy already contained in the linear \(\Delta\) terms; without it the ground-state energy \(E_0\) and hence the condensation energy would be wrong.
Derive the condensation energy \(E_s-E_n=-\tfrac12N(0)\Delta_0^2\) at \(T=0\).
Solution
The ground-state energy is \(E_0=\sum_{\mathbf{k}}(\xi_{\mathbf{k}}-E_{\mathbf{k}})+\Delta^2/V\). Converting to \(N(0)\int_{-\hbar\omega_D}^{\hbar\omega_D}(\xi-\sqrt{\xi^2+\Delta^2})\,d\xi+N(0)\Delta^2/(N(0)V)\), evaluating the integral for \(\Delta\ll\hbar\omega_D\) gives \(-\tfrac12N(0)\Delta^2[\,1-\ln(\cdots)\,]\); using the gap equation \(1/N(0)V=\ln(2\hbar\omega_D/\Delta_0)\) to cancel the log term leaves \(E_s-E_n=-\tfrac12N(0)\Delta_0^2\). Numerically for Al (\(N(0)\sim10^{47}\,\text{J}^{-1}\text{m}^{-3}\), \(\Delta_0\sim0.18\,\text{meV}\)) this is a tiny \(\sim1\,\mu\text{eV per atom}\), consistent with the fragility of superconductivity.
A superconductor has \(\Theta_D=275\,\text{K}\) and \(N(0)V=0.25\). Compute \(T_c\) and \(\Delta_0\).
Solution
\(T_c=1.13\times275\times e^{-1/0.25}=1.13\times275\times e^{-4}=1.13\times275\times0.0183=5.69\,\text{K}\). Then \(\Delta_0=1.76\,k_BT_c=1.76\times8.617\times10^{-5}\times5.69\,\text{eV}=8.6\times10^{-4}\,\text{eV}=0.86\,\text{meV}\).
Verify the universal ratio: from \(\Delta_0=2\hbar\omega_D e^{-1/N(0)V}\) and \(k_BT_c=(2e^{\gamma}/\pi)\hbar\omega_D e^{-1/N(0)V}\), obtain \(2\Delta_0/k_BT_c\) numerically (\(\gamma=0.5772\)).
Solution
\(\dfrac{2\Delta_0}{k_BT_c}=\dfrac{2\cdot2\hbar\omega_D e^{-1/N(0)V}}{(2e^{\gamma}/\pi)\hbar\omega_D e^{-1/N(0)V}}=\dfrac{4}{2e^{\gamma}/\pi}=\dfrac{2\pi}{e^{\gamma}}\). With \(e^{0.5772}=1.781\), this is \(2\pi/1.781=6.283/1.781=3.528\). The exponentials and all material parameters cancel, which is the origin of the universality.
Near \(T_c\), expand the gap equation to show \(\Delta(T)\propto\sqrt{1-T/T_c}\) and find the coefficient.
Solution
Subtract the \(T_c\) equation (\(\Delta=0\)) from the finite-\(\Delta\) equation and expand the integrand to \(O(\Delta^2)\). The temperature-derivative term gives \(\ln(T_c/T)\approx1-T/T_c\) on the left, while the \(\Delta^2\) correction produces the coefficient \(\frac{7\zeta(3)}{8\pi^2}(\Delta/k_BT_c)^2\) from \(\sum_n(2n+1)^{-3}=\tfrac{7}{8}\zeta(3)\). Equating gives \(\Delta(T)=k_BT_c\sqrt{\dfrac{8\pi^2}{7\zeta(3)}}\sqrt{1-T/T_c}\). With \(\zeta(3)=1.202\), the prefactor is \(\sqrt{8\pi^2/(7\times1.202)}=\sqrt{9.38}=3.06\), so \(\Delta(T)\approx3.06\,k_BT_c\sqrt{1-T/T_c}\), the classic mean-field square-root onset.