PHYSICS

D. N. ZUBAREV and Yu. A. TSERKOVNIKOV

THERMODYNAMICS OF SUPERCONDUCTORS

(Presented by Academician N. N. Bogolyubov, 9 VII 1958)

In papers \((^1,^2)\) the thermodynamics of superconductors was considered on the basis of a model Hamiltonian in which the electron-phonon interaction is replaced by a direct interaction between electrons with opposite momenta and spins. In this case it proved possible to construct an asymptotically exact solution \((^2)\). In the present paper we shall consider the thermodynamics of superconductors by means of the Fröhlich Hamiltonian, in which the electron-phonon interaction is taken explicitly into account. In doing so we shall use thermodynamic perturbation theory. Following N. N. Bogolyubov’s method, we shall also carry out a renormalization of the electron-phonon interaction constant, which will make it possible to improve the convergence of the expansions, i.e. will give the same advantages as in the case of zero temperatures \((^3)\). The initial Hamiltonian has the form

\[ H=\sum_{k,\sigma}\varepsilon_k a^+_{k\sigma}a_{k\sigma} +\sum_q \omega_q\left(b_q^+b_q+\frac12\right) +\sum_{k'-k=q,\sigma} g_q\sqrt{\frac{\omega_q}{2V}}\,a^+_{k\sigma}a_{k'\sigma}(b_q^+ + b_{-q}), \tag{1} \]

where \(\varepsilon_k=k^2/2m-\mu\); \(\mu\) is the chemical potential; \(g_q\) is the electron-phonon interaction constant; \(a^+_{k\sigma}, a_{k\sigma}, b_q^+\) and \(b_q\) are the creation and annihilation operators, respectively, of fermions and phonons. Following the papers of N. N. Bogolyubov \((^3,^4)\), we perform the canonical transformation of operators:

\[ \begin{gathered} a_{k,\,1/2}=u_k\alpha_{k0}+v_k\alpha^+_{k1},\qquad a_{-k,\,-1/2}=u_k\alpha_{k1}-v_k\alpha^+_{k0},\\ b_q=\lambda_q\beta_q+\mu_q\beta^+_{-q},\qquad u_k^2+v_k^2=1,\qquad \lambda_q^2-\mu_q^2=1, \end{gathered} \tag{2} \]

where \(u_k, v_k, \lambda_k\) and \(\mu_q\) are real parameters depending on \(|k|\) and \(|q|\).

The Hamiltonian (1) in the new operators is written in the form:

\[ H=H_0+H',\qquad H'=H_1+H_2+H_3+H_4+H_{\mathrm I}+H_{\mathrm{II}}, \tag{3} \]

where

\[ H_0=U+\sum_k \tilde{\varepsilon}_k(\alpha^+_{k0}\alpha_{k0}+\alpha^+_{k1}\alpha_{k1}) +\sum_q \tilde{\omega}_q\beta^+_q\beta_q,\qquad U=2\sum_k \varepsilon_k v_k^2+\frac12\sum_q\omega_q(\lambda_p^2+\mu_q^2); \tag{4} \]

\[ H_1=\sum_k\{\varepsilon_k(u_k^2-v_k^2)-\tilde{\varepsilon}_k\} (\alpha^+_{k0}\alpha_{k0}+\alpha^+_{k1}\alpha_{k1}),\qquad H_2=\sum_q\{\omega_q(\lambda_q^2+\mu_q^2)-\tilde{\omega}_q\}\beta^+_q\beta_q, \]

\[ H_3=\sum_k 2\varepsilon_k u_kv_k(\alpha^+_{k0}\alpha^+_{k1}+\alpha_{k1}\alpha_{k0}),\qquad H_4=\sum_q \omega_q\lambda_q\mu_q(\beta^+_q\beta^+_{-q}+\beta_{-q}\beta_q), \]

\[ H_{\mathrm I}= \sum_{k'-k=q} g_{\mathrm I}(k,k')(\beta_q+\beta^+_{-q}) (\alpha^+_{k'0}\alpha_{k0}+\alpha^+_{k1}\alpha_{k'1}), \tag{5} \]

\[ H_{\mathrm{II}}= \sum_{k'-k=q} g_{\mathrm{II}}(k,k')(\beta_q+\beta^+_{-q}) (\alpha^+_{k'0}\alpha^+_{k1}+\alpha_{k'1}\alpha_{k0}); \]

\[ \begin{aligned} g_{\mathrm I}(kk')&=g_q(\omega_q/2V)^{1/2}(\lambda_q+\mu_q)(u_{k'}u_k-v_{k'}v_k),\\ g_{\mathrm{II}}(kk')&=g_q(\omega_q/2V)^{1/2}(\lambda_q+\mu_q)(u_{k'}v_k+v_{k'}u_k). \end{aligned} \tag{6} \]

In the Hamiltonian (4), the quantities \(\widetilde{\varepsilon}_k\) and \(\widetilde{\omega}_q\) are “renormalized” energies of the fermion and boson excitations, which will be determined below.

To calculate the thermodynamic potential \(\Omega\), we shall use thermodynamic perturbation theory in the form given in [5]:

\[ \Omega=\Omega_0+\langle R\rangle_c, \tag{7} \]

where

\[ \langle R\rangle_c=-\frac{1}{\beta}\sum_{n\geqslant 1}\frac{(-1)^n}{n!} \int_0^\beta \cdots \int_0^\beta dt_1\ldots dt_n \langle T H'(t_1)\ldots H'(t_n)\rangle_c, \]

\[ H'(t)=e^{H_0t}H'e^{-H_0t},\qquad \beta=1/\theta. \tag{8} \]

The index \(c\) in formula (8) means that only connected diagrams are taken into account. The averaging is performed over the grand Gibbs ensemble with Hamiltonian \(H_0\). After integration over \(t_1,\ldots,t_n\), (8) takes the form

\[ \langle R\rangle_c=\sum_{n\geqslant 1}(1)^{n-1} \left\langle H'\frac{1}{E_1}H'\frac{1}{E_2}\cdots \frac{1}{E_{n-1}}H'\right\rangle_c, \tag{9} \]

where \(E_i\) are energy denominators (see [5]).

In the problem under consideration, the unperturbed thermodynamic potential is

\[ \Omega_0=U-2\theta\sum_k \ln\{1+e^{-\widetilde{\varepsilon}_k/\theta}\} +\theta\sum_q \ln\{1-e^{-\widetilde{\omega}_q/\theta}\}. \tag{10} \]

Substituting the operator \(H'\) ((3), (5)) into (9), it is not difficult to see that the operators \(H_3\) and \(H_4\) already in the second approximation lead to the appearance of “dangerous” denominators \(1/2\widetilde{\varepsilon}_k\) and \(1/2\widetilde{\omega}_q\) (see, for example, the diagrams in Fig. 1a), which lead to divergences upon integration over \(k\) and \(q\). The parameters \(u_k, v_k\) and \(\lambda_q,\mu_q\) must be chosen so that these dangerous terms do not appear. For this purpose, for example, in the case of the parameters \(u_k\) and \(v_k\) it is sufficient to set equal to zero, in the second and third approximations, the sum of the contributions to the thermodynamic potential from all diagrams with two lines entering (or leaving) the vertex corresponding to the operator \(H_3\) (Fig. 1a). Also excluding the terms with dangerous denominators arising from the operator \(H_4\), and choosing the elementary excitations \(\widetilde{\varepsilon}_k\) and \(\widetilde{\omega}_q\) so that \(H_1\) and \(H_2\) do not give a contribution with the factor \(\beta=1/\theta\) corresponding to the “zero” denominator in the Bloch series (9) (see below), in the second approximation we obtain

\[ \Omega=\Omega_0+\sum_k\{\varepsilon_k(u_k^2-v_k^2)-\widetilde{\varepsilon}_k\}\,2n_k +\sum_q\{\omega_q(\lambda_q^2+\mu_q^2)-\widetilde{\omega}_q\}\nu_q+ \]

\[ +\sum_{k'-k=q}\{g_{\mathrm I}^2(k,k')J_{\mathrm I}(k,k') +g_{\mathrm{II}}^2(k,k')J_{\mathrm{II}}(k,k')\}; \tag{11} \]

\[ J_{\mathrm I}(k,k')= \frac{(1+\nu_q)(1-n_k)n_{k'}}{-\widetilde{\omega}_q-\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'}} +\frac{\nu_q(1-n_k)n_{k'}}{\widetilde{\omega}_q-\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'}} +\frac{(1+\nu_q)n_k(1-n_{k'})}{-\widetilde{\omega}_q+\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'}} +\frac{\nu_q n_k(1-n_{k'})}{\widetilde{\omega}_q+\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'}}, \]

\[ J_{\mathrm{II}}(k,k')= \frac{(1+\nu_q)(1-n_k)(1-n_{k'})}{-\widetilde{\omega}_q-\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'}} +\frac{\nu_q(1-n_k)(1-n_{k'})}{\widetilde{\omega}_q-\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'}} +\frac{(1+\nu_q)n_kn_{k'}}{-\widetilde{\omega}_q+\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'}} + \]

\[ +\frac{\nu_q n_kn_{k'}}{\widetilde{\omega}_q+\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'}}, \qquad n_k=\{e^{\widetilde{\varepsilon}_k/\theta}+1\}^{-1}, \qquad \nu_q=\{e^{\widetilde{\omega}_q/\theta}-1\}^{-1}. \tag{12} \]

The equation obtained as a result of the compensation of the diagrams in Fig. 1a,

\[ \frac{1}{4}\frac{\delta \Omega}{\delta s_k} = -\varepsilon_k u_k v_k(1-2n_k)+ \]

\[ +\sum_{k'\,(k'-k=q)} g_{\mathrm I}(k,k')g_{\mathrm{II}}(k,k') \{J_{\mathrm I}(k,k')-J_{\mathrm{II}}(k',k')\}=0 \tag{13} \]

coincides with the equation obtained by setting equal to zero the variation of the thermodynamic potential (11) with respect to \(u_k\) and \(v_k\) \((\delta u_k=v_k\delta s_k,\ \delta v_k=-u_k\delta s_k)\).

Fig. 1

Arguing in an analogous manner, for \(\lambda_q\) and \(\mu_q\) we shall have \((\delta \lambda_q=\mu_q\delta s_q,\ \delta \mu_q=-\lambda_q\delta s_q)\):

\[ \frac{1}{2}\frac{\delta \Omega}{\delta s_q} = \omega_q\lambda_q\mu_q(1+2\nu_q)+ \]

\[ +\sum_{k'\,(k'-k=q)} \{g_{\mathrm I}^{\,2}(k,k')J_{\mathrm I}(k,k') + g_{\mathrm{II}}^{\,2}(k,k')J_{\mathrm{II}}(k,k')\}=0. \tag{14} \]

The energies of the elementary excitations of fermions and bosons are determined by excluding diagrams with insertions of the operators \(H_1\) and \(H_2\) at the upper (or lower) vertex of the diagram (see, for example, Fig. 1b). Such compensation ensures the disappearance, in the perturbation-theory series, of terms proportional to \(\beta\), and is equivalent to the conditions

\[ \frac{\delta}{\delta n_k}\langle R\rangle_C=0, \qquad \frac{\delta}{\delta \nu_q}\langle R\rangle_C=0. \tag{15} \]

Let us now write equations (13), (14) in a more convenient form. Put:

\[ C_k(1-2n_k) = \frac{1}{V} \sum_{k'\,(k'-k=q)} g_q^2\omega_q(\lambda_q+\mu_q)^2u_{k'}v_{k'} \{J_{\mathrm I}(k,k')-J_{\mathrm{II}}(k,k')\}, \]

\[ S_k(1-2n_k) = \frac{1}{2V} \sum_{k'\,(k'-k=q)} g_q^2\omega_q(\lambda_q+\mu_q)^2 (u_{k'}^2-v_{k'}^2) \{J_{\mathrm I}(k,k')-J_{\mathrm{II}}(k,k')\}, \]

\[ P_q(1+2\nu_q) = -\frac{2}{V} \sum_{k'\,(k'-k=q)} g_q^2 \{(u_ku_{k'}-v_kv_{k'})^2J_{\mathrm I}(k,k') + (u_kv_{k'}+u_{k'}v_k)^2J_{\mathrm{II}}(k,k')\}. \]

Then, as is not difficult to verify, we shall have \((\xi_k=\varepsilon_k-S_k)\)

\[ u_k^2=\frac{1}{2} \left\{1+\frac{\xi_k}{\sqrt{C_k^2+\xi_k^2}}\right\}, \qquad v_k^2=\frac{1}{2} \left\{1-\frac{\xi_k}{\sqrt{C_k^2+\xi_k^2}}\right\}, \qquad (\lambda_q+\mu_q)^2=(1-P_q)^{-1/2}. \tag{16} \]

The functions \(C_k\) and \(S_k\) are determined by the equations

\[ (1-2n_k)C_k=\frac{1}{2V}\sum_{k'(k'-k=q)} g_q^2\omega_q(\lambda_q+\mu_q)^2 \bigl(J_{\mathrm I}(k,k')-J_{\mathrm{II}}(k,k')\bigr) \frac{C_{k'}}{\sqrt{C_{k'}^2+\xi_{k'}^2}}, \tag{17} \]

\[ (1-2n_k)S_k=\frac{1}{2V}\sum_{k'(k'-k=q)} g_q^2\omega_q(\lambda_q+\mu_q)^2 \bigl(J_{\mathrm I}(k,k')-J_{\mathrm{II}}(k,k')\bigr) \frac{\xi_{k'}}{\sqrt{C_{k'}^2+\xi_{k'}^2}}. \tag{18} \]

Expanding equations (15), we find

\[ \widetilde{\varepsilon}_k \left\{ 1+2\sum_{k'(k'-k=q)} g_{\mathrm I}^2 \frac{ (\widetilde{\omega}_q^2-\widetilde{\varepsilon}_k^2+\widetilde{\varepsilon}_{k'}^2)(1+2\nu_q) -2\widetilde{\omega}_q\widetilde{\varepsilon}_k(1-2n_{k'}) }{ [\widetilde{\omega}_q^2-(\widetilde{\varepsilon}_q+\widetilde{\varepsilon}_{k'})^2] [\widetilde{\omega}_q^2-(\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'})^2] } \right\} = \]

\[ = \left\{ \varepsilon_k-\sum_{k'(k'-k=q)} g_q^2\frac{\omega_q(\lambda_q+\mu_q)^2}{2V} \frac{\widetilde{\omega}_q(1-2n_{k'}) -(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})(1+2\nu_q)} {\widetilde{\omega}_q^2-(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})^2} \,(u_{k'}^2-v_{k'}^2) \right\} \times \]

\[ \times (u_k^2-v_k^2) +4u_kv_k\sum_{k'(k'-k=q)} g_q^2\frac{\omega_q(\lambda_q+\mu_q)^2}{2V} \frac{\widetilde{\omega}_q(1-2n_{k'}) -(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})(1+2\nu_q)} {\widetilde{\omega}_q^2-(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})^2} u_{k'}v_{k'}, \tag{19} \]

\[ \widetilde{\omega}_q = \omega_q(\lambda_q^2+\mu_q^2) + \sum_{k(k'-k=q)} 2g_{\mathrm{II}}^2 \frac{(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})(1-2n_{k'})} {\omega_q^2-(\widetilde{\varepsilon}_k+\widetilde{\varepsilon}_{k'})^2} + \sum_{k(k'-k=q)} 4g_{\mathrm I}^2 \frac{(\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'})n_{k'}} {\omega_q^2-(\widetilde{\varepsilon}_k-\widetilde{\varepsilon}_{k'})^2}. \]

Equations (17) and (18), together with the expressions for the elementary excitations (19) and the condition for determining the chemical potential, completely determine the coefficients of the canonical transformations (2). For \(\theta=0\) we obtain the results of work \({}^{3}\). As usual, these equations have a “normal” solution \(C_k=0\) \((u_kv_k=0)\). The nontrivial solution determines the superconducting state of the system.

V. A. Moskalenko \({}^{6}\) considered the thermodynamics of the superconducting state for the Fröhlich Hamiltonian by means of an approximate variational method without taking into account the renormalization of the phonons; for the free energy he obtained an expression analogous to relation (12).

We now proceed to the calculation of thermodynamic functions. Taking into account equations (13)—(15), we see that, in calculating the entropy from the thermodynamic potential (12), differentiation with respect to temperature may be performed only with respect to the explicitly entering parameter \(\theta\) \((d\Omega/d\theta=\partial\Omega_0/\partial\theta)\). Calculating the heat capacity, for the jump of the heat capacity at the critical point we obtain (taking into account the dependence on \(\theta\) through \(C_k\))

\[ \Delta C = -\left\{ \frac{1}{\theta} \left[ \sum_k n_k(1-n_k)\frac{\partial\widetilde{\varepsilon}_k^2}{\partial\theta} + \frac{1}{2}\sum_q \nu_q(1+\nu_q)\frac{\partial\widetilde{\omega}_q^2}{\partial\theta} \right] \right\}_{\theta=\theta_{\mathrm{cr}}}, \tag{20} \]

where \(\theta=\theta_{\mathrm{cr}}\) is the critical temperature at which the nontrivial solution \((u_kv_k\ne0)\) of the system of equations (17), (18) and the energy gap in the spectrum \(\widetilde{\varepsilon}_k\) (19) appear. The first term in formula (20) coincides with the expression for the jump of the heat capacity obtained in \({}^{1,2}\), where as \(C_k\) one must take the solution of equation (17).

In conclusion, we express our deep gratitude to Academician N. N. Bogolyubov for valuable advice, to V. A. Moskalenko for discussion of the work, and to K. Bloch for sending a preprint.

Mathematical Institute named after V. A. Steklov
Academy of Sciences of the USSR

Received
20 VI 1958

CITED LITERATURE

\({}^{1}\) J. Bardeen, L. N. Cooper, T. R. Schrieffer, Phys. Rev., 108, 1175 (1957).
\({}^{2}\) N. N. Bogolyubov, D. N. Zubarev, Yu. A. Tserkovnikov, DAN, 117, 788 (1957).
\({}^{3}\) N. N. Bogolyubov, Scientific Reports of Higher School, No. 1 (1958).
\({}^{4}\) N. N. Bogolyubov, ZhETF, 34, 58 (1958).
\({}^{5}\) C. Bloch, C. de Dominicis, Nuch. Phys., 7, 459 (1958).
\({}^{6}\) V. A. Moskalenko, DAN, 123, No. 3 (1958).