# Migdal-Eliashberg Theory as a Classical Spin Chain

Emil A. Yuzbashyan  Department of Physics and Astronomy, Rutgers University,
Piscataway, NJ 08854, USA  Boris L. Altshuler  Physics Department, Columbia
University, 538 West 120th Street, New York, NY 10027, USA

###### Abstract

We formulate the Migdal-Eliashberg theory of electron-phonon interactions in
terms of classical spins by mapping the free energy to a Heisenberg spin chain
in a Zeeman magnetic field. Spin components are energy-integrated normal and
anomalous Green’s functions and sites of the chain are fermionic Matsubara
frequencies. The Zeeman field grows linearly with the spin coordinate and
competes with ferromagnetic spin-spin interaction that falls off as the square
of the inverse distance. The spin chain representation makes a range of
previously unknown properties plain to see. In particular, infinitely many new
solutions of the Eliashberg equations both in the normal and superconducting
states emerge at strong coupling. These saddle points of the free energy
functional correspond to spin flips. We argue that they are also fixed points
of kinetic equations and play an essential role in far from equilibrium
dynamics of strongly coupled superconductors. Up to an overall phase, the
frequency-dependent gap function that minimizes the free energy must be
nonnegative. There are strong parallels between our Eliashberg spins and
Anderson pseudospins, though the two sets of spins never coincide.

##  I Introduction

Electron-phonon interactions determine many properties of quantum metals, such
as the charge and heat transport, thermodynamics, and superconductivity. A
well-established approach to these phenomena is the Migdal-Eliashberg theory
[1, 2]. This theory operates with two functions of Matsubara frequency,
$\Sigma_{n}\equiv\Sigma(\omega_{n})$ and $\Phi_{n}\equiv\Phi(\omega_{n})$,
which have the meaning of the normal and anomalous self-energies. The
functions $\Sigma_{n}$ and $\Phi_{n}$ must satisfy two nonlinear algebraic
equations known as the Eliashberg equations, whose kernel is the phonon-
mediated electron-electron interaction $\lambda_{l}\equiv\lambda(\omega_{l})$.
The order parameter is the frequency-dependent gap function
$\Delta_{n}=\Phi_{n}/Z_{n}$, where $Z_{n}=1+\Sigma_{n}/\omega_{n}$.

The main dimensionless parameter in the Migdal-Eliashberg theory is the
electron-phonon coupling defined as the electron-electron interaction
evaluated at zero Matsubara frequency, $\lambda=\lambda(\omega_{l}=0)$. The
coupling constant $\lambda$ is inversely proportional to the square of the
characteristic phonon frequency $\omega_{\mathrm{ch}}$. In conventional
electron-phonon models, such as the Holstein or Frölich Hamiltonian,
interactions with electrons strongly renormalize phonon frequencies [1, 2, 3].
Because of this it is crucial to distinguish bare ($\lambda_{0}$) and
renormalized ($\lambda$) coupling constants [4].

The Migdal-Eliashberg theory enjoyed a great deal of success over many decades
as a quantitative theory of conventional superconductivity. With phonon
spectrum extracted from experiment, the theory makes accurate predictions for
the superconducting transition temperature $T_{c}$, zero temperature energy
gap, jump in the specific heat at $T_{c}$, density of states, etc. for a broad
range of superconductors, such as Al, V, Ta, Sn, Tl, In, Nb, Pb, their various
alloys, Hg, and MgB2 [6, 7, 8, 9, 10, 5]. The coupling constant in these
materials ranges between $\lambda\approx 0.4$ for Al and $\lambda\approx 1.6$
for Hg in simple elements and up to $\lambda\approx 3.0$ (Pb0.5Bi0.5) in
alloys.

Despite the success and maturity of the Migdal-Eliashberg theory, there is
also a great deal of controversy surrounding it. A number of publications [11,
12, 13, 14, 15, 16] rediscover Migdal and Eliashberg’s observation [1, 2] that
the theory is inapplicable to $\lambda_{0}\gtrsim 0.5,$ Migdal theorem
notwithstanding [17]. Others see no upper limit on $\lambda$ [4] and even
hypothesize that the strong coupling, $\lambda\to\infty$, limit of the Migdal-
Eliashberg theory is attainable [18]. A recent study claims Eliashberg
equations acquire a one-parameter family of solutions in the strong coupling
limit and as a consequence the superconducting transition temperature vanishes
[19]. Even the $\lambda\to 0$ limit is not free from controversy. While it is
generally accepted that the weak coupling limit of the Migdal-Eliashberg
theory is the more famous Bardeen-Cooper-Schrieffer (BCS) theory [20], some
argue that this is in fact untrue [21, 22].

Our goal is to resolve these issues in a series of papers [24, 25, 23]. Here
we establish several important properties of the free energy of the electron-
phonon system. We show that the gap function $\Delta_{n}$ must be nonnegative
at the global minimum point and argue that the minimum is unique by symmetry
at all temperatures and coupling strengths, up to an overall phase
$e^{i\phi}$. Besides the global minimum, we find an infinite set of new saddle
points at strong coupling. As usual, stationary points other than the global
minimum do not affect equilibrium properties in the thermodynamic limit.
However, we reason that these saddle points play a major role in the far from
equilibrium dynamics of the electronic subsystem. Besides, their emergence and
proliferation seem to be tied to the subsequent breakdown of the Migdal-
Eliashberg theory discussed in the next paragraph. For the sake of
completeness, we also provide an educational appendix where we explain within
the path integral framework that the BCS theory undoubtedly is the weak
coupling limit of the Migdal-Eliashberg theory when this limit is properly
taken.

In the next paper [23], we will show that the Migdal-Eliashberg theory looses
its validity at a finite value of the renormalized electron-phonon coupling
$\lambda_{c}$ due to a phase transition which breaks the lattice translational
symmetry. We will then construct a new theory – theory of lattice-fermionic
superfluidity – that works past $\lambda_{c}$. Prior to the transition it
reduces to the Migdal-Eliashberg theory. Afterwards, it describes the new
state of the system.

Our main tool to achieve these goals is a representation of the Migdal-
Eliashberg theory as a classical Heisenberg spin chain, which we describe
below. Positions of the Eliashberg spins are the fermionic Matsubara
frequencies $\omega_{n}=\pi T(2n+1)$. Spin components are momentum-integrated
normal and anomalous Green’s functions. The interaction between them is
ferromagnetic and falls off as $(\omega_{n}-\omega_{m})^{-2}$ at large
separation. In addition, the spins are subject to a position-dependent Zeeman
field $2\pi\omega_{n}$ along the $z$-axis. The energy of the spin chain is
proportional to the free energy density of the electron-phonon system.
Classical spins provide a simple and intuitive picture of the normal and
superconducting states and of the transition between them. In the normal
state, the spins are parallel to the $z$-axis,
$\bm{S}_{n}=\mbox{sgn}(\omega_{n})\hat{\bm{z}}$, where $\hat{\bm{z}}$ is a
unit vector along the $z$-axis. Below the superconducting transition
temperature, they acquire $xy$-components softening the sharp domain wall
between $\omega_{-1}$ and $\omega_{0}$ in the normal state as shown in Fig. 1.

Not surprisingly, Eliashberg equations emerge as stationary point equations
for the free energy functional. Their solutions therefore correspond to either
minima or maxima or saddle points of the free energy. The spin chain
formulation makes it straightforward to demonstrate several fundamental
properties of these stationary points. In particular, we use it to prove that
the Eliashberg gap function at the global minimum is of the form
$\Delta_{n}=e^{i\phi}|\Delta_{n}|$, where $|\Delta_{n}|$ is even in
$\omega_{n}$. Furthermore, an infinite discrete set of saddle points emerges
as we increase $\lambda$. In terms of the spin chain, they are equilibria with
a certain number of spins flipped against effective magnetic fields acting on
them. The minimum is the stable, lowest energy spin configuration. The new
saddle points begin to proliferate just before the Migdal-Eliashberg theory
breaks down. We argue that these saddle points play a special role in the far
from equilibrium collisionless dynamics of strongly coupled conventional
superconductors. Namely, they are the fixed points of the corresponding
kinetic equations, which are Hamilton’s equations of motion for Eliashberg
spins. When sufficiently many spins are flipped, these stationary points are
unstable and give rise to rich soliton-like dynamics.

The content of this paper is as follows. In Sect. II we map the Migdal-
Eliashberg theory to a classical spin chain building on the path integral
formulation we develop in Appendices A and B. We interpret the superconducting
transition in terms of spins in Sect. III, construct a divergence free form of
Eliashberg equations as well as new (spin flip) solutions for them in Sect.
IV, and relate spins and the Eliashberg gap function $\Delta_{n}$ in Sect. V.
These three sections lay the groundwork for the analysis of stationary points
of the free energy in Sect. VI. Sect. VII focuses on the strong coupling limit
$\lambda=\infty$ and in Sect. VIII we compare the two notions of spins in the
theory of superconductivity – Eliashberg spins introduced in this paper and
Anderson pseudospins [26] – and discuss the role of the spin flip solutions of
the Eliashberg equations in the collisionless dynamics of strongly coupled
superconductors.

##  II Mapping to a Heisenberg spin chain

We work with two models in this paper. One is the Holstein model
(dispersionless phonons) with arbitrary hopping matrix and onsite potential.
The other is a rather general electron-phonon Hamiltonian with arbitrary
phonon dispersion and momentum-dependent electron-phonon coupling. The
Holstein Hamiltonian reads,

$H=\sum_{\bm{i}\bm{j}\sigma}h_{\bm{i}\bm{j}}c^{\dagger}_{\bm{i}\sigma}c_{\bm{j}\sigma}+\sum_{\bm{i}}\left[\frac{p_{\bm{i}}^{2}}{2M}+\frac{K_{0}x_{\bm{i}}^{2}}{2}\right]+\alpha\sum_{\bm{i}}n_{\bm{i}}x_{\bm{i}},$
(1)

where ${\bm{i}}$ and ${\bm{j}}$ label lattice sites, $h_{{\bm{i}}{\bm{j}}}$
are the matrix elements of an arbitrary single-electron Hamiltonian $\hat{h}$
[32], $c_{\bm{i}\sigma}$ annihilates an electron on site ${\bm{i}}$ with spin
projection $\sigma$,
$n_{\bm{i}}=\sum_{\sigma}c^{\dagger}_{{\bm{i}}\sigma}c_{{\bm{i}}\sigma}$ is
the fermion occupation operator, and $p_{\bm{i}}$ and $x_{\bm{i}}$ are ion
momentum and position operators. The bare phonon frequency is
$\Omega_{0}=\sqrt{K_{0}/M}$. The lattice and its dimensionality are at this
point arbitrary.

The Hamiltonian for electrons interacting with dispersing phonons is,

$\begin{split}H=&\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}+\sum_{\bm{q}}{\omega_{0}({\bm{q}})}b^{\dagger}_{\bm{q}}b_{\bm{q}}+\\\
&\frac{1}{\sqrt{N}}\sum_{{\bm{p}}{\bm{q}}\sigma}\frac{\alpha_{{\bm{q}}}}{\sqrt{2M\omega_{0}({\bm{q}})}}c^{\dagger}_{{\bm{p}}+{\bm{q}}\sigma}c_{{\bm{p}}\sigma}\left[b^{\dagger}_{-{\bm{q}}}+b_{\bm{q}}\right],\end{split}$
(2)

where $M$ is the ion mass and $N$ is the number of lattice sites.

We proceed with a path integral formulation of the Migdal-Eliashberg theory.
The first step is to integrate out phonons, which leaves us with an effective
fermion-fermion interaction quartic in fermionic fields. The second step is to
decouple the quartic term with three Hubbard-Stratonovich fields $\Phi$,
$\Sigma_{\uparrow}$ and $\Sigma_{\downarrow}$ that are functions of two
imaginary times and, in general, two space points. In the next step, we
integrate out the fermions to obtain an effective action solely in terms of
the Eliashberg fields $\Phi$, $\Sigma_{\uparrow}$ and $\Sigma_{\downarrow}$.
The fourth and final step is to obtain the stationary point of this effective
action, where the Eliashberg fields depend on the time difference only. We
detailed the above procedure in Appendix A for the Holstein model and Appendix
B for dispersing phonons (2). The end result is an expression for the free
energy functional [31] of the system per site,

$\begin{split}f=\nu_{0}T^{2}\sum_{nl}\left[\Phi_{n+l}^{*}\Lambda_{l}\Phi_{n}+\Sigma_{n+l}\Lambda_{l}\Sigma_{n}\right]-\\\
2\pi\nu_{0}T\sum_{n}\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}},\end{split}$
(3)

where $\nu_{0}$ is the density of states at the Fermi energy per site per spin
projection. The field $\Phi_{n}\equiv\Phi(\omega_{n})$ is complex and
$\Sigma_{n}\equiv\Sigma(\omega_{n})$ is real. Both fields are functions of the
fermionic Matsubara frequency $\omega_{n}$. At the stationary point, these
fields equal the anomalous and normal self-energies.

The quantity $\Lambda_{l}$ in Eq. (3) is the Fourier transform of
$1/\lambda(\tau)$ at bosonic Matsubara frequency $\omega_{l}=2\pi Tl$. Here
$\lambda(\tau)$ is the effective fermion-fermion interaction in the imaginary
time domain. In the Matsubara frequency domain this interaction reads for the
Holstein model,

$\lambda(\omega_{l})=\frac{g^{2}}{\omega_{l}^{2}+\Omega^{2}},\quad
g^{2}=\nu_{0}\alpha^{2}M^{-1}.$ (4)

We also define the dimensionless electron-phonon coupling constant as,

$\begin{split}\lambda=\lambda(\omega_{l}=0)=\frac{g^{2}}{\Omega^{2}}=\frac{\nu_{0}\alpha^{2}}{K},\\\
\end{split}$ (5)

where $\Omega$ and $K$ are the renormalized phonon frequency and spring
constant. See Appendix B for $\lambda(\omega_{l})$ and
$\lambda=\lambda(\omega_{l}=0)$ for dispersing phonons. At times we will
consider the strong coupling limit $\lambda\to\infty$, which is equivalent to
$\Omega\to 0$ or $K\to 0$ for the Holstein model.

The stationary point equations for the free energy (3), are the Eliashberg
equations [2],

$\displaystyle\Phi_{n}=\pi
T\sum_{m}\lambda_{nm}\frac{\Phi_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}},$
(6a) $\displaystyle\Sigma_{n}=\pi
T\sum_{m}\lambda_{nm}\frac{\omega_{m}+\Sigma_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}},$
(6b)

where

$\lambda_{nm}=\lambda(\omega_{n}-\omega_{m}).$ (7)

Note also that

$\lambda_{nn}=\lambda,$ (8)

diverges in the strong coupling limit.

It is convenient to introduce new variables – complex $F(\tau)$ and real
$G(\tau)$, such that,

$\Phi(\tau)=\pi\lambda(\tau)F(\tau),\quad\Sigma(\tau)=\pi\lambda(\tau)G(\tau).$
(9)

In frequency representation we have,

$\begin{split}\Phi_{n}=\pi T\sum_{m}\lambda_{nm}F_{m},\\\ \Sigma_{n}=\pi
T\sum_{m}\lambda_{nm}G_{m},\end{split}$ (10)

and Eq. (6) becomes,

$\begin{split}F_{n}=\frac{\Phi_{n}}{\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}},\\\
G_{n}=\frac{\omega_{n}+\Sigma_{n}}{\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}}.\end{split}$
(11)

On the stationary point, the fields $F_{n}$ and $G_{n}$ correspond to the
anomalous and normal Green’s functions integrated over the single-particle
energy, see Appendix A4. We also show in Appendix A4 that assuming time
reversal symmetry, $\Sigma_{n}$ is real and odd in $\omega_{n}$ and
$|\Phi_{n}|$ is even. Therefore, $G_{n}$ is real and odd and $|F_{n}|$ is
even.

Observe that Eq. (11) implies a constraint on the variables $G_{n}$ and
$F_{n}$,

$G_{n}^{2}+|F_{n}|^{2}=1.$ (12)

This means that we can trade these variables for three components of a
classical spin $\bm{S}_{n}$ of unit length, $\bm{S}^{2}_{n}=1$,

$S_{n}^{z}=G_{n},\quad S_{n}^{x}=\mathrm{Re}(F_{n}),\quad
S_{n}^{y}=\mathrm{Im}(F_{n}).$ (13)

It follows from Eq. (11) that,

$F_{n}\Phi_{n}^{*}+G_{n}(\omega_{n}+\Sigma_{n})=\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}.$
(14)

This allows us to rewrite the free energy (3) as,

$H_{s}\equiv\frac{f}{\nu_{0}T}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}T\sum_{nm}\lambda_{nm}{\bm{S}_{n}}\cdot{\bm{S}_{m}}.$
(15)

We interpret the free energy (15) as a Hamiltonian, $H_{s}$, of an open
classical Heisenberg spin chain in an inhomogeneous Zeeman magnetic field. The
positions of the spins are fermionic Matsubara frequencies $\omega_{n}$. Spin-
spin interactions are ferromagnetic and fall off at large “distance” as
$\lambda_{nm}\propto(\omega_{n}-\omega_{m})^{-2}$. The “magnetic field” is
linear in the position of the spin and goes to $\pm\infty$ as
$\omega_{n}\to\pm\infty$. Note that the Boltzmann weight is
$e^{-fN/T}=e^{-\nu_{0}NH_{s}}=e^{-H_{s}/\delta}$, where
$\delta=(\nu_{0}N)^{-1}$ is the single-particle (electron) level spacing in
the original electron-phonon problem. Therefore, the spin chain is at an
effective temperature,

$T_{s}=\delta.$ (16)

Let us also write down the classical spin Hamiltonian for the Holstein model
as a visual [substitute Eq. (4) into Eq. (15)],

$H_{s}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}Tg^{2}\sum_{nm}\frac{{\bm{S}_{n}}\cdot{\bm{S}_{m}}}{(\omega_{n}-\omega_{m})^{2}+\Omega^{2}}.$
(17)

See the list below Eq. (45) for more properties of the spin chain
representation of the free energy.

Sometimes an on-site Hubbard repulsion is added [6] to the Eliashberg
equations by replacing $\lambda_{nm}$ in Eq. (6a) with $\lambda_{nm}-u$. This
adds a long-range $xy$ term to the spin Hamiltonian (15),

$H_{C}=\pi^{2}Tu\sum_{nm}\left(S_{n}^{x}S_{m}^{x}+S_{n}^{y}S_{m}^{y}\right),$
(18)

and the classical spin Hamiltonian becomes $H_{s}+H_{C}$.

Figure 1: Superconducting transition in terms of classical spins. As shown in
the text, Migdal-Eliashberg theory maps to a classical Heisenberg spin chain.
The positions of the spins are fermionic Matsubara frequencies $\omega_{n}$.
The spin-spin interactions are purely ferromagnetic and the spins are subject
to a Zeeman magnetic field $2\pi\omega_{n}$ along the $z$-axis. The figure
shows (a) the normal state and (b) a superconducting state. In the
superconducting state, spins acquire $x$-components, which implies nonzero
anomalous Green’s function. The sharp domain wall in the normal state is
smeared in the superconducting state. The normal-superconductor transition
therefore translates into softening of the domain wall in the spin language.

##  III Superconductivity in terms of spins

In the spin language, the superconducting transition translates into softening
of the domain wall as shown in Fig. 1. This is similar to the Anderson
pseudospin description of the BCS superconductivity [26] even though
Eliashberg and Anderson spins are not the same (see Sect. VIII).

To understand what happens as we lower the temperature, let us analyze the
spin texture minimizing the energy $H_{s}$ as a function of $T$. The Zeeman
magnetic field and ferromagnetic interaction compete in $H_{s}$. The spin
configuration minimizing the Zeeman term is
$\bm{S}_{n}=\mbox{sgn}(\omega_{n})\hat{\bm{z}}$, and the Zeeman field
inevitably prevails far from the origin, so that
$\bm{S}_{n}\to\pm\hat{\bm{z}}$ for $\omega_{n}\to\pm\infty$. These fixed
values at infinity serve as boundary conditions for the interaction term.

Above the superconducting $T_{c}$, the anomalous averages vanish, $F_{n}=0$.
According to the definition (13) of the classical spin, this means that all
spins are parallel to the $z$-axis. From the behavior of $\bm{S}_{n}$ at large
$\omega_{n}$ and by symmetry, it is then clear that the minimum energy spin
texture is $\bm{S}_{n}=\mbox{sgn}(\omega_{n})\hat{\bm{z}}$. This is the normal
state in the spin language. The characteristic feature of the normal state is
a sharp domain wall between $\omega_{-1}$ and $\omega_{0}$ with an abrupt
maximal jump of the $z$-component of spin from $S_{-1}^{z}=-1$ to
$S_{0}^{z}=+1$, see Fig. 1.

Below $T_{c}$ the anomalous averages are nonzero, i.e., the spins acquire
$x$-components ($F_{n}$ can be made real in the spin chain ground state). This
means softening of the domain wall. Now the change in $S_{n}^{z}$ from $-1$ at
$\omega_{n}\to-\infty$ to $+1$ at $\omega_{n}\to+\infty$ is gradual and the
jump $S_{0}^{z}-S_{-1}^{z}<2$. Spin configurations with nonzero
$xy$-components are superconducting states.

##  IV  New solutions of Eliashberg equations

We begin this section by verifying that solutions of the Eliashberg equations
are equilibrium points of the spin chain. This correspondence allows us to
identify new solutions of the Eliashberg equations both with zero and nonzero
anomalous self-energy $\Phi_{n}$, i.e., both normal and superconducting. These
solutions correspond to equilibria of the spin chain with a number of spins
flipped against their effective magnetic fields, as opposed to all spins being
along their fields as in the ground state.

###  1 Eliashberg equations as equilibrium condition for spins

Consider our spin Hamiltonian in its most general form,

$H_{s}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}T\sum_{nm}\lambda_{nm}{\bm{S}_{n}}\cdot{\bm{S}_{m}},$
(19)

Terms in the Hamiltonian that contain spin $\bm{S}_{n}$ are,

$h_{n}=-\bm{b}_{n}\cdot\bm{S}_{n},\quad\bm{b}_{n}=2\pi\omega_{n}\hat{\bm{z}}+2\pi^{2}T\sum_{m\neq
n}\lambda_{nm}\bm{S}_{m},$ (20)

where $\bm{b}_{n}$ is the effective magnetic spin acting on spin $\bm{S}_{n}$,
which is the Zeeman field plus the field from other spins.

Equilibrium is when each spin is collinear with its field (parallel or
antiparallel), i.e.,

$\bm{S}_{n}=e_{n}\frac{\bm{b}_{n}}{|\bm{b}_{n}|},\quad e_{n}=\pm 1,$ (21)

where $e_{n}=+1$ indicates that the spin $\bm{S}_{n}$ is parallel to its
effective field $\bm{b}_{n}$, while $e_{n}=-1$ – that the spin is flipped
(antiparallel to the field). However, not all spin equilibria correspond to
stationary points of the Eliashberg free energy (3). The reason is that we
extracted the square root in Eq. (14) to obtain the spin chain. Only such spin
equilibria are stationary points of the free energy for which the left hand
side of Eq. (14) is nonnegative.

To interpret this condition in terms of spins, let us introduce an additional
magnetic field, which includes the action of the spin on itself,

$\bm{B}_{n}=\bm{b}_{n}+2\pi^{2}T\lambda\bm{S}_{n}.$ (22)

Recall that $\lambda=\lambda(\omega_{l}=0)=\lambda_{nn}$. The definition of
$\bm{b}_{n}$ in Eq. (20), the definition of spins (13), and Eq. (10) imply,

$B_{n}^{+}=2\pi\Phi_{n},\quad B_{n}^{z}=2\pi(\omega_{n}+\Sigma_{n}).$ (23)

Here and below we use the notation $V^{+}\equiv V^{x}+iV^{y}$, where $V^{x}$
and $V^{y}$ are the $x$ and $y$ components of a vector $\bm{V}$. Since in
equilibrium $\bm{b}_{n}$ is collinear with $\bm{S}_{n}$, so is $\bm{B}_{n}$.
Then, the requirement that the expression
$F_{n}\Phi_{n}^{*}+G_{n}(\omega_{n}+\Sigma_{n})$ on the right hand side of Eq.
(3) be nonnegative is equivalent to,

$\bm{B}_{n}\cdot\bm{S}_{n}=\bm{b}_{n}\cdot\bm{S}_{n}+2\pi^{2}\lambda T\geq 0.$
(24)

This condition always holds when $e_{n}=+1$ for all $n$ (no spin flips).

We saw that in equilibrium $\bm{S}_{n}$ is either parallel or antiparallel to
$\bm{B}_{n}$. Equilibria of the spin chain with each $\bm{S}_{n}$ parallel to
$\bm{B}_{n}$ correspond to solutions of the Eliashberg equations and vice
versa. We consider only this type of equilibria in this paper. However, it is
important to keep in mind that spin flips are nevertheless allowed as
$\bm{S}_{n}$ can be parallel to $\bm{B}_{n}$, but antiparallel to
$\bm{b}_{n}$. Later in this section we will point out equilibrium spin
configurations of this type. To derive Eliashberg equations from
$\bm{S}_{n}\parallel\bm{B}_{n}$, note that this condition is equivalent to,

$S_{n}^{z}=\frac{B_{n}^{z}}{|\bm{B}_{n}|},\quad
S_{n}^{+}=\frac{B_{n}^{+}}{|\bm{B}_{n}|}.$ (25)

Eq. (22) now implies the following two self-consistency conditions:

$B_{n}^{z}=\sum_{m}\lambda_{nm}\frac{B_{m}^{z}}{|\bm{B}_{m}|},\quad
B_{n}^{+}=\sum_{m}\lambda_{nm}\frac{B_{m}^{+}}{|\bm{B}_{m}|}.$ (26)

Using Eq. (23), we see that these are precisely the Eliashberg equations (6).

The spin chain approach naturally leads to an alternative set of Eliashberg
equations that are free of divergencies in the strong coupling limit
$\lambda\to\infty$ both in the superconducting and normal states [18]. Indeed,
we see from Eqs. (6) and (8) that this limit presents a problem for the
Eliashberg equations as they contain a diverging $\lambda_{nn}=\lambda$ term.
As a result, both self-energies, $\Sigma_{n}$ and $\Phi_{n}$, are proportional
to $\lambda$ and diverge in this limit. Eq. (10) shows that the divergent
parts of $\Sigma_{n}$ and $\Phi_{n}$ are the $m=n$ terms, $\pi\lambda TG_{n}$
and $\pi\lambda TF_{n}$, respectively. Since $G_{n}$ and $|F_{n}|$ are both of
order 1 [see Eq. (12)], these terms diverge in the strong coupling limit.

All we have to do to avoid these divergences is to use Eq. (21) instead of Eq.
(25). Let us introduce reduced self-energies [cf. Eq. (10)],

$\begin{split}\omega_{n}+\Sigma^{\prime}_{n}&=(\omega_{n}+\Sigma_{n})-\pi
T\lambda G_{n}=\frac{b_{n}^{z}}{2\pi},\\\ \Phi^{\prime}_{n}&=\Phi_{n}-\pi
T\lambda F_{n}=\frac{b_{n}^{+}}{2\pi}.\end{split}$ (27)

$\Sigma^{\prime}_{n}$ and $\Phi^{\prime}_{n}$ are the same as $\Sigma_{n}$ and
$\Phi_{n}$ in Eq. (10) but without the $n=m$ terms. It follows from Eqs. (21)
and the definition (20) of $\bm{b}_{n}$ that,

$\begin{split}\Phi^{\prime}_{n}=\pi T\sum_{m\neq
n}\lambda_{nm}\frac{\Phi^{\prime}_{m}}{\sqrt{(\omega_{m}+\Sigma^{\prime}_{m})^{2}+|\Phi^{\prime}_{m}|^{2}}},\\\
\Sigma^{\prime}_{n}=\pi T\sum_{m\neq
n}\lambda_{nm}\frac{\omega_{m}+\Sigma^{\prime}_{m}}{\sqrt{(\omega_{m}+\Sigma^{\prime}_{m})^{2}+|\Phi^{\prime}_{m}|^{2}}},\end{split}$
(28)

where for simplicity we took $e_{n}=+1$ (no spin flips). As we discuss below,
this is always the case at not too strong coupling and only these $e_{n}=+1$
solutions have been considered in the literature until now. Then, Eqs. (28)
are equivalent to the original Eliashberg equations (6), but do not contain
$n=m$ terms. Note also that Eqs. (21) and (25) imply
$S_{n}^{z}/S_{n}^{+}=b_{n}^{z}/b_{n}^{+}=B_{n}^{z}/B_{n}^{+}$ and therefore,

$\frac{\Phi_{n}}{\omega_{n}+\Sigma_{n}}=\frac{\Phi^{\prime}_{n}}{\omega_{n}+\Sigma^{\prime}_{n}}.$
(29)

The $z$-component of spin, $S^{z}(\omega_{n})\equiv S_{n}^{z}=G_{n}$, is odd
in $\omega_{n}$, because $G_{n}$ is odd. Now consider an equilibrium
configuration such that $\mbox{sgn}(S_{n}^{z})=\mbox{sgn}(\omega_{n})$. Taking
into account that $S_{n}^{z}$ is odd, the expression (20) for $\bm{b}_{n}$ and
Eq. (22) imply
$\mbox{sgn}(b_{n}^{z})=\mbox{sgn}(B_{n}^{z})=\mbox{sgn}(\omega_{n})$. Since
$b_{n}^{z}$ and $S_{n}^{z}$ have the same sign, this must be a spin
configuration with no spin flips. Though less obvious, the converse is also
true, i.e.,
$\mbox{sgn}(b_{n}^{z})=\mbox{sgn}(B_{n}^{z})=\mbox{sgn}(\omega_{n})$ holds in
any configuration with no spin flips.

It is also not difficult to show that
$\mbox{sgn}(b_{n}^{z})=\mbox{sgn}(\omega_{n})$ continues to hold when only a
small number of spins are flipped out of this configuration (which we will
assume). On the other hand, we saw above that $\bm{S}_{n}$ must be parallel to
$\bm{B}_{n}$ in any spin equilibrium that corresponds to a solution of
Eliashberg equations. Therefore, for such equilibria with few spin flips we
have,

$\mbox{sgn}(B_{n}^{z})=\mbox{sgn}(S_{n}^{z})=e_{n}\mbox{sgn}(\omega_{n}).$
(30)

###  2 Spin flip solutions

Let us determine when spin flip solutions of Eliashberg equations exist. The
spin chain has equilibria with any number of spins at arbitrary positions
being antiparallel to their fields. However, to be also solutions of the
Eliashberg equations these configurations must satisfy the inequality (24).
With the help of Eq. (21) we rewrite this inequality as,

$e_{n}|\bm{b}_{n}|+2\pi^{2}\lambda T\geq 0.$ (31)

First, it is clear that $e_{n}=+1$ always works, i.e., stationary
configurations of spins (equilibria) with no spin flips are always solutions
of the Eliashberg equations.

Now suppose $\bm{S}_{n}$ is antiparallel to its effective field $\bm{b}_{n}$.
Then, $e_{n}=-1$ and Eq. (31) becomes,

$2\pi^{2}\lambda T\geq|\bm{b}_{n}|.$ (32)

We see that no spin flips are allowed at weak coupling, $\lambda\to 0$, and at
$T=0$ for any finite $\lambda$, because $2\pi^{2}\lambda T$ in this inequality
is either negligible or zero, while the right hand side is positive.
Conversely, arbitrary spin flips are permitted in the strong coupling limit
$\lambda\to\infty$ at any finite $T$, because this term is infinite and
positive. Therefore, we expect that at finite temperature spin flip solutions
first appear at $\lambda\sim 1$.

The existence of spin flip solutions in the double limit $\lambda\to\infty$
and $T\to 0$ depends on whether the product $\lambda T$ goes to zero (spin
flip solutions do not exist), infinity (solutions with any set of $\bm{S}_{n}$
antiparallel to $\bm{b}_{n}$ exist), or a finite value (spin flip solutions
may or may not exist). Recall that $\lambda=g^{2}/\Omega^{2}$ and that we are
increasing $\lambda$ by decreasing $\Omega$. Then, arbitrary spin flips are
allowed as long as $Tg/\Omega^{2}\to\infty$ in the double limit $T\to 0$ and
$\Omega\to 0$ and there are numerous such solutions in the regime
$T\gg\Omega^{2}/g=g/\lambda=\Omega/\sqrt{\lambda}$.

We mentioned in the previous subsection that $m=n$ contributions to the self-
energies $\Sigma_{n}$ and $\Phi_{n}$ in Eqs. (10) diverge in the strong
coupling limit. It is precisely these divergent terms that make spin flip
solutions possible. Indeed, Eqs. (22–24) and the definition of $\bm{b}_{n}$ in
Eq. (20) show that the second term in the inequality (31), without which the
spin flip solutions would not exist, arises from just these two terms.

Spin flip solutions occur already in the normal state where
$\bm{S}_{n}=\mbox{sgn}(\omega_{n})\hat{\bm{z}}$. Suppose we flip the spin at
Matsubara frequency $\omega_{k}>0$. Since $S_{n}^{z}$ is odd in $\omega_{n}$,
we also have to symmetrically flip the spin at $-\omega_{k}$. The new spin
configuration is: $\bm{S}_{n}=\mbox{sgn}(\omega_{n})\hat{\bm{z}}$ for
$|\omega_{n}|\neq\omega_{k}$ and
$\bm{S}_{n}=-\mbox{sgn}(\omega_{n})\hat{\bm{z}}$ for
$|\omega_{n}|=\omega_{k}$. Let us use the Holstein model (17) for simplicity.
The inequality (31) becomes,

$\begin{split}\lambda\geq
2k+1+\sum_{m=1}^{k+1}\frac{2g^{2}}{4\pi^{2}T^{2}m^{2}+\Omega^{2}}.\end{split}$
(33)

We see that the first spin flip solution of the Eliashberg equations appears
at $\lambda\geq 1$ for $k=0$ and $T\to\infty$, the next one at $\lambda\geq 3$
etc. A similar analysis in the normal and superconducting states at lower $T$
again shows that these solutions appear at $\lambda\gtrsim 1$ and proliferate
at larger $\lambda$.

Let us also determine the lowest temperature at which normal
($\Delta_{n}\equiv 0$) solutions with spins antiparallel to their fields
exist. It costs least to flip the spins at the first two Matsubara frequencies
$\omega_{-1}$ and $\omega_{0}$. Setting $k=0$ in Eq. (33) and keeping in mind
that $\lambda=g^{2}/\Omega^{2}$, we find,

$T\geq\frac{\Omega}{2\pi}\sqrt{\frac{\lambda+1}{\lambda-1}}.$ (34)

At strong coupling, this becomes $T\geq\Omega/(2\pi)$, which implies the
condition $T\gg\Omega/\sqrt{\lambda}$ discussed above. Note also that in the
strong coupling limit $\Omega\to 0$, the quantity $\lambda T$ diverges as
$g^{2}/\Omega$ or stronger when $T$ satisfies the inequality (34).

We argue in Sect. VIII2 that spin flip solutions are important for
understanding the short-time nonlinear dynamics of strongly coupled
superconductors – most of them are unstable fixed points of kinetic equations,
similar to an inverted pendulum, that generate soliton-like waves with a rich
structure. As mentioned in the Introduction, the Migdal-Eliashberg theory
breaks down at a certain finite $\lambda_{c}$. Interestingly, both the
emergence of the spin flip saddle points at $\lambda<\lambda_{c}$ and the
breakdown of the theory originate from the same source. We saw in this section
that it is the divergent $n=m$ terms in the self energies that bring about the
spin flip solutions. In a subsequent paper [23] we will see that these terms
are also among the main culprits responsible for the breakdown of the Migdal-
Eliashberg theory. In this connection, let us mention the mountain pass
theorem [33] according to which there must be a saddle point between two
minima. Therefore, it may be that the emergence of these saddle points
indicates that the free energy develops an additional minimum besides the one
described by the Eliashberg equations.

##  V Spins and superconducting gap function

In this section, we introduce the superconducting gap function $\Delta_{n}$,
write down a new, more general gap equation that accounts for spin flips, and
relate spin components to the gap function.

As usual, we introduce the Eliashberg gap function
$\Delta(\omega_{n})\equiv\Delta_{n}$ through equations,

$\omega_{n}+\Sigma_{n}=\omega_{n}Z_{n},\quad\Phi_{n}=\Delta_{n}Z_{n}.$ (35)

Eqs. (23) and (30) imply,

$\mbox{sgn}(Z_{n})=e_{n}.$ (36)

In other words, $Z_{n}$ is negative if spin $\bm{S}_{n}$ is flipped and
positive otherwise. As far as we are aware, only solutions with positive
$Z_{n}$, i.e., with no spin flips have been considered in the literature until
now.

Since $\Sigma_{n}$ is odd and $|\Phi_{n}|$ is even in frequency, $Z_{n}$ and
$|\Delta_{n}|$ are both even. With the substitution (35), the Eliashberg
equation for $\Sigma_{n}$ in (6) determines $Z_{n}$ for a given $\Delta_{n}$,

$Z_{n}=1+\frac{\pi
T}{\omega_{n}}\sum_{m}\lambda_{nm}\frac{e_{m}\omega_{m}}{\sqrt{\omega_{m}^{2}+|\Delta_{m}|^{2}}},$
(37)

while the equation for $\Phi_{n}$ becomes the gap equation,

$\omega_{n}\Delta_{n}=\pi
T\sum_{m}\lambda_{nm}\frac{\omega_{n}e_{m}\Delta_{m}-\Delta_{n}e_{m}\omega_{m}}{\sqrt{\omega_{m}^{2}+|\Delta_{m}|^{2}}},$
(38)

where we used Eq. (37) to eliminate $Z_{n}$. As soon as the sign of $Z_{n}$ is
fixed, the gap equation decouples from $Z_{n}$ and the Eliashberg equations
reduce to a single equation – the gap equation (38). The gap equation and the
gap function determine various thermodynamic properties of the system, for
example, the specific heat, superconducting $T_{c}$, and the condensation
energy. Continued to real frequencies, the gap function determines the density
of states and various response functions, such as the optical conductivity.

The $m=n$ term vanishes in the gap equation (38). This equation is therefore
divergence free in the limit $\lambda_{nn}=\lambda\to\infty$. The same is not
true for $Z_{n}$ which diverges in this limit, because the $m=n$ term does not
cancel from Eq. (37). To get rid of this divergence, we use the reduced self-
energies defined in Sect. IV and introduce $Z_{n}^{\prime}$ and $\Delta_{n}$
through,

$\omega_{n}+\Sigma^{\prime}_{n}=\omega_{n}Z^{\prime}_{n},\quad\Phi^{\prime}_{n}=\Delta_{n}Z^{\prime}_{n}.$
(39)

Eq. (29) guarantees that $\Delta_{n}$ here is the same as in Eq. (35).
Therefore, the gap equation remains the same, while the expression for
$Z^{\prime}_{n}$ is,

$Z^{\prime}_{n}=1+\frac{\pi T}{\omega_{n}}\sum_{m\neq
n}\lambda_{nm}\frac{e_{m}\omega_{m}}{\sqrt{\omega_{m}^{2}+|\Delta_{m}|^{2}}}.$
(40)

The only difference between $Z^{\prime}_{n}$ and $Z_{n}$ is that the $m=n$
term is absent in $Z^{\prime}_{n}$.

Spin components in terms of the gap function are,

$S_{n}^{z}=\frac{e_{n}\omega_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}},\quad
S_{n}^{+}=\frac{e_{n}\Delta_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}},$
(41)

as follows from Eqs. (11) and (13). Since the interaction is purely
ferromagnetic, it is clear that for fixed $S_{n}^{z}$ the free energy is
minimal when the $xy$-projections of all spins, $\bm{S}^{\perp}_{n}$, point in
the same direction. Without loss of generality we take this direction to be
the $x$-axis. The spin Hamiltonian (15) becomes,

$H_{s}=-2\pi\sum_{n}\omega_{n}\cos\theta_{n}-\pi^{2}T\sum_{nm}\lambda_{nm}\cos(\theta_{n}-\theta_{m}),$
(42)

where $\theta_{n}$ is the angle the spin $\bm{S}_{n}$ makes with the $z$-axis,
$S_{n}^{z}=\cos\theta_{n}$ and $S_{n}^{x}=\sin\theta_{n}$. The stationary
point of $H_{s}$ with respect to $\theta_{n}$ is,

$\omega_{n}\sin\theta_{n}=\pi
T\sum_{m}\lambda_{nm}\sin(\theta_{m}-\theta_{n}).$ (43)

With the help of Eq. (41), we see that this is nothing but the gap equation
(38) written in terms of the angular variable $\theta_{n}$ and that
$\Delta_{n}=\omega_{n}\tan\theta_{n}$.

Let us also mention here various expressions for the free energy of boson
mediated superconductors that have been proposed over the years. Eliashberg
derived one such expression in 1962 [34], see also Ref. 35. Notably, Bardeen
and Stephen used Eliashberg’s result to evaluate the corrections to the BCS
condensation energy due to finite ratio of the BCS gap to the Debye energy
[36]. We discovered the spin chain representation of the free energy (17) in
2002 announcing it afterwards at various venues [38, 37]. To reconcile our
result with Eliashberg’s, Haslinger and Chubukov obtained a somewhat less
general (with all $e_{n}=+1$) answer for the free energy of the Holstein model
in terms of $\Delta_{n}$ by integrating Eliashberg’s expression over the
electron momentum [39] . This answer follows from Eq. (17) if we substitute
Eq. (41) with $e_{n}=+1$ into it. Subsequently, this answer has been used in a
number of publications.

It is important to emphasize that our expression (3) for the free energy is
more general than Eq. (17) and therefore Eliashberg’s expression. Eq. (3)
takes into account certain fluctuations of the fields $\Sigma$ and $\Phi$ in
imaginary time which Eq. (17) does not. The two expressions are equivalent
only at $\lambda=\infty$, while away from this limit they coincide only at the
stationary points. We will discuss this in more detail below.

##  VI Stationary points of the free energy

This section contains several results regarding stationary points of the free
energy. One is that at the minimum of the free energy
$\Delta(\omega_{n})=e^{i\phi}|\Delta_{n}|$, where $|\Delta_{n}|$ is even in
$\omega_{n}$ and $e^{i\phi}$ is an overall phase. In other words, up to an
overall phase $\Delta(\omega_{n})$ is nonnegative and even. Another is that
any spin configuration where $\Delta(\omega_{n})$ is real, but changes sign
cannot be a local or global minimum and is higher in energy than the global
minimum by an amount proportional to the system size. In particular, spin
flips correspond to saddle points. Fixing the overall phase, we also show that
there cannot be a continuous family of stationary points connected to the
global minimum and argue based on symmetry that the global minimum is unique.

We have seen in Sect. IV that stationary points of the free energy are the
equilibria of the spin chain. At weak coupling only equilibria where each spin
is parallel to its magnetic field ($e_{n}=+1$) are stationary points. As
$\lambda$ increases, equilibrium configurations with arbitrary number of spins
antiparallel to their fields ($e_{n}=-1$) become stationary points as well.
These spin flip solutions of the Eliashberg equations are necessarily saddle
points. Indeed, suppose $\bm{S}_{n}$ is antiparallel to $\bm{b}_{n}$. Then,
its contribution $h_{n}=-\bm{b}_{n}\cdot\bm{S}_{n}$ to the spin Hamiltonian is
positive and decreases when we rotate $\bm{S}_{n}$ keeping all other spins
fixed. Given that there are also spins parallel to their fields, the
contribution of any of them increases when it deviates from its equilibrium
position. Therefore, such solutions cannot be local minima or maxima, but are
saddle points.

We see that at a minimum all spins must be parallel to their fields. This
condition is only necessary, but not sufficient as there can be a collective
mode involving many spins that lowers the free energy. This happens, for
example, when the normal state looses stability at $T=T_{c}$. In the normal
state all spins are parallel to their fields and yet it is a saddle point and
not a minimum below $T_{c}$. Hence for the minimum we take $e_{n}=+1$ for all
$n$ in Eq. (41),

$S_{n}^{z}=\frac{\omega_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}},\quad
S_{n}^{+}=\frac{\Delta_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}},$ (44)

and

$\omega_{n}\Delta_{n}=\pi
T\sum_{m}\lambda_{nm}\frac{\omega_{n}\Delta_{m}-\Delta_{n}\omega_{m}}{\sqrt{\omega_{m}^{2}+|\Delta_{m}|^{2}}}.$
(45)

is the corresponding gap equation.

Before we proceed to analyze the global minimum, it is helpful to summarize
the main features of the spin chain representation of the free energy $f$:

  1. (a)

The spin chain formula, $f=\nu_{0}TH_{s}$, for the free energy holds at
stationary points of $f$. It does not necessarily hold away from such points
because we used the Eliashberg equations (6) in deriving it. However, the
strong coupling limit, $\lambda\to\infty$, is an exception. As we show in
Sect. VII, in this limit $f=\nu_{0}TH_{s}$ holds at every point
$(G_{n},F_{n})$ of the configuration space of the system, not only at the
stationary points.

  2. (b)

Every stationary point of $f$ is also a stationary point of $H_{s}$ and every
stationary point of $H_{s}$ with $e_{n}=+1$ is a stationary point of $f$. Most
importantly, the global minimum of $f$ is the global minimum of $H_{s}$. We
also saw in Sect. IV that in $\lambda\to\infty$ limit $e_{n}$ are
unconstrained and therefore the correspondence between stationary points of
$f$ and $H_{s}$ is one to one.

For easier reference, let us write down the spin chain Hamiltonian (19) here
in a slightly different form,

$H_{s}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}T\sum_{nm}\lambda_{nm}({\bm{S}_{n}}\cdot{\bm{S}_{m}}-1),$
(46)

and remind that $\lambda_{nm}>0$ for all $n$ and $m$. In Eq. (46) we
regularized the otherwise divergent interaction term by subtracting a constant
from it [40]. This does not affect the stationary point (gap) equation in any
way.

###  1  Global minimum of the free energy

Here we establish the following properties of the free energy of the electron-
phonon system:

  1. (1)

At the global minimum $\Delta(\omega_{n})$ is a nonnegative and even function
of $\omega_{n}$ up to an overall phase. Moreover, either
$\Delta(\omega_{n})=0$ for all $n$ (normal state) or $\Delta(\omega_{n})>0$
for all $n$ (superconducting).

  2. (2)

For a fixed overall phase, there are no curves of stationary points that
contain the global minimum.

  3. (3)

The only symmetry of the spin chain is the symmetry with respect to rotations
around the $z$-axis.

  4. (4)

Any $\Delta(\omega_{n})$ that is real, but changes sign with $\omega_{n}$
cannot be a minimum (local or global) of $H_{s}$.

These statements remain true at arbitrary temperature, including $T=0$ [42].
Let us begin with the proof of property (1). Suppose $S_{n}^{z}$ are fixed.
Then, we have to minimize the $xy$-part of the spin Hamiltonian. Clearly,
since $\lambda_{nm}>0$ for all $n$ and $m$, this requires all
$\bm{S}^{\perp}_{n}$ being aligned in the same direction. Any deviation from
this alignment means a finite energy cost in $H_{s}$ and because the total
free energy is $Nf=\nu_{0}NTH_{s}$, the total energy cost is proportional to
the total number of sites $N$. As before we designate the direction in which
$\bm{S}^{\perp}_{n}$ align to be the positive $x$-direction [43]. Eq. (44)
then implies that $\Delta_{n}$ is real and, most importantly, it must be
nonnegative for all $n$ at the global minimum. Since $|\Delta_{n}|$ is even,
$\Delta_{n}=|\Delta_{n}|$ must also be an even function of frequency.
Therefore, $\Delta(\omega_{n})$ is a nonnegative and even function. Given that
$\Delta_{n}$ is nonnegative, the gap equation (45) further shows that either
all $\Delta_{n}$ are zero or none.

Clearly, if $\Delta_{n}$ is the global minimum, so is $e^{i\phi}\Delta_{n}$.
Note that $e^{i\phi}\Delta_{n}$ defines a closed curve (circle) with no end
points in the configuration space. Let us show that there can be no other
1-parameter families of solutions of the Eliashberg equations that contain the
global minimum [property (2)]. Suppose such a family exists. It then similarly
defines a curve in the configuration space and both $f$ and $H_{s}$ must be
constant along this curve. It follows that all their derivatives vanish on
this curve as well. We saw above that at the global minimum we can always fix
the global phase so that $S_{n}^{y}=0$ and $S_{n}^{x}\geq 0$. Eq. (44) shows
that $S_{n}^{z}\geq 0$ for $\omega_{n}\geq 0$. Therefore, the angle
$\theta_{n}$ that spin $\bm{S}_{n}$ makes with the $z$-axis is in the range
$[0,\pi/2]$. Then, the 1-parameter family in question must have end points. As
all derivatives of $H_{s}$ vanish along the curve of solutions before the end
point and at least some of the first derivatives must be nonzero after the end
point, $f$ cannot be an analytic function of all $\theta_{n}$. This
contradicts Eq. (42), which shows that $H_{s}$ is analytic in all $\theta_{n}$
[44].

The model (46) is an inhomogeneous Heisenberg spin chain. Conditions under
which this type of models are integrable for both quantum and classical spins
have been studied in the literature [45, 46, 47, 48, 49, 50]. In particular,
there is a class of integrable models that describe BCS-like pairing between
fermions, which become spin chains when written in terms of Anderson spins
[48, 49, 50]. These models, integrable for both quantum and classical spins,
are similar to the Hamiltonian (46) in that there is an inhomogeneous Zeeman
field linear in the spin coordinate and interactions between spins are
ferromagnetic. Models of this type are integrable only for special choices of
coupling constants $\lambda_{nm}$. Integrable cases are a subset of measure
zero among all Hamiltonians of this form, because we need only of order
$N_{s}$ parameters to specify $\lambda_{nm}$ in integrable cases as opposed to
$N_{s}^{2}$ in the generic case, where $N_{s}$ is the number of spins. It is
straightforward to verify that our $\lambda_{nm}$ do not correspond to any
known integrable model. It is safe to conjecture that the Hamiltonian (46) is
not integrable in the sense of Ref. 51, i.e., that it does not possess
parameter-dependent integrals of motion (parameters here are $\omega_{n}$ and
$\lambda_{nm}$). Then, any symmetry must Poisson-commute with each term in the
Hamiltonian individually, i.e., with each $S_{n}^{z}$ and
$\bm{S}_{n}\cdot\bm{S}_{m}$. The only such symmetry is the $z$-component of
the total spin $J_{z}=\sum_{n}S_{n}^{z}$. Thus, we arrive at property (3)
assuming non-integrability of the spin chain (46).

Note that conservation of $J_{z}$ is responsible for the degeneracy of the
global minimum with respect to rotations around the $z$-axis, i.e., with
respect to the overall phase of $\Delta_{n}$. Since there are no other
symmetries, we expect the global minimum to be generally unique (barring
accidental degeneracy) apart from the arbitrary overall phase of $\Delta_{n}$.

Now consider a spin configuration with real sign-alternating
$\Delta(\omega_{n})$. We already know from property (1) that such a
configuration cannot be a global minimum and is associated with a macroscopic
energy cost. Let us investigate if it can at least be a local minimum. Eq.
(44) implies that the projection of spin $\bm{S}_{n}$ onto the $xy$-plane is
$\bm{S}_{n}^{\perp}=\mbox{sgn}(\Delta_{n})|\bm{S}_{n}^{\perp}|\hat{\bm{x}}$,
where $\hat{\bm{x}}$ is the unit vector along the $x$-axis. Let us uniformly
rotate all spins with $\mbox{sgn}(\Delta_{n})<0$ around the $z$-axis towards
the positive $x$-axis by a small angle $\delta\phi$, i.e., for these spins
$S_{n}^{+}\to S_{n}^{+}e^{i\delta\phi}$. This does not change the Zeeman
energy and the interaction energies of spins with $\mbox{sgn}(\Delta_{n})<0$
among themselves and of spins with $\mbox{sgn}(\Delta_{n})\geq 0$ among
themselves, while the ferromagnetic interaction energy of
$\mbox{sgn}(\Delta_{n})<0$ spins with $\mbox{sgn}(\Delta_{n})>0$ decreases as
they become more aligned with each other. Since the energy decreases for an
infinitesimal deviation from this configuration, it cannot be a local minimum
of $H_{s}$. At $\lambda=\infty$, when $f=\nu_{0}TH_{s}$ holds in the entire
configuration space, spin configurations with sign-alternating
$\Delta(\omega_{n})$ cannot be local minima of $f$ as well. This proves
property (4).

In connection with the above, let us mention papers by Wu et. al. that claim
that at $T=0$ and $\lambda=\infty$ there is a 1-parameter family of sign-
alternating solutions of the Eliashberg equations [19, 52, 53]. One end point
of this purported family of solutions is said to be the global minimum and the
other – the normal state, i.e., the free energy must increase along this
continuous family of stationary points (see, e.g., Fig. 4 of Ref. 19 and Fig.
1 of Ref. 52). This is impossible because any curve of stationary points must
also be a curve of constant free energy, such as, e.g., the curve traced out
by changing the global phase $\phi$. It also contradicts property (2) above.

Moreover, we saw above that the energy cost of configurations with sign-
alternating $\Delta_{n}$ is proportional to the system volume. Their Boltzmann
weight relative to the global minimum is $e^{-N\Delta f/T}$, where $N$ is the
number of sites and $\Delta f$ is the free energy density difference (which is
finite in $N\to\infty$ limit). Naturally, such states cannot contribute to
thermal equilibrium properties of a bulk system. In the present case not only
$N\to\infty$, but also $T\to 0$, so the global minimum becomes the ground
state and entirely determines all equilibrium physics.

##  VII Strong coupling limit

There are two universal limits of the Migdal-Eliashberg theory – the weak and
strong coupling limits, $\lambda\to 0$ and $\lambda\to\infty$. ‘Universal’
here means that the theory becomes independent of the phonon dispersion and
momentum dependence of the electron-phonon coupling. In both limits there is
only one low energy scale. The weak coupling limit of the Migdal-Eliashberg
theory is the BCS theory, see Appendix B2.

We review the strong coupling limit [6, 8, 27, 28, 29, 30] in this section for
two related reasons. First, the derivation of the spin chain representation of
the free energy does not rely on Eliashberg equations in this limit and
therefore holds everywhere in the configuration space. Second, we will need it
in a later study. However, it is important to emphasize that the strong
coupling limit of the Migdal-Eliashberg theory is unphysical, because the
theory breaks down at finite $\lambda$ before reaching $\lambda=\infty$ as
discussed in Introduction. We will use the strong coupling limit as a
technical tool only, e.g., to show that the breakdown occurs universally,
irrespective of the model electron-phonon Hamiltonian.

Consider the Holstein model first. By definition,

$\lambda=\frac{g^{2}}{\Omega^{2}}=\frac{\nu_{0}\alpha^{2}M^{-1}}{\Omega^{2}},$
(47)

where $\alpha$ is the dimensionful electron-phonon coupling constant. Note
that $\alpha$ does not have units of energy, while $g$ does. We define the
strong coupling limit as $\alpha\to\infty$. This implies $\lambda\to\infty$
and $\alpha,g\propto\lambda^{1/2}$. As $\Omega$ is negligible compared to $g$,
the latter remains the only intrinsic low energy scale in the problem [55].
Other energies characterizing the superconductor are proportional to $g$, for
example, $T_{c}\approx 0.18g$ and the spectral gap is $1.16g$ [6, 30].
Measured in units of $g$ these quantities are finite, while the phonon
frequency goes to zero, $\Omega\propto\lambda^{-1/2}$. Therefore, an
equivalent, but more convenient way to obtain the strong coupling limit is to
keep $g$ fixed and send $\Omega$ to zero.

We already saw that $Z_{n}$ and therefore $\Phi_{n}$ and $\Sigma_{n}$ diverge
when $\lambda\to\infty$, while $F_{n}$ and $G_{n}$ stay finite. Separating the
divergent part in the variable change (10), $\Sigma_{n}=\lambda\pi
TG_{n}+\dots$, $\Phi_{n}=\lambda\pi TF_{n}+\dots$ and expanding the square
root in Eq. (3) in $\Omega^{2}$, we directly obtain the spin Hamiltonian (15)
up to terms of order $\Omega^{2}$ that vanish in the strong coupling limit.
Moreover, it turns out that the mass of fluctuations violating the constraint
(12) is infinite [25] and hence $\bm{S}_{n}^{2}=1$. Note that here we cast the
free energy into the spin chain form without ever using the stationary point
(Eliashberg) equations. Therefore, in the strong coupling limit the spin chain
representation is guaranteed at any point $(G_{n},F_{n})$ in the configuration
space with no stationary point constraints on $G_{n}$ and $F_{n}$. The spin
chain Hamiltonian (17) in the strong coupling limit, $\Omega\to 0$, becomes,

$H_{s}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}Tg^{2}\sum_{nm}\frac{{\bm{S}_{n}}\cdot{\bm{S}_{m}}-1}{(\omega_{n}-\omega_{m})^{2}},$
(48)

where we subtracted a constant from the Hamiltonian as discussed below Eq.
(19).

The procedure for dispersing phonons is the same. Now,

$\begin{split}\lambda(\omega_{l})&=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\\!\\!\\!\\!\frac{g_{q}^{2}qdq}{\omega_{l}^{2}+\omega_{q}^{2}},\\\
\lambda&=\lambda(\omega_{l}=0)=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\frac{g_{q}^{2}qdq}{\omega_{q}^{2}},\end{split}$
(49)

see Appendix B. Suppose $g_{q}\to\infty$ for some range of $q$. Then, the
energy $g$ defined as,

$g^{2}\equiv\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\\!\\!\\!\\!g_{q}^{2}qdq,$
(50)

diverges and all other low energies are either negligible or proportional to
it. This is equivalent to keeping $g$ fixed and sending the phonon frequencies
$\omega_{q}$ to zero. In this limit, the spin chain Hamiltonian (15) for
dispersing phonons turns into the Hamiltonian (48), same as for the Holstein
model. We see that the spin chain (48) provides a complete and universal
description of the the thermodynamics of the strong coupling limit – it
determines the Boltzmann weight of any field configuration and is independent
of the underlying electron-phonon model except for a single energy $g$.

##  VIII Comparison to Anderson pseudospins

We already mentioned in Sect. III that the description of the normal and
superconducting states and of the transition between them in terms of spins
introduced in this paper is very similar to the Anderson pseudospin
description of the BCS theory of superconductivity [26]. In particular, our
Fig. 1 is identical to Fig. 1 of Ref. 26, except that in the latter figure the
cites of the chain are single-fermion energies $\xi_{\bm{p}}$ instead of
fermionic Matsubara frequencies $\omega_{n}$. Indeed, we will see in this
section that equilibrium configurations of Eliashberg and Anderson spins map
into each other under the interchange of $\xi_{\bm{p}}$ with $\omega_{n}$ and
of BCS gap $\Delta$ with the Eliashberg gap function $\Delta_{n}$.

An overall principle common to both spin formulations of superconductivity is
that there are three real functions of energy or frequency – the normal
average and real and imaginary parts of the anomalous average – that admit an
interpretation as three components of a classical spin. Nevertheless,
Eliashberg and Anderson spins are distinct. Eliashberg spins are energy-
integrated normal and anomalous thermal Green’s functions. They exist at any
temperature and coupling $\lambda$. Anderson spins are frequency-integrated
normal and anomalous Keldysh Green’s functions (see below). They are well-
defined both in and out of equilibrium, but only at $T=0$. The two set of
spins do not coincide even where their domains of definition overlap, i.e., at
$T=0$ in the weak coupling (BCS) limit of the Migdal-Eliashberg theory.

###  1 In and out of equilibrium BCS superconductivity in terms of classical
Anderson pseudospins – brief review

Soon after the publication of the BCS theory, Anderson realized [26] that the
BCS Hamiltonian,

$H=\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}-\lambda\delta\sum_{{\bm{p}}{\bm{p}}^{\prime}}c^{\dagger}_{-{\bm{p}}\uparrow}c^{\dagger}_{{\bm{p}}\downarrow}c_{{\bm{p}}^{\prime}\downarrow}c_{-{\bm{p}}^{\prime}\uparrow},$
(51)

where $\delta=(\nu_{0}N)^{-1}$, maps to a spin-1/2 model. We will need quantum
averages of pseudospin-1/2 operators defined by Anderson [56],

$s_{\bm{p}}^{z}=\frac{\langle
c^{\dagger}_{{\bm{p}}\uparrow}c_{{\bm{p}}\uparrow}+c^{\dagger}_{-{\bm{p}}\downarrow}c_{-{\bm{p}}\downarrow}\rangle-1}{2},\quad
s_{\bm{p}}^{-}=\langle c_{{\bm{p}}\downarrow}c_{-{\bm{p}}\uparrow}\rangle,$
(52)

where $s_{\bm{p}}^{\pm}=s_{\bm{p}}^{x}\pm is_{\bm{p}}^{y}$. Note that
$\bm{s}_{\bm{p}}$ depend on time out of equilibrium, when the state of the
system in which we evaluate the averages (52) is itself time-dependent.

In the mean field treatment, which is exact for the BCS model in the
thermodynamic limit [57, 58, 59], the Hamiltonian (51) becomes [61],

$H=\sum_{{\bm{p}}}2\xi_{\bm{p}}s_{\bm{p}}^{z}-\lambda\delta\sum_{{\bm{p}}{\bm{p}}^{\prime}}s_{\bm{p}}^{+}s_{{\bm{p}}^{\prime}}^{-}.$
(53)

The eigenstates of the BCS Hamiltonian are product states of the form
$\prod_{\bm{p}}(u_{\bm{p}}+v_{\bm{p}}c^{\dagger}_{{\bm{p}}\uparrow}c^{\dagger}_{-{\bm{p}}\downarrow})|0\rangle$.
For any state of this form,

$s_{\bm{p}}^{z}=\frac{|v_{\bm{p}}|^{2}-|u_{\bm{p}}|^{2}}{2},\quad
s_{\bm{p}}^{-}=u_{\bm{p}}^{*}v_{\bm{p}}.$ (54)

Therefore, the normalization condition $|u_{\bm{p}}|^{2}+|v_{\bm{p}}|^{2}=1$
implies that $\bm{s}_{\bm{p}}=(s_{\bm{p}}^{x},s_{\bm{p}}^{y},s_{\bm{p}}^{z})$,
is a vector of length $1/2$. Given the spin configuration, we determine the
state of the system using Eq. (54).

Equilibria of the spin Hamiltonian (53) correspond to the eigenstates of the
BCS Hamiltonian [26, 62, 63, 60]. Spin $\bm{s}_{\bm{p}}$ experiences an
effective magnetic field
$\bm{B}_{\bm{p}}=(-2\Delta_{x},-2\Delta_{y},2\xi_{\bm{p}})$, where
$\Delta_{x}$ and $-\Delta_{y}$ are the real and imaginary parts of the BCS
order parameter defined by,

$\Delta=\lambda\delta\sum_{\bm{p}}s_{{\bm{p}}}^{-}.$ (55)

In equilibrium the spin is either antiparallel ($e_{\bm{p}}=+1$) or parallel
($e_{\bm{p}}=-1$) to its effective field. It follows that equilibrium spin
configurations are,

$2s_{\bm{p}}^{z}=-\frac{e_{\bm{p}}\xi_{\bm{p}}}{\sqrt{\smash{\xi_{\bm{p}}^{2}}+|\Delta|^{2}}},\quad
2s_{\bm{p}}^{-}=\frac{e_{\bm{p}}\Delta}{\sqrt{\smash{\xi_{\bm{p}}^{2}}+|\Delta|^{2}}}.$
(56)

The self-consistency condition (55) becomes,

$2\Delta=\lambda\delta\sum_{\bm{p}}\frac{e_{\bm{p}}\Delta}{\sqrt{\smash{\xi_{\bm{p}}^{2}}+|\Delta|^{2}}}.$
(57)

The BCS ground state is obtained by aligning all spins antiparallel to their
fields, i.e., $e_{\bm{p}}=+1$ for all ${\bm{p}}$. States where one of the
spins is flipped opposite to its ground state orientation correspond to
excited or “real” pairs in the terminology of Bardeen, Cooper, and Schrieffer
[20, 26]. Note that Eqs. (56) and (57) describe two types of equilibrium
configurations – configurations with $\Delta=0$ [which is always a solution of
Eq. (57)] and configutrations where $\Delta\neq 0$. In the terminology of Ref.
63, these are the normal and anomalous states, respectively.

Classical Anderson pseudospins play a central role in understanding
collisionless dynamics of BCS condensates in response to fast perturbations
[61, 63, 62, 65, 66, 67, 64, 68] when the sample size is of the order of the
superconducting coherence length or smaller [65]. Hamilton’s equations of
motion for the classical spin Hamiltonian (53) with the usual angular momentum
Poisson brackets for components of $\bm{s}_{\bm{p}}$ are equivalent to the
time-dependent Bogoliubov-de Gennes equations [62]. These classical spin
equations of motion have been used by Anderson to analyze the collective modes
of a BCS superconductor and more recently to study its time evolution after a
quantum quench [61, 63, 66, 67, 64, 68]. Normal and anomalous states mentioned
in the preceding paragraph play a special role in the dynamics [62]. Normal
states are dynamically unstable and anomalous states are unstable when
sufficiently many spins are flipped. These unstable equilibria give rise to
normal and anomalous multi-solitons that start in an unstable equilibrium at
$t=-\infty$ and return into it at $t=+\infty$, similar to an inverted
pendulum. Moreover, dynamics for many physical initial conditions (e.g., for
an interaction quench) can be described in terms of multi-solitons.

Let us also mention the relationship between Anderson pseudospins and equal
time, zero temperature Keldysh Green’s functions defined by,

$\begin{split}\mathsf{G}_{\bm{p}}(t,t^{\prime})=-i\langle[c_{{\bm{p}}\sigma}(t),c_{{\bm{p}}\sigma}^{\dagger}(t^{\prime})]\rangle,\\\
\mathsf{F}_{\bm{p}}(t,t^{\prime})=-i\langle[c_{{\bm{p}}\uparrow}(t),c_{-{\bm{p}}\downarrow}(t^{\prime})]\rangle,\end{split}$
(58)

where the square brackets denote the commutator and we assumed time reversal
symmetry, so that $\mathsf{G}_{\bm{p}}(t,t^{\prime})$ does not depend on
$\sigma$. We see that $G_{\bm{p}}(t,t)=2is_{\bm{p}}^{z}$ and
$F_{\bm{p}}(t,t)=2is_{\bm{p}}^{-}$. Fourier transforming the Green’s functions
with respect to $t^{\prime}-t$, we find

$2is_{\bm{p}}^{z}=\int_{-\infty}^{\infty}\\!\\!\\!\\!d\omega\mathsf{G}_{\bm{p}}(\omega),\quad
2is_{\bm{p}}^{-}=\int_{-\infty}^{\infty}\\!\\!\\!\\!d\omega\mathsf{F}_{\bm{p}}(\omega).$
(59)

We suppressed the dependence on $t$ in spin components,
$\mathsf{G}_{\bm{p}}(\omega)$ and $\mathsf{F}_{\bm{p}}(\omega)$. Thus, we see
that Anderson spins are frequency-integrated normal and anomalous Keldysh
Green’s functions (divided by $2i$).

###  2  Strongly coupled superconductors in and out of equilibrium in terms of
Eliashberg spins

Let us first compare equilibrium configurations of Eliashberg and Anderson
spins. In equilibrium, Eliashberg spins are [Eq. (41)],

$S_{n}^{z}=\frac{e_{n}\omega_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}},\quad
S_{n}^{+}=\frac{e_{n}\Delta_{n}}{\sqrt{\omega_{n}^{2}+|\Delta_{n}|^{2}}}.$
(60)

We see that they turn into Anderson spins (56), up to a factor of 2 and
opposite sign of $y$ and $z$ components, if we replace the Matsubara frequency
$\omega_{n}$ with the single-particle energy $\xi_{\bm{p}}$ and the Eliashberg
gap function $\Delta_{n}$ with the BCS gap $\Delta$. The factor of 2
difference comes from normalization – the arbitrary choice of spin length –
1/2 for Anderson spins and 1 for Eliashberg spins. Opposite signs of $y$ and
$z$ components of spins is a similar “gauge” degree of freedom that does not
affect the Poisson brackets between spin components. It is straightforward to
redefine Anderson or Eliashberg spins to eliminate these differences. However,
the gap equations (57) and (38) are rather different due to the retarded
nature of interactions and the presence of $zz$ interactions in the spin
formulation of the Migdal-Eliashberg theory.

There are also similarities between the two approaches in the overall physical
picture they provide. For example, in both cases the normal state is a state
where all spins are down (up) below a certain energy and up (down) above it as
shown in Fig. 1 and spins acquire $xy$-components in the superconducting state
thus softening the domain wall at zero frequency (energy).

Nevertheless, the two sets of spins never coincide. Anderson spins are ill-
defined away from $\lambda\to 0$, $T=0$ limit, i.e., in the Migdal-Eliashberg
theory at finite $\lambda$ or finite $T$, in the sense that the model cannot
be formulated entirely in terms of them and the length of the vector
$\bm{s}_{\bm{p}}$ depends on ${\bm{p}}$ and is not conserved by the dynamics.
In contrast, Eliashberg spins remain well-defined in the BCS ($\lambda\to 0$)
limit and, in particular, the classical spin Hamiltonian (46) minus its value
in the normal state gives the BCS condensation energy as a function of $T$.
Even in the BCS limit the two sets of spins are distinct and provide
alternative descriptions of the system – Eliashberg spins determine the free
energy as a function of $T$, while Anderson spins work at $T=0$ and determine
the ground state and excited states of the Hamiltonian.

While Anderson spins are frequency-integrated Keldysh Green’s functions [see
Eq. (59)], Eliashberg spins correspond to energy-integrated normal and
anomalous Matsubara Green’s functions,

$\begin{split}\mathcal{G}_{\sigma{\bm{p}}}(\tau-\tau^{\prime})=-\langle
T_{\tau}c_{{\bm{p}}\sigma}(\tau)c^{\dagger}_{{\bm{p}}\sigma}(\tau^{\prime})\rangle,\\\
\mathcal{F}_{\bm{p}}(\tau-\tau^{\prime})=\langle
T_{\tau}c_{-{\bm{p}}\downarrow}(\tau)c_{{\bm{p}}\uparrow}(\tau^{\prime})\rangle.\end{split}$
(61)

We have evaluated these Green’s functions in the Matsubara frequency domain in
Appendix A4,

$\begin{split}\mathcal{G}_{{\bm{p}}n}=-\frac{i(\omega_{n}+\Sigma_{n})+\xi_{\bm{p}}+\chi_{n}}{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\bm{p}})^{2}},\\\
\mathcal{F}_{{\bm{p}}n}=-\frac{\Phi_{n}}{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\bm{p}})^{2}}.\end{split}$
(62)

Integrating $\mathcal{G}_{{\bm{p}}n}$ and $\mathcal{F}_{{\bm{p}}n}$ and
comparing with the expressions (11) for Eliashberg spins on the stationary
point, we see that,

$S_{n}^{z}=-i\int_{-\infty}^{\infty}\\!\\!\\!\\!d\xi_{\bm{p}}\mathcal{G}_{{\bm{p}}n},\quad
S_{n}^{+}=-\int_{-\infty}^{\infty}\\!\\!\\!\\!d\xi_{\bm{p}}\mathcal{F}_{{\bm{p}}n}.$
(63)

Note that the single-particle energy $\xi_{\bm{p}}$ plays the same role for
Eliashberg spins as the frequency $\omega$ does for Anderson spins.

We explained in the previous subsection how Anderson spins describe the far
from equilibrium dynamics of a BCS condensate. Likewise, we expect the time
evolution of Eliashberg spins with the Hamiltonian (46) to describe the
dynamics of electronic degrees of freedom of a strongly coupled (Eliashberg)
superconductor in a similar regime. Hamiltonian (Bloch) equations of motion
for Eliashberg spins read,

$\begin{split}\frac{d\bm{S}_{n}}{dt}&=-\bm{b}_{n}\times\bm{S}_{n},\\\
\bm{b}_{n}&=2\pi\omega_{n}\hat{\bm{z}}+2\pi^{2}T\sum_{m\neq
n}\lambda_{nm}\bm{S}_{m}.\end{split}$ (64)

These equations or a certain version of them should be valid for the short-
time far from equilibrium dynamics, at least when the sample size is of the
order of the coherence length or smaller. In this case, dynamics are spatially
uniform [65], phonons are in thermal equilibrium with the outside environment
[69], and collision integrals are negligible [70]. Indeed, kinetic equations
very similar to Eq. (64) have already been derived by Eliashberg [69], see
also Ref. 70. In the weak coupling limit, $\lambda_{nm}\to\lambda$ at relevant
frequencies and, with the replacement of $\omega_{n}$ with $\xi_{\bm{p}}$, Eq.
(64) is equivalent to the Bloch equation for Anderson spins $\bm{s}_{\bm{p}}$
defined by Eq. (52).

The spin flip solutions we identified in Sect. IV2 play a special role in the
dynamics. These saddle points of free energy are by construction stationary
solutions of the Bloch Eq. (64). Just as with flipped Anderson spins discussed
below Eq. (57), these equilibria are dynamically unstable when sufficiently
many spins are flipped and there are two types of them – normal and anomalous
– depending on whether or not the $xy$-components of spins are zero.
Therefore, we anticipate that these saddle points should give rise to dynamics
analogous to normal and anomalous solitons in the far from equilibrium BCS
superconductivity.

##  IX Conclusion

In this paper, we mapped the Migdal-Eliashberg theory of superconductivity to
a classical Heisenberg spin chain in Zeeman magnetic field. Lattice sites are
fermionic Matsubara frequencies $\omega_{n}$. The interaction between spins is
purely ferromagnetic. It depends on the phonon spectrum and falls off as
$(\omega_{n}-\omega_{m})^{-2}$ at large separation between spins. The Zeeman
field is proportional to $\omega_{n}$ – the spin coordinate along the chain.
As an example, the spin Hamiltonian for the Holstein model is,

$H_{s}=-2\pi\sum_{n}\omega_{n}S_{n}^{z}-\pi^{2}Tg^{2}\sum_{nm}\frac{{\bm{S}_{n}}\cdot{\bm{S}_{m}}-1}{(\omega_{n}-\omega_{m})^{2}+\Omega^{2}},$
(65)

where $T$ is the temperature, $g$ is the electron-phonon coupling, $\Omega$ is
the renormalized Einstein phonon frequency, and $\bm{S}_{n}^{2}=1$. The free
energy density of the system of fermions interacting through phonons is
$f=\nu_{0}TH_{s}$.

The spin chain formulation made the analysis of the free energy functional
simple. In particular, we saw that, up to an overall phase $e^{i\phi}$, the
Eliashberg order parameter $\Delta(\omega_{n})$ is nonnegative and even at the
global minimum and there can be no continuous families of stationary points
that include the global minimum. We also discussed symmetry and non-
integrability arguments that further support our claim that the global minimum
is unique.

It became apparent that the free energy acquires infinitely many new saddle
points (which are new, additional solutions of the Eliashberg equations) at
strong coupling, which correspond to spin flips. We saw that these saddle
points most likely play a significant role in the collisionless far from
equilibrium dynamics of strongly coupled superconductors. Apparently, they are
the fixed points of the corresponding kinetic equations and give rise to rich
multi-soliton dynamics. The new saddle points emerge and proliferate just
before the breakdown of the Migdal-Eliashberg theory and, moreover, these two
phenomena have similar origin as we discussed in Sect. IV2. It is evident that
the spin chain has unstable equilibria where a number of spins are flipped
against their magnetic fields. However, without the spin chain representation
these solutions are hard to notice and indeed they have not been seen before.

Our spin flip solutions of Eliashberg equations are new and unrelated to the
continuous family of sign-alternating solutions for $\Delta_{n}$ conjectured
by Wu et. al. [19] at $\lambda=\infty$. We believe that this conjecture is
internally inconsistent in several ways. Therefore, the claims by the
aforementioned reference of vanishing of the superconducting $T_{c}$,
destruction of the superconducting phase coherence etc. are unsubstantiated.
Thinking in terms of the spin chain also made it clear that any sign-
alternating $\Delta_{n}$ costs a macroscopic amount of energy to create at all
$\lambda$ and therefore cannot contribute to thermal equilibrium properties of
a macroscopic system. Most importantly, because as we discussed above the
Migdal-Eliashberg theory breaks down at finite $\lambda$, its strong coupling
limit is unphysical altogether.

Classical spins in the Migdal-Eliashberg theory (Eliashberg spins) are in many
ways similar to Anderson spins in the BCS theory of superconductivity, though
the two sets of spins never coincide. Common to both spin notions is that
three real numbers – the normal Green’s function and real and imaginary parts
of the anomalous Green’s function integrated over single-fermion energy
(Eliashberg) or frequency (Anderson) – become the three components of a
classical spin vector. Eliashberg spins exist at any temperature and any
$\lambda$, while Anderson spins are well-defined in and out of equilibrium,
but only at $T=0$ and in the BCS limit of the Migdal-Eliashberg theory.
Equilibrium configurations of Anderson and Eliashberg spins map into each
other under the interchange of single-fermion energy $\xi_{\bm{p}}$ with
fermionic Matsubara frequency $\omega_{n}$ and of the BCS gap $\Delta$ with
Eliashberg gap function $\Delta_{n}$. An interesting open question is to study
the short-time far from equilibrium dynamics of strongly coupled conventional
superconductors with the help of Eliashberg spins.

###### Acknowledgements.

We thank A. Abanov, A. V. Chubukov, and M. Kiessling for discussions and
comments.

##  Appendix A Eliashberg free energy for the Holstein model with arbitrary
single-particle potential

In this appendix, we derive the Eliashberg free energy from the path integral
for the Holstein model, see also Refs. 71, 72. For the most part the
derivation follows a standard sequence of steps. One notable difference is
that we keep the single-electron potential arbitrary throughout the entire
calculation and rewrite the action in its eigenbasis. This will be especially
useful later, in our study of the breakdown of the Eliashberg theory at finite
$\lambda$ [23].

###  1 Effective action

The Lagrangian corresponding to the Holstein Hamiltonian (1) is,

$\begin{split}L=\sum_{\bm{i}\bm{j},\sigma}\psi^{*}_{\bm{i}\sigma}G_{0{\bm{i}}{\bm{j}}}^{-1}\psi_{\bm{j}\sigma}+\sum_{\bm{i}}\left[\frac{M\Omega_{0}^{2}\varphi^{2}_{\bm{i}}}{2}+\frac{M(\partial_{\tau}\varphi_{\bm{i}})^{2}}{2}\right]+\alpha\sum_{\bm{i}\sigma}\psi^{*}_{\bm{i}\sigma}\psi_{\bm{i}\sigma}\varphi_{\bm{i}}+\sum_{{\bm{i}}\sigma}(J^{*}_{{\bm{i}}\sigma}\psi_{{\bm{i}}\sigma}+\psi^{*}_{{\bm{i}}\sigma}J_{{\bm{i}}\sigma}),\end{split}$
(66)

where
$G_{0{\bm{i}}{\bm{j}}}^{-1}=\partial_{\tau}\delta_{{\bm{i}}{\bm{j}}}+h_{\bm{i}\bm{j}}-\mu\delta_{{\bm{i}}{\bm{j}}}$,
$\psi_{\bm{i}\sigma}$ and $\psi_{\bm{i}\sigma}^{*}$ are Grassmann variables,
and $\varphi_{\bm{i}}$ is a real field that corresponds to the ion
displacement operator $x_{\bm{i}}$. We also introduced a chemical potential
$\mu$ and Grassmann source fields $J^{*}_{{\bm{i}}\sigma}$ and
$J_{{\bm{i}}\sigma}$, which we will use to evaluate Green’s functions. All
fields in the Lagrangian depend on the imaginary time $\tau$.

Integrating out phonons, we obtain the following effective action:

$\begin{split}S_{\mathrm{eff}}=\sum_{{\bm{i}}{\bm{j}}\tau\sigma}\psi^{*}_{{\bm{i}}\tau\sigma}G_{0{\bm{i}}{\bm{j}}}^{-1}\psi_{{\bm{j}}\tau\sigma}-\frac{1}{2}\sum_{\bm{i}\tau\tau^{\prime}\sigma\sigma^{\prime}}\psi^{*}_{\bm{i}\tau^{\prime}\sigma^{\prime}}\psi_{\bm{i}\tau^{\prime}\sigma^{\prime}}D_{\tau^{\prime}\tau}\psi^{*}_{\bm{i}\tau\sigma}\psi_{\bm{i}\tau\sigma}+\sum_{{\bm{i}}\tau\sigma}\left(J^{*}_{{\bm{i}}\tau\sigma}\psi_{{\bm{i}}\tau\sigma}+\psi^{*}_{{\bm{i}}\tau\sigma}J_{{\bm{i}}\tau\sigma}\right)\\!.\end{split}$
(67)

Summations indicate summation over ${\bm{i}}$ and integration over $\tau$ and
$\tau^{\prime}$,

$D_{\tau^{\prime}\tau}=T\sum_{l}\frac{\alpha^{2}M^{-1}}{\omega_{l}^{2}+\Omega^{2}}e^{i\omega_{l}(\tau^{\prime}-\tau)},$
(68)

is the effective electron-electron interaction, and $\omega_{l}=2\pi Tl$ are
bosonic Matsubara frequencies. The interaction is proportional to the phonon
propagator as usual [3].

We replaced the bare phonon frequency $\Omega_{0}$ with the renormalized
frequency $\Omega$. To obtain an equation for $\Omega$ within the path
integral framework, one needs to introduce an additional Hubbard-Stratonovich
field $\Pi_{\bm{i}}(\tau^{\prime},\tau)$ for phonons [71]. On the stationary
point this leads to the usual phonon renormalization procedure, see Ref. 3, p.
181. We reject it here, because Holstein and other standard electron-phonon
Hamiltonians do not renormalize phonons correctly [73, 74]. Instead, we treat
$\Omega$ as a parameter of the model. Renormalized phonon frequencies are
generally momentum dependent even when the bare spectrum is dispersionless.
Then, the effective action is that for dispersing phonons, which we treat in
Appendix B. Here we disregard the momentum dependence of $\Omega$ for
simplicity – Eq. (67) is a legitimate model of the electron-phonon system in
its own right, no less legitimate than the Holstein model itself.

We decouple the four-fermion term in Eq. (67) with a Hubbard-Stratonovich
transformation in the particle-particle and particle-hole channels with three
fields $\Phi_{\bm{i}}(\tau^{\prime},\tau)$,
$\Sigma_{{\bm{i}}\uparrow}(\tau^{\prime},\tau)$, and
$\Sigma_{{\bm{i}}\downarrow}(\tau^{\prime},\tau)$,

$\begin{split}S_{\mathrm{eff}}=\sum_{\bm{i}}\int\\!\\!\\!\\!\int
d\tau^{\prime}d\tau\left[\frac{\Phi^{*}_{\bm{i}}(\tau^{\prime},\tau)\Phi_{\bm{i}}(\tau^{\prime},\tau)}{D(\tau^{\prime}-\tau)}+\sum_{\sigma}\frac{\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)\Sigma_{{\bm{i}}\sigma}(\tau,\tau^{\prime})}{2D(\tau^{\prime}-\tau)}\right]+\\\
+\sum_{x^{\prime}x}\Psi^{\dagger}_{x^{\prime}}M_{x^{\prime}x}\Psi_{x}+\sum_{x}\left(K_{x}^{\dagger}\Psi_{x}+\Psi_{x}^{\dagger}K_{x}\right)\\!,\end{split}$
(69)

where $\Psi_{x}$ and $K_{x}$ are two-component Nambu fields, $M_{x^{\prime}x}$
is a $2\times 2$ matrix in the Nambu space,

$\Psi_{x}=\begin{bmatrix}\psi_{\uparrow x}\\\ \psi^{*}_{\downarrow x}\\\
\end{bmatrix}\\!,\quad K_{x}=\begin{bmatrix}J_{{\bm{i}}\uparrow}(\tau)\\\
-J_{{\bm{i}}\downarrow}^{*}(\tau)\\\ \end{bmatrix}\\!,\quad
M_{x^{\prime}x}=\begin{bmatrix}G_{0{\bm{i}}{\bm{j}}}^{-1}\delta(\tau^{\prime}-\tau)-i\Sigma_{{\bm{i}}\uparrow}(\tau^{\prime},\tau)\delta_{{\bm{i}}{\bm{j}}}&\Phi_{{\bm{i}}}(\tau^{\prime},\tau)\delta_{{\bm{i}}{\bm{j}}}\\\
\Phi_{{\bm{i}}}^{*}(\tau,\tau^{\prime})\delta_{{\bm{i}}{\bm{j}}}&{\mkern
1.5mu\overline{\mkern-1.5muG\mkern-1.5mu}\mkern
1.5mu}_{0{\bm{i}}{\bm{j}}}^{-1}\delta(\tau^{\prime}-\tau)+i\Sigma_{{\bm{i}}\downarrow}(\tau,\tau^{\prime})\delta_{{\bm{i}}{\bm{j}}}\\\
\end{bmatrix}\\!,$ (70)

and ${\mkern 1.5mu\overline{\mkern-1.5muG\mkern-1.5mu}\mkern
1.5mu}_{0{\bm{i}}{\bm{j}}}^{-1}=\partial_{\tau}\delta_{{\bm{i}}{\bm{j}}}-h_{\bm{j}\bm{i}}+\mu\delta_{{\bm{i}}{\bm{j}}}$.
By definition the fields $\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$ are
Hermitian

$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)=\Sigma_{{\bm{i}}\sigma}^{*}(\tau,\tau^{\prime}).$
(71)

The statement here is that integrating $e^{-S_{\mathrm{eff}}}$ with
$S_{\mathrm{eff}}$ from Eq. (69) over the fields
$\Phi_{\bm{i}}(\tau^{\prime},\tau)$ and
$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$, we obtain Eq. (67). On the
stationary point, the fields $\Phi_{\bm{i}}(\tau^{\prime},\tau)$ and
$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$ depend only on the difference
$\tau^{\prime}-\tau$ and represent the anomalous and normal self-energies (see
below).

The action (69) is quadratic in fermion fields. Performing the Gaussian
integral over these fields, we obtain the effective action for the fields
$\Phi_{\bm{i}}(\tau^{\prime},\tau)$ and
$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$,

$\begin{split}S_{\mathrm{eff}}=\sum_{\bm{i}}\int\\!\\!\\!\\!\int
d\tau^{\prime}d\tau\left[\frac{\Phi^{*}_{\bm{i}}(\tau^{\prime},\tau)\Phi_{\bm{i}}(\tau^{\prime},\tau)}{D(\tau^{\prime}-\tau)}+\sum_{\sigma}\frac{\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)\Sigma_{{\bm{i}}\sigma}(\tau,\tau^{\prime})}{2D(\tau^{\prime}-\tau)}\right]-\mathrm{Tr}\ln
M-\sum_{x^{\prime}x}K^{\dagger}_{x^{\prime}}M^{-1}_{x^{\prime}x}K_{x},\end{split}$
(72)

The action is a scalar – matrix trace. It is convenient to transform the
matrices and vectors in $S_{\mathrm{eff}}$ to fermionic Matsubara frequencies
$\omega_{n}$ and the eigenbasis of the single-electron Hamiltonian $\hat{h}$
with a unitary matrix,

$U_{\tau
n,{\bm{i}}\alpha}=\sqrt{T}\begin{bmatrix}\pi_{{\bm{i}}\alpha}e^{-i\tau\omega_{n}}&0\\\
0&\pi_{{\bm{i}}\alpha}e^{-i\tau\omega_{n}}\\\ \end{bmatrix}\\!,$ (73)

which leaves $S_{\mathrm{eff}}$ invariant. Here $\pi_{\alpha}$ are the
eigenstates of $\hat{h}$,

$\sum_{\bm{j}}h_{{\bm{i}}{\bm{j}}}\pi_{{\bm{j}}\alpha}=\epsilon_{\alpha}\pi_{{\bm{i}}\alpha}.$
(74)

For example, the matrix $M$ transforms into
$(U^{\dagger}MU)_{nm,\alpha\beta}=\sum_{{\bm{i}}{\bm{j}}}\int\\!\\!\int
d\tau^{\prime}d\tau
U^{\dagger}_{n\tau,\alpha{\bm{i}}}M_{{\bm{i}}\tau,{\bm{i}}^{\prime}\tau^{\prime}}U_{\tau^{\prime}m,{\bm{i}}^{\prime}\beta}$
and similarly $K\to U^{\dagger}K$ and $\Psi\to U^{\dagger}\Psi$. Carrying out
this transformation in Eq. (72), we find,

$\begin{split}{S}_{\mathrm{eff}}=T\\!\\!\\!\\!\sum_{nml\alpha\beta}\left[\left(\Phi^{\alpha\beta}_{n+l,m+l}\right)^{*}D^{-1}_{l}\Phi_{nm}^{\alpha\beta}+\frac{1}{2}\sum_{\sigma}\Sigma_{\sigma,m+l,n+l}^{\alpha\beta}D^{-1}_{l}\Sigma_{\sigma
nm}^{\beta\alpha}\right]-\mathrm{Tr}\,\\!\ln
M-TK^{\dagger}M^{-1}K.\end{split}$ (75)

Here $D_{l}^{-1}$ is the bosonic Matsubara frequency $\omega_{l}$ component of
$1/D(\tau)$,

$\Phi_{nm}^{\alpha\beta}=T\sum_{\bm{i}}\\!\int\\!\\!\\!\int
d\tau^{\prime}d\tau
e^{i\omega_{n}\tau^{\prime}}\pi^{*}_{{\bm{i}}\alpha}\Phi_{\bm{i}}(\tau^{\prime},\tau)\pi_{{\bm{i}}\beta}e^{-i\omega_{m}\tau},$
(76)

and similarly for $\Sigma_{\sigma nm}$. Hermiticity of
$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$ implies Hermiticity of
$\Sigma_{\sigma nm}^{\alpha\beta}$, since the transformation is unitary. The
matrix $M$ in the new basis reads,

$M=\begin{bmatrix}(-i\omega_{n}+\xi_{\alpha})\delta_{nm}-i\Sigma_{\uparrow
nm}^{\alpha\beta}&\Phi_{nm}^{\alpha\beta}\\\
\left(\Phi_{mn}^{\beta\alpha}\right)^{*}&(-i\omega_{n}-\xi_{\alpha})\delta_{nm}+i\Sigma_{\downarrow,-m,-n}^{\alpha\beta}\\\
\end{bmatrix},\\\ $ (77)

and the source fields are,

$\begin{split}K_{\alpha m}=\begin{bmatrix}J_{\uparrow\alpha m}\\\
-J_{\downarrow,\alpha^{\prime},-m}^{*}\\\ \end{bmatrix},\quad
J_{\uparrow\alpha m}=\sum_{\bm{i}}\int d\tau
e^{i\omega_{m}\tau}\pi_{{\bm{i}}\alpha}^{*}J_{{\bm{i}}\uparrow}(\tau),\quad
J_{\downarrow\alpha^{\prime}m}=\sum_{\bm{i}}\int d\tau
e^{i\omega_{m}\tau}(\pi_{{\bm{i}}\alpha}^{*})^{*}J_{{\bm{i}}\downarrow}(\tau).\end{split}$
(78)

Note that $\alpha^{\prime}$ labels the state $\pi_{\alpha}^{*}$, which is
related to $\pi_{\alpha}$ by time reversal operation and $A_{-m,-n}$ stands
for $A(-\omega_{m},-\omega_{n})$.

###  2 Stationary point

The Eliashberg equations are a stationary point of the effective action, such
that the fields $\Phi_{\bm{i}}(\tau^{\prime},\tau)$ and
$\Sigma_{{\bm{i}}\sigma}(\tau^{\prime},\tau)$ are spatially uniform and depend
only on the time difference,

$\Phi_{\bm{i}}(\tau^{\prime},\tau)=\Phi(\tau^{\prime}-\tau),\quad\Sigma_{{\bm{i}}\sigma}(\tau^{\prime}-\tau)=\Sigma_{\sigma}(\tau^{\prime}-\tau).$
(79)

Eq. (76) then implies,

$\Phi_{nm}^{\alpha\beta}=\Phi_{n}\delta_{nm}\delta_{\alpha\beta},\quad\Sigma_{\sigma
nm}^{\alpha\beta}=\Sigma_{\sigma n}\delta_{nm}\delta_{\alpha\beta}.$ (80)

In words, offdiagonal matrix elements are zero and diagonal ones depend only
on the Matsubara frequency.

The derivative of the action (75) with respect to any offdiagonal matrix
element is a sum of terms each of which is proportional to one of the other
offdiagonal matrix elements. This means that setting these matrix elements to
zero automatically solves all stationary point equations for them. To
determine stationary point equations for the diagonal matrix elements, we
substitute Eq. (80) into Eq. (75). Setting also the sources to zero, we
obtain,

$\begin{split}{S}_{\mathrm{eff}}=TN\sum_{nl}\left[\Phi_{n+l}^{*}D_{l}^{-1}\Phi_{n}+\Sigma_{n+l}D_{l}^{-1}\Sigma_{n}-\chi_{n+l}D_{l}^{-1}\chi_{n}\right]-\sum_{n\alpha}\ln\left[(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\alpha})^{2}\right],\end{split}$
(81)

where $N$ is the number of sites (which is equal to the number of the single-
electron states $\pi_{\alpha}$) and

$\Sigma_{n}=\frac{\Sigma_{\uparrow n}-\Sigma_{\downarrow,-n}}{2},\quad
i\chi_{n}=\frac{\Sigma_{\uparrow n}+\Sigma_{\downarrow,-n}}{2},$ (82)

where $A_{-n}$ means $A(-\omega_{n})$. Minimizing the effective action with
respect to $\Phi_{n}^{*},\Sigma_{n}$ and $\chi_{n}$, we derive three
Eliashberg equations,

$\begin{split}\Phi_{n}=\frac{T}{N}\sum_{m\alpha}D_{n-m}\frac{\Phi_{m}}{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}+(\chi_{m}+\xi_{\alpha})^{2}},\\\
\Sigma_{n}=\frac{T}{N}\sum_{m\alpha}D_{n-m}\frac{\omega_{m}+\Sigma_{m}}{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}+(\chi_{m}+\xi_{\alpha})^{2}},\\\
\chi_{n}=-\frac{T}{N}\sum_{m\alpha}D_{n-m}\frac{\xi_{\alpha}+\chi_{m}}{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}+(\chi_{m}+\xi_{\alpha})^{2}},\end{split}$
(83)

where

$D_{n-m}=\frac{\alpha^{2}M^{-1}}{(\omega_{n}-\omega_{m})^{2}+\Omega^{2}}.$
(84)

When $\xi_{\alpha}$ are symmetric with respect to zero, $\chi_{n}\equiv 0$
solves the last equation in (83).

As is normally done in the Migdal-Eliashberg theory, we now send the Fermi
energy to infinity and take the density of states per site per spin
orientation, $\nu_{0}$, to be constant. We implement this limit in Eq. (81) by
integrating the log over $\xi_{\alpha}$ from $-\mathbb{\Lambda}$ to
$\mathbb{\Lambda}$, discarding the constant term that depends only on
$\mathbb{\Lambda}$, and then taking the limit $\mathbb{\Lambda}\to\infty$. The
result is

$\begin{split}{S}_{\mathrm{eff}}=TN\sum_{nl}\left[\Phi_{n+l}^{*}D_{l}^{-1}\Phi_{n}+\Sigma_{n+l}D_{l}^{-1}\Sigma_{n}\right]-2\pi\nu_{0}N\sum_{n}\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}.\end{split}$
(85)

Eliashberg equations become,

$\begin{split}\Phi_{n}=\pi
T\sum_{m}\lambda(\omega_{n}-\omega_{m})\frac{\Phi_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}},\\\
\Sigma_{n}=\pi
T\sum_{m}\lambda(\omega_{n}-\omega_{m})\frac{\omega_{m}+\Sigma_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}}.\end{split}$
(86)

with

$\lambda(\omega_{n}-\omega_{m})=\frac{g^{2}}{(\omega_{n}-\omega_{m})^{2}+\Omega^{2}},\quad
g^{2}=\nu_{0}\alpha^{2}M^{-1}.$ (87)

The fermionic part of the free energy functional per site is
$f=TS_{\mathrm{eff}}/N$,

$\begin{split}f=\nu_{0}T^{2}\sum_{nl}\left[\Phi_{n+l}^{*}\Lambda_{l}\Phi_{n}+\Sigma_{n+l}\Lambda_{l}\Sigma_{n}\right]-2\pi\nu_{0}T\sum_{n}\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}},\end{split}$
(88)

where $\Lambda_{l}$ is the discrete Fourier transform of $1/\lambda(\tau)$ at
bosonic Matsubara frequency $\omega_{l}$.

###  3 Spatially nonuniform stationary point

For future reference let us also consider spatially nonuniform solutions,

$\Phi_{\bm{i}}(\tau^{\prime},\tau)=\Phi_{\bm{i}}(\tau^{\prime}-\tau),\quad\Sigma_{{\bm{i}}\sigma}(\tau^{\prime}-\tau)=\Sigma_{{\bm{i}}\sigma}(\tau^{\prime}-\tau).$
(89)

Let $\Phi_{{\bm{i}}n}$ be the Fourier transform of
$\Phi_{\bm{i}}(\tau^{\prime}-\tau)$ with respect to $\tau^{\prime}-\tau$. Its
matrix elements in the eigenbasis of the single-particle Hamiltonian
$h_{{\bm{i}}{\bm{j}}}$ are
$\Phi_{n}^{\alpha\beta}=\sum_{\bm{i}}\pi^{*}_{{\bm{i}}\alpha}\Phi_{{\bm{i}}n}\pi_{{\bm{i}}\beta}$.
If the eigenstates $\pi_{{\bm{i}}\alpha}$ are localized or delocalized but
highly oscillatory, the offdiagonal, $\alpha\neq\beta$, matrix elements are
negligible. In this case, a suitable ansatz for the stationary point is,

$\Phi_{nm}^{\alpha\beta}=\Phi_{n}^{\alpha}\delta_{nm}\delta_{\alpha\beta},\quad\Sigma_{\sigma
nm}^{\alpha\beta}=\Sigma_{\sigma n}^{\alpha}\delta_{nm}\delta_{\alpha\beta}.$
(90)

The derivation of the stationary point equations is similar to the uniform
case. The effective action now reads,

$\begin{split}{S}_{\mathrm{eff}}=T\sum_{nl\alpha}\left[(\Phi_{n+l}^{\alpha})^{*}D_{l}^{-1}\Phi_{n}^{\alpha}+\Sigma_{n+l}^{\alpha}D_{l}^{-1}\Sigma_{n}^{\alpha}-\chi_{n+l}^{\alpha}D_{l}^{-1}\chi_{n}^{\alpha}\right]-\sum_{n\alpha}\ln\left[(\omega_{n}+\Sigma_{n}^{\alpha})^{2}+|\Phi_{n}^{\alpha}|^{2}+(\chi_{n}^{\alpha}+\xi_{\alpha})^{2}\right],\end{split}$
(91)

where

$\Sigma_{n}^{\alpha}=\frac{\Sigma_{\uparrow
n}^{\alpha}-\Sigma_{\downarrow,-n}^{\alpha}}{2},\quad
i\chi_{n}^{\alpha}=\frac{\Sigma_{\uparrow
n}^{\alpha}+\Sigma_{\downarrow,-n}^{\alpha}}{2}.$ (92)

It is incorrect to minimize ${S}_{\mathrm{eff}}$ with respect to
$\Phi_{n}^{\alpha}$, because
$\Phi_{n}^{\alpha}=\sum_{\bm{i}}|\pi_{{\bm{i}}\alpha}|^{2}\Phi_{{\bm{i}}n}$
are not independent. The original independent variables are $\Phi_{{\bm{i}}n}$
and the matrix $R_{{\bm{i}}\alpha}=|\pi_{{\bm{i}}\alpha}|^{2}$ that relates
$\Phi_{n}^{\alpha}$ to $\Phi_{{\bm{i}}n}$ is degenerate. For example, the
completeness relation
$\sum_{\alpha}\pi_{{\bm{i}}\alpha}\pi_{{\bm{j}}\alpha}=\delta_{{\bm{i}}{\bm{j}}}$
implies $\sum_{\alpha}R_{{\bm{i}}\alpha}=1$ meaning that the columns of the
matrix $R$ are linearly dependent. Similarly, the normalization condition
$\sum_{\bm{i}}|\pi_{{\bm{i}}\alpha}|^{2}=1$ says that the raws of $R$ are
linearly dependent too. The same applies to $\Sigma_{n}^{\alpha}$ and
$\chi_{n}^{\alpha}$.

We should instead minimize Eq. (91) with respect to $\Phi_{{\bm{i}}n}$,
$\Sigma_{{\bm{i}}n}$, and $\chi_{{\bm{i}}n}$ using

$\frac{\partial\Phi_{n}^{\alpha}}{\partial\Phi_{{\bm{i}}n}}=|\pi_{{\bm{i}}\alpha}|^{2},\mbox{
etc.}$ (93)

We obtain,

$\begin{split}\sum_{\alpha}\Phi_{n}^{\alpha}|\pi_{{\bm{i}}\alpha}|^{2}=T\sum_{m\alpha}D_{n-m}\frac{\Phi_{m}^{\alpha}|\pi_{{\bm{i}}\alpha}|^{2}}{(\omega_{m}+\Sigma_{m}^{\alpha})^{2}+|\Phi_{m}^{\alpha}|^{2}+(\chi_{m}^{\alpha}+\xi_{\alpha})^{2}},\\\
\sum_{\alpha}\Sigma_{n}^{\alpha}|\pi_{{\bm{i}}\alpha}|^{2}=T\sum_{m\alpha}D_{n-m}\frac{(\omega_{m}+\Sigma_{m}^{\alpha})|\pi_{{\bm{i}}\alpha}|^{2}}{(\omega_{m}+\Sigma_{m}^{\alpha})^{2}+|\Phi_{m}^{\alpha}|^{2}+(\chi_{m}^{\alpha}+\xi_{\alpha})^{2}},\\\
\sum_{\alpha}\chi_{n}^{\alpha}|\pi_{{\bm{i}}\alpha}|^{2}=-T\sum_{m\alpha}D_{n-m}\frac{(\xi_{\alpha}+\chi_{m}^{\alpha})|\pi_{{\bm{i}}\alpha}|^{2}}{(\omega_{m}+\Sigma_{m}^{\alpha})^{2}+|\Phi_{m}^{\alpha}|^{2}+(\chi_{m}^{\alpha}+\xi_{\alpha})^{2}}.\end{split}$
(94)

If the single-particle Hamiltonian $h_{{\bm{i}}{\bm{j}}}$ has the periodicity
of the lattice, $\pi_{{\bm{i}}\alpha}$ are plane waves and
$|\pi_{{\bm{i}}\alpha}|^{2}=N^{-1}$. Then, $\Phi_{n}^{\alpha}=\Phi_{n}$,
$\Sigma_{n}^{\alpha}=\Sigma_{n}$, $\chi_{n}^{\alpha}=\chi_{n}$, and we recover
the usual Eliashberg equations (83).

###  4 Green’s functions

Let us evaluate the Fourier transforms $\mathcal{G}_{\sigma\alpha n}$ and
$\mathcal{F}_{\alpha n}$ of the normal and anomalous thermal Green’s functions
at the stationary point. The latter are defined as,

$\mathcal{G}_{\sigma\alpha}(\tau-\tau^{\prime})=-\langle
T_{\tau}c_{\sigma\alpha}(\tau)c^{\dagger}_{\sigma\alpha}(\tau^{\prime})\rangle,\quad\mathcal{F}_{\alpha}(\tau-\tau^{\prime})=\langle
T_{\tau}c_{\downarrow\alpha^{\prime}}(\tau)c_{\uparrow\alpha}(\tau^{\prime})\rangle.$
(95)

We obtain $\mathcal{G}_{\sigma\alpha n}$ and $\mathcal{F}_{\alpha n}$ by
differentiating the partition function ${\cal Z}$ with respect to the source
fields,

$T\mathcal{G}_{\sigma\alpha n}=\frac{1}{\cal Z}\frac{\partial^{2}{\cal
Z}}{\partial J_{\sigma\alpha n}\partial J^{*}_{\sigma\alpha n}},\quad
T\mathcal{F}_{\alpha n}=\frac{1}{\cal Z}\frac{\partial^{2}{\cal Z}}{\partial
J^{*}_{\downarrow\alpha^{\prime}n}\partial J^{*}_{\uparrow\alpha,-n}}.$ (96)

The source fields enter through the $K^{\dagger}M^{-1}K$ term in the effective
action (75). At the stationary point this term becomes,

$\begin{split}K^{\dagger}M^{-1}K=\sum_{\alpha
n}K_{n\alpha}^{\dagger}M^{-1}_{\alpha n}K_{n\alpha}=\sum_{\alpha
n}\frac{1}{\Theta_{\alpha n}}\left[\hbox
to0.0pt{$\displaystyle\phantom{\prod}$\hss}a^{+}_{\alpha
n}J^{*}_{A}J_{A}-\Phi_{n}J^{*}_{A}J^{*}_{A^{\prime}}-(\Phi_{n})^{*}J_{A^{\prime}}J_{A}+a^{-}_{\alpha
n}J_{A^{\prime}}J^{*}_{A^{\prime}}\right],\end{split}$ (97)

where $a^{\pm}_{\alpha n}=i(\omega_{n}+\Sigma_{n})\pm(\xi_{\alpha}+\chi_{n})$,
the labels $A$ and $A^{\prime}$ stand for $\uparrow\\!\\!\alpha n$ and
$\downarrow\\!\alpha^{\prime},\\!-n$, respectively, and $\Theta_{\alpha
n}=(\omega_{n}+\Sigma_{n}^{\alpha})^{2}+|\Phi_{n}^{\alpha}|^{2}+(\chi_{n}^{\alpha}+\xi_{\alpha})^{2}=-\det
M_{\alpha n}$.

The source-dependent part of the partition function at the stationary point is
${\cal Z}_{s}=e^{-TK^{\dagger}M^{-1}K}$. Using Eq. (96), we find

$\begin{split}\mathcal{G}_{\uparrow\alpha
n}=-\frac{i(\omega_{n}+\Sigma_{n})+\xi_{\alpha}+\chi_{n}}{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\alpha})^{2}},\\\
\mathcal{G}_{\downarrow\alpha
n}=-\frac{i(\omega_{n}-\Sigma_{-n})+\xi_{\alpha}+\chi_{-n}}{(\omega_{n}-\Sigma_{-n})^{2}+|\Phi_{-n}|^{2}+(\chi_{-n}+\xi_{\alpha})^{2}},\\\
\mathcal{F}_{\alpha
n}=-\frac{\Phi_{-n}}{(\omega_{n}-\Sigma_{-n})^{2}+|\Phi_{-n}|^{2}+(\chi_{-n}+\xi_{\alpha})^{2}}.\end{split}$
(98)

Suppose there is a time reversal symmetry so that $\Sigma_{\uparrow
n}=\Sigma_{\downarrow n}$ and $\mathcal{G}_{\uparrow\alpha
n}=\mathcal{G}_{\downarrow\alpha n}$. Eq. (82) then implies that $\Sigma_{n}$
is odd in frequency and $\chi_{n}$ is even. It further follows from the above
formulas for $\mathcal{G}_{\uparrow\alpha n}$ and
$\mathcal{G}_{\downarrow\alpha n}$ that $|\Phi_{n}|$ is even and we arrive at
the usual expressions for the thermal Green’s functions in the Migdal-
Eliashberg theory,

$\begin{split}\mathcal{G}_{\uparrow\alpha n}=\mathcal{G}_{\downarrow\alpha
n}\equiv\mathcal{G}_{\alpha
n}=-\frac{i(\omega_{n}+\Sigma_{n})+\xi_{\alpha}+\chi_{n}}{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\alpha})^{2}},\\\
\mathcal{F}_{\alpha
n}=-\frac{\Phi_{n}}{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}+(\chi_{n}+\xi_{\alpha})^{2}}.\end{split}$
(99)

##  Appendix B Path integral formulation of the Migdal-Eliashberg theory for
dispersing phonons

Here we derive Eliashberg equations and free energy for arbitrary phonon
spectrum and momentum-dependent electron-phonon coupling. We will see that
these quantities have the same form as for the Holstein model, but with a more
general kernel which now involves an integral over the phonon spectrum.

Our starting point is the standard electron-phonon Hamiltonian,

$\begin{split}H=\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}+\sum_{\bm{q}}{\omega_{0}({\bm{q}})}b^{\dagger}_{\bm{q}}b_{\bm{q}}+\frac{1}{\sqrt{N}}\sum_{{\bm{p}}{\bm{q}}\sigma}\frac{\alpha_{{\bm{q}}}}{\sqrt{2M\omega_{0}({\bm{q}})}}c^{\dagger}_{{\bm{p}}+{\bm{q}}\sigma}c_{{\bm{p}}\sigma}\left[b^{\dagger}_{-{\bm{q}}}+b_{\bm{q}}\right],\end{split}$
(100)

Hermitian property requires $\alpha_{\bm{q}}=\alpha_{-{\bm{q}}}^{*}$. For
simplicity, we take $\alpha_{\bm{q}}$ and $\omega_{0}({\bm{q}})$ to depend
only on the magnitude of the momentum ${\bm{q}}$. Integrating out phonons, we
obtain an effective action for the fermionic fields,

$\begin{split}S_{\mathrm{eff}}=T\sum_{\mathsf{p}\sigma}\psi^{*}_{\mathsf{p}\sigma}G_{0\mathsf{p}}^{-1}\psi_{\mathsf{p}\sigma}-\frac{1}{2}\frac{T^{3}}{N}\sum_{\mathsf{p}\mathsf{p}^{\prime}\mathsf{q}\sigma\sigma^{\prime}}\frac{|\alpha_{q}|^{2}M^{-1}}{\omega_{m}^{2}+\omega_{q}^{2}}\psi^{*}_{\mathsf{p}+\mathsf{q}\sigma}\psi_{\mathsf{p}\sigma}\psi^{*}_{\mathsf{p}^{\prime}\\!-\mathsf{q}\sigma^{\prime}}\psi_{\mathsf{p}^{\prime}\sigma^{\prime}},\end{split}$
(101)

where $G_{0\mathsf{p}}^{-1}=-i\omega_{n}+\xi_{\bm{p}}$, and
$\mathsf{p}=(\omega_{n},{\bm{p}})$; $\mathsf{q}=(\omega_{m},{\bm{q}})$ are the
frequency-momentum 4-vectors. As in Eq. (68) we replace the bare phonon
frequencies $\omega_{0}(q)$ with the renormalized frequencies $\omega_{q}$ in
the denominator of the effective electron-electron interaction.

It is helpful to rewrite the interaction as,

$\begin{split}\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\frac{T^{2}}{N}\frac{|\alpha_{q}|^{2}M^{-1}}{\omega_{n}^{2}+\omega_{q}^{2}}\psi^{*}_{\mathsf{p}+\mathsf{q}\sigma}\psi_{\mathsf{p}\sigma}\psi^{*}_{\mathsf{p}^{\prime}-\mathsf{q}\sigma^{\prime}}\psi_{\mathsf{p}^{\prime}\sigma^{\prime}}=\psi^{*}_{\mathsf{p}_{1}^{\prime}\uparrow}\psi^{*}_{-\mathsf{p}_{2}^{\prime}\downarrow}D^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\psi_{-\mathsf{p}_{2}\downarrow}\psi_{\mathsf{p}_{1}\uparrow}-\\\
\frac{1}{2}\psi^{*}_{\mathsf{p}_{1}^{\prime}\uparrow}\psi_{\mathsf{p}_{2}^{\prime}\uparrow}D^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\psi^{*}_{\mathsf{p}_{2}\uparrow}\psi_{\mathsf{p}_{1}\uparrow}-\frac{1}{2}\psi^{*}_{\mathsf{p}_{1}^{\prime}\downarrow}\psi_{\mathsf{p}_{2}^{\prime}\downarrow}D^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\psi^{*}_{\mathsf{p}_{2}\downarrow}\psi_{\mathsf{p}_{1}\downarrow},\end{split}$
(102)

where we treat $D$ as a matrix with matrix elements,

$D^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}=D_{\mathsf{q}}\delta_{\mathsf{p}_{1}^{\prime},\mathsf{p}_{1}+q}\delta_{\mathsf{p}_{2}^{\prime},\mathsf{p}_{2}+\mathsf{q}},\quad
D_{\mathsf{q}}=\frac{T^{2}}{N}\frac{|\alpha_{q}|^{2}M^{-1}}{\omega_{n}^{2}+\omega_{q}^{2}}.$
(103)

As before, we decouple the interaction with fields $\Phi$ and
$\Sigma_{\sigma}$,

$\begin{split}S_{\mathrm{eff}}=T\\!\\!\\!\\!\\!\sum_{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}\mathsf{p}_{1}\mathsf{p}_{2}}\left\\{\Phi_{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}^{*}[D^{-1}]^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\Phi_{\mathsf{p}_{1}\mathsf{p}_{2}}+\frac{1}{2}\sum_{\sigma}\Sigma_{\sigma\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}^{*}[D^{-1}]^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\Sigma_{\sigma
p_{1}p_{2}}\right\\}+T\sum_{\mathsf{p}_{1}\mathsf{p}_{2}}\Psi_{\mathsf{p}_{1}}^{\dagger}M_{\mathsf{p}_{1}\mathsf{p}_{2}}\Psi_{\mathsf{p}_{2}},\end{split}$
(104)

where

$M_{\mathsf{p}_{1}\mathsf{p}_{2}}=\begin{pmatrix}[-i\omega_{n}+\xi_{{\bm{p}}_{1}}]\delta_{\mathsf{p}_{1}\mathsf{p}_{2}}-i\Sigma_{\uparrow\mathsf{p}_{1}\mathsf{p}_{2}}&\Phi_{\mathsf{p}_{1}\mathsf{p}_{2}}\\\
\Phi_{\mathsf{p}_{2}\mathsf{p}_{1}}^{*}&[-i\omega_{n}-\xi_{{\bm{p}}_{1}}]\delta_{\mathsf{p}_{1}\mathsf{p}_{2}}+i\Sigma_{\downarrow,-\mathsf{p}_{2},-\mathsf{p}_{1}}\\\
\end{pmatrix},\quad\Psi_{p}=\begin{pmatrix}\psi_{\mathsf{p}\uparrow}\\\
\psi^{*}_{-\mathsf{p}\downarrow}\\\ \end{pmatrix}.$ (105)

Integrating out the fermions, we find

$S_{\mathrm{eff}}=T\\!\\!\\!\\!\\!\sum_{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}\mathsf{p}_{1}\mathsf{p}_{2}}\left\\{\Phi_{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}^{*}[D^{-1}]^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\Phi_{\mathsf{p}_{1}\mathsf{p}_{2}}+\frac{1}{2}\sum_{\sigma}\Sigma_{\sigma\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}^{*}[D^{-1}]^{\mathsf{p}_{1}^{\prime}\mathsf{p}_{2}^{\prime}}_{\mathsf{p}_{1}\mathsf{p}_{2}}\Sigma_{\sigma
p_{1}p_{2}}\right\\}-\mathrm{Tr}\,\\!\ln M.$ (106)

###  1 Stationary point

It is reasonable to expect that in a translationally invariant system, the
stationary point is also translationally invariant. Translational invariance
means that matrices $\Phi_{\mathsf{p}_{1}\mathsf{p}_{2}}$ and
$\Sigma_{\sigma\mathsf{p}_{1}\mathsf{p}_{2}}$ are diagonal,

$\Phi_{\mathsf{p}_{1}\mathsf{p}_{2}}=\Phi_{\mathsf{p}_{1}}\delta_{\mathsf{p}_{1}\mathsf{p}_{2}},\quad\Sigma_{\sigma\mathsf{p}_{1}\mathsf{p}_{2}}=\Sigma_{\sigma\mathsf{p}_{1}}\delta_{\mathsf{p}_{1}\mathsf{p}_{2}},$
(107)

since $\mathsf{p}_{1}=\mathsf{p}_{2}$ ensures that
$e^{i\mathsf{p}_{1}x_{1}-i\mathsf{p}_{2}x_{2}}=e^{i\mathsf{p}_{1}(x_{1}-x_{2})}$
depends only on ${\bm{r}}_{1}-{\bm{r}}_{2}$ and $\tau_{1}-\tau_{2}$.

Setting the offdiagonal matrix elements to zero in the effective action (106),
we obtain

$\begin{split}{S}_{\mathrm{eff}}=T\sum_{\mathsf{p}^{\prime}\mathsf{p}}\left\\{\Phi_{\mathsf{p}^{\prime}}^{*}[D^{-1}]^{\mathsf{p}^{\prime}}_{\mathsf{p}}\Phi_{\mathsf{p}}+\Sigma_{\mathsf{p}^{\prime}}[D^{-1}]^{\mathsf{p}^{\prime}}_{\mathsf{p}}\Sigma_{\mathsf{p}}-\chi_{\mathsf{p}^{\prime}}[D^{-1}]^{\mathsf{p}^{\prime}}_{\mathsf{p}}\chi_{\mathsf{p}}\right\\}-\sum_{\mathsf{p}}\ln\left[(\omega_{n}+\Sigma_{\mathsf{p}})^{2}+|\Phi_{p}|^{2}+(\chi_{\mathsf{p}}+\xi_{\bm{p}})^{2}\right],\end{split}$
(108)

where

$\Sigma_{\mathsf{p}}=\frac{\Sigma_{\uparrow\mathsf{p}}-\Sigma_{\downarrow,-\mathsf{p}}}{2},\quad
i\chi_{\mathsf{p}}=\frac{\Sigma_{\uparrow\mathsf{p}}+\Sigma_{\downarrow,-\mathsf{p}}}{2}.$
(109)

Now we evaluate the stationary point of this effective action with respect to
$\Phi_{\mathsf{p}}^{*}$, $\Sigma_{\mathsf{p}}$ and $\chi_{\mathsf{p}}$ and
multiply the resulting expression by matrix $D$ on both sides. The result is,

$\begin{split}\Phi_{n{\bm{p}}}=\frac{T}{N}\sum_{m{\bm{q}}}\frac{|\alpha_{q}|^{2}M^{-1}}{(\omega_{n}-\omega_{m})^{2}+\omega_{q}^{2}}\frac{\Phi_{m{\bm{p}}^{\prime}}}{(\omega_{m}+\Sigma_{m{\bm{p}}^{\prime}})^{2}+|\Phi_{m{\bm{p}}^{\prime}}|^{2}+(\chi_{m{\bm{p}}^{\prime}}+\xi_{{\bm{p}}^{\prime}})^{2}},\\\
\Sigma_{n{\bm{p}}}=\frac{T}{N}\sum_{m{\bm{q}}}\frac{|\alpha_{q}|^{2}M^{-1}}{(\omega_{n}-\omega_{m})^{2}+\omega_{q}^{2}}\frac{\omega_{m}+\Sigma_{m{\bm{p}}^{\prime}}}{(\omega_{m}+\Sigma_{m{\bm{p}}^{\prime}})^{2}+|\Phi_{m{\bm{p}}^{\prime}}|^{2}+(\chi_{m{\bm{p}}^{\prime}}+\xi_{{\bm{p}}^{\prime}})^{2}},\\\
\chi_{n{\bm{p}}}=-\frac{T}{N}\sum_{m{\bm{q}}}\frac{|\alpha_{q}|^{2}M^{-1}}{(\omega_{n}-\omega_{m})^{2}+\omega_{q}^{2}}\frac{\chi_{m{\bm{p}}^{\prime}}+\xi_{{\bm{p}}^{\prime}}}{(\omega_{m}+\Sigma_{m{\bm{p}}^{\prime}})^{2}+|\Phi_{m{\bm{p}}^{\prime}}|^{2}+(\chi_{m{\bm{p}}^{\prime}}+\xi_{{\bm{p}}^{\prime}})^{2}},\end{split}$
(110)

where ${\bm{p}}^{\prime}={\bm{p}}+{\bm{q}}$. These are the momentum-dependent
Eliashberg equations for the Hamiltonian (100), cf. Eq. (83) for the Holstein
model.

We look for an isotropic solution of Eq. (110), so that $\Phi_{n{\bm{p}}}$,
$\Sigma_{n{\bm{p}}}$ and $\chi_{n{\bm{p}}}$ are independent of the direction
of ${\bm{p}}$. Then, it is useful to convert the summation over ${\bm{q}}$ in
Eq. (110) into an integral over $q$ and $p^{\prime}$, where $q$ and
$p^{\prime}$ are the magnitudes of vectors $\bm{q}$ and $\bm{p}^{\prime}$.
Using ${p^{\prime}}^{2}=p^{2}+q^{2}+2pq\cos\theta$, we evaluate the Jacobian
for the change of variables $(q,\cos\theta)\to(q,p^{\prime})$ and find,

$\sum_{\bm{q}}=\frac{V}{4\pi^{2}}\int\\!\\!\\!\int
q^{2}dqd\cos\theta=\frac{V}{4\pi^{2}}\frac{1}{p}\int qdq\int
p^{\prime}dp^{\prime}.$ (111)

Since the right hand side of Eq. (110) is a product of a function of $q$ only
and a function of $p^{\prime}$ only, the integration over ${\bm{q}}$
factorizes into a product of an integral over $q$ and an integral over
$p^{\prime}$.

In the Migdal-Eliashberg theory, the fluctuations around the stationary point
are negligible only when the Fermi energy is much larger than any other
characteristic energy [1, 2]. This means that the integration is confined to
the vicinity of the Fermi surface, which we take to be spherical, and
$p\approx p^{\prime}\approx p_{F}$, where $p_{F}$ is the Fermi momentum. Since
${\bm{q}}$ connects ${\bm{p}}$ and ${\bm{p}}^{\prime}$, which are both on the
Fermi surface, its magnitude varies between 0 and $2p_{F}$. This allows as to
rewrite Eq. (111) as,

$\sum_{\bm{q}}=\frac{V}{4\pi^{2}}\frac{1}{p_{F}^{2}}\int_{0}^{2p_{F}}\\!\\!\\!\\!qdq\int{p^{\prime}}^{2}dp^{\prime}=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\\!\\!\\!\\!qdq\sum_{{\bm{p}}^{\prime}}.$
(112)

Using this in Eq. (110), we observe $\Phi_{n{\bm{p}}}$, $\Sigma_{n{\bm{p}}}$
and $\chi_{n{\bm{p}}}$ are independent of the momentum ${\bm{p}}$ and
performing the summation over ${\bm{p}}^{\prime}$ we arrive at the following
two equations:

$\begin{split}\Phi_{n}=\pi
T\sum_{m}\lambda(\omega_{n}-\omega_{m})\frac{\Phi_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}},\\\
\Sigma_{n}=\pi
T\sum_{m}\lambda(\omega_{n}-\omega_{m})\frac{\omega_{m}+\Sigma_{m}}{\sqrt{(\omega_{m}+\Sigma_{m})^{2}+|\Phi_{m}|^{2}}},\end{split}$
(113)

where,

$\lambda(\omega_{l})=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\frac{g_{q}^{2}qdq}{\omega_{l}^{2}+\omega_{q}^{2}},\quad
g_{q}^{2}=|\alpha_{q}|^{2}M^{-1}.$ (114)

This expression is for a spherical Fermi surface in $d=3$ dimensions, but it
is straightforward to extend it to any $d\geq 2$. Eqs. (113) are of the same
form as Eqs. (86), but with a more general kernel
$\lambda(\omega_{n}-\omega_{m})$. For dispersionless phonons,
$\omega_{q}=\Omega$ and $g_{q}=g$, this kernel too coincides with that in Eq.
(86). By definition the dimensionless coupling constant for dispersing phonons
is,

$\lambda=\lambda(\omega_{l}=0)=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\frac{g_{q}^{2}qdq}{\omega^{2}_{q}}.$
(115)

The effective action evaluated at the stationary point is [cf. Eq. (85)],

$\begin{split}{S}_{\mathrm{eff}}=T\nu_{0}V\sum_{nl}\left[\Phi_{n+l}^{*}\Lambda_{l}\Phi_{n}+\Sigma_{n+l}\Lambda_{l}\Sigma_{n}\right]-2\pi\nu_{0}V\sum_{n}\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}.\end{split}$
(116)

Here $\Lambda_{l}$ is the Fourier transform of $1/\lambda(\tau)$ at frequency
$\omega_{l}$. The corresponding free energy $f=TS_{\mathrm{eff}}/N$ per
lattice site is,

$\begin{split}f=\nu_{0}T^{2}\sum_{nl}\left[\Phi_{n+l}^{*}\lambda_{l}^{-1}\Phi_{n}+\Sigma_{n+l}\lambda_{l}^{-1}\Sigma_{n}\right]-2\pi\nu_{0}T\sum_{n}\sqrt{(\omega_{n}+\Sigma_{n})^{2}+|\Phi_{n}|^{2}}.\end{split}$
(117)

###  2 BCS theory as the weak coupling limit of the Migdal-Eliashberg theory

Confusion exists in the literature to this day regarding whether or not the
BCS theory is the weak coupling limit of the Migdal-Eliashberg theory [21,
22]. This confusion arises from attempts to compare cutoff-dependent
quantities which are illegitimate within the BCS theory. The BCS theory is
only valid in the limit $\Delta_{\mathrm{BCS}}/\Xi\to 0$, where
$\Delta_{\mathrm{BCS}}$ is the BCS ground state gap and $\Xi$ is the cutoff of
the order of the Debye energy for acoustic phonons and $\Omega_{0}$ for
Einstein phonons (Holstein model).

In other words, one should take the limit $\Xi\to\infty$ and the coupling
constant $\lambda\to 0$ while keeping $\Delta_{\mathrm{BCS}}\propto\Xi
e^{-1/\lambda}$ fixed. Only quantities that survive this limit are legitimate
within the BCS theory. In particular, it is meaningless to compare
$\Delta_{\mathrm{BCS}}/\Xi$ or $T_{c}/\Xi$, but, for example,
$\Delta_{\mathrm{BCS}}/T_{c}$, the condensation energy, the normalized jump in
the specific heat are meaningful and their BCS values should agree with those
in the weak coupling limit of the Migdal-Eliashberg theory, which they indeed
do [36, 75, 76, 8].

It is well-understood how the BCS theory emerges from the Migdal-Eliashberg
theory in the weak coupling limit. Here we reproduce known arguments within
our path integral framework. Our aim in doing so is to be able to compare
later the weak and strong coupling limits of the Migdal-Eliashberg theory as
well as Anderson and Eliashberg spins. We define the weak coupling limit as
the limit where the ratio of the energy gap to the characteristic phonon
frequency $\omega_{\mathrm{ch}}$ goes to zero. This is equivalent to
$\omega_{\mathrm{ch}}\to\infty$, cf. $\omega_{\mathrm{ch}}=\Omega\to 0$ in the
strong coupling limit of the Holstein model, see Eq. (5).

Since relevant energies and frequencies are of the order of the gap,
$\omega_{n}$ in the denominator is negligible in the effective action (101).
Then, the interaction is instantaneous and the effective action corresponds to
the Hamiltonian,

$H=\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}-\frac{1}{2N}\sum_{{\bm{p}}{\bm{p}}^{\prime}{\bm{q}}\sigma\sigma^{\prime}}\frac{\alpha_{q}^{2}M^{-1}}{\omega_{q}^{2}}c^{\dagger}_{{\bm{p}}+{\bm{q}}\sigma}c_{{\bm{p}}\sigma}c^{\dagger}_{{\bm{p}}^{\prime}-{\bm{q}}\sigma^{\prime}}c_{{\bm{p}}^{\prime}\sigma^{\prime}}.$
(118)

The hierarchy of scales
$\varepsilon_{F}\gg\omega_{\mathrm{ch}}\gg\Delta_{\mathrm{BCS}}$ implies that
the pairing interaction is confined to the vicinity of the Fermi surface and
that momenta of phonons mediating the interaction are of the order of the
Fermi momentum. Therefore, the range of the electron-electron interaction is
of order of the Fermi wavelength $\lambda_{F}$. The characteristic length
scale (coherence length) in a weakly coupled superconductor is of the order of
$v_{F}/\Delta_{\mathrm{BCS}}$, which is orders of magnitude larger than
$\lambda_{F}$. We see that the dimensionless range of interactions goes to
zero in the weak coupling limit. This means the interaction potential is a
delta function in the position space and a constant in the momentum space.
Since the pairing occurs at the Fermi surface in this limit, we determine the
constant by averaging the interaction over this surface,

$\mathrm{const}=\frac{1}{2p_{F}^{2}}\int_{0}^{2p_{F}}\frac{\alpha_{q}^{2}M^{-1}qdq}{\omega_{q}^{2}}=\frac{\lambda}{\nu_{0}},$
(119)

where we used Eq. (5). This is very similar to how we obtained Eq. (114) in
the Migdal-Eliashberg theory. The interaction in the position basis takes the
form
$c^{\dagger}_{\sigma}({\bm{r}})c_{\sigma}({\bm{r}})c^{\dagger}_{\sigma^{\prime}}({\bm{r}})c_{\sigma^{\prime}}({\bm{r}})$
For $\sigma=\sigma^{\prime}$ this is equal to
$c^{\dagger}_{\sigma}({\bm{r}})c_{\sigma}({\bm{r}})$, which we absorb into the
chemical potential term. We are left with $\sigma\neq\sigma^{\prime}$ terms,
i.e., with the BCS Hamiltonian,

$H=\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}-\lambda\delta\sum_{{\bm{p}}{\bm{p}}^{\prime}{\bm{q}}}c^{\dagger}_{{\bm{p}}+{\bm{q}}\uparrow}c^{\dagger}_{{\bm{p}}^{\prime}-{\bm{q}}\downarrow}c_{{\bm{p}}^{\prime}\downarrow}c_{{\bm{p}}\uparrow},$
(120)

where $\delta=(\nu_{0}N)^{-1}$ is the single-particle level spacing at the
Fermi energy.

Both BCS and Eliashberg theories are mean-field theories. In mean-field
approach, only Cooper pairs with zero total momentum determine the ground
state and low lying excitation spectrum. Keeping only such pairs, we end up
with the reduced BCS Hamiltonian [20],

$H=\sum_{{\bm{p}}\sigma}\xi_{\bm{p}}c^{\dagger}_{{\bm{p}}\sigma}c_{{\bm{p}}\sigma}-\lambda\delta\sum_{{\bm{p}}{\bm{p}}^{\prime}}c^{\dagger}_{-{\bm{p}}\uparrow}c^{\dagger}_{{\bm{p}}\downarrow}c_{{\bm{p}}^{\prime}\downarrow}c_{-{\bm{p}}^{\prime}\uparrow}.$
(121)

We expect the above sequence of steps leading to Eq. (121) to be exact when we
take both the weak coupling and thermodynamic limits. The neglect of frequency
and momentum dependencies of the electron-electron interaction is valid at
energies much smaller than $\omega_{\mathrm{ch}}$. It resulted in infinite
range interaction in momentum space between Cooper pairs in Eq. (121) [or,
equivalently, spins in Eq. (53)]. Therefore, this infinite range interaction
has to be cutoff at $\Xi\sim\omega_{\mathrm{ch}}$. Any quantity to which
energies of order $\omega_{\mathrm{ch}}$ contribute is beyond the BCS theory.

The weak coupling limit of the Migdal-Eliashberg theory – the BCS theory – is
universal in the sense that it is independent of the phonon dispersion and
momentum dependence of the electron-phonon coupling. There is only one low
energy scale in this limit – the BCS gap $\Delta_{\mathrm{BCS}}$. The gap
itself cannot be determined accurately from within the BCS theory. However,
the ratio of any other energy to the gap, e.g., $T_{c}/\Delta_{\mathrm{BCS}}$
is well-defined. The strong coupling limit of the Migdal-Eliashberg theory is
similarly universal [30].

## References

  * [1] A. B. Migdal, Interaction between Electrons and Lattice Vibrations in a Normal Metal, Sov. Phys. JETP 7, 996 (1958).
  * [2] G. M. Eliashberg, Interactions between Electrons and Lattice Vibrations in a Superconductor, Sov. Phys. JETP 11, 696 (1960).
  * [3] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Dover Publications, Inc., New York, 1975).
  * [4] J. Bauer, J. E. Han, and O. Gunnarsson, Quantitative reliability study of the Migdal-Eliashberg theory for strong electron-phonon coupling in superconductors, Phys. Rev. B 84, 184531 (2011).
  * [5] W. L. McMillan and J. M. Rowell, Tunneling and Strong-Coupling Superconductivity, in Superconductivity, edited by R.D. Parks (Marcel Dekker, Inc., New York, 1969) p. 561.
  * [6] P. B. Allen and R. C. Dynes, Transition temperature of strong-coupled superconductors reanalyzed, Phys. Rev. B 12, 905 (1975).
  * [7] P. B. Allen and B. Mitrovic, Theory of superconducting $T_{c}$, in Solid State Physics, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic, New York, 1982), Vol. 37, p. 1.
  * [8] J.P. Carbotte, Properties of boson-exchange superconductors, Rev. Mod. Phys. 62, 1027 (1990).
  * [9] J. Franck, Experimental studies of the isotope effect in high temperature superconductors, in Physical Properties of High Temperature Superconductors IV, edited by D.M. Ginsberg (World Scientific, Singapore, 1994) p. 189.
  * [10] H. J. Choi, M. L. Cohen, S. G. Louie, Anisotropic Eliashberg theory of MgB2: $T_{c}$, isotope effects, superconducting energy gaps, quasiparticles, and specific heat, Physica C 385, 66 (2003).
  * [11] P. Benedetti and R. Zeyher, Holstein model in infinite dimensions at half-filling, Phys. Rev. B 58, 14320 (1998).
  * [12] A. S. Alexandrov, Breakdown of the Migdal-Eliashberg theory in the strong-coupling adiabatic regime, EPL, 56, 92 (2001).
  * [13] D. Meyer, A. C. Hewson, and R. Bulla, Gap Formation and Soft Phonon Mode in the Holstein Model, Phys. Rev. Lett. 89, 196401(2002).
  * [14] M. Capone and S. Ciuchi, Polaron Crossover and Bipolaronic Metal-Insulator Transition in the Half-Filled Holstein Model, Phys. Rev. Lett. 91, 186405 (2003).
  * [15] I. Esterlis et. al., Breakdown of the Migdal-Eliashberg theory: A determinant quantum Monte Carlo study, Phys. Rev. B 97, 140501(R) (2018).
  * [16] I. Esterlis, S. A. Kivelson, and D. J. Scalapino, Pseudogap crossover in the electron-phonon system, Phys. Rev. B 99, 174516 (2019).
  * [17] This theorem says that corrections to the Migdal-Eliashberg theory are negligible when $\lambda_{0}\omega_{\mathrm{ch}}/\varepsilon_{F}\to 0$, where $\varepsilon_{F}$ is the Fermi energy, and when $\lambda_{0}$ is sufficiently small [1, 2], namely, “not too close too $\frac{1}{2}$” [3].
  * [18] A. V. Chubukov, A. Abanov, I. Esterlis, S. A. Kivelson, Eliashberg theory of phonon-mediated superconductivity – When it is valid and how it breaks down, Ann. Phys. 417, 168190 (2020).
  * [19] Y.-M. Wu, S.-S. Zhang, A. Abanov, A. V. Chubukov, Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal: V. The $\gamma$ model and its phase diagram. The case $\gamma=2$, Phys. Rev. B 103, 184508 (2021).
  * [20] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of Superconductivity, Phys. Rev. 108, 1175 (1957).
  * [21] F. Marsiglio, Eliashberg theory in the weak-coupling limit, Phys. Rev. B 98, 024523 (2018).
  * [22] S. Mirabi, R. Boyack, and F. Marsiglio Eliashberg theory in the weak-coupling limit: Results on the real frequency axis, Phys. Rev. B 101, 064506 (2020).
  * [23] E. A. Yuzbashyan and B. L. Altshuler, Breakdown of Migdal-Eliashberg Theory and Theory of Lattice-Fermionic Superfluidity, unpublished.
  * [24] E. A. Yuzbashyan et. al., Superconductivity near a quantum critical point in an extreme retardation limit (tentative title), unpublished.
  * [25] E. A. Yuzbashyan and B. L. Altshuler, Fluctuations in boson mediated superconductors and the meaning of the Migdal theorem (tentative title), unpublished.
  * [26] P. W. Anderson, Random-Phase Approximation in the Theory of Superconductivity, Phys. Rev. 112, 1900 (1958).
  * [27] G. Bergmann and D. Rainer, The sensitivity of the transition temperature to changes in $\alpha^{2}F(\omega)$, Z. Phys. 59, 263 (1973).
  * [28] F. Marsiglio and J. P. Carbotte, Gap function and density of states in the strong-coupling limit for an electron-boson system, Phys. Rev. B 43, 5355 (1991).
  * [29] A. E. Karakozov, E. G. Maksimov, and A. A. Mikhailovsky, The investigation of Eliashberg equations for superconductors with strong electron-phonon interaction, Solid State Commun. 79, 329 (1991).
  * [30] R. Combescot, Strong-coupling limit of Eliashberg theory, Phys. Rev. B 51, 11625 (1995).
  * [31] We often refer to $f$ as the free energy or free energy density. To be more precise, taken at its minimum $f$ is the grand potential per site as we work at fixed chemical potential.
  * [32] The arbitrariness of $\hat{h}$ does not play an essential role in this paper, but will be important later [23].
  * [33] J. Bisgard, Mountain Passes and Saddle Points, SIAM Rev., 57, 275 (2015).
  * [34] G. M. Eliashberg, The low temperature specific heat of metals, Sov. Phys. JETP 16, 780 (1963) [Russian original: Zh. Eksperim. i Teor. Fiz. 43, 1005 (1962).]
  * [35] Y. Wada, Critical Field and Specific Heat of Strong Coupling Superconductors, Phys. Rev, 135, A1481 (1964).
  * [36] J. Bardeen and M. Stephen, Free-Energy Difference Between Normal and Superconducting States, Phys. Rev. 136, A1485 (1964).
  * [37] A. Abanov, B. Altshuler, A. V. Chubukov, E. Yuzbashyan, On a Strong Coupling Theory of Superconductivity, Bulletin of the American Physical Society, abstract N13.014, (2004).
  * [38] E. Tsoncheva and A. V. Chubukov, Erratum: Condensation Energy in Eliashberg theory: From weak to strong coupling, Phys. Rev. B 72, 149902 (2005).
  * [39] R. Haslinger and A. V. Chubukov, Condensation energy in strongly coupled superconductors, Phys. Rev. B 68, 214508 (2003).
  * [40] This constant is in fact zero by the zeta function regularization prescription [41]. Indeed, the summand here depends only on $p=n-m$. Summation over $p$ leaves us with a sum over $l=n+m$ of the form $\sum_{l=-\infty}^{\infty}1=1+2\zeta(0)=0$. See also Ref. 23.
  * [41] H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets (World Scientific Pub Co Inc., 2009).
  * [42] As usual, $T=0$ on Matsubara axis is understood as the limit $T\to 0$ and we assume $\Delta(\omega)$ [$T\to 0$ limit of $\Delta(\omega_{n})$] is at least piecewise continuous. It is clear then that all these properties commute with $T\to 0$ limit.
  * [43] This sets the arbitrary overall phase $e^{i\phi}$ factor of $\Delta(\omega_{n})$ to one.
  * [44] One must subtract a constant from $H_{s}$, replacing $\omega_{n}\cos\theta_{n}$ with $\omega_{n}\cos\theta_{n}-|\omega_{n}|$ and $\cos(\theta_{n}-\theta_{m})$ with $\cos(\theta_{n}-\theta_{m})-1$, to make the summations converge. The first replacements subtracts the free energy of noninteracting Fermi gas, the second is discussed below Eq. (19). This does not affect the stationary points in any way.
  * [45] M. Gaudin, The Bethe Wavefunction (Cambridge University Press, 2014).
  * [46] E K Sklyanin, Boundary conditions for integrable quantum systems, J. Phys. A, 21, 2375 (1988).
  * [47] E. K. Sklyanin, Separation of variables in the classical integrable $SL(3)$ magnetic chain, Commun.Math. Phys. 150, 181 (1992).
  * [48] E. A. Yuzbashyan, B. L. Altshuler, V. B. Kuznetsov, V. Z. Enolskii, Solution for the dynamics of the BCS and central spin problems, J. Phys. A 38, 7831, (2005).
  * [49] G. Ortiz, R. Somma, J. Dukelsky, S. Rombouts, Exactly-solvable models derived from a generalized Gaudin algebra, Nucl. Phys. B 707 421 (2005).
  * [50] E. Stouten, P. W. Claeys, J.-S. Caux, and V. Gritsev, Integrability and duality in spin chains, Phys. Rev. B 99, 075111 (2019).
  * [51] E. A. Yuzbashyan and B. S. Shastry, Quantum integrability in systems with finite number of levels, J. Stat. Phys. 150, 704 (2013).
  * [52] S.-S. Zhang, Y.-M. Wu, A. Abanov, A. V. Chubukov, Interplay between superconductivity and non-Fermi liquid at a quantum-critical point in a metal. VI. The $\gamma$ model and its phase diagram at $2<\gamma<3$, arXiv:2107.14340 (2021).
  * [53] These papers use the same expression for the free energy as we do [see Eq. (39) in Ref. 19], except they subtract the normal state free energy and accordingly call it the condensation energy. They also take $\lambda_{nm}=g^{2}/(\omega_{n}-\omega_{m})^{\gamma}$ and consider other values of $\gamma$ in addition to $\gamma=2$ (our case). However, all results of this section hold equally well for all values of $\gamma>0$.
  * [54] Indeed, $|\bm{S}_{n}\cdot\bm{S}_{m}-1|\leq 2$ and Eq. (114) shows that $\lambda_{nm}\leq g^{2}/(\omega_{n}-\omega_{m})^{2}$. It follows that the interaction (free) energy of any spin with other spins cannot exceed $\nu_{0}g^{2}\zeta(2)$ in magnitude. Because of this the interaction cannot compensate for a deviation of $S_{n}^{z}$ from $\mp 1$ at $\omega_{n}=\pm\infty$. Therefore, $\Delta(\omega_{n})$ must vanish at $\omega_{n}=\pm\infty$.
  * [55] The high energy scale is the Fermi energy $\varepsilon_{F}$. The hierarchy of scales in the strong coupling limit is $\varepsilon_{F}\gg g\gg\Omega$ and must be maintained when $\lambda\to\infty$. The analogous hierarchy in the weak coupling limit is $\varepsilon_{F}\gg\Omega\gg\Delta_{\mathrm{BCS}}$.
  * [56] We reversed the signs of the $z$ and $y$ components of spins with respect to Anderson’s original definition.
  * [57] R. W. Richardson, Pairing in the limit of a large number of particles, J. Math. Phys. 18, 1802 (1977).
  * [58] J. Roman, G. Sierra, and J. Dukelsky, Large N limit of the exactly solvable BCS model: Analytics versus numerics, Nucl. Phys. B 634, 483 (2002).
  * [59] E. A. Yuzbashyan, A. A. Baytin, and B. L. Altshuler, Finite-size corrections for the pairing Hamiltonian, Phys. Rev. B 71, 094505 (2005).
  * [60] More precisely, they correspond to those eigenstates that do not contain unpaired fermions, i.e., to the BCS ground state and excitations of the BCS condensate.
  * [61] E. A. Yuzbashyan, M. Dzero, V. Gurarie, and M. S. Foster, Quantum quench phase diagrams of an s-wave BCS-BEC condensate, Phys. Rev. A 91, 033628 (2015).
  * [62] E. A. Yuzbashyan, Normal and anomalous solitons in the theory of dynamical Cooper pairing, Phys. Rev. B 78, 184507 (2008).
  * [63] E. A. Yuzbashyan, O. Tsyplyatyev and B. L. Altshuler, Relaxation and Persistent Oscillations of the Order Parameter in Fermionic Condensates, Phys. Rev. Lett. 96, 097005 (2006).
  * [64] R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Collective Rabi Oscillations and Solitons in a Time-Dependent BCS Pairing Problem, Phys. Rev. Lett. 93, 160401 (2004).
  * [65] M. Dzero, E. A. Yuzbashyan and B. L. Altshuler, Cooper pair turbulence in atomic Fermi gases, EPL 85, 20004 (2009).
  * [66] R. A. Barankov and L. S. Levitov, Synchronization in the BCS Pairing Dynamics as a Critical Phenomenon, Phys. Rev. Lett. 96, 230403 (2006).
  * [67] E. A. Yuzbashyan and M. Dzero, Dynamical Vanishing of the Order Parameter in a Fermionic Condensate, Phys. Rev. Lett. 96, 230404 (2006).
  * [68] N. Tsuji and H. Aoki, Theory of Anderson pseudospin resonance with Higgs mode in superconductors, Phys. Rev. B 92, 064508 (2015).
  * [69] G.M. Eliashberg, Inelastic electron collisions and nonequilibrium stationary states in superconductors, Sov. Phys. JETP 34, 668 (1972).
  * [70] A. F. Volkov and S. M. Kogan, Collisionless relaxation of the energy gap in superconductors, Sov. Phys. JETP, 38, 1018 (1974).
  * [71] D. Hauck, M. J. Klug, I. Esterlis, J. Schmalian, Eliashberg equations for an electron-phonon version of the Sachdev-Ye-Kitaev model: Pair breaking in non-Fermi liquid superconductors, Ann. Phys. 417 168120, (2020).
  * [72] M. Protter, R. Boyack, and F. Marsiglio, Functional-integral approach to Gaussian fluctuations in Eliashberg theory, Phys. Rev. B 104, 014513 (2021).
  * [73] E.G. Brovman and Yu. Kagan, The phonon spectrum of metals, Sov. Phys. JETP 25, 365 (1967).
  * [74] B. T. Gelikman, Adiabatic perturbation theory for metals and the problem of lattice stability, Sov. Phys. Usp. 18, 190 (1975).
  * [75] B. Mitrovíc, H.G. Zarate, and J.P. Carbotte, The ratio $2\Delta_{0}/(k_{B}T_{c})$ within Eliashberg theory, Phys. Rev. B 29 184 (1984).
  * [76] F. Marsiglio and J.P. Carbotte, Strong-coupling corrections to Bardeen-Cooper-Schrieffer ratios, Phys. Rev. B 33 6141 (1986).