Superconductivity in twisted Graphene ${\rm NbSe_2}$ heterostructures

TL;DR

This paper investigates how twisting graphene on NbSe2 can induce and tune superconductivity in graphene, demonstrating robustness against magnetic fields and dependence on twist angle.

Contribution

It develops a continuous low-energy model incorporating twist angle, spin-orbit coupling, and superconductivity, with parameters derived from ab-initio calculations, revealing tunable induced superconductivity.

Findings

Superconductivity can be induced in graphene over a wide range of twist angles.

The induced superconducting gap is highly robust against in-plane magnetic fields.

The size of the superconducting gap can be tuned by varying the twist angle.

Abstract

We study the low-energy electronic structure of heterostructures formed by one sheet of graphene placed on a monolayer of ${\rm NbSe_2}$. We build a continuous low-energy effective model that takes into account the presence of a twist angle between graphene and ${\rm NbSe_2}$, and of spin-orbit coupling and superconducting pairing in ${\rm NbSe_2}$. We obtain the parameters entering the continuous model via ab-initio calculations. We show that despite the large mismatch between the graphene's and ${\rm NbSe_2}$'s lattice constants, due to the large size of the ${\rm NbSe_2}$'s Fermi pockets, there is a large range of values of twist angles for which a superconducting pairing can be induced into the graphene layer. In addition, we show that the superconducting gap induced into the graphene is extremely robust to an external in-plane magnetic field. Our results show that the size of the induced superconducting gap, and its robustness against in-plane magnetic fields, can be significantly tuned by varying the twist angle.

Tables

Table 1. Table 1: Parameters for graphene- NbSe 2 subscript NbSe 2 {\rm NbSe_{2}} commensurate structures.

TMD	($m_1,m_2,n_1,n_2$)	$a_s\left(\AA \right)$	$a_g\left(\AA \right)$	$\%\delta a_g$	$\theta$	$\left|A\right|\left(\AA \right)$	$d\left(\AA \right)$	${\mu }_G\left(eV\right)$
${\mathrm{NbSe}}_2$	$(-2,1,-4,-3)$	3.48 Ding et al. (2011)	2.55	$3.7\%$	$-{65.2}^0$	9.2	3.57	-0.40
${\mathrm{NbSe}}_2$	$(-1,2,1,4)$	3.48 Ding et al. (2011)	2.55	$3.7\%$	${33.0}^0$	9.2	3.57	-0.40

Table 2. Table 2: Values of the twist angle θ 𝜃 \theta for which the graphene’s FS overlap with NbSe 2 subscript NbSe 2 {\rm NbSe_{2}} ’s FS pocket around the K 𝐾 K point or Γ Γ \Gamma point. For θ m ​ ( K ) − δ ​ θ ​ ( K ) ≤ θ ≤ θ m ​ ( K ) + δ ​ θ ​ ( K ) subscript 𝜃 m K 𝛿 𝜃 K 𝜃 subscript 𝜃 m K 𝛿 𝜃 K \mathrm{\theta_{m}(K)-\delta\theta(K)}\leq\theta\leq\mathrm{\theta_{m}(K)+\delta\theta(K)} , θ m ​ ( Γ ) − δ ​ θ ​ ( Γ ) ≤ θ ≤ θ m ​ ( Γ ) + δ ​ θ ​ ( Γ ) subscript 𝜃 m Γ 𝛿 𝜃 Γ 𝜃 subscript 𝜃 m Γ 𝛿 𝜃 Γ \mathrm{\theta_{m}(\Gamma)-\delta\theta(\Gamma)}\leq\theta\leq\mathrm{\theta_{m}(\Gamma)+\delta\theta(\Gamma)} , graphene’s FS overlaps NbSe 2 subscript NbSe 2 {\rm NbSe_{2}} ’s K 𝐾 K pocket, Γ Γ \Gamma pocket, respectively. n 𝑛 n is an integer between 0 and 5.

TMD (1L)	${\theta }_m\left(K\right)$	$\delta \theta \left(K\right)$	${\theta }_m\left(\Gamma \right)$	$\delta \theta \left(\Gamma \right)$
${\mathrm{NbSe}}_2$	$0^0+n\ast {60}^0$	${7.2}^0$	${21.9}^0+n\ast {60}^0$	${3.9}^0$
			${37.5}^0+n\ast {60}^0$	${3.9}^0$

Equations

$$m_1{\mathbf{a}}_{1s}+m_2{\mathbf{a}}_{2s}=n_1{\mathbf{a}}_{1g}+n_2{\mathbf{a}}_{2g}$$

$$(m_1^2+m_2^2-m_1{m}_2)=\frac{a_g^2}{a_s^2}(n_1^2+n_2^2-n_1{n}_2).$$

$$H_{{\mathbf{K}}_g}$$

$$H_{{\mathbf{K}}_g^{\prime }}$$

$$H_{{\Gamma }_s}={ϵ}_{0\Gamma }(\mathbf{k}){\sigma }_0+{\lambda }_{\Gamma }(\mathbf{k}){\sigma }_z$$

$${ϵ}_{0\Gamma }(\mathbf{k})$$

$${\lambda }_{\Gamma }(\mathbf{k})$$

$${\eta }_{0\Gamma }$$

$${\eta }_{2\Gamma }$$

$$l_{3\Gamma }$$

$$H_{{\mathbf{K}}_s}$$

$$H_{{\mathbf{K}}_s^{\prime }}$$

$${ϵ}_0(\mathbf{k})$$

$${ϵ}_3(\mathbf{k})$$

$$\lambda (\mathbf{k})$$

$${\eta }_0$$

$${\eta }_2$$

$${\eta }_3$$

$$l_0$$

$$l_2$$

$${\mathbf{p}}_s+{\mathbf{G}}_s={\mathbf{p}}_g+{\mathbf{G}}_g,$$

$${\mathbf{k}}_s={\mathbf{k}}_g+({\mathbf{K}}_g-{\mathbf{K}}_s)+({\mathbf{G}}_g-{\mathbf{G}}_s)$$

$${\mathbf{k}}_s={\mathbf{k}}_g+{\mathbf{K}}_g+({\mathbf{G}}_g-{\mathbf{G}}_s)$$

$$\begin{gathered} {\hat{H}}_t=\sum _{{\mathbf{G}}_g{\mathbf{G}}_s \\ \tau \sigma }{\hat{T}}_{\tau \sigma {\sigma }^{\prime }}({\mathbf{p}}_g+{\mathbf{G}}_g)e^{-i{\mathbf{G}}_g\cdot {\mathbf{d}}_{\tau }}{c}_{{\mathbf{p}}_g\tau \sigma }^{†}{d}_{{\mathbf{p}}_g+({\mathbf{G}}_g-{\mathbf{G}}_s){\sigma }^{\prime }}+h.c. \end{gathered}$$

$${\hat{H}}_t=\sum _{i=1}^3{c}_{{\mathbf{k}}_g\tau \sigma }^{†}{T}_{{\mathbf{K}}_g,i,\tau \sigma {\sigma }^{\prime }}^{†}{d}_{{\mathbf{k}}_g+{\mathbf{q}}_i,{\sigma }^{\prime }}+h.c.$$

$$\displaystyle T_{{\bf K}_{g},1}^{\dagger}=\left[\begin{array}[]{cccc}t&0&t&0\\ 0&t&0&t\end{array}\right]$$

$$\displaystyle T_{{\bf K}_{g},2}^{\dagger}=\left[\begin{array}[]{cccc}t&0&te^{-i{\bf G}_{g2}\cdot{\bf d}_{B}}&0\\ 0&t&0&te^{-i{\bf G}_{g2}\cdot{\bf d}_{B}}\end{array}\right]$$

$$\displaystyle T_{{\bf K}_{g},3}^{\dagger}=\left[\begin{array}[]{cccc}t&0&te^{-i{\bf G}_{g3}\cdot{\bf d}_{B}}&0\\ 0&t&0&te^{-i{\bf G}_{g3}\cdot{\bf d}_{B}}\end{array}\right].$$

$$C_{\mathbf{k}}^{†}$$

$$D_{\Gamma \mathbf{k}}^{†}$$

$$D_{K,\mathbf{k}}^{†}$$

$${\Psi }_{K_g{\Gamma }_s\mathbf{k}}^{†}$$

$${\Psi }_{K_g{K}_s,\mathbf{k}}^{†}$$

$$\displaystyle H_{K_{g}K_{s}}({\bf k})=\left[\begin{array}[]{cccc}H_{{\bf K}_{g}}({\bf k})&T_{{\bf K}_{g},1}&T_{{\bf K}_{g},2}&T_{{\bf K}_{g},3}\\ T_{{\bf K}_{g},1}^{\dagger}&H_{{\bf K}_{s}+{\bf G}_{sK1}}({\bf k}+{\bf q}_{1K})&0&0\\ T_{{\bf K}_{g},2}^{\dagger}&0&H_{{\bf K}_{s}+{\bf G}_{sK2}}({\bf k}+{\bf q}_{2K})&0\\ T_{{\bf K}_{g},3}^{\dagger}&0&0&H_{{\bf K}_{s}+{\bf G}_{sK3}}^{S}({\bf k}+{\bf q}_{3K})\end{array}\right].$$

$${\hat{H}}_{{\Gamma }_s}^{(\mathrm{SC})}=\sum _{\mathbf{k}}{\Psi }_{\mathbf{k}s}^{†}{H}_{{\Gamma }_s}^{(\mathrm{SC})}{\Psi }_{\mathbf{k}s},$$

$$\displaystyle H_{\Gamma_{s}}^{(\rm SC)}=\left[\begin{array}[]{cc}H_{\Gamma_{s}}({\bf k})&i\Delta_{\Gamma}\sigma_{2}\\ -i\Delta_{\Gamma}\sigma^{*}_{2}&-H^{T}_{\Gamma_{s}}(-{\bf k})\end{array}\right],$$

$${\hat{H}}_{sK}^{(\mathrm{SC})}=\sum _{{\mathbf{k}}_n}{\Psi }_{{\mathbf{k}}_s}^{†}{H}_{sK}^{(\mathrm{SC})}{\Psi }_{{\mathbf{k}}_s},$$

$$\displaystyle H_{sK}^{(\rm SC)}=\left[\begin{array}[]{cc}H_{s{\bf K}_{s}}({\bf k})&i\Delta_{K}\sigma_{2}\\ -i\Delta_{K}\sigma^{*}_{2}&-H^{T}_{s{\bf K}^{\prime}_{s}}(-{\bf k})\end{array}\right],$$

Full text

I Introduction

Transition metal dichalcogenides (TMDs) are extremely interesting materials due to their unique electronic properties Ramasubramaniam (2012); Latzke et al. (2015); Cazalilla et al. (2014); Wang et al. (2012, 2015); Gmitra et al. (2016); Duerloo et al. (2014); Ugeda et al. (2015); Xi et al. (2015a, b); Xu et al. (2014) and the fact that in recent years experimentalists have been able to isolate and probe TMD films only few atoms thick, down to the monolayer limit. Some TMDs monolayers, like ${\rm MoSe_{2}}$ and ${\rm MoS_{2}}$, are insulators with gaps of the order of 1.5-2 eV. Other TMDs monolayers, such as ${\rm NbSe_{2}}$, ${\rm NbS_{2}}$, ${\rm TaSe_{2}}$, ${\rm TaS_{2}}$ are metallic at room temperature and superconducting at low temperature. One feature that all TMDs have in common is a strong spin-orbit coupling (SOC). In monolayer TMDs the strongest effect of the SOC is a spin-splitting of the conduction and valence bands around the $K$, and $K^{\prime}$, points of the Brillouin zone (BZ) Böker et al. (2001); Ding et al. (2011); Mak et al. (2010). For the TMDs that are superconducting at low temperature, such a spin splitting causes the superconducting pairing to be of the Ising type Xi et al. (2015a) and therefore extremely robust to external in-plane magnetic fields Dvir et al. (2018a); Lu et al. (2015); Saito et al. (2015); de la Barrera et al. (2018). The ability of metallic TMDs to exhibit superconductivity even in the limit in which they are only one-atom thick, and the robustness of such superconducting state to external magnetic fields make them very interesting systems both from a fundamental point of view and for possible applications.

Recent advances in fabrication techniques have made possible the realization of van der Waals (vdW) heterostructures obtained by stacking crystals that are only few atoms thick Dean et al. (2010); Geim and Grigorieva (2013) In these structures the different layers are held together by van der Waals forces. As a consequence the crystals that can be used to create the structures, and their stacking configuration, are not limited to the configurations allowed by chemical bonds. This makes possible the realization of systems with unique properties such as graphene–topological-insulator heterostructures in which graphene has a tunable spin-orbit coupling depending on the stacking configuration Kou et al. (2013); Zhang et al. (2016, 2014); Zalic et al. (2017); Rodriguez-Vega et al. (2017).

In graphene the conduction and valence bands touch at the corners ($K$ and $K^{\prime}$ points) of the hexagonal BZ, and around these points the electrons behave as massless Dirac Fermions Novoselov et al. (2005); Castro Neto et al. (2009). This fact makes graphene an ideal semimetal in which the polarity of the carriers can easily be tuned via external gates. In addition, graphene has a very high electron mobility due to its very low concentration of defects and the fact that electron-phonons scattering processes do not contribute significantly to the resistivity for temperatures as high as room temperature Rossi et al. (2009); Das Sarma et al. (2011); Rodriguez-Vega et al. (2014). All these features make graphene an ideal system to probe, via tunneling setups, other materials and to realize novel vdW heterostructures with tunable properties. In particular, the fact, that the low energy states of graphene, in momentum space, are located just at the $K$ points of its BZ, in vdW structures implies that by simply varying the twist angle, graphene can be used as a momentum selective probe of the electronic structure, and properties, of the substrate. The work that we present below is an example of such momentum-selective probing capability of graphene. In monolayer ${\rm NbSe_{2}}$ the Fermi surface (FS) is formed by a pocket around the $\Gamma$ point, and pockets around the $K$, and $K^{\prime}$ points. Contrary to bulk ${\rm NbSe_{2}}$, in monolayer ${\rm NbSe_{2}}$ there is no selenium-like FS pocket around the $\Gamma$ point. As a consequence monolayer ${\rm NbSe_{2}}$ is expected to be a single-gap superconductor with the same gap at the $\Gamma$ pocket as at the $K$ pockets Khestanova et al. (2018). However, the $\Gamma$ and $K$ pockets differ in the magnitude, and $k$ dependence around the pocket, of the spin-splitting induced by the spin-orbit coupling. The splitting is much larger for the $K$ pockets and therefore the superconducting gap for these pockets is much more robust to external in-plane magnetic fields than for the $\Gamma$ pocket. As we show below a graphene-${\rm NbSe_{2}}$ heterostructure allows to probe separately ${\rm NbSe_{2}}$’s states around the $\Gamma$ point, and ${\bf K}$ point simply by tuning the relative twist angle between graphene and ${\rm NbSe_{2}}$ and therefore to study the difference between pockets of the interplay between spin-orbit coupling and superconducting pairing.

In this work we study vdW heterostructures formed by graphene and monolayer ${\rm NbSe_{2}}$. Our results show that despite the large mismatch between the lattice constants of graphene and ${\rm NbSe_{2}}$ in these structures a large superconducting pairing can be induced into the graphene layer. In addition, we show how such pairing depends, both in nature and structure, on the stacking configuration. Our results are relevant also to other graphene-TMD heterostructures such as the ones that can be obtained by replacing the ${\rm NbSe_{2}}$ monolayer by a monolayer of ${\rm NbS_{2}}$, ${\rm TaSe_{2}}$, or ${\rm TaS_{2}}$ that have also been shown to be superconducting at low temperature Leroux et al. (2012); Kačmarčík et al. (2010); de la Barrera et al. (2018); Navarro-Moratalla et al. (2016) and show how graphene can be used to probe in these systems the momentum-dependent superconducting gap and in particular its multiband structure.

II Method

In graphene the carbon atoms are arranged in a 2D hexagonal structure formed by two triangular sublattices, $A$ and $B$, with lattice constant $a_{g}=\sqrt{3}a=2.46\mbox{\AA}$, with $a=1.42\mbox{\AA}$ the carbon-carbon atomic distance. The 2D structure of ${\rm NbSe_{2}}$ is also formed by two triangular sublattices. One of the sublattices is formed by the Nb atoms, the other by two Se atoms symmetrically displaced by a distance $u=1.679\mbox{\AA}$ above and below the plane formed by the Nb atoms. The lattice constant of ${\rm NbSe_{2}}$ is $a_{s}=3.48\mbox{\AA}$. Ding et al. (2011). Figure 1 shows the Brillouin zone of graphene and ${\rm NbSe_{2}}$. In this figure and in the remainder we take $k_{x}$ to be in the direction connecting the valley ${\bf K}$ with its time-reversed partner ${\bf K}^{\prime}$. Figure 1 (a) shows the relative orientation of the graphene’s and ${\rm NbSe_{2}}$’s BZs for the case when the twist angle $\theta$ is zero and Figure 1 (b) for a case when $\theta\neq 0$.

To obtain the electronic structure of the graphene-${\rm NbSe_{2}}$ structure for a generic twist angle and in the presence of superconducting pairing in the ${\rm NbSe_{2}}$, we first need to estimate the charge transfer between the graphene layer and ${\rm NbSe_{2}}$, and the strength of the tunneling $t$ between graphene and the ${\rm NbSe_{2}}$ monolayer. To this effect we first obtain via ab-initio the electronic structure of a commensurate graphene-${\rm NbSe_{2}}$ structure. Let ${\bf a}_{1s}=a_{s}[\cos(\pi/3-\theta){\bf x}-\sin(\pi/3-\theta){\bf y}]$, ${\bf a}_{2s}=a_{s}[\cos(\pi/3+\theta){\bf x}+\sin(\pi/3+\theta){\bf y}]$, be the primitive lattice vectors for ${\rm NbSe_{2}}$, and ${\bf a}_{1g}=a_{g}[\cos(\pi/3){\bf x}-\sin(\pi/3){\bf y}]$, ${\bf a}_{2g}=a_{g}[\cos(\pi/3){\bf x}+\sin(\pi/3){\bf y}]$, the primitive vectors for graphene, with ${\bf x}$ and ${\bf y}$ the unit vectors in the $x$ and $y$ direction, respectively. In a commensurate stacking configuration the primitive vectors satisfy the equation:

[TABLE]

where $(m_{1},m_{2},n_{1},n_{2})$ are four integers constrained by the following second order Diophantine equation:

[TABLE]

Given that the lattice constant of graphene and NbSe are highly incommensurate with respect to each other, Eq. 1 (or, equivalently, Eq. 2) can only be satisfied for structures with primitive cells comprising a very large number of atoms ($\sim$1000). It is computationally extremely expensive to study structures with such large primitive cells using ab-initio methods. For this reason we allow for few percents strain of the graphene’s lattice so that Eq. 1 (or, equivalently, Eq. 2) can be satisfied for structures with primitive cells comprising 100 atoms or less. In general the relative strain of the graphene’s and NbSe2’s lattices will depend on the specific structure considered. We did not perform an energy minimization analysis and chose to strain graphene rather than NbSe2 for convenience. This is justified considering that the amount of charge transfer between the graphene layer and NbSe2, and the magnitude of the graphene-NbSe2 tunneling strength, are not expected to be affected by a small change of the graphene’s or NbSe2’s lattice constant

The ab-initio calculation were performed using the Quantum-Espresso package Giannozzi et al. (2009, 2017). We use full-relativistic ultrasoft pseudopotentials with the wavefunction kinetic energy cutoff of 50 Ry. We adopted the Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996) as the exchange and correlation functional. We set the vacuum thickness equal to $25\mathrm{\AA}$ to isolate the heterostructure and avoid the interactions between the periodic layers along the direction, ($z$), perpendicular to the layers. The interlayer distance between graphene and ${\rm NbSe_{2}}$ was obtained by full relaxation in the z-direction. The total energy was calculated by using a $18\times 18\times 1$ Monkhorst-Pack scheme grid for the $k$ points.

After having obtained the amount of charge transfer and the strength of the tunneling between the graphene layer and ${\rm NbSe_{2}}$ via ab-initio, we use a continuum model Mele (2010, 2012); Bistritzer and MacDonald (2011); Zhang et al. (2014) to obtain the low-energy spectrum of the graphene-${\rm NbSe_{2}}$ heterostructure for different values of the twist angle $\theta$. In general, the Hamiltonian $\hat{H}$ describing the graphene-${\rm NbSe_{2}}$ heterostructure can be written as: $\hat{H}=\hat{H}_{g}+\hat{H}_{s}+\hat{H}_{t}$ where $\hat{H}_{g}$ is the Hamiltonian for graphene, $\hat{H}_{s}$ is the Hamiltonian for ${\rm NbSe_{2}}$ and $\hat{H}_{t}$ is the term describing tunneling processes between graphene and ${\rm NbSe_{2}}$.

In graphene the low energy states are located at the ${\bf K}_{g}$ and ${\bf K}^{\prime}_{g}$ points of the BZ: ${\bf K}_{g}=(4\pi/(3a_{g}),0)$, ${\bf K}^{\prime}_{g}=(-4\pi/(3a_{g}),0)$ (and equivalent points connected by reciprocal lattice wave vectors). Close the ${\bf K}_{g}$ and ${\bf K}^{\prime}_{g}$ points in graphene the electrons, at low energies, are well described as massless Dirac fermions with Hamiltonians $\hat{H}_{{\bf K}_{g}}=\sum\limits_{{\bf k},\tau\tau^{\prime}\sigma\sigma^{\prime}}c^{\dagger}_{{\bf K}_{g}+{\bf k},\tau^{\prime}\sigma^{\prime}}H_{{\bf K}_{g}}c_{{\bf K}_{g}+{\bf k},\tau\sigma}$, $\hat{H}_{{\bf K}^{\prime}_{g}}=\sum\limits_{{\bf k},\tau\tau^{\prime}\sigma\sigma^{\prime}}c^{\dagger}_{{\bf K}^{\prime}_{g}+{\bf k},\tau^{\prime}\sigma^{\prime}}H_{{\bf K}^{\prime}_{g}}c_{{\bf K}^{\prime}_{g}+{\bf k},\tau\sigma}$, where

[TABLE]

$c^{\dagger}_{{\bf p},\tau\sigma}$ ($c_{{\bf p},\tau\sigma}$) is the creation (annihilation) operator for an electron, in the graphene sheet, with spin $\sigma$ and two-dimensional momentum $\hbar{\bf p}=\hbar(p_{x},p_{y})$, ${\bf k}$ is a wave vector measured from ${\bf K}$ (${\bf K}^{\prime}$), $v_{F}=10^{6}$ m/s is graphene’s Fermi velocity, $\mu_{g}$ graphene’s chemical potential, and ${\tau_{i}}$, ${\sigma_{i}}$ ($i=0,1,2,3)$) are the $2\times 2$ Pauli matrices in sublattice and spin space, respectively. As a consequence, when considering the states of graphene close to the ${\bf K}_{g}$ (${\bf K}^{\prime}_{g}$) point we have $H_{g}=H_{{\bf K}_{g}}$ ($H_{g}=H_{{\bf K}^{\prime}_{g}}$).

In ${\rm NbSe_{2}}$ the low energy states are located close to the $\Gamma$, ${\bf K}$, and ${\bf K}^{\prime}$ points of the BZ: ${\bf K}_{s}=(4\pi/(3a_{s}),0)$, ${\bf K}^{\prime}_{s}=(-4\pi/(3a_{s}),0)$ (and equivalent points connected by reciprocal lattice wave vectors). Close the $\Gamma$ point the effective low-energy Hamiltonian for ${\rm NbSe_{2}}$ takes the form $H_{\Gamma_{s}}=\sum_{{\bf k}\sigma\sigma^{\prime}}d^{\dagger}_{{\bf k},\sigma}H_{\Gamma_{s}}d_{{\bf k},\sigma^{\prime}}$, where $d^{\dagger}_{{\bf k},\sigma}$ ($d_{{\bf k},\sigma}$) is the creation (annihilation) operator for an electron in ${\rm NbSe_{2}}$ with momentum ${\bf k}$ and spin $\sigma$, and $H_{\Gamma_{s}}$ is the effective low energy Hamiltonian matrix for the conduction band of ${\rm NbSe_{2}}$. By fitting the ab-initio results we obtain:

[TABLE]

where

[TABLE]

$k_{\pm}=k_{x}\pm ik_{y}$, and $\eta_{0\Gamma}$, $\eta_{2\Gamma}$, $\l_{3\Gamma}$ are constants:

[TABLE]

Close to the corners of the BZ of ${\rm NbSe_{2}}$, the ${\bf K}_{s}$ and ${\bf K}^{\prime}_{s}$ points, for ${\rm NbSe_{2}}$ we have $H_{{\bf K}_{s}}=\sum_{{\bf k}\sigma\sigma^{\prime}}d^{\dagger}_{{\bf k},\sigma}H_{{\bf K}_{s}}d_{{\bf k},\sigma^{\prime}}$, $H_{{\bf K}^{\prime}_{s}}=\sum_{{\bf k}\sigma\sigma^{\prime}}d^{\dagger}_{{\bf k},\sigma}H_{{\bf K}^{\prime}_{s}}d_{{\bf k},\sigma^{\prime}}$, where ${\bf k}$ is now a wave vector measured from the ${\bf K}_{s}$, ${\bf K}^{\prime}_{s}$ point, respectively, and

[TABLE]

where,

[TABLE]

and $\eta_{0}$, $\eta_{2}$, $\eta_{3}$, $l_{0}$, $l_{2}$, are constants that we extracted from the ab-initio results for an isolated monolayer of ${\rm NbSe_{2}}$:

[TABLE]

Let ${\bf p}_{g}$, ${\bf p}_{s}$, be the wave vector of an electron in graphene, ${\rm NbSe_{2}}$, respectively. In the remainder we consider only momentum and spin conserving tunneling processes. Conservation of crystal momentum requires

[TABLE]

where ${\bf G}_{g}$ and ${\bf G}_{s}$ are reciprocal lattice vectors for graphene and ${\rm NbSe_{2}}$ respectively. For the purpose of developing a continuum low energy model for a graphene-${\rm NbSe_{2}}$ heterostructure it is more convenient to consider the twist angle $\theta$ as relative twist between BZ’s, as shown in Fig. 1. For $\theta=0$ the ${\bf K}$ point of graphene’s and ${\rm NbSe_{2}}$’s BZs are on the same axis. Depending on the value of $\theta$ we can have two situations: the low energy states of graphene, in momentum space, are close to ${\rm NbSe_{2}}$’s Fermi pockets around the ${\bf K}$ and ${\bf K}^{\prime}$ points, or, considering ${\rm NbSe_{2}}$’s extended BZ, to ${\rm NbSe_{2}}$’s Fermi pocket around the $\Gamma$ point. In the first case the conservation of the crystal momentum, Eq. (13), takes the form:

[TABLE]

where ${\bf k}_{s}$ ${\bf k}_{g}$ are momentum wave vectors measured from ${\bf K}_{g}$ and ${\bf K}_{s}$, respectively. By replacing ${\bf K}_{g}$, ${\bf K}_{s}$, with ${\bf K}^{\prime}_{g}$, and ${\bf K}^{\prime}_{s}$ in Eq. (14) we obtain the momentum conservation equation valid for momenta taken around the ${\bf K}^{\prime}$ points. In the second case Eq. (13) takes the form:

[TABLE]

and similarly for momenta around ${\bf K}^{\prime}_{g}$.

The conservation of the crystal momentum implies that the tunneling term takes the form:

[TABLE]

where ${\bf d}_{\tau}$ is the position of the carbon atom on sublattice $\tau$ within the primitive cell of the graphene sheet. For sublattice $A$ ${\bf d}_{\tau}=(0,0)$, for sublattice $B$ ${\bf d}_{\tau}=(a_{0},0)$, with $a_{0}$ the carbon-carbon distance.

Considering that, as shown in table 1, the separation $d=3.57\AA$ between the graphene sheet and ${\rm NbSe_{2}}$ is much larger than the interatomic distance in each material, in momentum space, the tunneling amplitude $t({\bf p})$ decays very rapidly as a function of ${\bf p}$ Bistritzer and MacDonald (2011) and so in Eq.(16) we can just keep the terms for which $({\bf p}_{g}+{\bf G}_{g})$ is smallest, i.e., restrict the sum to ${\bf G}_{g}=0$ and the two ${\bf G}_{g}$ that map ${\bf K}$ (${\bf K}^{\prime}$) to the two other equivalent points in the BZ and set $t=t({\bf K})$. The sum over ${\bf G}_{s}$ is restricted by the fact that we only need to keep terms for which the graphene’s and ${\rm NbSe_{2}}$’s states have energy separated by an amount of the order of $t$.

Let ${\bf q}={\bf k}_{s}-{\bf k}_{g}$. The above considerations imply that for the case when the ${\bf K}_{g}$ and ${\bf K}_{s}$ are close we only need to keep the terms for which $|{\bf q}|=|{\bf K}_{g}-{\bf K}_{s}|$, given that these are the terms for which $({\bf p}_{g}+{\bf G}_{g})$ that satisfies Eq. (14) is smallest. Due to the $C_{3v}$ symmetry of the hexagonal structure there are three equivalent ${\bf K}$ points, ${\bf K}_{1}$, ${\bf K}_{2}$, ${\bf K}_{3}$, (and ${\bf K}^{\prime}$ points), i.e. two reciprocal lattice wave vectors ${\bf G}$ connecting equivalent corners of the BZ. There are three vectors ${\bf q}_{iK}=({\bf K}_{g}-{\bf K}_{s})+({\bf G}_{gi}-{\bf G}_{si})$ ($i=1,2,3$) such that $|{\bf q}_{i}|=|{\bf K}_{g}-{\bf K}_{s}|$. ${\bf q}_{1K}$ is obtained by taking ${\bf G}_{g1}=0$ and ${\bf G}_{s}={\bf G}_{sK1}\equiv 0$, ${\bf q}_{2K}$ by taking ${\bf G}_{g}={\bf G}_{g2}\equiv 4\pi/(\sqrt{3}a_{g})[\cos(5\pi/6),\sin(5\pi/6)]$, ${\bf G}_{s}={\bf G}_{sK2}\equiv 4\pi/(\sqrt{3}a_{s})[\cos(5\pi/6+\theta),\sin(5\pi/6+\theta)]$, and ${\bf q}_{3K}$ by taking ${\bf G}_{g}={\bf G}_{g3}\equiv 4\pi/(\sqrt{3}a_{g})[\cos(7\pi/6),\sin(7\pi/6)]$, ${\bf G}_{s}={\bf G}_{sK3}\equiv 4\pi/(\sqrt{3}a_{s})[\cos(7\pi/6+\theta),\sin(7\pi/6+\theta)]$.

When the graphene’s low energy states are close to the $\Gamma$ pocket of ${\rm NbSe_{2}}$’s second BZ the smallest possible value of $|{\bf q}|$ is $|{\bf K}_{g}-{\bf G}_{s}|$ with ${\bf G}_{s}=4\pi/(\sqrt{3}a_{s})[\cos(-\pi/6+\theta),\sin(-\pi/6+\theta)]$. As before, considering the $C_{3v}$ symmetry, there are three vectors ${\bf q}_{i\Gamma}$ with this magnitude: ${\bf q}_{1\Gamma}$ obtained by taking ${\bf G}_{g}=0$, ${\bf G}_{s}={\bf G}_{s\Gamma 1}\equiv 4\pi/(\sqrt{3}a_{s})[\cos(-\pi/6+\theta),\sin(-\pi/6+\theta)]$, ${\bf q}_{2\Gamma}$ obtained by taking ${\bf G}_{g}={\bf G}_{g2}$, ${\bf G}_{s}={\bf G}_{s\Gamma 2}\equiv 4\pi/(\sqrt{3}a_{s})[\cos(\pi/2+\theta),\sin(\pi/2+\theta)]$, and ${\bf q}_{3\Gamma}$ obtained by taking ${\bf G}_{g}={\bf G}_{g3}$, ${\bf G}_{s}={\bf G}_{s\Gamma 3}\equiv 4\pi/(\sqrt{3}a_{s})[\cos(7\pi/6+\theta),\sin(7\pi/6+\theta)]$,

By retaining only the tunneling terms for which $t({\bf p}_{g}+{\bf G}_{g})$ is largest, when considering the graphene states close to the ${\bf K}_{g}$ point so that $H_{g}=H_{{\bf K}_{g}}$, we can rewrite $\hat{H}_{t}$ as

[TABLE]

with:

[TABLE]

[TABLE]

[TABLE]

In the remainder, supported by DFT results, we take $t$ to be the same both when the graphene’s low energy states tunnel into states around the ${\bf K}$ (${\bf K}^{\prime}$) point and the $\Gamma$ point of ${\rm NbSe_{2}}$. Let $\gamma\equiv t/\hbar v_{F}|{\bf q}_{i}|$. When $\gamma<1$ we can develop a perturbative approach in which $\gamma$ is the small parameter Lopes dos Santos et al. (2007); Bistritzer and MacDonald (2011): terms of order $\gamma^{n}$ correspond n-tuple tunneling processes. For our situation, as we show in the following section, $\gamma\ll 1$ and so we can retain just the lowest order terms in $\gamma$.

It is convenient to define the following spinors:

[TABLE]

For the case when the graphene’s FS overlaps with the ${\rm NbSe_{2}}$’s pocket close to the $K$ point, we can then express the Hamiltonian for the graphene-${\rm NbSe_{2}}$ system as $\hat{H}_{K_{g}K_{s}}=\sum_{\bf k}\Psi^{\dagger}_{{\bf k},K_{g}K_{s}}H_{K_{g}K_{s}}({\bf k})\Psi_{{\bf k},K_{g}K_{s}}$ with

[TABLE]

For the case when we consider graphene states close to the ${\bf K}^{\prime}_{g}$ point, so that $H_{g}=H_{{\bf K}^{\prime}_{g}}$, the expression of the Hamiltonian matrix $H_{K^{\prime}_{g}K^{\prime}_{s}}({\bf k})$ for the graphene-${\rm NbSe_{2}}$ system, within the approximations described above, can be obtained from Eq. (31) by doing the following substituions: ${\bf K}_{s}\to{\bf K}^{\prime}_{s}$, ${\bf G}_{gi}\to-{\bf G}_{gi}$, ${\bf G}_{si}\to-{\bf G}_{si}$, ${\bf q}_{iK}\to-{\bf q}_{iK}$ and noticing that $T_{{\bf K}^{\prime}_{g},i}=T^{*}_{{\bf K}_{g},i}$. Similarly, when the low energy states of graphene are close to the $\Gamma$ point of ${\rm NbSe_{2}}$ the Hamiltonian $H_{K_{g}\Gamma}({\bf k})$ ($H_{K^{\prime}_{g}\Gamma}({\bf k})$) is obtained from the expression (31) for $H_{K_{g}K_{s}}({\bf k})$ via the substitutions ${\bf K}_{s}+{\bf G}_{sKi}\to{\bf G}_{s\Gamma i}$ (${\bf K}^{\prime}_{s}-{\bf G}_{sKi}\to-{\bf G}_{s\Gamma i}$), and ${\bf q}_{iK}\to{\bf q}_{i\Gamma}$ (${\bf q}^{\prime}_{iK}\to-{\bf q}_{i\Gamma})$.

Including the superconducting pairing, the effective low-energy Hamiltonian for ${\rm NbSe_{2}}$ for states close to the $\Gamma$ point takes the form

[TABLE]

where $\Psi^{\dagger}_{{\bf k}s}$ is the Nambu spinor $\Psi^{\dagger}_{{\bf k}s}=(D^{\dagger}_{{\bf k}},D_{-{\bf k}})$,

[TABLE]

$H_{\Gamma_{s}}({\bf k})$ is given by Eq. (5), and $\Delta_{\Gamma}$ is the size of the superconducting gap of ${\rm NbSe_{2}}$ close to the $\Gamma$ point.

For states close to ${\bf K}_{s}$, including the superconducting pairing, the Hamiltonian for ${\rm NbSe_{2}}$ becomes

[TABLE]

where now ${\bf k}$ ($-{\bf k}$) is understood to be measured from ${\bf K}_{s}$ (${\bf K}^{\prime}_{s}$), and

[TABLE]

$H_{s{\bf K}_{s}}({\bf k})$, $H_{s{\bf K}^{\prime}_{s}}({\bf k})$ are given by Eq. (8).

For monolayer ${\rm NbSe_{2}}$ the superconducting gap is expected to have the same value, $\Delta$, on the $\Gamma$ and $K$ pocket. In the remainder we conservatively assume $\Delta=0.5$ meV Khestanova et al. (2018).

The Hamiltonian for the graphene-${\rm NbSe_{2}}$ system including the superconducting pairing in ${\rm NbSe_{2}}$. For the case when ${\bf K}_{g}$ is close to ${\bf K}_{s}$ the Hamiltonian becomes $\hat{H}_{K_{g}K_{s}}^{(\rm SC)}=\sum_{{\bf k}}\Psi^{\dagger}_{K_{g}K_{s},SC,{\bf k}}H_{K_{g}K_{s}}^{(\rm SC)}({\bf k})\Psi_{K_{g}K_{s},SC,{\bf k}}$, with $\Psi^{\dagger}_{K_{g}K_{s},SC,{\bf k}}=(\Psi^{\dagger}_{K_{g}K_{s},{\bf k}},\Psi^{T}_{K^{\prime}_{g}K^{\prime}_{s},-{\bf k}}$),

[TABLE]

and

[TABLE]

where $0_{m\times n}$ is the zero matrix with $m$ rows and $n$ columns.

Similarly, for the case when the low energy states of graphene are close to the $\Gamma$ point of the extended BZ of ${\rm NbSe_{2}}$ the Hamiltonian for the whole system becomes $\hat{H}_{K_{g}\Gamma_{s}}^{(\rm SC)}=\sum_{{\bf k}}\Psi^{\dagger}_{K_{g}\Gamma_{s},SC,{\bf k}}H_{K_{g}\Gamma_{s}}^{(\rm SC)}({\bf k})\Psi_{K_{g}\Gamma_{s},SC,{\bf k}}$, with $\Psi^{\dagger}_{K_{g}\Gamma_{s},SC,{\bf k}}=(\Psi^{\dagger}_{K_{g}\Gamma_{s},{\bf k}},\Psi^{T}_{K^{\prime}_{g}\Gamma_{s},-{\bf k}}$),

[TABLE]

III Results

The large lattice mismatch between graphene and ${\rm NbSe_{2}}$ would suggest that even in the absence of any twist angle the electronic states of the two systems would not hybridize. However, this does not take into account the large size of ${\rm NbSe_{2}}$’s Fermi pockets. As shown in Fig. 2 there is a large set of values of $\theta$ for which the Dirac point of graphene intersects the ${\rm NbSe_{2}}$’s FS either around the $K$ points, or around the $\Gamma$ point in the repeated zone scheme. For these points the electronic states of graphene and ${\rm NbSe_{2}}$ are expected to hybridize.

From the results shown in Fig. 2 we see that for small values of $\theta$, we can expect that the graphene’s low energy states close to the Dirac point will hybridize with the ${\rm NbSe_{2}}$’s states close to the $K$ point. For values of $\theta$ close to $30\degree$ we see that graphene’s states will hybridize with ${\rm NbSe_{2}}$’s states close to the $\Gamma$ point. For this reason, to estimate the charge transfer and the strength of the graphene-${\rm NbSe_{2}}$ tunneling in the two situations, we performed ab-initio calculations for a commensurate heterostructure with $\theta=-65.2\degree$, and one with $\theta=33.0\degree$. The parameters identifying these commensurate structures are given in table 1 and the corresponding primitive cells are shown in Fig. 3.

The ab-initio calculations return the band structure shown in Fig. 4, 5. In these figures the dashed blue lines show the bands of isolated graphene. The left panels show the results obtained without including spin-orbit effects and the right panels the results obtained taking into account the presence of spin orbit coupling. Panels (c) and (d) show an enlargement at low energies of the results shown in panels (a) and (b).

The results of Fig. 4, 5 clearly show that there is a significant charge transfer between graphene and monolayer ${\rm NbSe_{2}}$ resulting in hole doping of the graphene sheet corresponding to a Fermi energy of about -0.4 eV. They also show that the amount of charge transfer does not depend on the value of the twist angle $\theta$. Considering the finite extension of the graphene’s FS due to the charge-transfer shown in Fig. 4 5 between ${\rm NbSe_{2}}$ and graphene, we obtain that there is a significant range of values of $\theta$ for which the graphene’s FS intersects the ${\rm NbSe_{2}}$ FS and for which we can then expect non-negligible hybridization of the graphene’s and ${\rm NbSe_{2}}$ states. This is shown in Fig. 2 in which the red circles delimit the boundaries of the graphene’s FS as $\theta$ is varied. Table 2 shows the range of values of $\theta$ extracted from Fig. 2 for which the graphene’s FS is expected to intersect either one of the ${\rm NbSe_{2}}$’s FS pockets around the $K$ ($K^{\prime}$) point, or around the $\Gamma$ point. In this table $\theta_{m}(K)$ ($\theta_{m}(\Gamma)$) is the angle in the middle of the range $2\delta\theta(K)$ ($2\delta\theta(\Gamma)$) of angles for which the graphene’s FS intersects the ${\rm NbSe_{2}}$’s FS.

The ab-initio results allow us also to estimate the strength of the tunneling between graphene and ${\rm NbSe_{2}}$. In Figs. 4 (c), (d), 5 (c), (d) we can see the avoided crossings close to the Fermi energy between the graphene’s and ${\rm NbSe_{2}}$’s bands. The amplitude of such crossings provides an estimate of the tunneling strength $t$ between the graphene sheet and the monolayer of ${\rm NbSe_{2}}$. We find that both for the case when the graphene’s FS intersects the ${\rm NbSe_{2}}$’s pocket around the ${\bf K}$ point and when it intersects the ${\rm NbSe_{2}}$’s FS pocket around the $\Gamma$ point, $t\approx 20$ meV and so in the remainder we set $t=20$ meV.

We first consider the case when graphene’s FS intersects the FS pocket of ${\rm NbSe_{2}}$ close to the $K$ point, i.e. $-7.2\degree<\theta<7.2\degree$, and $\Delta=0$. Figure 6 shows the results for the FS of the hybridized system in the limit when no superconducting pairing is present in ${\rm NbSe_{2}}$: the left (right) column shows the FS around the ${\bf K}$ (${\bf K}^{\prime}$) of graphene. Figure 6 (a), (b) show the relative position in momentum space of graphene’s FS and ${\rm NbSe_{2}}$’s FS for the case when $\theta=0$ and $t=0$, taking into account the “folding” of the ${\rm NbSe_{2}}$’s FS pockets due to the fact that the three ${\bf K}$ (${\bf K}^{\prime}$) corners of the BZ are equivalent. The graphene FS is shown in red and the spin splitted ${\rm NbSe_{2}}$’s FS in blue and green. We use this color-convention throughout this work. A zoom closer to the graphene’s $K$ point, Figs 6 (c), (d), clearly shows the overlap of the graphene’s FS with the ${\rm NbSe_{2}}$’ FS pockets. When $t\neq 0$ the graphene’s and ${\rm NbSe_{2}}$’s states hybridize giving rise to the reconstructed FSs shown in Fig. 6 (e), (f). Figures 6 (e), (f) show that the graphene’s FS, due to the hybridization with ${\rm NbSe_{2}}$, becomes spin split.

Figure 7 shows the results for the case when $\theta=2\degree$, left column, and $\theta=6\degree$, right columns. For these values of the twist angle the low energy states of graphene are still close to the low energy states of ${\rm NbSe_{2}}$ located around ${\rm NbSe_{2}}$’s $K$ points. For $\theta=2\degree$ the graphene’s and ${\rm NbSe_{2}}$’s low energy states are still close enough (in momentum and energy) that, for $t=20$ meV, the hybridization is strong enough to significantly modify the FS of the combined system, as shown in Fig 7 (c), obtained setting $\Delta=0$. For $\theta=6\degree$ the graphene’s and ${\rm NbSe_{2}}$’s FSs are tangent at isolated points as shown in Fig. 7 (b). As a consequence, when $t\neq 0$ the states at the FS of graphene and ${\rm NbSe_{2}}$ only hybridize around these “tangent-points”, as shown in Fig. 7 (d) obtained for $t=20$ meV and $\Delta=0$.

We now consider the case when a superconducting gap is present in ${\rm NbSe_{2}}$. We find that for $\theta=0$ the FS is completely gapped but the gap is not uniform. Figure 8 (a) shows the lowest positive electron energy, $E_{c}$, as a function of ${\bf k}$. The smallest value of $E_{c}({\bf k})$ corresponds to the induced superconducting gap $\Delta_{\rm ind}$. For $\theta=0$ we find $\Delta_{\rm ind}=0.05$ meV. By calculating the smallest value of $E_{c}({\bf k})$ for each angle $\phi_{k}=\arctan(k_{y}/k_{x})$ we obtain the angular dependence of $\Delta_{\rm ind}$. This is shown in Fig. 8 (b) for the case when the twist angle is zero. We see that $\Delta_{\rm ind}$ is strongly anisotropic, with a $C_{3v}$ symmetry, a reflection of the structure of the reconstructed FS, Fig. 6 (e), 7 (c).

As the twist angle $\theta$ increases $\Delta_{\rm ind}$ decreases becoming vanishing small for $\theta\gtrsim 9\degree$. Figure 8 (c) shows $E_{c}({\bf k})$ when $\theta=9\degree$. From this figure we see that the location where $E_{c}({\bf k})$ is minumum appears to correspond to the original graphene’s FS for which $|{\bf k}|=k_{F,g}$. A closer inspection however reveals small oscillations as a function of $\phi_{k}$, as shown in Fig. 8 (d) where $E_{c}({\bf k})$ is plotted as function of $\phi_{k}$ and $|{\bf k}|$ for a small range of $|{\bf k}|$ centered at $k_{F,g}$.

We now consider the case when the graphene’s FS touches, in the extended BZ, the ${\rm NbSe_{2}}$’s FS pocket around the $\Gamma$ point. Figure 9 shows the results when $\theta=20\degree$, situation for which the overlap between the graphene’s FS and the ${\rm NbSe_{2}}$’s pocket at the $\Gamma$ point is largest. The left row show the results for the ${\bf K}$ point, the right the ones for the ${\bf K}^{\prime}$ point. Figure 9 (a), (b), show, on a fairly large scale, the configuration of the graphene’s and ${\rm NbSe_{2}}$’ FSs, in the absence of any interlayer tunneling, and the corresponding ${{\bf q}_{i}}$ vectors. Figure 9 (c), (d) show a zoom, at small momenta, of Fig. 9 (a) and (b), respectively, from which we can see that the graphene’s FS and the ${\rm NbSe_{2}}$’s spin-split FS intersect at several points. At these intersections the graphene’s and ${\rm NbSe_{2}}$’s states strongly hybridize causing the FS of the system to take the form shown in Fig. 9 (e), (f), for the case when $t=20$ meV, and $\Delta_{\Gamma}=0$.

As $\theta$ moves away from $20\degree$ the overlap of the graphene’s and ${\rm NbSe_{2}}$’s FSs is reduced. For $\theta=18\degree$ the overlap is still significant, the graphene’s and ${\rm NbSe_{2}}$’s FS still intersect, Fig. 10 (a), resulting in a significantly modified FS for the graphene-${\rm NbSe_{2}}$ system, Fig. 10 (c). For $\theta=16\degree$ the graphene’s and ${\rm NbSe_{2}}$’s FSs merely touch, Fig. 10 (b). As a consequence the FS of the hybridized system, for $t=20$ and $\Delta_{\Gamma}$ =0, is quite similar to the FS of the two isolated systems.

The superconducting gap on the ${\rm NbSe_{2}}$’s Gamma pocket induces a gap in the graphene layer when $\theta$ is around $22\degree$. Figure 11 (a)-(c) show the profile of $E_{c}({\bf k})$ for $\theta=(20\degree,22\degree,16\degree)$, respectively. As $\theta$ moves away from $22\degree$ $\Delta_{\rm ind}$ decrease. Figure 11 (d) show $E_{c}({\bf k})$ as function of $\phi_{k}$ and $|{\bf k}|$ for a small range of $|{\bf k}|$ centered at $k_{F,g}$ for the case when $\theta=16\degree$ and the original FSs of graphene and ${\rm NbSe_{2}}$ barely touch. As for the case then $\theta=9\degree$ we see that also for $\theta=16\degree$ $\Delta_{\rm ind}$ is very small and oscillates as function of $\phi_{k}$ for $|{\bf k}|\approx k_{F,g}$.

Using tunneling experiments Steinberg et al. (2015); Dvir et al. (2018b) it is possible to obtain the density of states, DOS, of van der Waals systems like graphene-${\rm NbSe_{2}}$. From the DOS it is then straightforward to extract the value of the induced superconducting gap. Figure 12 (a) shows the total DOS as a function of energy on a linear-log scale. We observe the coherence peaks corresponding to the ${\rm NbSe_{2}}$’s superconducting gap. Below such coherence peaks the DOS remains finite, because of the graphene’s states, until the energy is equal to $\Delta_{\rm ind}$. When the energy is equal to $\Delta_{\rm ind}$ the DOS rapidly goes to zero given that at that energy also the graphene’s states become gapped. By analyzing the DOS at small energies we can find how it depends on the twist angle, as shown in Fig. 12 (b), and (c). Figure 12 (b) shows the low energy DOS for several values of $\theta$ close to zero, i.e., for the case when ${\bf K}_{g}$ is close to ${\bf K}_{s}$, and Figure 12 (b) shows it for several values of $\theta$ close to $20\degree$, i.e., for the case when the ${\bf K}_{g}$ is close to $\Gamma$ point of ${\rm NbSe_{2}}$’s extended BZ.

From results like the ones showed in Figs. 12 (b), (c), we can extract the size of the induced superconducting gap and in particular its dependence on the twist angle, Fig. 13. We see that $\Delta_{\rm ind}$ has a fairly sharp peak for $\theta=23\degree$ (we used a $0.5\degree$ resolution) where it reaches the value of 0.087 meV. This is due to the fact that for $\theta\approx 23\degree$ there is a very strong overlap of the graphene’s and ${\rm NbSe_{2}}$ Fermi surfaces. $\Delta_{\rm ind}$ rapidly decrease as $\theta$ deviates from $23\degree$ and becomes an order of magnitude smaller when $\theta=16\degree$. $\Delta_{\rm ind}(\theta)$ has a lower and broader peak for $\theta=0$, for wich $\Delta_{\rm ind}$ =0.05 meV, i.e., for the situation in which the graphene’s FS has the maximum overlap with the ${\rm NbSe_{2}}$ $K$ pockets. As $\theta$ increases from zero $\Delta_{\rm ind}$ smoothly decreases and becomes negligible for $\theta\approx 9\degree$. Due to the symmetry of the system the behavior of $\Delta_{\rm ind}(\theta)$ has a “mirror” symmetry around $\theta=30\degree$ and is periodic with period equal to $60\degree$, as exemplified by Fig. 13. We notice that the range of values of $\theta$ for which $\Delta_{\rm ind}$ is not vanishingly small is larger than what we can infer by simply looking at the overlaps of the graphene’s and ${\rm NbSe_{2}}$’s FSs, Fig. 2. The reason is that for finite $t$ graphene’s and ${\rm NbSe_{2}}$’s states that are within the energy window $|t|$ can still hybridize resulting in a nonzero $\Delta_{\rm ind}$.

Figure 13 shows that in a graphene-${\rm NbSe_{2}}$ structure the superconducting gap can be strongly tuned by varying the twist angle and that, counterintuitively, the maximum induced gap is achieved for a value of $\theta$ for which the graphene’s FS overlaps with the $\Gamma$ pocket of ${\rm NbSe_{2}}$ in the second BZ.

Due to the strong spin-orbit coupling in ${\rm NbSe_{2}}$ the in plane critical field is much larger than the field corresponding to the Pauli paramagnetic limit. Due to the fact that SOC is also induced into the graphene layer via proximity effect we find that also for graphene-${\rm NbSe_{2}}$ heterostructures the in plane upper critical field is much larger than the Pauli paramagnetic limit. This is shown in Fig. 14 in which we plot the evolution of $\Delta_{\rm ind}$ in the presence of a Zeeman term $V_{z}$ both for values of $\theta$ corresponding to the case when the graphene’s FS overlaps ${\rm NbSe_{2}}$’s $K$ pockets (solid lines and circles), and for values of $\theta$ corresponding to the case when the graphene’s FS overlaps ${\rm NbSe_{2}}$’s $\Gamma$ pocket (dashed lines and squares). We see that in both cases $\Delta_{\rm ind}$ remains finite for $V_{z}$ as large as 40 times the induced gap of the system at zero magnetic field. However, it is also evident that the suppression of $\Delta_{\rm ind}$ due to the magnetic field is weaker, and almost independent of $\theta$, for the case when graphene’s FS overlaps ${\rm NbSe_{2}}$’s $K$ pockets. This is a consequence of the fact that in ${\rm NbSe_{2}}$ the bands’ spin splitting due to SOC is much stronger for the $K$ pockets than for the $\Gamma$ pocket.

From Fig. 14 we notice that for $\theta=22\degree$ the dependence of $\Delta_{\rm ind}$ on the Zeeman term deviates from the dependence that we find for the other values of $\theta$: $\Delta_{\rm ind}$ suddenly decreases when $V_{z}\approx 15\Delta_{\rm ind}(V_{z}=0)$, and it exhibits oscillations for larger values of $V_{z}$. The reason is that for this value of $\theta$ there are several points in momentum space for which the induced gap is close to the minimum value and, as shown in Figs. 15 (a)-(c), as $V_{Z}$ increases the point, ${\bf k}^{*}$, in momentum space where the induced gap is minimum moves. This is in contrast to what happens for other value of $\theta$, for which the gap is minimum always around the same points in $k$ space, Figs. 15 (d), regardless of the value of $V_{z}$. This implies, for $\theta=22\degree$, depending on the value of $V_{z}$ the minimum gap will be located at points with significantly different amount of SOC-induced spin splitting of the original FSs, and therefore different robustness against an in-plane magnetic field.

IV Conclusions

In conclusion, we have shown that, despite the large lattice mismatch between graphene’s and monolayer ${\rm NbSe_{2}}$’s lattice constants, in graphene-${\rm NbSe_{2}}$ heterostructures graphene exhibit a significant proximity-induced superconducting gap for a large range of stacking configurations. This is due to the fact that ${\rm NbSe_{2}}$ has large FS pockets that overlap with the FS of graphene for most twist angles. Using ab-initio calculations we have obtained the amount of charge transfer between graphene and ${\rm NbSe_{2}}$ and estimated the strength of the interlayer tunneling. We have then obtained a continuum model to describe the low energy electronic structure valid in the limit of small interlayer tunneling, condition that the ab-initio results show is satisfied. The continuum model takes into account both the presence of SOC and superconducting pairing in ${\rm NbSe_{2}}$ and the fact that, depending on the twist angle, graphene’s FS overlaps either with ${\rm NbSe_{2}}$’s FS around the ${\bf K}$ point or the $\Gamma$ point. Using this model, and the value of the parameters from ab-initio calculations, we find that, assuming conservatively the gap in ${\rm NbSe_{2}}$ monolayer to be equal to 0.5 mev, and the graphene-${\rm NbSe_{2}}$ tunneling to be 20 meV, the maximum induced superconducting gap in graphene is $\sim 0.09$ meV, obtained for a situation when the graphene FS has maximum overlap the ${\rm NbSe_{2}}$’s FS around the $\Gamma$ point. We have shown that the superconducting gap induced into the graphene layer is very robust to external in plane magnetic fields: the superconducting gap remains finite for values of the Zeeman term more than 40 times larger then the value of the induced gap in the absence of magnetic fields. In addition, we have shown that such robustness strongly depends on the twist angle in the sense that if $\theta$ is such that the graphene’s FS overlaps with the ${\rm NbSe_{2}}$ pockets around the $K$ points the induced gap is much more robust to an external in-plane magnetic field than if $\theta$ is such that the graphene’s FS overlaps with the ${\rm NbSe_{2}}$ pocket around the $\Gamma$ pocket. This is a consequence of the fact that the spin-splitting of the ${\rm NbSe_{2}}$ bands due to SOC is much stronger at the $K$ point than at the $\Gamma$ point.

The strong dependence on the external magnetic fields of the superconducting gap induced into the graphene layer is a reflection of the fact that graphene can be used, by simply varying the twist angle, as a momentum-selective probe of the electronic structure, and properties, of the substrate. We can therefore envision that tunneling experiments on graphene-based heterostructures could provide very useful, momentum selective, information on the gap structure of systems with more complex gap profiles.

Considering the similarities between the Fermi surface structure of monolayer ${\rm NbSe_{2}}$ and other transition metal dichalcogenides our results are relevant also to other graphene-TMDs heterostructures. This also applies to the case in which, instead of a monolayer, a few atomic layers TMD is used. Our results suggests that in general, for a large range of stacking configurations, the graphene and TMD states, despite the large lattice mismatch, are expected to hybridize and, when the TMD is superconducting, induce a significant superconducting gap into the graphene layer. It would be interesting to study how such proximity affect can affect the ground state of twisted-bilayer graphene systems Lu et al. (2016); Cao et al. (2018a, b); Yankowitz et al. (2019); Huang et al. (2018).

V Acknowledgments

We thank Eric Walters for helpful discussions. This work is supported by BSF Grant 2016320. YSG and ER acknowledge support from NSF CAREER grant DMR-1350663. HS is supported by European Research Council Starting Grant (No. 637298, TUNNEL). ER also thanks ONR and ARO for support. The numerical calculations have been performed on computing facilities at William & Mary which were provided by contributions from the NSF, the Commonwealth of Virginia Equipment Trust Fund, and ONR.