Topological superconductivity from magnetic impurities on monolayer NbSe$_2$

TL;DR

This paper explores how magnetic impurity chains on monolayer NbSe2 can host topological superconducting phases with Majorana bound states, influenced by spin-orbit coupling and impurity configurations.

Contribution

It demonstrates the conditions for topological phases and Majorana states in magnetic impurity chains on monolayer NbSe2, highlighting the role of Ising spin-orbit coupling.

Findings

Topological phases are stabilized by strong Ising spin-orbit coupling.

Majorana bound states appear at the edges of impurity chains.

Conditions depend on impurity spacing, spin projection, and magnetic exchange strength.

Abstract

Recent experimental studies have found that magnetic impurities deposited on superconducting monolayer NbSe$_2$ generate coupled Yu-Shiba-Rusinov bound states. Here we consider ferromagnetic chains of impurities which induce a Yu-Shiba-Rusinov band and harbor Majorana bound states at the chain edges. We show that these topological phases are stabilized by strong Ising spin-orbit coupling in the monolayer and examine the conditions under which Majorana phases appear as a function of distance between impurities, impurity spin projection, orientation of chains on the surface of the monolayer, and strength of magnetic exchange energy between impurity and superconductor.

Tables

Table 1. Table 1: Parameters (in eV eV \mathrm{e}\mathrm{V} ) of the tight-binding Hamiltonian ( 2 ) from a fitting of ab initio data.

$\mu$	$t_1$	$t_2$	$t_3$	$t_4$	$t_5$	${\lambda }_1$	${\lambda }_2$
0.023	0.0134	0.097	0.0066	-0.0102	-0.0144	0.0163	0.0013

Table 2. Table 2: Parameters (in eV eV \mathrm{e}\mathrm{V} ) of the three-band tight-binding Hamiltonian ( 40 ) from a fitting of ab initio data.

$t_0$	$t_1$	$t_2$	$t_{11}$	$t_{12}$	$t_{22}$	${\epsilon }_1$	${\epsilon }_2$	$\lambda$
$u_0$	$u_1$	$u_2$	$u_{11}$	$u_{12}$	$u_{22}$	$r_0$	$r_1$	$r_2$	$r_{11}$	$r_{12}$
$-0.2399$	0.3245	0.3389	0.2608	0.2918	$-0.0511$	$1.5047$	1.9518	0.1372
0.0686	$-0.0421$	0.0105	0.0919	$-0.0254$	$-0.1127$	0.0098	$-0.0727$	0.0101	0.0365	0.0871

Equations

$${\mathcal{H}}_0=\frac{1}{2}\sum _{\mathbf{k}}{C}_{\mathbf{k}}^{†}{H}_0(\mathbf{k})C_{\mathbf{k}},$$

$$H_0(\mathbf{k})={\xi }_k{\tau }_z+{\Lambda }_{\mathbf{k}}{\tau }_z{\sigma }_z+\Delta {\tau }_x,$$

$$H^{\eta K}(\mathbf{q})$$

$$H^{\Gamma }(\mathbf{q})=\frac{{\hslash }^2{q}^2}{2m_{\Gamma }}-{\mu }_{\Gamma }+{\gamma }_3{q}^3\cos (3\chi ){\sigma }_z+\mathcal{O}(q^4),$$

$$G_{ij}(E)=\frac{1}{\Omega }\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot {\mathbf{r}}_{ij}}\frac{1}{E-H_0(\mathbf{k})+i{0}^{+}},$$

$$G_{ij}(E)\simeq G_{ij}^K(E)+G_{ij}^{\Gamma }(E),$$

$$G_{ij}^K(E)=\frac{1}{3\Omega }\sum _{n=0}^2\sum _{\eta ,\lambda =\pm }\sum _{\mathbf{q}}{e}^{i(\eta {\mathbf{K}}_n+\mathbf{q})\cdot {\mathbf{r}}_{ij}}{G}^{\eta \lambda }(E,\mathbf{q}),$$

$$G^{\eta \lambda }(E,\mathbf{q})=\frac{1}{E-{\xi }_{q,\eta \lambda }{\tau }_z-\Delta {\tau }_x}\frac{1+\lambda {\sigma }_z}{2}.$$

$$G_{ij}^{\Gamma }=\frac{1}{\Omega }\sum _{\mathbf{q}}{e}^{i\mathbf{q}\cdot {\mathbf{r}}_{ij}}\frac{1}{E-{\xi }_q^{\Gamma }{\tau }_z-\Delta {\tau }_x},$$

$$G_{ij}^K(E)$$

$$G_{ij}^{\Gamma }(E)$$

$$G_{ii}(E)\simeq -\frac{\pi \rho (0)}{2}\frac{E+\Delta {\tau }_x}{\sqrt{{\Delta }^2-E^2}},$$

$$H_1(\mathbf{r})=-JS\sum _{j=1}^N\mathbf{n}\cdot \mathbf{\sigma }\delta (\mathbf{r}-{\mathbf{r}}_j),$$

$$[H_0(\mathbf{r})+H_1(\mathbf{r})]{\psi }_{\mathbf{r}}=E{\psi }_{\mathbf{r}}.$$

$$H_0(\mathbf{q})\psi (\mathbf{q})-\frac{JS\mathbf{n}\cdot \mathbf{\sigma }}{\sqrt{\Omega }}\sum _j{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_j}{\psi }_j=E\psi (\mathbf{q}).$$

$$\psi (\mathbf{q})=-\frac{JS}{\sqrt{\Omega }}\sum _{j}G(E,\mathbf{q})e^{-i\mathbf{q}\cdot {\mathbf{r}}_j}\mathbf{n}\cdot \mathbf{\sigma }{\psi }_j.$$

$$G_{ij}(E)=\frac{1}{\Omega }\sum _{\mathbf{q}}G(E,\mathbf{q})e^{i\mathbf{q}\cdot {\mathbf{r}}_{ij}},$$

$${\psi }_i=-JS\sum _j{G}_{ij}(E)\mathbf{n}\cdot \mathbf{\sigma }{\psi }_j.$$

$$G_{ij}(E)\simeq G_{ij}(0)+E{\delta }_{ij}{\partial }_E{G}_{ij}(0),$$

$$\sum _{j=1}^N[{\delta }_{ij}+G_{ij}(0)JS\mathbf{n}\cdot \mathbf{\sigma }]{\psi }_j=-E{\partial }_E{G}_{ii}(0)JS\mathbf{n}\cdot \mathbf{\sigma }{\psi }_i.$$

$$\sum _{j=1}^N{H}_{ij}{\psi }_j=E{\psi }_i,$$

$$H_{ij}$$

$$E\psi =\frac{\Delta }{\alpha }(\mathbf{n}\cdot \mathbf{\sigma }-\alpha {\tau }_x)\psi .$$

$$\frac{E}{\Delta }=\pm \frac{1-{\alpha }^2}{1+{\alpha }^2}.$$

$${\psi }_{+}=\frac{1}{2}\left(\begin{array}{c}1\\ 1\end{array}\right)\otimes \left(\begin{array}{c}\cos \frac{\theta }{2}\\ \sin \frac{\theta }{2}{e}^{i\phi }\end{array}\right),{\psi }_{-}=\frac{1}{2}\left(\begin{array}{c}1\\ -1\end{array}\right)\otimes \left(\begin{array}{c}\sin \frac{\theta }{2}{e}^{-i\phi }\\ -\cos \frac{\theta }{2}\end{array}\right).$$

$${\bar{H}}_{ij}=\left(\begin{array}{cc}{h}_{ij}+b_{ij}& {\bar{\Delta }}_{ij}\\ -{\bar{\Delta }}_{ij}^{\ast }& -h_{ij}+b_{ij}\end{array}\right),$$

$$\frac{h_{ij}}{\Delta }$$

$$\frac{b_{ij}}{\Delta }$$

$$\frac{{\bar{\Delta }}_{ij}}{\Delta }$$

$$I_{0,1}(0,|{\mathbf{r}}_{ij}|)\sim \frac{1}{\sqrt{|{\mathbf{r}}_{ij}|}}{e}^{-|{\mathbf{r}}_{ij}|/{\xi }_0},$$

$$\bar{H}(k)=\sum _j{e}^{-ikjd}{\bar{H}}_j,$$

$$\bar{H}(k)=\left(\begin{array}{cc}h(k)+b(k)& \bar{\Delta }(k)\\ {\bar{\Delta }}^{\ast }(k)& -h(k)+b(k)\end{array}\right)$$

Full text

I Introduction

Topological superconductivity in one dimension is characterized by the presence of zero-energy Majorana bound states (MBS) which encode nonlocally electronic degrees of freedom. Kitaev (2001) This property, together with the non-Abelian statistics of their exchange, holds great promise to realize fault-tolerant topological quantum computation, where the coherence of qubits built from MBS is protected against local perturbations. Kitaev (2006); Nayak et al. (2008); Alicea (2012); Leijnse and Flensberg (2012); Beenakker (2013); Sarma et al. (2015); Elliott and Franz (2015); Sato and Fujimoto (2016); Aguado (2017); Lutchyn et al. (2018)

The first experimental signatures of MBS were recorded in conductance measurements of semiconducting nanowires proximitized by $s$-wave superconductors Mourik et al. (2012); Lutchyn et al. (2010); Oreg et al. (2010). Nevertheless, studies have proposed a wealth of additional platforms for the realization of one-dimensional $p$-wave topological superconductivity with MBS, including topological insulator-superconductor hybrids, Fu and Kane (2009); Cook and Franz (2011) cold atom quantum wires Jiang et al. (2011), proximitized half-metallic ferromagnets Duckheim and Brouwer (2011), nanomagnets on superconductors, Kjaergaard et al. (2012) etc. The present paper focuses on proposals which generate MBS through the interplay between magnetic impurities and superconductivity.

The magnetic impurities break time-reversal and spin-rotation symmetries in $s$-wave superconductors and generate localized in-gap bound states, called Yu-Shiba-Rusinov (YSR) states Yu (1965); Shiba (1968); Rusinov (1969); Balatsky et al. (2006), first detected experimentally for \ceMn and \ceGd adatoms on superconducting \ceNb. Yazdani et al. (1997) When the impurities are close enough, the YSR states hybridize and form a band inside the superconducting gap. Effective superconducting pairing may be induced in the band to generate stable topological phases. Two theoretical avenues were envisaged to achieve this. First, the superconductor substrate has no spin-orbit coupling, but helical ordering of impurity spins, opens a gap in the YSR band and generates effective topological $p$-wave superconductivity Choy et al. (2011); Martin and Morpurgo (2012); Braunecker and Simon (2013); Pientka et al. (2013); Klinovaja et al. (2013); Nadj-Perge et al. (2013); Pöyhönen et al. (2014); Kim et al. (2014). Alternatively, the superconductors have a strong spin-orbit coupling, which allows ferro- and antiferromagnetic chains of impurities to exhibit MBS at the edges Heimes et al. (2014, 2015); Brydon et al. (2015); Kim et al. (2015); Kaladzhyan et al. (2016). In this second scenario, the inversion symmetry is broken in the superconductor and the order parameter becomes a mixture of singlet and triplet pairing, Bauer and Sigrist (2012); Yip (2014); Samokhin (2015); Smidman et al. (2017) with the latter essential to generate effective $p$-wave superconductivity in the chains. Since the helical configuration of impurity spins proved harder to stabilize in experiments, the first detection of MBS from YSR states has followed the second avenue by employing chains of \ceFe atoms deposited on superconducting \cePb Nadj-Perge et al. (2014); Ruby et al. (2015); Pawlak et al. (2016); Feldman et al. (2016); Jeon et al. (2017). One experimental advantage of such setups is that scanning tunneling microscopy (STM) allowed direct imaging of the zero-energy local density of states characteristic for MBS, and, in principle, could be used to carry out braiding of MBS. Li et al. (2016) The spin-polarized STM is also used as a tool to detect MBS, Jeon et al. (2017); Jäck et al. (2019) since MBS wave functions have spin polarization sensitive to the impurity magnetization. Sticlet et al. (2012); He et al. (2014); Kotetes et al. (2015); Chirla and Moca (2016); Devillard et al. (2017); Li et al. (2018)

The search of topological superconducting phases supporting MBS has led lately to a surge in experimental investigations on YSR states Heinrich et al. (2018). In particular, our theoretical study is prompted by the recent experimental breakthroughs in detecting YSR states in superconducting \ceNbSe2. Ménard et al. (2015); Kezilebieke et al. (2018)

The bulk \ceNbSe2 is a layered transition metal dichalcogenide, with layers connected by weak van der Waals forces. Each layer is composed of three covalently bonded planar triangular lattices, two of \ceSe atoms, with one of \ceNb atoms intercalated between them, such that the unit cell has a trigonal prismatic structure. The bulk \ceNbSe2 is known to display a $3\times 3$ charge density wave ordering bellow $33\text{\,}\mathrm{K}$ and to become a superconductor bellow $7.2\text{\,}\mathrm{K}$. Recent technical advances have allowed to separate out monolayers of \ceNbSe2 from the bulk. Novoselov et al. (2005); Cao et al. (2015) It was shown that the critical temperature to the superconducting phase for monolayers is suppressed to approximately $1$ Wang et al. (2017) or $2\text{\,}\mathrm{K}$, Cao et al. (2015); Ugeda et al. (2016), while the critical temperature to realize charge density wave ordering increases to $145\text{\,}\mathrm{K}$ due to enhanced electron-phonon coupling. Xi et al. (2015); Ugeda et al. (2016) Since the monolayer has a thickness on the order of $1.2\text{\,}\mathrm{n}\mathrm{m}$, it may be considered quasi-two-dimensional, placing \ceNbSe2 in a select class of two-dimensional superconductors. Moreover, superconductivity in monolayer \ceNbSe2 was shown to have surprising features, most striking being its robustness to high in-plane magnetic fields. Xi et al. (2016); Xing et al. (2017); Sohn et al. (2018) This is due to the strong Ising spin-orbit coupling, which favors out-of-plane alignment of spins in opposite valleys. This results in Cooper pairs with pinned out-of-plane spins, which, together with suppressed orbital response in the two-dimensional limit, leads to an enhanced in-plane critical magnetic field, factors of magnitude above the Pauli paramagnetic limit. Xi et al. (2016); Xing et al. (2017); Sohn et al. (2018) These experimental findings are not unique to \ceNbSe2 and are shared by the group-VIB TMDs, with most notable representative being monolayer \ceMoS2. Xiao et al. (2012); Taniguchi et al. (2012); Yuan et al. (2014); Lu et al. (2015); Shi et al. (2015); Saito et al. (2015); Costanzo et al. (2016) While \ceNbSe2 is metallic in its normal state, the latter are usually semiconductors which are gated into the conduction band before becoming superconductors.

The presence of a strong spin-orbit coupling in the superconductor was shown above to be an essential ingredient for a class of proposals to realize topological superconductivity. This has motivated the recent experimental investigations of YSR states in \ceNbSe2. In-gap bound states were generated either from Fe adatoms Ménard et al. (2015) or magnetic molecules (cobalt phthalocyianine) Kezilebieke et al. (2018). There are two main reasons why this findings are of interest to us. First, \ceNbSe2 is a quasi-two-dimensional superconductor, which leads to a slower decay of YSR wave function in comparison with the three-dimensional case ($1/\sqrt{r}$ versus $1/r$, for distance smaller than the superconducting coherence length $\xi_{0}$). Ménard et al. (2015) This allows to couple the localized bound states easier as shown in Ref. Kezilebieke et al., 2018, where YSR dimers were first reported. This advances pave the way towards realization of MBS from coupled YSR states in \ceNbSe2. The second reason has to do with the exotic Ising spin-orbit coupling in \ceNbSe2. Due to broken in-plane reflection symmetry in the quasi-two-dimensional plane, but conserved out-of-plane mirror symmetry, the spin-orbit coupling favors a spin alignment normal to the monolayer for electrons at the Fermi surface in $K$ valleys, which leads to robust superconductivity to in-plane magnetic fields. This is different from most cases studied before Heimes et al. (2014, 2015); Brydon et al. (2015); Kim et al. (2015); Kaladzhyan et al. (2016) where spin-orbit coupling is Rashba type, which favors in-plane orientation of spins. For the latter, the depairing effect of in-plane magnetic fields is large, rendering problematic the manipulation of MBS with magnetic fields Li et al. (2016).

The theoretical studies on YSR states in superconducting \ceNbSe2 have for now limited themselves to consider only a few magnetic impurities. Besides the theoretical work accompanying Refs. Ménard et al., 2015; Kezilebieke et al., 2018, the interest in the effect of magnetic impurities on \ceNbSe2 goes back to early 2000 Flatté and Reynolds (2000), where a model without spin-orbit coupling was studied. More recently, the local density of states due to one to three impurities was investigated within the self-consistent formalism Ptok et al. (2017). Also the crossover from YSR states to Kondo impurity states came under focus Kezilebieke et al. (2019), as well as the interplay between charge density wave ordering and YSR states Liebhaber et al. (2019). Instead, our study is greatly indebted and extends the recent theoretical works Zhang and Aji (2016); Sharma and Tewari (2016) on monolayer superconducting TMDs like \ceMoS2, which shares with \ceNbSe2 the same $D_{3h}$ point-group symmetry. These works show that topological superconductivity is stabilized by ferromagnetic chains of impurities in superconductors with Ising spin-orbit coupling. The Fermi surface for \ceNbSe2 is however quite different from gated \ceMoS2 considered before, Zhang and Aji (2016); Sharma and Tewari (2016) with pockets not only at $K$, but also at $\Gamma$ valley [see Figs. 1(b) and 1(c)]. Moreover, the Fermi surface at $K$ valleys in \ceNbSe2 are spin split, rendering the problem complicated as more scattering channels are allowed. We show however that stable MBS are still possible, and provide a map of topological phases as a function of chain orientation on the two-dimensional superconducting substrate, the strength of magnetic exchange energy, and the distance between impurities.

The structure of the article is as follows. Section II introduces two effective models to describe superconductivity in monolayer \ceNbSe2. The first model is an effective tight-binding Hamiltonian for the conduction band. The model is mainly used in numerical simulations, and serves to derive the second, analytical model, which provides a better insight into generic features of topological superconductivity in \ceNbSe2. Section III investigates the addition of a periodic ferromagnetic chain of impurities in \ceNbSe2 and obtains an effective model for the chain in the limit of dilute impurities, energetically close to the middle of the superconducting gap (the deep-impurity limit). In the following, we obtain the topological phase diagram as a function of the distance between impurities and the orientation of impurity chains on the monolayer. Section IV displays topological phases with MBS using numerical methods, confirming the previous analytical calculations. The numerical simulations are done in the deep-impurity limit, and, beyond that, by using exact diagonalization of finite two-dimensional lattices. The last concluding Sec. V discusses our results and envisages possible extensions. Several appendices detail specific points in the paper regarding ab initio simulations, tight-binding models, indirect exchange interactions between impurity spins, additional numerical results, characteristic lengths from our calculations, etc.

II Model for superconducting monolayer \ceNbSe2

II.1 Tight-binding model

The band energy dispersion of \ceNbSe2 from experimentsXi et al. (2016); Bawden et al. (2016); Xing et al. (2017) and ab initio simulations (Appendix A) shows a single spin-split band crossing the Fermi level [Fig. 1(b)]. Since we aim to describe physics at the Fermi surface, at energy scales on the order of the superconducting gap $\Delta=$1\text{,}\mathrm{m}\mathrm{e}\mathrm{V}$$, far smaller than the bandwidth, it is necessary to devise a simple Hamiltonian to model electron dynamics in this band. From orbital-resolved density of states it follows that the largest contribution to this low-energy band comes from $d$ orbitals on the \ceNb atoms. Therefore in a first approximation the effect of \ceSe atoms is incorporated in an effective Hamiltonian for the triangular sublattice of \ceNb atoms Liu et al. (2013); He et al. (2018).

The Hamiltonian has to meet the symmetry requirements of the underlying triangular lattice [see Fig. 1(a)]. There is an in-plane threefold rotation symmetry represented by the operator $C_{3}=\exp(-\frac{2\pi i}{3}\sigma_{z})$. There is an out-of-plane mirror symmetry represented by $M_{z}=-i\sigma_{z}$. There are three in-plane reflection symmetry axes, the one about the $y$ axis, and rotations of it under the threefold symmetry. Characteristic for the TMDs and responsible for the Ising spin-orbit coupling is breaking three in-plane reflection symmetries (about $x$ axis and its threefold rotations). Breaking the in-plane reflection symmetry generates a spin-orbit coupling characterized by an out-of-plane spin orientation.

The Bogoliubov-de Gennes (BdG) Hamiltonian modeling superconductivity in the band crossing the Fermi level reads

[TABLE]

with the creation operators $C^{\dagger}_{\bm{k}}=(c^{\dagger}_{\bm{k}\uparrow},c^{\dagger}_{\bm{k}\downarrow},c_{\bm{k}\downarrow},-c_{\bm{k}\uparrow})$ and Bloch Hamiltonian:

[TABLE]

with the Pauli matrices in spin space denoted by $\bm{\sigma}$, and $\bm{\tau}$, in particle-hole space. The kinetic term $\xi_{k}$ follows from a Slater-Koster Slater and Koster (1954) expansion across the Brillouin zone. Variations of this term up to different orders in number of nearest-neighbor hoppings have been used before in describing triangular lattices for TMDs. Smith et al. (1985); Rossnagel et al. (2005); Inosov et al. (2008); Ménard et al. (2015) Here we have used up to fifth nearest-neighbor hopping terms (see Appendix B).

In addition, a spin-orbit coupling term $\Lambda_{\bm{k}}$ is added to the model due to broken in-plane reflection symmetries. This so-called Ising spin-orbit term favors an out-of-plane orientation of spins ($\sigma_{z}$) and it is responsible for the spin-splitting of the conduction band [Fig. 1(b)]. Since $\Lambda_{\bm{k}}$ is odd in momentum normal to the broken symmetry axes, it vanishes on the $\Gamma-M$ lines in the Brillouin zone (Appendix B). We will show that this leads to gapless superconducting phases when chains of magnetic impurities are oriented along the unbroken symmetry axes, identical to the theoretical prediction in \ceMoS2 TMDs Zhang and Aji (2016).

Details about the structure of the Bloch Hamiltonian and the fitting parameters to the ab initio data are relegated to Appendix B, together with a derivation of the spin-orbit term $\Lambda_{\bm{k}}$, in Appendix C. Further comparison with a three-band model for metallic \ceNbSe2 is contained in Appendix D.

II.2 Analytical model in the parabolic approximation

In order to study the superconducting properties of \ceNbSe2, we derive an effective model for \ceNbSe2 in its normal metallic phase. Since the Fermi surface of \ceNbSe2 shows electron pockets at $K$ and $\Gamma$ points, a valuable strategy is to obtain effective Hamiltonians around those points through a series expansion of the tight-binding Hamiltonian. In order to make analytical progress, the ensuing Hamiltonians are further mapped to effective Hamiltonians with a parabolic energy dispersion.

We start by showing the first terms at order $q^{3}$ near the high-symmetry points in the Brillouin zone, where $\bm{q}$ is momentum measured from either $\pm\bm{K}_{n}$ or $\Gamma$. At the two $K$ valleys, the Hamiltonian reads

[TABLE]

where $q=|\bm{q}|$, the effective mass at $K$ valleys is $m_{K}$, with chemical potential $\mu_{K}$. The valley index is denoted throughout by $\eta=\pm$, and the angle between $\bm{q}$ and the $x$ axis, by $\chi$. The cubic term parametrized by $w$ represents the warping of the Fermi surface. The spin-orbit coupling amplitude is parametrized by the set of $\kappa_{i}$, and it is consistent with low-energy continuum models developed for triangular lattices with an extremum at $K$ valleys as for \ceMoS2, MoSe2, or \ceMoTe2. Wang and Wu (2014); Wang et al. (2018)

For small momenta $\bm{q}$ near the $\Gamma$ point, the Hamiltonian reads

[TABLE]

recovering the continuum Hamiltonians derived in Ref. He et al., 2018.

In the parabolic approximation, we restrict ourselves only at quadratic terms in momentum. Therefore the \ceNbSe2 valleys have circular Fermi surfaces and important quantities to our analysis such as Fermi momentum and velocity simplify to scalars. In order to keep close to the experimental reality, we derive the effective Fermi momentum and velocity at each valley from the full tight-binding model with warped Fermi surfaces. More precisely, the effective $v_{F}$ and $k_{F}$ are obtained form the angular average at the Fermi surface $\langle*\rangle=\frac{1}{2\pi}\int\!*\,d\chi$ (see Appendix E).

The retarded Green’s function for the superconductor reads:

[TABLE]

with $\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}$ and system size $\Omega$. The sum carries over the total number of momentum states $\bm{k}$ in the first Brillouin zone. For energies $E\sim\Delta\sim$1\text{,}\mathrm{m}\mathrm{e}\mathrm{V}, much smaller than the normal metal bandwidth ($\sim$1\text{\,}\mathrm{e}\mathrm{V}), only the contribution from momenta at the Fermi surface is important. Therefore the Green’s function may be decomposed in a sum over $K$ and $\Gamma$ pockets:

[TABLE]

with $G_{ij}^{K}(E)$, the Green’s function for electrons at $\pm\bm{K}_{n}$ valleys, and $G_{ij}^{\Gamma}(E)$, the Green’s function at the $\Gamma$ valley.

The Green’s function for the six $K$ valleys is modeled as:

[TABLE]

where $\bm{q}$ is a small momentum near the $\pm\bm{K}_{n}$ valley. The factor $1/3$ comes from an average over the three equivalent valleys for a fixed valley index $\eta$. We have used the spectral decomposition of the Green’s function in the helicity subspace $\lambda$ Gor’kov and Rashba (2001) to obtain:

[TABLE]

The product is understood here as a Kronecker product between particle-hole space and spin space, with $1$ being interpreted when necessary as identity in the respective space. The normal dispersion $\xi_{q,\eta\lambda}$ at the $\pm\bm{K}_{n}$ valleys is indexed by valley $\eta$ and helicity $\lambda$, which models polarized electrons with spin up ($\lambda=+$) and down ($\lambda=-$). From now on, we will denote for notational simplicity the sums over valley and helicity as $\sum_{n\eta\lambda}$. The normal dispersion is approximated using the leading term $\kappa_{0}$ in spin-orbit coupling, $\xi_{q,\eta\lambda}=\hbar^{2}q^{2}/2m_{K}-\mu_{K}+\eta\lambda\kappa_{0}$. The second term ($\lambda\sigma_{z}$) in Eq. (8) is due to spin-orbit coupling and generates spin-triplet superconducting correlations, necessary to engineer topological superconductivity. Brydon et al. (2015)

At the $\Gamma$ pocket, there is no contribution from spin-orbit coupling to $\mathcal{O}(q^{2})$. Only the cubic term contributes in creating spin-triplet superconducting pairing and helps stabilizing topological nontrivial phases. Nevertheless, the physics is dominated by the contribution from $K$ valleys, since the spin-orbit coupling at $\pm\bm{K}_{n}$ is by an order of magnitude larger than at $\Gamma$. 111from the tight-binding model, $\kappa_{0}\simeq$0.078\text{,}\mathrm{e}\mathrm{V} and $\gamma_{3}/a^{3}\simeq$0.0027\text{\,}\mathrm{e}\mathrm{V}, with $a=$0.344\text{,}\mathrm{n}\mathrm{m}$$ the distance between \ceNb atoms. Therefore we neglect in the following the cubic contribution from the $\Gamma$ valley. Under this approximation, the Green’s function at the $\Gamma$ point reads:

[TABLE]

with $\xi_{q}^{\Gamma}=\hbar^{2}q^{2}/2m_{\Gamma}-\mu_{\Gamma}$.

All Green’s functions may be written in terms of integrals $I_{0}(E,\bm{r}_{ij})$ and $I_{1}(E,\bm{r}_{ij})$, which are analytically solvable for parabolic bands near the Fermi surface:

[TABLE]

The $I$ integrals are a sum of Bessel and Struve functions in agreement with other works on two-dimensional superconductors Brydon et al. (2015); Heimes et al. (2015); Zhang and Aji (2016) (see Appendix F). Due to time-reversal symmetry of the pure superconductor, the integrals at $K$ valleys obey $I_{0,1}^{\eta\lambda}=I_{0,1}^{-\eta,-\lambda}$. The equal-position total Green’s function is trivial:

[TABLE]

with $\rho(0)$ the density of states at the Fermi surface for both spin projections.

III YSR states from a ferromagnetic chain of impurities

III.1 General formulation

In this section, we review the method to obtain the equation for the YSR bound states for a ferromagnetic chain of impurities. Let us consider $N$ identical impurities coupled to the superconducting substrate and placed at equally distanced positions $\bm{r}_{j}$, with a spacing $d$. By convention, the monolayer sits in $xy$ plane, and the direction of the linear chain is parametrized by the angle $\phi$ between the chain and $x$ axis. Each impurity spin $\bm{S}$ is oriented along direction $\bm{n}$, specified by angles $\theta$ and $\varphi$, $\bm{n}=(\cos\varphi\sin\theta,\sin\varphi\sin\theta,\cos\theta)$. We assume that in the following that the impurities are exchange-coupled identically to the substrate and the interaction is purely local in space. Under these assumptions, the interaction Hamiltonian between the impurities and the substrate reads

[TABLE]

where $J>0$ is the exchange coupling constant between impurity and substrate, while $S$ is the impurity spin amplitude. Within a classical approximation, the spins are treated as vectors, with amplitude $S\to\infty$, and small coupling constant, $J\to 0$, with magnetic exchange energy $JS$ finite. Shiba (1968) Note that we neglect in the following the nonmagnetic interaction between impurities and substrate, since this potential scattering has typically a small influence on YSR state energy in $s$-wave superconductors. Balatsky et al. (2006)

The preference for the ferromagnetic ordering in the dilute impurity limit has to be derived from a minimization of the indirect exchange Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction energy between spins. Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957) In Appendix G, we discuss the classical RKKY Hamiltonian showing that the Ising spin-orbit coupling induces, besides Heisenberg and Ising spin-spin interactions, also a Dzyaloshinskii-–Moriya interaction, which favor noncollinear ordering of impurity spins (similar to the case of \ceMoS2). Parhizgar et al. (2013) The results show that the RKKY interaction is direction dependent due to the $K$ valley electrons. Therefore the ground state admits a variety of spin orders, depending on the distance between spins and chain orientation. The analysis in Appendix G neglects additional crystal field effects which favors collinear order of spins Heimes et al. (2014, 2015) and could stabilize the spin order in Eq. (12). Therefore, at the moment, Eq. (12) remains an ansatz whose soundness requires additional numerical works or experimental investigations using spin-polarized STM. Zhou et al. (2010); Khajetoorians et al. (2012); Nadj-Perge et al. (2014); Jeon et al. (2017); Jäck et al. (2019)

The stationary Schrödinger equation reads:

[TABLE]

Using the Fourier transform $\psi_{\bm{r}}=\frac{1}{\sqrt{\Omega}}\sum_{\bm{k}}\psi(\bm{k})e^{i\bm{k}\cdot\bm{r}}$ and integrating over both sides in Eq. (13) with $\frac{1}{\sqrt{\Omega}}\int d\bm{r}e^{-i\bm{q}\cdot\bm{r}}$ yields

[TABLE]

Using the free Green’s function in momentum space $G(E,\bm{q})^{-1}=E-H_{0}(\bm{q})$, it follows the relation to determine $\psi(\bm{q})$:

[TABLE]

Integrating over momentum with $\frac{1}{\sqrt{\Omega}}\int d\bm{q}e^{i\bm{q}\cdot\bm{r}_{i}}$ and using the Green’s function in real space, which for equidistant impurities depends only on the distance between impurities located at $\bm{r}_{i}$ and $\bm{r}_{j}$,

[TABLE]

we obtain the equation for the YSR states:

[TABLE]

The equation determines the wave function of YSR bound states localized at impurity positions, depending on the impurity spin orientations, the magnitude of the exchange coupling energy, and the substrate superconductor Green’s function. Due to the presence of the Green’s function, it is however a nonlinear eigenvalue equation in the energy $E$. To make analytical progress in determining the YSR wave functions and energies, it is necessary to perform additional approximations.

III.2 Deep-impurity limit

The YSR bound states may be obtained through a linearization of Eq. (17) in the so-called deep-impurity limit. This approximation is correct when the YSR band is close to the Fermi level, when $|E|\ll|\Delta|$. Such approximation is a priori known to be physically applicable in the case of dilute impurities, where only weak interactions exist between neighboring impurity sites. Under this condition, the picture from the single magnetic impurity physics is expected to hold, i.e., the YSR energy (and therefore the entire YSR band) may be brought at the Fermi level when free tuning the magnetic exchange energy $JS$ [see below Eq. (23)]. In the present section, we will therefore assume dilute weakly coupled impurities ($k_{F}d\gg 1$) (for a dense-impurity approach, see Ref. Peng et al., 2015) with energy $|E|\ll|\Delta|$. Therefore one expands Eq. (17) both in the energetic distance of the bound states to the Fermi level and in a small coupling between impurity sites. Keeping Eq. (17) at first order in the small parameters implies a linear dependence on energy in the same-position Green’s function, while the Green’s function controlling hopping is estimated at zero order in YSR energy: Pientka et al. (2013)

[TABLE]

where $\partial_{E}G_{ij}(0)$ denotes $\partial_{E}G_{ij}(E)|_{E=0}$ and $\delta_{ij}$ is the Kronecker delta. Therefore a generalized eigenvalue equation follows from Eq. (17):

[TABLE]

This result is further simplified with the use of Eqs. (6) and (11) to cast the problem to a simple tight-binding stationary Schrödinger equation:

[TABLE]

with the Hamiltonian

[TABLE]

with $\alpha=\pi\rho(0)JS/2$, and $\rho(0)$, the total density of states at the Fermi surface. The on-site terms $(\delta_{ij})$ in the equation are explicitly separated from the hopping terms $(1-\delta_{ij})$, which depend on the Green’s function for $\Gamma$ and $K$ valley electrons calculated before (II.2).

If there is only one impurity added to the system, then the tight-binding equation contains only the on-site term:

[TABLE]

The YSR bound states follow readily, yielding the classical solutions Balatsky et al. (2006):

[TABLE]

The two solutions cross each other at zero energy ($\alpha=1$) and a quantum phase transition takes place, resulting in a quasiparticle excited on the YSR state. The two solutions are:

[TABLE]

The tight-binding model simplifies when projecting it in the single impurity basis $\{\psi_{+},\psi_{-}\}$. The resulting Hamiltonian has the same block structure as in Ref. Zhang and Aji, 2016. Moreover, we recover their results if we conflate helicity and valley index and eliminate the contribution from the $\Gamma$ pocket (as it would be correct in gated \ceMoS2 near the lowest point in the conduction band). The effective Hamiltonian $\bar{H}$ reads:

[TABLE]

with

[TABLE]

The above Hamiltonian describes quasiparticle dynamics inside the YSR band, at subgap energies (note the scaling of all matrix elements with $\Delta$). The terms proportional to $\lambda$ are due to spin-triplet correlations in the Green’s function (8), and lead to antisymmetric “effective magnetic field” $b$ and effective superconducting pairing $\bar{\Delta}$, $b_{ij}=-b_{ji}$ and $\bar{\Delta}_{ij}=-\bar{\Delta}_{ji}$. The antisymmetric pairing $\bar{\Delta}$ is reminiscent of the $p$-wave pairing in the Kitaev chain Kitaev (2001) and is essential to obtain topological superconductivity. The integrals $I_{0,1}$ are expanded in the dilute limit $k_{F}d\gg 1$, yielding the following asymptotic behavior as a function of distance between impurities (see Appendix F):

[TABLE]

with superconducting coherence length $\xi_{0}=\hbar v_{F}/\Delta$. The slow power-law decay of $I_{0,1}(0)$, up to distances on the order of $\xi_{0}\sim 700a$ in our system (see Appendix E), results in a long-range hopping tight-binding model with direct hopping between sites at distance of $\mathcal{O}(\xi_{0})$.

Two particular cases where there is no topological superconductivity may be directly inferred from the effective tight-binding Hamiltonian (25). First, if the spins are entirely out-of-plane, then $\sin\theta=0$, and the effective gap $\bar{\Delta}$ closes in the model, making it impossible to have topological phases. Second, the spin-orbit coupling which provides spin-triplet pairing is odd in momentum normal to the broken reflection symmetry axes ($x$ axis and its threefold rotations). Therefore it must vanish along $y$ (and its threefold rotations), which leads to vanishing spin-orbit splitting along $\Gamma-M$ lines in the Brillouin zone [see Fig. 1(b)]. When the chain with direction given by $\bm{r}_{ij}$ sits along the $y$ axis or its threefold rotations, there is no contribution from triplet pairing and both $b$ and the effective gap $\bar{\Delta}$ vanish. This leads again to a gapless trivial phase. The above argument may be checked by writing $\bm{K}_{n}\cdot\bm{r}_{ij}=K|\bm{r}_{ij}|\cos(\frac{2\pi n}{3}-\phi)$, with the $K=4\pi/3a$, the distance from $\Gamma$ to $\pm\bm{K}_{n}$. Hence, using the time-reversal symmetry property of $I$ integrals, $I^{\eta\lambda}_{0,1}=I^{-\eta,-\lambda}_{0,1}$, it follows that the effective magnetic field and superconducting gap are proportional to $\sum_{n=0}^{2}\sin(K|\bm{r}_{ij}|\cos(\frac{2\pi n}{3}-\phi))$. Therefore they indeed vanish on the unbroken symmetry axes $\phi=(2m+1)\pi/6$, with $m$ an arbitrary integer.

The topological phases are determined by studying the Bloch Hamiltonian, obtained after a Fourier transform for momenta $k$ along the chain:

[TABLE]

with $jd$ the distance between impurities, with spacing $d$. The effective Bloch Hamiltonian for the chain becomes:

[TABLE]

with

[TABLE]

and

[TABLE]

The expressions make use of the polylogarithm function $\textrm{Li}_{\nu}(z)=\sum_{j=1}^{\infty}z^{j}/j^{\nu}$. The superconducting coherence length at pocket $x$ ($\Gamma$ or $K$) reads $\xi_{0}^{x}=\hbar v_{F}^{x}/\Delta$. The values of Fermi momenta $k_{F}^{x}$, Fermi velocities $v_{F}^{x}$, and densities of states, at respective pockets, are explicit in Appendix E. Note that $\rho_{\Gamma}(0)$ is the total density of states at $\Gamma$, for both spin projections, while $\rho_{\eta\lambda}(0)$ is the density of states for valley $\eta$ and helicity $\lambda$. Due to time-reversal invariance in pure \ceNbSe2, the following relation holds: $\rho_{\eta\lambda}(0)=\rho_{-\eta,-\lambda}(0)$.

III.3 Topological phases

Due to $C_{3}$ symmetry, all information about the Hamiltonian is contained for chain orientation angle $\phi$ spanning the interval $[0,\pi/3]$. Moreover, due to unbroken in-plane reflection symmetry about $\pi/6$, there is the same physics at $\pi/6+\phi^{\prime}$ and $\pi/6-\phi^{\prime}$. Therefore it is sufficient to study the effective Hamiltonian in the restricted window $[0,\pi/6]$, where we compute in the following the topological phase diagrams.

The effective Hamiltonian $\bar{H}$ has particle-hole symmetry represented by the operator $\tau_{x}\mathcal{K}$. The presence of the magnetic field term $b(k)$ breaks time-reversal symmetry, and therefore the one-dimensional Hamiltonian is generally in D class, Schnyder et al. (2008); Ryu et al. (2010) described by the topological number:

[TABLE]

where $\textrm{Pf}(A)$ is the Pfaffian of matrix $A$. In numerics, we usually choose an in-plane ordering of impurity spins ($\theta=\pi/2$), since $b(k)$ vanishes in this case, ensuring that the YSR band is gapped. Closing of the gap occurs at possible topological phase transitions or when the impurity chain sits along the in-plane unbroken symmetry lines ($y$ axis and its threefold rotations). Note however that the case of in-plane spins is rather particular Pientka et al. (2013, 2014) because an effective time-reversal symmetry gets restored. The Hamiltonian becomes real under a unitary transformation $U$ which gauges out the superconducting phase $e^{i\varphi}$:

[TABLE]

Consequently, the Hamiltonian is time-reversal symmetric under the action of $\tau_{z}\mathcal{K}$ and is classified in the one-dimensional BDI class. In this case, the topological invariant is $\mathbb{Z}$, Schnyder et al. (2008); Ryu et al. (2010) and multiple MBS are allowed in the model. The presence of BDI phases in two-dimensional superconductors with in-gap YSR states have been studied in several other works where the substrate has Rashba spin-orbit coupling. Pientka et al. (2013, 2014); Heimes et al. (2014); Pöyhönen et al. (2014)

Such behavior is not surprising and is apparent from symmetry considerations of the original Hamiltonian (13). In that case, the pure superconductor Hamiltonian $H_{0}$ is time-reversal invariant under the action of $T=i\sigma_{y}\mathcal{K}$. The impurity spins contributing to $H_{1}$ act as a magnetic field breaking TRS. If the impurity spins are constrained to point in plane, then they also break the out-of-plane mirror symmetry $M_{z}=-i\sigma_{z}$. Nevertheless, the combination $M_{z}T=-i\sigma_{x}\mathcal{K}$ restores the correct transformation of spins under the time-reversal symmetry, placing the total Hamiltonian $H$ in the two-dimensional chiral-symmetric BDI class. This conclusion holds more generally for a large class of TMDs including gated \ceMoS2, MoSe2, etc. Wang et al. (2018) Since the above argument is valid for any completely in-plane orientation of spins, all effective one-dimensional Hamiltonians derived from above and living in the plane, remain in BDI class. Therefore, indeed, under this constraint on spins, the effective impurity chain Hamiltonian is in the one-dimensional BDI class.

In the present paper, we have limited ourselves to cases where only two MBS are present in the system, such that the Majorana number (32) remains useful in characterizing the topological phases also for the chiral-symmetric BDI system. For more generic cases of purely in-plane spins, the topological characterization requires the measurement of a winding number. Tewari and Sau (2012) It is beyond the scope of the present study to inventory possible high-winding number phases of the purely planar orientation of spins. It is expected that generally there would be a nonvanishing out-of-plane component to spin projections, which renders high-winding number phases unlikely to be observed in experiments.

Figure 2 shows the resulting topological phase diagram as a function of distance between impurities and chain orientation $\varphi$. The magnetic exchange energy $JS$ is fixed at $2/\pi\rho(0)$, the point where the single-impurity YSR crosses zero, where we expect that the deep-impurity limit holds. Despite the complicated structure of the phase diagram, several extended connected regions of nontrivial topological phases are present. In particular, we will focus in the following numerical investigations on the region centered at $\phi=0$ and $d=3a$, and show that the topological phases are also robust when varying the magnetic exchange energy. Some additional results are contained in Appendix H, when the distance between impurities is varied in a range relevant for the case of magnetic molecular chains with, e.g., phthalocyanine.

IV Numerical results

In the following, we study the impurity problem using numerical methods, as a direct proof of MBS presence. Section IV.1 investigates the deep-impurity limit in order to obtain the energy and wave functions of YSR states without resorting to the parabolic approximation of energy dispersion at the \ceNbSe2 valleys. Section IV.2 dispenses with both the parabolic approximation and the deep-impurity limit by using exact diagonalization of a finite-size two-dimensional \ceNbSe2 superconductor with impurities deposited on system sites. All simulations manifest Majorana states at values predicted by the analytical study.

The numerical simulations come with one technical inconvenience. Due to the small superconducting gap $\Delta=$1\text{,}\mathrm{m}\mathrm{e}\mathrm{V} in \ceNbSe2, the coherence lengths are on the order of $\xi_{0}\sim 700a$ (see Appendix [E](#A5)). To avoid unphysical finite size effects due to level spacing being larger that the superconducting gap, it is necessary for the tight-binding lattices to have linear lengths exceeding $\xi_{0}$. This imposes a strain of numerical simulations, which is usually remedied by artificially increasing the superconducting pairing amplitude. This increase does no modify qualitatively single-impurity physics since YSR energies are scaled by $\Delta$. However, in the case of a chain of impurities, this has the effect to decrease the number of hopping terms involved in the effective Hamiltonian, and tends to favor the growth of topological phases size in parameter space (see Appendix [H](#A8)). Here we use in numerics $\Delta=$10\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}, i.e., a tenfold decrease of the coherence length.

IV.1 Effective 1D models in the deep-impurity limit

Here we obtain the YSR bands and MBS in the deep-impurity limit, without resorting to the parabolic approximation to the energy dispersion near the valleys.

The equal-position Green’s function (11) is trivial and its energy derivative follows analytically. Therefore only left-hand side in Eq. (19) needs to be computed numerically. After denoting $\Xi_{j}=\bm{\sigma}\cdot\bm{n}\psi_{j}$, and rearranging the terms, the equation to determine YSR states becomes:

[TABLE]

The pure superconductor Green’s function $G_{ij}(0)$ from Eq. (5) is calculated by numerical integration over the two-dimensional Brillouin zone. We consider a periodic patch of sites defined by vectors $m\bm{a}_{1}$ and $n\bm{a}_{2}$, with $\bm{a}_{1}$ and $\bm{a}_{2}$ primitive vectors of the lattice, and $m,n$ integers. This induces a discretization of the Brillouin zone, which contains $m\times n$ lattice momenta. In numerical experiments, we usually choose $m=n=600$.

It was sufficient to consider small chains of 60 impurities, in order to have well-localized Majorana states. The impurities sit at lattice positions and form a chain in the middle of the periodic patch. For definiteness, we have considered a chain of impurities along $x$ axis ($\phi=0$) with spins pointing in the positive $x$ direction, spaced by $d=3a$. The results are shown in Fig. 3 (upper panel), which explores the behavior of the YSR states and MBS in a range of magnetic exchange energies. There are two transition points where the YSR band closes, with an extended nontrivial topological phase in between, starting around $JS\sim$0.2\text{,}\mathrm{e}\mathrm{V}$$. The topological phases are marked by the presence of oscillating MBS, separated by a gap from the rest of the YSR band. The oscillating behavior of MBS is due to finite size effects, since oscillations decrease with the length of impurity chain and the size of the periodic substrate. For a fixed value of $JS$, we present in Fig. 3 (lower panel) the two MBS eigenvalues, and the amplitude of their functions, which is localized at the edges of the impurity chain.

The numerical simulations confirm the analytical prediction that the YSR gap closes for completely out-of-plane orientation of impurity spins Fig. 4(a) and for impurity chains aligned along the conserved in-plane reflection symmetry axes Fig. 4(b). In both cases, no MBS develop in the system. In contrast, Fig. 4(c) shows a generic case with MBS for an arbitrary chain orientation, away from any symmetry axis of the triangular lattice, with an arbitrary orientation of impurity spins (with some nonzero out-of-plane component), and impurity positions incommensurate with the underlying lattice.

IV.2 Exact diagonalization of 2D tight-binding models

In this section, MBS presence is proved through exact diagonalization of a finite triangular \ceNbSe2 lattice. To render simulations comparable to those of the previous section, we define impurities along $x$ axis, with spins polarized in the $x$ direction, and with distance between impurities $d=3a$. The lattice is a rectangle of size approximately $140a\times 70a$, such that it is twice longer than the estimated superconducting coherence length. On the lattice, we fit a chain of 40 impurities arranged as in the lower panel of Fig. 5, which shows a cutout (limited in $y$ direction) of the system. A scan over different magnetic exchange energies shows MBS over a region comparable to that in the previous section. Nevertheless, the finite size effects are more noticeable, compared to those in the larger lattices previously explored. This leads to larger oscillations of MBS at zero energy and to a higher $JS$ where transition between trivial and nontrivial phases takes place. We have checked by scaling up from smaller lattices, that the topological phases in $JS$ decreases, presumably towards the results shown in Fig. 3. In the latter case of the previous section, we could check that the positions of topological phase transition do not change for larger lattices, but MBS oscillations are, as expected, reduced.

The inset of Fig. 5 shows the MBS eigenvalues, while the lower panel, the corresponding amplitudes of the wave functions for a fixed value of $JS$ in the nontrivial region. The Majorana wave functions show strong localization at the edges of the impurity chains (lower panel in Fig. 5). Finite size effects are present in the slight overlap of MBS, leading to a small energy split near zero energy, as well as in the reflection of the MBS wave function from the vertical edges of the lattice (Fig. 5, lower panel). Simulations on the finite tight-binding systems have been performed using the Kwant package Groth et al. (2014).

V Discussion

In this paper, we have studied topological phases supporting Majorana bound states due to ferromagnetic chains of impurities on monolayer \ceNbSe2 superconductors. Our analysis is based on realistic tight-binding modeling for the conduction band of the superconducting substrate using parameters extracted from fits to ab initio band structure calculations. Our results indicate that the strength of Ising spin-orbit coupling in the \ceNbSe2 superconductor is sufficient to stabilize topological phases supporting Majorana bound states localized at the end of the impurity chains. Similar to previous studies in transition metal dichalcogenides, such as \ceMoS2, Zhang and Aji (2016); Sharma and Tewari (2016) we also showed that Majorana bound states are expected for generic orientation of impurity chains on the substrate and impurity spin projections. Nevertheless, the Yu-Shiba-Rusinov bands close when the chains are oriented on the lines of unbroken in-plane reflection symmetry, rendering Majorana phases impossible in such cases. Moreover, it is necessary for the existence of topological phases that impurity spins have a finite in-plane component. Since topological phases in \ceNbSe2 exist for spin-split bands at the $K$ valleys, it is reasonable to expect that topological phases would develop also in \ceMoS2, \ceMoSe2, etc. for generic gating in the conduction band, when two Fermi surfaces develop at each valley.

Our findings are encouraging for the recent experimental forays in generating YSR states in \ceNbSe2 Ménard et al. (2015); Kezilebieke et al. (2018, 2019). In particular, a valid experimental work program would be to investigate whether self-assembled chains of magnetic molecules offers a robust alternative to designing one-dimensional topological superconductivity. Nevertheless, additional studies are needed to determine the magnitude of magnetic exchange energies between the superconducting substrate and, e.g., magnetic porphyrin or phtalocyanine molecules. Our work has determined that a necessary condition to generate nontrivial topological phases consists in having magnetic exchange energies of at least $JS\sim$0.2\text{,}\mathrm{e}\mathrm{V}$$.

Finally, a remaining open problem is to evaluate the effect of charge density wave ordering on the formation of the Yu-Shiba-Rusinov band and on the topological phases in \ceNbSe2. Recent experiments Liebhaber et al. (2019) demonstrate strong dependence of an YSR state energy on modulations in the local environment of the impurity. Therefore we expect that the YSR band will be modified, and it remains to quantify the impact on topological superconducting phases.

Acknowledgements.

We thank C.P. Moca, F. von Oppen, and L.P. Zârbo, for illuminating discussions, and to M.P. Nowak for going over the final manuscript. This research was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CNCS-UEFISCDI, with project No. PN-III-P1-1.2-PCCDI-2017-0338.

Appendix A Ab initio calculation

The ab initio simulations were performed using the Siesta code Ordejón et al. (1996); Soler et al. (2002) which uses norm-conserving pseudopotentials Troullier and Martins (1991) and expands the wave functions of valence electrons using linear combinations of atomic orbitals (LCAO). We used GGA/PBE. Perdew et al. (1996) as exchange-correlation function. The LCAO basis set is double-zeta polarized with an energy shift of $100\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. The lattice constant is $0.344\text{\,}\mathrm{n}\mathrm{m}$ and the integration in the Brillouin zone uses a $40\times 40$ Monkhorst-Pack grid.

The resulting band structure along the high-symmetry lines is included in Fig. 1(b). The density of states for the conduction band is due almost entirely to $4p$ \ceSe and $4d$ \ceNb orbitals. The orbital-resolved density of states for the conduction band is given in Fig. 6(a). Since the dominant contribution is due to \ceNb atoms, it is justified to devise an effective model where the effect of \ceSe atoms is included in a Hamiltonian (2) for hopping only in the triangular lattice of \ceNb atoms.

Appendix B Low-energy tight-binding model

The explicit form for the Bloch-BdG Hamiltonian modeling monolayer \ceNbSe2 from Eq. (2) reads

[TABLE]

where $\alpha=k_{x}a/2$ and $\beta=\sqrt{3}k_{y}a/2$.

The Hamiltonian parameters are obtained through a fit to the ab initio data along the high-symmetry lines in BZ. The resulting parameters are presented in Table 1.

Appendix C Derivation of the spin-orbit coupling term

An out-of-plane spin-orbit term is allowed due to breaking of the in-plane reflection symmetry $M_{x}=-i\sigma_{x}$. The spin-orbit coupling term $\Lambda_{\bm{k}}$ derives from the spin-orbit Hamiltonian for electrons on a lattice:

[TABLE]

where $m$ is the effective electron mass, $c$ the speed of light and $V$, the lattice potential. The electronic cloud centered at $\ce{Se}$ atoms on top and at the bottom of the sheet of $\ce{Nb}$ atoms exert a Coulomb repulsion on electrons hopping in plane between the $\ce{Nb}$ atoms. The effect is maximal for hopping along $\bm{d}_{j}$, when the force normal to the hopping is largest. For hopping along the other directions, the projection of $\ce{Se}$ atoms on the $\ce{Nb}$ plane passes close to the hopping path, such that there is compensation between the effect of the charges above and below the plane. For example, the compensation is perfect for next nearest-neighbor hopping, leading to a vanishing force [see Fig. 6(b)]. Under these condition, we investigate first the effect due to nearest-neighbor hopping.

Since the out-of-plane reflection symmetry $M_{z}=-i\sigma_{z}$ is preserved, the spin-orbit term may have only a component along $\sigma_{z}$. Therefore the effect of spin-orbit term to nearest-neighbor hopping becomes in the second quantization:

[TABLE]

where the sum over $i$ runs over all lattice sites and sum over spin indices is implied. We used the symmetry to factor out the equal magnitude of forces $|\bm{F}_{\perp}|=|\bm{F}_{j\perp}|$ for any $j$. Denoting the constant term $\lambda_{1}=\hbar|\bm{F}_{\perp}|/4m^{2}c^{2}a$ and Fourier transforming into momentum space, we find a contribution from spin-orbit coupling:

[TABLE]

Such spin-orbit coupling was used before Xi et al. (2016); Sharma and Tewari (2016) to describe \ceMoS2 type compounds modeled on a triangular lattice and may be derived generically for noncentrosymmetric superconductors Samokhin (2007). As noted, the contribution to next nearest-neighbor hopping vanishes. For fitting purposes, we also add the contribution of spin-orbit coupling to third nearest-neighbor hoppings, leading to the form that we used in our tight-binding Hamiltonian:

[TABLE]

with $\alpha=k_{x}a/2$ and $\beta=\sqrt{3}k_{y}a/2$. As expected, time-reversal symmetry of pure \ceNbSe2, together with the conserved reflection symmetry $M_{y}$, leads to vanishing spin-orbit coupling along $y$-axis (and its in-plane rotations under the threefold symmetry).

Appendix D Comparison with the three-band model

More elaborate models have been devised in order to reproduce the low-energy band structure of \ceNbSe2 [as well as for other monolayer TMDs \ceMX2, with metal \ceM in {\ceMo, W} and chalcogen \ceX in {\ceS, Se,Te}], the most used being a model for hopping only between three $d$ orbitals of \ceNb in a triangular lattice. Liu et al. (2013); He et al. (2018) This model has more fitting parameters, allowing it to better reproduce the energy dispersion from the ab initio data. In addition, it also provides a fit for the first two bands above the Fermi level-crossing band. In this appendix, we compare results between one-band (2) and three-band models in case of a single magnetic impurity in the system.

The three-band Bloch Hamiltonian for superconducting \ceNbSe2 reads

[TABLE]

The product of matrices is a Kronecker product, with Pauli matrices $\bm{\tau}$ in particle-hole space, and $\bm{\sigma}$ in spin space, with the 3-by-3 matrices in the orbital space. Absence of one matrix in a product is understood as presence of identity matrix in the respective space. The matrix elements of the kinetic term are:

[TABLE]

where $\alpha=k_{x}a/2$ and $\beta=\sqrt{3}k_{y}a/2$. With the three-band model, the fit reproduces the conduction band together with the first two higher bands coming from the ab initio calculations. Note that the model may be used to fit perfectly only the conduction band, at the expense of a poor fit for the higher energy bands. Instead, we fit all the three above mentioned bands, which leads to small deviations of the fitted model in the conduction band, mostly at the minimum between $\Gamma$ and $K$ valleys. The resulting fitting parameters are shown in Table 2, together with the energy dispersion along the high-symmetry lines, in Fig. 7(a).

To obtain YSR states for a single impurity, we first connect a site in the substrate to one impurity. The classical spin orientation is irrelevant to bound state energy in case of a single impurity and we fix it here normal to the monolayer. From the three-band model, we extract the density of states at the Fermi surface ($\rho(0)\simeq 2.468$ states$/$\mathrm{e}\mathrm{V}) and use Eq. ([23](#S3.E23)) to determine the YSR energies as a function of $JS$ in Fig. [7](#A4.F7)(b). In the one-band case we perform the calculation of YSR energies numerically, using the methods described in Sec. [IV.1](#S4.SS1). The results are remarkably similar due to only a small variation in the Fermi surface density of states between the two models. They both identify the crossing of YSR states at zero energy at $JS\approx$0.25\text{\,}\mathrm{e}\mathrm{V}. These results are encouraging for the focus in our paper on the simpler one-band model, which is amenable to analytical treatment in the parabolic approximation, and provides a less CPU time-consuming numerical implementation.

Appendix E Parameters in the analytical calculation

The only free parameters in the analytical calculation are the chain orientation on top of the two-dimensional superconductor, given by the angle $\phi$, the orientation of the classical impurity spins, the distance $d$ between adjacent impurities, and the amplitude of magnetic exchange energy $JS$. The other constants are fixed from the tight-binding model and is extracted in an investigation of the energy dispersion at the $\pm K$ and $\Gamma$ pockets.

Near the $\Gamma$ pocket, the Fermi surface is close to spherical. From a fit of the analytical model to the ab initio data, we obtain:

[TABLE]

Near the $\pm K$ pockets, the surface is hexagonally warped even in the absence of spin-orbit coupling. To obtain constants appropriate for a spherically symmetric effective Hamiltonian at the pockets, we perform the angular average of the Fermi surface $\bar{x}_{F}=\frac{1}{2\pi}\int d\chi x_{F}(\chi)$, with $x$ either Fermi velocity or momentum. This yields effective Fermi momenta and velocities (omitting the bar):

[TABLE]

From the split of the Fermi sheets at the $K$ point, we extract the spin-orbit momentum:

[TABLE]

The above results imply that the superconducting coherence length is on the order of $\xi_{0}\sim\hbar v_{F}/\Delta\sim 700a$ with $a=$0.344\text{,}\mathrm{n}\mathrm{m}$$.

Using numerical integration at the Fermi surface with a broadening of the Green’s function $\epsilon=10^{-3}$, we obtain the densities of states:

[TABLE]

in states per electronvolt and primitive unit cell area $\Omega_{0}=\sqrt{3}a^{2}/2$. Note that $\rho_{\Gamma}(0)$ and $\rho(0)$ are total densities of states for both spin projections, and $\rho_{\eta\lambda}(0)$, is the density of states at $K$ valley $\eta=\pm$ for helicity $\lambda=\pm$.

Appendix F Integrals in the analytical model

The integrals in Eq. (II.2) are

[TABLE]

[TABLE]

in the limit of many momenta $\bm{q}$ at the Fermi surface for the $K$ pockets. At energies $|E|<|\Delta|$ much smaller than the normal metal bandwidth, we assume that the density of states varies slowly, such that the following approximation holds $\frac{1}{\Omega}\sum_{\bm{q}}=\frac{\rho_{\eta\lambda}(0)}{2\pi}\int_{0}^{2\pi}d\vartheta\int d\xi$, with $\rho_{\eta\lambda}(0)$ the density of states at the Fermi level for the Fermi surface at valley $\eta$ and helicity $\lambda$. In this approximation the momenta are linearized at the Fermi surface, $\bm{q}\cdot\bm{r}_{ij}=q(\xi)|\bm{r}_{ij}|\cos\vartheta$, with

[TABLE]

The last approximation assumes a parabolic band instead of the warped surface seen in experiments or ab initio simulations. In order to make our calculations realistic, we extract an effective (parabolic band) Fermi momentum $k_{F}^{K}$ and effective velocity $v_{F}^{K}$ from the tight-binding model by performing an angular average over the warped Fermi surface (see Appendix E).

Under these approximations, the $I_{0}$ integral is well-behaved at the Fermi surface and readily yields the result:

[TABLE]

with $J_{0}$, the zero order Bessel function and $H_{0}$, the zero order Struve function. Their argument is

[TABLE]

with superconducting coherence length $\xi_{E}^{K}=\hbar v_{F}^{K}/\sqrt{\Delta^{2}-E^{2}}$.

The $I_{1}$ integral has logarithmic divergences at high energies, which are resolved by introducing a regularization, which ensures convergence of the integral at energies higher than Debye energy $\omega_{D}$ Pientka et al. (2013):

[TABLE]

In the limit of dilute impurities $k_{F}d\gg 1$, the asymptotic expansion of Bessel and Struve functions in $I_{0,1}$ yields:

[TABLE]

The integrals at the $\Gamma$ point are solved identically, and the results follow from above through the replacements $k_{\eta\lambda}\to k_{F}^{\Gamma}$, $v_{F}^{K}\to v_{F}^{\Gamma}$, $\xi_{E}^{K}\to\xi_{E}^{\Gamma}$, and $\rho_{\eta\lambda}(0)\to\rho_{\Gamma}(0)/2$.

Appendix G RKKY Hamiltonian and competing spin-spin interactions

In this Appendix, we analyze in the parabolic approximation of Sec. II the RKKY Hamiltonian which is responsible for the ordering of impurity spins on the metallic \ceNbSe2. To obtain the Hamiltonian, we will work in the limit of dilute impurities $k_{F}r\gg 1$ and we assume that the ordering of impurity spins is set prior to \ceNbSe2 becoming superconducting. Our analysis follows similar works Heimes et al. (2015) by using Matsubara Green’s function formalism to single out the dominant contribution from conduction electrons near the Fermi surface.

Let us assume each impurity coupled to the surface through a Hamiltonian:

[TABLE]

where in contrast to Eq. (12), there is no preference for a particular ordering of spins $\bm{S}_{j}$. Since $\bm{\sigma}$ are the Pauli spin matrices, and $\bm{S}_{j}$ are classical spins at position $\bm{r}_{j}$, the coupling constant $J$ has units (length)${}^{2}/$(time). The spins order through indirect interactions mediated by the conduction electrons, i.e. the RKKY mechanism.

A second order perturbation theory yields the RKKY Hamiltonian, Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957) which in the zero temperature limit reads:

[TABLE]

The prefactor $1/2$ is added to avoid double counting $(jl)$ pairs. The Matsubara Green’s function $g_{jl}(i\omega)$ for normal-state metallic \ceNbSe2 is evaluated at the Fermi surface where it can be decomposed into separate contributions from $\Gamma$ and $K$ electrons. In direct parallel to Eqs. (5-9) in the main text, we find the Green’s functions:

[TABLE]

The integrals over momenta are performed formally following the same route as in the previous appendix. The summation of momenta in the limit of many $\bm{q}$ becomes an integral. Focusing on an small energy interval near the Fermi surface, the integrals over momenta may be approximated using the constant density of states at the Fermi surface, $\frac{1}{\Omega}\sum_{\bm{q}}\to\frac{\rho(0)}{2\pi}\int d\vartheta\int d\xi$. To simplify the ensuing expressions, we assume only in this appendix that the densities of states for different helicities are identical $\rho_{\eta\lambda}(0)\simeq\rho_{\eta,-\lambda}(0)\simeq\rho_{K}(0)/2$, with $\rho_{K}(0)\simeq 0.78/$\mathrm{e}\mathrm{V}$\Omega_{0}$ the total density of states at a given valley. This approximation is reasonable since $\rho_{++}$ is only $5\%$ larger than $\rho_{+-}$ [see Eq. (E)]. The qualitative effects due to spin-orbit coupling enter through the oscillating exponential factors in the Green’s functions. Under these conditions, the angular integrals follow immediately, yielding a first kind zero order Bessel function $J_{0}(q(\xi)r)$. In the limit of dilute impurities $k_{F}r\gg 1$, we use the asymptotic expansion of the Bessel function for large argument:

[TABLE]

The prefactor in the asymptotic expansion of $J_{0}$ is evaluated at Fermi surface, neglecting also a small spin-orbit coupling $k_{F}\gg k_{so}$. Then performing the last integral over $\xi$ yields the following Matsubara Green’s functions at $K$ and $\Gamma$:

[TABLE]

We have used the same notations as before for Fermi velocities, Fermi momenta, densities of states, and $k_{\eta\lambda}=k_{F}^{K}-\eta\lambda k_{so}$. Note that $g_{jl}^{\Gamma}(i\omega)$ and $I_{jl}^{\eta\lambda}(i\omega)$ depend only on the absolute distance $r_{jl}=|\bm{r}_{jl}|$, so they are even under permutation of real-space indices.

Let us return into the RKKY Hamiltonian (54), using the Green’s function decomposition at the Fermi surface. The ensuing Hamiltonian contains three types of contributions. The interaction between the spin at $j$ and $l$ is mediated entirely through conduction electrons at $\Gamma$ valley, or at $K$ valley, or through both $\Gamma$ and $K$ valleys. Therefore the Hamiltonian comprises three possible integrals over Green’s functions:

[TABLE]

Subsequently, the trace over spins is performed using the identities:

[TABLE]

Combining Eqs. (G) and (G) in the RKKY Hamiltonian (54) yields:

[TABLE]

The first term in the $H_{\rm RKKY}$ describes the indirect exchange interaction between spins mediated by electrons in the $\Gamma$ valley. Since in the parabolic approximation the electrons at $\Gamma$ valley are free, the first term recovers the asymptotic form $k_{F}r\gg 1$ of a classical result for the two-dimensional electron gas. Fischer and Klein (1975); Béal-Monod (1987); Litvinov and Dugaev (1998) The second term contains the mixed contribution to the spin susceptibility from both $K$ and $\Gamma$ electrons. The third term represents the contribution from $K$ valley alone and has the same type of spin-spin interactions as discussed for \ceMoS2. Parhizgar et al. (2013)

The Hamiltonian reads after grouping the terms by the spin-spin interactions and summing over $\lambda$ and $\eta$:

[TABLE]

The final Hamiltonian contain three types of spin-spin interactions: Heisenberg $\bm{S}_{j}\cdot\bm{S}_{l}$, Ising $S_{j}^{z}S_{l}^{z}$, and Dzyaloshinskii-–Moriya (DM) $(\bm{S}_{j}\times\bm{S}_{l})_{z}$ terms. Some conclusions may be drawn immediately from the classical Hamiltonian (62). The first two terms favor collinear ordering of spins, while the DM term favors orders where spins rotate in the $xy$ plane. The classical Hamiltonian also shows that the spin-spin interactions decay with square distance between the spins, as in the two-dimensional electron gas, but, in contrast, the interactions are anisotropic in space due to the $K$ valley electrons from \ceNbSe2, similar to \ceMoS2. Parhizgar et al. (2013) The oscillations are roughly with a period $r=\pi/k_{F}$. As seen in Fig. 8, the interactions obey the underlying $C_{3}$ symmetry in \ceNbSe2.

The DM and Ising terms exist only in presence of Ising spin-orbit coupling. Only the spin-rotational invariant Heisenberg interactions survive, when $k_{so}\to 0$ in Eq. (62). Due to vanishing spin-orbit splitting on the conserved in-plane reflection symmetry axes (see App. B), the DM and Ising terms also vanish along them (see Fig. 8).

In general the three terms in $H_{\rm RKKY}$ compete to set the ground state’s ordering of spins. Let us start with the case of two classical spins to get a sense of the competing orders in our system. The numerical results are displayed in Fig. 8 showing that the amplitudes of the three terms are on the same order of magnitude (especially going to even larger distances that shown here). It is immediately apparent that interactions obey the $C_{3}$ symmetry of \ceNbSe2 and that Ising and DM interactions vanish on the conserved in-plane reflection symmetry axes. For a particular choice of distance between spins and their position on the \ceNbSe2 surface, there are situations where DM term dominates, so in those cases the ansatz in Eq. (12) would not be correct. It is interesting to note that for the case $r=3a$ and $\phi=0$, which is the focus of our numerical study in Sec. IV, the DM and Ising interactions almost vanish, leaving a ferromagnetic Heisenberg interaction $J^{\rm H}>0$ in accordance to the ansatz of Eq. (12).

An exhaustive study of all possible spin orderings is outside the scope of the present paper. As we have seen the RKKY interaction is direction dependent, so the ground state spin orientation will depend both on the distance between spins and the chain orientation. In the following, we will minimize the RKKY Hamiltonian starting from a classical ansatz for spin ordering: a transverse conical magnetic structure which allows to investigate the competition between the Heisenberg and the DM interactions. Note that cycloidal, with spins rotating in a plane normal to the \ceNbSe2 plane, or helical ordering, with spins rotating in a plane normal to the impurity chain axis, will result in vanishing DM interactions, since $(\bm{S}_{j}\times\bm{S}_{l})_{z}=0$ in both cases.

The spins precess around the $\hat{\bm{z}}$ axis with a polar angle $\varphi$, at a fixed azimuthal angle $\theta$. The transverse conical structure is generated by spins:

[TABLE]

where $\hat{\bm{r}}$ points in $xy$ plane along the impurity chain. Initial numerical study shows in this case a preference for $\theta=\pi/2$ orientation. Therefore we simplify the problem by fixing $\theta=\pi/2$ and finding the minimum energy under a variation of $\varphi$ in the interval $[0,\pi]$. Note that ferro- and antiferromagnetic order is realized for $\varphi=0$, and, respectively, $\varphi=\pi$, with in-plane $xy$ helical ordering for intermediate values. The RKKY Hamiltonian becomes in this case, for $N\to\infty$ impurities:

[TABLE]

The sum is over distances $jd$ between the spins, with terms decaying as $j^{-2}$. Therefore in the minimization process it is sufficient to truncate the sum to a finite number of terms (here 50). The results are shown in Fig. 9. There are large areas in parameter space where collinear order is preferred. Note that, in contrast, a two-dimensional metal with Rashba spin-orbit coupling does not stabilize collinear order of spin. Instead additional crystal field effects, which at first order favor an Ising coupling of spins, are invoked to produce a collinear order. Heimes et al. (2015) More dedicated numerical studies are required to include these effects, to decide whether they are enough to stabilize ferromagnetic order for any chain orientation, and whether the classical ground state is stable to quantum and thermal fluctuations Heimes et al. (2015); Kim et al. (2014).

Appendix H Effect of superconducting coherence length on the topological phase diagram

This appendix investigates the behavior of the phase diagram as the coherence length is gradually decreased. The coherence length controls the number of distant-neighbor hopping terms relevant in the Hamiltonian. Reducing the coherence length, or equivalently increasing $\Delta$, simplifies the phase diagram, which is seen to gradually losing its fine structure as Hamiltonian matrix elements are less oscillating in momentum space [Figs. 2, 10(a) and 10(b), where $\Delta$ is, respectively, $1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, $10\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, and $100\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$]. The topological phases tend to grow at smaller coherence length. A particular case is in the unphysical short coherence length limit $\xi_{0}\ll d$, such that the effective Hamiltonian has only nearest-neighbor hoppings. The Bloch Hamiltonian follows directly from Eq. (29) under the approximation of polylogarithmic functions with their argument. This leads at $\alpha=1$ or $JS=2/\pi\rho(0)$ to a topological phase extension to the entire allowed [where $h(k)$ is not vanishing] parameter space in $d$ and $\phi$, since $h(k)$ changes sign between $k=0$ and $k=\pi/d$ due to a $\cos(kd)$ behavior:

[TABLE]

Additionally, Fig. 10(c) shows the topological phase diagram for larger distances between magnetic impurities, on the order expected for molecular chains of magnetic porphyrins and phthalocyanines.