Ising superconductivity and magnetism in NbSe$_2$

Recent studies on superconductivity in NbSe$_2$ have demonstrated a large anisotropy in the superconducting critical field when the material is reduced to a single monolayer. Motivated by this recent discovery, we use density functional theory (DFT) calculations to quantitatively address the superconducting properties of bulk and monolayer NbSe$_2$. We demonstrate that NbSe$_2$ is close to a ferromagnetic instability, and analyze our results in the context of experimental measurements of the spin susceptibility in NbSe$_2$. We show how this magnetic instability, which is pronounced in a single monolayer, can enable sizeable singlet-triplet mixing of the superconducting order parameter, contrary to contemporary considerations of the pairing symmetry in monolayer NbSe$_2$, and discuss approaches as to how this degree of mixing can be addressed quantitatively within our DFT framework. Our calculations also enable a quantitative description of the large anisotropy of the superconducting critical field, using DFT calculations of monolayer NbSe$_2$ in the normal state


I. INTRODUCTION
The transition metal dichalcognides (TMD) exhibit an astonishing variety of phenomena and phase transitions, which includes charge-density waves (CDW) [1], superconductivity [2], and magnetism [3]. Bulk 2H-NbSe 2 , which is one of the canonical transition metal dichalcogenides, exhibits a rich phase diagram, which includes a superconducting and a CDW phase [1]. Superconductivity in bulk NbSe 2 has been studied extensively both experimentally [2,[4][5][6] and theoretically [7,8], and the superconducting transition temperature, T c , has been experimentally identified as ∼ 7 K [5]. While coupling between the superconducting and CDW order parameters is certainly possible, it was found to be weak in NbSe 2 : the superconducting phase remains robust while the CDW phase collapses as a function of increasing pressure [1]. Since NbSe 2 is also a layered van der Waals material, this has inspired several studies of superconductivity in monolayer NbSe 2 [9][10][11][12] and several intriguing proposals that seek to exploit proximity induced effects at interfaces between monolayer NbSe 2 and magnetic materials [13][14][15][16].
Monolayer NbSe 2 , unlike the bulk structure, lacks inversion symmetry, which leads to a large spin-orbit (SO) splitting of the states at the K, and its inversion partner, K ′ , points [17] (there is an additional Fermi surface with states around Γ, which we will also discuss later in this study). The magnitude of the SO-splitting is larger than the superconducting order parameter. The zero-magnetic field T c of monolayer NbSe 2 is ∼ 3 K, [9,10], which is lower than the bulk T c . The combination of SOcoupling and broken inversion symmetry locks the pseudospins near K and K ′ to be parallel to the c-axis of the monolayer. Due to time-reversal symmetry, pseudospins at the K and K ′ points are antiparallel, and their energies are degenerate. Hence, the Cooper pairs that form completely break their rotational invariance in spin space. This leads to a novel phenomenon, aptly named Ising superconductivity. One key consequence of this unique pairing is that the superconducting phase survives in the presence of in-plane magnetic fields that considerably exceeds the Pauli limit [9,10].
Thus far, theoretical analyses of the superconducting pairing mechanism in monolayer NbSe 2 [18,19] have relied on model descriptions of superconductivity in materials that lack inversion symmetry [20][21][22], loosely based on the band structure calculated from first principles [19] However, a quantitative description of superconductivity in monolayer is lacking. There is also a lack of consistency between first-principles descriptions of superconductivity in bulk NbSe 2 and experimental results. State-of-the-art first-principles calculations that are usually very accurate for superconductors where the pairing is entirely due to electron-phonon coupling, overestimate T c in bulk NbSe 2 and isostructural NbS 2 [23,24] by a factor of ∼ 3 and the zero-temperature gap by a factor of ∼ 4. Furthermore, experimental measurements of the spin susceptibility, χ s , in bulk NbSe 2 [25] find a χ s ∼ 3×10 −4 emu/mole, which, as we will show later, considerably exceeds the bare bulk Pauli susceptibility, χ 0 .
Two plausible mechanisms that can be invoked to explain this discrepancy between theory and experiment are the potential role of strong electron-electron interactions and strong spin fluctuations. In Ref. 24, the authors suggested the overestimation of T c in their first-principles calculations can be corrected by accounting for a reduction in the effective mass induced by electron-electron interactions, which they described within the GW approximation. However, reducing the effective mass by a factor of (m * /m) reduces the density of states (DOS) by a factor of (m * /m) and increases the magnitude of the electronphonon matrix element by a factor of (m * /m). Indeed, the DOS is proportional to the one-electron Green function and the electron-phonon matrix element is proportional to the derivative of the inverse Green function. Since the electron-phonon coupling constant depends linearly on the DOS and quadratically on the matrix element, reducing the effective mass increases the strength of the electron-phonon coupling. Hence, strong electronelectron interaction effects alone do not provide a solution to this discrepancy.
The latter mechanism, which is the role of strong spin fluctuations in NbSe 2 , has thus far remained unaddressed. Fluctuations in the magnetic moment and magnetic order have been shown to be a source of pairing, or pair-breaking, of Cooper pairs in a number of other materials [26,27]. Furthermore, strong spin fluctuations can also lead to a sizeable Stoner renormalization of χ 0 . To our knowledge, all theoretical studies of bulk and monolayer NbSe 2 at their equilibrium lattice parameters have found the material to be non-magnetic. However, calculations of monolayer NbSe 2 subject to tensile strain exceeding 2% have predicted a ferromagnetic ground state [28].
This seeming lack of consistency between the various theoretical and experimental results on bulk NbSe 2 reported in the literature indicates there is still a need to explore the fundamental properties of NbSe 2 . For example, if spin fluctuations are operative in NbSe 2 it is unclear how this may impact arguably the most interesting aspect of Ising superconductivity in monolayer NbSe 2 , which is the possibility of a singlet-triplet mixed state. While Ising superconductivity is well understood at the phenomenological level, it has never been described on a quantitative level using first-principles calculations. This also precludes a quantitative description of Ising superconductivity, which has been observed in monolayers of several other transition metal dichalcogenides beyond NbSe 2 . [10,29].
In this paper we will demonstrate that, indeed, bulk and especially monolayer NbSe 2 are close to a magnetic instability. We will "translate" the existing model theory that has been developed to analyze superconductivity in materials that lack inversion symmetry and bands split by SO-interaction [30], such as monolayer NbSe 2 , into DFT parlance, which allows us to develop a quantitative theory of the critical field anisotropy in this material. We use nonrelativistic and relativistic fixed-spin moment calculations to determine the spin susceptibility of bulk and monolayer NbSe 2 . Finally, we will use the insights obtained from our calculations to discuss possible ramifications on the superconducting order parameter in NbSe 2 , in particular, the factors that control the scale of the triplet admixture to the order parameter.

II. RESULTS AND DISCUSSION
A. Electronic structure The bulk unit cell of 2H-NbSe 2 consists of two monolayers of NbSe 2 ( Fig.1(a)) where in a single monolayer the Nb atoms are in a trigonal prismatic coordination with the Se atoms. The Nb atoms in one of the monolayers are vertically above the Se atoms of the second monolayer (cf. Fig. 1(b)), which leads to a center of inversion that is between the two monolayers of the unit cell. The calculated bulk lattice constants, a=3.449Å , and c=12.550Å, are in agreement with experimental reports of the lattice constants of bulk NbSe 2 [31].
The trigonal crystal field splits the 4d-states of Nb 4+ into three different groups: d z 2 , [d x 2 −y 2 , d xy ] and [d xz ,d yz ]. The bulk band structure of NbSe 2 , which has been studied extensively [32,33], has three bands that cross the Fermi level. Two of the bands are derived from Nb d-states and the third band is derived from Se p z states. Spin-orbit interaction leads to a mixing of the Se p and Nb d-states along the Γ-K-M path of the Brillouin zone, but does not lead to SO-splitting. The density of states at the Fermi level, N (E F ) is 2.7 states/eV ·cell, which leads to a bulk Pauli susceptibility, χ 0 = 0.87×10 −4 emu/mole, which is a factor of ∼ 3.5 lower than the experimentally reported spin susceptibility of bulk NbSe 2 [25], which suggests spin fluctuations are operative in NbSe 2 .
For the case of the monolayer, there is one band that crosses the Fermi level several times, leading to three Fermi contours, one contour around the Γ-point and two contours around K and K ′ (these are related by inversion symmetry). At Γ, the band character is Nb d z 2 , with a minor admixture of Se p z . As the band progresses toward K or K ′ , this leads to a larger admixture of Nb d x 2 −y 2 and d xy orbitals in addition to a minor contribution of Se p xy states. The states at the K and K ′ contours are comprised entirely of Nb [d x 2 −y 2 , d xy ] states. In the absence of SOinteraction, this band is spin-degenerate. When we allow for SO-interaction, the lack of a center of inversion in the monolayer leads to SO-splitting everywhere except along the Γ-M line (cf. Fig .2(a)).
To understand why the pseudospin state does not have an in-plane component and why the splitting is small near Γ, it is instructive to rationalize this splitting from the band structure point of view. If we neglect the minor admixture of Se p-states to the bands that cross the Fermi level in monolayer NbSe 2 , a state at a given wave vector k can be defined as follows: where η 2 + β 2 + γ 2 = 1. Note that d z 2 corresponds to |l, m = |2, 0 , d x 2 −y 2 to (|2, 2 + |2, −2 ) √ 2, and d xy to (|2, 2 − |2, −2 )/i √ 2. Accounting for spin, the Hamiltonian at each k point is a (2 × 2) matrix, and, by virtue of the z/ − z mirror symmetry, does not include contributions from the |2, ±1 orbitals. Thus, the nondiagonal matrix elements L ± are zero. However, it is easy to show that the diagonal element L z = 2(η Im β − β Im η). One phase can always be selected as real, for instance, η, then L z = 2η Im β.
In the centrosymmetric bulk 2H-NbSe 2 , β can also be chosen to be real, and there is no SO-induced spinsplitting (but there is splitting due to doubling of the unit cell). In the monolayer, β is complex everywhere except the Γ−M direction, and therefore the diagonal elements of this (2 × 2) matrix have opposite signs, ±λη Im β, where λ measures the strength of the SO-coupling. Consequently, the splitting is small around the Γ pocket (the maximum splitting at this Fermi-surface is ∼70 meV, which occurs where it intersects with the Γ−K and Γ−K ′ lines), where |γ 2 | ≫ |ηβ|, but sizeable (∼ 150 meV) on the K and K ′ contours, where |γ 2 | ≪ |ηβ|. Due to the absence of nondiagonal coupling, the pseudospin-split states are also pure S z = ± 1 2 spin states and the direction of the pseudospin flips between the K and K ′ valleys as illustrated in Fig. 2(b). As we will discuss later, pseudospin states are no longer pure S z states when an external inplane magnetic field is applied.

B. Magnetism without spin-orbit coupling
To quantitatively address the role of spin fluctuations we calculate the spin susceptibility, χ, by considering the effect that a uniform external magnetic field, H, has on the magnetic moment, m, in NbSe 2 , via the Zeeman interaction, where χ = ∂m/∂H. In practice, one way to apply a magnetic field within a first-principles calculation, is to apply a constraint on the magnetization and compute the total energy, E, as a function of the magnetic moment, m. This "fixed-spin moment" (FSM) approach allows us to define χ as χ = (∂ 2 E/∂m 2 ) −1 . Similarly, relativistic FSM calculations that include SO-interaction let us determine the change in total energy and in turn χ for directions parallel to the c-axis ( 001 ) and perpendicular to the c-axis ( 100 ) of NbSe 2 .
The non-relativistic and relativistic (along 001 and 100 ) FSM calculations result in the same qualitative trends; the total energy increases monotonically with respect to the total energy of the non-magnetic state as a function of increasing magnetic moment. As an example, we illustrate the results of our non-relativistic FSM calculations for bulk and monolayer NbSe 2 in Fig. 3. If we express the expansion of the DFT total energy as we can determine the static spin-susceptibility, χ, for low values of m by using the coefficient a 1 obtained by fitting the data in Fig. 3 to Eq. 2. The results are summarized in Table I.
Based on our FSM calculations, we can draw the following conclusions. First, DFT repoduces the experimentally observed bulk susceptibility reasonably well, only slightly overestimating it compared to the experimental value. This overestimation in the calculated value of χ is known to occur in itinerant metals close to a magnetic instability, and is due to a fluctuational reduction of the mean-field DFT moment [34,35]. Applying Moriya's theory [34] to NbSe 2 , we can estimate the average magnitude of spin fluctuations as χ ∼ 0.28 µ B . To put this into context, the average magnitude of spin fluctuations in palladium (Pd), a known superparamagnet (which at some point was considered a candidate for triplet superconductivity [27]) was calculated (in LDA, as opposed to our GGA calculation) to be χ ∼ 0.15 µ B [35]. This is also consistent with our disordered local moment calculations (DLM), where the energy cost of creating a local spin fluctuation with an amplitude of ∼ 0.2 µ B is nearly twice higher in Pd then in NbSe 2 . This does not mean that NbSe 2 is closer to ferromagnetism compared to Pd. The molar spin susceptibility of NbSe 2 is a factor of two lower compared to Pd. However, it does mean spin fluctuations in NbSe 2 are soft over a large part of the Brillouin zone. Furthermore, it is important to note the susceptibility of the monolayer structure q (2π/a) is ∼50% larger than that of the bulk, which indicates that spin fluctuations are stronger in a monolayer. Indeed, this is consistent with monolayer NbSe 2 having a lower superconducting transition temperature compared to bulk NbSe 2 [9].
To verify that NbSe 2 is indeed close to a ferromagnetic instability, we also calculated the exchange coupling between fluctuating moments within the DLM formalism (the calculation details are similar to methods used in Ref. [36]) for bulk NbSe 2 . The exchange coupling is largely dominated by the nearest-neighbor coupling, which we find to be ferromagnetic. In order to transform the exchange interactions into reciprocal space, J 0 (q), we have defined where c 0 is a constant of the order of N (0). In Figure 4 we plot the renormalization factor 1/[1 − c 0 J(q)], using c 0 = 5.3 eV −1 , which was chosen so as to have the renormalization at q = 0 be approximately consistent with a renormalization factor of 3.3. The peak near q = 0 in Figure  4 indicates that the system is close to a ferromagnetic instability.

C. Magnetism with spin-orbit coupling
Having demonstrated that bulk and monolayer NbSe 2 are indeed close to a ferromagnetic instability, we now focus on the effect of SO-coupling, and specifically on the response to an external magnetic field applied parallel to the c-axis and perpendicular to the c-axis. First, it is evident the values of χ reported in Table I from our relativistic calculations are isotropic along 001 and 100 , within a few percent, for both the bulk and monolayer structures, as opposed to the susceptibility in the superconducting state.
To understand this, let us analyze how the bands that cross the Fermi level evolve as a function of the magnitude and direction of an applied magnetic field. In the absence of SO-interaction, the states at Γ, K and K ′ are degenerate. The Zeeman interaction, regardless of the direction of the field, splits the bands by approximately the same magnitude, ±H, where H is the Stonerenhanced external field (we have absorbed the Bohr magneton in the units of H). Indeed, in our calculations, the splitting at Γ, K and K ′ increases linearly and by approximately the same amount as a function of increasing magnetic moment. The magnitude of the Fermi surface splitting in reciprocal space is δk F (k) =2H/v F (k). The area between the spin-split contours will determine the total magnetization for a given magnitude of the magnetic field, and will be 2N ↑ H = N ↑↓ H, where N ↑ is the total number of states with pseudospins alongẑ and N ↑↓ is the total number of states with pseudospins alongẑ and −ẑ.
We now define the following generic Hamiltonian with SO-interaction for a given point on the Fermi surface of monolayer NbSe 2 subject to an external magnetic field H z alongẑ (parallel to the c-axis), where the spinquantization axis is alongẑ.
Based on Eq. 1, λ k = λη k Im β k , and the matrix is in the S z spin space. Inversion in the momentum space changes the upper sign to the lower sign for ±λ and ∓λ. The splitting between the Fermi contours for a given pseudospin along +ẑ will increase by approximately the same magnitude as the Fermi contours of the majority spin channel in the case of the nonrelativistic calculations. The splitting of the Fermi contours for pseudospins along −ẑ will decrease by the same magnitude. For instance, if the Fermi contour splitting around K increases by 2H z , the splitting around K ′ will decrease by the same amount. Hence, the total magnetization that is induced in terms of the population of pseudospins, to lowest order in λ, will be exactly the same as determined by our nonrelativistic calculations.
In the relativistic DFT calculations, a pseudospin state is formally a combination of both spins, so we introduce an effective g factor, which describes the difference between the pure spin susceptibility and the "pseudospin" susceptibility. For a magnetic field alongẑ, the g factor is exactly 2, so we expect χ along 001 to be very similar to χ obtained from nonrelativistic calculations -which is consistent with our calculations in Table I. The splitting of the states at Γ, K and K ′ for magnetic moments alongẑ is illustrated in Fig. 5(a) and they indeed change linearly with respect to the magnetic moment.
For an in-plane magnetic field alongx (perpendicular to the c-axis), H x , the Zeeman interaction, S x H x = (S + + S − )H x /2, couples to the off-diagonal components of the spin-orbit Hamiltonian as follows: To linear order in H x , the splitting between the Fermi contours at K and K ′ will not change. However, the −H£¤ wave functions change, and thus the effective g-factor will deviate linearly from 2. For example, applying standard perturbation theory to the pseudospin |+ states, gives: Pseudospin |− states will acquire the opposite magnetization, and their g factor will be reduced by the same amount. We now observe that the total pseudomoment around the K point will be proportional to the area between the split concentric Fermi contours, ±N K λ K , where N K , is the density of states for this contour at K and λ K is the average splitting at this contour. Around the K ′ contour, the area between the split concentric Fermi contours is ∓N K λ K , and this total pseudomoment does not depend on H x . Multiplying it by the difference in the g factors of Eq. 4, which deviate from 2 by the same amount, but in the opposite directions, we get a spin susceptibilty of χ 100 ≈ N tot = χ Pauli , where N tot is the total density of states. Thus, no anisotropy in the spin susceptibility appears in the lowest order of the SO-coupling. DFT calculations fully conform with this description: the splitting of the one-electron energies at K and K ′ is quadratic with respect to magnetic moments oriented alongx (cf. Fig. 5(b)).
Within our considerations of the Zeeman interaction, H is the total magnetic field, which includes the Stoner renormalization. Within DFT, the RPA approximation is exact, since one can write the total DFT exchangecorrelation energy, E xc , in an external magnetic field as [37,38]: where I = δ 2 E xc /dm 2 is the DFT Stoner factor, which in the DFT language combines the diagonal (Hubbard U ) and off-diagonal (Hund's J) interactions. Indeed, where χ 0 is the bare Pauli susceptibility.

D. Magnetism and superconductivity
Within this framework, it is especially easy to address the effect of superconductivity on the spin susceptibility. Indeed, the opening of the superconducting gap, ∆, only affects states that are close to the Fermi surface. Since the spin susceptibility parallel to the c-axis, χ 001 , is determined entirely by the shift of the Fermi contours as a function of an increasing magnetic field (Fig. 5(a)), the spin susceptibility is suppressed by superconductivity in exactly the same way as without SO-coupling. In contrast, the spin susceptibility perpendicular to the c-axis, χ 100 , as we just saw, is defined by the states removed from the Fermi level by ∼ λ k ≫ ∆, and as a result is not affected by ∆. Thus, the thermodynamic critical field, H C0 , which is determined by the free energies in the normal and superconducting state as behaves conventionally for magnetic fields parallel to the c-axis, but is essentially infinite for magnetic fields perpendicular to the c-axis. However, if one examines the SO-splitting for the Fermi contour around the Γ-point, we find that it is nodal along the Γ-M and the Γ-M ′ line, which makes H C0 finite, but, greatly enhanced for magnetic fields perpendicular to the c-axis. Figure 2(b), illustrates the calculated splitting due to SO-coupling along the Γ contour is finite and has nodes along all Γ− M and Γ−M ′ directions. The magnitude of this splitting is low, but, along the antinodal line Γ−K, the maximum splitting is still larger than the superconducting gap. In the Appendix (Sec. C), we will derive an analytical expression, which generalizes the considerations presented above onto a SO-coupled nodal case of the Γ-centered Fermi contour.

E. Singlet and triplet superconductivity: DFT point of view
We now use the above considerations to determine the symmetry of possible pairing interactions. The fact that NbSe 2 lacks inversion symmetry formally allows for parity mixing, but it has been argued [12] that the triplet component must be vanishingly small. As we will discuss below, this is not necessarily true. Strong spinfluctuations, which we have demonstrated to be operative in NbSe 2 , and/or a particular structure of the electronphonon coupling have the ability to generate a sizeable triplet component.  Fig 2(b)). Solid circles represent pseudospin |↑ states and dotted circles represent pseudospin |↓ states. The possible pairing interactions due to phonons or spin fluctuations between the pseudospin states are denoted with red dotted arrows. The relative signs of these interactions are summarized in Table II. Subscript o refers to the outer contour and the subscript i refers to the inner contour at a given valley.
In the spirit of band theory, we consider the oneelectron Hamiltonian to be fully diagonalized before we consider superconducting pairing. First, we only consider the K and K ′ contours. The Cooper pairs at these contours are comprised of states that reside on either the inner or outer contours at K and K ′ . We assume that the outer contour around K has spin-up, and the inner contour has spin-down states, which we denote as [|K, o, ↑ ,|K, i, ↓ ] The contours at K ′ are degenerate in energy with the contours at K (the states at K and K ′ are related by time-reversal symmetry) and is given by [|K ′ , o, ↓ ,|K ′ , i, ↑ ]. No other combinations are allowed. Schematically, these four contours can be represented as two pairs of concentric rings as depicted in Fig. 6. In this basis, the anomalous averages that appear in the problem are Notice that we introduced a minus sign in Eq. 8 for d i,k , this is to ensure that a usual spin-singlet state is given by d i,k = d o,k and that the interactions we discuss later reduce to the usual interactions when no spin-orbit coupling is included. While one has the freedom to define the relative phase between the superconducting order parameters at different momenta, it is rarely used and appreciated. However, in some cases it is important to distinguish between the relative phase to ensure the pairing interactions lead to meaningful physical observable quantities [39].
Since a singlet pair is defined as (|↑↓ − |↓↑ )/ √ 2 and the triplet pair is defined as (|↑↓ + |↓↑ )/ √ 2, the order parameter on the outer contour, ∆ o , derived from the anomalous average while the order parameter on the inner contour, ∆ i , is Within this definition, ∆ S , is the order parameter for a singlet pair and ∆ T is the order parameter for a triplet pair. This picture implies four types of pairing interactions, which corresponds to the following scattering processes of Cooper pairs: then in most (albeit not necessarily all) experiments the triplet component cancels out. For example, such a situation can arise following the considerations of Shaffer et al. [12] where they take the intraband scattering within the same valley, d o,K ⇐⇒ d o,K or d i,K ⇐⇒ d i,K to be the same (denoted as g 2 in Ref. [12]), which differs from their consideration of intraband scattering between the K and [12]).
If the pairing is due to phonons, then the matrix element for intraband pairing interactions, g p oo (or g p ii ), is defined as, where the phonon Green function is included in g, The nondiagonal electron-phonon matrix element that couples the outer and inner contours, g p oi is defined as, Since scattering by phonons does not flip spin, electron-phonon coupling is only relevant for |K, ν ⇐⇒ |K, ν ′ scattering if ν = ν ′ or for |K, ν ⇐⇒ |K ′ , ν ′ if ν = ν ′ . The indices ν and ν ′ can be either o or i to denote a state on an outer (o) or inner (i) contour.
On the other hand, ferromagnetic spin-fluctuations, even though they can only couple states withing the same valley, can work within both the inter-and intraband channels. To determine the relative sign and strength of these interactions, we define the amplitude of a spin fluctuation as S, and to a reasonable approximation their correlators can be assumed to be spin-rotationally invari- , leading to an isotropic spin susceptibility, χ. Within this definition, we can now describe the interaction of an electronic state that crosses the Fermi level with a fluctuating spin moment as σ · S = σ z S z +(σ + S − +σ − S + )/2. Hence, just as we did for the intraband electron-phonon interaction (Eq. 9), we can write down the following expression for the intraband interaction due to spin fluctuations, g s oo (or g s ii ): ≈0   TABLE II. Sign of the pairing interaction that involves either phonons, g p , or ferromagnetic spin fluctuations, g s , between two states on the outer and/or inner contour at K and K ′ . The states involved in pairing are denoted as |k, ν where k is a state either at K or K ′ , ν can take on an index o or i to denote a state on the outer (o) or inner (i). The pseudospin (↑ or ↓) associated with states on each contour is depicted schematically in Figure 6. Attractive pairing interactions are positive and repulsive interactions are negative.
Note that the sign of g s oo is negative, indicating intraband spin-fluctuation mediated interactions are repulsive, as would occur in a singlet channel. We also define an expression for spin-fluctuation mediated interactions that couple the outer and inner contours, g s oi , as where the minus sign in this matrix element appears because of our phase convention for ∆ i . Notice the prefactor g s oi , is a factor of two larger than g s oo . Hence, intraband and interband interactions around a fixed valley (K or K ′ ) have the same sign as expected for singlet pairs. However, they also have distinct prefactors, which is as if the standard singlet rotational factor of 3 has distributed itself in a ratio of 2:1 between the interband and intraband contributions within the SO-split bands. This combination of pairing interactions due to phonons and spin fluctuations is summarized in Table  II. We have provided an alternative and more complete derivation of these interactions in Appendix B.
The line of reasoning that leads to ∆ o = ∆ i , is based on the fact that the bands in question are two-dimensional and nearly parabolic, so their DOS is essentially the same (which DFT calculations confirm), while the direction of the spins is anti-aligned. If the strength of the pairing interaction is similar, this hypothesis of ∆ o = ∆ i is confirmed, and the net superconducting order parameter exhibits a singlet character (see also Appendix B). However, does the strength of the pairing interaction really have to be similar on the outer and inner contours, keeping in mind that, while the difference in the DOS is vanishingly small, the k F splitting is definitely not negligible in NbSe 2 (and, as a side note, even larger in another candidate Ising superconductor, TaS 2 [10]? First, we discuss intraband interactions. They have a pairing component due to phonons that is determined by the Eliashberg function, α 2 F (q, ω), which will have a characteristic momentum, q, and energy, ω, dependence, and a pair-breaking component due to ferromagnetic spin fluctuations, which is determined by the qdependent spin susceptibility, χ(q, ω). The q-dependent spin susceptibility is sharply peaked at q ∼ 0 (Fig. 4). The q-dependence of the electron-phonon coupling likely has a non-negligible q-dependence as well.
If the intraband interaction due to phonons is the same for states on the outer and inner contours, g p oo = g p ii (an approximation adopted within Ref. [12]), then the ratio of the superconducting gap values, |∆ o |/|∆ i |, is inversely proportional to the square root of the ratio of the density of states [40] of the outer and inner contour, N o /N i , which is, as we know, essentially 1. The phase, however, depends on the net sign of the interaction: if |g p oi | > |g s oi |, the phase is the same, and the net order parameter is essentially singlet. However, if |g p oi | < |g s oi |, which is feasible, the order parameter is net-triplet.
Let us now discuss the potential for parity mixing due to a given structure of the intraband coupling. We assume that the net intraband interaction due to electronphonon, g p (q), and spin-fluctuation mediated interactions, g s (q), is peaked at small q. To be specific, we take this net interaction to have a Lorentzian dependence on q (it may as well have a sharp minimum at q=0, or have some other comparable structure within momentum space): The net intraband coupling constants, g oo and g ii , can be obtained by averaging Eq. 13 over all k and k ′ on a circular Fermi contour with a radius k F . If we now consider how g(q) changes if k F changes by δk F , where δk F is the SO-coupling induced splitting of the Fermi contours, we find: Unlike the DOS, which for a parabolic two-dimensional band does not depend on k F , this expression for the net pairing interaction due to intraband interactions depends linearly on δk F and thus on the magnitude of the SO-coupling. If ξ ≫ 2k F , it vanishes, but if ξ ∼ 2k F , it leads to a non-negligible correction. Note that in NbSe 2 δk F /k F ∼ 1/3, and in TaS 2 δk F /k F ∼ 1/2! While momentum-resolved calculations of the electronphonon coupling in monolayer NbSe 2 are underway [41], the qualitative considerations we have presented above challenges the current notion that the superconducting order parameter in monolayer NbSe 2 is purely singlet, and demonstrates the order parameter can indeed host a measurable admixture of triplet character. Finally, while this is not the main subject of our paper, we briefly comment on the ramifications and plausible experimental probes of the singlet-triplet mixing of the order parameter of monolayer NbSe 2 . Indeed, a number of recent studies have alluded to the possibility of singlet-triplet mixing of the order parameter in monolayer NbSe 2 [11,15,18] by invoking extrinsic mechanisms such as impurities and strain. The discussion we have presented above suggests this parity mixing of the order parameter can have an intrinsic origin depending on the interplay between momentum-dependent phonon and spin-fluctuation induced couplings. Experiments that attempt to elucidate this parity mixing need to access the parity-dependent coherence factors. One possibility is quasiparticle interference, where the main challenge is to separate the intraband (o−o and i−i) scattering from the interband scattering processes. Magneto-optical spectroscopy using microwaves in the deep infrared region of the spectrum is another potential experimental probe. Finally, in the spirit of Ref. [18], one expects that impurities may affect superconductivity differently, depending on the parity. All of these probes require quantitative theories, that are beyond the scope of this present paper. Our primary goal was to demonstrate that mixed parity Ising superconductivity is possible in the transition metal dichalcogenides, and we hope this will encourage further theoretical and experimental research into its manifestation.

III. CONCLUSIONS
We have developed a formalism that adapts the model theory for Ising superconductivity into first-principles DFT calculations. We demonstrated that bulk and monolayer NbSe 2 are close to a magnetic instability, and spin-fluctuation induced interactions cannot be neglected when addressing superconductivity in NbSe 2 . Finally, we outlined two parametrically admissible situations where superconductivity in monolayer NbSe 2 may be partially triplet or even predominantly triplet without invoking an external magnetic field or exchange bias, and point to the need to reexamine the symmetry of the order parameter in monolayer NbSe 2 . This perspective on the role of magnetism in monolayer NbSe 2 will also be crucial to understand and control the superconducting properties of monolayer NbSe 2 in the presence of an external magnetic field or with heterostructures between monolayer NbSe 2 and magnetic materials.
where n q = k,s c † k+q/2,s c k−q/2,s and S i,q = q,s,s ′ c † k+q/2,s σ i,s,s ′ c k−q/2,s , and σ i is a Pauli matrix. For clarity, we have allowed the spin-interaction J i (q) to depend upon spin direction i, and will later impose isotropy J i (q) = J(q). Eq. B1 assumes interactions take the same form as when inversion symmetry is present, implying we only consider inversion symmetry breaking through single particle interactions. Noting that for sufficiently large Ising spin-orbit coupling, pairing will only occur between states of opposite spin, the contribution of the above interaction towards superconductivity can be written as where we have used ρ(k) = ρ(−k) and J i (k) = J i (−k). We now examine Cooper pairs formed from Fermions near the K and K ′ points. To this end we define operators where δk o (δk i ) denote a wavevector on the outer (inner) Fermi pocket at the K point. Here we have introduced the same sign convention for d † i as in Eq. 8 in the main text. For these operators, we find intraband, g ii and g oo , and interband, g oi and g io , interactions due to electronphonon and spin can be defined as: where Q = 2K and we have imposed spin-isotropy J i (k) = J(k). From this expression, and taking J(Q + δk) ≈ 0, the coupling constants found in Table II can be readily deduced. It is instructive to consider the limit δk i = δk o → 0, then, when J(Q) ≈ 0, g ii = g oo = −g p (0) + g s (0) and g io = g oi = −g p (Q) + 2g s (0) where the constants g p (0) and g p (Q)are defined to be positive, corresponding to attractive electron-phonon interactions, and g s (0) is negative corresponding to repulsive ferromagnetic interactions. In this case, a pure singlet state corresponds to the operator [d i (k) + d o (k)]/ √ 2 for which the interaction is v s = −g p (0) − g p (Q) − 3g s (0).
A pure triplet case corresponds to the operator [d i (k) − d o (k)]/ √ 2, for which the interaction is v t = −g p (0) + g p (Q) + g s (0). (B11) These expressions reveal how spin-fluctuations strongly suppress the spin-singlet state and enhance the spintriplet state. Notice that once the spin fluctuations become sufficiently strong, that is |g p (Q)| < 2|g s (0)|, the triplet solution will have a higher T c than the singlet solution.
Appendix C: Critical field anisotropy for a nodal Fermi surface As we discuss in Sec. II D the third Fermi pocket, around Γ, has zero SOC splitting along the Γ−M and Γ−M ′ directions. Here, we rederive the expression for the spin susceptibility for this band topology that accounts for the nodes along these directions.
Assuming that the SO-splitting varies angularly as λ cos(3ϕ), we derive the the susceptibility, dm/dH x , for an in-plane magnetic field applied alongx.
where λ is the maximal SOC splitting on this Fermi contour and ζ = H/λ. This gives the same Pauli expression as before, but with a logarithmic correction: dm dH x = χ P auli (1 + 3 4 ζ 2 log ζ) (C2) In the superconducting state, λ cos(3ϕ) is replaced by λ 2 cos 2 (3ϕ) + ∆ 2 , where ∆ is the average superconducting gap along the Γ contour and ζ defined above is replaced with ζ = √ ∆ 2 + H 2 /λ. Then in the superconducting state m = 1 2π where ζ 0 = ∆/λ. That is to say, the anisotropy of the thermodynamic critical field is not infinite, but, roughly, where N Γ(K) is the DOS around the Γ contour and N K is the total DOS around the K and K ′ pockets. While this factor is formally finite, it is a very large number of the order of 10 3 . It is evident, other factors that limit the anisotropy of the critical field are more important.