# Néel proximity effect at antiferromagnet/superconductor interfaces

G. A. Bobkov  Moscow Institute of Physics and Technology, Dolgoprudny, 141700
Moscow region, Russia  I. V. Bobkova  Institute of Solid State Physics,
Chernogolovka, 142432 Moscow region, Russia  Moscow Institute of Physics and
Technology, Dolgoprudny, 141700 Moscow region, Russia  National Research
University Higher School of Economics, 101000 Moscow, Russia  A. M. Bobkov
Institute of Solid State Physics, Chernogolovka, 142432 Moscow region, Russia
Moscow Institute of Physics and Technology, Dolgoprudny, 141700 Moscow region,
Russia  Akashdeep Kamra  Condensed Matter Physics Center (IFIMAC) and
Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma
de Madrid, E-28049 Madrid, Spain

###### Abstract

Spin-splitting induced in a conventional superconductor weakens
superconductivity by destroying spin-singlet and creating spin-triplet Cooper
pairs. We demonstrate theoretically that such an effect is also caused by an
adjacent compensated antiferromagnet, which yields no net spin-splitting. We
find that the antiferromagnet produces Néel triplet Cooper pairs, whose
pairing amplitude oscillates rapidly in space similar to the antiferromagnet’s
spin. The emergence of these unconventional Cooper pairs reduces the singlet
pairs’ amplitude, thereby lowering the superconducting critical temperature.
We develop a quasiclassical Green’s functions description of the system
employing a two-sublattice framework. It successfully captures the rapid
oscillations in the Cooper pairs’ amplitude at the lattice spacing scale as
well as their smooth variation on the larger coherence length scale. Employing
the theoretical framework thus developed, we investigate this Néel proximity
effect in a superconductor/antiferromagnet bilayer as a function of
interfacial exchange, disorder, and chemical potential, finding rich physics.
Our findings also offer insights into experiments which have found a larger
than expected suppression of superconductivity by an adjacent antiferromagnet.

##  I Introduction

Conventional superconductors are formed by spin-singlet Cooper pairs [1, 2].
Exposing them to a net spin-splitting, such as via an applied magnetic field
or due to interfacial exchange interaction with a ferromagnet (F), causes
spin-singlet pairs to be converted into their spin-triplet counterparts [3, 4,
5, 6, 7]. This weakens the conventional superconducting state and lowers the
critical temperature [8, 9, 10]. On the other hand, since the net
magnetization in an antiferromagnet (AF) vanishes, an adjacent superconductor
(S) interfaced to the former via a compensated interface is expected to
experience no net spin-splitting or reduction in critical temperature [11,
12]. Nevertheless, unconventional Andreev reflection and bound states at such
S/AF interfaces have been predicted [13, 14]. The rich Josephson physics in
S/AF/S hybrids has also been investigated, theoretically [15, 16, 17, 18] and
experimentally [19, 20, 21].

Several experiments have found that AFs lower the critical temperature of an S
layer [19, 22, 23, 24], despite the no net spin-splitting argument above. In
some cases, the effect has been comparable to or even larger than that induced
by a ferromagnet layer [23]. To understand this, several potential
consequences of the AF layer have been considered. First, an AF doubles the
spatial period of the lattice due to its antiparallel spins on the two
sublattices. This can open a bandgap in the adjacent conductor, which may
reduce the normal-state density of states in S and thus, superconductivity
[25, 26]. Second, it has been shown than an uncompensated interface, which
seems to be common in experiments [27, 11, 28, 29], to an AF insulator does
induce a net spin-splitting [11], just like a ferromagnet. Furthermore, the
interfacial disorder was found to cause spin-flip scattering and reduce
superconductivity [11, 30]. While these offer potential mechanisms for
affecting the S, they do not account for the phenomena that underlie the
previously considered unconventional Andreev reflection [13, 14]. Also, the
question of how can AFs affect adjacent S more than ferromagnets remains
unanswered. Furthermore, a recent Bogoliubov-de Gennes numerical analysis of a
hybrid comprising a compensated AF interfaced with an S suggested the
interface to be spin-active [31]. Hence, a key piece of the puzzle in
understanding S/AF bilayers appears to be missing.

In this article, we undertake a detailed theoretical investigation of an S/AF
bilayer with a compensated interface. For simplicity, we assume the AF to be
an insulator. We first analyze this system numerically solving the Bogoliubov-
de Gennes equation on a two-dimensional lattice. In this analysis, we find
that the Néel order of the AF induces spin-triplet correlations in the S.
Their amplitude flips sign from one lattice site to the next, just like the
Néel spin order in the AF. Thus, we call these Néel triplet Cooper pairs. The
AF with its compensated interface has induced rapidly oscillating triplet
correlations, that escape the conventional quasiclassical description of the S
since it cannot resolve variations on such short length scales [32, 33].

In order to adequately account for such rapid oscillations, we develop a two-
sublattice quasiclassical Green’s functions description of the S. Employing
it, we obtain analytic results for the Néel pairing amplitude and the critical
temperature of the ensuing superconducting state, finding results consistent
with our Bogoliubov-de Gennes numerics. Furthermore, the developed framework
allows us to address the role of chemical potential in the normal state,
disorder, and the strength of interfacial exchange semi-analytically thereby
providing valuable insights.

We find that the Néel triplets are formed due to interband pairing within the
two-sublattice model. Thus, they are formed when the chemical potential is
such that there are two bands around the Fermi surface within the energy
$\Delta_{0}$, the superconducting gap of S without the adjacent AF. The
formation of such Néel spin-triplets comes at the expense of destroying the
spin-singlet correlation which reduces the S critical temperature. On the
other hand, disorder destroys this interband pairing and diminishes the
proximity effect of the AF on the S. Thus, our investigated mechanism weakens
superconductivity strongly for clean systems, yielding an opposite trend with
respect to the role of disorder-mediated spin-flip scattering. Thus, the role
of disorder in proximity effect with an AF can be positive or negative on the
superconductivity. These competing effects further allow maxima to be feasible
in experiments [23, 22, 24].

##  II Bogoliubov-de Gennes analysis

We begin by numerically examining the system of interest - an
antiferromagnetic insulator interfaced via a compensated interface to a thin
superconductor [Fig. 1(a)]. To this end, we set up a tight-binding Bogoliubov-
de Gennes Hamiltonian for the two materials [34] and assume on-site local
s-wave correlations in the superconductor, assumed conventional. For
simplicity, we consider a two-dimensional system ($12\times 100$ spatial
cluster was considered) and employ periodic boundary conditions along the
interfacial direction. The latter emulates an infinite-length interface. The
alternating Néel spin order in the AF induces a correspondingly oscillating
spin-splitting at the interfacial lattice sites in the S via exchange
interaction between the AF spins and the S electrons. Describing the
Hamiltonian and methodology details in the appendix, we numerically
diagonalize the Hamiltonian and evaluate the superconducting state self-
consistently.

Figure 1: (a) Sketch of the antiferromagnetic insulator interfaced via a
compensated interface to the thin superconductor, considered in the framework
of a tight-binding Bogoliubov-de Gennes Hamiltonian in Sec. II. The system
represents a two-dimensional 12 × 100 spatial cluster. Blue points are S
sites, red and green points correspond to AF sites with opposite directions of
the on-site magnetization $\pm\bm{m}$. (b) Spatial variation of the triplet
correlations amplitude $F_{\bm{i}}^{t}$ in the investigated AF/S bilayer. Each
colored square codes the value of $F_{\bm{i}}^{t}$ at a given site. Only small
vertical part of the bilayer is shown, which is marked with a black rectangle
in panel (a). An alternating sign of the correlations in S commensurate with
the Néel order in the AF can be seen along the interfacial direction. The
triplet amplitude is normalized to the hopping amplitude, see appendix.

The spin-triplet correlations amplitude $F_{\bm{i}}^{t}$ at each lattice site
with the radius-vector $\bm{i}$ is evaluated by summing the anomalous Green’s
function over the positive Matsubara frequencies, as detailed in the appendix.
Figure 1(b) plots the spatially resolved spin-triplet pairing amplitude in the
investigated bilayer. A clear imprinting of the AF Néel order is seen on the
triplet pairing amplitude in the direction parallel to the interface. This
perfectly commensurate variation of the triplet correlations is disturbed in
the direction perpendicular to the interface by Friedel oscillations [34] and
a general lack of out-of-plane momentum conservation due to the interface. A
finite $F_{\bm{i}}^{t}$ in the first few lattice sites of the AF is due to a
small leakage of the electron wavefunctions into the insulating AF.
Furthermore, although not explicitly shown, the superconducting critical
temperature is found to be reduced substantially by the AF.

Our numerical analysis clearly demonstrates a large proximity effect of the AF
on the S, despite a compensated interface resulting in no net spin-splitting
on spatial averaging over the S coherence length scale [11]. It also shows
that the interesting physics is taking place on the lattice constant length
scale, which is beyond the resolution of the conventional quasiclassical
Green’s functions description of the superconductor [32]. Johnsen and
coworkers also recently found the compensated interface between AF and S to be
spin-active using a similar Bogoliubov-de Gennes numerics [31]. However, they
considered trilayers involving an additional ferromagnet, which prevents a
clear understanding of the phenomena taking place at the AF/S interface.

##  III Quasiclassical Green’s functions description of a two-sublattice
system

The description of superconductors in terms of quasiclassical Green’s
functions has proven immensely powerful in understanding all kinds of hybrids
involving magnets and superconductors [33, 32, 6, 5, 7]. This framework is
made tractable by averaging over rapid spatial variations on the Fermi
wavelength length scale, which is comparable to or larger than the lattice
spacing. Such a procedure adequately captures the properties of the
superconducting state while ignoring some small details the underlying normal-
metal state.

Motivated and guided by our numerical results based on solving the Bogoliubov-
de Gennes equation (Sec. II), we wish to develop a quasiclassical Green’s
function description capable of capturing these effects semi-analytically. In
doing so, we notice that the rapid oscillations in the pairing amplitude on
the lattice constant scale are merely an expression of the two-sublattice
nature of the AF/S system. Thus, these oscillations can be adequately captured
by working with a two-sublattice framework and the spatial variations on an
individual sublattice are expected to remain slow, as compared to the lattice
spacing. This method is directly analogous to capturing the spin order in an
antiferromagnet [35], that shows rapid oscillation between the two sublattices
but slow spatial variations within a single sublattice.

###  III.1 Two-sublattice Eilenberger equation

Figure 2: Schematic depiction of the setup under consideration. A Néel ordered
antiferromagnet (AF) is interfaced via a compensated interface to a
superconductor (S). The interface is in the $y-z$ plane and the first
superconducting layer is at $i_{x}=0$. The lattice in both parts of the
structure is divided into two sublattices A and B. The red arrows depict
localized spins in the AF. The basis vectors of the original (prior to the
introduction of A and B sublattices) lattice in the superconductor are
$\bm{a}_{x}$, $\bm{a}_{y}$, $\bm{a}_{z}$.

The unit cell with two sublattices - A and B - is introduced as shown in Fig.
2. In the framework of this two-sublattice approach the unit cells as a whole
are marked by radius-vector $\bm{i}$. Then the staggered magnetism is
described by $\bm{m}_{\bm{i},A(B)}=+(-)\bm{m}_{\bm{i}}$, where
$\bm{m}_{\bm{i}}$ is the local magnetic moment at site A of the unit cell with
the radius vector $\bm{i}$ in the AF. This allows us to consider
$\bm{m}_{\bm{i}}$ as a slow function of the spatial coordinate. We consider
the homogeneously ordered Néel state of the AF here, such that
$\bm{m}_{\bm{i}}$ does not depend on the position $\bm{i}$. The superconductor
S is described by the Hamiltonian:

$\displaystyle\hat{H}=$
$\displaystyle-t\sum\limits_{\langle\bm{i}\bm{j}\nu\bar{\nu}\rangle,\sigma}\hat{\psi}_{\bm{i}\sigma}^{\nu\dagger}\hat{\psi}_{\bm{j}\sigma}^{\bar{\nu}}+\sum\limits_{\bm{i},\nu}(\Delta_{\bm{i}}^{\nu}\hat{\psi}_{\bm{i}\uparrow}^{\nu\dagger}\hat{\psi}_{\bm{i}\downarrow}^{\nu\dagger}+H.c.)-\mu\sum\limits_{\bm{i}\nu,\sigma}\hat{n}_{\bm{i}\sigma}^{\nu}$
$\displaystyle+\sum\limits_{\bm{i}\nu,\alpha\beta}\hat{\psi}_{\bm{i}\alpha}^{\nu\dagger}(\bm{h}_{\bm{i}}^{\nu}\bm{\sigma})_{\alpha\beta}\hat{\psi}_{\bm{i}\beta}^{\nu}+\sum\limits_{\bm{i}\nu,\sigma}V_{i}^{\nu}\hat{n}_{\bm{i}\sigma}^{\nu},$
(1)

where $\nu=A,B$ is the sublattice index, $\bar{\nu}=A(B)$ if $\nu=B(A)$ means
that the corresponding quantity belongs to the opposite sublattice,
$\langle\bm{i}\bm{j}\nu\bar{\nu}\rangle$ means summation over the nearest
neighbors,
$\hat{\psi}_{\bm{i}\sigma}^{\nu\dagger}(\hat{\psi}_{\bm{i}\sigma}^{\nu})$ is
the creation (annihilation) operator for an electron with spin $\sigma$ at the
sublattice $\nu$ of the unit cell $\bm{i}$. $t$ parameterizes the hopping
between adjacent sites, $\Delta_{\bm{i}}^{\nu}$ accounts for on-site s-wave
pairing, $\mu$ is the electron chemical potential, and $V_{i}^{\nu}$ is the
local on-site potential that is later employed to capture the effect of
impurities and disorder.
$\hat{n}_{\bm{i}\sigma}^{\nu}=\hat{\psi}_{\bm{i}\sigma}^{\nu\dagger}\hat{\psi}_{\bm{i}\sigma}^{\nu}$
is the particle number operator at the site belonging to sublattice $\nu$ in
unit cell $\bm{i}$. The influence of the antiferromagnetic insulator on the
superconductor is described by the exchange field
$\bm{h}_{\bm{i}}^{A(B)}=\bm{h}(-\bm{h})\delta_{i_{x},0}$, where the
$\delta$-symbol means that the exchange field is only nonzero at the AF/S
interface sites corresponding to $i_{x}=0$ in the superconductor [11]. We
assume that the interface is fully compensated, that is the interface exchange
field is staggered with zero average value.

The Matsubara Green’s function in the two-sublattice formalism is $8\times 8$
matrix in the direct product of spin, particle-hole and sublattice spaces.
Introducing the two-sublattice Nambu spinor
$\check{\psi}_{\bm{i}}=(\hat{\psi}_{{\bm{i}},\uparrow}^{A},\hat{\psi}_{\bm{i},\downarrow}^{A},\hat{\psi}_{\bm{i},\uparrow}^{B},\hat{\psi}_{\bm{i},\downarrow}^{B},\hat{\psi}_{\bm{i},\uparrow}^{A\dagger},\hat{\psi}_{\bm{i},\downarrow}^{A\dagger},\hat{\psi}_{\bm{i},\uparrow}^{B\dagger},\hat{\psi}_{\bm{i},\downarrow}^{B\dagger})^{T}$
we define the Green’s function as follows:

$\displaystyle\check{G}_{\bm{i}\bm{j}}(\tau_{1},\tau_{2})=-\langle
T_{\tau}\check{\psi}_{\bm{i}}(\tau_{1})\check{\psi}_{\bm{j}}^{\dagger}(\tau_{2})\rangle,$
(2)

where $\langle T_{\tau}...\rangle$ means imaginary time-ordered thermal
averaging. At first, we consider the clean case corresponding to
$V_{\bm{i}}^{\nu}=0$. For the system described by Hamiltonian (1) the
Heisenberg equation of motion for spinor $\check{\psi}_{\bm{i}}$ takes the
form:

$\displaystyle\frac{d\check{\psi}_{\bm{i}}}{d\tau}=[\hat{H},\check{\psi}_{\bm{i}}]=\tau_{z}\Bigl{(}t\hat{j}+\mu-\check{\Delta}_{\bm{i}}\sigma_{y}-\bm{h}\check{\bm{\sigma}}\delta_{i_{x},0}\rho_{z}\Bigr{)}\check{\psi}_{\bm{i}},$
(3)

where $\sigma_{k}$, $\tau_{k}$ and $\rho_{k}$ ($k=x,y,z$) are Pauli matrices
in spin, particle-hole and sublattice spaces, respectively.
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})^{T}$ is the vector of Pauli
matrices in spin space,
$\check{\bm{\sigma}}=\bm{\sigma}(1+\tau_{z})/2+\bm{\sigma}^{*}(1-\tau_{z})/2$
is the quasiparticle spin operator and
$\check{\Delta}_{\bm{i}}=\Delta_{\bm{i}}\tau_{+}+\Delta_{\bm{i}}^{*}\tau_{-}$
with $\tau_{\pm}=(\tau_{x}\pm i\tau_{y})/2$. Here we assume that
$\Delta_{\bm{i}}^{A}=\Delta_{\bm{i}}^{B}=\Delta_{\bm{i}}$, that is the order
parameter values are equal for $A$ and $B$ sites of a unit cell. For this
reason $\check{\Delta}_{\bm{i}}$ is proportional to the unit matrix in the
sublattice space. It follows from physical considerations since the change
$A\to B$ is equivalent to $\bm{h}\to-\bm{h}$ and the order parameter should
not depend on odd powers of the exchange field for the problem under
consideration. Also this physical assumption is confirmed by subsequent self-
consistent calculations. In general, for a case of non-equivalent sublattices
matrix $\check{\Delta}_{\bm{i}}$ is also diagonal in the sublattice space, but
the diagonal components can be not equal
$\check{\Delta}_{\bm{i}}^{A}\neq\check{\Delta}_{\bm{i}}^{B}$.

$\displaystyle\hat{j}\check{\psi}_{\bm{i}}=\rho_{+}\sum\limits_{\bm{a}}\check{\psi}_{\bm{i}+\bm{a}-\bm{a}_{z}}+\rho_{-}\sum\limits_{\bm{a}}\check{\psi}_{\bm{i}+\bm{a}+\bm{a}_{z}}.$
(4)

Here $\bm{a}\in\\{\pm\bm{a}_{x},\pm\bm{a}_{y},\pm\bm{a}_{z}\\}$, see Fig. 2
for definition of these vectors. $\rho_{\pm}=(\rho_{x}\pm i\rho_{y})/2$. The
Green’s function Eq. (2) obeys equation:

$\displaystyle\frac{d\check{G}_{\bm{i}\bm{j}}}{d\tau_{1}}=-\delta(\tau_{1}-\tau_{2})\delta_{\bm{i}\bm{j}}-\langle
T_{\tau}\frac{d\check{\psi}_{\bm{i}}(\tau_{1})}{d\tau_{1}}\check{\psi}_{\bm{j}}^{\dagger}(\tau_{2})\rangle.$
(5)

Substituting Eq. (3) into Eq. (5) and expressing the Green’s function in terms
of the Matsubara frequencies $\omega_{m}=\pi T(2m+1)$ we obtain:

$\displaystyle
G_{\bm{i},l}^{-1}\check{G}_{\bm{i}\bm{j}}(\omega_{m})=\delta_{\bm{i}\bm{j}},$
(6) $\displaystyle
G_{\bm{i},l}^{-1}=\tau_{z}\left(t\hat{j}+\mu-\check{\Delta}_{\bm{i}}\i\sigma_{y}-\bm{h}\check{\bm{\sigma}}\delta_{i_{x},0}\rho_{z}\right)+i\omega_{m}.$
(7)

Analogously one obtains the equation, where the operator $G^{-1}$ acts on the
Green’s function from the right:

$\displaystyle\check{G}_{\bm{i}\bm{j}}(\omega_{m})G_{\bm{j},r}^{-1}=\delta_{\bm{i}\bm{j}},$
(8) $\displaystyle
G_{\bm{j},r}^{-1}=\left(t\hat{j}+\mu+\check{\Delta}_{\bm{i}}\i\sigma_{y}-\bm{h}\check{\bm{\sigma}}\delta_{i_{x},0}\rho_{z}\right)\tau_{z}+i\omega_{m}.$
(9)

where the operator $\hat{j}$ acts on the Green’s function from the right as
follows

$\check{G}_{\bm{i}\bm{j}}\hat{j}=\sum\limits_{\bm{a}}\check{G}_{\bm{i},\bm{j}+\bm{a}-\bm{a}_{z}}\rho_{-}+\sum\limits_{\bm{a}}\check{G}_{\bm{i},\bm{j}+\bm{a}+\bm{a}_{z}}\rho_{+}.$
(10)

Further we consider the Green’s function in the mixed representation:

$\displaystyle\check{G}(\bm{R},\bm{p})=F(\check{G}_{\bm{i}\bm{j}})=\int
d^{3}re^{-i\bm{p}(\bm{i}-\bm{j})}\check{G}_{\bm{i}\bm{j}},$ (11)

where $\bm{R}=(\bm{i}+\bm{j})/2$ and the integration is over $\bm{i}-\bm{j}$.
The term $\hat{j}\check{G}_{ij}$ in the mixed representation takes the form:

$\displaystyle
F(\hat{j}\check{G}_{ij})=te^{ip_{z}a_{z}}\rho_{+}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{G}(\bm{R}+\frac{\bm{a}}{2}-\frac{\bm{a}_{z}}{2})+te^{-ip_{z}a_{z}}\rho_{-}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{G}(\bm{R}+\frac{\bm{a}}{2}+\frac{\bm{a}_{z}}{2}),~{}~{}~{}~{}~{}~{}$
(12)

$a_{i}$ is the absolute value of $\bm{a}_{i}$ for $i=x,y,z$. Now we define the
following transformed Green’s function:

$\displaystyle\check{\tilde{G}}(\bm{R},\bm{p})=\left(\begin{array}[]{cc}1&0\\\
0&-i\sigma_{y}\end{array}\right)_{\tau}\rho_{x}e^{\frac{-\displaystyle
ip_{z}a_{z}\rho_{z}}{\displaystyle
2}}\check{G}(\bm{R},\bm{p})e^{\frac{\displaystyle
ip_{z}a_{z}\rho_{z}}{\displaystyle 2}}\left(\begin{array}[]{cc}1&0\\\
0&-i\sigma_{y}\end{array}\right)_{\tau},$ (17)

where subscript $\tau$ means that the explicit matrix structure corresponds to
the particle-hole space. Then taking into account Eq. (12) from Eq. (6) we
obtain:

$\displaystyle\left[i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\right]\rho_{x}\check{\tilde{G}}(\bm{R},\bm{p})+t\frac{\rho_{0}+\rho_{z}}{2}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}-\frac{\bm{a}_{z}}{2},\bm{p})$
$\displaystyle+t\frac{\rho_{0}-\rho_{z}}{2}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}+\frac{\bm{a}_{z}}{2},\bm{p})=1,~{}~{}$
(18)

where $\rho_{0}$ is the unit matrix in the sublattice space. The dependence of
the Green’s function on $\bm{R}$ is assumed to be slow, therefore from now we
can consider $\bm{R}$ as a continuous variable. We also generalized our
consideration to the case of an arbitrary spatial profile of the exchange
field $\bm{h}(\bm{R})$.

Analogously, from Eq. (8) it follows:

$\displaystyle\check{\tilde{G}}(\bm{R},\bm{p})\left[i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\right]\rho_{x}+t\sum\limits_{\bm{a}}e^{i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}+\frac{\bm{a}_{z}}{2},\bm{p})\frac{\rho_{0}+\rho_{z}}{2}$
$\displaystyle+t\sum\limits_{\bm{a}}e^{i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}-\frac{\bm{a}_{z}}{2},\bm{p})\frac{\rho_{0}-\rho_{z}}{2}=1.~{}~{}$
(19)

Subtracting Eqs. (18) and (19) we obtain:

$\displaystyle\left[\left(i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\right)\rho_{x},\check{\tilde{G}}(\bm{R},\bm{p})\right]+\check{A}=0,$
(20)

where,

$\displaystyle\check{A}=$ $\displaystyle
t\frac{\displaystyle\rho_{0}+\rho_{z}}{\displaystyle
2}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}-\frac{\bm{a}_{z}}{2},\bm{p})+t\frac{\displaystyle\rho_{0}-\rho_{z}}{\displaystyle
2}\sum\limits_{\bm{a}}e^{-i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}+\frac{\bm{a}_{z}}{2},\bm{p})$
$\displaystyle-t\sum\limits_{\bm{a}}e^{i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}+\frac{\bm{a}_{z}}{2},\bm{p})\frac{\displaystyle\rho_{0}+\rho_{z}}{\displaystyle
2}-t\sum\limits_{\bm{a}}e^{i\bm{p}\bm{a}}\check{\tilde{G}}(\bm{R}+\frac{\bm{a}}{2}-\frac{\bm{a}_{z}}{2},\bm{p})\frac{\displaystyle\rho_{0}-\rho_{z}}{\displaystyle
2}.$ (21)

Now we introduce quasiclassical $\xi$-integrated Green’s function:

$\displaystyle\check{g}(\bm{R},\bm{p}_{F})=-\frac{1}{i\pi}\int\check{\tilde{G}}(\bm{R},\bm{p})d\xi,$
(22)

where $\xi(\bm{p})=-2t(\cos p_{x}a_{x}+\cos p_{y}a_{y}+\cos p_{z}a_{z})$ is
the normal state electron dispersion counted from the Fermi energy. In
general, in the framework of the two-sublattice formalism the Brillouin Zone
(BZ) is reduced along the $p_{z}$ direction, that is
$p_{z}\in[-\pi/2a_{z},\pi/2a_{z}]$. This results in appearance of two normal
state quasiparticle dispersion branches $\xi_{\pm}=\mp 2t(\cos p_{x}a_{x}+\cos
p_{y}a_{y}+\cos p_{z}a_{z})-\mu$ in the reduced BZ. Only one of them has a
Fermi surface in the first BZ. Let it be $\xi_{+}$, then the Fermi momentum
$\bm{p}_{F}$ is determined by the equation $\xi_{+}(\bm{p}_{F})=0$. The
quasiclassical Green’s function (22) depends only on the direction of the
momentum on the Fermi surface, as usual. Please note that in the framework of
the developed formalism $\mu$ is assumed to be small with respect to $t$ and,
therefore, can be neglected in $\xi_{+}(\bm{p})$, which is reduced to
$\xi(\bm{p})$ in this case. Moreover, if $\mu$ is of the order of $t$, the
quasiclassical formalism developed here fails to work. Technically it is
because $\mu\sim t$ still enters Eq. (20) in contrast to the well-known one-
sublattice formalism, and the quasiclassical formalism does not allow such a
high-energy term. Physically it means suppression of all the correlations
between electrons belonging to different normal state dispersion branches and
Umklapp scattering processes[36, 37, 35], see Sec. III.4 for further details.
In the framework of the quasiclassical approximation
$\check{\tilde{G}}(\bm{R}+\bm{a}/2\pm\bm{a}_{z}/2,\bm{p})\approx\check{\tilde{G}}(\bm{R},\bm{p})+(\bm{a}\pm\bm{a}_{z})\bm{\nabla}\check{\tilde{G}}(\bm{R},\bm{p})/2$
and we obtain

$\displaystyle-\frac{1}{i\pi}\int\check{A}d\xi=it\sum\limits_{\bm{a}}\sin\bm{p}_{F}\bm{a}(\bm{\nabla}\check{g}(\bm{R},\bm{p}_{F})\bm{a})=i\bm{v}_{F}\bm{\nabla}\check{g}(\bm{R},\bm{p}_{F}),$
(23)

where
$\bm{v}_{F}=d\xi/d\bm{p}|_{\bm{p}=\bm{p}_{F}}=2t(\bm{a}_{x}\sin[p_{x}a_{x}]+\bm{a}_{y}\sin[p_{y}a_{y}]+\bm{a}_{z}\sin[p_{z}a_{z}])$.
After the $\xi$-integration the resulting Eilenberger equation for the
quasiclassical Green’s function takes the form:

$\displaystyle\left[\left(i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\right)\rho_{x},\check{g}(\bm{R},\bm{p}_{F})\right]+i\bm{v}_{F}\bm{\nabla}\check{g}(\bm{R},\bm{p}_{F})=0.$
(24)

Eq. (24) should be supplemented by the normalization condition

$\displaystyle\check{g}^{2}(\bm{R},\bm{p}_{F})=1,$ (25)

which is typical for the quasiclassical theory. In order to verify its
validity in our case one can multiply Eq. (24) by
$\check{g}(\bm{R},\bm{p}_{F})$ from the left, then from the right and add the
resulting equations. This procedure leads to the conclusion that
$\check{g}^{2}(\bm{R},\bm{p}_{F})$ obeys the same Eq. (24) and, therefore,
$\check{g}^{2}(\bm{R},\bm{p}_{F})=1$ is a solution of this equation for
arbitrary spatial profiles of $\check{\Delta}(\bm{R})$ and $\bm{h}(\bm{R})$.
In particular, one can assume that $\check{\Delta}(\bm{R})$ and
$\bm{h}(\bm{R})$ evolve smoothly to zero values, that is the system transforms
to a normal metal. Given that Eq. (25) holds in the normal metal one can
conclude that $\check{g}(\bm{R},\bm{p}_{F})$ obeying Eqs. (24) and (25) is a
true solution of the problem under consideration. Below we directly show that
Eq. (25) is valid in the limit of the normal metal $\check{\Delta}(\bm{R})\to
0$ and $\bm{h}(\bm{R})\to 0$.

Equations (24) and (25) are the desired two-sublattice Eilenberger equation
together with the normalization condition on the quasiclassical Green’s
function matrix. These constitute the main result of this subsection. The
sublattice degree of freedom adds the $\rho_{x,y,z}$ Pauli matrices to the
framework.

###  III.2 Inclusion of disorder

So far, we have considered a homogeneous and clean superconductor. A key
strength of the quasiclassical framework is that it can account for several
nonidealities, such as impurity scattering or disorder. These are especially
important in modeling and understanding experiments which employ
superconductors with a varying degree of disorder.

Considering impurity potential in the superconductor in the framework of the
Born approximation [32], it can be obtained that Eq. (20) is modified as
follows:

$\displaystyle\Bigl{[}\Bigl{(}i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\Bigr{)}\rho_{x}-\check{\Sigma}_{\mathrm{imp}}(\bm{R}),\check{\tilde{G}}(\bm{R},\bm{p})\Bigr{]}+\check{A}=0,~{}~{}~{}~{}~{}~{}~{}$
(26)

where the impurity self-energy $\check{\Sigma}_{\mathrm{imp}}$ takes the form:

$\displaystyle\check{\Sigma}_{\mathrm{imp}}=\frac{1}{\pi
N_{F}\tau}\int\frac{d^{3}p}{(2\pi)^{3}}[\rho_{+}\check{\tilde{G}}(\bm{R},\bm{p})\rho_{+}+\rho_{-}\check{\tilde{G}}(\bm{R},\bm{p})\rho_{-}]$
$\displaystyle=\frac{1}{\pi N_{F}\tau}\int d\xi\langle
N_{F}(\bm{p}_{F})[\rho_{+}\check{\tilde{G}}(\bm{R},\bm{p})\rho_{+}+\rho_{-}\check{\tilde{G}}(\bm{R},\bm{p})\rho_{-}]\rangle_{\bm{p}_{F}},~{}~{}~{}~{}~{}~{}$
(27)

where $\langle...\rangle_{\bm{p}_{F}}$ means Fermi-surface averaging,
$N_{F}(\bm{p}_{F})$ is the momentum direction-dependent density of states
(DOS) at the Fermi surface. $N_{F}$ is the momentum - averaged DOS at the
Fermi-surface. $\tau$ is the quasiparticle mean free time, which is connected
to the on-site impurity potential $V_{i}$ as $\langle V_{i}V_{j}\rangle=(1/\pi
N_{F}\tau)\delta_{ij}$. In this case one can express the impurity self-energy
Eq. (27) via the quasiclassical Green’s function:

$\displaystyle\check{\Sigma}_{\mathrm{imp}}(\bm{R})=\frac{i}{\tau}[\rho_{+}\check{\bar{g}}(\bm{R})\rho_{+}+\rho_{-}\check{\bar{g}}(\bm{R})\rho_{-}],$
(28)

where
$\check{\bar{g}}(\bm{R})=\langle\frac{N_{F}(\bm{p}_{F})}{N_{F}}g(\bm{R},\bm{p}_{F})\rangle_{\bm{p}_{F}}$.
Then Eq. (24) in the presence of impurities takes the form:

$\displaystyle\Bigl{[}\Bigl{(}i\omega_{m}\tau_{z}+\mu+\tau_{z}\check{\Delta}(\bm{R})-\bm{h}(\bm{R})\bm{\sigma}\tau_{z}\rho_{z}\Bigr{)}\rho_{x}-\check{\Sigma}_{\mathrm{imp}}(\bm{R}),\check{g}(\bm{R},\bm{p}_{F})\Bigr{]}+i\bm{v}_{F}\bm{\nabla}\check{g}(\bm{R},\bm{p}_{F})=0.~{}~{}~{}~{}$
(29)

###  III.3 Green’s function of a homogeneous normal metal at an arbitrary
disorder

Now let’s find the quasiclassical Green’s function in the spatially
homogeneous normal metal in the presence of spatially homogeneous, but
staggered, exchange field $\bm{h}$ as produced by a compensated
antiferromagnet interface. Then the Green’s function
$\check{\tilde{G}}(\bm{R},\bm{p})$ does not depend on $\bm{R}$ and Eq. (18) is
reduced to:

$\displaystyle\Bigl{[}\Bigl{(}i\omega_{m}\tau_{z}+\mu-\bm{h}\bm{\sigma}\tau_{z}\rho_{z}\Bigr{)}\rho_{x}-\xi(\bm{p})-\check{\Sigma}_{\mathrm{imp}}\Bigr{]}\check{\tilde{G}}(\bm{p})=1,~{}~{}~{}~{}$
(30)

where we have also added the impurity self-energy. In the clean limit
$\tau\to\infty$, i.e., if we disregard impurity scattering, Eq. (30) can be
easily solved resulting in the following answer for the $\xi$-integrated
Green’s function Eq. (22):

$\displaystyle\check{g}=\tau_{z}\frac{\bm{h}\bm{\sigma}\rho_{y}+(\omega_{m}-i\tau_{z}\mu)\rho_{x}}{\sqrt{h^{2}+(\omega_{m}-i\mu\tau_{z})^{2}}}.$
(31)

It is seen that the normalization condition Eq. (25) immediately follows from
Eq. (31). It is important that the Néel structure of the exchange field
radically changes the well-known answer for the quasiclassical Green’s
function in the normal metal $\check{g}_{N}=1$, which is valid for the case of
ferromagnet-type exchange field $\bm{h}$. This is primarily because the
staggered spin-splitting due to the antiferromagnet opens a gap in the normal-
state local density of states (LDOS) of the superconductor:

$\displaystyle N(\varepsilon)=\frac{1}{4}{\rm
Tr}_{8}\Bigl{[}\frac{1+\tau_{z}}{2}\rho_{x}{\rm
Re}[\check{g}(i\omega_{m}\to\varepsilon+i\delta)]\Bigr{]}$ (32)

at $\varepsilon=-\mu$. The LDOS at $\mu=0$ and in the absence of impurities is
shown in Fig. 3 by the red line. The LDOS at finite $\mu$ is obtained by
shifting the corresponding plots along the energy axis $N_{\mu\neq
0}(\varepsilon)=N_{\mu=0}(\varepsilon+\mu)$.

Figure 3: Normal state LDOS for different values of impurity scattering time
$\tau$, which is measured in units of the inverse bulk superconducting
critical temperature $T_{c0}^{-1}$. $\varepsilon$ is measured in units of
$T_{c0}$. We consider $\mu=0$ and $h=0.3T_{c0}$ in this figure.

In the presence of impurities solution of Eq. (30) can not be written
explicitly, but the quasiclassical Green’s function still has the off-diagonal
structure in the sublattice space, that is

$\displaystyle\check{g}=\left(\begin{array}[]{cc}0&\check{g}^{AB}\\\
\check{g}^{BA}&0\end{array}\right),$ (35)

where

$\displaystyle\check{g}^{AB}=\sqrt{[\check{\beta}+\frac{i}{\tau}\check{g}^{BA}]^{-1}[\check{\alpha}+\frac{i}{\tau}\check{g}^{AB}]},$
$\displaystyle\check{g}^{BA}=\sqrt{[\check{\alpha}+\frac{i}{\tau}\check{g}^{AB}]^{-1}[\check{\beta}+\frac{i}{\tau}\check{g}^{BA}]}$
(36)

with

$\displaystyle\check{\alpha}(\check{\beta})=\mp\bm{h}\bm{\sigma}\tau_{z}+i(\omega_{m}\tau_{z}-i\mu).$
(37)

From Eqs. (35) and (36) it follows that in the presence of impurities the
normalization condition Eq. (25) also holds.

Irrespective of the impurity strength the general answer for normal
quasiclassical Green’s function takes the form:

$\displaystyle\check{g}=\left(\begin{array}[]{cc}\hat{g}_{N}&0\\\
0&\hat{\tilde{g}}_{N}\end{array}\right)_{\tau},$ (40)

where the electron and hole components of the Green’s function take the
following structure in the sublattice space:

$\displaystyle\hat{g}_{N}=A\rho_{x}+B\rho_{y}\bm{n}_{h}\bm{\sigma},$
$\displaystyle\hat{\tilde{g}}_{N}=\tilde{A}\rho_{x}+\tilde{B}\rho_{y}\bm{n}_{h}\bm{\sigma},$
(41)

$\bm{n}_{h}=\bm{h}/h$. The difference between two sublattices is contained in
the second terms in Eq. (41). It is seen that it is directly determined by the
presence of the staggered exchange field. The normal state LDOS at finite
impurity concentration is presented in Fig. 3. It is seen that the gap in the
LDOS is graduallly suppressed with impurity strength.

###  III.4 Relating single-sublattice and two-sublattice pictures

Figure 4: Electron dispersion of the normal-state S in the reduced Brillouin
zone (BZ) $pa\in[-\pi/2,\pi/2]$ considering a 1D system with two sites in the
unit cell $\xi_{\pm}(p)=\mp 2t\cos pa-\mu$. The reciprocal lattice vector due
to the periodicity enforced by the AF is $Q_{1D}=\pi/a$. The spectrum branches
are doubled in the BZ due to the reduction of the BZ volume. The blue line
indicates ordinary pairing between $(p_{0},\xi_{1}=0)$(1) and
$(-p_{0},\xi_{2}=0)$(2) electrons corresponding to the zero total pair
momentum. The green line indicates Néel pairing between $p_{0}$(1) and
$-p_{0}+Q_{1D}$ (3) corresponding to the total pair momentum $Q_{1D}$. From
the point of view of the first BZ it is an interband pairing between electrons
(1) and (3’). Taking into account that $p_{0}$ is defined from the condition
$-2t\cos p_{0}a-\mu=0$ one immediately obtains that
$\xi_{1}-\xi_{3^{\prime}}=2\mu$. That is, the energy difference between (1)
and (3’) electrons grows with $\mu$ thus reducing the efficiency of pairing.
The antiferromagnetic gap opening as discussed in Fig. 3 has been disregarded
in the present simplified figure. $\Delta_{0}$ is the zero-temperature gap of
the bulk S.

We now discuss a physical picture for the Néel Cooper pairs. The essential
physics is captured conveniently within a one-dimensional (1D) model and we
consider that here for clarity [38]. In this model, a 1D AFM is interfaced to
a 1D superconductor running along the AFM. Therefore, the electron wavevector
bears only one component which is along the interface.

Now, if we disregard the AFM, the normal-state electronic dispersion of S can
be depicted as in Fig. 4 with a BZ $pa\in[-\pi,\pi]$, where $a$ is the lattice
constant. Within this single-sublattice dispersion for the normal-state of S,
the AF Néel order imprinted on the S via a staggered spin-splitting field
causes scattering between electronic states that differ by the wavenumber
$Q_{1D}=\pi/a$, which is the reciprocal state unit vector for the AF. This has
sometimes been termed Umklapp scattering in the literature [36, 37, 35].
However, we should keep in mind that the momentum change is the reciprocal
unit vector of the AF, and not the S. Thus, the AF converts conventional spin-
singlet pairing between $+p$ and $-p$ electronic states into Néel spin-triplet
pairing between, for example, $+p$ and $-p+Q_{1D}$ (see Fig. 4). Such a
pairing oscillates rapidly in space similar to the Néel order with the
wavenumber $Q_{1D}$.

In the discussion above, we have disregarded the gap opening (see Fig. 3)
caused by the imposition of staggered spin-splitting on the S. This gap
opening seems natural when we recognize that the Néel ordered AF reduces the
periodicity of the adjacent S by imposing a staggered spin-splitting on it and
reducing the BZ to $pa\in[-\pi/2,\pi/2]$. Thus, within this adequately
rigorous two-sublattice picture, the gap opening is natural and happens at the
BZ boundary [38, 26]. We continue to ignore it in the ongoing discussion and
in Fig. 4, returning to the gap in understanding further physical phenomena
below. Furthermore, within this two-sublattice picture, we now have two bands
in the electronic dispersion. What appeared as pairing between $+p_{0}$ and
$-p_{0}+Q_{1D}$ states in the single-sublattice picture is actually pairing
between the $+p_{0}$ state from one band with the $-p_{0}$ state from the
other band, as depicted in Fig. 4.

We also see that the Néel or interband pairing is taking place between two
states which generally have different energies. They becomes degenerate only
under the special condition when the chemical potential $\mu$ is zero so that
we have perfect half filling in the system, and when we disregard the gap
opening due to the staggered spin-splitting. As Cooper pairing takes place
between electrons within an energy $\Delta_{0}$ (the superconducting gap of
the bulk S) from the Fermi surface, the strength of interband pairing depends
on the LDOS in the normal-state. This understanding will be employed to
qualitatively understand several of our results below.

##  IV Critical temperature of the superconductor/antiferromagnet bilayer

We now employ the quasiclassical formalism developed in the previous section
to study the proximity effect and the superconducting critical temperature in
the AF/S bilayer. The S layer is assumed to be thin with respect to the
coherence length i.e., $d_{S}\ll\xi_{S}$. In this case the Green’s function
$\check{g}$ is approximately constant in the superconductor except for the
trajectories nearly parallel to the AF/S interface [39, 40] and Eq. (24) can
be integrated over the S width $d_{S}$.

We work at temperatures close to the critical temperature and linearize the
Eilenberger equation with respect to $\Delta/T_{c}$. In this approximation the
Green’s function takes the form:

$\check{g}=\left(\begin{array}[]{cc}\hat{g}_{N}&\hat{f}\\\
\hat{\tilde{f}}&\hat{\tilde{g}}_{N}\end{array}\right)_{\tau},$ (42)

where all the components are $4\times 4$ matrices in the direct product of
spin and sublattice spaces. The diagonal components $\hat{g}_{N}$ and
$\hat{\tilde{g}}_{N}$ are to be calculated in the normal state of the
superconductor, what is done in the previous section. The anomalous components
$\hat{f}$ and $\hat{\tilde{f}}$ are of the first order with respect to
$\Delta/T_{c}$ and contain singlet and Néel triplet correlations. The
resulting equation for the anomalous Green’s function $\hat{f}$ takes the
form:

$\displaystyle\left\\{i\omega_{m}\rho_{x}-\bm{h}_{eff}\bm{\sigma}\rho_{z}\rho_{x},\hat{f}\right\\}+\left[\mu\rho_{x},\hat{f}\right]+\Delta\left(\rho_{x}\hat{\tilde{g}}_{N}-\hat{g}_{N}\rho_{x}\right)+$
$\displaystyle\frac{i}{\tau}\left([\rho_{+}\hat{g}_{N}\rho_{+}+\rho_{-}\hat{g}_{N}\rho_{-}]\hat{f}-\hat{f}[\rho_{+}\hat{\tilde{g}}_{N}\rho_{+}+\rho_{-}\hat{\tilde{g}}_{N}\rho_{-}]\right)=0,~{}~{}~{}~{}$
(43)

where $h_{eff}=ha_{x}/d_{S}$ is the absolute value of the effective staggered
exchange field, induced in the superconductor due to proximity with the AF
insulator. Now let $\bm{h}_{eff}=h_{eff}\bm{e}_{z}$. Then the solution of Eq.
(43) takes the form:

$\displaystyle\hat{f}=f_{s}\rho_{x}+f_{t}\rho_{y}\sigma_{z},$ (44)

where

$\displaystyle f_{s}=$
$\displaystyle\Delta\frac{-i\left(A^{2}-\tilde{A}^{2}\right)+\tau\left[(\tilde{A}-A)\mu+ih_{eff}(B+\tilde{B})\right]}{2\left[(B+\tilde{B})h_{eff}+\omega_{m}(A+\tilde{A}-i\mu\tau)\right]},~{}~{}~{}~{}$
(45) $\displaystyle f_{t}=$
$\displaystyle\Delta\frac{i(B+\tilde{B})(A-\tilde{A}+\omega_{m}\tau)}{2\left[(B+\tilde{B})h_{eff}+\omega_{m}(A+\tilde{A}-i\mu\tau)\right]}.~{}~{}~{}~{}$
(46)

In the clean limit, Eqs. (45) and (46) are reduced to

$\displaystyle f_{s}=$
$\displaystyle\frac{\Delta}{2i\omega_{m}}\Bigl{(}\frac{\omega_{m}-i\mu}{\sqrt{h_{eff}^{2}+(\omega_{m}-i\mu)^{2}}}+\frac{\omega_{m}+i\mu}{\sqrt{h_{eff}^{2}+(\omega_{m}+i\mu)^{2}}}\Bigr{)}$
$\displaystyle+\frac{\Delta
h_{eff}^{2}}{2\mu\omega_{m}}\Bigl{(}\frac{1}{\sqrt{h_{eff}^{2}+(\omega_{m}-i\mu)^{2}}}-\frac{1}{\sqrt{h_{eff}^{2}+(\omega_{m}+i\mu)^{2}}}\Bigr{)},$
(47) $\displaystyle f_{t}=$
$\displaystyle\frac{h_{eff}\Delta}{2\mu}\Bigl{(}\frac{1}{\sqrt{h_{eff}^{2}+(\omega_{m}-i\mu)^{2}}}-\frac{1}{\sqrt{h_{eff}^{2}+(\omega_{m}+i\mu)^{2}}}\Bigr{)}.$
(48)

From Eq. (44), it is seen that the triplet part of the anomalous Green’s
function is opposite at two sublattices (due to the $\rho_{y}$ term). Thus,
the quasiclassical approach adequately captures and demonstrates the Néel
triplets previously seen in the numerical solution (Fig. 1).

Figure 5: Anomalous Green’s function of the Néel triplet correlations summed
over positive Matsubara frequencies
$F_{A}^{t}=\sum\limits_{\omega_{m}>0}f_{t}$ as a function of $h_{eff}$ for
different values of the mean free time $\tau$. Red line represents the same
quantity for an S/F interface with a ferromagnetic insulator producing the
same value of the effective exchange field (but homogeneous, not staggered) in
the superconductor. The S/F interface is not sensitive to impurities, for this
reason only one line is shown for the ferromagnetic case. Each line ends at
the critical value of $h_{eff}$ corresponding to the full suppression of
superconductivity. We consider $\mu=0$ here.

Eq. (48) demonstrates that the Néel triplets are odd-frequency correlations
[5, 41]. However, in contrast to the usual odd-frequency triplets in S/F
heterostructures [5, 7], they are suppressed by impurities as demonstrated in
Fig. 5, which presents the dependence of the Néel triplets on $h_{eff}$ at
different impurity strengths. The suppression of the Néel triplets by
impurities is because they are constituted by interband pairing and scattering
from impurities mixes the states in the two bands. Furthermore, it is seen
that in the absence of impurities the absolute value of the Neel triplet
correlations is very close to the value of triplets at the S/F interface with
the same parameters. The difference between them becomes essential at higher
value of $h_{eff}$, where the suppression of superconductivity at the S/AF
interface due to the antiferromagnetic gap also plays an important role.

Figure 6: Critical temperature of the AF/S bilayer as a function of $h_{eff}$
for different values of mean free time $\tau$. $\mu=0.01T_{c0}$. Orange line
represents $T_{c}(h_{eff})$ for an S/F interface with a ferromagnetic
insulator producing the same value of the effective exchange field (but
homogeneous, not staggered) in the superconductor. The dashed parts of the
curves correspond to the unstable solutions, see text for further details.

The critical temperature of the AF/S bilayer is calculated from the self-
consistency equation

$\displaystyle 1=i\pi\lambda
T_{c}\sum\limits_{\omega_{m}}\frac{f_{s}}{\Delta},$ (49)

where $\lambda$ is the dimensionless coupling constant. The dependence of the
critical temperature of the AF/S bilayer on the effective exchange field is
presented in Fig. 6 for different values of the impurity strength. We can
observe two important physical phenomena in this figure.

First, the critical temperature is indeed suppressed by the staggered
exchange. In the clean case the efficiency of suppression by the staggered
field is of the same order, and even higher, as the suppression by the
ferromagnet with the same absolute value of the exchange field. The stronger
suppression of the superconductivity by the staggered exchange as compared to
the uniform ferromagnetic exchange field is explained by the presence of the
antiferromagnetic gap at the Fermi surface [see Fig. 3], which prevents
electronic states inside this gap from superconducting pairing.

Second, the critical temperature for a given $h_{eff}$ grows with impurity
concentration. That is, the impurities reduce the efficiency of suppression.
It is in sharp contrast with the behavior of a F/S bilayer, where the degree
of suppression is not sensitive to the impurity concentration. This highly
unusual behavior results from two facts, which act in parallel. First of all,
as it is discussed above, the odd-frequency Néel triplets are suppressed by
impurities. Second, the antiferromagnetic gap at the Fermi surface is also
suppressed by impurities. It makes the corresponding electronic states
available for superconducting pairing.

The dashed parts of the curves in Fig. 6 represent the unstable branches of
$T_{c}(h_{eff})$. It really means that starting from some value of $h_{eff}$
the self-consistency equation for the superconducting order parameter (beyond
the linearized limit) has two different solutions corresponding to two
different critical temperatures. One of them (with higher critical
temperature) represents the stable solution and the other (with lower critical
temperature) corresponds to the maximum of the system free energy, that is
absolutely unstable. The situation is similar to the physics of S/F bilayers
[42, 8, 9, 7].

Figure 7: Anomalous Green’s function of the Neel triplet correlations as a
function of $h_{eff}$ for different values of $\mu$. Each line ends at the
critical value of $h_{eff}$ corresponding to the full suppression of
superconductivity. We consider the clean limit $\tau^{-1}=0$ here.

Néel triplets are also suppressed by finite values of $\mu$ for a given
$h_{eff}$. This is demonstrated in Fig. 7. The physical reason for the
suppression is understood from the discussion in Sec. III.4 and the
corresponding Fig. 4. Electrons can effectively pair if their energies are
within the bulk superconducting gap $\Delta_{0}$ from the Fermi surface. Néel-
type or interband pairing involves electrons with different energies with the
energy difference $2\mu$ between them (see Fig. 4). When this energy
difference exceeds $\Delta_{0}$, the Néel-type pairing is not effective. This
simple picture is valid for small values of $h_{eff}$ when the
antiferromagnetic gap in the electron dispersion can be ignored. As it is seen
from Fig. 7, at $h_{eff}\gtrsim\mu$ the suppression of Néel triplets by $\mu$
is greatly weakened. The suppression of the Neel triplets by nonzero $\mu$
naturally results in the reduced sensitivity of the critical temperature to
the exchange field at larger $\mu$, as has been demonstrated in Fig. 8.

Figure 8: Critical temperature of the AF/S bilayer as a function of $h_{eff}$
for different values of $\mu$. Unstable branches are not shown. The results
correspond to the clean limit $\tau^{-1}=0$.

##  V Discussion

A theoretical accounting of the bandgap created in the normal-state electron
dispersion by the staggered spin-splitting induced by an adjacent Néel ordered
AF has previously been considered [26], without disorder. Furthermore,
emergence of a finite spin-splitting in the S due to the commonly found
uncompensated interface and spin-flip scattering offer other known mechanisms
of proximity effect in an AF/S bilayer [11]. Our analysis in this article has
revealed a previously unexplored mechanism - converting spin-singlet into Néel
spin-triplet Cooper pairs - of suppressing superconductivity in the S by an
adjacent compensated AF. Furthermore, employing the two-sublattice
quasiclassical Green’s functions theory developed, our analysis reveals the
effect of disorder on all these mechanisms. This seemingly complete picture
allows us to envisage the various competing effects taking place in a
realistic AF/S bilayer [19, 22, 24, 23].

Let us understand the qualitative dependencies of these mechanisms. The role
of finite spin-splitting in S due to an uncompensated interface with an AF
together with spin-flip scattering due to interfacial disorder is nearly
identical to that of proximity effect due to an adjacent ferromagnet [30, 11].
For the case of a compensated interface that we have considered here, both the
bandgap opening mechanism and the induction of Néel triplets are most active
when the Fermi wavevector in the normal state of the S (without AF) is close
to the reciprocal space unit vector of the AF. In that case, the bandgap
opening tends to make the normal-state of S into an insulator thereby
diminishing density of states and weakening superconductivity. At the same
time, under those conditions, the two bands (Fig. 4) are energetically close
resulting in conversion of spin-singlet into Néel (or interband) spin-triplet
Cooper pairs. Notably, disorder suppresses both these effects. Furthermore, if
the S normal state does not have the required Fermi surface, both of these
mechanisms become inactive and one recovers the Werthammer treatment of the
proximity effect in the AF/S bilayer [12], where the AF, if metallic, merely
acts as a normal metal.

Hence, to some extent, the experimental observation of a finite proximity
effect and reduction of superconductivity in AF/S bilayers can be understood
as due to the uncompensated interface that causes spin-splitting and spin-flip
scattering. However, a proximity effect that turns out to be stronger than
that due to a ferromagnet and that shows a maximum at specific AF thickness
cannot be explained in this manner [19, 22, 23]. These observations hint at
the role of bandgap opening and/or the Néel triplets mechanism investigated
here. A decrease in the average disorder with increasing thickness of the AF
should allow for a stronger suppression of the superconductivity as per these
two mechanisms. With further increase in the AF thickness, the S Cooper pairs
do not leak into the metallic AF far enough and the proximity effect gets
saturated. Furthermore, with reduced disorder, the weakening of
superconductivity due to spin-flip scattering is decreased that may result in
a recovery of the S critical temperature. Nevertheless, a fully detailed
analysis of such a system should adequately consider a metallic AF, in
contrast with our consideration of proximity effect due to an insulating AF.

##  VI Summary

We have demonstrated theoretically the emergence of unconventional spin-
triplet Cooper pairs in a superconductor exposed via a compensated interface
to an antiferromagnetic insulator. The pairing amplitude for these alternates
sign on the lattice length scale, similar to the Néel spin order in the
antiferromagnet. Thus, we call these Néel triplets and show that these are
formed from pairing between electronic states from two different bands i.e.,
interband pairing. Formulating a two-sublattice quasiclassical Green’s
functions description of the bilayer, we investigate the effect of this
interband pairing on the superconducting critical temperature. The ensuing
powerful formulation provides a semi-analytic description and a clear physical
understanding of the key phenomena as a function of disorder and chemical
potential. Such Néel or interband triplets may already be playing an important
role in experiments on proximity effect in superconductor/antiferromagnet
hybrids.

## Acknowledgments

We thank Lina Johnsen for valuable discussions. We acknowledge financial
support from the Spanish Ministry for Science and Innovation – AEI Grant
CEX2018-000805-M (through the “Maria de Maeztu” Programme for Units of
Excellence in R&D) and from the Russian Science Foundation via the RSF project
No.22-22-00522.

##  Appendix A Bogoliubov-de Gennes analysis: formulation and solution

In this appendix, we formulate and solve the Bogoliubov-de Gennes equation
[34] which has been employed for the numerical analysis presented in Sec. II.
The system of our interest is described via the following Hamiltonian:

$\displaystyle
H=-t\sum\limits_{\langle\bm{i}\bm{j},\rangle\sigma}\hat{c}_{\bm{i}\sigma}^{\dagger}\hat{c}_{\bm{j}\sigma}+\sum\limits_{\bm{i}}(\Delta_{\bm{i}}\hat{c}_{\bm{i}\uparrow}^{\dagger}\hat{c}_{\bm{i}\downarrow}^{\dagger}+H.c.)-\mu\sum\limits_{\bm{i}\sigma}\hat{n}_{\bm{i}\sigma}+\sum\limits_{\bm{i},\alpha\beta}\hat{c}_{\bm{i}\alpha}^{\dagger}(\bm{m}_{\bm{i}}\bm{\sigma})_{\alpha\beta}\hat{c}_{\bm{i}\beta},$
(50)

where $\bm{i}=(i_{x},i_{y})^{T}$ is the radius vector of the site and greek
letters correspond to the spin indices. $\langle\bm{i}\bm{j}\rangle$ means
summation over the nearest neighbors. $\Delta_{\bm{i}}$ and $\bm{m}_{\bm{i}}$
denote the on-site superconducting and magnetic order parameters at site
$\bm{i}$, respectively.
$\hat{c}_{\bm{i}\sigma}^{\dagger}(\hat{c}_{\bm{i}\sigma})$ creates
(annihilates) an electron of spin $\sigma=\uparrow,\downarrow$ on the site
$\bm{i}$, $t$ denotes the nearest-neighbor hopping integral, $\mu$ is the
filling factor.
$\hat{n}_{\bm{i}\sigma}=\hat{c}_{\bm{i}\sigma}^{\dagger}\hat{c}_{\bm{i}\sigma}$
is the particle number operator at site $\bm{i}$. We also define the vector of
Pauli matrices in spin space
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})^{T}$. We assume that the
antiferromagnet is of G-type. Then the magnetic order parameter can be taken
in the form $\bm{m}_{\bm{i}}=(-1)^{i_{x}+i_{y}}\bm{m}$ inside the
antifferomagnet. $x$\- and $y$-axes are taken normal to the AF/S interface and
parallel to it, respectively. It is assumed that the antiferromagnetism is due
to the localized electrons and the amplitude of the on-site magnetization is
not influenced by the adjacent metal. Therefore, we do not calculate $\bm{m}$
self-consistently and it is assumed to be constant inside the AF region. We
diagonalize the Hamiltonian (50) by the Bogolubov transformation:

$\displaystyle\hat{c}_{\bm{i}\sigma}=\sum\limits_{n}u^{\bm{i}}_{n\sigma}\hat{b}_{n}+v^{\bm{i}*}_{n\sigma}\hat{b}_{n}^{\dagger},$
(51)

where $\hat{b}_{n}^{\dagger}(\hat{b}_{n})$ are the creation (annihilation)
operators of Bogolubov quasiparticles. Then the resulting Bogolubov-de Gennes
equations take the form:

$\displaystyle-\mu
u^{\bm{i}}_{n,\sigma}-t\sum\limits_{\bm{j}\in\langle\bm{i}\rangle}u_{n,\sigma}^{\bm{j}}+\sigma\Delta_{\bm{i}}v^{\bm{i}}_{n,-\sigma}+(\bm{m}_{\bm{i}}\bm{\sigma})_{\sigma\alpha}u_{n,\alpha}^{\bm{i}}$
$\displaystyle=\varepsilon_{n}u_{n,\sigma}^{\bm{i}}$ $\displaystyle-\mu
v^{\bm{i}}_{n,\sigma}-t\sum\limits_{\bm{j}\in\langle\bm{i}\rangle}v_{n,\sigma}^{\bm{j}}+\sigma\Delta_{\bm{i}}^{*}u^{\bm{i}}_{n,-\sigma}+(\bm{m}_{\bm{i}}\bm{\sigma}*)_{\sigma\alpha}v_{n,\alpha}^{\bm{i}}$
$\displaystyle=-\varepsilon_{n}v_{n,\sigma}^{\bm{i}},$ (52)

where $\bm{j}\in\langle\bm{i}\rangle$ means summation over nearest neighbors
of site $\bm{i}$ and $\varepsilon_{n}$ are the eigen-state energies of the
Bogolubov quasiparticles.

The superconducting order parameter in the S layer is calculated self-
consistently:

$\displaystyle\Delta_{\bm{i}}=g\langle\hat{c}_{\bm{i}\downarrow}\hat{c}_{\bm{i}\uparrow}\rangle=g\sum\limits_{n}(u_{n,\downarrow}^{\bm{i}}v_{n,\uparrow}^{\bm{i}*}(1-f_{n})+u_{n,\uparrow}^{\bm{i}}v_{n,\downarrow}^{\bm{i}*}f_{n}),$
(53)

where $g$ is the coupling constant. The quasiparticle distribution function is
assumed to be the equilibrium Fermi distribution $f_{n}=\langle
b_{n}^{\dagger}b_{n}\rangle=1/(1+e^{\varepsilon_{n}/T})$.

Further we investigate the structure of superconducting correlations at the
AF/S interface with infinite superconducting layer. The anomalous Green’s
function in Matsubara representation can be calculated as
$F_{\bm{i},\alpha\beta}=-\langle\hat{c}_{\bm{i}\alpha}(\tau)\hat{c}_{\bm{i}\beta}(0)\rangle$,
where $\tau$ is the imaginary time. The component of this anomalous Green’s
function for a given Matsubara frequency $\omega_{m}=\pi T(2m+1)$ is
calculated as follows:

$\displaystyle
F_{\bm{i},\alpha\beta}(\omega_{m})=\sum\limits_{n}(\frac{u_{n,\alpha}^{\bm{i}}v_{n,\beta}^{\bm{i}*}}{i\omega_{m}-\varepsilon_{n}}+\frac{u_{n,\beta}^{\bm{i}}v_{n,\alpha}^{\bm{i}*}}{i\omega_{m}+\varepsilon_{n}})$
(54)

Only off-diagonal in spin space components, corresponding to opposite-spin
pairs, are nonzero for the case under consideration. The singlet (triplet)
correlations are described by
$F_{\bm{i}}^{s,t}(\omega_{m})=F_{\bm{i},\uparrow\downarrow}(\omega_{m})\mp
F_{\bm{i},\downarrow\uparrow}(\omega_{m})$. Please note that the on-site
triplet correlations are odd in Matsubara frequency, as it should be according
to the general fermionic symmetry. Therefore we calculate the sum over the
positive Matsubara frequencies
$F_{\bm{i}}^{t}=\sum\limits_{\omega_{m}>0}F_{\bm{i}}^{t}(\omega_{m})$.

## References

  * Cooper [1956] L. N. Cooper, Bound electron pairs in a degenerate fermi gas, Phys. Rev. 104, 1189 (1956).
  * Bardeen  _et al._ [1957] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Microscopic theory of superconductivity, Phys. Rev. 106, 162 (1957).
  * Maki and Tsuneto [1964] K. Maki and T. Tsuneto, Pauli paramagnetism and superconducting state, Progress of Theoretical Physics 31, 945 (1964).
  * Maki [1969] K. Maki, Gapless superconductivity, in  _Superconductivity: Part 2_ , Superconductivity, edited by R. Parks (Taylor & Francis, 1969).
  * Bergeret  _et al._ [2005] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Odd triplet superconductivity and related phenomena in superconductor-ferromagnet structures, Rev. Mod. Phys. 77, 1321 (2005).
  * Buzdin [2005] A. I. Buzdin, Proximity effects in superconductor-ferromagnet heterostructures, Rev. Mod. Phys. 77, 935 (2005).
  * Bergeret  _et al._ [2018] F. S. Bergeret, M. Silaev, P. Virtanen, and T. T. Heikkilä, Colloquium: Nonequilibrium effects in superconductors with a spin-splitting field, Rev. Mod. Phys. 90, 041001 (2018).
  * Chandrasekhar [1962] B. S. Chandrasekhar, A note on the maximum critical field of high‐field superconductors, Applied Physics Letters 1, 7 (1962).
  * Clogston [1962] A. M. Clogston, Upper limit for the critical field in hard superconductors, Phys. Rev. Lett. 9, 266 (1962).
  * Sarma [1963] G. Sarma, On the influence of a uniform exchange field acting on the spins of the conduction electrons in a superconductor, Journal of Physics and Chemistry of Solids 24, 1029 (1963).
  * Kamra  _et al._ [2018] A. Kamra, A. Rezaei, and W. Belzig, Spin splitting induced in a superconductor by an antiferromagnetic insulator, Phys. Rev. Lett. 121, 247702 (2018).
  * Hauser  _et al._ [1966] J. J. Hauser, H. C. Theuerer, and N. R. Werthamer, Proximity effects between superconducting and magnetic films, Phys. Rev. 142, 118 (1966).
  * Bobkova  _et al._ [2005] I. V. Bobkova, P. J. Hirschfeld, and Y. S. Barash, Spin-dependent quasiparticle reflection and bound states at interfaces with itinerant antiferromagnets, Phys. Rev. Lett. 94, 037005 (2005).
  * Andersen  _et al._ [2005] B. M. Andersen, I. V. Bobkova, P. J. Hirschfeld, and Y. S. Barash, Bound states at the interface between antiferromagnets and superconductors, Phys. Rev. B 72, 184510 (2005).
  * Andersen  _et al._ [2006] B. M. Andersen, I. V. Bobkova, P. J. Hirschfeld, and Y. S. Barash, $0-\pi$ transitions in josephson junctions with antiferromagnetic interlayers, Phys. Rev. Lett. 96, 117005 (2006).
  * Bulaevskii  _et al._ [2017] L. Bulaevskii, R. Eneias, and A. Ferraz, Superconductor-antiferromagnet-superconductor $\pi$ josephson junction based on an antiferromagnetic barrier, Phys. Rev. B 95, 104513 (2017).
  * Enoksen  _et al._ [2013] H. Enoksen, J. Linder, and A. Sudbø, Pressure-induced 0-$\pi$ transitions and supercurrent crossover in antiferromagnetic weak links, Phys. Rev. B 88, 214512 (2013).
  * Rabinovich  _et al._ [2019] D. S. Rabinovich, I. V. Bobkova, and A. M. Bobkov, Anomalous phase shift in a josephson junction via an antiferromagnetic interlayer, Phys. Rev. Research 1, 033095 (2019).
  * Bell  _et al._ [2003] C. Bell, E. J. Tarte, G. Burnell, C. W. Leung, D.-J. Kang, and M. G. Blamire, Proximity and josephson effects in superconductor/antiferromagnetic $\mathrm{Nb}/\gamma-{\mathrm{fe}}_{50}{\mathrm{mn}}_{50}$ heterostructures, Phys. Rev. B 68, 144517 (2003).
  * Komissinskiy  _et al._ [2007] P. Komissinskiy, G. A. Ovsyannikov, I. V. Borisenko, Y. V. Kislinskii, K. Y. Constantinian, A. V. Zaitsev, and D. Winkler, Josephson effect in hybrid oxide heterostructures with an antiferromagnetic layer, Phys. Rev. Lett. 99, 017004 (2007).
  * Zhou  _et al._ [2019] X. Zhou, M. Lan, Y. Ye, Y. Feng, X. Zhai, L. Gong, H. Wang, J. Zhao, and Y. Xu, Andreev reflection and $0-\pi$ transition in graphene-based antiferromagnetic superconducting junctions, EPL (Europhysics Letters) 125, 37001 (2019).
  * Hübener  _et al._ [2002] M. Hübener, D. Tikhonov, I. A. Garifullin, K. Westerholt, and H. Zabel, The antiferromagnet/superconductor proximity effect in Cr/V/Cr trilayers, Journal Of Physics-Condensed Matter 14, 8687 (2002).
  * Wu  _et al._ [2013] B. L. Wu, Y. M. Yang, Z. B. Guo, Y. H. Wu, and J. J. Qiu, Suppression of superconductivity in nb by irmn in irmn/nb bilayers, Applied Physics Letters 103, 152602 (2013).
  * Seeger  _et al._ [2021] R. L. Seeger, G. Forestier, O. Gladii, M. Leiviskä, S. Auffret, I. Joumard, C. Gomez, M. Rubio-Roy, A. I. Buzdin, M. Houzet, and V. Baltz, Penetration depth of cooper pairs in the irmn antiferromagnet, Phys. Rev. B 104, 054413 (2021).
  * Buzdin and Bulaevskiĭ [1986] A. I. Buzdin and L. N. Bulaevskiĭ, Antiferromagnetic superconductors, Soviet Physics Uspekhi 29, 412 (1986).
  * Krivoruchko [1996] V. N. Krivoruchko, Upper critical fields of the superconducting state of a superconductor-antiferromagnetic metal superlattice, Zh. Eksp.Teor. Fiz. 109, 649 (1996).
  * Belashchenko [2010] K. D. Belashchenko, Equilibrium magnetization at the boundary of a magnetoelectric antiferromagnet, Phys. Rev. Lett. 105, 147204 (2010).
  * He  _et al._ [2010] X. He, Y. Wang, N. Wu, A. N. Caruso, E. Vescovo, K. D. Belashchenko, P. A. Dowben, and C. Binek, Robust isothermal electric control of exchange bias at room temperature, Nature Materials 9, 579 (2010).
  * Manna and Yusuf [2014] P. Manna and S. Yusuf, Two interface effects: Exchange bias and magnetic proximity, Physics Reports 535, 61 (2014).
  * Abrikosov and Gor’kov [1961] A. A. Abrikosov and L. P. Gor’kov, Contribution to the theory of superconducting alloys with paramagnetic impurities, Sov. Phys. JETP 12, 1243 (1961).
  * Johnsen  _et al._ [2021] L. G. Johnsen, S. H. Jacobsen, and J. Linder, Magnetic control of superconducting heterostructures using compensated antiferromagnets, Phys. Rev. B 103, L060505 (2021).
  * Kopnin [2001] N. Kopnin,  _Theory of Nonequilibrium Superconductivity_ , International Series of Monographs (Clarendon Press, 2001).
  * Belzig  _et al._ [1999] W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, Quasiclassical green’s function approach to mesoscopic superconductivity, Superlattices and Microstructures 25, 1251 (1999).
  * Zhu [2016] J. Zhu,  _Bogoliubov-de Gennes Method and Its Applications_, Lecture Notes in Physics (Springer International Publishing, 2016).
  * Baltz  _et al._ [2018] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Antiferromagnetic spintronics, Rev. Mod. Phys. 90, 015005 (2018).
  * Cheng  _et al._ [2014] R. Cheng, J. Xiao, Q. Niu, and A. Brataas, Spin pumping and spin-transfer torques in antiferromagnets, Phys. Rev. Lett. 113, 057601 (2014).
  * Takei  _et al._ [2014] S. Takei, B. I. Halperin, A. Yacoby, and Y. Tserkovnyak, Superfluid spin transport through antiferromagnetic insulators, Phys. Rev. B 90, 094408 (2014).
  * Simensen  _et al._ [2020] H. T. Simensen, A. Kamra, R. E. Troncoso, and A. Brataas, Magnon decay theory of gilbert damping in metallic antiferromagnets, Phys. Rev. B 101, 020403 (2020).
  * Eschrig  _et al._ [2015] M. Eschrig, A. Cottet, W. Belzig, and J. Linder, General boundary conditions for quasiclassical theory of superconductivity in the diffusive limit: application to strongly spin-polarized systems, New Journal Of Physics 17, 083037 (2015).
  * Eschrig [2015] M. Eschrig, Spin-polarized supercurrents for spintronics: a review of current progress, Reports On Progress In Physics 78, 104501 (2015).
  * Linder and Balatsky [2019] J. Linder and A. V. Balatsky, Odd-frequency superconductivity, Rev. Mod. Phys. 91, 045005 (2019).
  * Bobkova and Bobkov [2014] I. V. Bobkova and A. M. Bobkov, Bistable state in superconductor/ferromagnet heterostructures, Phys. Rev. B 89, 224501 (2014).