Anisotropic Andreev Reflection and Josephson Effect in Ballistic Phosphorene

TL;DR

This paper investigates how anisotropic band structure in ballistic phosphorene affects Andreev reflection and Josephson supercurrent, revealing orientation-dependent behavior and control via gate voltage.

Contribution

It demonstrates the anisotropic effects on supercurrent and Andreev reflection in phosphorene, highlighting the role of transverse modes and orientation dependence.

Findings

Supercurrent varies significantly with direction in phosphorene.

Gate voltage controls Andreev reflection probability.

Supercurrent oscillates with junction length and chemical potential.

Abstract

We study Andreev reflection and the Josephson effect in a ballistic monolayer of black phosphorous, known as phosphorene. Due to the anisotropic band structure of this system, the supercurrent changes with an order of magnitude when comparing tunneling along two perpendicular directions in the monolayer. We show that the main reason for this effect is a large difference in the number of transverse modes in Andreev bound states. The oscillatory behavior of the supercurrent as a function of the length and chemical potential of the junction also differs substantially depending on the orientation of the superconducting electrodes deposited on the phosphorene sheet. For Andreev reflection, we show that a gate voltage controls the probability of this process and that the anisotropic behavior found in the supercurrent case is also present for conductance spectra.

Equations

$$H$$

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

$${\psi }_{\mathbf{k}}^{†}=[c_{\mathbf{k},1}^{†},c_{\mathbf{k},2}^{†}]$$

$$f_{\mathbf{k}}$$

$$g_{\mathbf{k}}$$

$$\displaystyle+2t_{3}\cos\Big{(}\frac{5k_{x}}{2\sqrt{3}}\Big{)}\cos\Big{(}\frac{k_{y}}{2}\Big{)}+t_{5}\cos\Big{(}\frac{2k_{x}}{\sqrt{3}}\Big{)},$$

$$h_{\mathbf{k}}$$

$$\displaystyle+2t_{3}\sin\Big{(}\frac{5k_{x}}{2\sqrt{3}}\Big{)}\cos\Big{(}\frac{k_{y}}{2}\Big{)}-t_{5}\sin\Big{(}\frac{2k_{x}}{\sqrt{3}}\Big{)}.$$

$${ϵ}_{\pm ,\mathbf{k}}=f_{\mathbf{k}}\pm \sqrt{g_{\mathbf{k}}^2+h_{\mathbf{k}}^2},$$

$$\displaystyle E_{\text{gap}}=2\sqrt{g_{\bm{k}}^{2}+h_{\bm{k}}^{2}}\Big{|}_{{\bm{k}}=(0,0)}\simeq 1.5\text{ eV}$$

$$f_{\mathbf{k}}$$

$$g_{\mathbf{k}}$$

$$h_{\mathbf{k}}$$

$$t_1$$

$$t_5$$

$$H_{\text{sc}}=\sum _{\mathbf{k}}(\Delta c_{1\mathbf{k}}^{†}{c}_{1,-\mathbf{k}}^{†}+\Delta c_{2\mathbf{k}}^{†}{c}_{2,-\mathbf{k}}^{†}+\text{h.c.}).$$

$$\left(\begin{array}{c}{\eta }_{+}\\ {\eta }_{-}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& \frac{g-\mathrm{i}h}{R}\\ 1& \frac{-g+\mathrm{i}h}{R}\end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\end{array}\right)$$

$$H=\sum _{\mathbf{k}}{ϵ}_{+,\mathbf{k}}{\eta }_{+,\mathbf{k}}^{†}{\eta }_{+,\mathbf{k}}+\sum _{\mathbf{k}}(\Delta {\eta }_{+,\mathbf{k}}^{†}{\eta }_{+,-\mathbf{k}}^{†}+\text{h.c.}).$$

$${\psi }_{\text{left}}={\psi }_{\text{right}},{\hat{v}}_n{\psi }_{\text{left}}={\hat{v}}_n{\psi }_{\text{right}}$$

$$\mathcal{H}=\left(\begin{array}{cc}{ϵ}_{+,\mathbf{k}}& \Delta \\ {\Delta }^{\ast }& -{ϵ}_{+,\mathbf{k}}\end{array}\right)$$

$$\left(\begin{array}{c}{k}_x\\ k_y\end{array}\right)=\left(\begin{array}{cc}\cos \delta & -\sin \delta \\ \sin \delta & \cos \delta \end{array}\right)\left(\begin{array}{c}{k}_n\\ k_{\perp }\end{array}\right).$$

$${ϵ}_{\mathbf{k}}\equiv {ϵ}_{+,\mathbf{k}}=f_{\mathbf{k}}+\sqrt{g_{\mathbf{k}}^2+h_{\mathbf{k}}^2}$$

$$f_{\mathbf{k}}$$

$$g_{\mathbf{k}}$$

$$h_{\mathbf{k}}$$

$${\psi }_L$$

$${\psi }_N$$

$${\psi }_R$$

$$\displaystyle F(k)\equiv\frac{\partial\epsilon_{\bm{k}}}{\partial k_{n}}\Bigg{|}_{k_{n}=k},$$

$$L_e+L_h{\mathrm{e}}^{\mathrm{i}\zeta }$$

$$L_e{\mathrm{e}}^{\mathrm{i}\zeta }+L_h$$

$$R_e{\mathrm{e}}^{\mathrm{i}{k}_{n}L}+R_h{\mathrm{e}}^{\mathrm{i}\zeta +\mathrm{i}\varphi }{\mathrm{e}}^{-\mathrm{i}{k}_{n}L}$$

$$R_e{\mathrm{e}}^{\mathrm{i}\zeta -\mathrm{i}\varphi }{\mathrm{e}}^{\mathrm{i}{k}_{n}L}+R_h{\mathrm{e}}^{-\mathrm{i}{k}_{n}L}$$

$$F(-k_n)L_e+F(k_n)L_h{\mathrm{e}}^{\mathrm{i}\zeta }$$

$$F(-k_n)L_e{\mathrm{e}}^{\mathrm{i}\zeta }+F(k_n)L_h$$

$$F(k_n)R_e{\mathrm{e}}^{\mathrm{i}{k}_{n}L}$$

$$=F(k_n^{\prime })$$

$$F(k_n)R_e{\mathrm{e}}^{\mathrm{i}{k}_{n}L}{\mathrm{e}}^{\mathrm{i}\zeta -\mathrm{i}\varphi }$$

$$=F(k_n^{\prime })$$

$$I=\frac{2e}{\hslash }\sum _{k_{\perp }}\sum _{\pm }\frac{\partial E_{\pm }}{\partial \varphi }f(E_{\pm }),$$

Full text

I Introduction

There is currently much research addressing the physics of two-dimensional materials consisting of a single or very few atomic layers. Besides their interest from a fundamental viewpoint, such materials typically also feature unusual electronic properties which has spurred efforts to identify possible technological applications. Notable two-dimensional systems where the electrons display Dirac physics despite moving at non-relativistic velocities include graphene castroneto_rmp_09 , silicene kara_ssr_12 , and transition metal dichalcogenides duan_csr_15 . Phosphorene is a two-dimensional material which may be regarded as a single layer of black phosphorous, analogously to how graphene is a single layer of graphite. It was recently isolated by mechanical exfoliation li_natnano_14 ; koenig_apl_14 ; liu_acsnano_14 and has attracted attention due to the existence of an intrinsic band gap, in contrast to what occurs in graphene. This thickness-dependent band gap substantially increases the potential for using phosphorene in semiconductor-based technology such as transistors and solar cells. In addition to its promising electronic properties, phosphorene also displays high mechanical flexibility.

The proximity effect in two-dimensional materials was studied early on in graphene and various features of superconducting graphene have been clarified beenakker_prl_06 ; sengupta_prl_06 ; ludwig_prb_07 ; linder_prl_07 (more recently also in silicene linder_prb_14 ; li_prb_16 ; li_prb_16b ; zhou_prb_16 ; li_prb_16c ; kuzmanovski_prb_16 ; paul_prb_16 ). It has been shown that Andreev reflection in superconducting graphene junction is specularbeenakker_prl_06 . Josephson current in ballistic graphene at the Dirac point is formally identical to that in a disordered normal metalTitov . The tunneling conductancesengupta_prl_06 ; linder_prl_07 and Josephson currentMaiti in graphene junctions are oscillatory functions of a width and height of the barrier at the interface.

In this paper, we study how superconducting correlations are manifested in phosphorene. In practice, this is accomplished via the proximity effect to a host superconducting material where a tunnel coupling allows Cooper pairs to penetrate a finite distance into the phosphorene sheet. We address two of the most fundamental superconducting transport phenomena, the Josephson effect josephson_pl_62 and Andreev reflection andreev_jetp_64 ; btk , in a phosphorene sheet. In contrast to superconducting transport in e.g. graphene and silicene, the supercurrent and conductance spectra for phosphorene are strongly anisotropic depending on how the superconducting electrodes are placed on the phosphorene sheet, differing with an order of magnitude depending on the orientation. Moreover, the supercurrent and conductance display oscillations as a function of the distance between the superconducting/normal electrodes with a period that also depends on the orientation of the electrodes. These results highlight that there may exists interesting opportunities with regard to tailoring anisotropic superconducting transport due to the geometry of the setup by using phosphorene.

II Theory

We here describe the derivation of the Andreev bound state (ABS) energies and the resulting supercurrent transport. The starting point is the two-band model of phosphorene ezawa_njp_14 ; sarkar_arxiv_16

[TABLE]

where we use a basis vector of operators

[TABLE]

and the subscript $i$ on the fermionic $c_{{\bm{k}},i}^{\dagger}$ creation operators denote the two atoms in the reduced unit cell. The quantities $f,g,h$ are defined as ezawa_njp_14 ; sarkar_arxiv_16 :

[TABLE]

We have not shown the lattice constants $a_{x}$ and $a_{y}$ explicitly for brevity of notation. The normal-state eigenvalues are

[TABLE]

where $\epsilon_{+,{\bm{k}}}$ has a minimum at the $\Gamma$ point ${\bm{k}}=(0,0)$. The band gap is

[TABLE]

which is much larger than the superconducting gap. Therefore, specular Andreev reflection is disregarded and we focus on the conduction band $\epsilon_{+,{\bm{k}}}$. Performing a low-energy expansion around the $\Gamma$ point yields:

[TABLE]

Above, $t_{4}<0$ and all the other parameters are defined via the hopping terms $t_{i}$ between various atoms ezawa_njp_14 ; sarkar_arxiv_16 (all in units of eV):

[TABLE]

Inserting these low-energy expansions into the expression for the normal-state eigenvalues, one obtains that the spectrum is parabolic in $k_{y}$ when $k_{x}=0$ whereas it is close to linear (but not strictly linear) in $k_{x}$ at $k_{y}=0$ since $\gamma\gg|\alpha|$ and $\gamma\gg|t_{4}|$. To describe proximity-induced superconductivity, we now add a standard superconducting term:

[TABLE]

Any coupling between the conduction and valence bands due to Andreev reflection is irrelevant due to the large band gap $E_{\text{gap}}$, and we therefore simply project out the fermion operators belonging to the $\epsilon_{-,{\bm{k}}}$ valence band. This can be done after rewriting the original fermion operators $c_{i,{\bm{k}}}$ in terms of the band-basis operators which we denote $\eta_{\pm,{\bm{k}}}$. We find (dropping momentarily the ${\bm{k}}$ index for brevity)

[TABLE]

where $R\equiv\sqrt{g^{2}+h^{2}}$. Inserting this into the superconducting Hamiltonian, using that $f$ and $g$ are symmetric in ${\bm{k}}$ whereas $h$ is antisymmetric, and discarding all final terms that contain $\eta_{-}$, we end up with:

[TABLE]

In other words, the effective model is comprised of the Bardeen-Cooper-Schrieffer Hamiltonian, but with a new, anisotropic normal-state dispersion $\epsilon_{+,{\bm{k}}}$. In the band basis, the order parameter remains of the conventional $s$-wave type.

We will consider the formation of ABS in an superconductor/normal/superconductor (SNS) junction which does not necessarily extend along the $x$-axis. In this way, we will be able to probe the effect of the anisotropic band-structure of phosphorene on Andreev reflection and the supercurrent. We use an extended version of the Blonder-Tinkham-Klapwijk formalism btk adapted to materials with a dispersion relation deviating from that of free electrons, using a similar procedure and notation as in linder_prl_08 ; linder_prb_14 . Let $\hat{n}$ be the interface normal to the SN interfaces. Usually, one considers $\hat{n}=\hat{x}$. Now, we want to consider any $\hat{n}$ between $\hat{x}$ and $\hat{y}$ to probe how the supercurrent changes due to the anisotropic normal-state dispersion which is linear in momentum in one direction and quadratic in momentum in another. The setup is shown in Fig. 1 where the angle $\delta$ defines the orientation of the electrodes, so that $\hat{n}=(\cos\delta,\sin\delta)$. For a given orientation, transverse modes will be fully taken into account. The dispersion of the upper band $\epsilon_{+}$ is shown in Fig. 2.

Realistically, there will be a Fermi vector-mismatch between the S and N regions due to charge-transfer between the SC electrode and the phosphorene layer under it. Thus, our calculations include the possibility for different chemical potentials $\mu$ in the two regions, i.e. $\mu_{S}$ and $\mu_{N}$. We will consider doped phosphorene where the chemical potential lies in the conduction band, allowing for a finite density of states that thus supports the formation of a proximity-induced superconducting state. The resulting boundary conditions at each interface of the structures that we consider, both in the supercurrent setup shown in Fig. 1(a) and the conductance spectroscopy setup shown in (b), are:

[TABLE]

where $\psi_{\text{left/right}}$ refers to the wavefunction on the left/right side of a given interface. The velocity operator is defined as $\hat{v}_{n}=\partial\mathcal{H}/\partial k_{n}$ and

[TABLE]

where the normal-state dispersion may be expressed via the momentum parallell to $\hat{n}$ and perpendicularly to it (instead of $k_{x}$ and $k_{y}$) by using that

[TABLE]

Note that when applying the velocity operator, care should be taken to use the correct value of the momentum perpendicular to the interface which is determined by the chemical potential $\mu_{S}$ or $\mu_{N}$ (the S or N region). Let $c\equiv\cos\delta$ and $s\equiv\sin\delta$. In this case, we get

[TABLE]

where we defined

[TABLE]

To determine the properties of the Josephson effect in an SNS phosphorene junction, we set up the wavefunctions in each region (left S, normal, right S):

[TABLE]

Here, we have allowed for a different $\mu$ in the N and S regions by distinguishing the wavevectors $k_{n}^{\prime}$ and $k_{n}$ in these regions while $n$ and $n_{\perp}$ denote the coordinates parallell to and perpendicular to the interface normal $\hat{n}$. Moreover, $\phi$ is the SC phase difference and $\zeta=\text{acos}(E/\Delta_{0})$ where $E$ is the quasiparticle energy. The wavevectors $\{k_{n},k_{n}^{\prime}\}$ have to be obtained from the dispersion relation numerically for a fixed value of $k_{\perp}$. For a given direction of the electrodes $\hat{n}$, we only consider the contribution from propagating modes in the N region, i.e. only contributions from the $k_{\perp}$ values that give a real $k_{n}$ and $k_{n}^{\prime}$. This is required to be consistent with the diagonalization of the Hamiltonian. Unlike Ref. sarkar_arxiv_16 , we will not make any approximations in the band-structure (such as setting $\alpha=0$) and we thus keep all terms in the dispersion relation.

Defining the quantity

[TABLE]

we can write down the system of equations that determine the ABS energies:

[TABLE]

stem from the continuity of the wavefunction, whereas

[TABLE]

and also

[TABLE]

stem from the continuity of particle flux. By writing this system of equations as $Ax=0$ where $A$ is an $8\times 8$ matrix and $x$ is a vector containing all the scattering coefficients, one determines the ABS energies from the requirement det($A$)=0. The supercurrent is then obtained from:

[TABLE]

where the ABS energies have the form $E_{\pm}=\pm\Delta_{0}\sqrt{R(\phi)}$. Here, $f$ is the Fermi-Dirac distribution function. This expression is derived from the fundamental thermodynamical relation between the free energy $F$ of the junction and the supercurrent $I$, namely $(2e/\hbar)dF/d\phi=I$. It is valid beenakker_prl_91 in the short-junction limit $L/\xi\ll 1$ which we will consider in this paper. Here,$L$ and $\xi$ are the length of the junction and superconducting coherence length, respectively. The expression for $R(\phi)$ is too lengthy to be particularly useful, but may be numerically implemented. Using Eq. (21), we can then study how the supercurrent differs for propagation along the $x$-axis ($\delta=0$) and the $y$-axis ($\delta=\pi/2$), as well as its dependence on the length $L$ of the junction. For the transverse modes, we use an effective width $W=500$ nm for the junction in both cases with a spacing of $\pi/W$ between the transverse modes, and the temperature is set to $\Delta_{0}/k_{B}T=100$ (corresponding to $T\ll T_{c}$). Our choice of width $W\gg L$ is large enough to effectively mimic a semi-infinite junction in the transverse direction, as the results are quantitatively indistinguishable upon further increasing $W$.

To study Andreev reflection and conductance spectra in a phosphorene junction, a NN’S junction is described by the following wavefunctions:

[TABLE]

With the boundary conditions, we obtain the scattering coefficients. The expressions for the coefficients $a$ and $b$ are given in the Appendix.

The normalized conductance at zero temperature is then calculated as

[TABLE]

We have included a normal region N’ separating the N and S regions in order to study how the transport properties are affected by gating on N’, which controls the local chemical potential.

III Results and Discussion

III.1 Josephson effect

We begin by considering the direction dependence for a fixed length $L=15$ nm, corresponding to $L/W=0.03$. Due to the charge-transfer between the host superconductor and the region of the phosphorene sheet directly underneath it, we set $\mu_{S}>\mu_{N}$ in order to use experimentally relevant values. Specifically, we set $\mu_{S}=4t_{4}+1.5m$. As a reference value, we note that $\mu=4t_{4}+m$ corresponds to the bottom of the conduction band. Fig. 3 shows the supercurrent vs. phase. It is clear that the supercurrent is highly anisotropic and differs by an order of magnitude when comparing propagation in the $x$- and $y$-directions ($\delta=0$ and $\delta=\pi/2$, respectively) for e.g. $\mu_{N}=4t_{4}+1.05m$ (meaning that the N region is assumed to be only slightly doped). To understand this, one should note that not only the effective dispersion relation is different along these directions, but the number of allowed propagating modes also differs greatly due to the anisotropic Fermi surface. As a result, the supercurrent magnitude becomes larger in the $x$-direction along which the dispersion is closer to being linear, hosting more propagating modes characterized by $k_{\perp}$. The maximum of the supercurrent is shifted slightly away from $\phi=\pi/2$ due to the presence of higher harmonics which are typically present in high-transparency ballistic Josephson junctions.

Next, we consider the length and chemical potential dependence of the supercurrent. The supercurrent shows oscillation due to Klein tunneling pertaining to Dirac fermions.Maiti Besides the difference in magnitude, a qualitative difference emerges between the two directions. In the $\delta=\pi/2$ direction, the oscillation period is much smaller than in the $\delta=0$ direction. This is physically reasonable upon considering the different dispersion relations which affects the wavevector magnitude and thus the oscillation period, since the supercurrent depends on the product $k_{n}L$. For a fixed value of the Fermi level $\mu_{N}$, the Fermi wavevector is much smaller in the $\delta=0$ case for normal incidence $k_{y}=0$ in Fig. 2 than the Fermi wavevector in the $\delta=\pi/2$ case for normal incidence $k_{x}=0$. As a result, the oscillations occur on a shorter length-scale in the $\delta=\pi/2$ case as shown in Fig. 4(b).

III.2 Andreev reflection

We next consider how Andreev reflection is manifested in phosphorene. Experimentally, gating is a commonly used way to influence the transport properties. For this reason, we take into account a phosphorene region covered by a gate electrode that separates the normal and superconducting parts. In effect, we are considering an NN’S phosphorene junction where the chemical potentials are denoted by $\mu_{N}$, $\mu_{N^{\prime}}$, and $\mu_{S}$. In the superconducting region, we fix $\mu_{S}=4t_{4}+1.5m$, and we proceed to determine how the conductance spectra of the system under consideration depends on applied bias voltage $eV$ and the chemical potentials $\mu_{N}$ and $\mu_{N^{\prime}}$. We fix the length $L$ of the N’ layer to $L=$10 nm.

In Fig. 5, we show the conductance for $\mu_{N}=4t_{4}+1.01m$ with (a) $\delta=0$ and (b) $\delta=\pi/2$ and several values of $\mu_{N^{\prime}}$. A clear anisotropy is seen for the two directions $\delta=0$ and $\delta=\pi/2$. Due to a large Fermi wave-vector mismatch, Andreev reflection is suppressed and the conductance shows a gap-like structure. Figure 6 depicts the conductance for $\mu_{N}=4t_{4}+1.1m$ with (a) $\delta=0$ and (b) $\delta=\pi/2$ and several values of $\mu_{N^{\prime}}$. The conductance obtained for the two directions $\delta=0$ and $\delta=\pi/2$ remains noticeably different. Since the Fermi wavevector mismatch is small, the gap-like structures become shallower. In Fig. 7, we show the conductance at zero bias as a function of $\mu_{N^{\prime}}$ for (a) $\mu_{N}=4t_{4}+1.01m$ and (b) $\mu_{N}=4t_{4}+1.1m$. The conductance oscillates as a function of $\mu_{N^{\prime}}$ with different periods for two directions. These oscillations stem from wavenature of the wavefunction for Dirac fermions in the N’ region. As explained in the previous subsection, the oscillation period is smaller in the $\delta=\pi/2$ case due to the difference in Fermi wavevector magnitude when comparing the two directions of propagation. Moreover, the result indicates that by gating, one can tune the anisotropy of the Andreev reflection probability and conductance. Note that the non-normalized conductances for $\delta=0$ and $\pi/2$ differ in magnitude with about one order due to the different number of the transverse modes $k_{\perp}$ contributing to the transport.

IV Conclusion

In conclusion, we have presented a study of the anisotropic superconducting transport properties of phosphorene, a single layer of black phosphorous. The system setup consisting of superconducting or normal electrodes deposited at different locations of a phosphorene sheet should be experimentally feasible in light of the recent experimental reports of stable, isolated phosphorene via exfoliation. Due to the anisotropic band structure of this system, the supercurrent magnitude changes with an order of magnitude when comparing tunneling along two perpendicular directions in the monolayer. The oscillatory behavior of the supercurrent as a function of the length and chemical potential of the junction is different when modifying the orientation of the superconducting electrodes deposited on the phosphorene sheet. For Andreev reflection, we show that gate voltaging controls the probability of this process and that the anisotropic behavior found in the supercurrent case is also present for conductance spectra. The oscillatory behaviors of the supercurrent and conductance found here are manifestation of Dirac dispersions, as seen in graphene junctions.sengupta_prl_06 ; Maiti

Interesting future directions to explore include non-local transport in multiterminal geometries and in particular the crossed Andreev reflection process, as well as the inclusion of magnetic elements in the system setup.

Acknowledgements.

J.L. acknowledges funding via the Outstanding Academic Fellows program at NTNU, the NV-Faculty, and the Research Council of Norway Grant numbers 216700 and 240806. This work was supported by a Grant-in-Aid for Scientific Research on Innovative Areas ”Topological Materials Science” (KAKENHI Grant No. JP16H00988) from JSPS of Japan.

Appendix A Scattering coefficients

The coefficients $a$ and $b$ are given by

[TABLE]

where

[TABLE]

and $\alpha={e^{i{k_{n}}^{\prime\prime}L}}$.