BCS Mean-Field Theory

BCS mean-field theory is the self-consistent quadratic approximation to a fermionic pairing Hamiltonian. It replaces pair scattering among many momentum states by a complex pair field, diagonalizes the resulting particle–hole Hamiltonian with a fermionic Bogoliubov transformation, and determines that field from the state of the quasiparticles it creates.

The logical chain is

$$\text{Fermi surface} \longrightarrow \text{Cooper instability} \longrightarrow \text{pairing saddle} \longrightarrow \text{gapped quasiparticles}.$$

Each arrow carries an assumption. A Fermi surface supplies a large phase space of nearly degenerate pair states. An attractive eigenvalue in a pairing channel produces a logarithmic normal-state instability. Mean-field factorization selects a saddle in that channel. The quadratic saddle then supports Bogoliubov quasiparticles. None of these steps proves that every interacting fermion system is a weak-coupling BCS state.

The exact reduced pairing Hamiltonian conserves particle number. Its familiar mean-field representative does not: it mixes sectors differing by two particles while preserving fermion parity. This distinction is essential in finite systems and whenever the phrase “broken U(1)” is used.

Helpful background. Fermionic operators fixes pair signs and parity; the Fermi surface explains the low-energy shell. The exact reduced BCS model treats finite number sectors.

Time-reversed pairs and the pairing window

We use the grand Hamiltonian

$$K = H-\mu\widehat N,$$

and define

$$\xi_{\mathbf k} = \epsilon_{\mathbf k}-\mu.$$

The elementary model assumes:

• two fermionic internal states denoted $\uparrow$ and $\downarrow$;
• time-reversal- and inversion-symmetric dispersion, so $\xi_{-\mathbf k}=\xi_{\mathbf k}$;
• zero-center-of-mass pairing of $\mathbf k\uparrow$ with $-\mathbf k\downarrow$;
• a real, nonnegative, even form factor $w_{\mathbf k}=w_{-\mathbf k}$ in the elementary model;
• an attractive coupling $g>0$ in the sign convention used below;
• system volume $\mathcal V$;
• a single-spin density of states per unit volume, denoted $\nu(\xi)$.

For the elementary isotropic model,

$$w_{\mathbf k} = \begin{cases} 1, & |\xi_{\mathbf k}|<\omega_{\mathrm c}, \\ 0, & |\xi_{\mathbf k}|>\omega_{\mathrm c}. \end{cases}$$

Here $\omega_{\mathrm c}$ denotes an energy, despite the frequency-like letter; for a phonon frequency $\omega_{\mathrm D}$ the corresponding energy scale is $\hbar\omega_{\mathrm D}$. It is an interaction cutoff, not a universal constant. The shell shorthand is centered on a chosen reference chemical potential. When taking thermodynamic derivatives in $\mu$, hold the microscopic energies, $g$ and the selected form factors fixed. Moving a sharp shell with $\mu$ differentiates the effective Hamiltonian itself and adds boundary terms; neglecting those terms is a separate weak-coupling approximation. In the simplest phonon-window model it is of order a characteristic phonon energy. In a continuum contact model, a bare cutoff and coupling must instead be eliminated in favor of a physical scattering parameter.

Momentum sums run over one label $\mathbf k$ for each ordered pair of spin orbitals $(\mathbf k\uparrow,-\mathbf k\downarrow)$. With this convention, every spin orbital appears once across the product state and quasiparticle sum.

From an Interaction to the Reduced Pairing Model

A general two-body interaction contains many scattering channels. The reduced BCS model retains scattering of time-reversed pairs with zero total momentum. Define

$$b_{\mathbf k}^\dagger = c_{\mathbf k\uparrow}^\dagger c_{-\mathbf k\downarrow}^\dagger, \qquad b_{\mathbf k} = c_{-\mathbf k\downarrow} c_{\mathbf k\uparrow}.$$

The collective pair operator is

$$B^\dagger = \sum_{\mathbf k} w_{\mathbf k}b_{\mathbf k}^\dagger.$$

The reduced grand Hamiltonian is

$$K_{\mathrm{red}} ={} \sum_{\mathbf k,\sigma} \xi_{\mathbf k} c_{\mathbf k\sigma}^\dagger c_{\mathbf k\sigma} - \frac{g}{\mathcal V} B^\dagger B.$$

Equivalently,

$$K_{\mathrm{red}} ={} \sum_{\mathbf k,\sigma} \xi_{\mathbf k} c_{\mathbf k\sigma}^\dagger c_{\mathbf k\sigma} - \frac{g}{\mathcal V} \sum_{\mathbf k,\mathbf k'} w_{\mathbf k}w_{\mathbf k'} b_{\mathbf k}^\dagger b_{\mathbf k'}.$$

The factor $1/\mathcal V$ makes the interaction energy extensive when $g$ is held fixed in the thermodynamic limit.

What has been discarded

The reduction is physical input, not an operator identity. It suppresses or absorbs:

• density and exchange terms;
• finite-center-of-mass pairing;
• retardation of the interaction;
• detailed frequency and momentum dependence;
• Coulomb pseudopotential effects;
• competing spin, charge, and magnetic channels;
• self-energy corrections to the normal-state dispersion.

In a conventional weak-coupling application, some of those effects motivate the attractive low-energy interaction and renormalize $\xi_{\mathbf k}$ and $g$. They are not derived by solving the reduced Hamiltonian.

Exact Number Symmetry

The pair operators carry number charges

$$[ \widehat N, b_{\mathbf k}^\dagger ] = 2b_{\mathbf k}^\dagger, \qquad [ \widehat N, b_{\mathbf k} ] = -2b_{\mathbf k}.$$

The pair-scattering product is neutral:

$$[ \widehat N, b_{\mathbf k}^\dagger b_{\mathbf k'} ] = 0.$$

Therefore

$$[ \widehat N, K_{\mathrm{red}} ] = 0.$$

The exact model moves a pair from one time-reversed orbital pair to another. The exact Hilbert space decomposes into fixed-$N$ sectors, and the reduced Hamiltonian preserves each one.

The Hamiltonian also preserves fermion parity,

$$P_F = (-1)^{\widehat N}.$$

Parity will remain a symmetry even after the usual mean-field approximation mixes sectors differing by two particles.

Pair Operators Are Not Elementary Bosons

For one momentum label,

$$[ b_{\mathbf k}, b_{\mathbf k}^\dagger ] = 1 - n_{\mathbf k\uparrow} - n_{-\mathbf k\downarrow}.$$

The right-hand side is state dependent. Moreover,

$$\left( b_{\mathbf k}^\dagger \right)^2 = 0.$$

One pair orbital can be empty, paired, or blocked by a singly occupied fermionic mode. Pair operators for distinct labels commute because they have even fermion parity, but each label retains a hard-core constraint.

This is why a weak-coupling Cooper pair should not be pictured as a pointlike elementary boson. Its constituents occupy a broad, Pauli-constrained superposition of momentum states.

Pseudospin representation

Define Anderson pseudospins

$$S_{\mathbf k}^+ = b_{\mathbf k}^\dagger, \qquad S_{\mathbf k}^- = b_{\mathbf k},$$

and

$$S_{\mathbf k}^z = \frac{1}{2} \left( n_{\mathbf k\uparrow} + n_{-\mathbf k\downarrow} - 1 \right).$$

On the empty and paired subspace, these obey a spin-$1/2$ algebra. Apart from a constant and blocked levels,

$$K_{\mathrm{red}} = 2\sum_{\mathbf k} \xi_{\mathbf k}S_{\mathbf k}^z - \frac{g}{\mathcal V} \sum_{\mathbf k,\mathbf k'} w_{\mathbf k}w_{\mathbf k'} S_{\mathbf k}^+S_{\mathbf k'}^-.$$

This representation exposes the competition between energy, which favors pseudospins pointing according to the sign of $\xi_{\mathbf k}$, and pairing, which favors transverse alignment. It also connects the reduced model to the exact Richardson solution. Mean field is not the only possible treatment of the reduced Hamiltonian.

Cooper Instability Preview

The Cooper problem asks what happens when two additional fermions interact attractively above an inert filled Fermi sea. Use the trial state

$$|\Psi\rangle = \sum_{\mathbf k}^{0<\xi_{\mathbf k}<\omega_{\mathrm c}} \phi_{\mathbf k} b_{\mathbf k}^\dagger |\mathrm{FS}\rangle.$$

Let the pair energy relative to $2\mu$ be $-E_{\mathrm b}$, with $E_{\mathrm b}>0$. For a constant attraction in the shell, the Schrödinger equation gives

$$\left( 2\xi_{\mathbf k} + E_{\mathrm b} \right) \phi_{\mathbf k} = \frac{g}{\mathcal V} \sum_{\mathbf k'} \phi_{\mathbf k'}.$$

Hence

$$\phi_{\mathbf k} \propto \frac{1}{ 2\xi_{\mathbf k}+E_{\mathrm b} }.$$

The consistency condition is

$$1 = \frac{g}{\mathcal V} \sum_{\mathbf k}^{0<\xi_{\mathbf k}<\omega_{\mathrm c}} \frac{1}{ 2\xi_{\mathbf k}+E_{\mathrm b} }.$$

For an approximately constant single-spin density of states $\nu(0)$,

$$1 = \frac{g\nu(0)}{2} \ln \left( \frac{ 2\omega_{\mathrm c}+E_{\mathrm b} }{ E_{\mathrm b} } \right).$$

Solving gives

$$E_{\mathrm b} = \frac{ 2\omega_{\mathrm c} }{ \exp[ 2/(g\nu(0)) ] - 1 }.$$

At weak coupling,

$$E_{\mathrm b} \simeq 2\omega_{\mathrm c} \exp \left[ - \frac{2}{g\nu(0)} \right].$$

The logarithm diverges as $E_{\mathrm b}\to0^+$. In the idealized continuum problem, every $g>0$ therefore produces a bound pair above the Fermi sea.

What the Cooper problem establishes

It establishes that the normal Fermi sea is singularly susceptible to attraction in the pair channel. A large set of pair states close to the Fermi surface are nearly degenerate, so repeated scattering cannot be represented by a regular finite-order expansion in $g$. Perturbation Theory in Many-Body Systems places this logarithm in the broader map of infrared failures and resummations.

What it does not establish

The Cooper problem:

• freezes the Fermi sea;
• treats only one added pair;
• uses only unoccupied states above the surface;
• does not determine a many-body order parameter;
• does not produce the BCS number equation or thermodynamics.

The exponent differs from the BCS gap exponent because the many-body gap equation uses coherent particle and hole states on both sides of the Fermi surface.

Within the paired shell, the normal branches $\pm\xi$ cross at $\xi=0$, whereas a nonzero pair field gives the Nambu eigenvalues $\pm\sqrt{\xi^2+|\Delta|^2}$. The positive branch has minimum $|\Delta|$; the vertical separation between branches there is $2|\Delta|$. The negative branch represents the particle–hole partner.

Anomalous Mean-Field Decoupling

Write

$$B = \langle B\rangle + \delta B.$$

Then

$$B^\dagger B ={} \langle B^\dagger\rangle B + B^\dagger\langle B\rangle - \langle B^\dagger\rangle \langle B\rangle + \delta B^\dagger\delta B.$$

BCS mean field discards the connected quadratic fluctuation

$$\delta B^\dagger\delta B.$$

Define the pair field

$$\Delta = \frac{g}{\mathcal V} \langle B\rangle = \frac{g}{\mathcal V} \sum_{\mathbf k} w_{\mathbf k} \langle b_{\mathbf k}\rangle.$$

The interaction becomes

$$- \frac{g}{\mathcal V} B^\dagger B \;\longrightarrow\; - \Delta B^\dagger - \Delta^*B + \frac{\mathcal V|\Delta|^2}{g}.$$

Thus

$$K_{\mathrm{MF}} ={} \sum_{\mathbf k,\sigma} \xi_{\mathbf k} c_{\mathbf k\sigma}^\dagger c_{\mathbf k\sigma} - \sum_{\mathbf k} w_{\mathbf k} \left( \Delta b_{\mathbf k}^\dagger + \Delta^*b_{\mathbf k} \right) + \frac{\mathcal V|\Delta|^2}{g}.$$

The final term is the mean-field subtraction constant. Omitting it leaves the quasiparticle spectrum unchanged but gives the wrong grand potential, gap equation, and condensation energy.

Mean field is self-consistent

The quadratic Hamiltonian depends on $\Delta$, but $\Delta$ is not a freely chosen external source. The state or thermal ensemble of $K_{\mathrm{MF}}(\Delta)$ must reproduce the expectation value used to define it:

$$\Delta = \frac{g}{\mathcal V} \sum_{\mathbf k} w_{\mathbf k} \langle c_{-\mathbf k\downarrow} c_{\mathbf k\uparrow} \rangle_{\Delta}.$$

This fixed-point condition is the gap equation.

Relation to Hartree–Fock

The decoupling is anomalous because it uses an expectation value of two annihilation operators. A more general Hartree–Fock–Bogoliubov treatment may retain normal densities and exchange fields as well. In the reduced model, normal-state shifts are usually omitted or absorbed into $\xi_{\mathbf k}$ and $\mu$ so that the pairing structure remains visible.

Fermionic Nambu Hamiltonian

Introduce the Nambu spinor

$$\Psi_{\mathbf k} = \begin{pmatrix} c_{\mathbf k\uparrow} \\ c_{-\mathbf k\downarrow}^\dagger \end{pmatrix}.$$

Then

$$K_{\mathrm{MF}} = \frac{\mathcal V|\Delta|^2}{g} + \sum_{\mathbf k} \left[ \xi_{\mathbf k} + \Psi_{\mathbf k}^\dagger h_{\mathbf k} \Psi_{\mathbf k} \right],$$

where

$$h_{\mathbf k} = \begin{pmatrix} \xi_{\mathbf k} & -w_{\mathbf k}\Delta \\ -w_{\mathbf k}\Delta^* & -\xi_{\mathbf k} \end{pmatrix}.$$

The added $\xi_{\mathbf k}$ is the constant generated when the hole component is reordered:

$$c_{-\mathbf k\downarrow} c_{-\mathbf k\downarrow}^\dagger = 1 - c_{-\mathbf k\downarrow}^\dagger c_{-\mathbf k\downarrow}.$$

Keeping this constant is as important as retaining the mean-field subtraction term.

The Nambu matrix is Hermitian. Its eigenvalues are

$$\lambda_{\mathbf k,\pm} = \pm E_{\mathbf k},$$

with

$$E_{\mathbf k} = \sqrt{ \xi_{\mathbf k}^2 + w_{\mathbf k}^2|\Delta|^2 }.$$

The negative branch is a particle–hole partner in the doubled Nambu description. It is not an independent negative-energy excitation.

Bogoliubov–Valatin Transformation

Choose $u_{\mathbf k}\ge0$ and write the canonical transformation as

$$\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}.$$

Preserving fermionic anticommutators requires

$$u_{\mathbf k}^2 + |v_{\mathbf k}|^2 = 1.$$

This is an ordinary unitary rotation in each Nambu block, unlike the hyperbolic normalization of bosonic Bogoliubov theory.

The coherence factors may be chosen as

$$u_{\mathbf k}^2 = \frac{1}{2} \left( 1 + \frac{\xi_{\mathbf k}}{E_{\mathbf k}} \right),$$

and

$$|v_{\mathbf k}|^2 = \frac{1}{2} \left( 1 - \frac{\xi_{\mathbf k}}{E_{\mathbf k}} \right).$$

Their relative phase is fixed by

$$u_{\mathbf k}v_{\mathbf k} = \frac{ w_{\mathbf k}\Delta }{ 2E_{\mathbf k} }.$$

The diagonal Hamiltonian is

$$K_{\mathrm{MF}} ={} K_{\mathrm{vac}} + \sum_{\mathbf k} E_{\mathbf k} \left( \gamma_{\mathbf k\uparrow}^\dagger \gamma_{\mathbf k\uparrow} + \gamma_{-\mathbf k\downarrow}^\dagger \gamma_{-\mathbf k\downarrow} \right),$$

where

$$K_{\mathrm{vac}} = \frac{\mathcal V|\Delta|^2}{g} + \sum_{\mathbf k} \left( \xi_{\mathbf k} - E_{\mathbf k} \right).$$

Every quasiparticle excitation has nonnegative energy $E_{\mathbf k}$.

Physical Meaning of the Coherence Factors

Far above the Fermi surface,

$$\xi_{\mathbf k} \gg |\Delta|,$$

so

$$u_{\mathbf k}^2 \to 1, \qquad |v_{\mathbf k}|^2 \to 0.$$

The positive-energy quasiparticle is predominantly a particle.

Far below the Fermi surface,

$$\xi_{\mathbf k} \ll -|\Delta|,$$

so

$$u_{\mathbf k}^2 \to 0, \qquad |v_{\mathbf k}|^2 \to 1.$$

The same positive-energy quasiparticle is predominantly a hole.

At the Fermi surface,

$$\xi_{\mathbf k} = 0,$$

and

$$u_{\mathbf k}^2 = |v_{\mathbf k}|^2 = \frac{1}{2}.$$

Particle and hole character are maximally mixed. Pairing turns their normal-state crossing into an avoided crossing.

Bogoliubov Quasiparticles develops the excitation-centered interpretation of these amplitudes, including addition and removal weight, expected microscopic charge, parity, and the contrast with bosonic squeezing.

The BCS Variational State

The quasiparticle vacuum satisfies

$$\gamma_{\mathbf k\uparrow} |\mathrm{BCS}\rangle = 0, \qquad \gamma_{-\mathbf k\downarrow} |\mathrm{BCS}\rangle = 0.$$

It is

$$|\mathrm{BCS}\rangle = \prod_{\mathbf k} \left( u_{\mathbf k} + v_{\mathbf k} b_{\mathbf k}^\dagger \right) |0\rangle.$$

Each factor is normalized because

$$u_{\mathbf k}^2 + |v_{\mathbf k}|^2 = 1.$$

The state is a coherent product of empty and paired configurations in every time-reversed orbital pair. It is not a tensor product of distinguishable real-space molecules.

Where $u_{\mathbf k}\ne0$, it can be written as

$$|\mathrm{BCS}\rangle = \left( \prod_{\mathbf k}u_{\mathbf k} \right) \exp \left[ \sum_{\mathbf k} \frac{v_{\mathbf k}}{u_{\mathbf k}} b_{\mathbf k}^\dagger \right] |0\rangle.$$

The exponential truncates automatically at each $\mathbf k$ because $(b_{\mathbf k}^\dagger)^2=0$.

Occupation and anomalous amplitude

At zero temperature,

$$\langle c_{\mathbf k\uparrow}^\dagger c_{\mathbf k\uparrow} \rangle = |v_{\mathbf k}|^2,$$

and

$$\langle c_{-\mathbf k\downarrow} c_{\mathbf k\uparrow} \rangle = u_{\mathbf k}v_{\mathbf k} = \frac{ w_{\mathbf k}\Delta }{ 2E_{\mathbf k} }.$$

The sharp Fermi step is rounded over an energy scale of order $|\Delta|$. This is a coherent ground-state redistribution, not thermal smearing.

Finite-Temperature Averages

For the thermal state of the quadratic mean-field Hamiltonian,

$$\langle \gamma_{\mathbf k\sigma}^\dagger \gamma_{\mathbf k\sigma} \rangle = f(E_{\mathbf k}),$$

where

$$f(E) = \frac{1}{ e^{\beta E}+1 }.$$

The anomalous average becomes

$$\begin{aligned} \langle c_{-\mathbf k\downarrow} c_{\mathbf k\uparrow} \rangle ={}& u_{\mathbf k}v_{\mathbf k} \left[ 1 - 2f(E_{\mathbf k}) \right] \\ ={}& \frac{ w_{\mathbf k}\Delta }{ 2E_{\mathbf k} } \tanh \left( \frac{\beta E_{\mathbf k}}{2} \right). \end{aligned}$$

Thermally excited quasiparticles therefore reduce the pair amplitude.

The Gap Equation

Insert the anomalous average into the definition of $\Delta$:

$$\Delta = \frac{g}{\mathcal V} \sum_{\mathbf k} \frac{ w_{\mathbf k}^2\Delta }{ 2E_{\mathbf k} } \tanh \left( \frac{\beta E_{\mathbf k}}{2} \right).$$

The normal solution

$$\Delta = 0$$

always exists. For a nonzero solution, one may divide by $\Delta$:

$$1 = \frac{g}{\mathcal V} \sum_{\mathbf k} \frac{ w_{\mathbf k}^2 }{ 2E_{\mathbf k} } \tanh \left( \frac{\beta E_{\mathbf k}}{2} \right).$$

Dividing too early hides the normal branch and can obscure bifurcations or first-order behavior in generalized models.

General momentum-dependent pairing

For an interaction kernel $V_{\mathbf k\mathbf k'}$ with attraction represented by a positive pairing eigenvalue, define

$$\Delta_{\mathbf k} = \frac{1}{\mathcal V} \sum_{\mathbf k'} V_{\mathbf k\mathbf k'} \langle b_{\mathbf k'} \rangle.$$

The self-consistency equation is

$$\Delta_{\mathbf k} = \frac{1}{\mathcal V} \sum_{\mathbf k'} V_{\mathbf k\mathbf k'} \frac{ \Delta_{\mathbf k'} }{ 2E_{\mathbf k'} } \tanh \left( \frac{\beta E_{\mathbf k'}}{2} \right),$$

with

$$E_{\mathbf k} = \sqrt{ \xi_{\mathbf k}^2 + |\Delta_{\mathbf k}|^2 }.$$

This is an eigenvalue problem when linearized near the transition. Different eigenfunctions encode different pairing symmetries. Their microscopic selection and material consequences are beyond the reduced isotropic model.

The Number Equation

The mean occupation of a paired block is

$$\begin{aligned} \langle n_{\mathbf k\uparrow} + n_{-\mathbf k\downarrow} \rangle ={}& 1 - \frac{\xi_{\mathbf k}}{E_{\mathbf k}} \left[ 1 - 2f(E_{\mathbf k}) \right] \\ ={}& 1 - \frac{\xi_{\mathbf k}}{E_{\mathbf k}} \tanh \left( \frac{\beta E_{\mathbf k}}{2} \right). \end{aligned}$$

At fixed density,

$$n = \frac{N}{\mathcal V} = \frac{1}{\mathcal V} \sum_{\mathbf k} \left[ 1 - \frac{\xi_{\mathbf k}}{E_{\mathbf k}} \tanh \left( \frac{\beta E_{\mathbf k}}{2} \right) \right].$$

The gap and number equations must be solved together for $\Delta$ and $\mu$ when $N$ rather than $\mu$ is fixed.

In the normal limit $\Delta\to0$, the bracket reduces to

$$2f(\xi_{\mathbf k}),$$

including both spins. This is a useful sign and factor-of-two check.

Why the number equation matters

In weak-coupling metals with an approximately particle–hole-symmetric band, the shift of $\mu$ can be tiny. That does not make the number equation optional in general. It becomes important for:

• asymmetric densities of states;
• shallow bands;
• finite systems;
• low densities;
• strong pairing;
• multiband models;
• BCS–BEC crossover calculations.

Grand Potential and Stationarity

The mean-field grand potential is

$$\Omega_{\mathrm{MF}}(\Delta,\mu,T) ={} \frac{\mathcal V|\Delta|^2}{g} + \sum_{\mathbf k} \left( \xi_{\mathbf k} - E_{\mathbf k} \right) - \frac{2}{\beta} \sum_{\mathbf k} \ln \left( 1+e^{-\beta E_{\mathbf k}} \right).$$

An equivalent form is

$$\Omega_{\mathrm{MF}} = \frac{\mathcal V|\Delta|^2}{g} + \sum_{\mathbf k} \xi_{\mathbf k} - \frac{2}{\beta} \sum_{\mathbf k} \ln \left[ 2\cosh \left( \frac{\beta E_{\mathbf k}}{2} \right) \right].$$

Stationarity gives

$$\frac{ \partial\Omega_{\mathrm{MF}} }{ \partial\Delta^* } = 0,$$

which reproduces the gap equation. Number follows from

$$N = - \frac{ \partial\Omega_{\mathrm{MF}} }{ \partial\mu }.$$

Take these derivatives at fixed microscopic energies and form factors, as specified above. For the optimized potential the implicit gap derivative vanishes by stationarity. The identities fail if the subtraction or Nambu reordering constant is dropped.

A root is not automatically stable

Solving the gap equation finds stationary points. One must still:

1. compare $\Omega_{\mathrm{MF}}$ among all solutions at fixed $\mu$ and $T$;
2. compare the Helmholtz free energy at fixed $N$;
3. examine the Hessian in amplitude and competing-channel directions;
4. verify that the assumed homogeneous, zero-momentum saddle is not unstable.

For the elementary weak-coupling model, the nonzero solution is the lower grand-potential branch below $T_{\mathrm c}$.

Zero-Temperature Weak-Coupling Gap

Assume:

• isotropic $w_{\mathbf k}=1$ in $|\xi|<\omega_{\mathrm c}$;
• constant single-spin density of states $\nu(0)$;
• $\omega_{\mathrm c}$ small compared with the scale over which $\nu(\xi)$ changes;
• $\mu$ fixed at its normal-state value to leading weak-coupling order.

At $T=0$,

$$\tanh \left( \frac{\beta E}{2} \right) \to 1.$$

The gap equation becomes

$$1 = g\nu(0) \int_0^{\omega_{\mathrm c}} \frac{ d\xi }{ \sqrt{ \xi^2+\Delta_0^2 } }.$$

Therefore

$$\frac{1}{g\nu(0)} = \operatorname{arsinh} \left( \frac{\omega_{\mathrm c}}{\Delta_0} \right).$$

Within the constant-density cutoff model,

$$\Delta_0 = \frac{ \omega_{\mathrm c} }{ \sinh[ 1/(g\nu(0)) ] }.$$

For

$$g\nu(0) \ll 1,$$

this reduces to

$$\Delta_0 \simeq 2\omega_{\mathrm c} \exp \left[ - \frac{1}{g\nu(0)} \right].$$

The gap is nonanalytic at $g=0$. No finite power series in the attractive coupling can produce the exponential scale.

Linearized Gap Equation and Transition Temperature

At a continuous transition, let $\Delta\to0$. Then

$$E_{\mathbf k} \to |\xi_{\mathbf k}|.$$

The linearized equation is

$$1 = \frac{g}{\mathcal V} \sum_{\mathbf k} \frac{ w_{\mathbf k}^2 }{ 2|\xi_{\mathbf k}| } \tanh \left( \frac{ |\xi_{\mathbf k}| }{ 2k_{\mathrm B}T_{\mathrm c} } \right).$$

For the constant-density cutoff model,

$$1 = g\nu(0) \int_0^{\omega_{\mathrm c}} \frac{d\xi}{\xi} \tanh \left( \frac{\xi}{2k_{\mathrm B}T_{\mathrm c}} \right).$$

In weak coupling,

$$k_{\mathrm B}T_{\mathrm c} \simeq \frac{ 2e^{\gamma} }{ \pi } \omega_{\mathrm c} \exp \left[ - \frac{1}{g\nu(0)} \right],$$

where $\gamma$ is the Euler–Mascheroni constant.

Combining this with the zero-temperature gap gives

$$\frac{ 2\Delta_0 }{ k_{\mathrm B}T_{\mathrm c} } = \frac{2\pi}{e^\gamma} \simeq 3.528.$$

This ratio is universal only within the weak-coupling, isotropic, constant-density BCS assumptions. Strong coupling, anisotropy, multiple bands, pair breaking, and unconventional nodes can change it.

Near $T_{\mathrm c}$, the same model gives

$$\Delta(T) \simeq 3.06 k_{\mathrm B}T_{\mathrm c} \sqrt{ 1-\frac{T}{T_{\mathrm c}} }.$$

This square-root behavior is the mean-field critical law. Critical fluctuations can modify the asymptotic behavior sufficiently close to a real transition.

Pair Susceptibility View

Apply an external pair source $K^\prime=-\eta B^\dagger-\eta^*B$ and define the physical response by $\delta\langle B\rangle/\mathcal V=\chi_{\mathrm{pair}}\eta$ in the normal state. This fixes the sign and volume normalization; susceptibilities develops the general complex-source convention.

Define the normal-state static pair susceptibility for the chosen form factor:

$$\chi_{\mathrm{pair}}(T) = \frac{1}{\mathcal V} \sum_{\mathbf k} \frac{ w_{\mathbf k}^2 }{ 2\xi_{\mathbf k} } \tanh \left( \frac{\beta\xi_{\mathbf k}}{2} \right).$$

The factor $\tanh(\beta\xi/2)/(2\xi)$ is even and has the finite limit $\beta/4$ at $\xi=0$. A nonsymmetric density of states need not make the entire momentum integral particle–hole symmetric. The transition criterion is

$$1 = g\chi_{\mathrm{pair}}(T_{\mathrm c}).$$

For a Fermi surface with nonzero $\nu(0)$,

$$\chi_{\mathrm{pair}}(T) \sim \nu(0) \ln \left( \frac{\omega_{\mathrm c}}{k_{\mathrm B}T} \right)$$

at low temperature in the idealized normal state. This logarithm is the many-body form of the Cooper instability.

The criterion identifies a normal-state instability. The nonlinear gap equation determines the ordered saddle below it. Collective response and fluctuation corrections require a conserving treatment beyond the static saddle.

Condensation Energy

At $T=0$, compare the paired and normal grand potentials at the same $\mu$. For the elementary constant-density model,

$$\frac{ \Omega_{\mathrm s}-\Omega_{\mathrm n} }{ \mathcal V } ={} \frac{\Delta_0^2}{g} - 2\nu(0) \int_0^{\omega_{\mathrm c}} d\xi \left[ \sqrt{ \xi^2+\Delta_0^2 } - \xi \right].$$

Using the gap equation and expanding at weak coupling gives

$$\frac{ \Omega_{\mathrm s}-\Omega_{\mathrm n} }{ \mathcal V } \simeq - \frac{1}{2} \nu(0)\Delta_0^2.$$

Here $\nu(0)$ is the single-spin density of states. If a two-spin density of states is used instead, the prefactor must be adjusted. Many apparent factor-of-two disagreements come from changing this convention silently.

The negative sign shows that the paired saddle is energetically favored. The gain comes from coherent mixing near the Fermi surface and is balanced against the cost $\mathcal V|\Delta|^2/g$.

What the Excitation Gap Means

For an isotropic paired shell containing a continuum Fermi surface, in the weak-coupling regime $|\Delta|\ll\omega_{\mathrm c}$,

$$\min_{\mathbf k}E_{\mathbf k} = |\Delta|.$$

This is the minimum energy of one mean-field quasiparticle in the grand-canonical description. Outside the sharp pairing window the model has $w=0$ and energies $|\xi|\ge\omega_{\mathrm c}$; if one extrapolates to $|\Delta|>\omega_{\mathrm c}$, the global gap is no longer $|\Delta|$. Finite discrete levels need not contain a level exactly at $\xi=0$ either.

An isolated even-parity system cannot change fermion parity by a number-conserving perturbation. The simplest pair-breaking excitation creates two quasiparticles, with threshold

$$2|\Delta|$$

in the ideal clean isotropic model.

These are different statements:

• $|\Delta|$ is the one-quasiparticle spectral gap;
• $2|\Delta|$ is the elementary two-quasiparticle threshold under parity-preserving excitation;
• collective modes need not begin at either value;
• anisotropic $\Delta_{\mathbf k}$ can have smaller minima or nodes;
• disorder, interactions, and lifetime effects can alter observed thresholds.

The order parameter and an experimentally inferred spectral gap coincide only under the assumptions of the simple model.

Number Fluctuations in the BCS State

The unprojected variational state is not an eigenstate of $\widehat N$. Its mean number is

$$\langle\widehat N\rangle = 2\sum_{\mathbf k} |v_{\mathbf k}|^2.$$

Because each pair block contains either zero or two particles,

$$\operatorname{Var}(\widehat N) = 4\sum_{\mathbf k} u_{\mathbf k}^2 |v_{\mathbf k}|^2.$$

The variance is extensive in the number of active pair orbitals, while the relative fluctuation scales to zero as

$$\frac{ \sqrt{ \operatorname{Var}(\widehat N) } }{ \langle\widehat N\rangle } \sim \mathcal V^{-1/2}$$

under ordinary thermodynamic scaling.

Vanishing relative fluctuations help explain why the unprojected state gives accurate bulk observables. They do not turn it into an exact number eigenstate.

Number projection

A fixed-$N$ state can be formed with

$$\mathcal P_N = \frac{1}{2\pi} \int_0^{2\pi} d\theta\, e^{i\theta(\widehat N-N)}.$$

Then

$$|\mathrm{BCS};N\rangle \propto \mathcal P_N |\mathrm{BCS}\rangle.$$

Projection onto a sector with nonzero weight restores sharp particle number while retaining pair correlations. This even-parity product has no odd-$N$ component; an odd system requires a blocked unpaired orbital. Variational Many-Body States places number-projected and Gutzwiller-projected BCS states alongside determinant, Jastrow, MPS, and neural ansätze. In finite systems, projected variation and exact Richardson methods can distinguish effects hidden by the thermodynamic saddle.

Number rotations and the phase of a paired state

Use the active state rotation $|\psi_\theta\rangle=U(\theta)|\psi\rangle$ with $U(\theta)=e^{-i\theta\widehat N}$. The operator conjugation obeys

$$U c U^\dagger=e^{i\theta}c,\qquad U b U^\dagger=e^{2i\theta}b.$$

The expectation in the rotated state instead contains $U^\dagger bU$:

$$\langle b\rangle_{\psi_\theta} =e^{-2i\theta}\langle b\rangle_\psi, \qquad \Delta_\theta=e^{-2i\theta}\Delta.$$

Accordingly,

$$U K_{\mathrm{MF}}(\Delta)U^\dagger =K_{\mathrm{MF}}(e^{-2i\theta}\Delta).$$

Keeping the state transformation and operator conjugation distinct avoids a phase-sign error. The exact reduced Hamiltonian commutes with $\widehat N$, while a fixed nonzero mean-field representative does not; choosing one phase selects a branch of its covariant family.

Finite systems

In an exact number eigenstate,

$$\langle b_{\mathbf k} \rangle = 0$$

by the number selection rule. Pair order can instead be diagnosed through number-neutral quantities, such as the pair density matrix

$$M_{\mathbf k\mathbf k'} = \langle b_{\mathbf k}^\dagger b_{\mathbf k'} \rangle.$$

A macroscopic eigenvalue of this matrix is a number-conserving signature of pair condensation. Off-diagonal long-range order explains pair-space normalization and distinguishes condensation from pair binding and the excitation gap.

Thermodynamic order of limits

A symmetry-breaking construction may add an infinitesimal source

$$K_\eta = K_{\mathrm{red}} - \eta B^\dagger - \eta^*B,$$

take $\mathcal V\to\infty$, and only then let $\eta\to0$. At fixed $T>0$, removing the source first restores the symmetric finite-volume Gibbs state. At $T=0$, the same conclusion holds when the finite-volume ground state is unique; an exact ground-state degeneracy requires separate analysis. The general symmetry-breaking treatment distinguishes these conditions from bulk order (Koma and Tasaki, 1994, §1.1).

For an explicit finite check, retain one pair orbital with $w=1$, set $\lambda=g/\mathcal V>0$ and $\xi=\lambda/2$, and take a real source $\eta>0$. Then

$$K_\eta=\frac{\lambda}{2}(n_\uparrow+n_\downarrow) -\lambda n_\uparrow n_\downarrow-\eta(b^\dagger+b), \qquad b=c_\downarrow c_\uparrow.$$

At zero source, $|0\rangle$ and $|2\rangle=b^\dagger|0\rangle$ are degenerate ground states, while each singly occupied state costs $\lambda/2$. In the even sector the source gives the matrix $-\eta\sigma_x$. Its ground state and anomalous average are

$$|\psi_\eta\rangle=\frac{|0\rangle+|2\rangle}{\sqrt2}, \qquad \langle b\rangle_{\psi_\eta}=\frac12.$$

This average survives $\eta\to0^+$ at $T=0$ without increasing the system size. It reflects selection within an isolated finite degeneracy. It does not establish a thermodynamic ordered phase, and a number eigenstate or the symmetric zero-source mixture still has $\langle b\rangle=0$.

What remains unbroken

Pair terms change particle number by two, so fermion parity remains:

$$[ P_F, K_{\mathrm{MF}} ] = 0.$$

The continuous number symmetry is reduced to its $\mathbb Z_2$ subgroup in the mean-field representative.

Global symmetry and electromagnetic gauge redundancy

For a neutral fermionic superfluid, the broken global number U(1) language directly identifies a phase mode. For a charged superconductor, electromagnetic gauge transformations are redundancies of description, and physical statements must be gauge invariant. The mean-field phase, electromagnetic field, and vertex corrections must be treated together.

It is therefore safer to say that BCS mean field uses a particle-number-symmetry-breaking representation than to say that a gauge redundancy itself is an observable broken symmetry. Random Phase Approximation explains the generic self-consistent response structure and why a paired saddle requires matching phase vertices. Collective Modes supplies the general phase–amplitude, response-matrix, hybridization, and damping framework; London theory develops the local Meissner response, while a superconducting collective-mode derivation remains beyond this static treatment.

Phase and Amplitude Fluctuations

Write a slowly varying pair field as

$$\Delta(\mathbf r,t) = \left[ \Delta_0 + h(\mathbf r,t) \right] e^{i\phi(\mathbf r,t)}.$$

The fields have different roles:

• $h$ changes the magnitude of pairing;
• $\phi$ changes its phase;
• in a neutral system, long-wavelength phase fluctuations contain the Goldstone sound mode;
• in a charged system, coupling to electromagnetism reorganizes the phase response;
• the amplitude response often lies near the pair-breaking continuum and need not be a sharp mode.

The static gap equation fixes the saddle but not the full propagators of $h$ and $\phi$. Those require Gaussian fluctuations or response theory with the conservation-law constraints kept intact.

Variational Derivation

The BCS state can also be used directly as a variational ansatz. For real $u_{\mathbf k}$ and normalized factors,

$$\mathcal E = \langle\mathrm{BCS}| K_{\mathrm{red}} |\mathrm{BCS}\rangle.$$

At leading thermodynamic order,

$$\mathcal E ={} 2\sum_{\mathbf k} \xi_{\mathbf k} |v_{\mathbf k}|^2 - \frac{g}{\mathcal V} \left| \sum_{\mathbf k} w_{\mathbf k} u_{\mathbf k}v_{\mathbf k} \right|^2.$$

For the displayed Hamiltonian including the diagonal $\mathbf k=\mathbf k'$ terms, the exact expectation adds

$$-\frac{g}{\mathcal V}\sum_{\mathbf k}w_{\mathbf k}^2|v_{\mathbf k}|^4.$$

Indeed, $\langle b^\dagger_{\mathbf k}b_{\mathbf k}\rangle=|v_{\mathbf k}|^2$, whereas its anomalous factorization is $u_{\mathbf k}^2|v_{\mathbf k}|^2$. With a bounded form factor and an extensive number of active levels this correction is order one, so it does not change the leading bulk energy density. It matters in a finite-level variational calculation and cannot be silently discarded while claiming the exact expectation.

Parameterize

$$u_{\mathbf k} = \cos \left( \frac{\theta_{\mathbf k}}{2} \right),$$

and

$$v_{\mathbf k} = e^{i\varphi} \sin \left( \frac{\theta_{\mathbf k}}{2} \right).$$

Variation aligns each pseudospin with its effective field and gives

$$\cos\theta_{\mathbf k} = \frac{\xi_{\mathbf k}}{E_{\mathbf k}},$$

and

$$\sin\theta_{\mathbf k} = \frac{ w_{\mathbf k}|\Delta| }{ E_{\mathbf k} }.$$

These are exactly the coherence-factor relations. After the stated subextensive correction is neglected, variation and anomalous decoupling give the same leading thermodynamic saddle.

Self-Consistency as a Pseudospin Alignment

In the mean-field Hamiltonian, each pseudospin sees

$$\mathbf h_{\mathbf k} = \left( - w_{\mathbf k}\operatorname{Re}\Delta, \; w_{\mathbf k}\operatorname{Im}\Delta, \; \xi_{\mathbf k} \right).$$

At zero temperature, the ground-state pseudospin aligns opposite to the appropriate Hamiltonian field. Deep below the Fermi surface it is nearly paired; high above it is nearly empty; near $\xi_{\mathbf k}=0$ it tilts into the transverse plane.

Self-consistency requires the transverse components of all pseudospins to regenerate the same $\Delta$. This picture makes three facts immediate:

1. only states near the Fermi surface rotate strongly;
2. all active pair amplitudes share a phase in the simple separable model;
3. the gap equation is collective rather than a separate two-body equation for each pair.

Controlled Regime

The elementary BCS mean-field solution is most reliable when:

• the dimensionless attraction $g\nu(0)$ is small;
• the system has many active states within the pairing shell;
• the coherence length is large compared with microscopic spacing;
• order-parameter fluctuations are weak;
• competing channels are separated in scale;
• disorder and pair breaking are absent or perturbative;
• the normal state has reasonably long-lived fermionic excitations;
• the pairing interaction and cutoff are defined consistently.

In the conventional weak-coupling regime, the small ratio

$$\frac{\Delta_0}{E_{\mathrm F}} \ll 1$$

implies that many overlapping pairs contribute to a smooth collective field.

What Mean Field Omits

Pair-field fluctuations

Amplitude and phase fluctuations are discarded at the saddle. They matter near critical points, in low dimensions, and in small or disordered systems.

Quasiparticle interactions

The diagonal Hamiltonian treats quasiparticles as independent. Residual interactions determine lifetimes, collective modes, and corrections to response.

Retardation and strong coupling

A static constant $g$ does not reproduce the frequency dependence of a microscopic boson-mediated interaction. Strong-coupling superconductivity requires a frequency-dependent self-energy and pairing kernel rather than the elementary BCS equations.

Competing orders

Magnetism, density waves, nematicity, and other pair channels can alter or preempt the simple saddle. Choosing one decoupling channel in advance can miss the actual minimum.

Exact finite-size structure

Level spacing, parity effects, blocked orbitals, and number projection become important when the number of active levels is not large. The exact reduced Hamiltonian remains meaningful even where spontaneous-symmetry-breaking language is not.

Non-BCS normal states

The standard logarithm assumes a Fermi surface with suitable coherent low-energy fermions. Strongly incoherent or non-Fermi-liquid normal states require a more careful pairing analysis.

Exercises

Normal-state limit of the number equation

Show that

$$1 - \frac{\xi}{E} \tanh \left( \frac{\beta E}{2} \right)$$

reduces to $2f(\xi)$ as $\Delta\to0$.

Solution

When $\Delta\to0$,

$$E \to |\xi|.$$

For $\xi>0$, the expression becomes

$$1 - \tanh \left( \frac{\beta\xi}{2} \right) = \frac{2}{ e^{\beta\xi}+1 } = 2f(\xi).$$

For $\xi<0$, use $|\xi|=-\xi$:

$$1 + \tanh \left( \frac{\beta|\xi|}{2} \right) = \frac{2}{ e^{\beta\xi}+1 } = 2f(\xi).$$

The two cases agree and include the two spin occupations in one paired block.

Number variance

For

$$|\mathrm{BCS}\rangle = \prod_{\mathbf k} \left( u_{\mathbf k} + v_{\mathbf k}b_{\mathbf k}^\dagger \right) |0\rangle,$$

derive $\langle\widehat N\rangle$ and $\operatorname{Var}(\widehat N)$.

Solution

For one $\mathbf k$ block, the particle number is $0$ with probability $u_{\mathbf k}^2$ and $2$ with probability $|v_{\mathbf k}|^2$. Therefore

$$\langle N_{\mathbf k}\rangle = 2|v_{\mathbf k}|^2,$$

and

$$\langle N_{\mathbf k}^2\rangle = 4|v_{\mathbf k}|^2.$$

Its variance is

$$\begin{aligned} \operatorname{Var}(N_{\mathbf k}) &= 4|v_{\mathbf k}|^2 - 4|v_{\mathbf k}|^4 \\ &= 4u_{\mathbf k}^2 |v_{\mathbf k}|^2. \end{aligned}$$

Distinct product factors are statistically independent, so the means and variances add:

$$\langle\widehat N\rangle = 2\sum_{\mathbf k}|v_{\mathbf k}|^2,$$ $$\operatorname{Var}(\widehat N) = 4\sum_{\mathbf k} u_{\mathbf k}^2|v_{\mathbf k}|^2.$$

Condensation-energy check

At $T=0$, start from

$$\frac{ \Omega_{\mathrm s}-\Omega_{\mathrm n} }{ \mathcal V } = \frac{\Delta_0^2}{g} - 2\nu(0) \int_0^{\omega_{\mathrm c}} d\xi \left( \sqrt{ \xi^2+\Delta_0^2 } - \xi \right).$$

Show that its weak-coupling limit is $-\nu(0)\Delta_0^2/2$.

Solution

The integral is

$$I ={} \frac{1}{2} \omega_{\mathrm c} \sqrt{ \omega_{\mathrm c}^2+\Delta_0^2 } - \frac{\omega_{\mathrm c}^2}{2} + \frac{\Delta_0^2}{2} \operatorname{arsinh} \left( \frac{\omega_{\mathrm c}}{\Delta_0} \right).$$

The first line simplifies at weak coupling:

$$\frac{1}{2} \omega_{\mathrm c} \sqrt{ \omega_{\mathrm c}^2+\Delta_0^2 } - \frac{\omega_{\mathrm c}^2}{2} \simeq \frac{\Delta_0^2}{4}.$$

The gap equation gives

$$\frac{1}{g} = \nu(0) \operatorname{arsinh} \left( \frac{\omega_{\mathrm c}}{\Delta_0} \right).$$

Substitute both results:

$$\begin{aligned} \frac{ \Omega_{\mathrm s}-\Omega_{\mathrm n} }{ \mathcal V } &\simeq \nu(0)\Delta_0^2 \operatorname{arsinh} \left( \frac{\omega_{\mathrm c}}{\Delta_0} \right) - 2\nu(0) \left[ \frac{\Delta_0^2}{4} + \frac{\Delta_0^2}{2} \operatorname{arsinh} \left( \frac{\omega_{\mathrm c}}{\Delta_0} \right) \right] \\ &= - \frac{1}{2} \nu(0)\Delta_0^2. \end{aligned}$$

The logarithmic terms cancel only when the subtraction constant and gap equation are both retained.

A finite-level pairing correction

For one active pair orbital, $|\psi\rangle=u|0\rangle+v|2\rangle$ with real $u$ and $u^2+|v|^2=1$, compare the exact interaction expectation $-\lambda\langle b^\dagger b\rangle$ with $-\lambda|\langle b\rangle|^2$.

Solution

The exact value is $-\lambda|v|^2$, while anomalous factorization gives $-\lambda u^2|v|^2$. Their difference is $-\lambda|v|^4$. For a fully occupied pair, $\langle b\rangle=0$ even though the exact interaction energy is $-\lambda$. An anomalous average alone does not reproduce every finite-level diagonal contribution.

Finite-temperature source removal

For the one-pair model $K_\eta$ above, compute $\langle b\rangle$ in the full four-state Gibbs ensemble. Compare removing $\eta$ at fixed $T>0$ with taking $T\to0$ first.

Solution

The even eigenstates have energies $-\eta$ and $+\eta$, with $\langle b\rangle=+1/2$ and $-1/2$. The two odd states have energy $\lambda/2$ and zero anomalous average. Therefore

$$Z=2\cosh(\beta\eta)+2e^{-\beta\lambda/2}, \qquad \langle b\rangle_{\beta,\eta} =\frac{\sinh(\beta\eta)} {2\cosh(\beta\eta)+2e^{-\beta\lambda/2}}.$$

At every finite $\beta$, this tends to zero as $\eta\to0$. At fixed $\eta>0$, taking $\beta\to\infty$ instead gives $1/2$, which remains after source removal. These noncommuting temperature and source limits concern the stated finite crossing; no thermodynamic limit has been taken.

References

• P. W. Anderson, “Random-Phase Approximation in the Theory of Superconductivity,” Physical Review 112, 1900–1916 (1958), doi:10.1103/PhysRev.112.1900.
• J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of Superconductivity,” Physical Review 108, 1175–1204 (1957), doi:10.1103/PhysRev.108.1175.
• N. N. Bogoliubov, “A New Method in the Theory of Superconductivity. I,” Soviet Physics JETP 7, 41–46 (1958), official JETP archive.
• N. N. Bogoljubov, “On a New Method in the Theory of Superconductivity,” Il Nuovo Cimento 7, 794–805 (1958), doi:10.1007/BF02745585.
• L. N. Cooper, “Bound Electron Pairs in a Degenerate Fermi Gas,” Physical Review 104, 1189–1190 (1956), doi:10.1103/PhysRev.104.1189.
• P. G. de Gennes, Superconductivity of Metals and Alloys (Westview Press, 1999).
• A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (Dover, 2003).
• T. Koma and H. Tasaki, “Symmetry Breaking and Finite Size Effects in Quantum Many-Body Systems,” Journal of Statistical Physics 76, 745–803 (1994), author manuscript.
• Y. Nambu, “Quasi-Particles and Gauge Invariance in the Theory of Superconductivity,” Physical Review 117, 648–663 (1960), doi:10.1103/PhysRev.117.648.
• R. W. Richardson, “A Restricted Class of Exact Eigenstates of the Pairing-Force Hamiltonian,” Physics Letters 3, 277–279 (1963), doi:10.1016/0031-9163(63)90259-2.
• J. R. Schrieffer, Theory of Superconductivity, revised ed. (Perseus Books, 1999).
• M. Tinkham, Introduction to Superconductivity, 2nd ed. (Dover, 2004).
• J. G. Valatin, “Comments on the Theory of Superconductivity,” Il Nuovo Cimento 7, 843–857 (1958), doi:10.1007/BF02745589.