Magnetic Field Driven Nodal Topological Superconductivity in Monolayer Transition Metal Dichalcogenides

Recently, Ising superconductors which possess in-plane upper critical fields 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. In this work, 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 $H_{c2}$, a 2H-structure monolayer NbSe$_2$ or simiarly TaS$_2$ becomes a nodal topological superconductor. The bulk nodal points appear on the $\Gamma- M$ lines of the Brillouin zone where the Ising SOC vanishes. The nodal points are connected by Majorana flat bands, similar to the Weyl points being connected by surface Fermi arcs in Weyl semimetals. The Majorana flat bands are associated with a large number of zero energy Majorana fermion edge modes which induce spin-triplet Cooper pairs. This work demonstrates an experimentally feasible way to realise Majorana fermions in nodal topological superconductor, without any fining tuning of experimental parameters.


I. INTRODUCTION
In recent experiments, a new type of superconductors called Ising superconductors were discovered. It was shown that monolayer superconducting transition metal dichalcogenides (TMD) such as MoS 2 [1,2] and NbSe 2 [3] 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 a special type of SOC called Ising SOC [1,3,4]. Ising SOC is caused by the lattice structure of monolayer TMD which breaks an in-plane mirror symmetry and it pins electrons spins to the out-of-plane directions. This is in sharp contrast to the Rashba SOC which pins electron spins to in-plane directions [5,6]. While it is intriguing that Ising SOC can enhance the upper critical fields of superconductors, it is not 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 the in-plane H c2 of superconducting TMDs is strongly enhanced, a superconducting monolayer TMD can be easily driven to the nodal topological phase by an in-plane magnetic field. As shown in Fig.1d, using monolayer NbSe 2 as an example, that the nodal topological phase occupies a large part of the phase diagram in Fig.1d and the predicted nodal topological regime has been achieved in recent experiments [3,7,8]. As shown in Fig.2a, nodal points on Γ-M lines are created by an in-plane magnetic field and these nodal points are connected by Majorana flat bands as shown in Fig.2b. We further show that the Majorana flat bands can be detected by tunnelling experiments as illustrated in The H-T phase diagram of the monolayer NbSe2. The red dots represent the experimental data from Ref. [3]. The calculated Hc2 is denoted by the 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 Hc2, the system can go through a phase transition from a fully gapped superconducting phase to a nodal topological phase.
on semiconductor nanowire/superconductor heterostructures [9] and quantum anomalous Hall/superconductor heterostructures [10] have provided evidence of Majorana fermions in 1D and 2D fully gapped topological superconductors [11][12][13][14]. In these systems, fine tuning of experimental parameters is needed to achieve the topological regime as the topological regimes are extremely narrow. This work demonstrates an experimentally feasible way to realise Majorana fermions in nodal topological superconductor, without engineering heterostructures and without any fining tuning of experimental parameters.
Moreover, the appearance of the nodal topological phase predicted does not depend on the detailed band structure of the system. It only depends on the fact that the Fermi surfaces near the Γ point of the superconductors are split by Ising SOC. Our study applies to a whole class of TMDs such as monolayer superconducting TaS 2 [15,16] and NbSe 2 [3] which are under intense experimental studies recently due to the coexistence of superconductivity with charge density way and strong Ising SOC.
In the following sections, we first present a realistic band structure of a monolayer NbSe2. Secondly, using the tight-binding model, we explain the origin of the enhanced in-plane H c2 in NbSe 2 . Thirdly, we discuss how the nodal topological phase can be induced by an inplane magnetic field resulting in nodal points which are connected by Majorana flat bands. Finally, we study how to detect this topological nodal phase by tunnelling experiments.

A. 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 y direction is broken and this gives rise to the Ising SOC which pins electron spins to the out-ofplane directions and split the energy bands [17]. The band structure of a monolayer NbSe 2 is obtained through first-principle calculations taking into account SOC using ABINIT package [18]. 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 [2]. However, unlike Mo and W based materials which are insulating intrinsically, a 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 150meV. Importantly, the band splitting is significant at the Fermi energy even 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 which consists of the d z 2 , d xy and d x 2 −y 2 orbitals [1,2] of the Nb atoms. Including the third-nearest-neighbor hoppings, this band structure of the DFT calculations can be well explained. The details of the tight-binding model H N (k) and the fitting parameters can be found in the Appendix and the Supplementary Material [21]. The Fermi surfaces from the tight-binding model are shown in Fig.1c. At the Fermi level, there are Fermi pockets around Γ, K and -K points. In general, the bands are split by Ising SOC except for the states lying along the Γ − M lines. It is important to note that the Fermi surfaces of NbSe 2 as depicted in Fig.1c have been observed through angle-resolved photoemission spectroscopy [4,22].

B. Ising superconductivity in monolayer NbSe2
It is well-known that bulk NbSe 2 is superconducting with T c of 7K [24]. However, monolayer NbSe 2 was shown to be superconducting only very recently with T c onset at about 5K and zero resistance appears at about 3K [3,7,22,23]. In bulk NbSe 2 , magnetic field can create vortices which eventually destroy superconductivity with H c2 of about 14.5T [25]. 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 35T 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 tight-binding model introduced in the Appendix, 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 [21] between the electrons is introduced. U is determined by the self-consistent gap equation ∆ = U/V k,o c k,↓,o c −k,↑,o at zero magnetic field. Here, ∆ = 1.76k 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 Bogoliubovde Gennes Hamiltonian H BDG as shown in detail in the Supplementary Material [21]. The phase diagram for this Ising superconductor is determined and shown in Fig.1(d), with H p ≈ 1.84T c the Pauli limit [26]. As shown in Fig.1(d), the calculated H c2 (solid line) is strongly enhanced. The theoretical values of H c2 , using the parameters from the tight-binding model without any tuning parameters, are higher than the experimental values (the red dots). Such discrepancy can be due to the fact that the NbSe 2 sample has ripples [27], and the ripples introduce Rashba type SOC [28] 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 [1]. 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.1(d). It is evident that the Ising SOC enhances H c2 so that it becomes much higher than the Pauli limit.

C. Nodal topological phase
From Fig.1d, 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 A(E, R) = − 1 π Tr [ImG (E, R)] of the a semi-infinite NbSe 2 stripe is shown in Fig.2b. 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 Material [21]. 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 . It can be seen that there are four sections of Majorana flat bands in Fig.2a. The Majorana flat bands connect the four pairs of bulk nodal points. These four pairs of nodal points are projections of the six pairs of bulk nodal points in the bulk energy spectrum as shown in Fig.2a.
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 Ising SOC is weaker.
Around the Γ pocket, the d z 2 orbital dominates [1,2]. In the basis of [c k,↑ , c k,↓ ] and in the absence of external magnetic fields, the Hamiltonian has to satisfy the point 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 † −k,↑ , c † −k,↓ ], the Hamiltonian in the 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 ydirections respectively, µ B is the Bohr magneton and g is the electron's gyromagnetic ratio. The pairing potential is denoted by ∆. 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, for fixed k y , H s respects the chiral symmetry CH ky (k x )C −1 = −H ky (k x ) where C = σ x τ y . As a result, for any fixed k y , H s is in the BDI class and the Hamiltonian can be topologically non-trivial [29][30][31][32]. For the range of k y where H s is non-trivial, there are Majorana zero energy modes on the edge of the system as shown in Fig.2b. By tuning k y as a parameter, the system can undergo a topological phase transition from a trivial regime to a non-trivial regime by closing the bulk gap. These topological phase transition k y points are the nodal points in Fig.2b. The nodal points and the Majorana flat bands shown in Fig.2 can be reproduced by the simple Hamiltonian H s in Eq.3.

D. 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| vxvy , where E is the energy and v x , v y denote the Fermi velocity in the x, ydirections 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 tightbinding model H BDG defined in the Appendix and the Supplementary Material [21], 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.3a. 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.3b. 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 strongly enhanced (red curve). Both the bulk density of states and the edge density of states can be detected by tunnelling spectroscopy experiments.
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 [6,33,34]: 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 tripletpairing correlation and ψ characterizes the spin-singlet pairing correlation. In Fig.4a, it is shown that the triplet pairing correlation at zero energy on the armchair edge of the sample is non-zero only when the Zeeman energy V z = 1 2 µ B g|H| is larger than ∆. Fig.4b shows the spatial variation of |d| across the sample for 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 [35] as shown in Fig.4c. This indicates that the Cooper pairs injected into the superconductor by the normal lead are spin-polarized with spin pointing to s-direction [32,36]. As a result, even though the lead attached to the superconductor is non-magnetic, the current from the lead to the superconductor is spinpolarized. 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 [37,38].

III. 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 disorder. This is in sharp contrast to the d x 2 −y 2wave 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 is not protected by the bulk gap (since the system is nodal) but by the chiral symmetry C = σ x τ y of H s in Eq.3 [39]. This chiral symmetry is not broken by on-site disorders. This is similar to the toy model studied in Ref. [31], where an in-plane magnetic field creates Majorana flat bands in a time-reversal invariant p±ip topological superconductor.
Second, we emphasis that the Γ pocket is essential for creating the nodal topological phase. In the case of gated MoS 2 or WS 2 , superconductivity appear 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 non-trivial regimes for all k x parallel to the zigzag 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.

Appendix: Tight-binding Hamiltonians
In the basis of [c k,d z 2 ,↑ , c k,dxy,↑ , c k,d x 2 −y 2 ,↑ , c k,d z 2 ,↓ , c k,dxy,↓ , c k,d x 2 −y 2 ,↓ ], the six-band model used has the form [2] where I n means the n × n identity matrix. The details about these matrix elements can be found in supplementary material [21]. Then the Bogliubov-de Gennes Hamiltonian in Nambu spinor representation can be written as   t0  t1  t2  t11  t12  t22  r0  r1  r2  r11  r12  u0  u1  u2  u11  u12 In this way the Bogliubov-de Gennes Hamiltonian for the superconducting monolayer NbSe 2 ribbon can be obtained as With this Hamiltonian, the spectral function A = Tr [ImG (E, R)] for the semi-infinite NbSe 2 ribbon with armchair edge is calculated in Fig. 2. Here R is the last slab of the NbSe 2 ribbon and the Green's function G = E + i0 + − H BdG ribbon −1 can be obtained through the surface Green function method [S3]. We set ∆ = 0.02eV in