# Density functional Bogoliubov-de Gennes theory for superconductors
implemented in the SIESTA code

R. Reho1, N. Wittemeier2, A. H. Kole1, P. Ordejón2, Z. Zanolli  Chemistry
Department, Debye Institute for Nanomaterials Science and ETSF, Condensed
Matter and Interfaces, Utrecht University, PO Box 80.000, 3508 TA Utrecht, The
Netherlands.  Catalan Institute of Nanoscience and Nanotechnology - ICN2 (CSIC
and BIST), Campus UAB, Bellaterra, 08193 Barcelona, Spain

###### Abstract

We present SIESTA-BdG, an implementation of the simultaneous solution of the
Bogoliubov-de Gennes (BdG) and Density Functional Theory (DFT) problem in
SIESTA, a first-principles method and code for material simulations which uses
pseudopotentials and a localized basis set. This unified approach describes
both conventional and unconventional superconducting states, and enables a
description of inhomogeneous superconductors and heterostructures. We
demonstrate the validity, accuracy, and efficiency of SIESTA-BdG by computing
physically relevant quantities (superconducting charge density, band
structure, superconducting gap features, density of states) for conventional
singlet (Nb, Pb) and unconventional (FeSe) superconductors. We find excellent
agreement with experiments and results obtained within the KKR-BdG
computational framework. SIESTA-BdG forms the basis for modelling quantum
transport in superconducting devices and including - in an approximate fashion
- the superconducting DFT potential of Oliveira, Gross, and Kohn.

##  I Introduction

Superconductivity is an intrinsically quantum phenomenon, defying the
classical picture of dissipative current flow. It is one of the striking cases
where quantum effects manifest at the macroscopic scale, as demonstrated by
the Meissner effect [1]. The microscopic theory developed by Bardeen, Cooper,
and Schrieffer (BCS) in 1957 [2, 3] was instrumental in understanding the
underlying mechanism in terms of attractive interaction between electrons
mediated by crystal lattice vibrations (phonons), forming Cooper pairs below a
(very low) critical temperature. However, there are other complex realizations
of superconductivity which cannot be explained within the BCS theory. These
are referred to as unconventional or exotic forms of superconductivity.
Examples include high-temperature superconductivity and multiband
superconductivity, where multiple bands cross the Fermi surface, leading to
several superconducting instabilities [4, 5, 6, 7]. These and other exotic
forms of superconductivity can arise due to strong interactions between
electrons [8], from the interplay between disorder and dimensionality [9], or
emerge from the proximity interaction of a conventional superconductor and a
topological insulator (topological superconductivity) [10]. In particular,
topological superconductors [11, 12, 13, 14] could revolutionize the field of
quantum computing by hosting quantum bits based on Majorana zero modes [15].

Modelling and understanding the origins of effects like proximity-induced
superconductivity and unconventional superconductivity requires covering one
or more of the following points: (a) a realistic description of
heterostructures and interfaces between normal and superconducting materials
(for example, to model topological superconductors) (b) a unified treatment of
normal and superconducting electronic properties to predict band structure,
density of states, normal and anomalous charge density, and (c) the ability to
model unconventional pairing mechanisms between electrons (e.g. $p-$wave spin-
triplet), beyond BCS theory. Superconducting density functional theory (SCDFT)
[16] satisfies these requirements: it combines density functional theory (DFT)
[17, 18], which successfully predicts structural and electronic properties of
solids, with the Bogoliubov-de Gennes (BdG) formalism [19, 20, 21], which is a
Hartree-Fock-like mean-field approximation of the many-body Hamiltonian
generalizing the BCS theory [2] to treat inhomogeneous superconductors. SCDFT
allows for a general form of the pairing interaction and, hence, can describe
normal/superconductor interfaces and unconventional superconductors [22, 23],
in addition to conventional ones [24, 25, 26].

In SCDFT, a non-local superconductive effective potential (depending on both
normal and anomalous charge densities) replaces the standard DFT one (Section
II.1 and Appendix A). Then, the energy functional’s minimization with respect
to both normal and anomalous densities leads to the BdG equations. SCDFT is
formally an exact theory. However, it does not give a prescription for the
construction of the superconducting coupling (kernel) $\lambda$. In this work,
we follow the semi-phenomenological approach of Suvasini et al [27, 28], where
the superconducting kernel is assumed to be local and the explicit functional
dependence on the normal ($\rho$) and anomalous ($\chi$) densities is replaced
by site-dependent constants (Section II.1), namely we neglect the dynamical
update of the superconducting kernel due to quasiparticle interactions (such
as magnons or phonons).

Semi-phenomenological SCDFT has been recently implemented within multiple
scattering theory in the Korringa-Kohn-Rostoker (KKR) approximation (KKR-BdG)
[29, 30] to study conventional superconductors [30], proximity-induced
superconductivity [31, 32, 33], and the formation of Yu-Shiba-Rusinov states
(YSR) [34, 35, 36]. A different approach, combining dynamical mean field
theory (DMFT) and the quasiparticle self-consistent GW approximation, has been
recently developed and implemented in the Questaal code [37] and successfully
used to model (un)conventional superconductors [37, 38, 39]. In Questaal, the
superconducting properties are extracted from pairing susceptibilities in the
particle-particle channel, computed with the Bethe-Salpeter equation.

In this paper, we present the theoretical description and software
implementation of semi-empirical SCDFT in the SIESTA DFT method and code [40,
41], which we will refer to as SIESTA-BdG. SIESTA is a well-documented, open
source code with an exhaustive set of tutorials 111https://github.com/siesta-
project and https://www.youtube.com/@siestaproject, and a large community of
developers and users. SIESTA uses a basis of strictly localized orbitals,
which leads to sparse matrices which can be diagonalized using highly
efficient linear algebra algorithms: SIESTA has a native interface for the
ScaLAPACK, MRRR and ELPA solvers. These solvers and recently developed GPU-
accelerated interface in ELPA can also be used in SIESTA-BdG. In addition, the
ELSI library allows simultaneous parallelization over k-points, orbitals, and
spins, which is particularly useful given the extremely dense sampling of the
Brillouin Zone required to capture the (sub)meV energy scales characterizing
the superconducting gap. All these features make SIESTA-BdG an ideal choice
for studying superconductivity in systems consisting of several hundreds of
atoms such as interfaces, or heterostructures.

In SIESTA-BdG, the specific value of the superconducting coupling can be taken
from experiment, simulations (electron-phonon DFT, strong correlation model,
…), or chosen to explore the properties of the system as a function of the
input parameters. Then, one can perform a self-consistent calculation of both
normal and superconducting (anomalous) charge densities to compute the pairing
potential and superconducting properties, such as local density of states
(directly comparable with scanning tunnelling microscopy experiments), shape
and size of the superconducting gap, quasiparticle spectrum, electron and hole
wave functions, Fermi surface, and electron momentum distribution function
(Section II.2). SIESTA-BdG can be employed to analyse the role that different
atomic orbitals play in the formation of Cooper pairs by selectively
“activating” the pairing between specific orbitals. This feature is useful to
investigate the mechanism behind exotic forms of superconductivity. SIESTA-BdG
is implemented using the Spin-Orbit coupling (SOC) formalism present in SIESTA
[43], so that SOC and related effects (such as topology and topological
superconductivity) can be fully described in synergy with normal state
properties. Hence, within SIESTA, one can compute the full spectrum of
structural, electronic, magnetic, and topological properties of both normal
and superconducting materials.

Quantum electron transport can also be computed using the same software
package, ensuring a coherent description of all these properties for the whole
simulation workflow. The SIESTA-BdG method is an ideal starting point for
superconducting quantum transport calculations using non-equilibrium Green’s
functions techniques (NEGF), because the strictly localized basis sets used in
SIESTA make it easy to define local Hamiltonians and local Green’s functions.
Indeed, the SIESTA package comes with a set of routines and utilities that
implement methods to compute quantum transport properties in normal state
device, namely TranSiesta [44, 45]. In this work, we show how the SIESTA-BdG
method can be applied to study superconducting quantum transport in a (5,5)
armchair carbon nanotube (CNT) using NEGF methods in a one-particle framework
[46, 47, 48, 49, 50]. We compute the normal and Andreev transmission and
reflection functions for an artificial N-S junction.

The manuscript is organized as follows. In Section II we adapt the Bogoliubov-
de Gennes formalism [20] to linear combination of atomic orbitals (LCAO)
methods, present the SIESTA-BdG implementation, three different self-
consistency schemes to solve the SIESTA-BdG equations, and two methods for the
initialization of the superconducting pairing potential and coupling. In
Section III we present SIESTA-BdG simulations for conventional (Nb) and
unconventional superconductors (FeSe), and validate them against experimental
and/or theoretical simulations with KKR-BdG. In Section IV, we formulate
electron transport in the presence of superconductivity within the SIESTA-BdG
method and discuss the foreseen impact on the field. In the appendices we
discuss the spin generalized superconducting DFT formalism (Appendix A),
computational strategies to perform SIESTA-BdG simulations (Appendix B)
including efficiency, accuracy, SIESTA-BdG solution methods, initialization of
the superconducting potential, and computational details. Adaptation of
scattering theory for transport in superconductors to the SIESTA-BdG formalism
is covered in Appendix C.

##  II SIESTA-BdG

###  II.1 Bogoliubov-de Gennes equations in a localized basis set

In the SIESTA-BdG method, the spin generalized Kohn-Sham BdG equations \-
reviewed in Appendix A \- are written in terms of the SIESTA linear
combination of atomic orbitals (LCAO) basis set [40]. An atomic orbital
$\phi_{\mu}$ is identified by an atom label ($I$) and its quantum numbers
($l,m,n$), collected in a super-index $\mu=\\{I,l,m,n\\}$ which labels the
basis orbitals in the unit cell:

$\varphi_{\mu}^{\alpha}(\mathbf{r})=\varphi_{Ilmn\alpha}(\mathbf{r})=\varphi_{Ilmn}(\mathbf{r})\otimes|\alpha\rangle$
(1)

with $\alpha$ the spin index. The electron $u_{i}^{\alpha}(\mathbf{r})$ and
hole $v_{i}^{\alpha}(\mathbf{r})$ eigenfunctions in the LCAO basis set can
still be grouped into four-component Nambu spinors
$\Psi_{i}(\mathbf{r})=({u}_{i}^{\uparrow}\left(\mathbf{r}\right)u_{i}^{\downarrow}\left(\mathbf{r}\right)v_{i}^{\uparrow}\left(\mathbf{r}\right)v_{i}^{\downarrow}\left(\mathbf{r}\right))^{T}$,
and read:

$u_{i}^{\alpha}(\mathbf{r})=\sum_{\mu}\varphi_{\mu}(\mathbf{r})u_{i\mu}^{\alpha}~{},\quad
v_{i}^{\alpha}(\mathbf{r})=\sum_{\mu}\varphi_{\mu}^{*}(\mathbf{r})v_{i\mu}^{\alpha}$
(2)

where $u_{i\mu}^{\alpha}$ and $v_{i\mu}^{\alpha}$ are the expansion
coefficients, and $i$ the index of the superconducting bands (Bogoliubon
excitation state).

SIESTA solves the DFT problem using periodic boundary conditions. Non-periodic
systems (molecules, tubes, slabs, interfaces, …) are modelled using the
supercell approach. Assuming a lattice periodic pairing potential, it is
convenient to write and solve the BdG eigenvalue problem in the reciprocal
space of the simulation cell. Exploiting the Bloch theorem, the electron and
hole wavefunction can be expanded in Bloch waves

$u^{\alpha}_{i}(\mathbf{k},\mathbf{r})={}\sum_{\mu,\mathbf{R}}e^{i\mathbf{k}(\mathbf{R}+\bm{\tau_{\mu}})}\varphi_{\mu}(\mathbf{r})u_{i\mu}^{\alpha}(\mathbf{k})~{},\quad
v^{\alpha}_{i}(\mathbf{k},\mathbf{r})={}\sum_{\mu,\mathbf{R}}e^{i\mathbf{k}(\mathbf{R}+\bm{\tau_{\mu}})}\varphi^{*}_{\mu}(\mathbf{r})v_{i\mu}^{\alpha}(\mathbf{k})$
(3)

where the sum runs through all orbitals in space, given by the lattice vectors
$\mathbf{R}$ and orbitals positions $\bm{\tau}_{\mu}$ in the unit cell, and
$\mathbf{k}$ is the crystal momentum of the supercell. Writing the BdG
equations in the LCAO basis is crucial to reveal the contribution of specific
orbitals to superconductivity, explain complex interactions between different
bands (multi-band superconductors), and analyse competing superconducting
phases as those breaking time-reversal symmetry or related to different spin
configurations [51]. The relevant equations for the SIESTA-BdG method in
$\mathbf{k}$-space are:

  1. 1.

the BdG Hamiltonian $\mathbf{\widehat{H}_{BdG}}$

$\widehat{H}_{BdG,\mu\nu}^{\alpha,\beta}(\mathbf{k})=\sum_{\nu,\mathbf{R}}e^{i\mathbf{k}(\mathbf{R}+\bm{\tau}_{\nu}-\bm{\tau}_{\mu})}\widehat{H}^{\alpha,\beta}_{BdG,\mu\nu}$
(4)

  2. 2.

the eigenvalue problem

$\displaystyle\sum_{\mu\alpha}\left(\begin{array}[]{cc}\mathbf{h^{KS}}^{\beta\alpha}_{\nu\mu}(\mathbf{k})-\mu
S_{\nu\mu}(\mathbf{k})\delta_{\beta\alpha}&\Delta^{\beta\alpha}_{\nu\mu}(\mathbf{k})\\\
-\Delta^{*\beta\alpha}_{\nu\mu}(-\mathbf{k})&-\mathbf{h^{KS}}^{*\beta\alpha}_{\nu\mu}(-\mathbf{k})+\mu
S_{\nu\mu}(\mathbf{k})\delta_{\beta\alpha}\end{array}\right)\left(\begin{array}[]{c}u^{\alpha}_{i\mu}(\mathbf{k})\\\
v_{i\mu}^{\alpha}(\mathbf{k})\end{array}\right)$ (9)
$\displaystyle=\varepsilon_{i}\sum_{\mu\alpha}S_{\nu\mu}(\mathbf{k})\delta_{\beta\alpha}\left(\begin{array}[]{c}u^{\alpha}_{i\mu}(\mathbf{k})\\\
v_{i\mu}^{\alpha}(\mathbf{k})\end{array}\right),$ (12)

where $\bf{h^{KS}}$ is the normal state Kohn-Sham Hamiltonian,
$\Delta^{\beta\alpha}_{\nu\mu}(\mathbf{k})$ is the superconducting pairing
potential, $S_{\nu\mu}(\mathbf{k})$ the overlap matrix at $\mathbf{k}$ and
$\mu$ the chemical potential,

  3. 3.

the normal and anomalous densities

$\displaystyle\rho^{\alpha\beta}_{\mu\nu}(\mathbf{k})=\sum_{i}f\left(\varepsilon_{i}(\mathbf{k})\right)u_{i\mu}^{*\alpha}(\mathbf{k})u_{i\nu}^{\beta}(\mathbf{k})+\sum_{i}\bar{f}(\varepsilon_{i}(\mathbf{-k}))v_{i\mu}^{\alpha}(\mathbf{-k})v_{i\nu}^{*\beta}(\mathbf{-k})$
(13)
$\displaystyle\chi^{\alpha\beta}_{\mu\nu}(\mathbf{k})=\sum_{i}f\left(\varepsilon_{i}(\mathbf{k})\right)u_{i\nu}^{\beta}(\mathbf{k})v_{n\mu}^{*\alpha}(\mathbf{k})+\sum_{i}\bar{f}\left(\varepsilon_{i}(\mathbf{-k})\right)u_{i\mu}^{\alpha}(\mathbf{-k})v_{i\nu}^{*\beta}(\mathbf{-k})$
(14)

where $f(\varepsilon_{i})$ and $\bar{f}(\varepsilon_{i})=1-f(\varepsilon_{i})$
are the occupation functions,

  4. 4.

the pairing potential $\Delta^{\alpha\beta}_{\mu\nu}(\mathbf{k})$:

$\Delta^{\alpha\beta}_{\mu\nu}(\mathbf{k})=-\sum_{l=0,1,2,3}\sum_{\gamma}i\sigma_{l}^{\alpha\gamma}\sigma_{2}^{\gamma\beta}\int_{\mathrm{BZ}}\mathrm{d}\mathbf{k^{\prime}}\lambda_{l}^{\mu\nu}(\mathbf{k}-\mathbf{k^{\prime}})\chi^{l}_{\mu\nu}(\mathbf{k^{\prime}})$
(15)

where $\sigma_{l}$ are the Pauli matrices, $\lambda_{l}^{\mu\nu}$ are the
orbital projected superconducting couplings. We denote as $\bm{\lambda}$ the
matrix with elements $\lambda_{l}^{\mu\nu}$.

The above equations are only valid under the assumption that the pairing
potential is lattice periodic. Superconductors with spatially non-uniform
pairing potential need to be modelled using a supercell approach. This is the
case of Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase of a superconductor
[52], which is characterized by Cooper pairs with finite momentum.

Even though the superconducting coupling $\bm{\lambda}$ is local, it can have
an explicit dependence on the atomic orbitals in real-space, which can be used
to initialize the simulations (Sections II.3 and B.4). The pairing potential
$\Delta$ is updated in the self-consistent steps via 15 (Section II.2) and it
acquires a non-trivial wave vector $\mathbf{k}$-dependence. Hence, the SIESTA-
BdG formalism allows one to compute non-isotropic superconducting band gap and
superconductors with multiple $\bf{k}$-dependent gaps (multiband
superconductors).

In the following, the quantities computed in real and reciprocal space are
denoted with a $\mathbf{R}$ or $\mathbf{k}$ superscript, respectively (e.g.
$\mathbf{H}^{\mathbf{R}}$, $\mathbf{H}^{\mathbf{k}}$). The matrix elements of
any function in real space $g(\mathbf{r})$ can be written in the LCAO basis as
$g^{\mathbf{R}}_{\mu,\nu}=\sum_{\mathbf{R}^{\prime}}\int_{\text{UC}+\mathbf{R}^{\prime}}\mathrm{d^{3}}\mathbf{r}\phi^{*}_{\mu}(\mathbf{r})g(\mathbf{r})\phi_{\nu}(\mathbf{r}-\mathbf{R})$.

###  II.2 SIESTA-BdG SCF Solution Methods

Preliminary to SIESTA-BdG simulations, we perform a normal state DFT
calculation to compute the electronic Hamiltonian $\bf{h^{KS}}$ and normal
state density $\rho$. The converged $\bf{h^{KS}}$ is used to construct the BdG
Hamiltonian $\mathbf{H}_{BdG,\mu\nu}^{\mathbf{R}}$ by introducing a pairing
potential $\bm{\Delta}$ to break charge conservation and enable the
superconducting phase. In SIESTA-BdG, $\bm{\Delta}$ and $\bm{\lambda}$ can be
initialized in the orbital or real space representation, as discussed in
Section II.3. Then, the SIESTA-BdG Hamiltonian and (normal and anomalous)
densities can be computed by solving the eigenvalue problem (Eq. 9) either
directly (non SCF-BdG) or by performing a self-consistent field cycle (SCF) on
only the normal density (fixed-$\Delta$) or both the normal and anomalous
densities (full SCF-BdG).

For a given choice of $\bm{\Delta}$ or $\bm{\lambda}$ full self-consistency
should also include the chemical potential in addition to the densities. The
number of electrons is given by a weighted sum over all states, using a
(chemical potential dependent) occupation function as the weighting factor.
For normal-state calculations, only the occupation function depends on the
chemical potential ($\mu$) while for SIESTA-BdG calculations, $\mu$ explicitly
enters the eigenvalue problem. This means that a change in $\mu$ requires one
to solve the eigenvalue problem again to find new solutions, making the
process more complex. To compute $\mu$ we follow the approach by Suvasini et
al. [28] where $\mu$ is updated in an outer loop that repeats the SCF cycle
until the correct number of particles is reached [28]. The user can skip this
loop to reduce computational costs. In this case, the total number of
electrons might not match the number of electrons in the normal state, but
experience shows that the deviation is small for typical values of the pairing
potential [28]. The SIESTA-BdG solution methods are illustrated in Figure 1
and described below.

STARTDFT SCF: Find self-consistent solution for the normal-state electronic
density $\rho_{0}$Construct $h^{\mathrm{KS}}_{\mu\nu}[\rho_{0}]$non SCF-
BdGmethod?Solve eigenvalue problem:
$\left(\begin{array}[]{cc}h^{KS}_{\mu\nu}&\Delta_{\mu\nu}\\\
-\Delta^{*}_{\mu\nu}&-h^{KS}*_{\mu\nu}\end{array}\right)\left(\begin{array}[]{c}u_{\nu}\\\
v_{\nu}\end{array}\right)=\varepsilon\left(\begin{array}[]{cc}\mathbf{S}_{\mu\nu}&0\\\
0&\mathbf{S}_{\mu\nu}\end{array}\right)\left(\begin{array}[]{c}u_{\nu}\\\
v_{\nu}\end{array}\right)$Construct $\rho$ and $\chi$ matrices from
$\left(\begin{array}[]{c}u_{\nu}\\\ v_{\nu}\end{array}\right)$Fixed-$\Delta$
or full SCF-BdG?Fixed-$\Delta$Construct $h^{KS}_{\mu\nu}$Self-consistent$\rho$
and $h^{KS}_{\mu\nu}$?full SCF-BdGConstruct $h^{KS}_{\mu\nu}$ and
$\Delta_{\mu\nu}$Self-consistent$\rho$, $h^{KS}_{\mu\nu}$, $\chi$,
$\Delta$?Charge conserved?Update Fermi level?Solve eigenvalue problem:
$\left(\begin{array}[]{cc}h^{KS}_{\mu\nu}&\Delta_{\mu\nu}\\\
-\Delta^{*}_{\mu\nu}&-h^{KS}*_{\mu\nu}\end{array}\right)\left(\begin{array}[]{c}u_{\nu}\\\
v_{\nu}\end{array}\right)=\varepsilon\left(\begin{array}[]{cc}\mathbf{S}_{\mu\nu}&0\\\
0&\mathbf{S}_{\mu\nu}\end{array}\right)\left(\begin{array}[]{c}u_{\nu}\\\
v_{\nu}\end{array}\right)$ for bands, PDOS, $\chi$, etc.Input
$\Delta_{\mu\nu}$Yesskip SIESTA-BdGSCF loopNoSIESTA-BdG SCF loopYesNoInput
$\lambda$YesNoNoYesYesNo

Figure 1: Schematic of the iterative process to compute self-consistent
solutions of the SIESTA-BdG secular equation. During the update step within
the SIESTA-BdG SCF loop the density or the Hamiltonian is mixed.

Non SCF-BdG: The $\mathbf{H}_{BdG}^{\mathbf{R}}$ Hamiltonian is constructed
and diagonalized (Eq. 4). No SIESTA-BdG SCF loop is performed; therefore,
$\bf{h^{KS}}$ is not updated in the presence of $\bm{\Delta}$. The resulting
eigenvalues and eigenvectors are post-processed to calculate electronic band
structure, density of states, normal and anomalous charge density. The
underlying assumption is that superconductivity is a small perturbation of
$\bf{h^{KS}}$ and the normal state density. The non SCF-BdG method has a
significantly lower computational cost compared to fixed-$\Delta$ and full
SCF-BdG, as the latter two methods perform an additional SCF cycle using the
four-component Nambu spinors $\Psi_{i}(\mathbf{r})$.

Fixed-$\bf{\Delta}$: The $\mathbf{H}_{BdG}^{\mathbf{R}}$ Hamiltonian is
diagonalized in the presence of a fixed pairing potential $\bm{\Delta}$. At
every step of the SIESTA-BdG SCF cycle, the Kohn-Sham Hamiltonian
($\bf{h^{KS}}$), the normal ($\rho$) and the anomalous ($\chi$) densities are
updated. Self-consistency is ensured on the normal-state part ($\bf{h^{KS}}$,
$\rho$) of the SIESTA-BdG Hamiltonian. The underlying assumption is that the
perturbations induced on the superconducting properties due to the anomalous
charge are negligible, namely the initial guess for $\bm{\Delta}$ is good
enough to model the superconducting system and does not need to be updated in
the SCF cycle via Eq. 15. In particular, for bulk Nb the difference between
fixed-$\Delta$ and non SCF-BdG in the density of states are negligible for
values of the pairing potential $\mathbf{\Delta}$ smaller than $\sim$ 100 meV
(Appendix B.3).

full SCF-BdG: is an implementation of the full self-consistent field SIESTA-
BdG scheme. The superconducting couplings $\bm{\lambda}$ and the external
pairing potential $d^{i}_{\mathrm{ext}}(\mathbf{r},\mathbf{r^{\prime}})$ (Eq.
39) must be specified in the input. The external pairing potential is needed
to ensure a non-zero anomalous density in the first SCF step, and is set to
zero thereafter (see discussion after Eq. 39). At each SCF step,
$\bf{h^{KS}}$, $\rho$, $\chi$ and $\bm{\Delta}=\sum_{i}\left(d_{\text{ext
}}^{i}\left(\mathbf{r}\right)-\lambda^{i}\left(\mathbf{r}\right)\chi^{i}\left(\mathbf{r}\right)\right)\sigma_{i}i\sigma_{2}$
are updated. Analogously to normal-state simulations, $\bm{\Delta}$ and $\chi$
can be mixed using generalized SIESTA algorithms. Convergence can be achieved
on $\bm{\Delta}$, on $\chi$, or on both of them simultaneously.

###  II.3 Initialization of the pairing potential and superconducting coupling

Depending on the solution method, either the superconducting coupling
$\bm{\lambda}$ (full SCF-BdG) or the pairing potential $\bm{\Delta}$ (non SCF-
BdG and fixed-$\Delta$) is a semi-empirical parameter. Their initial values
can be expressed either in terms of the LCAO atomic orbitals (orbital
representation) or as a modulation in real space of a superconducting strength
parameter $\bar{\lambda}$ or $\bar{\Delta}$ (superconducting strength
representation) by a function $g(\mathbf{r})$, i.e.
$\Delta(\mathbf{r})=\bar{\bm{\Delta}}g(\mathbf{r})$. Note that, when working
in this representation, we use a bold symbol to indicate real-space matrices
with orbital projected matrix elements (as $\bm{\lambda}$) and italics for
real space functions $\lambda(\mathbf{r})$. The orbital representation is
ideal for investigating the nature of the superconducting interaction by
selectively coupling specific orbitals. The user can directly specify the
matrix elements of the pairing potential ($d^{i\mathbf{R}}_{\mu\nu}$) for any
superconducting channel ($i$). The full pairing potential matrix is
constructed using the orbital representation of Eq. 40:

$\displaystyle\bm{\Delta}^{\mathbf{R}}_{\mu\nu}=\sum_{i=0,1,2,3}d^{i\mathbf{R}}_{\mu\nu}\sigma_{i}i\sigma_{2}=\sum_{i=0,1,2,3}\left(d_{\text{ext},\mu\nu}^{i\mathbf{R}}-\lambda^{\mathbf{R}\mu\nu}_{i}\chi^{i\mathbf{R}}_{\mu\nu}\right)\sigma_{i}i\sigma_{2}$
(16)

If the symmetries of the gap function are known, then the initial guess for
$\mathbf{\Delta}$ can be constructed as a matrix respecting those symmetries.
For instance, one can specify a pairing exclusively between electrons and
holes occupying the same orbital on the same atom (intra-orbital, diagonal
$\bm{\Delta}$), or between different electron and hole orbitals belonging to
the same atom (on-site, block-diagonal $\bm{\Delta}$), or couple all the
orbitals of all atoms (full overlap, full $\bm{\Delta}$ matrix), as
illustrated in Figure 2. The flexibility of the orbital representation has
been used to identify the orbitals responsible for unconventional
superconductivity in FeSe (Section III.2).

Figure 2: Initialization of the superconducting pairing in the orbital
representation: Schematic representation of coupling between orbitals in the
(a) intra-orbital, (b) on-site, and (c) full overlap pairing.

The superconducting strength representation has a more intuitive physical
meaning, as the superconducting pairing is defined on the SIESTA real-space
grid. In this case, the user can choose among the different options depending
on the knowledge of the physical properties and symmetries of the materials to
be studied. The computed physical observables (for instance, $\chi$) mainly
depends on the spatial basis overlap $S_{\mu\nu}^{\mathbf{R}}$ and not on
$g(\mathbf{r})$, as illustrated in Appendix B.4 for bulk Pb. The current
implementation includes the following options for the modulation function
$g(\mathbf{r})$:

  * •

Gaussians:
$g(\mathbf{r};\mathbf{x},\sigma)=e^{-|\mathbf{r}-\mathbf{x}|^{2}/(2\sigma^{2})}$,
with center $\mathbf{x}$ and width $\sigma$.

  * •

Spherical wells:
$g(\mathbf{r};\mathbf{x},\sigma)=\theta(|\mathbf{r}-\mathbf{x}|;\sigma,\eta)$
with center $\mathbf{x}$ and radius $\sigma$. The shape $\theta$ can be a
Heaviside step function (hard well) or a Fermi-Dirac function with broadening
$\eta$ (soft well).

  * •

Plane waves: $g(\mathbf{r};\mathbf{k},\phi)=\cos(\mathbf{k}\mathbf{r}+\phi)$,
where $\mathbf{k}$ is the wave vector and $\phi$ is the phase of the wave at
$\mathbf{x}=0$.

  * •

any combination of the above

Gaussian and spherical wells can be centred at any point within the unit cell
and have arbitrary widths, potentially extending beyond the unit cell
boundary. The lattice periodicity is automatically imposed. Plane waves can be
used to represent fully delocalized interactions. Superpositions of multiple
waves can be used to create pairing potentials with different symmetries. The
versatility of the initial choice can be used to model complex systems with
inhomogeneous superconducting potential.

##  III Results

###  III.1 Bulk Niobium

Figure 3: Bulk Nb: (a) Density of States (DOS) comparison among full SCF-BdG
SIESTA, JuKKR [30, 53] and STM measurements [54], showing good agreement
between simulations and experiment. For SIESTA-BdG we set $\lambda=0.188$
eV$\cdot V_{\mathrm{sphere}}$ with a spherical ($r=1.4136$ Å) square well
potential and a smearing temperature $T=0.1$ K. In the JuKKR code a value of
$\lambda=1.11$ eV was used [30]. (b) Superconducting band structure with
visible particle-hole symmetry.

Bulk Niobium is a good system for validating the SIESTA-BdG method, as several
theoretical [28, 30] and experimental [55, 56, 54] results are available in
the literature for comparison. We computed the superconducting state DOS using
the full SCF-BdG method and compared the SIESTA-BdG result with the JuKKR code
[30] and experimental data from STM measurements [54] (Figure 3). There is
good agreement between the two codes: both predict a DOS with two peaks at
$\sim\pm 2$ meV (main coherence peaks), characteristic of conventional
superconductors, and denoted U-shaped gap. In addition, the two codes predict
visible gap anisotropy, as revealed by the multiple peaks in the DOS for low
enough smearing temperatures. The two approaches differ in the position of the
innermost peaks, the relative shift is about $\pm 0.2$ meV, which can be
explained by differences in the assumptions underlying the two codes. The
overall gap shape also matches the one from the experiment. The missing gap
anisotropy in the experimental data could be due to temperature (smearing),
SIESTA-BdG predicts an additional gap anisotropy (i.e. each main coherence
peak splits in two) compared to JuKKR, possibly due to the finer sampling of
the energy axis in our calculation.

###  III.2 Bulk FeSe

Figure 4: DOS of FeSe, an unconventional superconductor, for different pairing
potentials calculated with the SIESTA-BdG method: (a) isotropic pairing
($\Delta(\mathbf{r})$ = 3 meV), (b-f) pairing potential limited to a different
Fe-$d$ orbital ($\Delta_{d_{\nu}d_{\nu}}=3$ meV). The simulations with
isotropic pairing (a) reproduce the size and shape of the experimentally
observed superconducting gap [57]. The superconducting gap emerges only if the
pairing potential is applied to the $d_{xz}$ (c) or $d_{yz}$ (d) orbitals.

The origin and character of high-temperature superconductivity in FeSe (and
other iron-based superconductors) are currently debated [58, 59, 60]. It has
been suggested that FeSe is an unconventional $d$-wave superconductor [61, 62,
63] driven by pairing interaction in the $3d_{xz}$ and $3d_{yz}$ orbitals of
iron [64]. Experimental evidence suggests that FeSe does not possess long-
range magnetic order. Even at low temperatures (4 K), spin fluctuations are
significant [51]. These fluctuations occur primarily between two magnetic
phases: Néel (or checkerboard) and staggered dimer. In the Néel magnetic
structure, the magnetic moments of neighbouring iron atoms are antiparallel.
In the following, we report results for the Néel magnetic structure. The
simulation of the nematic phase would require very large supercells in
periodic boundary conditions, which is beyond the scope of this work.

The superconducting gap of FeSe in the Néel phase is V shaped and of the order
of $\sim 1-10$ meV at the coherence peaks, depending on the experimental
conditions [57]. To model the superconducting state we use the non SCF-BdG
method and initialize the pairing potential with a constant
($\Delta(\mathbf{r})=\bar{\Delta}=3$ meV). The correct V shape of the
superconducting gap emerges from our calculations (Figure 4.a). This indicates
that orbital symmetries are sufficient to predict the conventional or
unconventional nature of the superconducting gap.

To analyse the role of different orbitals in enabling superconductivity, we
restrict the pairing interaction to specific orbitals: We use the orbital
representation method (Section II.3) and assign non-zero values to matrix
elements corresponding to a specific Fe-$3d$ orbital (Fe-$3d_{xz}$,
Fe-$3d_{yz}$, …) at a time. We observe that imposing a pairing potential
exclusively between $3d_{xz}$ or $3d_{yz}$ orbitals opens a V-shaped gap
(Figure 4 b-f). We can therefore identify those orbitals as responsible for
unconventional superconductivity in FeSe. The superconducting gap is smaller
than the one computed with isotropic pairing among all orbitals because the
pairing interaction is limited to a subset of all orbitals, reducing the
pairing potential’s average strength. No superconducting gap emerges if the
pairing interaction is restricted to any of the other Fe-$3d$ orbitals.

##  IV Quantum Transport within the SIESTA-BdG formalism

###  IV.1 The formalism

The SCDFT formalism maps the many-body problem, including the pair interaction
giving rise to Cooper pairs, onto an equivalent non-interacting problem
described by one-particle wave functions. In this auxiliary system, many
aspects of quantum transport can be described in one-particle frameworks in a
Landauer-based conductivity approach. As is commonly done in DFT approaches,
the KS wavefunctions are assumed to be a good approximation to the one-
particle wavefunction of the interacting many-body system [65, 48, 49, 50].
However, any current described in such a one-particle framework is
fundamentally a one-particle current. As such, Cooper pair currents are not
described.

In a non-superconducting device, the current through any electrode
$\mathfrak{e}$ at temperature $T_{\mathfrak{e}}$ and chemical potential
$\mu_{\mathfrak{e}}$ is given by [66]

$\displaystyle
I_{\mathfrak{e}}=\frac{e}{h}\int_{-\infty}^{\infty}d\varepsilon\,\int_{\mathrm{BZ}}d\mathbf{k}\,\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T_{\mathfrak{e}\mathfrak{e}^{\prime}}[f_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)]$
(17)

where $T_{\mathfrak{e}\mathfrak{e}^{\prime}}$ is the transmission function
between two electrodes $\mathfrak{e}$ and $\mathfrak{e}^{\prime}$, and
$f_{\mathfrak{e}}(\varepsilon)$ is the occupation function for the electrode
$\mathfrak{e}$. The energy and $\mathbf{k}$ dependence of the quantities is
implicit.

In the case of superconducting electrodes or devices, the single-particle
current can be carried either by the electrons or holes. The transmission
($T_{\mathfrak{e}\mathfrak{e}^{\prime}}$) and reflection probabilities
($R_{\mathfrak{e}}$) become 2x2 matrices

$\displaystyle R_{\mathfrak{e}}$
$\displaystyle\rightarrow\begin{pmatrix}R_{\mathfrak{e}}^{pp}&R_{\mathfrak{e}}^{hp}\\\
R_{\mathfrak{e}}^{ph}&R_{\mathfrak{e}}^{hh}\end{pmatrix}$ $\displaystyle
T_{\mathfrak{e}\mathfrak{e}^{\prime}}$
$\displaystyle\rightarrow\begin{pmatrix}T_{\mathfrak{e}\mathfrak{e}^{\prime}}^{pp}&T_{\mathfrak{e}\mathfrak{e}^{\prime}}^{hp}\\\
T_{\mathfrak{e}\mathfrak{e}^{\prime}}^{ph}&T_{\mathfrak{e}\mathfrak{e}^{\prime}}^{hh}\end{pmatrix}$
(18)

with $R_{\mathfrak{e}}^{pp}$ ($R_{\mathfrak{e}}^{hh}$) the reflection
probability for electrons (holes) in lead $\mathfrak{e}$ to scatter back into
electronic (hole) states on the same electrode, and $R_{\mathfrak{e}}^{ph}$
($R_{\mathfrak{e}}^{hp}$) the reflection probability for an electron (hole) in
lead $\mathfrak{e}$ to scatter into hole (electron) states in the same lead
(Andreev reflection). Similarly, the four components of
$T_{\mathfrak{e}\mathfrak{e}^{\prime}}$ describe the scattering probabilities
between electron and hole states in lead $\mathfrak{e}$ and electron and hole
states in $\mathfrak{e}^{\prime}$ (normal and anomalous transmission). The
transmission and reflection functions can be calculated using the Green’s
function formulation of quantum transport [66, 67]. Using all four functions
(see Appendix C for additional details), the current through any (normal)
electrode can be written as [68]

$\displaystyle
I_{\mathfrak{e}}=\frac{2e}{h}\int_{-\infty}^{\infty}d\varepsilon\,\int_{\mathrm{BZ}}d\mathbf{k}\,\Big{\\{}$
$\displaystyle\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}[f_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)]-T^{hh}_{\mathfrak{e}\mathfrak{e}^{\prime}}[\bar{f}_{\mathfrak{e}}(\varepsilon)-\bar{f}_{\mathfrak{e}^{\prime}}(\varepsilon)]$
$\displaystyle+$
$\displaystyle\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}[f_{\mathfrak{e}}(\varepsilon)-\bar{f}_{\mathfrak{e}^{\prime}}(\varepsilon)]-T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}[\bar{f}_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)]$
$\displaystyle+$
$\displaystyle[R^{hp}_{\mathfrak{e}}+R^{ph}_{\mathfrak{e}}][f_{\mathfrak{e}}(\varepsilon)-\bar{f}_{\mathfrak{e}}(\varepsilon)]\Big{\\}}$
(19)

where $\bar{f}_{\mathfrak{e}}(\varepsilon)$ is the occupation of hole states
in the electrode $\mathfrak{e}$. Moving to the low-bias limit
$f_{\mathfrak{e}}(\varepsilon)=f(\varepsilon)-\frac{\partial{f}}{\partial{\varepsilon}}(\varepsilon)[\mu_{\mathfrak{e}}-\mu]$
($\mu$ the reference chemical potential w.r.t to which the electrode
potentials are shifted) and invoking particle-hole symmetry, the expression
for the current can be rewritten as:

$\displaystyle I_{\mathfrak{e}}\approx$
$\displaystyle\frac{2e}{h}\int_{\mathrm{BZ}}d\mathbf{k}\,\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\mu)[\mu_{\mathfrak{e}^{\prime}}-\mu_{\mathfrak{e}}]+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\mu)[2\mu-\mu_{\mathfrak{e}^{\prime}}-\mu_{\mathfrak{e}}]+2R^{hp}_{\mathfrak{e}}(\mu)[\mu-\mu_{\mathfrak{e}}]$
(20)

In this formula three different components of the current can be identified:
(i) a normal component which depends on electron-electron (hole-hole)
transmission ($T^{pp}$, $T^{hh}$), (ii) an anomalous component which depends
on electron-hole (hole-electron) transmission ($T^{ph}$, $T^{hp}$), and (iii)
an Andreev contribution which depends solely on the Andreev reflections
($R^{hp}$, $R^{ph}$). Particle conservation requires that the currents flowing
in all electrodes sum up to zero (Kirchhoff’s current law [69]):
$\sum_{\mathfrak{e}}I_{\mathfrak{e}}=0$ [48, 70]. This condition uniquely
determines the chemical potential $\mu$: The current is not proportional to
the pairwise differences in the electrode chemical potentials, and explicitly
depends on the electrode chemical potentials. This is an important difference
with respect to a normal-state device.

Eq. IV.1 is valid for computing the current of any normal electrode with a
superconducting device. However, it is not suitable for a superconducting
electrode $\mathfrak{e}$, as it does not account for the contributions of
Cooper pair currents. Arguments on how to properly compute currents under non-
equilibrium conditions and in superconducting electrodes are given in [50, 48,
70] and references therein.

A non-zero conductance at zero applied voltage (zero-bias conductance) can
only arise due to Andreev (or anomalous) terms in this formalism. In a bulk
superconductor, no reflections occur, and the anomalous transmission is equal
to zero. Thus, the zero-bias conductance in any bulk superconductor is zero.
Only in the presence of a scattering region can Andreev reflections (or
anomalous transmission) occur and lead to a non-zero conductance. This, of
course, does not imply a finite resistance of the bulk superconductor, since
Cooper pair currents are not included in the description.

The SIESTA-BdG method, due to the strictly localized basis set, constitute an
ideal starting point for quantum transport calculations using non-equilibrium
Green’s functions techniques (NEGF) on open quantum systems with
superconducting components, such as superconductor junctions and topological
superconductor devices. The chemical potential in the leads does not change as
a consequence of the interaction with the scattering region [66].

Figure 5:  Normal-Superconductor (N-S) junction in an (5,5) carbon nanotube.
The normal ($\Delta=0$ meV) and superconducting ($\Delta=5$ meV) regions are
indicated in blue and orange, respectively. (a) Simulation cell for transport
simulations displaying scattering region and leads (shaded area) (b)
Electronic band structure of the scattering region of the N-S (5, 5) carbon
nanotube junction. (c) Transmission and (d) reflection in quantum conductance
units $G_{0}$. Within the superconductor gap, all electrons in the normal
electrode participate in Andreev reflections ($R^{eh}$). For energies larger
than the superconducting gap, the N-S junction behaves like a normal metal. At
the superconducting gap there is a smooth transition between the two regimes.

To showcase the possible application of the SIESTA-BdG method in the quantum
transport formalism, we consider a (5,5) armchair carbon nanotube (CNT) as a
transport system. Signatures of superconductivity have been observed in
different types of carbon nanotubes: (a) 4 Å single-wall carbon nanotubes at
$T_{C}\approx 15$ K (Meissner effect) [71], (b) multi-wall nanotubes with 10
nm outer tube at $T_{C}\approx 12$ K (sudden drop of DOS below $T_{C}$) [72],
and (c) ropes of single-wall carbon nanotubes at 0.55 K (drop of resistance)
[73]. However, it is difficult to reliably create superconductivity in one-
dimensional nanostructures due to low DOS and quantum fluctuations [74]. If
realized and stabilized, superconducting carbon nanotubes could be ideal one-
dimensional leads for nanoelectronics.

We construct an N-S junction within a (5,5) nanotube by selectively applying
the pairing potential to part of the scattering region and right electrode.
Our transport device consists of 3 rings of the (5,5) nanotube (scattering
region) which are contacted by the same nanotube on each side (electrodes,
each 3 rings wide), thus forming an infinite tube. We use the Hamiltonian
obtained from a periodic calculation neglecting the coupling with the infinite
leads. We then use the Green’s function formulation of quantum transport [66]
to obtain the transmission and reflection probabilities. We apply the non SCF-
BdG method and neglect any changes in the electronic structure of the leads
due to the superconducting pairing (Landauer assumption). Therefore, the
resulting N-S junction possesses a negligible barrier for the electrons and
electron-electron back-scattering should be close to non-existent. Thus, our
N-S junction is expected to behave like an ideal N-S junction, allowing for
Andreev reflections that can lead to non-vanishing reflection coefficients.
Furthermore, we use an isotropic superconducting pairing with a strength of 5
meV.

In the normal state, the band structure of a (5,5) armchair CNT near the Fermi
level (blue iti7lines in Figure 5.a) is described by four linear bands (two
per spin channel). In the superconducting state (orange lines, Figure 5.a), a
band gap opens up at the Fermi level. The band dispersion becomes parabolic,
and the number of bands is doubled. In addition, Figure 5 shows the different
transmission and reflection functions for the N-S junction between a normal
state CNT and superconducting CNT. The Andreev reflections $R^{ph}$ and normal
transmission $T^{pp}$ (Figures 5.b and 5.c) become significant at the edge of
the superconducting gap where electron and hole states hybridize in the
scattering region because of the pairing potential. Deep within the gap, the
anomalous transmission $T^{ph}$ and the normal reflection $R^{pp}$ are zero,
and all incoming states are perfectly reflected (Figures 5.b and 5.d). These
results exactly replicate the expected behaviour for a barrier-free N-S
junction [75].

## Conclusions and Outlook

In this work we present SIESTA-BdG, a method to calculate the superconducting
properties of materials and its implementation in the open-source SIESTA code
[41]. SIESTA-BdG implements the semi-phenomenological SCDFT formalism of
Suvasini [28] in localized basis sets (equations 9 and 15). We show that
SIESTA-BdG can accurately predict both conventional and unconventional
superconductors.

We propose three solution methods to solve the DFT-BdG secular eqution 9: non
SCF-BdG, fixed-$\Delta$ and full SCF-BdG, which mainly differ in the treatment
of self-consistency of the superconducting density and Hamiltonian. We present
different options for initializing the superconducting coupling/pairing, which
is a parameter in this method. We sustained our discussion with three test
cases: Nb, Pb and FeSe bulk, which demonstrate the robustness and accuracy of
our implementation. We critically compared our results with experiments and
KKR-BdG simulations, finding excellent quantitative agreement (Section III).
We introduced and discussed optimization schemes such as the adaptive
Brillouin Zone sampling, the tetrahedron method, parallelization over
$\mathbf{k}$-points and real space initialization of the pairing potential
that are indispensable for an accurate description of the superconducting
state. We conclude the manuscript with a perspective on the simulation of
quantum transport for superconductors by combining the SIESTA-BdG with the
NEGF formalism of Ref. [44, 45]. As a proof of concept, we compute
transmission and reflection probabilities of a normal/superconducting junction
in a carbon nanotube.

The limitations of SIESTA-BdG are related to the ignorance of the exact form
of the superconducting coupling ($\bm{\lambda}$) or pairing ($\bm{\Delta}$),
namely the magnitude and phase of the superconducting interaction is unknown.
One can make an estimate using full SCDFT for conventional superconductors, or
refer to theories that model the superconducting coupling strength itself [76,
77]. As a guideline to initialize the pairing potential/coupling, we suggest
using first the superconducting strength representation with a constant
modulation function to acquire familiarity with the system and, then, test
more specific real-space initializations. When the symmetries of the system
are known, it might be easier to impose them directly in the matrix elements
of the pairing potential (orbital representation, Section II.2). Exploring the
consequences of the initialization of the superconducting parameter (different
solution methods, new or more complex real space initialization, etc.) needs
to be addressed case-by-case, as it is material specific.

The applicability of our approach is vast, heterostructures and complex device
configurations being prime candidates that could benefit from this novel
implementation, given SIESTA’s inherent optimal scaling. In perspective,
SIESTA-BdG framework can be interfaced with existing SCDFT codes, mapping the
non-local SCDFT potential into an approximate local one and/or extended to a
non-local superconducting potential to allow for SCDFT calculations. SIESTA-
BdG is the basis for the description of transport in superconducting systems:
In general, the SCDFT formalism maps the many-body problem, including the pair
interaction giving rise to Cooper pairs, onto an equivalent non-interacting
problem described by one-particle wavefunctions. In this auxiliary system,
transport is described by a one-particle framework (e.g. scattering matrix
formalism or Green’s function formalisms [78]), as opposed to other Many-Body
formalism which focuses on the two-particle superconducting response functions
[38, 39, 79]. The current implementation will provide a useful tool to further
investigate and understand exotic forms of superconductivity and related
emerging physics such as unconventional superconductivity, proximity induced
superconductivity, topological superconductivity, interplay between magnetism
and superconductivity.

## Acknowledgments

The authors acknowledge the fruitful discussion with Philipp Rüßmann and
Andrés R. Botello-Mendéz, Miguel Pruneda, and Gábor Csire, Antonio Sanna,
Hervé Ness, and Enrico Perfetto. ZZ acknowledges the research program
“Materials for the Quantum Age” (QuMat) for financial support. This program
(registration number 024.005.006) is part of the Gravitation program financed
by the Dutch Ministry of Education, Culture and Science (OCW). RR and AK
financial support from Sector Plan Program 2019-2023. This work was sponsored
by NWO-Domain Science for the use of supercomputer facilities. We also
acknowledge that the results of this research have been achieved using the
Tier-0 PRACE Research Infrastructure resource Discoverer based in Sofia,
Bulgaria (OptoSpin project id. 2020225411). This project has received funding
from the European Union’s Horizon Europe research and innovation programme
under grant agreement No 101130384 (QUONDENSATE). PO and NW acknowledge
support from the EU MaX CoE (Grant No. 101093374) and grants PCI2022-134972-2
and PID2022-139776NB-C62 funded by the Spanish MCIN/AEI/10.13039/501100011033
and by the ERDF, A way of making Europe, and Grant No. NW further acknowledges
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska-Curie Grant Agreement No. 754558 (PREBIST
– COFUND). ICN2 is supported by the Severo Ochoa programme from Spanish MINECO
(grant no. CEX2021-001214-S) and by Generalitat de Catalunya (CERCA
programme).

## Data availability statement

The data that support the findings of this study are openly available on
Materials Cloud.

##  Appendix A Spin generalized Superconducting DFT

The SCDFT Hamiltonian
$\widehat{H}_{\mathrm{BdG}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ [16]
consists of four blocks describing particles (top-left), antiparticles
(bottom-right), and particle/antiparticle coupling $\Delta$ (off-diagonal)
222We interchangeably use the term electrons (holes) and particles (anti-
particles).:

$\displaystyle\widehat{H}_{\mathrm{BdG}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=$
$\displaystyle\left(\begin{array}[]{cc}\widehat{h}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r^{\prime}})&\Delta(\mathbf{r},\mathbf{r^{\prime}})\\\
\Delta^{*}(\mathbf{r},\mathbf{r^{\prime}})&-\widehat{h}^{*}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r^{\prime}})\end{array}\right)$
(23)

To extend the $2\times 2$ SCFDT Hamiltonian to the $4\times 4$ fully
relativistic case (i.e. including SOC), we follow the SIESTA formalism [41]:

$\displaystyle\widehat{H}_{\mathrm{BdG},\alpha\beta}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=$
$\displaystyle\left(\begin{array}[]{cc}\widehat{h}^{KS}_{\alpha\beta}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r^{\prime}})&\Delta_{\alpha\beta}(\mathbf{r},\mathbf{r^{\prime}})\\\
-\Delta_{\alpha\beta}^{*}(\mathbf{r},\mathbf{r^{\prime}})&-\widehat{h}^{*KS}_{\alpha\beta}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r^{\prime}})\end{array}\right)$
(26)

with spin indices $\alpha$ and $\beta$, and

$\widehat{h}^{KS}(\mathbf{r})=\left(-\nabla^{2}+\widehat{V}^{H}-\mu\right)\mathds{1}_{2}+\widehat{V}^{\mathrm{KB}}+\widehat{V}^{\mathrm{SO}}+\widehat{V}^{XC}$
(27)

the $2\times 2$ fully relativistic electronic Kohn-Sham Hamiltonian. The
latter includes the kinetic energy $\nabla^{2}$, the Hartree potential
$\widehat{V}^{H}$, the chemical potential $\mu$, the Kleinman–Bylander scalar-
relativistic pseudo-potential $\hat{V}^{KB}$, the spin-orbit term
$\hat{V}^{SO}$, and the exchange-correlation potential $\widehat{V}^{XC}$. The
relativistic contributions to the pairing potential $\bm{\Delta}$ are
neglected. The minus sign in the bottom-left block of Eq. (26) arises from the
anticommutation rules of the spin indices:
$-\Delta^{*}_{\alpha,\beta}=\Delta^{*}_{\beta,\alpha}$. Writing explicitly the
spin dependence, it reads:

$\widehat{H}_{BdG}(\mathbf{r},\mathbf{r^{\prime}})=\left(\begin{array}[]{cccc}\widehat{h}^{KS}_{\uparrow\uparrow}(\mathbf{r})&\widehat{h}^{KS}_{\uparrow\downarrow}(\mathbf{r})&\Delta_{\uparrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&\Delta_{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\\\
\widehat{h}^{KS}_{\downarrow\uparrow}(\mathbf{r})&\widehat{h}^{KS}_{\downarrow\downarrow}(\mathbf{r})&\Delta_{\downarrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&\Delta_{\downarrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\\\
-\Delta_{\uparrow\uparrow}^{*}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&-\Delta_{\uparrow\downarrow}^{*}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&-\widehat{h}^{*KS}_{\uparrow\uparrow}(\mathbf{r})&-\widehat{h}^{*KS}_{\uparrow\downarrow}(\mathbf{r})\\\
-\Delta_{\downarrow\uparrow}^{*}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&-\Delta_{\downarrow\downarrow}^{*}\left(\mathbf{r},\mathbf{r}^{\prime}\right)&-\widehat{h}^{*KS}_{\downarrow\uparrow}(\mathbf{r})&-\widehat{h}^{*KS}_{\downarrow\downarrow}(\mathbf{r})\end{array}\right)$
(28)

The spin generalized Kohn-Sham BdG equations are, then:

$\int\mathrm{d}\mathbf{r}^{\prime}\widehat{H}_{\mathrm{BdG}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\left(\begin{array}[]{c}u_{i}^{\uparrow}\left(\mathbf{r}^{\prime}\right)\\\
u_{i}^{\downarrow}\left(\mathbf{r}^{\prime}\right)\\\
v_{i}^{\uparrow}\left(\mathbf{r}^{\prime}\right)\\\
v_{i}^{\downarrow}\left(\mathbf{r}^{\prime}\right)\end{array}\right)=\varepsilon_{i}\left(\begin{array}[]{c}u_{i}^{\uparrow}(\mathbf{r})\\\
u_{i}^{\downarrow}(\mathbf{r})\\\ v_{i}^{\uparrow}(\mathbf{r})\\\
v_{i}^{\downarrow}(\mathbf{r})\end{array}\right).$ (29)

The eigenfunctions of the BdG Hamiltonian are four-component Nambu spinors
$\Psi_{i}(\mathbf{r})=({u}_{i}^{\uparrow}\left(\mathbf{r}\right)u_{i}^{\downarrow}\left(\mathbf{r}\right)v_{i}^{\uparrow}\left(\mathbf{r}\right)v_{i}^{\downarrow}\left(\mathbf{r}\right))^{T}$,
where $u^{\alpha}_{i}(r)$ and $v^{\alpha}_{i}(r)$ are the electron and hole
components. In the superconducting ground state ($\varepsilon_{0}$) all Cooper
pairs are in the same energy level. The eigenenergies $\varepsilon_{i}$ above
the ground state correspond to Bogoliubons, namely quasiparticles consisting
of a linear combination of excited electrons and holes. The BdG equations (26)
satisfy particle-hole symmetry [15]. Namely, if
$(u_{i}^{\uparrow}\left(\mathbf{r}\right)u_{i}^{\downarrow}\left(\mathbf{r}\right)v_{i}^{\uparrow}\left(\mathbf{r}\right)v_{i}^{\downarrow}\left(\mathbf{r}\right))^{T}$
is an eigenvector with solution $\varepsilon_{i}$, then
$(v_{i}^{\uparrow}\left(\mathbf{r}\right)v_{i}^{\downarrow}\left(\mathbf{r}\right)u_{i}^{\uparrow}\left(\mathbf{r}\right)u_{i}^{\downarrow}\left(\mathbf{r}\right))^{*T}$
leads to an equivalent solution with eigenvalue $-\varepsilon_{i}$.

In analogy with Kohn-Sham DFT, one can define single particle ($V$) and
pairing ($\bf{\Delta}$) effective potentials in terms of a superconducting
exchange-correlation functional $\Omega_{xc}[\rho,\bm{\chi}]$, which is a
functional of both the normal $\rho(\mathbf{r})$ and the anomalous
$\bm{\chi}(\mathbf{r},\mathbf{r^{\prime}})$ charge densities :

$\displaystyle V[\rho,\bm{\chi}](\mathbf{r})$
$\displaystyle=V^{H}(\mathbf{r})+V^{KB}(\mathbf{r})+V^{SO}(\mathbf{r})+\frac{\delta\Omega_{xc}[\rho,\bm{\chi}]}{\delta\rho(\mathbf{r})}$
(30)
$\displaystyle\bm{\Delta}[\rho,\bm{\chi}]\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
$\displaystyle=d^{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\mathrm{i}\hat{\sigma}_{2}+\sum_{j=1,2,3}d^{j}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\hat{\sigma}_{j}\mathrm{i}\hat{\sigma}_{2}$
(31) $\displaystyle
d^{i}[\rho,\bm{\chi}]\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
$\displaystyle=d_{\mathrm{ext}}^{i}\left(\mathbf{r},\mathbf{r}^{\prime}\right)-\frac{\delta\Omega_{xc}[\rho,\bm{\chi}]}{\delta\chi^{*i}\left(\mathbf{r},\mathbf{r}^{\prime}\right)}~{}~{}~{}\mathrm{with}~{}i=0,1,2,3$
(32)

where $\hat{\sigma}_{i}$ are the Pauli matrices, $d^{i}$ are the singlet
($i=0$) and triplet ($i=1,2,3$) spin components in the Balian-Werthamer (BW)
representations [81], and
$d^{i}_{\mathrm{ext}}(\mathbf{r},\mathbf{r^{\prime}})$ is an external pairing
potential introduced, formally, to break the charge symmetry of the system. In
practice, the SIESTA-BdG implementation builds the BdG Hamiltonian from
$d^{i}_{\mathrm{ext}}(\mathbf{r},\mathbf{r^{\prime}})$ in the first self-
consistent field iteration. The resulting anomalous density is used to build
the pairing potential in the successive step and
$d^{i}_{\mathrm{ext}}(\mathbf{r},\mathbf{r^{\prime}})$ is set to zero
thereafter. For this reason,
$d^{i}_{\mathrm{ext}}(\mathbf{r},\mathbf{r^{\prime}})$ is missing in
equation(15).

In the BW representation, the singlet and triplet spin components of the
anomalous charge density are

$\displaystyle\chi^{0}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{\chi^{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)-\chi^{\downarrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)}{2}=\frac{\chi^{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)+\chi^{\uparrow\downarrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)}{2},$
(33)
$\displaystyle\chi^{1}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{\chi^{\downarrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)-\chi^{\uparrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)}{2}=\frac{\chi^{\uparrow\uparrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)-\chi^{\downarrow\downarrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)}{2},$
(34)
$\displaystyle\chi^{2}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{\chi^{\downarrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)+\chi^{\uparrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)}{2\mathrm{i}}=\frac{\chi^{\downarrow\downarrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)+\chi^{\uparrow\uparrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)}{-2\mathrm{i}},$
(35)
$\displaystyle\chi^{3}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\frac{\chi^{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)+\chi^{\downarrow\uparrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)}{2}=\frac{\chi^{\uparrow\downarrow}\left(\mathbf{r},\mathbf{r}^{\prime}\right)-\chi^{\uparrow\downarrow}\left(\mathbf{r}^{\prime},\mathbf{r}\right)}{2}.$
(36)

The BW parametrization is useful to study the properties of $\bm{\chi}$ under
parity transformation $\mathbf{r}\rightarrow\mathbf{r^{\prime}}$. $\bm{\chi}$
is a fermionic object and, hence, must satisfy the Pauli exclusion principle,
namely it must be antisymmetric under exchange of the two electrons forming
the Cooper pair:
$\chi^{\beta\alpha}(\mathbf{r}^{\prime},\mathbf{r})=-\chi^{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})$.
This implies that if $\chi^{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})$ is
even (odd) under the parity transformation, then it has to be odd (even) under
spin exchange leading to a $S=0$ singlet state ($S=1$, triplet state).

The normal and anomalous densities that minimize the energy functional
$E[\rho,\chi]$ are [16]

$\begin{gathered}\rho(\mathbf{r})=\sum_{\alpha,i}f\left(\varepsilon_{i}\right)\left|u_{i}^{\alpha}(\mathbf{r})\right|^{2}+\sum_{\alpha,i}\bar{f}\left(\varepsilon_{i}\right)\left|v_{i}^{\alpha}(\mathbf{r})\right|^{2}\\\
\chi^{\alpha\beta}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=\sum_{i}f\left(\varepsilon_{i}\right)u_{i}^{\beta}\left(\mathbf{r}^{\prime}\right)v_{i}^{*\alpha}(\mathbf{r})+\sum_{i}\bar{f}\left(\varepsilon_{i}\right)u_{i}^{\alpha}(\mathbf{r})v_{i}^{*\beta}\left(\mathbf{r}^{\prime}\right)\end{gathered}$
(37)

where $f(\varepsilon_{i})$ is the occupation function of the $i^{th}$
Bogoliubon, and $\bar{f}(\varepsilon_{i})=1-f(\varepsilon_{i})$. In contrast
to conventional DFT, the sum over $\alpha$, and $i$ in Eq. 37 extends over all
states and is not limited by the Fermi level. Furthermore, the contribution of
the $uv^{*}$ terms in $\chi$ is expected to be non-zero only in an energy
region close to the Fermi surface [21].

In SCDFT, $\Omega_{xc}[\rho,\bm{\chi}]$ is written in terms of a
superconducting kernel which is, in general, a non-local pairing interaction
dependent on four spatial (and spins) variables:
$\lambda[\rho,\bm{\chi}]\left(\mathbf{r},\mathbf{r}^{\prime};\mathbf{r_{1}},\mathbf{r_{1}}^{\prime}\right)$.
In this work, we use Suvasini [28] assumption of a local superconducting
kernel $\lambda(\mathbf{r})$. As a consequence, the exchange-correlation
functional reduces to

$\Omega_{xc}[\rho,\bm{\chi}]=\Omega_{0,xc}[\rho]+\int\sum_{i}\mathrm{d}\mathbf{r}\chi^{i*}(\mathbf{r})\lambda^{i}(\mathbf{r})\chi^{i}(\mathbf{r})$
(38)

with $\Omega_{0,xc}$ the Kohn-Sham exchange-correlation potential, and

$d^{i}\left(\mathbf{r}\right)=d_{\text{ext
}}^{i}\left(\mathbf{r}\right)-\lambda^{i}\left(\mathbf{r}\right)\chi^{i}\left(\mathbf{r}\right)$
(39)

which is expanded in its tensor form in the BW representation as

$\bm{\Delta(\mathbf{r})}=\sum_{i=0,1,2,3}d^{i}(\mathbf{r})\sigma_{i}i\sigma_{2}=\sum_{i=0,1,2,3}[d_{\text{ext
}}^{i}\left(\mathbf{r}\right)-\lambda^{i}\left(\mathbf{r}\right)\chi^{i}\left(\mathbf{r}\right)]\sigma_{i}i\sigma_{2}$
(40)

##  Appendix B SIESTA-BdG simulations in practice

###  B.1 Convergence of SIESTA-BdG simulations

Figure 6: Convergence of SIESTA-BdG calculations for bulk Niobium. a)
Dependence of the DOS on smearing temperature at fixed $800\times 800\times
800$ $\mathbf{k}$-grid. Gap anisotropy is only resolved when the smearing
temperature is on the order of $0.1$ K or lower. b) DOS dependence on
$\mathbf{k}$-grid size ($n_{k}\times n_{k}\times n_{k}$) with a smearing
temperature $T=0.1$ K: an extremely dense sampling ($n_{k}=800$) is required
to converge the tails of the gap and resolve the coherence peaks.

A rigorous convergence study of simulation parameters is required for the
SIESTA-BdG method. This concerns the real-space grid cutoff, the SCF cycle
tolerances, the $\mathbf{k}-$grid for the SCF cycle and DOS/PDOS, and the
smearing of electron/hole occupations.

The cutoff value of the real-space grid is typically the same order of
magnitude as for normal-state SIESTA simulations, as it mostly depend on the
choice of the pseudopotentials and on the desired accuracy of the
representation of physical quantities (density, potential, charge) on the
SIESTA real-space grid.

SCF tolerances for non SCF-BdG and fixed-$\Delta$ are similar to fully-
relativistic normal state calculations, namely $10^{-5}$ for the density
matrix and $10^{-4}$ eV for the Hamiltonian). Full SCF-BdG calculations
require stricter SCF tolerance compared to normal state calculations by one to
two orders of magnitude.

The $\mathbf{k}$-grid for the SCF cycles can be converged in steps: first, for
the normal metal phase, then the superconducting-state DOS is converged (using
the normal-state SCF $\mathbf{k}$-grid), and finally the SCF $\mathbf{k}$-grid
is further refined against the shape and size of the superconducting gap.

The smearing temperature and DOS/PDOS sampling of the Brillouin zone need to
be pushed beyond the usual converged values for normal-state SIESTA
simulations. The smearing temperature is related to the Fermi-Dirac occupation
function for electrons and holes, as usually defined in first-principles DFT
simulations. The smaller the smearing temperature, the higher the energy
resolution in the DOS. The DOS can vary greatly with the smearing temperature,
as illustrated in Figure 6.a. Therefore, the DOS is a good quantity to be used
for converging the smearing temperature. As a rule of thumb, setting the
temperature to about one-tenth of the superconducting gap usually ensures
accurate simulations.

In order to obtain smooth and converged DOS/PDOS it is necessary to use
hundreds to thousands of $\mathbf{k}$-points in each direction. An under-
converged sampling results in noisy DOS, which manifests most clearly in the
tails that border the superconducting gap, as illustrated in Figure 6.b for
the DOS of bulk Nb computed with different $\mathbf{k}$-grids. Performing
SIESTA-BdG DOS/PDOS simulations at such extremely dense sampling of the
Brillouin Zone requires the use of the strategies for efficient sampling
described in Section B.2.

We extended these two strategies to also work for the SIESTA-BdG calculations
with Nambu spinors, structuring our algorithm analogously to the one already
present in SIESTA [41].

Uniform gridUniform grid + 1 zoomUniform grid + 2 zoomzoomzoom

Figure 7: Schematic illustration of the adaptive Brillouin Zone sampling
scheme. Starting with a coarse grid, the algorithm determines the points that
contribute to the DOS in the region of interest (points circled in orange). To
enhance the sampling, the contributing points are replaced by a denser grid
(we ”zoom” in on the interesting region). This procedure can be repeated to
iteratively improve the sampling.

###  B.2 Strategies for efficient sampling of the Brillouin Zone.

The parallelization techniques allow one to deal with numerous
$\mathbf{k}$-points by distributing the work among more processors. Another
approach is reducing the number of $\mathbf{k}$-points needed. For normal
state calculations, the tetrahedron method has been shown to be effective in
accurately calculating integrals over the Brillouin Zone with relatively few
points [82]. However, in a typical SIESTA-BdG calculation, the region of
interest is limited to a few meV around the Fermi energy. This can be
exploited by using an adaptive sampling scheme where a dense sampling grid is
only used in the regions of the Brillouin Zone close to the Fermi surface and
a coarse grid everywhere else. We implemented this as an iterative procedure,
schematically illustrated in Figure 7, using the following steps:

  1. 1.

diagonalize the BdG Hamiltonian $\mathbf{H}^{k}$ and find all the eigenvalues
using an initial uniform grid in $\bf{k}$-space,

  2. 2.

identify which points have at least one eigenvalue that lies in a user-
specified energy range,

  3. 3.

replace all these points by denser sub-grids centred on the points,

  4. 4.

diagonalize the Hamiltonian and find the eigenvalues using these new points,

  5. 5.

repeat steps 2 to 4 until the desired accuracy in the density of states is
reached.

As an illustration of the effectiveness of this procedure, we compute the DOS
of bulk Nb using the adaptive sampling and compare with the DOS computed with
a uniform $800\times 800\times 800$ $k$-grid (Figure 6.b). The results can be
seen in Figure 8. The adaptive procedure starts with a $100\times 100\times
100$ grid and at each iteration (repeated three times) contributing points
were replaced by a $2\times 2\times 2$ sub-grid. This procedure allows one to
effectively reach the uniform $800\times 800\times 800$ sampling but only in a
small energy window around the Fermi energy. The adaptive grid reproduces the
smooth DOS at a fraction of the computational cost.

Figure 8: The DOS of bulk Niobium computed with a uniform $800\times 800\times
800$ $\mathbf{k}$-grid and an adaptive grid. $\Delta$ has been set to 5 meV.
The adaptive grid was generated starting from a $100\times 100\times 100$
uniform grid and using, for three iterations, a $2\times 2\times 2$ sub-grid

###  B.3 Limit of applicability of non-SCF BdG solution method

To investigate the differences in results between using the non SCF-BdG and
fixed-$\Delta$ approaches, the density of states (DOS) of bulk Niobium was
computed for different magnitudes of the pairing potential $\Delta$ using both
methods. For each magnitude a real-space initialization of touching hard
spheres was used. The results of the two approaches were compared by
integrating the absolute differences between the curves. This difference was
divided by the integral of the DOS of fixed-$\Delta$ over the same energy
range to get a relative integrated difference. The results can be seen in
Figure 9. The differences in the DOS remain small (less than 2%) up to a
$\Delta$ of 100 meV. This is consistent with earlier statements that the
superconducting pairing potential is a small perturbation. For large values of
$\Delta$ ($\geq 1$eV) the difference become significant. This however lies
beyond the reasonable strength for Niobium and most known superconductors.
Nevertheless, we stress that one should verify which of the approaches is
sufficient for the system under study.

Figure 9: The relative integrated difference in the DOS of bulk Niobium for
non SCF-BdG and fixed-$\Delta$ calculations with different magnitudes of the
pairing potential $\Delta$.

###  B.4 Comparing SCF-BdG solution methods and real space initializations

In this Section, we use bulk Pb to compare the three different SIESTA-BdG SCF
solution methods introduced in Section II.2. We work in the superconducting
strength representation of Section II.3. For the full SCF-BdG method, we also
compare different ways of specifying the superconducting coupling
$\lambda(\mathbf{r})$.

Bulk Pb is a conventional $s$-wave superconductor, where the pairing potential
arises from electron-phonon interactions [83, 84, 85, 86, 87]. The
superconducting gap of lead is U-shaped, the density of states inside the
superconducting gap is flat and rises sharply at the edges of the gap
(coherence peaks) [84]. The size of the superconducting gap has been measured
between 2.3 and 2.5 meV, depending on the crystal orientation [83, 84].

Based on this knowledge, we use a singlet-pairing potential with a strength of
$\bar{\Delta}=$1.25 meV (half of the gap size). We assume this potential to be
uniform and isotropic in real space, i.e., in matrix form, the pairing
potential can be written as

$\Delta_{\mu\nu}(\mathbf{R})=\bar{\Delta}\left(\begin{array}[]{cc}0&S_{\mu\nu}(\mathbf{R})\\\
-S_{\mu\nu}(\mathbf{R})&0\end{array}\right)$ (41)

where $S_{\mu\nu}(\mathbf{R})$ is the overlap matrix.

We perform calculations with all three solution methods (non SCF-BdG,
fixed-$\Delta$, and full SCF-BdG) to obtain the DOS for the superconducting
phase of Pb. For the full SCF-BdG method in Figure 10.a, we define $\lambda$
isotropic in space and tune its strength to obtain a gap with the same width
as for the two other methods finding $\bar{\lambda}=0.315$ eV $\cdot V_{UC}$,
where $V_{UC}$ is the unit cell volume. The resulting superconducting gap
confirms the experiment in shape and magnitude (Figure 10.a) irrespective of
the solution method used. The coherence peaks at the edge of the gap are
sharp, and the density of states forms a small plateau between these peaks.
Simulations with a smearing temperatures of $0.002$ meV allow us to predict a
splitting of the coherence peak of $0.07$ meV, which agrees well with the
experimentally observed one ($0.1$ meV) (Figure 11).

Figure 10: DOS of bulk Pb calculated with the SIESTA-BdG method with smearing
temperature $0.01$ meV. (a) Comparison of solution methods non SCF-BdG,
fixed-$\Delta$, and full SCF-BdG, assuming isotropic pairing
$\Delta(\mathbf{r})=\bar{\Delta}=1.25$ meV;
$\lambda(\mathbf{r})=\bar{\lambda}=0.315$ eV$\cdot V_{UC}$. (b) Comparison
between real-space initialization methods within the full SCF-BdG approach
(constant, touching spherical hardwells, atom-centered Gaussian ($\sigma=0.3$
Å)). Figure 11: DOS of bulk Pb calculated with full SCF-BdG with smearing
temperature $0.01$ meV (red) and $0.002$ meV (blue) using the constant real-
space initialization method as in Figure 10.c. The plot shows that the
computed coherence peak splitting ($0.07$ meV) agrees well with the
experimental one $0.1$ meV. Figure 12:  Real-space visualization of (a – c)
the superconducting coupling parameter $\lambda(\mathbf{r})$, (d – f) the
pairing potential $\Delta(\mathbf{r})$, and (g – i) the anomalous density
$\chi(\mathbf{r})$ for bulk Pb. These quantities are computed for different
initializations of $\lambda(\mathbf{r})$ in the full SCF-BdG solution method.
The columns correspond to different initialization: (a) constant, (b) touching
spherical hardwell ($R=1.704$ Å), and (c) Gaussian shape ($\sigma=0.3$ Å)).
The intensity of each quantity is represented as a colour gradient. The dashed
parallelogram indicates the unit cell of Pb with atoms centred on the corners.

To assess the influence of the real-space initialization of the
superconducting coupling $\lambda\mathbf{r})$, we calculate the pairing
potential $\Delta(\mathbf{r})$ and the anomalous density $\chi(\mathbf{r})$ of
bulk lead using the full SCF-BdG approach and the following initializations:
(i) constant $\lambda(\mathbf{r})=\bar{\lambda}$, (ii) touching spherical
hardwell $\lambda(\mathbf{r})=\bar{\lambda}$ if
$|\mathbf{r}-\mathbf{r}_{\mathrm{Pb}}|<1.704$ Å, and (iii) Gaussian shape
$\lambda(\mathbf{r})=\bar{\lambda}e^{-|\mathbf{r}-\mathbf{r}_{\mathrm{Pb}}|^{2}/(2\sigma^{2})}$;
$\sigma=0.3$ Å , as illustrated in Figure 12. In each case, we tune the
parameter $\bar{\lambda}$ to obtain a gap with approximately $2.56$ meV width
(Figure 10.b), resulting in $\bar{\lambda}/V_{UC}=$0.315 eV, 0.350 eV, and
3.69 eV for (i), (ii), and (iii), respectively.

Despite the differences in the functional form of the superconducting
coupling, the computed physical properties (DOS and anomalous charge) are
qualitatively similar, indicating that the physically relevant information is
included in the overlap matrix, namely in the projection of $\lambda$ and/or
$\Delta$ on the orbitals close to the Fermi level. The width of the
superconducting gap is the same by construction. The extension of the
anomalous charge and the spatial representation of $\Delta$, instead, depend
on the choice of $\bar{\lambda}$ and the specific radius of the modulation
function. The case that differ most is the Gaussian initialization, as those
functions are more sharply localized on the atoms.

###  B.5 Computational Details

SIESTA and SIESTA-BdG: The simulations for bulk Nb, Pb, and FeSe were
performed with PBE [88] functional, norm-conserving pseudopotentials from the
Pseudo-Dojo database [89], and the default SIESTA double-$\zeta$ polarized
(DZP) basis set [40]. The crystal structure of each system was relaxed with a
maximum stress tolerance of $0.01$ GPa.

For bulk Nb, the KS equations were solved using $T=0.1$ K, a $120\times
120\times 120$ $\mathbf{k}$-grid, and a mesh cutoff of $600$ Ry. Structural
relaxation was performed with $T=1$ meV , and $40\times 40\times 40$
$\mathbf{k}$-grid. The DOS was sampled on a $400\times 400\times 400$
$\mathbf{k}$-grid and then refined iterating the Adaptive Brillouin Zone
sampling three times with a $2\times 2\times 2$ sub-grid. A spherical square
well potential with a constant value of $\lambda=0.188$ eV$\cdot
V_{\mathrm{sphere}}$ inside the sphere was used (radius $r=1.4136$ Å) to solve
the SIESTA-BdG equations.

For bulk Pb, the KS equations were solved using $T=1$K, a $60\times 60\times
60$ $\mathbf{k}$-grid, and a mesh cutoff of $2000$ Ry. The DOS was sampled on
a $300\times 300\times 300$ $\mathbf{k}$-grid and then refined with a $3\times
3\times 3$ sub-grid.

The FeSe simulations are converged for a smearing temperature of 0.5 meV, a
$29\times 29\times 29$ k-grid for KS-equations, and a mesh cutoff of 1300 Ry.
The DOS was sampled on a 480$\times$480$\times$120 k-grid.

Quantum transport: We solve the KS equations using a smearing temperature of 1
meV, a 1$\times$1$\times$151 k-grid, and a mesh cutoff of 800 Ry. The
simulations are performed with the LDA functional as parameterized by Perdew
and Zunger [90] and norm-conserving fully relativistic pseudopotentials from
the Pseudo-Dojo database [89]. The atomic positions were calculated according
to the structure of an ideal nanotube with lattice constant 2.46 Å. We perform
a SIESTA-BdG calculation using the non SCF-BdG method with an pairing of 5 meV
in the superconducting section of the scattering region and the right
electrode, and set the super conducting pairing to 0 in the normal conductor
section. While it is well known that the superconducting order parameter
decays as a power law in low-dimensional nanostructures
[GonzalezCooperpairPropagation2007, GonzalezCriticalCurrents2008], we enforce
a sharp decay here for a proof of concept. In order to calculate the quantum
conductance, we use a custom Python script (Supplementary Material).

##  Appendix C Scattering theory for quantum transport with the SIESTA-BdG
formalism

The relationship between current and voltage in superconductors is
fundamentally different from that encountered in normal conductors. In a non-
superconducting device, the current through any electrode is given by

$\displaystyle
I_{\mathfrak{e}}=\int_{-\infty}^{\infty}\\!\\!d\varepsilon[M_{\mathfrak{e}}(\varepsilon)-R_{\mathfrak{e}}(\varepsilon)]f_{\mathfrak{e}}(\varepsilon)-\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)f_{\mathfrak{e}^{\prime}}(\varepsilon)$
(42)

where $\mathbf{M}(\varepsilon)$ is the number of channels in the electrode
$\mathfrak{e}$ at energy $\varepsilon$, $\mathbf{R}(\varepsilon)$ is the
reflection function,
$\mathbf{T}_{\mathfrak{e},\mathfrak{e}^{\prime}}(\varepsilon)$ is the
transmission function between two electrodes $\mathfrak{e}$ and
$\mathfrak{e}^{\prime}$, and $f_{\mathfrak{e}}(\varepsilon)$ is the occupation
function for the electrode $\mathfrak{e}$. Each electrode has its own chemical
potential ($\mu_{\mathfrak{e}}$) and temperature ($T_{\mathfrak{e}}$). The
$\mathbf{k}$ dependence of all quantities here is implicit.

The transmission and reflection functions are directly related to the
scattering matrix
$\textbf{s}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)$:

$\displaystyle T_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)$
$\displaystyle=\mathrm{Tr}\left\\{\textbf{s}_{\mathfrak{e}\mathfrak{e}^{\prime}}^{\dagger}(\varepsilon)\textbf{s}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\right\\}$
$\displaystyle R_{\mathfrak{e}}(\varepsilon)$
$\displaystyle=M_{\mathfrak{e}}(\varepsilon)-\mathrm{Tr}\left\\{\textbf{s}_{\mathfrak{e}\mathfrak{e}}^{\dagger}(\varepsilon)\textbf{s}_{\mathfrak{e}\mathfrak{e}}(\varepsilon)\right\\}$
(43)

Furthermore, the unitarity of the scattering matrix implies
$M_{\mathfrak{e}}=R_{\mathfrak{e}}+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T_{\mathfrak{e}\mathfrak{e}^{\prime}}$
[66]. This relationship allows us to rewrite the current as a sum of pairwise
currents.

$\displaystyle I_{\mathfrak{e}}$
$\displaystyle=\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}I_{\mathfrak{e},\mathfrak{e}^{\prime}}$
$\displaystyle I_{\mathfrak{e},\mathfrak{e}^{\prime}}$
$\displaystyle=\int_{-\infty}^{\infty}\\!\\!d\varepsilon\,T_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\left[f_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)\right]$
(44)

In the low-bias limit the difference in the Fermi-occupation function of the
electrodes is proportional to the difference in their chemical potential and
the pairwise currents my be written as.

$\displaystyle I_{\mathfrak{e},\mathfrak{e}^{\prime}}$ $\displaystyle\approx
T_{\mathfrak{e}\mathfrak{e}^{\prime}}(\mu)(\mu_{\mathfrak{e}^{\prime}}-\mu_{\mathfrak{e}})$
(45)

In other words, the current is proportional to the applied bias.

In the case of superconducting electrodes or devices, the transmission and
reflection probabilities become 2$\times$2 matrices (Eq. 18) [48].

We obtain the reflection and transmission probabilities using Green’s function
techniques, adopting an approach similar to the generalized Fisher-Lee
relation [boumrar2020equivalence].

$\displaystyle
T_{\mathfrak{e}\mathfrak{e}^{\prime}}^{\alpha\beta}(\varepsilon)$
$\displaystyle=\mathrm{Tr}\left\\{\mathbf{s}^{\alpha\beta}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\mathbf{s}^{\alpha\beta}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)^{\dagger}\right\\}$
(46) $\displaystyle
R_{\mathfrak{e}\mathfrak{e}^{\prime}}^{\alpha\beta}(\varepsilon)$
$\displaystyle=\mathrm{Tr}\left\\{\mathbf{s}^{\alpha\beta}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\mathbf{s}^{\alpha\beta}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)^{\dagger}\right\\}$
(47) $\displaystyle s^{\alpha\beta}_{\mathfrak{e}\mathfrak{e}^{\prime}}$
$\displaystyle=-\delta_{\alpha\beta}\delta_{\mathfrak{ee}^{\prime}}M^{\alpha}_{\mathfrak{e}}(\varepsilon)+i\mathbf{\Gamma}_{\mathfrak{e}}^{1/2}(\varepsilon)\mathbf{G}(\varepsilon+i0^{+})\mathbf{\Gamma}_{\mathfrak{e}^{\prime}}^{1/2}(\varepsilon)$
(48) $\displaystyle\mathbf{G}(\varepsilon+i0^{+})$
$\displaystyle=(\varepsilon+i0^{+})\mathbf{S}-\mathbf{H}-\sum_{\mathfrak{e}}\mathbf{\Sigma}_{\mathfrak{e}}(\varepsilon)$
(49)

where
$\Gamma_{\mathfrak{e}}(\varepsilon)=i\left[\mathbf{\Sigma}_{\mathfrak{e}}(\varepsilon)-\mathbf{\Sigma}^{\dagger}_{\mathfrak{e}}(\varepsilon)\right]$
is the electrode broadening matrix,
$\mathbf{\Sigma}_{\mathfrak{e}}(\varepsilon)$ is the electrode self-energy,
$\mathbf{G}(\varepsilon+i0^{+})$ is the retarded Green’s function of the open
quantum system, and $M^{\alpha}_{\mathfrak{e}}(\varepsilon)$ is the number of
electrode channels.

The electrode self-energies are calculated using the Lopez-Sancho/Lopez-Sancho
algorithm [sancho1984quick]. The number of channels
$M^{\alpha}_{\mathfrak{e}}(\varepsilon)$ is determined in a separate
calculation of the bulk electrodes.

These additional scattering probabilities have to be taken into account when
calculating the current through any electrode. Taking into account the
opposite charge of electrons and holes, the current through the electrode
$\mathfrak{e}$ is given by:

$\displaystyle I_{\mathfrak{e}}=\frac{e}{h}\int_{-\infty}^{\infty}\Big{\\{}$
$\displaystyle\underbrace{\left[M_{\mathfrak{e}}(\varepsilon)-R^{pp}_{\mathfrak{e}}(\varepsilon)+R^{hp}_{\mathfrak{e}}(\varepsilon)\right]{f}_{\mathfrak{e}}(\varepsilon)-\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}\left[T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)-T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\right]{f}_{\mathfrak{e}^{\prime}}(\varepsilon)}_{\text{electronic
contribution}}$ $\displaystyle-$
$\displaystyle\underbrace{\left[M_{\mathfrak{e}}(\varepsilon)-R^{hh}_{\mathfrak{e}}(\varepsilon)+R^{ph}_{\mathfrak{e}}(\varepsilon)\right]{\bar{f}}_{\mathfrak{e}}(\varepsilon)-\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}\left[-T^{hh}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)+T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)]\right]{\bar{f}}_{\mathfrak{e}^{\prime}}(\varepsilon)}_{\text{hole
contribution}}\Big{\\}}d\varepsilon.$ (50) $\displaystyle
I_{\mathfrak{e}}=\frac{e}{h}\int_{-\infty}^{\infty}\Big{\\{}$
$\displaystyle\underbrace{\left[M_{\mathfrak{e}}(\varepsilon)-R^{pp}_{\mathfrak{e}}(\varepsilon)+R^{hp}_{\mathfrak{e}}(\varepsilon)\right]{f}_{\mathfrak{e}}(\varepsilon)-\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}\left[T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)-T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)\right]{f}_{\mathfrak{e}^{\prime}}(\varepsilon)}_{\text{electronic
contribution}}$ $\displaystyle-$
$\displaystyle\underbrace{\left[M_{\mathfrak{e}}(\varepsilon)-R^{hh}_{\mathfrak{e}}(\varepsilon)+R^{ph}_{\mathfrak{e}}(\varepsilon)\right]{\bar{f}}_{\mathfrak{e}}(\varepsilon)-\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}\left[-T^{hh}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)+T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)]\right]{\bar{f}}_{\mathfrak{e}^{\prime}}(\varepsilon)}_{\text{hole
contribution}}\Big{\\}}d\varepsilon.$ (51)

Here, $\bar{f}(\varepsilon)=1-f(\varepsilon)$ is the occupation function of
the hole-like states in the electrode. As in a normal metal, the unitarity of
the scattering matrix links the number of electrode channels to the reflection
and transmission functions [68]:

$\displaystyle M_{\mathfrak{e}}(\varepsilon)$
$\displaystyle=R^{pp}_{\mathfrak{e}}(\varepsilon)+R^{ph}_{\mathfrak{e}}(\varepsilon)+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}[T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)+T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)]$
$\displaystyle=R^{hh}_{\mathfrak{e}}(\varepsilon)+R^{hp}_{\mathfrak{e}}(\varepsilon)+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}[T^{hh}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)+T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)].$
(52)

Using these identities, we can rewrite the current through electrode $I_{e}$
(Eq. IV.1) as

$\displaystyle I_{\mathfrak{e}}=\frac{e}{h}\int_{-\infty}^{\infty}\Big{\\{}$
$\displaystyle[R^{hp}_{\mathfrak{e}}(\varepsilon)+R^{ph}_{\mathfrak{e}}(\varepsilon)][f_{\mathfrak{e}}(\varepsilon)-\bar{f}_{\mathfrak{e}}(\varepsilon)]$
$\displaystyle+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}[T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)+T^{hh}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)][f_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)]$
$\displaystyle+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)[\bar{f}_{\mathfrak{e}}(\varepsilon)-f_{\mathfrak{e}^{\prime}}(\varepsilon)]+T^{ph}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\varepsilon)[f_{\mathfrak{e}}(\varepsilon)-\bar{f}_{\mathfrak{e}^{\prime}}(\varepsilon)]\Big{\\}}d\varepsilon$
(53)

where $\bar{f}_{\mathfrak{e}}(\varepsilon)$ is the occupation of hole states
in the electrode $\mathfrak{e}$. Moving to the low-bias limit
($f_{\mathfrak{e}}(\varepsilon)=f(\varepsilon)-\frac{\partial{f}}{\partial{\varepsilon}}(\varepsilon)[\mu_{\mathfrak{e}}-\mu]$),
and invoking particle-hole symmetry allows us to rewrite the expression of the
current as (Eq. 20):

$\displaystyle I_{\mathfrak{e}}\approx$
$\displaystyle\frac{e}{h}\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{pp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\mu)[\mu_{\mathfrak{e}^{\prime}}-\mu_{\mathfrak{e}}]+\sum_{\mathfrak{e}^{\prime}\neq\mathfrak{e}}T^{hp}_{\mathfrak{e}\mathfrak{e}^{\prime}}(\mu)[\mu_{\mathfrak{e}^{\prime}}+\mu_{\mathfrak{e}}-2\mu]+2R^{hp}_{\mathfrak{e}}(\mu)(\mu_{\mathfrak{e}}-\mu)$
(54)

## References

  * [1] Walther Meissner and Robert Ochsenfeld.  Ein neuer effekt bei eintritt der supraleitfähigkeit.  Naturwissenschaften, 21(44):787–788, 1933.
  * [2] John Bardeen, Leon N Cooper, and John Robert Schrieffer.  Theory of superconductivity.  Physical review, 108(5):1175, 1957.
  * [3] John Bardeen, Leon N Cooper, and J Robert Schrieffer.  Microscopic theory of superconductivity.  Physical Review, 106(1):162, 1957.
  * [4] F. Steglich, J. Aarts, C. D. Bredl, W. Lieke, D. Meschede, W. Franz, and H. Schäfer.  Superconductivity in the Presence of Strong Pauli Paramagnetism: Ce Cu 2 Si 2.  Physical Review Letters, 43(25):1892–1896, December 1979.
  * [5] D. Jérome, A. Mazaud, M. Ribault, and K. Bechgaard.  Superconductivity in a synthetic organic conductor (TMTSF)2PF 6.  Journal de Physique Lettres, 41(4):95–98, 1980.
  * [6] C. C. Tsuei and J. R. Kirtley.  Pairing symmetry in cuprate superconductors.  Reviews of Modern Physics, 72(4):969–1016, October 2000.
  * [7] G. R. Stewart.  Superconductivity in iron compounds.  Reviews of Modern Physics, 83(4):1589–1652, December 2011.
  * [8] Leon Balents, Cory R Dean, Dmitri K Efetov, and Andrea F Young.  Superconductivity and strong correlations in moiré flat bands.  Nature Physics, 16(7):725–733, 2020.
  * [9] Benjamin Sacépé, Mikhail Feigel’man, and Teunis M Klapwijk.  Quantum breakdown of superconductivity in low-dimensional materials.  Nature Physics, 16(7):734–746, 2020.
  * [10] SM Frolov, MJ Manfra, and JD Sau.  Topological superconductivity in hybrid devices.  Nature Physics, 16(7):718–724, 2020.
  * [11] A Yu Kitaev.  Unpaired Majorana fermions in quantum wires.  Physics-Uspekhi, 44(10S):131–136, October 2001.
  * [12] Wen-Yu He, Benjamin T. Zhou, James J. He, Noah F. Q. Yuan, Ting Zhang, and K. T. Law.  Magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides.  Communications Physics, 1(1):40, July 2018.
  * [13] Joon Sue Lee, Borzoyeh Shojaei, Mihir Pendharkar, Anthony P. McFadden, Younghyun Kim, Henri J. Suominen, Morten Kjaergaard, Fabrizio Nichele, Hao Zhang, Charles M. Marcus, and Chris J. Palmstrøm.  Transport Studies of Epi-Al/InAs Two-Dimensional Electron Gas Systems for Required Building-Blocks in Topological Superconductor Networks.  Nano Letters, 19(5):3083–3090, May 2019.
  * [14] Elsa Prada, Pablo San-Jose, Michiel W. A. de Moor, Attila Geresdi, Eduardo J. H. Lee, Jelena Klinovaja, Daniel Loss, Jesper Nygård, Ramón Aguado, and Leo P. Kouwenhoven.  From Andreev to Majorana bound states in hybrid superconductor–semiconductor nanowires.  Nature Reviews Physics, 2(10):575–594, October 2020.
  * [15] Masatoshi Sato and Yoichi Ando.  Topological superconductors: a review.  Reports on Progress in Physics, 80(7):076501, 2017.
  * [16] L. N. Oliveira, E. K. U. Gross, and W. Kohn.  Density-Functional Theory for Superconductors.  Physical Review Letters, 60(23):2430–2433, June 1988.
  * [17] Pierre Hohenberg and Walter Kohn.  Inhomogeneous electron gas.  Physical review, 136(3B):B864, 1964.
  * [18] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Physical review, 140(4A):A1133, 1965.
  * [19] NN Bogoljubov, Vladimir Veniaminovic Tolmachov, and DV Širkov.  A new method in the theory of superconductivity.  Fortschritte der physik, 6(11-12):605–682, 1958.
  * [20] N. N. Bogoljubov.  On a new method in the theory of superconductivity.  Il Nuovo Cimento, 7(6):794–805, March 1958.
  * [21] P G DeGennes.  Superconductivity of Metals and Alloys.  (W. A. Benjamin, New York). 1966.
  * [22] W. M. Temmerman, Z. Szotek, B. L. Gyorffy, O. K. Andersen, and O. Jepsen.  Gap Anisotropy in the Layered High Temperature Superconductors.  Physical Review Letters, 76(2):307–310, January 1996.
  * [23] B. L. Gyorffy, Z. Szotek, W. M. Temmerman, O. K. Andersen, and O. Jepsen.  Quasiparticle spectra of high-temperature superconductors.  Physical Review B, 58(2):1025–1042, July 1998.
  * [24] P. Cudazzo, G. Profeta, A. Sanna, A. Floris, A. Continenza, S. Massidda, and E. K. U. Gross.  _Ab_ _Initio_ Description of High-Temperature Superconductivity in Dense Molecular Hydrogen.  Physical Review Letters, 100(25):257001, June 2008.
  * [25] Lilia Boeri and Giovanni B Bachelet.  Viewpoint: The road to room-temperature conventional superconductivity.  Journal of Physics: Condensed Matter, 31(23):234002, June 2019.
  * [26] Lilia Boeri.  Understanding Novel Superconductors with Ab Initio Calculations.  In Wanda Andreoni and Sidney Yip, editors, Handbook of Materials Modeling, pages 73–112. Springer International Publishing, Cham, 2020.
  * [27] M.B. Suvasini and B.L. Gyorffy.  A multiple scattering method for solving Bogoliubov-de Gennes equations of superconductivity.  Physica C: Superconductivity, 195(1-2):109–126, May 1992.
  * [28] M. B. Suvasini, W. M. Temmerman, and B. L. Gyorffy.  Computational aspects of density-functional theories of superconductors.  Phys. Rev. B, 48:1202–1210, Jul 1993.
  * [29] Gábor Csire, Balázs Újfalussy, József Cserti, and Balázs Győrffy.  Multiple scattering theory for superconducting heterostructures.  Physical Review B, 91(16):165142, 2015.
  * [30] Philipp Rüßmann and Stefan Blügel.  Density functional Bogoliubov-de Gennes analysis of superconducting Nb and Nb(110) surfaces.  Physical Review B, 105(12):125143, March 2022.
  * [31] Philipp Rüßmann and Stefan Blügel.  Proximity induced superconductivity in a topological insulator.  arXiv preprint arXiv:2208.14289, 2022.
  * [32] Yalin Zhang, Jingfeng Huang, Chunchen Zhang, Peng Wang, Zhihe Wang, Tong Wang, Zhongwen Xing, and DY Xing.  Proximity-effect-induced superconductivity in bi 2 te 3/fe se 0.5 te 0.5 thin-film heterostructures with different interface conditions.  Physical Review B, 102(6):064503, 2020.
  * [33] Cequn Li, Yi-Fan Zhao, Alexander Vera, Omri Lesser, Hemian Yi, Shalini Kumari, Zijie Yan, Chengye Dong, Timothy Bowen, Ke Wang, et al.  Proximity-induced superconductivity in epitaxial topological insulator/graphene/gallium heterostructures.  Nature Materials, 22(5):570–575, 2023.
  * [34] Bendegúz Nyári, András Lászlóffy, László Szunyogh, Gábor Csire, Kyungwha Park, and Balázs Ujfalussy.  Relativistic first-principles theory of yu-shiba-rusinov states applied to mn adatoms and mn dimers on nb (110).  Physical Review B, 104(23):235426, 2021.
  * [35] Kyungwha Park, Bendeguz Nyari, Andras Laszloffy, Laszlo Szunyogh, and Balazs Ujfalussy.  Effects of strong spin-orbit coupling on shiba states from magnetic adatoms using first-principles theory.  New Journal of Physics, 25(3):033022, 2023.
  * [36] András Lászlóffy, Bendegúz Nyári, László Szunyogh, Kyungwha Park, and Balázs Ujfalussy.  Classification of yu-shiba-rusinov states of magnetic impurities on superconducting surfaces: A fully relativistic first-principles study.  Physical Review B, 107(22):224515, 2023.
  * [37] Dimitar Pashov, Swagata Acharya, Walter RL Lambrecht, Jerome Jackson, Kirill D Belashchenko, Athanasios Chantis, Francois Jamet, and Mark van Schilfgaarde.  Questaal: A package of electronic structure methods based on the linear muffin-tin orbital technique.  Computer Physics Communications, 249:107065, 2020.
  * [38] Swagata Acharya, Dimitar Pashov, and Mark van Schilfgaarde.  Role of nematicity in controlling spin fluctuations and superconducting t c in bulk fese.  Physical Review B, 105(14):144507, 2022.
  * [39] Swagata Acharya, Mikhail I Katsnelson, and Mark van Schilfgaarde.  Vertex dominated superconductivity in intercalated fese.  npj Quantum Materials, 8(1):24, 2023.
  * [40] José M Soler, Emilio Artacho, Julian D Gale, Alberto García, Javier Junquera, Pablo Ordejón, and Daniel Sánchez-Portal.  The SIESTA method for ab initio order- N materials simulation.  Journal of Physics: Condensed Matter, 14(11):2745–2779, March 2002.
  * [41] Alberto Garcia, Nick Papior, Arsalan Akhtar, Emilio Artacho, Volker Blum, Emanuele Bosoni, Pedro Brandimarte, Mads Brandbyge, Jorge I Cerdá, Fabiano Corsetti, et al.  Siesta: Recent developments and applications.  The Journal of chemical physics, 152(20):204108, 2020.
  * [42] https://github.com/siesta-project and https://www.youtube.com/@siestaproject.
  * [43] R Cuadrado and JI Cerdá.  Fully relativistic pseudopotential formalism under an atomic orbital basis: spin–orbit splittings and magnetic anisotropies.  Journal of Physics: Condensed Matter, 24(8):086005, 2012.
  * [44] Mads Brandbyge, José-Luis Mozos, Pablo Ordejón, Jeremy Taylor, and Kurt Stokbro.  Density-functional method for nonequilibrium electron transport.  Physical Review B, 65(16):165401, March 2002.
  * [45] Nick Papior, Nicolás Lorente, Thomas Frederiksen, Alberto García, and Mads Brandbyge.  Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta.  Computer Physics Communications, 212:8–24, 2017.
  * [46] Hervé Ness, Ivan A Sadovskyy, Andrey E Antipov, Mark van Schilfgaarde, and Roman M Lutchyn.  Supercurrent decay in ballistic magnetic josephson junctions.  npj Computational Materials, 8(1):23, 2022.
  * [47] Herve Ness and Mark Van Schilfgaarde.  Ab initio transport calculations: from normal to superconducting current.  arXiv preprint arXiv:2309.06984, 2023.
  * [48] C J Lambert.  Generalized landauer formulae for quasi-particle transport in disordered superconductors.  Journal of Physics: Condensed Matter, 3(34):6579–6587, August 1991.
  * [49] F. Taddei, S. Sanvito, J. H. Jefferson, and C. J. Lambert.  Suppression of Giant Magnetoresistance by a Superconducting Contact.  Physical Review Letters, 82(24):4938–4941, June 1999.
  * [50] Gianluca Stefanucci, Enrico Perfetto, and Michele Cini.  Time-dependent quantum transport with superconducting leads: A discrete-basis Kohn-Sham formulation and propagation scheme.  Physical Review B, 81(11):115446, March 2010.
  * [51] Qisi Wang, Yao Shen, Bingying Pan, Xiaowen Zhang, K. Ikeuchi, K. Iida, A. D. Christianson, H. C. Walker, D. T. Adroja, M. Abdel-Hafiez, Xiaojia Chen, D. A. Chareev, A. N. Vasiliev, and Jun Zhao.  Magnetic ground state of FeSe.  Nature Communications, 7(1):12182, July 2016.
  * [52] Roberto Casalbuoni and Giuseppe Nardulli.  Inhomogeneous superconductivity in condensed matter and qcd.  Reviews of Modern Physics, 76(1):263, 2004.
  * [53] Rüßmann Philipp and Blügel Stefan.  Density functional Bogoliubov-de Gennes analysis of superconducting {Nb} and {Nb}(110) surfaces, October 2021.
  * [54] Lucas Schneider.  Realizing topological boundary modes in bottom-up constructed magnetic nanostructures on superconducting surfaces.  PhD thesis, Universität Hamburg, Hamburg, 2021.
  * [55] D. K. Finnemore, T. F. Stromberg, and C. A. Swenson.  Superconducting properties of high-purity niobium.  Phys. Rev., 149:231–243, Sep 1966.
  * [56] VR Karasik and I Yu Shebalin.  Superconducting properties of pure niobium.  SOV PHYS JETP, 30(6):1068–1075, 1970.
  * [57] Shigeru Kasahara, Tatsuya Watashige, Tetsuo Hanaguri, Yuhki Kohsaka, Takuya Yamashita, Yusuke Shimoyama, Yuta Mizukami, Ryota Endo, Hiroaki Ikeda, Kazushi Aoyama, Taichi Terashima, Shinya Uji, Thomas Wolf, Hilbert von Löhneysen, Takasada Shibauchi, and Yuji Matsuda.  Field-induced superconducting phase of FeSe in the BCS-BEC cross-over.  Proceedings of the National Academy of Sciences, 111(46):16309–16313, November 2014.
  * [58] J. K. Glasbrenner, I. I. Mazin, Harald O. Jeschke, P. J. Hirschfeld, R. M. Fernandes, and Roser Valentí.  Effect of magnetic frustration on nematicity and superconductivity in iron chalcogenides.  Nature Physics, 11(11):953–958, November 2015.
  * [59] Can-Li Song, Yi-Lin Wang, Peng Cheng, Ye-Ping Jiang, Wei Li, Tong Zhang, Zhi Li, Ke He, Lili Wang, Jin-Feng Jia, et al.  Direct observation of nodes and twofold symmetry in fese superconductor.  Science, 332(6036):1410–1413, 2011.
  * [60] Peter J Hirschfeld.  Using gap symmetry and structure to reveal the pairing mechanism in fe-based superconductors.  Comptes Rendus Physique, 17(1-2):197–231, 2016.
  * [61] Johnpierre Paglione and Richard L. Greene.  High-temperature superconductivity in iron-based materials.  Nature Physics, 6(9):645–658, September 2010.
  * [62] S-H. Baek, D. V. Efremov, J. M. Ok, J. S. Kim, Jeroen van den Brink, and B. Büchner.  Orbital-driven nematicity in FeSe.  Nature Materials, 14(2):210–214, February 2015.
  * [63] Pierre Massat, Donato Farina, Indranil Paul, Sandra Karlsson, Pierre Strobel, Pierre Toulemonde, Marie-Aude Méasson, Maximilien Cazayous, Alain Sacuto, Shigeru Kasahara, Takasada Shibauchi, Yuji Matsuda, and Yann Gallais.  Charge-induced nematicity in FeSe.  Proceedings of the National Academy of Sciences, 113(33):9177–9181, August 2016.
  * [64] P. Zhang, T. Qian, P. Richard, X. P. Wang, H. Miao, B. Q. Lv, B. B. Fu, T. Wolf, C. Meingast, X. X. Wu, Z. Q. Wang, J. P. Hu, and H. Ding.  Observation of two distinct d x z / d y z band splittings in FeSe.  Physical Review B, 91(21):214503, June 2015.
  * [65] H. Ness and M. van Schilfgaarde.  Ab initio transport calculations: From normal to superconducting current.  2023\.
  * [66] Supriyo Datta.  Electronic Transport in Mesoscopic Systems.  Cambridge University Press, 1 edition, September 1995.
  * [67] J Ferrer, C J Lambert, V M García-Suárez, D Zs Manrique, D Visontai, L Oroszlany, R Rodríguez-Ferradás, I Grace, S W D Bailey, K Gillemot, Hatef Sadeghi, and L A Algharagholy.  GOLLUM: A next-generation simulation tool for electron, thermal and spin transport.  New Journal of Physics, 16(9):093029, September 2014.
  * [68] C J Lambert and R Raimondi.  Phase-coherent transport in hybrid superconducting nanostructures.  Journal of Physics: Condensed Matter, 10(5):901–941, February 1998.
  * [69] Studiosus Kirchhoff.  Ueber den Durchgang eines elektrischen Stromes durch eine Ebene, insbesondere durch eine kreisförmige.  Annalen der Physik, 140(4):497–514, January 1845.
  * [70] Supriyo Datta, Philip F Bagwell, and Manjeri P Anantram.  Scattering theory of transport for mesoscopic superconductors.  1996\.
  * [71] Z. K. Tang, Lingyun Zhang, N. Wang, X. X. Zhang, G. H. Wen, G. D. Li, J. N. Wang, C. T. Chan, and Ping Sheng.  Superconductivity in 4 Angstrom Single-Walled Carbon Nanotubes.  Science, 292(5526):2462–2465, June 2001.
  * [72] I. Takesue, J. Haruyama, N. Kobayashi, S. Chiashi, S. Maruyama, T. Sugai, and H. Shinohara.  Superconductivity in Entirely End-Bonded Multiwalled Carbon Nanotubes.  Physical Review Letters, 96(5):057001, February 2006.
  * [73] M. Kociak, A. Yu. Kasumov, S. Guéron, B. Reulet, I. I. Khodos, Yu. B. Gorbatov, V. T. Volkov, L. Vaccarini, and H. Bouchiat.  Superconductivity in Ropes of Single-Walled Carbon Nanotubes.  Physical Review Letters, 86(11):2416–2419, March 2001.
  * [74] Ken-ichi Sasaki, Jie Jiang, Riichiro Saito, Seiichiro Onari, and Yukio Tanaka.  Theory of Superconductivity of Carbon Nanotubes and Graphene.  Journal of the Physical Society of Japan, 76(3):033702, March 2007.
  * [75] G. E. Blonder, M. Tinkham, and T. M. Klapwijk.  Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion.  Physical Review B, 25(7):4515–4532, April 1982.
  * [76] Yann in’t Veld, Mikhail I Katsnelson, Andrew J Millis, and Malte Rösner.  Screening induced crossover between phonon-and plasmon-mediated pairing in layered superconductors.  2D Materials, 2023.
  * [77] M Simonato, MI Katsnelson, and M Rösner.  Revised tolmachev-morel-anderson pseudopotential for layered conventional superconductors with nonlocal coulomb interaction.  Physical Review B, 108(6):064513, 2023.
  * [78] Gianluca Stefanucci, Enrico Perfetto, and Michele Cini.  Time-dependent quantum transport with superconducting leads: A discrete-basis kohn-sham formulation and propagation scheme.  Physical Review B, 81(11):115446, 2010.
  * [79] Swagata Acharya, Dimitar Pashov, Francois Jamet, and Mark van Schilfgaarde.  Electronic origin of t c in bulk and monolayer fese.  Symmetry, 13(2):169, 2021.
  * [80] We interchangeably use the term electrons (holes) and particles (anti-particles).
  * [81] R Balian and NR Werthamer.  Superconductivity with pairs in a relative p wave.  Physical review, 131(4):1553, 1963.
  * [82] Peter E. Blöchl, O. Jepsen, and O. K. Andersen.  Improved tetrahedron method for Brillouin-zone integrations.  Physical Review B, 49(23):16223–16233, June 1994.  Publisher: American Physical Society.
  * [83] B. L. Blackford and R. H. March.  Tunneling Investigation of Energy-Gap Anisotropy in Superconducting Bulk Pb.  Physical Review, 186(2):397–399, October 1969.
  * [84] G. I. Lykken, A. L. Geiger, K. S. Dy, and E. N. Mitchell.  Measurement of the Superconducting Energy Gap and Fermi Velocity in Single-Crystal Lead Films by Electron Tunneling.  Physical Review B, 4(5):1523–1530, September 1971.
  * [85] R. Heid, K.-P. Bohnen, I. Yu Sklyadneva, and E. V. Chulkov.  Effect of spin-orbit coupling on the electron-phonon interaction of the superconductors Pb and Tl.  Physical Review B, 81(17):174527, May 2010.
  * [86] I. Yu Sklyadneva, R. Heid, P. M. Echenique, K.-B. Bohnen, and E. V. Chulkov.  Electron-phonon interaction in bulk Pb: Beyond the Fermi surface.  Physical Review B, 85(15):155115, April 2012.
  * [87] Anshu Sirohi, Preetha Saha, Sirshendu Gayen, Avtar Singh, and Goutam Sheet.  Transport spectroscopy on trapped superconducting nano-islands of Pb: Signature of unconventional pairing.  Nanotechnology, 27(28):285701, July 2016.
  * [88] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized Gradient Approximation Made Simple.  Physical Review Letters, 77(18):3865–3868, October 1996.
  * [89] M.J. van Setten, M. Giantomassi, E. Bousquet, M.J. Verstraete, D.R. Hamann, X. Gonze, and G.-M. Rignanese.  The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table.  Computer Physics Communications, 226:39–54, May 2018.
  * [90] J. P. Perdew and Alex Zunger.  Self-interaction correction to density-functional approximations for many-electron systems.  Physical Review B, 23(10):5048–f5079, May 1981.