# Andreev and normal reflections in gapped bilayer graphene-superconductor
junctions

Panch Ram<EMAIL_ADDRESS>Fachbereich Physik, Universität Konstanz,
D-78457 Konstanz, Germany  Detlef Beckmann  Institute for Quantum Materials
and Technologies, Karlsruhe Institute of Technology, Karlsruhe D-76021,
Germany  Romain Danneau  Institute for Quantum Materials and Technologies,
Karlsruhe Institute of Technology, Karlsruhe D-76021, Germany  Wolfgang Belzig
<EMAIL_ADDRESS>Fachbereich Physik, Universität Konstanz,
D-78457 Konstanz, Germany

###### Abstract

We study the Andreev and normal reflection processes—retro as well as
specular—in a bilayer graphene-superconductor junction where equal and
opposite displacement fields are applied for the top and bottom layers to
induce a band gap. By employing the Dirac-Bogoliubov-de Gennes equation for
the gapped bilayer graphene-superconductor junction, we calculate the
reflections probabilities within the scattering theory approach. The subgap
conductance, calculated in the framework of Blonder-Tinkham-Klapwijk
formalism, shows the contribution from the Andreev retro-reflection (specular-
reflection) when the applied bias voltage is below (above) the Fermi energy.
Notably, both retro and specular reflections are modified in the presence of
the displacement field, and the retro-to-specular crossover gets amplified
when the displacement field is relatively small. They can be further tuned to
either specular or retro Andreev reflection by adjusting the Fermi energy.
Furthermore, our study reveals the simultaneous existence of double Andreev
reflections and double normal reflections when the displacement field becomes
comparable to the interlayer coupling strength. The existence of the normal
retro-reflection process in a bilayer graphene-superconductor junction is a
new finding which shows a distinctive feature in the conductance that can be
experimentally verified.

##  I Introduction

Andreev reflection (AR) is a scattering process that occurs at a normal-
superconductor (NS) junction and is solely responsible for converting a
dissipative normal current into a dissipationless supercurrent Andreev (1964);
Tinkham (2004). In this process, when an electron from the N side is incident
at the junction with excitation energy ($\varepsilon$) less than the
superconducting gap ($\Delta$), it is reflected back as a hole in a retro-
reflection manner and a charge $2e$ is transferred on S side as a Cooper pair
Bardeen  _et al._ (1957). The effect of retro Andreev reflection (RAR) on the
current-voltage relation has been studied in the seminal paper of the so-
called BTK theory Blonder  _et al._ (1982); Shelankov (1982). One consequence
of it is the existence of subgap conductance for bias voltage $eV<\Delta$
which can attain twice the value of the normal state conductance for a perfect
transparent junction. Experimentally, a direct effect of RAR has also been
observed as a sign change of velocity upon reflection Benistant  _et al._
(1983) and zero-bias anomaly Kastalsky  _et al._ (1991).

Graphene is an interesting system which has attracted the attention of the
condensed matter community, not only due to its distinctive transport
properties Novoselov  _et al._ (2005); Zhang  _et al._ (2005); Novoselov  _et
al._ (2006); Katsnelson  _et al._ (2006), owing to the relativistic linear
dispersion and vanishing density of states at the Dirac points, but also
because of a new phenomenon of electron-hole conversion at the graphene-
superconductor junction—specular Andreev reflection (SAR) Beenakker (2006,
2008). The nature of AR at the junction is expected to change from retro to
specular since the Fermi energy ($E_{F}$) in graphene is significantly lower
than the conventional metals. More precisely, when the Fermi energy is close
to the charge neutrality point (CNP), i.e. $E_{F}\to 0$, the incident electron
from the conduction band (above the CNP) is reflected (in a specular manner)
as a hole from the valence band (below the CNP). This interband electron-hole
conversion phenomenon is known as the SAR. The study of SAR Beenakker (2006,
2008) has prompted a surge of interest in the graphene-superconductor
junction, and thereafter various theoretical works have been put forward
Bhattacharjee and Sengupta (2006); Linder and Sudbø (2007, 2008); Zhang  _et
al._ (2008); Benjamin and Pachos (2008); Majidi and Zareyan (2012). But, the
experimental studies of Andreev processes have been limited so far Popinciuc
_et al._ (2012); Han  _et al._ (2018); Pandey  _et al._ (2019); Bhandari  _et
al._ (2020); Pandey  _et al._ (2021, 2022); Jois  _et al._ (2023) and the
observation of SAR in graphene has not been succeeded yet Sahu  _et al._
(2016) as the Fermi-energy fluctuation is $\delta E_{F}>\Delta$ which prevents
attaining the $E_{F}\to 0$ limit.

In contrast, the bilayer graphene (BLG) is regarded as a more suitable system
compared to the monolayer graphene for observing the SAR, due to its low
$\delta E_{F}$ fluctuation Efetov  _et al._ (2016); Efetov and Efetov (2016);
Ludwig (2007); Takane  _et al._ (2017). Indeed, an experiment has been
performed Efetov  _et al._ (2016) and the measured subgap differential
conductance shows a small dip around the bias $eV=E_{F}$ when $E_{F}$ is tuned
below $\Delta$. This characteristic feature has been accredited to a crossover
from retro-to-specular Andreev reflection. However, the retro-to-specular
crossover color plots show a very weak consent between the experiment and
theory that are presented in Figs. 2 and 3(a) in Ref. Efetov  _et al._ (2016),
see especially the different colorbar scales used for the experimental and
theoretical data plots. This raises speculation regarding the observation of
SAR in BLG. Moreover, a theoretical study to enhance the SAR contribution in
subgap conductance has been suggested by introducing a Zeeman field on the
normal side of the bilayer graphene NS junction Soori  _et al._ (2018).

In this paper, we propose an alternate study to enhance the retro-to-specular
crossover in BLG by applying two different displacement fields. Without loss
of generality, we consider that the generated electrostatic potentials, due to
the displacement fields, are $\lambda$ for the top layer and $-\lambda$ for
the bottom layer which creates a gap $\sim 2\lambda$ in the BLG band structure
when $\lambda\ll t_{\perp}$, where $t_{\perp}$ is being the interlayer
coupling. Consequently, the subgap differential conductance across the NS
junction vanishes and widens the retro-to-specular crossover, from a point at
bias $eV=E_{F}$ to a whole range of bias $|E_{F}-\lambda|<eV<|E_{F}+\lambda|$
with width $2\lambda$. This enhancement due to the displacement field would
possibly help to observe the crossover experimentally as the boundary becomes
wider. In addition, we also investigate the reflection processes and the
corresponding differential conductance when the displacement field becomes
large and comparable to the interlayer coupling, i.e., $\lambda\lesssim
t_{\perp}$. Interestingly, in this regime, four reflection processes exist
simultaneously: specular normal reflection (SNR), retro normal reflection
(RNR), specular Andreev reflection (SAR), and retro normal reflection (RAR).
This occurs due to the partial inversion of the lower conduction/valence band
upon the applied displacement field. The existence of RNR is a new finding in
the system of bilayer graphene-superconductor junction which exhibits a
distinctive characteristic feature in the differential conductance.

The rest of the paper is organized as follows. In Sec. II, we present the
effective low-energy Hamiltonian of the gapped bilayer graphene and the
corresponding band structure. The model and formalism are illustrated in Sec.
III, introducing the Dirac-Bogoliubov-de Gennes (DBdG) equation and the
excitation energy eigenstates for the normal and superconducting sides of the
junction. Sec. IV provides the results for the reflection probabilities and
normalized conductance. They are discussed in detail separately for the small
and large displacement fields in subsections IV.1, IV.2, and IV.3,
respectively. Sec. V concludes the paper.

##  II Gapped bilayer graphene

The bilayer graphene comprises two layers of graphene, wherein two non-
equivalent $A_{1}$ and $B_{1}$ carbon atoms in one layer and $A_{2}$ and
$B_{2}$ in the other layer. The two layers are stacked either in
$A_{1}$-$A_{2}$ or in $B_{1}$-$A_{2}$ structure Yan  _et al._ (2011); Brown
_et al._ (2012). A quantum Monte Carlo simulation suggests that the latter
structure is more stable Mostaani  _et al._ (2015), so we consider the
$B_{1}$-$A_{2}$ stacking which is commonly known as the Bernal stacking Yan
_et al._ (2011), see in Fig. 1(i). Within the tight-binding approximation and
considering only the first nearest-neighbour hopping in each layer and
$B_{1}$-$A_{2}$ interlayer hopping, the relevant low-energy Hamiltonian can be
deduced in momentum-space McCann and Fal’ko (2006); McCann (2006); Barbier
_et al._ (2009); Castro Neto  _et al._ (2009); McCann and Koshino (2013);
Kraft  _et al._ (2018a); Wu  _et al._ (2020) which reads in basis
$(\psi_{A_{1}}~{}\psi_{B_{1}}~{}\psi_{A_{2}}~{}\psi_{B_{2}})^{T}$ as

$\displaystyle H=\begin{pmatrix}\lambda&\hbar v_{F}k_{-}&0&0\\\\[1.99997pt]
\hbar v_{F}k_{+}&\lambda&-t_{\perp}&0\\\\[1.99997pt]
0&-t_{\perp}&-\lambda&\hbar v_{F}k_{-}\\\\[1.99997pt] 0&0&\hbar
v_{F}k_{+}&-\lambda\end{pmatrix}$ (1)

Here, $k_{\pm}=k_{x}\pm i\eta k_{y}\equiv-i(\partial_{x}\pm
i\eta\partial_{y})$ are the wave-vectors that are measured from the corners of
$2$D hexagonal Brillouin zone’s Dirac points $K(K^{\prime})$ for different
valley $\eta=+(-)$; $v_{F}\simeq 10^{6}$ m/s is the Fermi-velocity and $\hbar$
is the reduced Planck’s constant. $\lambda$ is added to account for the on-
site potential (equal in magnitude but opposite in sign for the two layers)
which can be tuned by different displacement fields. A schematic depiction to
induce $\lambda$ is shown in Fig. 1(ii). The non-zero $\lambda$ creates a band
gap in BLG band structure, and the gap can be tuned up to $250$ meV Ohta  _et
al._ (2006); Oostinga  _et al._ (2008); Zhang  _et al._ (2009); Taychatanapat
and Jarillo-Herrero (2010); Varlet  _et al._ (2014); Kraft  _et al._ (2018b);
Du  _et al._ (2018). We set $\hbar v_{F}=1$ for the calculation hereafter.

Figure 1: Bernal-stacked bilayer graphene structure and a schematic depiction
to induce the displacement field are in (i) and (ii), respectively. In (iii),
the low-energy band structure of bilayer graphene, $E_{\nu,\pm}$ for $\nu=\pm$
in Eq. (2), with respect to $k$. The behavior of conduction bands
($E_{+,\pm}$) and valence bands ($E_{-,\pm}$) for the displacement field (a)
$\lambda=0$, (b) $\lambda=0.06t_{\perp}$, and (c) $\lambda=0.6t_{\perp}$. We
observe a band gap opening and the Mexican-hat-like shape formation in
$E_{+,-}$ and $E_{-,-}$ as $\lambda$ is increased.

The eigenvalues of $H$ in Eq. (1) provide valley degenerated low-energy bands

$\displaystyle
E_{\nu,\pm}=\nu\sqrt{k^{2}+\frac{t^{2}_{\perp}}{2}+\lambda^{2}\pm\sqrt{\Big{(}\frac{t^{2}_{\perp}}{2}\Big{)}^{2}+k^{2}\left(4\lambda^{2}+t^{2}_{\perp}\right)}}$
(2)

where $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and the index $\nu=+(-)$ labels the
conduction (valence) bands. In Fig. 1(iii), we plot the conduction bands
$E_{+,\pm}$ (solid lines) and the valence bands $E_{-,\pm}$ (dashed lines) as
a function of $k$. The nature of BLG dispersion is quadratic, unlike the
linear dispersion for monolayer graphene; the former is due to hopping between
the layers, $t_{\perp}$. In absence of the displacement field, $\lambda=0$,
the lower conduction band $E_{+,-}$ and the upper valence band $E_{-,-}$ touch
at the zero energy at $k=0$ [see in Fig. 1(iii)(a)]; however, the bands
$E_{+,+}$ and $E_{-,+}$ show a gap of $\pm t_{\perp}$. When $\lambda$ is
present, we notice a gap opening and Mexican-hat-like shape in $E_{+,-}$ and
$E_{-,-}$. The results are shown in Figs. 1(iii)(b)-(c). Analytical evaluation
of the bands $E_{\pm,-}$ provides three extremal points $k_{I}=0$ and
$k_{II}^{\pm}=\pm\sqrt{\frac{4\lambda^{4}+2(\lambda
t_{\perp})^{2}}{4\lambda^{2}+t_{\perp}^{2}}}$ which correspond to energies
$E_{\pm,-}=\pm\lambda$ and $E_{\pm,-}=\pm\lambda
t_{\perp}/\sqrt{4\lambda^{2}+t_{\perp}^{2}}$, respectively. This leads to the
band gaps $2\lambda$ at $k_{I}$ and $2\lambda
t_{\perp}/\sqrt{4\lambda^{2}+t_{\perp}^{2}}$ at $k_{II}^{\pm}$. Notably, when
$\lambda\ll t_{\perp}$, the two gaps approximately become equal for $k_{I}$
and $k_{II}^{\pm}$ as can be seen in Fig. 1(iii)(b). In contrast, they are
different when $\lambda$ is comparable to $t_{\perp}$, i.e. $\lambda\lesssim
t_{\perp}$, see in Fig. 1(iii)(c). The minimum and maximum of $E_{+,-}$ are
respectively $E_{\rm min}=\lambda t_{\perp}/\sqrt{4\lambda^{2}+t_{\perp}^{2}}$
and $E_{\rm max}=\lambda$.

##  III Model and Formalism

In order to investigate the scattering processes and the transport properties,
we consider an NS junction on the gapped BLG sheet formed at $x=0$ in $x$-$y$
plane and assume that $x<0$ is the region N, while $x>0$ occupies the region
S. The superconductivity in the S region can be induced through the proximity
effect by covering an external $s$-wave superconducting electrode Volkov  _et
al._ (1995); Beenakker (2006). We employ the Dirac-Bogoliubov-de Gennes (DBdG)
equation de Gennes (1966); Beenakker (2006) which couples electron with the
time-reversed hole excitation wavefunctions via the superconducting (SC) pair
potential. It reads

$\displaystyle\begin{pmatrix}{\cal H}-E_{F}&\Delta({\bf r})\\\\[3.00003pt]
\Delta({\bf r})&E_{F}-{\cal T}{\cal H}{\cal
T}^{-1}\end{pmatrix}\begin{pmatrix}u_{e}\\\\[3.00003pt]
v_{h}\end{pmatrix}=\varepsilon\begin{pmatrix}u_{e}\\\\[3.00003pt]
v_{h}\end{pmatrix}$ (3)

where $\varepsilon\geq 0$ is the excitation energy, $u_{e}(v_{h})$ is the
electron (hole) spinor wavefunction, and ${\cal T}$ denotes time-reversal
operator. We consider the SC gap to be uniform with $\Delta({\bf
r})=\Delta\Theta(x)$ where ${\bf r}=(x,y)$ and $\Theta(x)$ is the Heaviside
step function. Also, ${\cal H}$ = diag(${\cal H}_{+},{\cal H}_{-}$) with
${\cal H}_{\pm}=H_{\pm}+U({\bf r})$ where $H_{\pm}$ are given in Eq. (1) for
the valley $\eta=\pm$ and $U({\bf r})=-U_{0}\Theta(x)$ is a doping potential
in S region. For $U_{0}\gg E_{F}$, the Fermi wavelength on the N side is much
larger compared to the S side which enables $\Delta({\bf r})$ to attain its
bulk value $\Delta$. In the absence of magnetic field, the Eq. (3) preserves
the time-reversal symmetry such that ${\cal T}{\cal H}{\cal T}^{-1}={\cal H}$.
However, this operation interchanges the valley $K$ to $K^{\prime}$ and vice-
versa. So, we can decouple the Eq. (3) into two sets of eight equations,

$\displaystyle\begin{pmatrix}{\cal H}_{\eta}-E_{F}&\Delta({\bf
r})\\\\[3.00003pt] \Delta({\bf r})&E_{F}-{\cal
H}_{\eta}\end{pmatrix}\begin{pmatrix}u_{e}\\\\[3.00003pt]
v_{h}\end{pmatrix}=\varepsilon\begin{pmatrix}u_{e}\\\\[3.00003pt]
v_{h}\end{pmatrix}.$ (4)

Notice that the subscript $\eta$ in the above equation serves another purpose;
the DBdG Eq. (4) for $\eta=+(-)$ couples electron excitation from the valley
$K(K^{\prime})$ to hole excitation at the valley $K^{\prime}(K)$.

To find the excitation energy spectrum for Eq. (4) on both sides of the
junction, we follow the similar calculation procedures as given in Ref.
Beenakker (2006) and consider a plane-wave solution
$(u_{e}~{}v_{h})^{T}e^{ik_{x}x+ik_{y}y}$. For the N side, the excitation
energies for electron ($e$) and hole ($h$) are

$\displaystyle\varepsilon_{\nu,\pm}^{e}=E_{\nu,\pm}-E_{F}\quad\mathrm{and}\quad\varepsilon_{\nu,\pm}^{h}=E_{F}-E_{\nu,\pm}$
(5)

Similarly, for the S side, the excitation energies for electron-like and hole-
like quasiparticles are

$\displaystyle\varepsilon_{\pm,\pm}^{(S,1)}$
$\displaystyle=\pm\sqrt{[(U_{0}+E_{F})\pm\gamma_{1}]^{2}+\Delta^{2}}$ (6a)
$\displaystyle\varepsilon_{\pm,\pm}^{(S,2)}$
$\displaystyle=\pm\sqrt{[(U_{0}+E_{F})\pm\gamma_{2}]^{2}+\Delta^{2}}$ (6b)

where
$\gamma_{i}=\sqrt{k^{2}+\frac{t_{\perp}^{2}}{2}+\lambda^{2}+(-1)^{i}\sqrt{\big{(}\frac{t_{\perp}^{2}}{2}\big{)}^{2}+k^{2}\left(4\lambda^{2}+t_{\perp}^{2}\right)}}$
for $i=1,2$.

Since the system is translationally invariant along the $y$-direction, the
transverse momentum $k_{y}$ is conserved during the scattering process.
Therefore, for a given $\varepsilon$ and $k_{y}$, we solve the Eq. (4) on the
N side to obtain the state vector and longitudinal momentum $k_{x}$. The
$\eta$-dependent states for electron and hole are

$\displaystyle u_{e}^{\eta}(\varepsilon,k_{x})$
$\displaystyle=\frac{1}{N_{e}}\begin{pmatrix}-t_{\perp}k_{-}(\lambda+\varepsilon_{+})\\\\[1.99997pt]
t_{\perp}(\lambda^{2}-\varepsilon_{+}^{2})\\\\[1.99997pt]
[(\lambda-\varepsilon_{+})^{2}-k^{2}](\lambda+\varepsilon_{+})\\\\[1.99997pt]
k_{+}[(\lambda-\varepsilon_{+})^{2}-k^{2}]\end{pmatrix}$ (7) $\displaystyle
v_{h}^{\eta}(\varepsilon,k_{x})$
$\displaystyle=\frac{1}{N_{h}}\begin{pmatrix}-t_{\perp}k_{-}(\lambda+\varepsilon_{-})\\\\[1.99997pt]
t_{\perp}(\lambda^{2}-\varepsilon_{-}^{2})\\\\[1.99997pt]
[(\lambda-\varepsilon_{-})^{2}-k^{2}](\lambda+\varepsilon_{-})\\\\[1.99997pt]
k_{+}[(\lambda-\varepsilon_{-})^{2}-k^{2}]\end{pmatrix}$ (8)

where $\varepsilon_{\pm}=(E_{F}\pm\varepsilon)$, $N_{e}(N_{h})$ is
normalization constant fulfilling the condition
${u_{e}^{\eta}}^{\dagger}u_{e}^{\eta}=1({v_{h}^{\eta}}^{\dagger}v_{h}^{\eta}=1)$,
and the corresponding longitudinal momenta for $e$ and $h$ are

$\displaystyle\pm k_{x}^{e\tau}$
$\displaystyle=\pm\sqrt{\varepsilon_{+}^{2}+\lambda^{2}+\tau\Sigma_{e}-k_{y}^{2}}$
(9) $\displaystyle\pm k_{x}^{h\tau}$
$\displaystyle=\pm\sqrt{\varepsilon_{-}^{2}+\lambda^{2}+\tau\Sigma_{h}-k_{y}^{2}}$
(10)

with index $\tau=\pm$ (denoting, number of incident modes) and
$\Sigma_{e(h)}=\sqrt{(4\lambda^{2}+t_{\perp}^{2})\varepsilon_{+(-)}^{2}-(\lambda
t_{\perp})^{2}}$. However, for S region, finding an explicit analytical
expression for the state vector is a substantially difficult task as the non-
zero $\Delta$ couples $u_{e}$ and $v_{h}$. Additionally, the matrix size of
the DBdG Eq. (4) for the gapped BLG is double in comparison to the monolayer
graphene Beenakker (2006). Therefore, the state vector in this region is
calculated numerically, say, $u_{S}^{\eta}(\varepsilon,k_{x})$. Nevertheless,
we deduce an analytical form for longitudinal momenta: $\\{\pm
k_{x,\pm}^{(S,1)},\pm k_{x,\pm}^{(S,2)}\\}$ with
$k_{x,\pm}^{(S,i)}=\sqrt{\Gamma_{i}^{2}+\lambda^{2}-k_{y}^{2}\pm\sqrt{(4\lambda^{2}+t_{\perp}^{2})\Gamma_{i}^{2}-(\lambda
t_{\perp})^{2}}}$ where
$\Gamma_{i}=(U_{0}+E_{F})-(-1)^{i}\sqrt{\varepsilon^{2}-\Delta^{2}}$ for
$i=1,2$.

Using these states, we construct the scattering wavefunctions for both sides
which consist Andreev and normal reflection coefficients and transmission
coefficients, and that are obtained by demanding the continuity of the
wavefunctions at the junction $x=0$. By utilizing these coefficients, we
calculate the Andreev and normal reflections probabilities and the
differential conductance in the framework of Blonder-Tinkham-Klapwijk (BTK)
formalism Blonder  _et al._ (1982). See the Appendix A and C for the technical
details that are given for the two regimes $\lambda\ll t_{\perp}$ and
$\lambda\lesssim t_{\perp}$, respectively.

Figure 2: Schematic depiction of the gapped bilayer graphene-superconductor
junction and the possible reflection processes when an electron $e_{i}$ (solid
arrow with filled green circle) is incident at an angle $\alpha$: (a)-(b)
$\lambda\ll t_{\perp}$ and (c)-(d) $\lambda\lesssim t_{\perp}$. The other
solid arrows and dotted arrows represent respectively the electron reflections
and hole reflections. In (b) and (d), the solid and dotted concentric circles
in $k_{x}$-$k_{y}$ plane denote isoenergy contours for electron and hole
excitations, see Eq. (5), and arrows at the horizontal dashed-purple-line are
the corresponding group velocities. For the incident electron $e_{i}$ in
(a)-(b), the reflected hole is RAR when $\varepsilon<E_{F}$ (dotted blue
arrow) and is SAR when $\varepsilon>E_{F}$ (dotted red arrow). The electron-
hole conversion for RAR is intraband (only the conduction band) however SAR is
due to the interband conversion (conduction band to valence band). In
contrast, double Andreev reflections (RAR and SAR) and double normal
reflections (RNR and SNR) exist simultaneously for (c)-(d) and the electron-
hole conversion is always intraband type.

In Fig. 2, we have schematically shown the possible reflection processes
involved in both regimes. For $\lambda\ll t_{\perp}$ regime, the Figs.
2(a)-(b) show the usual RAR, SAR, and SNR reflections and the involved
excitation energy contours for electron and hole that participate in the
scattering process. As $\lambda$ is small, it only opens a band gap $\sim
2\lambda$ [see in Fig. 1(iii)(b)] and shows the usual reflections Ludwig
(2007); Efetov  _et al._ (2016). However, for $\lambda\lesssim t_{\perp}$
regime in Figs. 2(c)-(d), four reflections happen simultaneously which are the
RAR, SAR, SNR, and RNR Cheng  _et al._ (2020). The extra RNR process occurs
because of the inversion of the lower conduction band [see in Fig. 1(iii)(c)]
for large $\lambda$. This supports two isoenergy-contours for electron
excitation $\varepsilon_{+,-}^{e}$ and two for hole excitation
$\varepsilon_{+,-}^{h}$ when Fermi energy $E_{F}$ is set in-between the
maximum and minimum of the lower conduction band. Consequently, the SNR and
RNR from $\varepsilon_{+,-}^{e}$ and the RAR and SAR from
$\varepsilon_{+,-}^{h}$ occur, which are shown in Fig. 2(d).

In the next section, we present the results for the reflection probabilities
and differential conductance at zero temperature and discuss them in detail.
They are calculated by using the Eqs. (12), (13), (15), and (16). For
numerical calculation, we fix $\Delta=1$ and set all energy parameters in
units of $\Delta$. Since the bulk SC pair potential can be achieved
$\Delta\sim 1.2$ meV (by depositing NbSe2 on the BLG sheet Efetov  _et al._
(2016)) and the interlayer coupling in BLG is roughly $t_{\perp}\sim 0.39$ eV,
we present all results for $t_{\perp}=400\Delta$, except in Figs. 4(c)-(d),
5(b), and 9(c)-(d).

##  IV Results

###  IV.1 Small displacement field, $\lambda\ll t_{\perp}$

In this subsection, we work in regime $(\lambda,E_{F},\varepsilon,\Delta)\ll
t_{\perp}$ which set the longitudinal momenta for $\tau=-$ mode $k_{x}^{e-}$
and $k_{x}^{h-}$ imaginary, and consequently the corresponding state vectors
in Eqs. (7) and (8) become evanescent type. So, we consider only $\tau=+$ mode
for the incident electron which governs the scattering mechanism. The possible
reflection processes are depicted in Figs. 2(a) and 2(b), and the reflection
probabilities and differential conductance formula are given in Appendix A.

####  IV.1.1 Reflection probabilities

We present the normal and Andreev reflection probabilities, $R_{n,+}^{\eta}$
and $R_{a,+}^{\eta}$, with respect to the excitation energy $\varepsilon$ and
incident angle $\alpha$. Figure 3 shows $R_{n,+}^{\eta}$ and $R_{a,+}^{\eta}$
plots at valley $\eta=K(K^{\prime})$ for parameters $E_{F}=0.5\Delta$ and
$U_{0}=10\Delta$. In Figs. 3(a)-(b), we set the displacement field
$\lambda=0$, while it is set $\lambda=0.1\Delta$ in Figs. 3(c)-(f). For the
excitation range $0<\varepsilon<1.5\Delta$, the scattering process involves
bands $\varepsilon^{e}_{+,-}$, $\varepsilon^{h}_{+,-}$, and
$\varepsilon^{h}_{-,-}$ from the N side; however, only band
$\varepsilon_{+,-}^{(S,1)}$ is involved from the S side. The incident electron
always come from $\varepsilon^{e}_{+,-}$ and the reflected hole belongs to
$\varepsilon^{h}_{+,-}$ for $0<\varepsilon<(E_{F}-\lambda)$ whereas it comes
from $\varepsilon^{h}_{-,-}$ for $\varepsilon>(E_{F}+\lambda)$. When
$\varepsilon>\Delta$, the transmitted quasiparticles (electron-like as well as
hole-like) belong to $\varepsilon_{+,-}^{(S,1)}$.

Figure 3: Heatmap plots for the normal reflection
$R_{n,+}^{\eta}(\varepsilon,\alpha)$ and Andreev reflection
$R_{a,+}^{\eta}(\varepsilon,\alpha)$ at valley $\eta\equiv K,K^{\prime}$ for
(a)-(b) $\lambda=0$ and (c)-(f) $\lambda=0.1\Delta$. The other parameters are
$E_{F}=0.5\Delta$ and $U_{0}=10\Delta$. The horizontal dashed-line is at
$\varepsilon=E_{F}$ where $R_{a,+}^{\eta}$ vanishes and $R_{n,+}^{\eta}$
becomes $1$ when $\lambda=0$, and it widens to the range
$(E_{F}-\lambda)<\varepsilon<(E_{F}+\lambda)$ of width $2\lambda$ when
$\lambda\neq 0$.

In Figs. 3(a)-(b), $R_{n,+}^{K(K^{\prime})}$ and $R_{a,+}^{K(K^{\prime})}$ are
valley degenerated and symmetric with respect to $\alpha$. When
$0<\varepsilon<E_{F}$, at a given $\varepsilon$, the $R_{n,+}^{\eta}$ is
enhanced while $R_{a,+}^{\eta}$ is weakened as $\alpha$ increases from $0$ to
$\pi/2$. The $R_{n,+}^{\eta}(R_{a,+}^{\eta})$ for $E_{F}<\varepsilon<\Delta$
starts approximately equal to $1(0)$ [see the colorbar of Fig. 3] at
$\alpha=0$ and decreases (increases) with $\alpha/\pi\lesssim 0.25$ and again
reaches to $1(0)$ for further increasing. We notice that
$R_{n,+}^{\eta}+R_{a,+}^{\eta}=1$ as no transmissions are allowed in subgap
($\varepsilon<\Delta$) region, while it is weakened
$R_{n,+}^{\eta}+R_{a,+}^{\eta}\neq 1$ in $\varepsilon>\Delta$ region because
the quasiparticles transmission also happen. When $\lambda=0.1\Delta$ in Figs.
3(c)-(f), both $R_{n,+}^{\eta}$ and $R_{a,+}^{\eta}$ follow similar behaviour
as in Figs. 3(a)-(b), but they are now slightly asymmetric about $\alpha$ for
each valley $\eta=K,K^{\prime}$. Changing the valley from $K\to K^{\prime}$,
this asymmetry is reversed, i.e.
$R_{n,+}^{K}(\alpha)=R_{n,+}^{K^{\prime}}(-\alpha)$ and
$R_{a,+}^{K}(\alpha)=R_{a,+}^{K^{\prime}}(-\alpha)$, which is a direct
consequence of the layer asymmetry as $\lambda$ is nonzero. The presence of
$\lambda$ opens a gap $\sim 2\lambda$ at $\varepsilon=E_{F}$ (cyan dashed-
line) for the incident angle $\alpha=0$ as neither $\varepsilon^{h}_{+,-}$ nor
$\varepsilon^{h}_{-,-}$ is available for the excitation range
$(E_{F}-\lambda)<\varepsilon<(E_{F}+\lambda)$, and as a result, Andreev
reflection vanishes $R_{a,+}^{\eta}=0$ and normal reflection reaches to
$R_{n,+}^{\eta}=1$.

####  IV.1.2 Differential conductance

Figure 4 shows the normalized conductance $G/G_{0}$ versus applied bias
voltage $\varepsilon=eV$ when $\lambda=0$ which is obtained by using Eq. (13).
First, we discuss the result in the limit $t_{\perp}\gg U_{0}$ in Figs.
4(a)-(b) for $t_{\perp}=400\Delta$ and $U_{0}=10\Delta$. This limit
corresponds to the work in Ref. Efetov and Efetov (2016) where authors
calculated the subgap conductance approximately. As discussed in the previous
subsection, the scattering process always involves $\varepsilon^{e}_{+,-}$ for
electron, and $\varepsilon^{h}_{+,-}$ for the reflected hole when $eV<E_{F}$
(RAR) else $\varepsilon^{h}_{-,-}$ when $eV>E_{F}$ (SAR), whereas
$\varepsilon_{+,-}^{(S,1)}$ accounts for the quasiparticles transmission when
$eV>\Delta$.

Figure 4: Normalized conductance $G/G_{0}$ vs. $eV/\Delta$, calculated from
Eq. (13), in the absence of displacement field $\lambda=0$. (a) The $G/G_{0}$
for several fixed $E_{F}$ in limit $t_{\perp}\gg U_{0}$. (b) The effect of
$U_{0}$ on $G/G_{0}$ for $E_{F}=0.5\Delta$. (c) In opposite limit
$t_{\perp}\ll U_{0}$, a dip in $G/G_{0}$ is observed near the gap
$eV\sim\Delta$ when $E_{F}<\Delta$. (d) For $E_{F}=0$, the dip near the SC gap
moves closure to $\Delta$ as the interlayer coupling $t_{\perp}$ increases.

In Fig. 4(a), we plot $G/G_{0}$ for several fixed $E_{F}$. For $E_{F}<\Delta$,
it begins with a finite value and starts decreasing as $eV$ increases and
vanishes at $eV=E_{F}$ since no Andreev reflection happens at any incident
angle. On further increasing $eV$, it again rises, exhibiting a singularity at
$eV=\Delta$ similar to the ordinary NS junction Blonder  _et al._ (1982). The
$G/G_{0}$ becomes weak when $eV>\Delta$, due to the quasiparticles
transmission. However, for $E_{F}\geq\Delta$, only RAR contributes to the
subgap conductance and reaches $G/G_{0}\to 2$ for large $E_{F}$, see the curve
for $E_{F}=10\Delta$. The Fermi energy mismatch at the junction, due to the
finite doping potential $U_{0}=10\Delta$, prevents achieving the maximum value
of $G/G_{0}=2$. In Fig. 4(b), we show $U_{0}$-dependent $G/G_{0}$ at
$E_{F}=0.5\Delta$. It shows a strong suppression on increasing $U_{0}$ as the
Fermi energy mismatch at the junction increases with $U_{0}$ which reduces the
AR, leading to the suppression. Our results are numerically exact and are in
agreement with the findings in Ref. Efetov and Efetov (2016).

Next, we set the parameters $t_{\perp}=10\Delta$ and $U_{0}=400\Delta$ and
work in the opposite limit $t_{\perp}\ll U_{0}$ to observe the dip in
$G/G_{0}$ near the gap $eV\sim\Delta$ which was attributed to the
pseudospin-$1$ effect in Ref. Ludwig (2007). Similar to the case in Figs.
4(a)-(b), the scattering process for $t_{\perp}\ll U_{0}$ involves the bands
$\varepsilon^{e}_{+,-}$ and $\varepsilon^{h}_{+,-}$ ($\varepsilon^{h}_{-,-}$)
for electron and reflected hole, but the quasiparticles now participate from
bands $\varepsilon_{+,-}^{(S,1)}$ and $\varepsilon_{+,-}^{(S,2)}$ for the
transmission. Consequently, we observe a dip in the conductance near $\Delta$,
see in Figs. 4(c)-(d). The $G/G_{0}$ plots for $E_{F}=0$, $0.5\Delta$, and
$10\Delta$ in Fig. 4(c) and for $t_{\perp}=3\Delta$ and $10\Delta$ in Fig.
4(d) are generated for the same values of parameters taken in Ludwig (2007).
We obtain exactly the same results for $G/G_{0}$ as obtained in Ref. Ludwig
(2007). This shows that our numerical method for the bilayer graphene-
superconductor junction produces exact results, and also works in the both
limits $t_{\perp}\gg U_{0}$ and $t_{\perp}\ll U_{0}$.

Figure 5: Normalized conductance $G/G_{0}$ with respect to $eV/\Delta$ at
fixed $E_{F}$ for different $\lambda$. The top panel (a) is for the limit
$t_{\perp}\gg U_{0}$ whereas the bottom panel (b) corresponds to the limit
$t_{\perp}\ll U_{0}$. The $G/G_{0}$ vanishes for the bias voltage range
$(E_{F}-\lambda)<eV<(E_{F}+\lambda)$ with gap width $2\lambda$ and shows a
strong modification in the presence of displacement field.

Now, we present the normalized conductance $G/G_{0}$ behaviour in the presence
of displacement field $\lambda$. The results are plotted in Fig. 5(a) for the
limit $t_{\perp}\gg U_{0}$ and in Fig. 5(b) for the limit $t_{\perp}\ll
U_{0}$. They are shown for different $\lambda$ at $E_{F}=0.5\Delta$. We
restrict $\lambda\leq E_{F}$. As can be seen that the $G/G_{0}$ vanishes
around $E_{F}$ for the bias range $(E_{F}-\lambda)<eV<(E_{F}+\lambda)$ at any
finite $\lambda$, creating a gap of width $\sim 2\lambda$ (solid lines) along
the applied bias voltage axis. It happens because no Andreev reflections
(neither RAR nor SAR) take place in this range as both the hole bands
$\varepsilon^{h}_{+,-}$ and $\varepsilon^{h}_{-,-}$ are absent. This strong
modification due to the finite displacement field $\lambda$ suggests that we
can tune the subgap conductance and amplify the retro-to-specular crossover
which would possibly help to realise the crossover boundary experimentally
Efetov  _et al._ (2016), since the crossover-boundary becomes wider.

###  IV.2 Large displacement field, $\lambda\lesssim t_{\perp}$

This section is devoted to studying the transport properties at the junction
in large displacement field regime, wherein the band structure illustrated in
Fig. 1(iii)(c) participates in the scattering process. In this regime, both
the $\tau=\pm$ modes with incident electron longitudinal wave-vectors
$k_{x}^{e\pm}$ are present. The possible reflection processes are
schematically shown in Figs. 2(c) and 2(d). We fix $E_{F}$ comparable to
$\lambda$, i.e. $E_{F}\simeq\lambda$, and set the parameters
$t_{\perp}=400\Delta$, $E_{F}=80\Delta$, and $U_{0}=10\Delta$ for the
calculation. We first calculate the reflection probabilities at a given
$\lambda$ with respect to $\alpha$ and $\varepsilon$. Using these
probabilities, the conductance is calculated by integrating over $\alpha$, see
the formulation and technical details in Appendix C. For brevity, we only
present the results for differential conductance and discuss the key findings
(the $\alpha$-dependent and $\varepsilon$-dependent reflection probabilities
are presented in Appendix C for completeness).

Figure 6(a) presents the normalized conductance $G/G_{0}$ versus $eV/\Delta$
for several values of $\lambda$ taken around $E_{F}$. The scattering process
involves only the conduction band excitations $\varepsilon^{e}_{+,-}$ and
$\varepsilon^{h}_{+,-}$ for electron and hole, respectively. So, the electron-
hole conversion for the Andreev reflections is always intraband in nature. To
better understand the behavior of $G/G_{0}$, we also plot quantities
$\delta_{\rm min}=E_{F}-E_{\rm min}$ and $\delta_{\rm max}=E_{\rm max}-E_{F}$
with respect to $\lambda$ in Fig. 6(b), which are depicted in the bottom panel
of Fig. 6(c). The filled and open circles, in Fig. 6(b), are at the selected
$\lambda$ points for which $G/G_{0}$ is plotted. The top panel of Fig. 6(c)
shows the situation when $E_{F}$ lies above $E_{\rm max}$ in the conduction
band, i.e, $E_{F}>E_{\rm max}$.

Figure 6: (a) Normalized conductance $G/G_{0}$, calculated by using the Eq.
(16), with $eV/\Delta$ for several $\lambda$’s. (b) Variation of $\delta_{\rm
min}=E_{F}-E_{\rm min}$ and $\delta_{\rm max}=E_{\rm max}-E_{F}$ with
$\lambda/\Delta$. The filled and empty circles correspond to the points
$(\delta_{\rm min}/\Delta,\delta_{\rm max}/\Delta)$ = $(7.33,-2)$, $(5.72,0)$,
$(5.32,0.5)$, $(4.92,1)$, $(4.13,2)$, $(2.55,4)$, $(1,6)$, $(0.61,6.5)$,
$(0.22,7)$ for the selected $\lambda$’s in (a). The schematic plot of the
conduction band in (c) is shown for the situations when $\lambda<E_{F}$ (top
panel) and $\lambda\geq E_{F}$ (bottom panel).

For $\lambda=78\Delta$, in Fig. 6(a), only the $\tau=+$ mode wave-vectors are
available for electron and hole excitations for bias $0<eV<1.6\Delta$, so the
SAR and RNR vanish as the $\tau=-$ mode is absent due to the condition
$E_{F}>E_{\rm max}$. Therefore, only the RAR and SNR contribute to the
conductance. $G/G_{0}\to 2$ in the subgap region remains almost constant
because $E_{F}\gg\Delta$ and provides a similar result as shown in Fig. 4(a)
for $E_{F}=10\Delta$ and behaves as an ordinary NS junction Blonder  _et al._
(1982). When $\lambda=80\Delta(=E_{F})$, the scattering process is still
governed by $\tau=+$ mode wave-vectors for electron excitation and thus
similar $G/G_{0}$ behavior, but now both Andreev reflections (RAR and SAR)
exist along with SNR due to the presence of $\tau=\pm$ wave-vectors for hole
excitation. Next, we notice that the $G/G_{0}$ starts with a lower value and
rises again to the previous value at the critical bias $(eV)_{c}=0.5\Delta$
for $\lambda=80.5\Delta$ (dotted green curve, $\delta_{\rm max}=0.5\Delta$)
and $(eV)_{c}=\Delta$ for $\lambda=81\Delta$ (dot-dashed purple curve,
$\delta_{\rm max}=\Delta$) because the $\tau=\pm$ modes participate for
electron and hole in range $0<eV<(eV)_{c}$. This leads to the appearance of a
new type of normal reflection—retro normal reflection (RNR). Correspondingly,
it weakens $G/G_{0}$ in the region $0<eV<(eV)_{c}$.

In contrast, when $\lambda>(E_{F}+\Delta)$, both $\tau=\pm$ modes wave-vectors
from $\varepsilon^{e}_{+,-}$ are always available for electron excitation, but
depending on the strength of $\lambda$ the hole excitation
$\varepsilon^{h}_{+,-}$ does not always participate in the scattering process
for the whole range of $eV$. For instance, $\varepsilon^{h}_{+,-}$ is always
involved at $\lambda=82\Delta$ and $84\Delta$, and thus all four reflections
contribute. The further decrease in $G/G_{0}$ within the subgap region is due
to the presence of RNR along with SNR, whereas the quasiparticle transmissions
are responsible for the decrease in the $eV>\Delta$ region. However, the hole
excitation $\varepsilon^{h}_{+,-}$ is absent when bias voltage exceeds the
critical bias $(eV)_{c}=\Delta$, $0.61\Delta$, and $0.22\Delta$ for the
displacement fields $\lambda=86\Delta$ ($\delta_{\rm min}=\Delta$),
$\lambda=86.5\Delta$ ($\delta_{\rm min}=0.61\Delta$), and $\lambda=87\Delta$
($\delta_{\rm min}=0.22\Delta$), respectively. As a result, the double Andreev
reflections SAR and RAR become zero, and the non-zero double normal
reflections SNR and RNR cause the zero conductance inside the gap for
$(eV)_{c}<eV<\Delta$. These results are consistent with the discussion of the
reflection probabilities, given in Appendix C. Overall, we observe that
$G/G_{0}$ clearly shows a distinct feature which could be experimentally
verified by varying the displacement field close to the Fermi energy of the
normal side.

###  IV.3 Experimentally-feasible results: $\lambda=0$ (in S region)

In the preceding two subsections, we examined the conductance characteristics
in the small and large displacement field regimes while keeping $\lambda$ non-
zero and equal on both the N and S sides of the junction. Nevertheless,
implementing a non-zero displacement field on the S side in an experimental
setup would likely pose considerable challenges, if not impossibilities.
Therefore, in this subsection, we focus on calculating the conductance in the
small and large displacement field regimes when $\lambda$ is turned off on the
S side; however, it is still applied on the N side.

Figure 7: The normalized conductance $G/G_{0}$ vs. $eV/\Delta$ for different
non-zero $\lambda$’s in normal region. (a) Conductance behaviour is almost
exactly same as obtained in the Fig. 5(a) for $U_{0}=10\Delta$. However, it
deviates slightly in (b) for $U_{0}=\Delta$. The other parameters are
$t_{\perp}=400\Delta$ and $E_{F}=0.5\Delta$. The subscript in $\lambda_{N(S)}$
is added to specify the displacement field $\lambda$ in N (S) region.

In Fig. 7, we show the normalized conductance $G/G_{0}$ with respect to
$eV/\Delta$ for different values of $\lambda$ in the regime $\lambda\ll
t_{\perp}$. The other fixed parameters are given in the figure panels and
caption. Notice that the conductance curves in Fig. 7(a) exhibit almost
exactly the same result as in Fig. 5(a) even though $\lambda=0$ in the S
region. This happens mainly because the Fermi energy mismatch at the junction
does not change significantly since $\lambda\ll U_{0}$ as $U_{0}=10\Delta$ in
the both calculation. In contrast, this difference is visible when we set
$U_{0}=\Delta$ in Fig. 7(a) as $\lambda\sim U_{0}$. However, the qualitative
behaviour is still similar.

Now, we present the $G/G_{0}$ behaviour in the large displacement field limit,
see Fig. 8, for the same set of parameters values as in Fig. 6, except here
the displacement field $\lambda$ is non-zero only on N side, i.e.,
$\lambda_{N}=\lambda$ and $\lambda_{S}=0$. We see that the qualitative
characteristics of $G/G_{0}$ is similar to the Fig. 6(a), but now the subgap
conductance has dropped significantly as the Fermi-energy mismatch is enlarged
because of $\lambda_{S}=0$. Consequently, the Andreev reflections (normal
reflections) contribution to the conductance would be decreased (increased),
reducing the subgap conductance.

Figure 8: Normalized conductance $G/G_{0}$ vs. $eV/\Delta$, in large
displacement field regime $\lambda\lesssim t_{\perp}$, for several
$\lambda_{N}$ values while $\lambda_{S}=0$. The other parameters are
$t_{\perp}=400\Delta$, $E_{F}=80\Delta$, and $U_{0}=10\Delta$.

##  V Conclusion

We have studied the scattering reflection processes at the gapped bilayer
graphene-superconducting junction by employing the DBdG equation within the
scattering theory approach. Since $\lambda$ provides an external tool to tune
the BLG band structure, we have thoroughly investigated its effect on the
normal reflections, the Andreev reflections, and the experimentally accessible
quantity the differential conductance in two regimes: $\lambda\ll t_{\perp}$
and $\lambda\lesssim t_{\perp}$. In the former regime, we have revisited the
$\lambda=0$ case first and obtained the expected specular (retro) Andreev
reflection for the excitation energy $\varepsilon>E_{F}$
($\varepsilon<E_{F}$); the retro-to-specular crossing point at bias $eV=E_{F}$
in the subgap conductance Ludwig (2007); Efetov and Efetov (2016). As our
numerical calculation provides the exact result, we also observe the dip in
conductance near $eV\sim\Delta$ when $t_{\perp}\ll U_{0}$, which was
attributed to the pseudospin-$1$ feature in Ref. Ludwig (2007). However, this
dip does not appear when $t_{\perp}\gg U_{0}$, as is also the case in Ref.
Efetov and Efetov (2016). Furthermore, introducing a non-zero small $\lambda$
broadens the retro-to-specular crossing for the bias range
$|E_{F}-\lambda|<eV<|E_{F}+\lambda|$ in the subgap conductance which is a
direct manifestation of the band gap between conduction and valence bands due
to finite displacement field. This suggests that tuning $\lambda$ the
crossover region can be easily modified and possibly obtained experimentally
as the authors in Efetov  _et al._ (2016) struggle to observe it.

However, for the latter regime $\lambda\lesssim t_{\perp}$, we have shown that
apart from SNR, SAR, and RAR, there exists also RNR due to the presence of the
Mexican-hat-shape band structure in the gapped BLG. This is a new finding in a
bilayer graphene NS junction which has never been explored as per our
knowledge. The normalized conductance shows a very distinct characteristic
feature when all four reflections contribute, see Fig. 6, and it can be
distinguished from the result when the Mexican-hat structure does not affect
the scattering process.

###### Acknowledgements.

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - 467596333. This work was partly supported by the
Helmholtz Association through program NACIP.

##  Appendix A Small displacement field ($\lambda\ll t_{\perp}$)

In this regime, we have only the $\tau=+$ mode with the incident electron
wave-vector $(k_{x}^{e+},k_{y}^{+})$; therefore, the scattering wavefunctions
can be expressed as $\Psi^{\eta}(x)e^{ik^{+}_{y}y}$ where

$\displaystyle\Psi^{\eta}(x)=$
$\displaystyle\big{[}u_{e}^{\eta}(\varepsilon,k_{x}^{e+})e^{ik_{x}^{e+}x}+r_{n,+}^{\eta}u_{e}^{\eta}(\varepsilon,-k_{x}^{e+})e^{-ik_{x}^{e+}x}$
$\displaystyle~{}+r_{n,-}^{\eta}u_{e}^{\eta}(\varepsilon,-i\kappa_{x}^{e-})e^{\kappa_{x}^{e-}x}\big{]}\begin{pmatrix}1\\\\[0.50003pt]
0\end{pmatrix}$
$\displaystyle~{}+\big{[}r_{a,+}^{\eta}v_{h}^{\eta}(\varepsilon,\beta
k_{x}^{h+})e^{i\beta k_{x}^{h+}x}$
$\displaystyle~{}+r_{a,-}^{\eta}v_{h}^{\eta}(\varepsilon,-i\kappa_{x}^{h-})e^{\kappa_{x}^{h-}x}\big{]}\begin{pmatrix}0\\\\[0.50003pt]
1\end{pmatrix},\quad\mathrm{for}~{}x\leq 0$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{4}t_{j}^{\eta}u_{S}^{\eta}(\varepsilon,k_{x,j}^{S})e^{ik_{x,j}^{S}x},\quad\mathrm{for}~{}x\geq
0$ (11)

Here, $k_{x}^{e-}(k_{x}^{h-})=-i\kappa_{x}^{e-}(-i\kappa_{x}^{h-})$ and
$k_{x,j}^{S}$ for $j=1,...,4$, are the wave-vectors of side S and are chosen
appropriately from $\\{\pm k_{x,\pm}^{(S,1)},\pm k_{x,\pm}^{(S,2)}\\}$. The
coefficients $r_{n,\pm}^{\eta}$, $r_{a,\pm}^{\eta}$, and $t_{j=1,..,4}$ are
normal reflections, Andreev reflections, and transmissions, respectively.
$r_{n,-}^{\eta}$ and $r_{a,-}^{\eta}$ are always zero because of evanescent
solutions for the $\tau=-$ mode. It is worth mentioning that $r_{n,+}^{\eta}$
is always SNR, but $r_{a,+}^{\eta}$ is RAR for $\varepsilon<(E_{F}-\lambda)$,
$\beta=1$ and SAR for $\varepsilon>(E_{F}+\lambda)$, $\beta=-1$. By using the
continuity condition, $\Psi^{\eta}|_{x=0^{-}}=\Psi^{\eta}|_{x=0^{+}}$, we
obtain these coefficients for a given $\varepsilon$ and incident angle
$\alpha$ defined from
$k^{+}_{y}=\sqrt{(E_{F}+\varepsilon)^{2}+\lambda^{2}+\Sigma_{e}}\sin\alpha$
where $\alpha\in[-\pi/2,\pi/2]$. So, the reflection probabilities are obtained
as:

$\displaystyle R_{n,+}^{\eta}$
$\displaystyle=|r_{n,+}^{\eta}(\varepsilon,\alpha)|^{2},~{}~{}R_{a,+}^{\eta}=\frac{V_{h+}^{\eta}}{V_{e+}^{\eta}}|r_{a,+}^{\eta}(\varepsilon,\alpha)|^{2}$
(12)

with velocities $V_{e+}^{\eta}=\langle
u_{e}^{\eta}(\varepsilon,k_{x}^{e+})|\frac{\partial H_{\eta}}{\partial
k_{x}}|u_{e}^{\eta}(\varepsilon,k_{x}^{e+})\rangle$ and $V_{h+}^{\eta}=\langle
v_{h}^{\eta}(\varepsilon,k_{x}^{h+})|\frac{\partial H_{\eta}}{\partial
k_{x}}|v_{h}^{\eta}(\varepsilon,k_{x}^{h+})\rangle$. The differential
conductance within the BTK framework Blonder  _et al._ (1982) at zero
temperature can be expressed as follows:

$\displaystyle
G(\varepsilon)=\sum_{\eta}g_{0}(\varepsilon)\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\big{(}1-R_{n,+}^{\eta}+R_{a,+}^{\eta}\big{)}\cos\alpha~{}d\alpha$
(13)

with $\varepsilon=eV$ and
$g_{0}(\varepsilon)=\frac{2e^{2}W}{h\pi}\sqrt{(E_{F}+\varepsilon)^{2}+\lambda^{2}+\Sigma_{e}}$
where $W$ is the width of BLG sheet and factor $2$ accounts for the spin
degeneracy. Also, the differential conductance for a normal-to-normal BLG
junction is $G_{0}(\varepsilon)=4g_{0}(\varepsilon)$ where the factor 4 is due
to sum of valley index $\eta=\pm$ and incident angle $\alpha\in[-\pi/2,\pi/2]$
instead of $\alpha\in[0,\pi/2]$ Beenakker (2006); Ludwig (2007); Efetov and
Efetov (2016).

Figure 9: Behavior of the normalized conductance, $G/G_{0}$, with respect to
$eV/\Delta$ for three values of $\lambda$ at $E_{F}=0.2\Delta$ and
$0.8\Delta$. In (a)-(b) $t_{\perp}\gg U_{0}$ for $t_{\perp}=400\Delta$ and
$U_{0}=10\Delta$ whereas in (c)-(d) $t_{\perp}\ll U_{0}$ for
$t_{\perp}=10\Delta$ and $U_{0}=400\Delta$.

##  Appendix B Effect of $\lambda$ on $G/G_{0}$ at different $E_{F}$

In order to be consistent with the main text in subsection IV.1, we illustrate
the effect of the displacement field $\lambda$ on the normalized conductance
$G/G_{0}$ for two values of Fermi energy $E_{F}=0.2\Delta$ and $0.8\Delta$ to
study the contribution of RAR and SAR in subgap region $eV<\Delta$. The
results are shown in Figs. 9(a)-(b) for the limit $t_{\perp}\gg U_{0}$ and in
Figs. 9(c)-(d) for the limit $t_{\perp}\ll U_{0}$. In both limits, when Fermi
energy is set at $E_{F}=0.2\Delta$, the contribution of the RAR and SAR
equally reduced about $E_{F}$ and eventually only SAR contributes to $G/G_{0}$
as $\lambda$ is increased to the value $\lambda=E_{F}$, see in Figs. 9(a) and
9(c). While for $E_{F}=0.8\Delta$ in Figs. 9(b) and 9(d), only the RAR
contributes. Hence, varying the displacement field and setting the Fermi
energy appropriately, the Andreev reflection contribution to the subgap
conductance can be tuned to either SAR or RAR.

##  Appendix C Large displacement field ($\lambda\lesssim t_{\perp}$)

In this regime, we have both the $\tau=\pm$ modes for the incident electron
with the wave-vectors $(k_{x}^{e+},k^{+}_{y})$ and $(-k_{x}^{e-},k^{-}_{y})$
where
$k_{y}^{\pm}=\sqrt{(E_{F}+\varepsilon)^{2}+\lambda^{2}\pm\Sigma_{e}}\sin\alpha$.
So, the scattering wavefunctions on both sides of the junction have the form
$\Psi_{\tau}^{\eta}(x)e^{ik^{\tau}_{y}y}$ with

$\displaystyle\Psi_{\tau}^{\eta}(x)=$
$\displaystyle~{}\big{[}u_{e}^{\eta}(\varepsilon,\tau k_{x}^{e\tau})e^{i\tau
k_{x}^{e\tau}x}+r_{n1,\tau}^{\eta}u_{e}^{\eta}(\varepsilon,-k_{x}^{e+})e^{-ik_{x}^{e+}x}$
$\displaystyle~{}+r_{n2,\tau}^{\eta}u_{e}^{\eta}(\varepsilon,k_{x}^{e-})e^{ik_{x}^{e-}x}\big{]}\begin{pmatrix}1\\\\[0.50003pt]
0\end{pmatrix}$
$\displaystyle~{}+\big{[}r_{a1,\tau}^{\eta}v_{h}^{\eta}(\varepsilon,-k_{x}^{h-})e^{-ik_{x}^{h-}x}$
$\displaystyle~{}+r_{a2,\tau}^{\eta}v_{h}^{\eta}(\varepsilon,k_{x}^{h+})e^{ik_{x}^{h+}x}\big{]}\begin{pmatrix}0\\\\[0.50003pt]
1\end{pmatrix},\quad\mathrm{for}~{}x\leq 0$ $\displaystyle=$
$\displaystyle~{}\sum_{j=1}^{4}t_{j,\tau}^{\eta}u_{S}^{\eta}(\varepsilon,k_{x,j,\tau}^{S})e^{ik_{x,j,\tau}^{S}x},\quad\mathrm{for}~{}x\geq
0$ (14)

where $r_{n1,\tau}^{\eta},r_{n2,\tau}^{\eta},r_{a1,\tau}^{\eta}$, and
$r_{a2,\tau}^{\eta}$ denote reflection coefficients, namely, SNR, RNR, SAR,
and RAR for $\tau=+$ mode; however, they become RNR, SNR, RAR, and SAR for
$\tau=-$ mode. The $t^{\eta}_{j,\tau}$ for $j=1,...,4$ are transmission
coefficients. Again, using the continuity condition, we obtain the reflection
probabilities as:

$\displaystyle R_{n1,\tau}^{\eta}$
$\displaystyle=\frac{V_{e+}^{\eta}}{V_{e\tau}^{\eta}}|r_{n1,\tau}^{\eta}(\varepsilon,\alpha)|^{2},~{}~{}R_{n2,\tau}^{\eta}=\frac{V_{e-}^{\eta}}{V_{e\tau}^{\eta}}|r_{n2,\tau}^{\eta}(\varepsilon,\alpha)|^{2}$
$\displaystyle R_{a1,\tau}^{\eta}$
$\displaystyle=\frac{V_{h-}^{\eta}}{V_{e\tau}^{\eta}}|r_{a1,\tau}^{\eta}(\varepsilon,\alpha)|^{2},~{}~{}R_{a2,\tau}^{\eta}=\frac{V_{h+}^{\eta}}{V_{e\tau}^{\eta}}|r_{a2,\tau}^{\eta}(\varepsilon,\alpha)|^{2}$
(15)

with velocities $V_{e\tau}^{\eta}=\langle u_{e}^{\eta}(\varepsilon,\tau
k_{x}^{e\tau})|\frac{\partial H_{\eta}}{\partial
k_{x}}|u_{e}^{\eta}(\varepsilon,\tau k_{x}^{e\tau})\rangle$ and
$V_{h\tau}^{\eta}=\langle v_{h}^{\eta}(\varepsilon,\tau
k_{x}^{h\tau})|\frac{\partial H_{\eta}}{\partial
k_{x}}|v_{h}^{\eta}(\varepsilon,\tau k_{x}^{h\tau})\rangle$. The differential
conductance formula in this regime becomes

$\displaystyle
G(\varepsilon)=\sum_{\eta,\tau}g_{0}^{\tau}(\varepsilon)\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\Big{(}1$
$\displaystyle-R_{n1,\tau}^{\eta}-R_{n2,\tau}^{\eta}+R_{a1,\tau}^{\eta}$
$\displaystyle+R_{a2,\tau}^{\eta}\Big{)}\cos\alpha~{}d\alpha$ (16)

where $\varepsilon=eV$ and
$g^{\tau}_{0}(\varepsilon)=\frac{2e^{2}W}{h\pi}\sqrt{(E_{F}+\varepsilon)^{2}+\lambda^{2}+\tau\Sigma_{e}}$
which leads to $G_{0}(\varepsilon)=4\sum_{\tau}g^{\tau}_{0}(\varepsilon)$ for
a gapped BLG normal-to-normal junction.  The results for the differential
conductance have been presented in the subsections. IV.2 and IV.3 in this
regime. Here, we briefly discuss the reflection probabilities for the
displacement field values $\lambda=84\Delta$ and $86.5\Delta$.

Figure 10: Reflection probabilities (at valley $K$ and $\tau=+$) vs. $\alpha$
for $E_{F}=80\Delta$, $t_{\perp}=400\Delta$, and $U_{0}=10\Delta$ with
excitation energy $\varepsilon=0.3\Delta$ in (a) and (d),
$\varepsilon=0.8\Delta$ in (b) and (e), and $\varepsilon=1.2\Delta$ in (c) and
(f). The displacement fields for left and right panels are $\lambda=84\Delta$
and $86.5\Delta$. Figure 11: Reflection probabilities (at valley $K$ and
$\tau=+$) vs. $\varepsilon$ for $E_{F}=80\Delta$, $t_{\perp}=400\Delta$, and
$U_{0}=10\Delta$ with incident angle $\alpha=0$ in (a) and (d),
$\alpha=0.25\pi$ in (b) and (e), and $\alpha=0.4\pi$ in (c) and (f). The
displacement fields for left and right panels are $\lambda=84\Delta$ and
$86.5\Delta$.

###  C.1 $\alpha$-dependent normal and Andreev reflections

We show the reflection probabilities as a function of $\alpha$ at three
$\varepsilon=0.3\Delta$, $0.8\Delta$, and $1.2\Delta$ for $\lambda=84\Delta$
and $86.5\Delta$ in Figs. 10(a)-(c) and 10(d)-(f), respectively. Since
$\lambda$ is non-zero, these probabilities are slightly asymmetric about
$\alpha$. When $\lambda=84\Delta$, the electron and hole excitations span
$\varepsilon^{e}_{+,-}\in[-2.55\Delta,\infty]$ and
$\varepsilon^{h}_{+,-}\in[-\infty,2.55\Delta]$; the double normal reflections
(SNR and RNR) and double Andreev reflections (SAR and RAR) exist for the
intersection of $-2.55\Delta\lesssim\varepsilon^{e}_{+,-}<4\Delta$ and
$-4\Delta<\varepsilon^{h}_{+,-}\lesssim 2.55\Delta$. All four probabilities
are non-zero only around $\alpha=0$ at $\varepsilon=0.3\Delta$ and
$0.8\Delta$. The RAR $R_{a2,+}^{K}$ is dominating the scattering process. As
$\alpha$ is increased from $0$ to $\pi/2$, the $R_{a2,+}^{K}$ enhances first,
reaching to $R_{a2,+}^{K}\approx 1$, and then starts decreasing and eventually
goes to zero, while RNR $R_{n2,+}^{K}$ and SAR $R_{a1,+}^{K}$ weaken and
slowly vanish, see the results in Figs. 10(a) and 10(b). Moreover, the SNR
$R_{n1,+}^{K}$ steadily decreases and almost vanishes before it strikingly
rises to the value $1$. In Fig. 10(c) at $\varepsilon=1.2\Delta$, the
magnitude of all reflections is weakened because of the quasiparticles
transmission for $\varepsilon>\Delta$. Notice that the sum of all reflection
probabilities equals $1$ for $\varepsilon<\Delta$ and becomes less than $1$
for $\varepsilon>\Delta$.

For $\lambda=86.5\Delta$, in Figs. 10(d)-(f), electron and hole excitations
now span $\varepsilon^{e}_{+,-}\in[-0.61\Delta,\infty]$ and
$\varepsilon^{h}_{+,-}\in[-\infty,0.61\Delta]$, and all the four probabilities
exist only for the intersection of
$-0.61\Delta\lesssim\varepsilon^{e}_{+,-}<6.5\Delta$ and
$-6.5\Delta<\varepsilon^{h}_{+,-}\lesssim 0.61\Delta$, see the results in Fig.
10(d) at $\varepsilon=0.3\Delta$. Interestingly, RAR $R_{a2,+}^{K}$ is
suppressed while RNR $R_{n2,+}^{K}$ is now dominating the scattering process
around $\alpha=0$. For $\varepsilon>0.61\Delta$ in Figs. 10(e) and 10(f), the
hole excitation $\varepsilon^{h}_{+,-}$ is absent so the double Andreev
reflections SAR $R_{a1,+}^{K}$ and RAR $R_{a2,+}^{K}$ vanish, but the double
normal reflections SNR $R_{n1,+}^{K}$ and RNR $R_{n2,+}^{K}$ still exist as
the electron excitation $\varepsilon^{e}_{+,-}$ is always present.

###  C.2 $\varepsilon$-dependent normal and Andreev reflections

Here, we discuss these probabilities versus $\varepsilon$, see in Fig. 11, at
three incident angles $\alpha=0$, $0.25\pi$, and $0.4\pi$. When
$\lambda=84\Delta$, both $\varepsilon^{e}_{+,-}$ and $\varepsilon^{h}_{+,-}$
are present for whole range $0<\varepsilon<1.5\Delta$ therefore all four
reflection probabilities are present at $\alpha=0$ [see Fig. 11(a)]. The SNR
$R_{n1,+}^{K}$, SAR $R_{a1,+}^{K}$, and RAR $R_{a2,+}^{K}$ are slightly
enhanced on increasing $\varepsilon$, and then they are suppressed after
$\varepsilon>\Delta$ due to the quasiparticles transmission, whereas the RNR
$R_{n2,+}^{K}$ is opposite in nature. At $\alpha=0.25\pi$ in Fig. 11(b), the
RAR starts with $R_{a2,+}^{K}=1$ as other reflections are zero, and it
gradually decreases due to small rise in SNR around $\varepsilon=\Delta$ and
then falls rapidly. However, at $\alpha=0.4\pi$, RNR and SAR are always zero
while RAR (SNR) decreases (increases) with $\varepsilon$ but RAR vanishes for
$\varepsilon<\Delta$ and SNR attains unity, see Fig. 11(c). The critical value
of $\varepsilon$ at which RAR vanishes is $\varepsilon_{c}\approx 0.58\Delta$
as the critical angle for retro-reflected hole becomes
$\alpha_{c}\equiv\arcsin{(k^{h+}/k^{e+})}=0.4\pi$ where $k^{h+}$ and $k^{e+}$
are the magnitude of hole and electron momenta.

When we set $\lambda=86.5\Delta$, the hole excitation $\varepsilon^{h}_{+,-}$
is absent for $\varepsilon>0.61\Delta$, and consequently, both Andreev
reflections SAR and RAR are zero. However, they are non-zero for
$\varepsilon<0.61\Delta$ and are decreasing slowly as $\varepsilon$ is
increased, see in Fig. 11(d) for $\alpha=0$. Concerning the double normal
reflections, they are non-zero and show non-trivial behaviour. The Fig. 11(e)
for $\alpha=0.25\pi$ shows similar results except RNR vanishes before Andreev
reflections become zero and SNR attain unity for
$0.61\Delta<\varepsilon<\Delta$. For $\alpha=0.4\pi$ in Fig. 11(f), the RAR
becomes zero at critical $\varepsilon_{c}\approx 0.24\Delta$, and SNR
$R_{n1,+}=1$ for $\varepsilon_{c}<\varepsilon<\Delta$ as the other two
reflections RNR and SAR are always zero. We find that the total probability is
equal to $1$ for $\varepsilon<\Delta$, while it is less than $1$ for
$\varepsilon>\Delta$.

## References

  * Andreev (1964) A F Andreev, “Thermal conductivity of the intermediate state of superconductors,” J. Exp. Theor. Phys. 19, 1228 (1964).
  * Tinkham (2004) Michael Tinkham,  _Introduction to superconductivity_ (Courier Corporation, 2004).
  * Bardeen  _et al._ (1957) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of superconductivity,” Phys. Rev. 108, 1175–1204 (1957).
  * Blonder  _et al._ (1982) G. E. Blonder, M. Tinkham, and T. M. Klapwijk, “Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion,” Phys. Rev. B 25, 4515–4532 (1982).
  * Shelankov (1982) A. L. Shelankov, “Thermal conductivity of the intermediate state of superconductors,” JETP Lett. 32, 111 (1982).
  * Benistant  _et al._ (1983) P. A. M. Benistant, H. van Kempen, and P. Wyder, “Direct observation of Andreev reflection,” Phys. Rev. Lett. 51, 817–820 (1983).
  * Kastalsky  _et al._ (1991) A. Kastalsky, A. W. Kleinsasser, L. H. Greene, R. Bhat, F. P. Milliken, and J. P. Harbison, “Observation of pair currents in superconductor-semiconductor contacts,” Phys. Rev. Lett. 67, 3026–3029 (1991).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, “Two-dimensional gas of massless dirac fermions in graphene,” Nature 438, 197–200 (2005).
  * Zhang  _et al._ (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, “Experimental observation of the quantum Hall effect and Berry’s phase in graphene,” Nature 438, 201–204 (2005).
  * Novoselov  _et al._ (2006) K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A. K. Geim, “Unconventional quantum Hall effect and Berry’s phase of 2$\pi$ in bilayer graphene,” Nature Physics 2, 177–180 (2006).
  * Katsnelson  _et al._ (2006) M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, “Chiral tunnelling and the Klein paradox in graphene,” Nature Physics 2, 620–625 (2006).
  * Beenakker (2006) C. W. J. Beenakker, “Specular Andreev reflection in graphene,” Phys. Rev. Lett. 97, 067007 (2006).
  * Beenakker (2008) C. W. J. Beenakker, “Colloquium: Andreev reflection and klein tunneling in graphene,” Rev. Mod. Phys. 80, 1337–1354 (2008).
  * Bhattacharjee and Sengupta (2006) S. Bhattacharjee and K. Sengupta, “Tunneling conductance of graphene nis junctions,” Phys. Rev. Lett. 97, 217001 (2006).
  * Linder and Sudbø (2007) J. Linder and A. Sudbø, “Dirac fermions and conductance oscillations in $s$\- and $d$-wave superconductor-graphene junctions,” Phys. Rev. Lett. 99, 147001 (2007).
  * Linder and Sudbø (2008) J. Linder and A. Sudbø, “Tunneling conductance in $s$\- and $d$-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism,” Phys. Rev. B 77, 064507 (2008).
  * Zhang  _et al._ (2008) Q. Zhang, D. Fu, B. Wang, R. Zhang, and D. Y. Xing, “Signals for specular Andreev reflection,” Phys. Rev. Lett. 101, 047005 (2008).
  * Benjamin and Pachos (2008) C. Benjamin and J. K. Pachos, “Detecting entangled states in graphene via crossed Andreev reflection,” Phys. Rev. B 78, 235403 (2008).
  * Majidi and Zareyan (2012) L. Majidi and M. Zareyan, “Enhanced Andreev reflection in gapped graphene,” Phys. Rev. B 86, 075443 (2012).
  * Popinciuc  _et al._ (2012) M. Popinciuc, V. E. Calado, X. L. Liu, A. R. Akhmerov, T. M. Klapwijk, and L. M. K. Vandersypen, “Zero-bias conductance peak and Josephson effect in graphene-nbtin junctions,” Phys. Rev. B 85, 205404 (2012).
  * Han  _et al._ (2018) T. Han, J. Shen, N. F. Q. Yuan, J. Lin, Z. Wu, Y. Wu, S. Xu, L. An, G. Long, Y. Wang, R. Lortz, and N. Wang, “Investigation of the two-gap superconductivity in a few-layer NbSe2-graphene heterojunction,” Phys. Rev. B 97, 060505 (2018).
  * Pandey  _et al._ (2019) P. Pandey, R. Kraft, R. Krupke, D. Beckmann, and R. Danneau, “Andreev reflection in ballistic normal metal/graphene/superconductor junctions,” Phys. Rev. B 100, 165416 (2019).
  * Bhandari  _et al._ (2020) S. Bhandari, G.-H. Lee, K. Watanabe, T. Taniguchi, P. Kim, and R.M. Westervelt, “Imaging Andreev reflection in graphene,” Nano Letters 20, 4890–4894 (2020).
  * Pandey  _et al._ (2021) P. Pandey, R. Danneau, and D. Beckmann, “Ballistic graphene Cooper pair splitter,” Phys. Rev. Lett. 126, 147701 (2021).
  * Pandey  _et al._ (2022) P. Pandey, D. Beckmann, and R. Danneau, “Energy distribution controlled ballistic josephson junction,” Phys. Rev. B 106, 214503 (2022).
  * Jois  _et al._ (2023) S. Jois, J.L. Lado, G. Gu, Q. Li, and J.U. Lee, “Andreev reflection and Klein tunneling in high-temperature superconductor-graphene junctions,” Phys. Rev. Lett. 130, 156201 (2023).
  * Sahu  _et al._ (2016) M. R. Sahu, P. Raychaudhuri, and A. Das, “Andreev reflection near the Dirac point at the graphene-NbSe2 junction,” Phys. Rev. B 94, 235451 (2016).
  * Efetov  _et al._ (2016) D. K. Efetov, L. Wang, C. Handschin, K. B. Efetov, J. Shuang, R. Cava, T. Taniguchi, K. Watanabe, J. Hone, C. R. Dean, and P. Kim, “Specular interband Andreev reflections at van der Waals interfaces between graphene and NbSe2,” Nature Physics 12, 328–332 (2016).
  * Efetov and Efetov (2016) D. K. Efetov and K. B. Efetov, “Crossover from retro to specular Andreev reflections in bilayer graphene,” Phys. Rev. B 94, 075403 (2016).
  * Ludwig (2007) T. Ludwig, “Andreev reflection in bilayer graphene,” Phys. Rev. B 75, 195322 (2007).
  * Takane  _et al._ (2017) Y. Takane, K. Yarimizu, and A. Kanda, “Andreev reflection in a bilayer graphene junction: Role of spatial variation of the charge neutrality point,” Journal of the Physical Society of Japan 86, 064707 (2017).
  * Soori  _et al._ (2018) A. Soori, M. R. Sahu, A. Das, and S. Mukerjee, “Enhanced specular Andreev reflection in bilayer graphene,” Phys. Rev. B 98, 075301 (2018).
  * Yan  _et al._ (2011) K. Yan, H. Peng, Y. Zhou, H. Li, and Z. Liu, “Formation of bilayer bernal graphene: Layer-by-layer epitaxy via chemical vapor deposition,” Nano Letters 11, 1106–1110 (2011).
  * Brown  _et al._ (2012) L. Brown, R.t Hovden, P. Huang, M. Wojcik, D. A. Muller, and J. Park, “Twinning and twisting of tri- and bilayer graphene,” Nano Letters 12, 1609–1615 (2012).
  * Mostaani  _et al._ (2015) E. Mostaani, N. D. Drummond, and V. I. Fal’ko, “Quantum Monte Carlo calculation of the binding energy of bilayer graphene,” Phys. Rev. Lett. 115, 115501 (2015).
  * McCann and Fal’ko (2006) E. McCann and V. I. Fal’ko, “Landau-level degeneracy and quantum Hall effect in a graphite bilayer,” Phys. Rev. Lett. 96, 086805 (2006).
  * McCann (2006) Edward McCann, “Asymmetry gap in the electronic band structure of bilayer graphene,” Phys. Rev. B 74, 161403 (2006).
  * Barbier  _et al._ (2009) M. Barbier, P. Vasilopoulos, F. M. Peeters, and J. M. Pereira, “Bilayer graphene with single and multiple electrostatic barriers: Band structure and transmission,” Phys. Rev. B 79, 155402 (2009).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. 81, 109–162 (2009).
  * McCann and Koshino (2013) E. McCann and M. Koshino, “The electronic properties of bilayer graphene,” Reports on Progress in Physics 76, 056503 (2013).
  * Kraft  _et al._ (2018a) R. Kraft, I. V. Krainov, V. Gall, A. P. Dmitriev, R. Krupke, I. V. Gornyi, and R. Danneau, “Valley subband splitting in bilayer graphene quantum point contacts,” Phys. Rev. Lett. 121, 257703 (2018a).
  * Wu  _et al._ (2020) Xiuqiang Wu, Hao Meng, Fanjie Kong, Haiyang Zhang, Yujie Bai, and Ning Xu, “Tunable nonlocal valley-entangled cooper pair splitter realized in bilayer-graphene van der waals spin valves,” Phys. Rev. B 101, 125406 (2020).
  * Ohta  _et al._ (2006) T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Rotenberg, “Controlling the electronic structure of bilayer graphene,” Science 313, 951–954 (2006).
  * Oostinga  _et al._ (2008) Jeroen B. Oostinga, Hubert B. Heersche, Xinglan Liu, Alberto F. Morpurgo, and Lieven M. K. Vandersypen, “Gate-induced insulating state in bilayer graphene devices,” Nature Materials 7, 151–157 (2008).
  * Zhang  _et al._ (2009) Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, “Direct observation of a widely tunable bandgap in bilayer graphene,” Nature 459, 820–823 (2009).
  * Taychatanapat and Jarillo-Herrero (2010) T. Taychatanapat and P. Jarillo-Herrero, “Electronic transport in dual-gated bilayer graphene at large displacement fields,” Phys. Rev. Lett. 105, 166601 (2010).
  * Varlet  _et al._ (2014) A. Varlet, M.-H. Liu, V. Krueckl, D. Bischoff, P. Simonet, K. Watanabe, T. Taniguchi, K. Richter, K. Ensslin, and T. Ihn, “Fabry-Pérot interference in gapped bilayer graphene with broken anti-Klein tunneling,” Phys. Rev. Lett. 113, 116601 (2014).
  * Kraft  _et al._ (2018b) R. Kraft, J. Mohrmann, R. Du, P. B. Selvasundaram, M. Irfan, U. N. Kanilmaz, F. Wu, D. Beckmann, H. von Löhneysen, R. Krupke, A. Akhmerov, I. Gornyi, and R. Danneau, “Tailoring supercurrent confinement in graphene bilayer weak links,” Nat. Commun. 9, 1722 (2018b).
  * Du  _et al._ (2018) R. Du, M.-H. Liu, J. Mohrmann, F. Wu, R. Krupke, H. von Löhneysen, K. Richter, and R. Danneau, “Tuning anti-Klein to Klein tunneling in bilayer graphene,” Phys. Rev. Lett. 121, 127706 (2018).
  * Volkov  _et al._ (1995) A.F. Volkov, P.H.C. Magnée, B.J. van Wees, and T.M. Klapwijk, “Proximity and josephson effects in superconductor-two-dimensional electron gas planar junctions,” Physica C: Superconductivity 242, 261–266 (1995).
  * de Gennes (1966) P. G. de Gennes,  _Superconductivity of Metals and Alloys_ (Benjamin, 1966).
  * Cheng  _et al._ (2020) Qiang Cheng, Zhe Hou, and Qing-Feng Sun, “Double andreev reflections and double normal reflections in nodal-line semimetal-superconductor junctions,” Phys. Rev. B 101, 094508 (2020).