#  Proximity effects and tunneling in an Ising superconductor-magnetic
insulator heterostructure

Darshana Wickramaratne  Center for Computational Materials Science, U.S. Naval
Research Laboratory, Washington, DC 20375, USA
<EMAIL_ADDRESS>Menashe Haim  The Racah Institute of
Physics, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel  Maxim
Khodas  The Racah Institute of Physics, The Hebrew University of Jerusalem,
Jerusalem 9190401, Israel  I.I. Mazin  Department of Physics and Astronomy,
George Mason University, Fairfax, VA 22030, USA  Quantum Science and
Engineering Center, George Mason University, Fairfax, VA 22030, USA

###### Abstract

Hybrid Ising superconductor-ferromagnetic insulator heterostructures provide a
unique opportunity to explore the interplay between proximity-induced
magnetism, spin-orbit coupling and superconductivity. Here we use a
combination of first-principles calculations of NbSe2/CrBr3 heterostructures
and an analytical theory of Ising superconductivity in the presence of
defects, and propose a number of nontrivial effects are at play in such
junctions. In particular, we address spin and momentum filtering in such
junctions and show that tunneling of the states at the K-point are suppressed
and the often overlooked states near the $\Gamma$-point in monolayer NbSe2
contribute to tunneling. We then show that the interplay of the proximity-
induced exchange splitting and scattering off paramagnetic defects leads to
the following nontrivial effects; an increase in the magnitude of the
superconducting gap, broadening of the tunneling peaks and $C_{2}$ symmetry
breaking. Our results provide a unifying microscopic description of tunneling
spectroscopy studies performed on these two-dimensional Josephson junction
heterostructures which form the basis for novel spintronic and superconducting
devices.

##  I Introduction

The advent of two-dimensional (2D) superconductors such as monolayers of the
transition metal dichalcogenides (TMD), NbSe2, TaS2, and TaSe2, [1, 2, 3, 4,
5] and 2D magnetic insulators offers a route to design novel devices such as
2D Josephson junctions [6]. Monolayers of the superconducting TMDs are Ising
superconductors [2, 3, 4, 5] and some of these materials have been suggested
to be on the verge of a magnetic instability[3, 7, 8]. If the insulating
material of the junction is ferromagnetic, this raises the unique opportunity
to couple superconductivity with magnetism via the proximity effect [9, 10,
11]. Superconductor/ferromagnetic insulator junctions have been used to
elucidate the fundamental properties of the superconducting contacts and are
also pursued for potential for applications in spintronics [12] or hosting
topologically protected states [13, 14].

The prospects for such devices may soon by realized given the recent discovery
of ferromagnetism in single monolayers of the chromium trihalides [15, 16] and
preliminary reports of tunneling microscopy measurements performed on bilayer
and trilayer heterostructures comprised of NbSe2 and these magnetic insulators
[17, 18, 19, 20, 21, 22, 23]. We highlight the observations of one such
experiment [22].

In the study by Kang et al. [22], tunneling conductance measurements were
performed on a vertical heterostructure comprised of CrBr3 sandwiched between
NbSe2, which serve as superconducting contacts. An in-plane magnetic field
(lower than the superconducting critical field of NbSe2), rotates the CrBr3
spins from the magnetic easy axis along $\hat{z}$ to being in-plane. This
leads to a modest increase in the conductance gap, $\Delta$, by $\sim$2% and a
concomitant increase in the broadening of the tunneling peaks by $\sim$50%.
This is counterintuitive. One may expect that when $\Delta$ increases the
width of the tunneling peaks should decrease. The same experiments find a
puzzling hysteresis in the tunneling conductance, that has an onset $\sim$ 2K
below $T_{c}$. A separate tunneling study on NbSe2/CrBr3 [19] also identified
evidence of a two-fold rotation symmetry of the superconducting state, which
is unexpected given the six-fold symmetry of the hexagonal lattice of NbSe2.

These experimental observations contain a lot of interpretative power and form
a three-pronged puzzle that we will provide microscopic insight into in this
study using a combination of first-principles calculations and analytical
calculations based on a theory of Ising superconductivity. Obtaining
microscopic insight into each of these phenomena is predicated on
understanding the subtle interplay between Ising superconductivity, magnetism,
the atomic and electronic structure of the heterostructure, the proximity-
induced magnetic field in the superconductor, and the impact of all of this on
tunneling processes.

We show the modest change in $\Delta$ and the large increase in the broadening
of the tunneling peaks are related to the combination of Ising spin-orbit
coupling [3] and spin-conserving scattering due to paramagnetic point defects
(which we refer to generically as defects henceforth) [24, 25, 26, 27, 28].
Based on this theoretical framework we present a physically intuitive
explanation for the two-fold rotational symmetry [29, 19] and the hysteresis
observed in tunneling spectroscopy by considering the role of extended defects
[25] that order in a manner that breaks the six-fold rotational symmetry of
the hexagonal lattice and thus generates anisotropic pair-breaking scattering.

##  II First-principles calculations

We begin by examining the electronic structure of the NbSe2/CrBr3/NbSe2
trilayer heterostructure using first-principles calculations (see Methods and
Supplementary Material). The atomic structure of the heterostructure with the
lowest energy is illustrated in Fig. 1(a).

Figure 1:  First-principles calculation of the electronic structure, band
alignments and interlayer coupling of the NbSe2-CrBr3-NbSe2 heterostructure.
(a) Schematic illustration of the trilayer heterostructure showing the side
view and top view. The $x,y$ and $z$ axes denote the cartesian axes. (b)
Alignment of the energy levels of NbSe2 with respect to monolayer CrBr3 at the
K-point. (c) Spin-polarized band structure of the trilayer heterostructure
around the K point. The purple bands denote contributions from CrBr3 and green
bands denote contributions from NbSe2. (d) Interlayer coupling, 2t⟂, of the
NbSe2/CrBr3/NbSe2 trilayer heterostructure as a function of momentum.
$\Lambda$ corresponds to the midpoint along the $\Gamma$-K path.

We schematically illustrate the alignment of the NbSe2 states at the Fermi
level at K with respect to the CrBr3 spin up and spin down states obtained
with respect to the vacuum level from our first-principles calculations, in
Fig. 1(b). The spin polarized band structure of the trilayer heterostructure
is shown in Fig. 1(c). The NbSe2 states reside within the spin-up gap of
CrBr3, close to the spin up conduction band states of CrBr3.

One striking change in the electronic structure of the heterostructure is the
large exchange splitting, $\lambda_{\mathrm{ex}}$, of the NbSe2 derived states
(the bias field $\mu_{B}B=\lambda_{ex}/2)$. For the heterostructure in Fig.
1(a), $\lambda_{\mathrm{ex}}$ is 121 meV between the spin up and spin down
states. The origin of this splitting can be understood as follows. In bilayer
NbSe2 [30], there are two pairs of spin degenerate bands, one pair corresponds
to states from the top monolayer while the second pair is contributed by the
bottom monolayer. Both pairs of bands are split due to interlayer coupling,
$t_{\perp}$. In the heterostructure calculations we find the Nb atoms in the
top and bottom monolayers aquire a magnetic moment, $m_{\rm Nb}$, of $\sim$
0.10 to 0.13 $\mu_{B}$ due to proximity induced coupling between the Cr and Nb
states. This manifests in a proximity induced exchange splitting,
$\lambda_{\rm ex}$, illustrated in Fig. 1(c) that breaks the spin degeneracy
of these bands.

The magnitude of $\lambda_{\mathrm{ex}}$ reflects the magnitude of orbital
overlap between the Nb and Cr $d-$electrons. Hence, it crucially depends on
the overlap between the least overlapping orbitals in the magnetic exchange
path, namely Se and Br $p$-states, which in turn is sensitive to the vertical
separation distance between the Nb and Cr atoms, $d_{\mathrm{Cr-Nb}}$. For the
different stacking configurations we explored, $d_{\mathrm{Cr-Nb}}$ varies due
to the steric interaction between the Br and Se atoms at the interface [30].

To form a commensurate trilayer heterostructure, we assume CrBr3 layer is
under biaxial tensile strain with respect to NbSe2. In reality, the strain due
to lattice mismatch will likely be relieved by a rotation of one layer with
respect to another, which leads to spatially varying stacking of NbSe2 with
respect to CrBr3. Hence, the effective overlap and $\lambda_{\mathrm{ex}}$
should be averaged over all possible mutual orientations between the two
layers, while the equilibrium distance corresponds (even ideally, in the
absence of any defects in the van-der-Waals gap) to the sterically least
favorable geometry, i.e., when Br and Se ions are aligned vertically. The
effect of this averaging [30] is that $\lambda_{\mathrm{ex}}$ at
$d_{\mathrm{Cr-Nb}}\approx 6.88$Å which, according to our calculations is the
maximal possible separation distance between NbSe2 and CrBr3, becomes
$\left\langle\lambda_{\mathrm{ex}}\right\rangle\approx 0.04$ meV, which is
equivalent to a magnetic exchange field, $B\approx$ 0.7 T.

Inserting a single monolayer of CrBr3 increases the interlayer separation
between the NbSe2 layers, which changes the interlayer coupling, $t_{\perp}$.
We determine the magnitude of 2$t_{\perp}$ (see Methods) along $\Gamma$-M and
along $\Gamma-\Lambda$ (where $\Lambda$ is the midpoint along the $\Gamma$-K
path). These results are illustrated in Fig. 1(d). Similar to the case of the
NbSe2 monolayers separated by vacuum [30], $t_{\perp}^{\Gamma}$ $\gg$
$t_{\perp}^{\mathrm{K}}$. We find $t_{\perp}^{\Gamma}$ is 31.2 meV for the
spin-up states and 18.1 meV for the spin down states at $\Gamma$ while
$t_{\perp}^{\mathrm{K}}$ is 1.8 meV for the spin-up states and 1.1 meV for the
spin down states at K. With two monolayers or more of CrBr3 (as used by Kang
et al. [22]), $t_{\perp}$ at K is suppressed significantly compared to
$t_{\perp}$ at $\Gamma$ [30].

Moving away from $\Gamma$ our calculations in Fig. 1(d) shows that
2$t_{\perp}$ along the diagonal path ($\Gamma$-K), is lower compared to
2$t_{\perp}$ along the $\Gamma$-M path. Note that the magnitude of the spin-
orbit coupling is finite along $\Gamma$-K and reaches a maximum at $\Lambda$,
while it is zero along $\Gamma$-M [3]. Hence, the orbitals that contribute the
least to $2t_{\perp}$ leads to the largest $\Delta_{\rm SOC}$. This has
crucial implications on interpreting tunneling measurements in these
heterostructures, which we discuss next.

##  III Role of magnetic exchange field and defects on the tunneling
conductance

Armed with this quantitative understanding of $\lambda_{\mathrm{ex}}$,
$t_{\perp}^{\Gamma}$, and $t_{\perp}^{\mathrm{K}}$ we proceed to perform model
Hamiltonian calculations to describe tunneling across an Ising superconductor
- ferromagnetic insulator - Ising superconductor junction [30]. We first
consider at a heuristic level, the impact of a magnetic exchange field that is
out-of-plane (parallel to $\hat{z}$ in Fig. 1(a)) and in-plane (parallel to
$\hat{x}$ in Fig. 1(a)), on the position of the conductance peak and the
broadening (full-width half maximum) of the conductance peak.

For a magnetic exchange field, $B$, that is out-of-plane ($B\parallel c$), the
spin-orbit coupling (SOC) field, $\Delta_{\mathrm{SOC}}$, and $B$ polarize the
electron spins along $\hat{z}$ regardless of the in-plane momentum. Hence, the
magnetic exchange interaction reduces the energy of the singlet Cooper pairs,
while SOC plays no role. This leads to a familiar Pauli limited
superconductivity where the critical magnetic field is of the order of the
gap, $\Delta$ [31]. Furthermore, the Cooper pairs retain their singlet
identity and are immune to the disorder scattering in accordance with the
Anderson theorem [32].

In contrast when $B\perp c$, the impact on the order parameter is weak due to
the strong $\Delta_{\mathrm{SOC}}$. However, the broadening of the conductance
peak is sensitive to $B\perp c$ and grows with the magnitude of $B$. This is
due to the fact that when $B\perp c$, the spins at $\mathbf{k}$ and
$\mathbf{k}^{\prime}$ will acquire a finite in-plane component and a finite
probability for spin-flip scattering processes. The degree of the tilt angle
is determined by the strength of $\Delta_{\mathrm{SOC}}$. Hence, the
paramagnetic defects behave as magnetic defects due to a finite in-plane $B$
[26, 28].

Figure 2:  Differential conductance $dI/dV$ as a function of the bias voltage,
$|e|V$. Results are shown for an (a) out-of-plane magnetic field that is along
$\hat{z}$ ($B\parallel c$) and an (b) in-plane magnetic field that is along
$\hat{x}$ ($B\perp c$). We use four values of the magnetic exchange field, $B$
that corresponds to $B$ = 0 $T_{c}$ (blue), 0.225 $T_{c}$ (red), 0.45 $T_{c}$
(green), and 0.67 $T_{c}$ (grey) to determine the change in differential
conductance as a function of $B$. We use T=0.5$T_{c}$, $\Delta_{\rm
SOC}$=20$T_{c}$ and a scattering rate, $\eta$=$T_{c}$ for all of our
differential conductance calculations. The inset in each panel (a) and (b)
shows the suppression of the order parameter $\Delta_{o}$ as a function of
$B$. $\Delta_{o}$ in panel (a) is calculated with T=0.5Tc and in panel (b) is
computed using T=0.74Tc. (c) Change in the peak position of the differential
conductance as when $B\parallel c$ and $B\perp c$ is illustrated on the left
vertical axis and the change in the FWHM of the differential conductance as a
function of $B\perp c$ is illustrated on the right vertical axis.

These qualitative considerations are summarized in Table 1.

Table 1: Parameters that control the position and broadening of the conductance peak as a function of the magnetic exchange field that is out-of-plane, $B\parallel c$ (for moderately low temperatures) and in-plane, $B\perp c$. $B$ denotes the magnitude of the exchange field, $\Delta_{o}$ is the order parameter, $\Delta_{\mathrm{SOC}}$ is the magnitude of spin-orbit coupling, and $\eta$ is the disorder scattering rate, which is lower than $\Delta_{\mathrm{SOC}}$.  | $B\parallel c$ | $B\perp c$
---|---|---
peak shift | $(B/\Delta_{o})^{2}$ | $(B/\Delta_{\mathrm{SO}})^{2}$
peak broadening | $0$ | $\eta B^{2}/(B^{2}+\Delta_{\mathrm{SO}}^{2})$

To put these qualitative considerations on a firm theoretical footing we use a
model band dispersion, $\xi(\mathbf{k})$, to describe the $\Gamma$-valley of
monolayer NbSe2 calculated with respect to $E_{F}$. We use this to calculate
the order parameter for $B\parallel c$ and $B\perp c$ and combine this with
the information from our first-principles calculations to determine the spin-
dependent tunneling conductance (see Methods and Supplementary Material).

The calculated $dI/dV$ for an out-of-plane $B$ is shown in Fig. 2(a). Note
that for each value of $B\parallel c$, we only find one $dI/dV$ peak at
$|e|V=2\Delta_{o}$ which is not split by Zeeman coupling. This is due to the
fact that the top and bottom NbSe2 layers undergo the same amount of proximity
induced exchange splitting, $\lambda_{\mathrm{ex}}$. This is because the
magnitude of $m_{\rm Nb}$ is the same in the top and bottom NbSe2 layers
(Sec.II). Hence, the superconducting density of states of NbSe2 is split by
the same amount and spin is conserved during tunneling which leads to the
single peak [33].

From Figure 2(a) it is evident that as the magnitude of $B$ increases, the
position of the $dI/dV$ peak decreases. This is due to the suppression of the
order parameter, $\Delta_{o}$, which is proportional to $B^{2}$. This is in
the contrast to the Zeeman split peaks in the density of states which shifts
linearly with the magnitude of the exchange field. We also find that the full-
width half maximum (FWHM) of the conductance peaks remains unchanged and is
insensitive to the amount of disorder that we consider.

In Fig. 2(b) we illustrate the calculated $dI/dV$ when $B\perp c$. We find a
number of striking changes compared to Figure 2(a). The peak position of the
$dI/dV$ decreases and is weakly dependent on the magnitude of $B$. Secondly,
the FWHM of the $dI/dV$ increases as the magnitude of the in-plane $B$
increases. This is consistent with the spin-flip scattering rate increasing
quadratically as $\eta B^{2}/2\Delta_{\mathrm{SOC}}^{2}$, where $\eta$ is the
scattering rate due to paramagnetic defects [26, 28].

In Figure 2(c) we summarize our calculations on how the peak position and FWHM
in a tunneling measurement are impacted by the magnitude and direction of $B$.
These results, which are consistent with our qualitative considerations in
Table 1, provides a physically intuitive explanation for the modest increase
in $\Delta$ and the broadening of the differential conductance that has been
observed in tunneling measurements [22].

In addition to paramagnetic point defects, as-grown NbSe2 is known to exhibit
grain boundaries [25], which breaks the $C_{6}$ rotational symmetry in a
manner determined by the orientation of the grain boundary with respect to the
basal plane. Some tunneling studies have found the a change in symmetry from
$C_{6}$ to $C_{2}$ with respect to the direction of the applied in-plane
magnetic field [29, 19]. The question we would like to ask is: can such
symmetry-breaking extended defects manifest in a tunneling conductance that
has $C_{2}$ symmetry with respect to the direction of the external magnetic
field without impacting the symmetry of the superconducting order parameter?

Within our theoretical framework the defect-induced broadening of the
tunneling conductance peaks will depend on the angle between the direction of
the applied magnetic field and the orientation of the extended defects within
the basal plane of NbSe2. This will broaden the superconducting density of
states near the conductance peak, break the $C_{6}$ rotational symmetry, but
not well above or below the conductance peak, in agreement with the existing
experimental observations [29, 19]. Within this interpretation, the symmetry
of the superconducting state is not affected. In contrast, the prevailing
interpretations invoke the emergence of a symmetry-broken superconducting
phase that is distinct from the $s$-wave Ising superconducting phase [29].
However, this would suggest that this lower-symmetry state should become more
prevalent as the temperature is lowered instead of being maximum near $T_{c}$,
which is inconsistent with what is seen in experiments.

The final aspect of this three-pronged puzzle is the hysteresis in the
tunneling conductance and why it occurs at $T\lesssim(T_{c}-2K)$ [22]. Recall
that NbSe2 is close to a magnetic instability [3, 7]. Hence, it is likely that
defects in NbSe2 may have a finite magnetic moment with a spin axis that is
along a given orientation. To determine the ramifications of this, we
calculated the scattering rate off of anisotropic magnetic impurities, within
the Born approximation. If the magnetic moment of the defect lies along the
basal plane of NbSe2, we expect it to lead to some anisotropy in the in-plane
critical field [27]. In Fig. 3(a) we illustrate the critical field as a
function of the orientation of $B\perp c$.

Figure 3:  Magnetic defects and their impact on scattering in NbSe2. (a)The
in-plane critical field, $B_{c}$ as a function of the field orientation,
specified by the angle $\phi_{B}$ formed by the magnetic field with respect to
$x$-direction. We consider the magnetic easy axis of the defect spin along
$\hat{x}$ (green), $\hat{y}$ (red) and $\hat{z}$ (grey) with spin-flip
scattering rates $\eta_{1}$,$\eta_{2}$ and $\eta_{3}$ equal to 0.25 Tc. We set
T = 0.2Tc and $\Delta_{\mathrm{SOC}}$ = 20 Tc. (b)Spin density of a single
selenium vacancy within a 10$\times$10$\times$1 supercell of monolayer NbSe2.
The different colors correspond to different signs of the magnetization. The
net magnetization is $\sim$ 0.6 $\mu_{B}$. The position of the missing
selenium atom is denoted with the black dotted circle.

From Fig. 3(a) it is clear that the suppression of the critical magnetic field
is maximum when the easy axis of the spin of the defect is orthogonal to the
applied magnetic field. When the defect spin is along $\hat{z}$, the easy axis
is perpendicular to the applied in-plane magnetic field for all orientations
of the field. Although such defects are a source of efficient pair breaking,
their effect is independent of the orientation of the applied magnetic field.
However, if the magnetic easy axis of the defect is along the basal plane of
NbSe2, this gives rise to a two-fold oscillation of the critical field. Within
this configuration, the spin flip scattering is a pair breaking process that
adds to the pair breaking effect with the applied magnetic field. Moreover,
for point defects with a perpendicular in-plane easy axis, the critical field
oscillates out-of-phase.

One candidate defect that can give rise to this phenomenon are selenium
vacancies ($V_{\rm Se}$), which have been identified to be magnetic and
present in as-grown NbSe2 [24]. The spin density of $V_{\rm Se}$ from our
first-principles calculations is illustrated in Fig. 3(b). We find a sizeable
magnetization ($\approx 0.6\ \mu_{B}$, within our 300 atom supercell).
Furthermore, the induced magnetization has a finite length scale and is
commensurate with the size of our large supercell. This places a lower bound
on the spatial extent of the magnetization, $R_{d}$, which in our calculations
is $\sim$ 15 Å.

Now consider that as the temperature is lowered below $T_{c}$, the
superconducting coherence length, $\xi$, decreases and at some point may
become lower than $R_{d}$. When $\xi<R_{d}$, scattering would occur within the
unitary limit which would result in superconductivity being suppressed at a
length scale $R_{d}$ from the site of the vacancy. This suppression only
occurs when the magnetic moment of the defect is oriented along $\hat{z}$.
When the pairing energy of the resulting “puddle” of finite magnetization,
which is $\sim\Delta^{2}N(0)R_{d}^{2}$, becomes larger than the magnetic
anisotropy energy (typically on the order of $\mu$eV for point defects), the
magnetic moment of the point defect would flop to be in-plane. We expect this
behavior to be hysteretic, as associated with any magnetic transition.

Hence, we find that both within the Born limit and unitary limit, there is a
reduction in the superconducting energy, which leads to pair breaking, albeit
for different reasons. This pair breaking effect can be avoided if the
magnetization of the point defect, is flipped to be in-plane. This change in
the direction of the magnetization of the point defect, which is likely to be
hysteretic as illustrated by the two-fold oscillation of the critical field
(Fig. 3(a)), provides a natural explanation for the hysteresis that has been
observed in the tunneling conductance at a critical temperature below $T_{c}$
which we surmize is when the coherence length of the Cooper pairs becomes
commensurate with or is lower than $R_{d}$ [22].

##  IV Conclusions

We have presented a detailed analysis of tunneling in NbSe2/CrBr3/NbSe2
heterostructures. We find CrBr3 leads to a proximity-induced magnetic moment
in NbSe2 that manifests in exchange splitting of the NbSe2 states. We show the
states at the $\Gamma$-valley contribute the most to the tunneling
conductance, as opposed to the states centered around the K-point that are
often considered in such analyses. The states at $\Gamma$ that tunnel across
the heterostructure scatter off of paramagnetic point defects, which leads to
a pronounced broadening of the tunneling peaks and a modest enhancement of the
superconducting gap when the magnetic exchange field is in-plane.

This nontrivial result offers a transparent means to interpret tunneling
spectra of Ising superconductors that are interfaced with magnetic insulators.
Our theoretical framework offers a unifying and natural explanation for
several experimental observations. This includes the $C_{2}$ “nematic”
symmetry breaking near $T_{c}$, which can be understood if one considers the
interaction of extended defects, with paramagnetic point defects and tunneling
processes. This explanation does not require one to invoke exotic symmetry-
breaking superconducting states. Finally, we show the proximity of NbSe2 to a
magnetic instability leads to defects that have a finite spin polarization
that spans a finite length scale which can manifest in a finite hysteresis in
the conductance, as has been observed in tunneling measurements at a critical
temperature below $T_{c}$.

##  V Methods

###  V.1 First-principles calculations methodology

Our first-principle calculations are based on density functional theory within
the projector-augmented wave method [34] as implemented in the VASP code [35,
36] using the generalized gradient approximation defined by the Perdew-Burke-
Ernzerhof (PBE) functional [37]. We use our DFT calculations to determine the
equilibrium atomic structure of the NbSe2/CrBr3 heterostructure, the magnitude
of the proximity-induced exchange splitting and the interlayer coupling for
different stacking configurations. For NbSe2 we found it is essential that Nb
$5s^{1},4s^{2},4p^{6},4d^{4}$ electrons and Se $4s^{2},4p^{4}$ electrons are
treated as valence. For CrBr3, Cr $4s^{1},3p^{6},3d^{5}$ and Br
$4s^{2},4p^{5}$ electrons are treated as valence. For the heterostructure
calculations we use a ($9\times 9\times 1$) $\Gamma$-centered $k$-point grid
when performing structural optimization and a ($16\times 16\times 1$)
$\Gamma$-centered $k$-point grid with the tetrahedron method to obtain
converged total energies and magnetic moments on Nb. All of the calculations
of the trilayer structures, monolayer CrBr3 and bilayer NbSe2 use a vacuum
spacing of 20Å along the $c$-axis to avoid spurious interactions with
periodically repeated surfaces.

We calculated the value of the interlayer coupling, $t_{\perp}$ for our
trilayer heterostructure by examining the magnitude of the splitting of the
NbSe2 states with the same spin character on different monolayers (i.e., the
bonding anti-bonding splitting $2t_{\perp})$. To ensure that our calculations
of the interlayer coupling are converged we use a planewave energy cutoff of
600 eV, which is significantly larger than the maximum recommended plane wave
cutoff for each of the elements in our heterostructure.

The atomic positions of the heterostructure were optimized with the Grimme-D3
van der Waals correction [38] using a force convergence criteria of 5 meV/Å.
We also investigated the effect of different van-der-Waals correction schemes
such as Grimme-D2 [39], Grimme-D3 with damping [40] and the vdW-DF functional
[41] and found they lead to minor differences in $d_{\mathrm{Nb-Nb}}$ on the
order of $\sim$0.04Å.

We performed collinear DFT calculations of a selenium vacancy ($V_{\rm Se}$)
starting from a 10$\times$10$\times$1 supercell (300 atoms) of monolayer
NbSe2. We use the $\Gamma$-point to relax the atomic coordinates. We verified
that the results are converged using a (2$\times$2$\times$1) $k$-point grid.

###  V.2 Analytical model of Ising superconductivity and tunneling

To describe the dispersion of NbSe2 we use the following model Hamiltonian.

$H_{0}(\mathbf{k})=\xi(\mathbf{k})\sigma_{0}+\left[\boldsymbol{\gamma}(\mathbf{k})-\mathbf{B}\right]\cdot\boldsymbol{\sigma}$
(1)

The SOC, $\Delta_{\mathrm{SOC}}(\mathbf{k})$, and the proximity-induced
magnetic exchange field, $B$, both couple to the spin, which is parametrized
by the set of Pauli matrices
$\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$. $\sigma_{0}$ stands
for the unit matrix in spin space. The effect of SOC at the $\Gamma$-valley is
described using
$\boldsymbol{\gamma}(\mathbf{k})=\hat{z}\sqrt{2}\Delta_{\mathrm{SOC}}\cos(3\phi_{\mathbf{k}})$,
where $\phi_{\mathbf{k}}=\arctan(k_{y}/k_{x})$ [3]. While the maximum value of
$\Delta_{\mathrm{SOC}}$ at $\Gamma$ is $\sim$70 meV, which occurs at $\Lambda$
[3], our calculations in Fig. 1(d) show that $t_{\perp}$ is suppressed at this
point. To account for the fact that $t_{\perp}$ is larger along $\Gamma$-M
while $\Delta_{\rm SOC}$ is suppressed along this path, we use a lower value
of $\Delta_{\mathrm{SOC}}$ of 20$T_{c}$ for all of our calculations.

The spin dependent density of states is expressed via the disorder-averaged
Green function, $\hat{G},$ in the standard way,

$N_{s}(\omega)=-\frac{1}{\pi}\sum_{\mathbf{k}}\mathrm{Im}\hat{G}_{ss}(\mathbf{k},\omega)\,.$
(2)

The Green function, $\hat{G}(\mathbf{k},\omega),$ and the anomalous Green
function, $\hat{F}(\mathbf{k},\omega)$ define the Nambu Green function,

$\bar{G}(\mathbf{k},\omega)=\begin{bmatrix}\hat{G}(\mathbf{k},\omega)&\hat{F}(\mathbf{k},\omega)\\\
-\hat{F}^{\ast}(-\mathbf{k},\omega)&-\hat{G}^{\ast}(-\mathbf{k},\omega)\end{bmatrix}$
(3)

which is a direct product of the particle hole and spin spaces. It satisfies
the Gor’kov equation,

$\left[(\omega+i\delta)\bar{1}-\bar{H}_{\mathrm{BdG}}(\mathbf{k})-\bar{\Sigma}(\mathbf{k})\right]\bar{G}(\mathbf{k},\omega)=\bar{1}\,,$
(4)

where $\bar{1}$ is a unit matrix in the Nambu space,
$\hat{H}_{\mathrm{BdG}}(\mathbf{k})$ is the Bogoliubov-de Gennes (BdG)
Hamiltonian, $\bar{\Sigma}(\mathbf{k})$ is the self energy accounting for the
effect of the disorder, and $\delta$ is a positive infinitesimal.

For isotropic scattering this self energy is defined within the self-
consistent Born approximation using Fermi’s golden rule as:

$\bar{\Sigma}(\mathbf{k})=\eta\int\frac{d\xi_{\mathbf{k}^{\prime}}}{\pi}\\!\\!\int\frac{d\phi_{\mathbf{k}^{\prime}}}{2\pi}\bar{\sigma}_{z}\bar{G}(\mathbf{k}^{\prime},\omega)\bar{\sigma}_{z}\,,$
(5)

where $\bar{\sigma}_{z}=\mathrm{diag}(\sigma_{0},-\sigma_{0})$, and $2\eta$ is
the scattering rate due to disorder as defined within Fermi’s Golden rule. In
practice we set $\delta$ as a small but finite number for regularization. The
BdG Hamiltonian takes the form,

$\bar{H}_{\mathrm{BdG}}(\mathbf{k})=\begin{bmatrix}H_{0}(\mathbf{k})&\Delta_{o}(\mathbf{k})\\\
\Delta_{o}^{\dagger}(\mathbf{k})&-H_{0}^{\ast}(-\mathbf{k})\end{bmatrix}\ $
(6)

where the normal state Hamiltonian,

$H_{0}(\mathbf{k})=\xi(\mathbf{k})\sigma_{0}+\left[\boldsymbol{\gamma}(\mathbf{k})-\mathbf{B}\right]\cdot\boldsymbol{\sigma}$
(7)

includes the dispersion $\xi(\mathbf{k})$ calculated with respect to $E_{F}$;
the spin orbit coupling (SOC), ${\bf\gamma}(\bf k)$, and exchange field,
$\mathbf{B}$, both of which couple to the spin parametrized by the set of
Pauli matrices $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$,
whereas $\sigma_{0}$ stands for the unit matrix in spin space.

For the order parameter, $\Delta_{o}(\mathbf{k})$ within the BdG formulation
in Eq. 6, we consider only spin singlet pairing, and assume the order
parameter to be isotropic,

$\Delta_{o}(\mathbf{k})=(i\sigma_{y})_{ss^{\prime}}\Delta_{o}\,.$ (8)

The derivation of the singlet order parameter for a magnetic field that is
out-of-plane follow Ref. [31]. For the magnetic field that is in-plane, the
order parameter is computed along the lines of [28] and is elaborated on
further in the Supplementary Material [30].

We use different approaches to determine the order parameter for these two
configurations of $B$ in the presence of disorder. When $B\parallel c$, SOC
does not affect the order parameter, $\Delta_{o}$. Indeed, in this case the
SOC field commutes with $B$ (since they are both along $\hat{z}$), and the
energy penalty of the singlet Cooper pairs due to the Zeeman field is the same
as without SOC. This allows us to compute $\Delta_{o}$ using the low-
temperature expansion of the self-consistency equations in terms of the Bessel
functions [31]. The zero-temperature and zero-field order parameter used in
the expansion relies on the BCS value, $\Delta_{\rm o}(T=0)=1.76T_{c}$. For
the range of temperatures that we investigate here, a small number of terms is
sufficient to achieve converged results. At low temperatures, $T\ll\Delta_{\rm
o}(T=0)$, we find the $\Delta_{o}$ has a rather weak field dependence that
progressively becomes more pronounced when $T\lesssim$ 0.5$T_{c}$. When
$B\perp c$, we use a Landau expansion of the free energy near $T_{c}$ to
determine $\Delta_{o}$ [30]. To ensure that $\Delta_{o}$ has the same value at
zero field for both magnetic field orientations we use $T$=0.74 $T_{c}$ for
the in-plane exchange field calculations.

We compute the conductance using the tunneling Hamiltonian approach, and only
retain the contribution of the quasi-particles to the tunnel current, $I$. The
magnitude of the conductance depends in part on the magnitude of the tunneling
matrix element, which in turn depends on the momenta and spins of the initial
and final states involved in tunneling. Since $t_{\perp}^{\Gamma}\gg
t_{\perp}^{\mathrm{K}}$ we limit our calculations to the $\Gamma$ valley,
where the spin-orbit splitting is small and angle-dependent.

The magnitude of $dI/dV$ due to tunneling across the CrBr3 band gap is
characterized by the tunneling matrix element
$\hat{t}_{\mathbf{k}s,\mathbf{k}^{\prime}s^{\prime}}$ which is determined by
the momenta and spins of the initial state, $\mathbf{k}s$, and final state,
$\mathbf{k}^{\prime}s^{\prime}$ respectively. Since the dispersion around
$\Gamma$ is approximately isotropic, we can neglect the momentum dependence of
the tunneling matrix element and define
$\hat{t}_{\mathbf{k}s,\mathbf{k}^{\prime}s^{\prime}}=t_{\perp}^{\Gamma}\delta_{ss^{\prime}}$.
A voltage, $V$, that is applied along $\hat{z}$ gives rise to an out-of-plane
tunneling current, $I$, which we define as :

$\begin{split}I(V)=e|t_{\perp}^{\Gamma}|^{2}\sum_{s}\int\frac{d\omega}{2\pi}&\left[N_{s}(\omega)N_{s}(\omega+eV)\right]\\\
&\times[n_{F}(\omega)-n_{F}(\omega+eV)]\ \end{split}$ (9)

where $e$ is the electron charge, $N_{s}(\omega)$ is the density of states for
quasiparticles with spin $s$, and $n_{F}(\omega)$ is the Fermi-Dirac function.
We include the effect of paramagnetic defects by defining a disorder self
energy, $\bar{\Sigma}(\mathbf{k})$, within the self-consistent Born
approximation.

We calculate $dI/dV$ as follows. At a given temperature, $T$ we calculate the
order parameter, $\Delta_{o}$ for $B$=0 Tc, $B$=0.225 Tc, $B$=0.45 Tc and
$B$=0.675 Tc, where Tc is the critical temperature at zero field [30]. For a
given $\Delta_{o}$ we then solve the Gor’kov equation, and find the density of
states. This allows us to compute $dI/dV$ using Eq. 9. We perform these
calculations for $B\parallel c$ and $B\perp c$.

## Data availability

The data that supports the findings of this study are available from the
corresponding author upon reasonable request.

## Competing interests

The authors declare no competing interests.

###### Acknowledgements.

We thank David Möckli, Kin Fai Mak and Jie Shan for helpful discussions. D.W
was supported by the Office of Naval Research through the Naval Research
Laboratory’s Basic Research Program. I.I.M. was supported by ONR through grant
N00014-20-1-2345. M.H. and M.K. acknowledge the financial support from the
Israel Science Foundation, Grant No. 2665/20. Calculations by D.W. and I.I.M.
were performed at the DoD Major Shared Resource Center at AFRL.

## References

  * Xi  _et al._ [2016] X. Xi, Z. Wang, W. Zhao, J.-H. Park, K. T. Law, H. Berger, L. Forró, J. Shan, and K. F. Mak, Nat. Phys. 12, 139 (2016).
  * Sergio  _et al._ [2018] C. Sergio, M. R. Sinko, D. P. Gopalan, N. Sivadas, K. L. Seyler, K. Watanabe, T. Taniguchi, A. W. Tsen, X. Xu, D. Xiao, and B. Hunt, Nat. Comm. 9, 1427 (2018).
  * Wickramaratne  _et al._ [2020] D. Wickramaratne, S. Khmelevskyi, D. F. Agterberg, and I. I. Mazin, Phys. Rev. X 10, 041003 (2020).
  * Shaffer  _et al._ [2020] D. Shaffer, J. Kang, F. Burnell, and R. M. Fernandes, Phys. Rev. B 101, 224503 (2020).
  * Möckli and Khodas [2020] D. Möckli and M. Khodas, Phys. Rev. B 101, 014510 (2020).
  * Dvir  _et al._ [2018] T. Dvir, F. Massee, L. Attias, M. Khodas, M. Aprili, C. H. Quay, and H. Steinberg, Nature communications 9, 1 (2018).
  * Divilov  _et al._ [2020] S. Divilov, W. Wan, P. Dreher, M. M. Ugeda, and F. Ynduráin, arXiv preprint arXiv:2005.06210 (2020).
  * [8] W. Wan, P. Dreher, R. Harsh, F. Guinea, and M. M. Ugeda, arXiv preprint arXiv:2101.04050 .
  * Tokuyasu  _et al._ [1988] T. Tokuyasu, J. A. Sauls, and D. Rainer, Phys. Rev. B 38, 8823 (1988).
  * Tedrow  _et al._ [1986] P. Tedrow, J. Tkaczyk, and A. Kumar, Phys. Rev. Lett. 56, 1746 (1986).
  * Heikkilä  _et al._ [2019] T. T. Heikkilä, M. Silaev, P. Virtanen, and F. S. Bergeret, Progress in Surface Science 94, 100540 (2019).
  * Eschrig [2015] M. Eschrig, Reports on Progress in Physics 78, 104501 (2015).
  * Sau  _et al._ [2012] J. D. Sau, S. Tewari, and S. D. Sarma, Phys. Rev. B 85, 064512 (2012).
  * Głodzik and Ojanen [2020] S. Głodzik and T. Ojanen, New Journal of Physics 22, 013022 (2020).
  * McGuire [2017] M. A. McGuire, Crystals 7, 121 (2017).
  * Kim  _et al._ [2019a] H. H. Kim, B. Yang, S. Li, S. Jiang, C. Jin, Z. Tao, G. Nichols, F. Sfigakis, S. Zhong, C. Li, S. Tian, D. Cory, G.-X. Miao, J. Shan, K. F. Mak, H. Lei, K. Sun, L. Zhao, and A. Tsen, Proc. Natl. Acad. Sci. U.S.A 116, 11131 (2019a).
  * Kezilebieke  _et al._ [2020a] S. Kezilebieke, M. N. Huda, V. Vaňo, M. Aapro, S. C. Ganguli, O. J. Silveira, S. Głodzik, A. S. Foster, T. Ojanen, and P. Liljeroth, Nature 588, 424 (2020a).
  * Kezilebieke  _et al._ [2020b] S. Kezilebieke, O. J. Silveira, M. N. Huda, V. Vaňo, M. Aapro, S. C. Ganguli, J. Lahtinen, R. Mansell, S. van Dijken, A. S. Foster, _et al._ , arXiv preprint arXiv:2009.13465 (2020b).
  * Hamill  _et al._ [2020] A. Hamill, B. Heischmidt, E. Sohn, D. Shaffer, K.-T. Tsai, X. Zhang, X. Xi, A. Suslov, H. Berger, L. Forró, F. Burnell, J. Shan, K. F. Mak, R. Fernandes, K. Wang, and V. Pribiag, arXiv preprint arXiv:2004.02999 (2020).
  * Kim  _et al._ [2019b] H. H. Kim, B. Yang, S. Tian, C. Li, G.-X. Miao, H. Lei, and A. W. Tsen, Nano Lett. 19, 5739 (2019b).
  * Idzuchi  _et al._ [2020] H. Idzuchi, F. Pientka, K.-F. Huang, K. Harada, Ö. Gül, Y. Shin, L. Nguyen, N. Jo, D. Shindo, R. Cava, _et al._ , arXiv preprint arXiv:2012.14969 (2020).
  * Kang  _et al._ [2021] K. Kang, S. Jiang, H. Berger, K. Watanabe, T. Taniguchi, L. Forró, J. Shan, and K. F. Mak, Giant anisotropic magnetoresistance in ising superconductor-magnetic insulator tunnel junctions (2021), arXiv:2101.01327 [cond-mat.supr-con] .
  * [23] L. Ai, E. Zhang, C. Huang, X. Xie, Y. Yang, Z. Jia, Y. Zhang, S. Liu, Z. Li, P. Leng, _et al._ , arXiv preprint arXiv:2101.04323 .
  * Nguyen  _et al._ [2017] L. Nguyen, H.-P. Komsa, E. Khestanova, R. J. Kashtiban, J. J. Peters, S. Lawlor, A. M. Sanchez, J. Sloan, R. V. Gorbachev, I. V. Grigorieva, _et al._ , ACS Nano 11, 2894 (2017).
  * Wang  _et al._ [2017] H. Wang, X. Huang, J. Lin, J. Cui, Y. Chen, C. Zhu, F. Liu, Q. Zeng, J. Zhou, P. Yu, _et al._ , Nat. Comm. 8, 1 (2017).
  * Sosenko  _et al._ [2017] E. Sosenko, J. Zhang, and V. Aji, Phys. Rev. B 95, 144508 (2017).
  * Möckli  _et al._ [2020] D. Möckli, M. Haim, and M. Khodas, Journal of Applied Physics 128, 053903 (2020).
  * Haim  _et al._ [2020] M. Haim, D. Möckli, and M. Khodas, Physical Review B 102, 214513 (2020).
  * woo Cho  _et al._ [2020] C. woo Cho, J. Lyu, T. Han, C. Y. Ng, Y. Gao, G. Li, M. Huang, N. Wang, J. Schmalian, and R. Lortz, arXiv preprint arXiv:2003.12467 (2020).
  * [30] See Supplemental Material at … for additional electronic structure calculations of the trilayer heterostructure and the derivation of the free energy of the order parameter for a magnetic exchange field that is applied parallel to the plane of the heterostructure.
  * Maki and Tsuneto [1964] K. Maki and T. Tsuneto, Progress of Theoretical Physics 31, 945 (1964).
  * Anderson [1959] P. W. Anderson, Journal of Physics and Chemistry of Solids 11, 26 (1959).
  * Meservey and Tedrow [1994] R. Meservey and P. Tedrow, Physics reports 238, 173 (1994).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Grimme [2006] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
  * Grimme  _et al._ [2011] S. Grimme, S. Ehrlich, and L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
  * Klimeš  _et al._ [2011] J. Klimeš, D. R. Bowler, and A. Michaelides, Phys. Rev. B 83, 195131 (2011).