Superconductivity in Nb: Impact of Temperature, Dimensionality and Cooper-Pairing

The ability to realistically simulate the electronic structure of superconducting materials is important to understand and predict various properties emerging in both the superconducting topological and spintronics realms. We introduce a tight-binding implementation of the Bogoliubov–de Gennes method, parameterized from density functional theory, which we utilize to explore the bulk and thin films of Nb, known to host a significant superconducting gap. The latter is useful for various applications such as the exploration of trivial and topological in-gap states. Here, we focus on the simulation’s aspects of superconductivity and study the impact of temperature, Cooper-pair coupling and dimensionality on the value of the superconducting pairing interactions and gaps.


I. INTRODUCTION
Superconductivity is a basic fundamental phenomenon in solid-state physics.Although being discovered more than a century ago, it still generates a plethora of research activities.
Recently, its importance flourished again in the context of topological superconductivity [1][2][3], which is an essential field of research for the realization of quantum computing.Certainly, the intricate interplay of standard spintronics and quantum information concepts can give rise to unforeseen breakthroughs in basic research and in future technologies for storing, computing and transmitting information.This calls for the development of theoretical frameworks based on a realistic description of the electronic structure of superconducting materials interfaced with magnetic systems in order to understand and predict various related phenomena.
Interesting topics emerged in superconducting spintronics [4][5][6][7], where the interplay of Cooper pairs and other electronic degrees of freedom (such as spin, charge and orbital) give rise to fascinating phenomena, ranging from supercurrents, triplet superconductivity, and torques.Simultaneously, the potential design of topological superconductivity and the creation of Majorana states [8][9][10] is driving a lot of research in the context of topological qubits [11].
The microscopic theory to describe conventional superconductivity goes back to the work of John Bardeen, Leon Cooper, and Robert Schrieffer (BCS) [12].In this context, the Bogoliubov-de Gennes (BdG) method [13][14][15][16] is an elegant mean-field approximation that relies upon Bogoliubov-Valatin transformations that take the Hamiltonian from a particle space into a particle-hole one-a framework that has been frequently used, see for example Refs.17-27.Nowadays, several methods based on a realistic description of the electronic structure of materials are capable of treating fundamental aspects related to superconductivity.For instance, powerful theoretical frameworks have been proposed and utilized within density functional theory (DFT) combined with the BdG formalism [28][29][30][31][32][33], or by proposing an extension of DFT to account for superconductivity [34][35][36][37].Tight-binding schemes are very appealing since they allow a versatile description of superconductivity-related physics with a potential control of various levels of approximations, which could enable the treatment of problems not attainable with conventional DFT.
The goal of this article is to introduce our recently developed method that allows the treatment of superconductivity on the basis of a tight-binding approach as implemented in TITAN [49][50][51][52].The framework enables a realistic description of the electronic structure, with parameters obtained from DFT [53].Examples of applications of TITAN ranges from dynamical magnetic responses [54], dynamical transport and torques [50,52], magnetic damping [51], as well as ultrafast magnetization dynamics triggered by laser pulses [55,56].
In the context of superconductivity, the method was already introduced in Ref. 57, with a focus on magnetic exchange interactions at superconducting-magnetic interfaces.Here, we pay attention to Nb in bulk and thin films with the aim of highlighting fundamental aspects in the simulations procedure.Aspects related to the choice of the electron-phonon coupling to open the appropriate band gap are discussed together with the impact of temperature and self-consistency on the underlying physics.Dimensionality affects superconductivity, as we unveil, which can be utilized to tune the superconducting gap for different applications in modern investigations of superconductivity-related problems.
The article is organized as follows.We first introduce the method in Sec.II, which is based on multi-orbital tight-binding theory.Then we present in Sec.III the results of our simulations in Nb bulk, (001) and (110) monolayers, which are followed by Nb(110) thin films of different thicknesses.Finally, we present our concluding remarks in Sec.IV.

II. TIGHT-BINDING AND THE BOGOLIUBOV-DE GENNES METHOD
To simulate the superconducting properties of Nb across dimensions, we introduce here the tight-binding methodology together with the BdG method that we utilized, as implemented in our code TITAN [49,51,52,[55][56][57].The hopping parameters were obtained from first-principles calculations [53].More details on superconducting related aspects are elaborated in Ref. 57.
The fundamental tight-binding Hamiltonian reads: where c † αµσ (k) and c βνη (k) consist of the creation and annihilation operators associated to electrons having a wave vector k and spin σ in the orbitals µ and spin η in the orbitals ν , respectively.α and β are specific atoms in the bulk unit cell, or layers in the film geometry.k is a reciprocal vector, while N is the number of wave vectors in the three-dimensional (for bulk) or two-dimensional (for films) Brillouin zone.To enable the formation of Cooper pairs, which ultimately gives rise to superconductivity, we introduce λ αµ in the second term of the previous equation, as induced by the electron-phonon coupling.
In the mean-field approximation, Eq. (1) simplifies to with ∆ αµ is known as the superconducting gap parameter, which corresponds to the energy required to scatter Cooper pairs [58].
The term H µν αβ,ση (k) corresponds to the regular Hamiltonian where superconductivity is not accounted for: with H 0µν αβ being the spin-independent tight-binding term, the second term is responsible for the intra-atomic exchange interaction (originating from a Hubbard-like contribution [50,52]), which can lead to magnetism, and the third term corresponds to the spin-orbit interaction.
The BdG methodology consists in utilizing the Bogoliubov-Valatin transformation [16] to rewrite Eq. ( 2) in the electron-hole representation instead of the electron one, which leads to the BdG Hamiltonian: with E F being the Fermi energy while the associated eigenvalue problem to solve reads: The transformation is canonical and enabled by rewriting with where the prime indicates that the sums run only over the states with positive energy [16,59].
This restriction in the sum is done to counteract the doubling of the degrees of freedom originated from the change of basis.
The eigenvector of Eq. ( 6) is Our numerical procedure to evaluate the various properties of a given material is based on the self-consistency of the charge density, which is also pursued once the electron-phonon coupling is included.In practice, we start the simulations considering the normal metallic phase then we incorporate the Cooper pairing and repeat the computational runs until convergence.

III. RESULTS
In this section we explore the results of simulations done for various Nb structures.Our choice on Nb is motivated by its relatively large superconducting gap, which facilitates the exploration of in-gap states.Simulating low temperatures phenomena requires finer computational meshes, thus for the first investigations we choose high temperatures and large values for the superconducting coupling λ.We start with bulk Nb, where we recover a superconducting gap opening in the electronic structure around the Fermi energy.Then we explore the case of thin films down to the single monolayer along the [001] or [110] directions.
We monitor how the superconducting gap is affected by the self-consistency of the electronic structure calculations.

A. Bulk superconducting Nb
Nb has a superconducting critical temperature of T C = 9.3 K [60], which leads to a significant superconducting gap 2∆ = 3.8 meV.A rough approximation relates T c with the size of the gap-at T=0 K-through [58] ∆ Large values of ∆ are desirable when investigating states that appear inside the gap.In This holds true for the DOS as well.We expect that by allowing for finite coupling λ, degeneracies occurring at crossing points will be lifted.
In subsection III B, we describe a rather optimized search method to tune this parameter.
From our explorations, we found that: (i) large values of λ do not always mean large gaps and (ii) the size of the gap is not a linear function of λ.From tabulated values of the electron-phonon λ e-p coupling constant we can evaluate the coupling parameter λ utilizing the following relation [62] where D(E F ) is the density of states at the Fermi energy.For example, for bulk Nb with the data from Refs.63 we obtain λ ≈ 1.3 eV.In Ref. 62, the electron-phonon coupling constant is calculated for slabs of Nb(100) and Nb(110) of 3, 6, 9, 12, and 15 layers.For the Nb(110), the results for the system of 9 layers-the slab size closer to the one we will study laterlead to coupling parameters of λ ≈ 5.28 eV at the surfaces and λ ≈ 1.3 eV in the middle layer.When plotting the curve associated to ∆ vs. λ, we obtain a quadratic function, in agreement with Ref. 64, where the self-consistent Korringa-Kohn-Rostoker method is used.
There, the authors found that in order to obtain ∆ ≈ 1 meV the coupling parameter must be λ ≈ 1.1 eV for bulk Nb.
As a first example of our investigations, we set λ µ = 5.44 eV.The occupancies are ⟨n s ⟩ = 0.70317397 , ⟨n p ⟩ = 0.66860504 , ⟨n d ⟩ = 3.62822099, which experience a minor change of the order of 10 −2 with respect to the case where λ is set to zero.The self-consistent value of the superconducting gap parameter is ∆ = 0.497 eV, which is obviously gigantic.The resulting band structure and LDOS plots are illustrated in Fig. 1(c).The red dotted lines delimit the range E ∈ [−∆, ∆].As a sanity check, it is reassuring that the gap in the band structure matches the self-consistent output value of ∆.The density of states does not show a perfect gap, but decays at the gap region.This effect is due to the thermal broadening induced artificially, which is in this case set to 43 meV.
In subsections III B and III D, we will systematically explore the impact of temperature and coupling strengths, searching for an optimal λ that yields the experimental value of ∆.
B. Superconducting gap tuning: single bcc (001) Nb monolayer Explorations to find the appropriate coupling λ that results in the desired gap parameter ∆ involve systematic repeated simulations with different parameters.In this section, we show the typical curves for ∆ as a function of the temperature (via the broadening parameter) and the Cooper pair coupling parameter.To exemplify, we use a 2D system, namely a (001) monolayer of bcc Nb.As in our previous examples, we set the same λ for all orbitals, with eight values between 4 and 9 eV.We consider forty different temperature points between 300 K and 14 000 K. Simulations at lower temperatures are computationally highly expensive since they require an increased number of k-points to be able to capture the narrow peaky structures in the Brillouin zone with lower broadening.All the simulations in this section use the same number of k-points, ≈ 2 × 10 6 .
Fig. 2(a) displays a heatmap with the results of the gap parameter for the different temperatures and couplings.We can see that the gap parameter drops to zero after reaching the corresponding critical temperature for every line.Directing our attention for fixed temperatures (i.e. the "columns" in the figure), it is clear that there is also a critical λ c below which ∆ is zero.Throughout the heatmap, ∆ ranges from 0.0 to 1.76 eV.We can analyze the same data from a different perspective, as shown in Fig. 2(b).There we see that the ∆ vs. T curve has indeed the shape predicted by the BCS theory.For sufficiently low temperatures, ∆ approaches a plateau.Around the critical temperature T c , the selfconsistency procedure may become unstable and difficult to achieve a precise value, which lead some curves to miss points in that range.We can proceed to an extrapolation to evaluate the missing points on the basis of the asymptotic behavior of ∆ around the critical temperature [58] ∆ We used the precedent equation to go from the data in Fig. 2(b) to the one in Fig. 2(a).The BCS theory predicts that, at zero Kelvin, the relation between ∆ and T c is the one in Eq. (10).
Experimentally the reported values range from 2∆ = 3.2k b T c to 2∆ = 4.6k b T c , we signaled this range with the colored bands in Fig. 2(b).This plot establishes the aforementioned statement regarding the non-linearity of ∆ with respect to λ.
For a wide range of simulations, it is desirable and recommendable to place the superconducting system well below T c , at a point in the plateau of the ∆ vs. T curve.The superconducting state is frailer at the vicinity of the critical temperature, and internal changes in the material can perturb the properties of interest.
C. Renormalization after self-consistency: single bcc (001) Nb monolayer A common problem in simulating superconductors interfaced with normal conductors is the artificially increased size of the gap in order to facilitate the exploration of emerging in-gap states and it is rather common to proceed to single-iteration simulations, which assumes the electronic structure is not significantly altered in comparison to the case, where instead one proceeds to self-consistency.The importance of self-consistency was highlighted in Ref. 65.In this subsection, we address this aspect of simulation of superconducting materials by comparing the size of the gap obtained after one iteration (starting from the converged electronic structure of the metallic phase) with respect to the one obtained after full self-consistency.In Fig. 2(c), we single out the value of the superconducting gap parameter through several simulations with different coupling parameters at a temperature of 43 meV.Our algorithm, based on the Powell method [66] can halt the simulations when no improvement is detected in the last five steps.All shown cases converged without a problem.
In Fig. 2(d) we portray the percentage change of the gap parameter (( between the first-and last-iteration calculated ∆ for different λ.We note that the change is bigger for smaller values of λ.
So far, we have been discussing the numerics of the simulations assuming both large superconducting gaps and large electronic temperatures.In section III D, we drop T to realistic values associated to a superconducting Nb layer.Additionally, we explore the impact on the orbital-dependent gaps ∆ µ .
the evolution of the superconducting gap parameter across the material.

E. Superconducting slab
It is difficult to predict if the coupling strength λ that works well for one given system still opens a reasonable gap for a different one.We started with the value found in the previous section, λ = 1.37 eV, which opened a gap ∆ = 1.4 meV for a single monolayer.It turns out, however, that as soon as the Nb film thickness is larger than 5 layers, a large coupling is required to open a gap.In the remained study, we assume λ = 2.72 eV at a temperature of 4.56 K.
In Fig. 4(a), we show the magnitude of the gap in the central layer of the films as a function of the slab thickness.The orange line corresponds to the gap parameter of the bulk Nb for λ = 2.72 eV at the same temperature.One notices that the thicker the film gets, the closer the gap is to the one obtained for the bulk, a behavior comparable to what has been reported in Ref. 67.
The superconducting gap is not homogeneous across the Nb films.An example is shown in Fig. 4(b) for a film with a thickness of 9 layers.The gap in the middle of the film is larger than at the surface and it seems to oscillate as a function of position, as we would expect for wave-like particles in a confined system.A lower gap at the surface is expected since the electrons there have fewer opportunities to couple in Cooper-pairs due to their lack of neighbors from one side.We must also consider surface effects; for example, we have around 1.5 electrons less than the inner parts.Electrons move from the surfaces to the middle layers.Interestingly some of the layers in the film can have a larger gap than the bulk one.This indicates that the dimensionality of a material, can be utilized to control the superconducting gap due to confinement effects.

IV. CONCLUSIONS
After briefly introducing the tight-binding method that we developed to tackle superconductivity within the Bogoliubov-de Gennes formalism, we presented the results of our simulations on a material that is currently heavily utilized to explore the physics of superconductivity.We addressed the case of bulk and thin films of Nb along different directions.
A major issue in this field from the computational point of view is to choose the right Cooper pair coupling that opens the desired superconducting gap since the calculations are extremely expensive to resolve gaps of a few meV or even smaller.We analysed how the superconducting order parameter correlates with both the Cooper pair coupling and temperature, and explored the impact of self-consistency on the emerging orbital-dependent superconducting gap.In the thin film geometry, we unveiled the possibility of engineering the magnitude of the gap via confinement effects.Indeed, the superconducting gap is found to be layer-dependent, either smaller or larger than the values found in the bulk phase.
our simulation procedure, we start with bulk Nb at high temperatures.In TITAN, the thermal effects enter through the parameter η = k b T /π, which broadens the Fermi-Dirac function; in this case, we will use k b T = 43 meV (η = 1 mRy).For a reference, this is higher than room temperature, which is at 26 meV, but nevertheless, considerably lower than the melting point of Nb 237 meV.Nb has a body-centered cubic crystal structure.We make use of the Slater-Koster (SK) parameters from Ref.53, the spin-orbit coupling strength comes from Ref. 61.In Fig. 1(a), we show the density of states of both spin channels (left and right panels), and in the middle the band structure of bulk Nb.The charge occupancy per orbital is ⟨n s ⟩ = 0.70 , ⟨n p ⟩ = 0.68 , ⟨n d ⟩ = 3.61.The density of states (DOS) and the band structure agree with the results displayed in Ref. 53.Repeating the calculations by solving the BdG equations but setting λ = 0 on all orbitals, we naturally recover the metallic band-structure shown in Fig.1(b), where both the electron and hole parts are illustrated.Both parts are symmetric with respect to the Fermi level.

FIG. 1 .
FIG. 1. Electronic structure of bulk Nb with and without superconductivity.Three cases are illustrated: a without the BdG formalism only the electron part is shown; b with the BdG formalism but setting the Cooper pair coupling λ to zero; c with BdG and λ = 5.44 eV in all orbitals.The red colored line evinces the location of the gap as predicted by the self-consistent value ∆ = 0.497 eV and thermal broadening η = 43 meV.The band structure is shown in the middle.To the right (left), the DOS of the majority (minority) spin channel is presented.The legend box depicts the colors used to signal the total, the band resolved DOS and their electron versus hole parts, i.e. the contributions from respectively the u and v components of the eigenvectors.

FIG. 2 .
FIG. 2. Dependence of the superconducting gap of Nb(001) monolayer on various parameters.a Heat map of the superconducting gap as function of the coupling λ and the thermal broadening parameter.All curves drop to zero after reaching their corresponding critical temperature.b Another version of a in order to better follow how the superconducting gap evolves against temperature.The bands in color reflect an rough extrapolation to the critical temperature.c Evolution of the superconducting gap across the self-consistency for several coupling strengths λ. d Relative change of the superconducting gap when comparing the values obtained after one single iteration to those after self-consistency ((∆ initial − ∆ final )/∆ initial ) for different values of λ.

FIG. 3 .FIG. 4 .
FIG. 3. Orbital-resolved superconducting gap for Nb(110) monolayer.a ∆ µ grows considerably for all the orbital except p z .b Two close ups to the region where the material transitions from the non-superconducting to the superconducting phase.Simulations performed at 4.56 K.