Magnetic field driven nodal topological superconductivity in monolayer transition metal dichalcogenides

Recently, Ising superconductors that possess in-plane upper critical fields Hc2 much larger than the Pauli limit field are under intense experimental study. Many monolayer or few layer transition metal dichalcogenides are shown to be Ising superconductors. Here we show that in a wide range of experimentally accessible regimes where the in-plane magnetic field is higher than the Pauli limit field but lower than Hc2, a 2H-structure monolayer NbSe2 or similarly TaS2 becomes a nodal topological superconductor. The bulk nodal points appear on the Γ−M lines of the Brillouin zone where the Ising SOC vanishes. The nodal points are connected by Majorana flat bands, and the flat bands are associated with a large number of Majorana zero energy edge modes that induce spin-triplet Cooper pairs. This work demonstrates an experimentally feasible way to realize Majorana fermions in nodal topological superconductor, without any fine-tuning of experimental parameters. There are many different types of superconducting phases each with unique properties and mechanisms behind their superconductivity. The authors investigate a type of superconductor called an Ising superconductor and demonstrate that by application of a magnetic field they can be driven into a nodal superconducting phase.

A monolayer transition metal dichalcogenide (TMD) is formed by a layer of transition metal atoms sandwiched between two layers of chalcogen atoms in the trigonal prismatic structure 1 . Because of its two-dimensional nature, high electron mobility, and massive Dirac energy spectrum, TMD materials are considered as candidates for next-generation transistors [2][3][4][5][6] . Moreover, since the crystal structure of TMD has no inversion symmetry but respects the out-of-plane mirror symmetry, the strong atomic spin−orbit coupling (SOC) field in TMD materials gives rise to a special type of SOC that we referred to as Ising SOC [7][8][9][10][11][12] . In sharp contrast to Rashba SOC that pins electrons spins to in-plane directions 13,14 , the Ising SOC acts as an effective Zeeman field that pins electron spins to opposite outof-plane directions for electrons with opposite momentum. As TMD materials possess strong Ising SOC, the novel optical and transport phenomena related to the Ising SOC are under intense studies [15][16][17] .
Recently, several TMD materials such as MoS 2 , MoSe 2 , WS 2 , NbSe 2 are found to be superconducting [18][19][20] even down to atomic thickness. More importantly, it was shown that monolayer superconducting MoS 2 7,8 and NbSe 2 9,10 possess extremely high in-plane upper critical field H c2 , more than six times higher than the Pauli limit. It was explained that the enhancement of H c2 is due to the presence of the strong Ising SOC [7][8][9][10] . The Ising SOC pins electron spins of the Cooper pairs to the out-of-plane directions, preventing the electron spins from being aligned with external in-plane magnetic fields and hence protects the superconductivity. While it is intriguing that Ising SOC can enhance the upper critical fields of superconductors, it is not yet understood whether Ising SOC can induce novel superconducting phases.
In this work, we show that, due to the strong Ising SOC and the fact that H c2 of superconducting TMDs is strongly enhanced, a superconducting monolayer TMD can be easily driven to a nodal topological phase by an in-plane magnetic field. Using monolayer NbSe 2 as an example, we demonstrate that the nodal topological phase occupies a large part of the phase diagram, and the predicted nodal topological regime has been achieved in recent experiments 9,21,22 . In the nodal topological phase, bulk nodal points located on Γ−M lines are created by an in-plane magnetic field and these nodal points are connected by Majorana flat bands. We further show that the Majorana flat bands can be detected by tunneling experiments. The realization of Majorana fermions in monolayer NbSe 2 is experimentally more accessible than other current schemes in semiconductor nanowire/superconductor heterostructures 23 and quantum anomalous Hall/ superconductor heterostructures 24 , which have provided evidence of Majorana fermions in 1D and 2D fully gapped topological superconductors [25][26][27][28] . This is due to the fact that driving the system to the nodal topological phase by an external magnetic field does not require the engineering of heterostructures nor the fine-tuning of experimental parameters. Furthermore, our study applies to other TMDs such as monolayer superconducting TaS 2 29,30 .

Results
Band structure in normal state. A monolayer NbSe 2 is formed by a layer of Nb atoms with triangular lattice sandwiched by two layers of Se atoms with triangular lattices. The material has a hexagonal lattice structure when viewed from the out-of-plane direction but with broken A−B sublattice symmetry as shown in Fig. 1a. Therefore, an in-plane mirror symmetry along the ydirection is broken and this gives rise to the Ising SOC that pins electron spins to the out-of-plane directions and splits the energy bands 31 . The band structure of a monolayer NbSe 2 is obtained through first-principle calculations taking into account SOC using ABINIT package (http://www.abinit.org). The results are presented in Fig. 1b. As expected, the shape of the band structure is almost identical to the band structures of monolayer MoS 2 , MoSe 2 , WSe 2 and WTe 2 found in previous works 32 . However, unlike Mo-and W-based materials which are insulating intrinsically, an Nb atom has one less d-electron in the outer most shell than Mo and W atoms and the chemical potential of NbSe 2 lies in the valence band. The band splitting at the K points due to Ising SOC is about 150 meV. Importantly, the band splitting is significant at the Fermi energy even though the Fermi surface is far away from the K points. As we show in the next section, this Ising SOC at the Fermi energy plays a crucial role in protecting the superconductivity from the paramagnetic effects of in-plane magnetic fields.
In order to study the superconducting properties of the material, we construct a six-band tight-binding model that consists of the d z 2 , d xy and d x 2 À y 2 orbitals 32,33 of the Nb atoms. Including the third-nearest-neighbor hoppings, the band structure of the DFT calculations can be well explained. The details of the tight-binding model H N (k) can be found in Supplementary Note 1 with the fitting parameters shown in Supplementary  Table 1 Fig. 1 The crystal structure and electronic bands for monolayer NbSe 2 . a Lattice structure of monolayer NbSe 2 (left) and its top view (right). M x denotes the in-plane mirror symmetry. b The band structure for the normal state of monolayer NbSe 2 from DFT calculations. c The energy contour at the Fermi level from tight-binding model important to note that the Fermi surfaces of NbSe 2 as depicted in Fig. 1c have been observed through angle-resolved photoemission spectroscopy 12,34 .
Ising superconductivity in monolayer NbSe 2 . It is well-known that bulk NbSe 2 is superconducting with T c of 7 K 35 . However, monolayer NbSe 2 was shown to be superconducting only very recently with T c on-set at about 5 K and zero resistance appears at about 3 K 9,21,34,36 . In bulk NbSe 2 , magnetic field can create vortices that eventually destroy superconductivity with H c2 of about 14.5 T 37 . In the monolayer case, the orbital effects of in-plane magnetic fields are completely suppressed and the magnetic field can suppress superconductivity only through Zeeman coupling with electron spins.
Similar to the gated superconducting MoS 2 , NbSe 2 was shown to have strongly enhanced H c2 at 35 T even though the T c is much lower than the bulk T c . While the experimental data in MoS 2 can be well explained by solving the self-consistence gap equations by taking into account the Ising SOC and Rashba SOC induced by gating, a microscopic theory for explaining the H c2 data of monolayer NbSe 2 is still lacking due to the more complicated band structures in the valence bands. In this work, with the tightbinding model introduced in the method section, we can take into account all the three Fermi pockets shown in Fig. 1c to explain the experimental data. An intra-orbital attractive interaction U between the electrons is introduced. U is determined by the selfconsistent gap equation at zero magnetic field. Here, Δ = 1.76 k B T c is the pairing gap at zero temperature, V is the area of the sample and c k ,↑/↓,o is the electron annihilation operator for orbital o.
After finding U at zero magnetic field, we can determine Δ as a function of magnetic field by minimizing the free energy density with E k,n the eigenvalue adopted from the Bogoliubov-de Gennes Hamiltonian H BDG as shown in detail in Supplementary Note 2. The phase diagram for this Ising superconductor is determined and shown in Fig. 2, with H P ≈ 1.84T c the Pauli limit 38 . As shown in Fig. 2, the calculated H c2 (solid line) is strongly enhanced. The theoretical values of H c2 , using the parameters from the tightbinding model without any tuning parameters, are higher than the experimental values (the red dots). Such discrepancy can be due to the disorder 39 and the fact that the NbSe 2 sample has ripples, and the ripples introduce Rashba-type SOC 40 which competes with Ising SOC. This Rashba-type SOC can lower the in-plane H c2 as in the case of gated superconducting MoS 2 7 , and a tiny small Rashba SOC is enough to match the in-plane H c2 with the experimental data, shown as the red line in Fig. 2. We ignore the effects of ripples in this work. To highlight the importance of Ising SOC, the H c2 in the absence of Ising SOC is denoted by the dashed line in Fig. 2. It is evident that the Ising SOC enhances H c2 so that it becomes much higher than the Pauli limit.
Nodal topological phase. From Fig. 2, one can see that superconductivity survives in the regime where the applied magnetic field is higher than the Pauli limit field. The question is: what are the properties of the superconducting phase under strong magnetic fields? To answer this question, the spectral function of the a semi-infinite NbSe 2 stripe is calculated. Here, E is energy and R is a point on the armchair edge and G is the Green's function of the BdG Hamiltonian defined in Supplementary Note 2. The armchair edges are parallel to the y-direction and subject to periodic boundary conditions. The in-plane magnetic field applied is chosen to be H = 3H P . In the bulk Brillouin zone, it can be seen that there are six pairs of nodal points lying in the Γ−M line, as is shown in Fig. 3a. After projecting onto the edge Brillouin zone along k y , there are four pairs of projected nodal points. The outer two pairs of projected nodal points are connected by two sections of Majorana flat bands respectively. The inner two pairs of projected nodal points are from the middle four pairs of nodal points in Fig. 3a, so they are connected by two sections of doubly degenerated Majorana flat bands. In the edge Brillouin zone, the four projected nodal points as well as the four sections of Majorana flat bands are shown in Fig. 3b. Other finite energy modes emerging below the bulk continuum spectrum in Fig. 3b are the trivial edge modes, which can be removed by local potentials.
To understand this nodal topological phase, we use an effective Hamiltonian approach and understand the system from a symmetry point of view. Before writing down the effective Hamiltonian near the Fermi energy, we note that the Ising SOC near the K points is very strong as shown in Fig. 1b. Therefore, the superconducting states near K points are hardly affected by the in-plane magnetic field at H ¼ 3H p ( H c2 . Therefore, to understand the gap closing effect of the in-plane magnetic field, we only need to focus on the states in the Γ pocket where the Ising SOC is weaker.
Around the Γ pocket, the d z 2 orbital dominates 32,33 . In the basis of [c k,↑ , c k,↓ ] and in the absence of external magnetic fields, the Hamiltonian has to satisfy the point group symmetries M z = −iσ z , M x = −iσ x , C 3 ¼ e Àiðπ=3Þσ z and time-reversal symmetry T = iσ y K. These symmetries restrict the effective Hamiltonian, up to third order in k, to have the form: Here, μ is the chemical potential and k ± = k x ± ik y . It is important to note that λ SOC term vanishes along the Γ − M lines as dictated by symmetries, consistent with the calculations in Fig. 1c.
Including the in-plane magnetic field, spin-singlet pairing Δ, and in the basis of c k;" ; c k;# ; c y Àk;" ; c y Àk;# h i , the Hamiltonian in the  Fig. 2 The magnetic field-temperature phase diagram for the superconducting monolayer NbSe 2 . The red dots represent the experimental data from ref. 9 . The calculated H c2 is denoted by the black solid line and it is much higher than the Pauli limit which is denoted by the dashed line. When the applied in-plane magnetic field is higher than the Pauli limit but lower than H c2 , the system can go through a phase transition from a fully gapped superconducting phase to a nodal topological phase. The sample ripples-induced Rashba spin−orbit coupling is further considered in the H c2 calculation denoted by red solid line, which agrees well with the experimental data superconducting phase can be written as: Here, τ i denotes Pauli matrices in the particle-hole basis. H x and H y are magnetic fields in the x-and y-directions respectively, μ B is the Bohr magneton and g is the electron's gyromagnetic ratio. Along the Γ−M line where the Ising SOC vanishes, the electrons' spins are all in-plane polarized by the in-plane magnetic field, so the spin-singlet Cooper pairs from the time-reversal partners are strongly suppressed. Once the Zeeman energy 1 2 μ B gH exceeds Δ the pairing gap closes along the Γ−M line. However, away from the Γ−M line, the Ising SOC pins the time-reversal partners' spin to the opposite direction along z and thus the spin-singlet Cooper pair will condensate so that the pairing gap remains finite. As a result, there will be six pairs of point nodes split out along the Γ −M line as the in-plane magnetic field is gradually applied, as is shown in Fig. 3a.
The six pairs of point nodes are protected by the chiral symmetry. Even though time-reversal is broken by the external magnetic field, H s respects a time-reversal like symmetry U T H s k x ; k y U À1 T ¼ H s Àk x ; k y with U T = M z Tτ z and k y unchanged under the symmetry operation. Moreover, H s respects a 1D particle-hole like symmetry U P H s k x ; k y U À1 P ¼ ÀH s Àk x ; k y with U P = τ x K. Therefore, H s respects the chiral symmetry CH s k x ; k y C À1 ¼ ÀH s k x ; k y where C = τ y σ x . As a result, for any fixed k y , H s is in the BDI class and the Hamiltonian can be topologically nontrivial [41][42][43][44] . For the range of k y where H s is nontrivial, there are Majorana zero energy modes on the edge of the system as shown in Fig. 3b. By tuning k y as a parameter, the system can undergo a topological phase transition from a trivial regime to a nontrivial regime by closing the bulk gap. These topological phase transition k y points are the nodal points in Fig. 3b. The nodal points and the Majorana flat bands shown in Fig. 2 can be reproduced by the simple Hamiltonian H s in Eq. 3.
Detection. In this section, we discuss the experimental detection of the nodal topological superconducting phase in NbSe 2 . As discussed above, superconducting NbSe 2 can be driven from a fully gapped superconducting phase to a nodal superconducting phase by an in-plane magnetic field. As in the case of d x 2 Ày 2 -wave superconductors, we expect the density of states of the system as a function of energy to be V-shape near zero energy when the system is nodal. From the low-energy bulk spectrum of the effective Hamiltonian H s , the density of states near zero energy can be found to be Nð2π E j jÞ=ðv x v y Þ, where E is the energy and v x , v y denote the Fermi velocity in the x,y-directions near the nodal points and N is the number of nodal points. Thus, the bulk density of states near zero energy is indeed linearly proportional to |E|.
Using the real space version of the six-band tight-binding model H BDG defined in Supplementary Note 1, including pairing and the in-plane magnetic field, we calculate the local density of states in the bulk of the sample.
here G R (E) = (E + iη − H BdG ) −1 is the retarded Green's function. The local density of states for a point in the bulk is shown in Fig. 4a. When the applied in-plane field is lower than the Pauli limit, the density of states as a function of energy is U-shape (blue curve) near zero energy indicating a fully gapped superconducting state. On the other hand, the density of states is changed to Vshape (red curve) when the system is driven to the nodal phase as expected. Even though the V-shape density of states in the bulk is only the signature of a nodal phase, not necessarily a nodal topological phase, the observation of the fully gapped to nodal phase transition can be important for identifying the topological phase transition point. Moreover, in the nodal topological phase, there are a large number of zero energy Majorana modes residing on the edges of the sample. These zero energy modes strongly enhance the local density of states on the edge of the superconductor. The local density of states at the edge of the sample, before and after the topological phase transition, are shown in Fig. 4b. It is evident that, in the nodal topological phase with in-plane field stronger than the Pauli limit field, the zero energy density of states are Interestingly, due to the zero energy Majorana modes, the pairing correlations induced on the edge of the superconductor are indeed spin-triplet. To show this, we calculated the pairing correlation for a site on the armchair edge of the sample. The pairing correlation can be parametrized as 11,14,45 : where σ and o are the spin and orbital indices respectively. In the matrix form, the pairing correlation is written as: Here, the so-called d-vector characterizes the triplet pairing correlation and ψ characterizes the spin-singlet pairing correlation. In Fig. 5a, it is shown that the triplet pairing correlation at zero energy on the armchair edge of the sample is nonzero only when the Zeeman energy V z ¼ ð1=2Þμ B gjHj is larger than Δ. Figure 5b shows the spatial variation of |d| across the sample for a strip of NbSe 2 which has a width of 200-sites in the x-direction (perpendicular to the armchair edge). It is evident that |d| is significant only near the two edges of the sample. On the other hand, the singlet pairing correlation amplitude ψ at zero energy is negligible. Furthermore, we calculated s = i(d × d * )/|d| 2 which gives the spin-polarization direction of the Cooper pairs 46 as shown in Fig. 5c. This indicates that the Cooper pairs injected into the superconductor by the normal lead are spin-polarized with spin pointing to s-direction 44,47 . As a result, even though the lead attached to the superconductor is nonmagnetic, the current from the lead to the superconductor is spin-polarized. In other words, the Majorana modes can act as spin filters and only allow electrons with spin pointing to s-direction to tunnel into the superconductor. We believe this property of monolayer NbSe 2 may lead to applications in superconducting spintronics 48,49 .

Discussion
We would like to discuss some important points about the nodal topological phase mentioned above. First, the large number of Majorana zero energy modes associated with the Majorana flat band in Fig. 2b are protected against on-site disorder. This is in sharp contrast to the d x 2 Ày 2 -wave superconductor in which the zero energy fermionic modes can be lifted to finite energy by disorder. One can show that the zero energy Majorana modes are not protected by the bulk gap (since the system is nodal) but by the chiral symmetry C = τ y σ x of H s in Eq. 3 [50][51][52][53][54][55] . In the presence of the chiral symmetry, the band topology can be characterized by net chirality number, which is defined as the difference between number of states with chirality +1 and −1 50 . For the finite energy modes, states with opposite chirality number ±1 always come in pairs so the net chirality is zero. In contrast, Majorana zero modes along the same edge have a nonzero net chirality number seen in Supplementary Note 3. Since the on-site disorder respects the chiral symmetry, perturbations due to the on-site disorder can only excite pairs of zero energy states with chirality number ±1 to finite energy modes. Thus, a finite number of Majorana edge modes always exist as long as the bulk topology is not changed. This is similar to the toy model studied in ref. 43 , where an in-plane magnetic field creates Majorana flat bands in a time-reversal invariant p ± ip topological superconductor.
Second, we emphasize that the Γ pocket is essential for creating the nodal topological phase. In the case of gated MoS 2 or WS 2 , superconductivity appears in the conduction bands near the K pockets only. Therefore, the nodal phase discussed above cannot be found in those superconductors. On the other hand, intrinsic superconductors NbSe 2 , NbS 2 , TaSe 2 and TaS 2 are candidate materials for realizing the nodal topological superconducting phase found in this work.
Third, the armchair edge is not essential for the observation of the Majorana flat bands as long as the edge is not parallel to the zig-zag edge (k x direction in Fig. 2a). On the zig-zag edge, the projection of the bulk nodal points can cancel each other such that there are no topologically nontrivial regimes for all k x parallel to the zig-zag edge. For all other edges, one can find finite sections of Majorana flat bands along the edges when the system is driven to the nodal topological phase.

Methods
Tight-binding Hamiltonians. In the basis of c k;d z 2 ;" ; c k;d xy ;" ; c k;d x 2 Ày 2 ;" ; c k;d z 2 ;# ; c k;d xy ;# ; c k;d x 2 Ày 2 ;# h i , the six-band model used has the form 32 with where I n means the n × n identity matrix. The details about these matrix elements can be found in Supplementary Note 1. Then the Bogliubov-de Gennes The real space version of H TNN and H BdG can be found in Supplementary Note 2.
Data availability. The data supporting the findings of this study are available within the paper, its Supplementary Information files, and from the corresponding author upon reasonable request.