# Anisotropic transport of normal metal-barrier-normal metal junctions in monolayer phosphorene

###### Abstract

We study transport properties of a phosphorene monolayer in the presence of single and multiple potential barriers of height and width , using both continuum and microscopic lattice models, and show that the nature of electron transport along its armchair edge ( direction) is qualitatively different from its counterpart in both conventional two-dimensional electron gas with Schrödinger-like quasiparticles and graphene or surfaces of topological insulators hosting massless Dirac quasiparticles. We show that the transport, mediated by massive Dirac electrons, allows one to achieve collimated quasiparticle motion along and thus makes monolayer phosphorene an ideal experimental platform for studying Klein paradox. We study the dependence of the tunneling conductance as a function of and , and demonstrate that for a given applied voltage its behavior changes from oscillatory to decaying function of for a range of with finite non-zero upper and lower bounds, and provide analytical expression for these bounds within which decays with . We contrast such behavior of with that of massless Dirac electrons in graphene and also with that along the zigzag edge ( direction) in phosphorene where the quasiparticles obey an effective Schrödinger equation at low energy. We also study transport through multiple barriers along and demonstrate that these properties hold for transport through multiple barrier as well. Finally, we suggest concrete experiments which may verify our theoretical predictions.

## I Introduction

Two dimensional crystals composed of single or few atomic layers are the focus of intense research currently, on account of their remarkable optical, electronic and mechanical properties. Moreover such materials, as shown in recent years revgr1 (); revti1 (), serve as test bed for Dirac physics. This property of these materials arises from the fact that their low-energy quasiparticles obey an effective Dirac-like equation. Such quasiparticles lead to a host of unconventional thermodynamic and transport properties. Examples of such unconventional properties seen in the context of graphene and topological insulators include unconventional quantum Hall effect hallref1 (), unusual Kondo effect kondoref1 (), and unconventional transport properties transref1 (); transref2 (). Indeed, the latter property serves as one of the key aspects of Dirac (relativistic) materials which distinguishes them from conventional materials whose quasiparticles obey Schrödinger equation. For example, the former class of materials display oscillatory behavior of the transmission through a potential barrier as a function of the barrier height or width; this behavior is in complete contrast to an exponentially decaying transmission function found for the latter class History_klein (); Klein_nature (); Robinson ().

Phosphorene ph-1 (); ph-review () is another such Dirac material which is being actively investigated for electronic and other applications ph-tr1 (); ph-tr2 (); ph-tr3 () on account of its anisotropic electronic ph-bs1-dft (), thermal ph-thermal1 (); ph-thermal2 () and optical properties ph-bs1-dft (). A hallmark of monolayer phosphorene is its anisotropic band-structure, which displays an almost flat parabolic Schrödinger like dispersion along the zigzag edge ( direction) and a predominantly Dirac like dispersion along the armchair edge ( direction)ph-bs1-dft (); ph-bs2-tb (); ph-bs2-tb-2 (); Doh (); ph-bs3-kp (); ph-bs3-tb (); ph-bs4-tb (). Such anisotropic bandstructure distinguishes phosphorene from other Dirac materials such as graphene or topological insulator surface where the effective Dirac theory is massless and isotropic. Accordingly, the transport properties across single or multiple potential barriers in phosphorene are expected to be different depending on the orientation of the barrier relative to the plane. The anisotropy of the low energy dispersion in phosphorene thus offers an opportunity to simultaneously probe the relativistic as well as the non-relativistic nature of electrons in a transport experiment based on a normal-barrier-normal (NBN) junctions. Depending on the orientation of the phosphorene monolayer with respect to the measuring electrodes, it is expected to display an oscillatory behavior in the transmission function as a function of the barrier strength as well as the barrier height along the armchair direction, and an exponential decaying transmission function along the zigzag direction. However, a detailed theoretical investigation of such properties has not been yet carried out.

It is the aim of this paper to highlight these anisotropic transport signatures of phosphorene across a NBN junction with different orientations. More specifically, we study the transport properties of quasiparticles in monolayer phosphorene in the presence of single and multiple potential barriers oriented along (subsequently referred to as the direction). The main results that we obtain from such a study are as follows. First, we show that the transport shows several signatures consistent with a gapped Dirac quasiparticle; this is in sharp contrast to the case when the barrier is along ( direction) for which the transport properties conform to those due to Schröndinger quasiparticle. Second, we show that for a barrier along the direction, the dominant contribution to the transport comes from the quasiparticles which impinge on the barrier at near-normal incidence provided the applied voltage is close to the bottom of the conduction band, leading to collimated transport of electrons. The degree of such collimation can be tuned by the external applied voltage . Third, we find that this property, which is experimentally inaccessible in gapless Dirac materials such as graphene and topological insulators, allows us to tune to a regime where the conductance mimics the behavior of normal transmission amplitude. This leads to the possibility of observing signature of Klein paradox through measurement of . Fourth, we find that both the normal transmission and displays oscillatory or decaying behavior as a function of the barrier width depending on the relative strength of the dimensionless barrier potential and the applied voltage , where is the mass gap; the behavior of and changes from oscillatory to a monotonically decreasing function of for . We also show analytically that and for using a continuum approximation to the lattice model of phosphorene and demonstrate that the phenomenon described above persists beyond the continuum approximation used to obtain the expression of and . We point out that such transport behavior is qualitatively different from their counterparts in both conventional Schrödinger and gapless Dirac materials. Fifth, we study transport of phosphorene quasiparticles through multiple barriers (each of height and width ) along and demonstrate that all of the above-mentioned features derived for transport through a single barrier holds in the multiple barrier case. The main difference between the two manifests itself in additional peaks in normal transmission or conductance as a function of in the oscillatory regime () for the multiple barrier case; we provide an analytical explanation for this phenomenon for barriers. Finally, we discuss experiments which can test our theory.

The manuscript is organized as follows. In Sec. II, we outline the band structure of phosphorene and chart out the continuum Hamiltonian used for transport calculation. Such calculations for the single potential barrier along is discussed in Sec. III, while transport through multiple barriers along is discussed in Sec. IV. Finally, we discuss our results, chart out possible experiments which can test them, and conclude in Sec. V.

## Ii Low energy effective Hamiltonian of Phosphorene

The band-structure of monolayer phosphorene, is well known from ab-initio calculations ph-bs1-dft (); ph-bs2-tb (), and it has been used to construct effective low energy Hamiltonian using several approaches such as method ph-bs3-kp (), the tight-binding approach ph-bs2-tb (); ph-bs2-tb-2 (); ph-bs3-tb (); ph-bs4-tb (); Doh (), and the methods of invariants ph-bs-inv1 (). All of these methods yield qualitatively similar band structure. Thus for this paper we use the two band Hamiltonian which is obtained from the four band tight-binding Hamiltonian on a discrete lattice, making use of the symmetry ph-bs3-tb () and expanding around the point. The origin of such symmetry can be understood in the following manner. The unit cell of phosphorene, shown in Fig. 1(a), contains four phosphorus atoms, such that the upper and the lower layers each contain two of these atoms. This leads to the point group invariance which constitutes invariance under a shift along the plane of the layer combined with exchange of the layer indices. It is therefore sufficient to consider a reduced unit cell of two atoms either in the upper or in the lower layer since the transfer energy of atoms in these layers are identical.

Our starting point is the effective two-band tight-binding model based on Refs. [ph-bs2-tb, ; ph-bs3-tb, ]. This tight-binding model is given by

(1) |

where the summation runs over the lattice sites, is the hopping matrix element between and sites, and is the annihilation operator of an electrons at site . The corresponding two-band Hamiltonian can be written in momentum space as ph-bs3-tb ()

(2) | |||||

where for are the different transfer matrix elements indicated in Fig. 1(a) with eV, eV, eV, eV and eV, are the Pauli matrices in the band basis, is the two component fermion field with being the annihilation operators corresponding to electrons of the two phosphorene atoms from the unit cell, denotes the identity matrix, and we have scaled wavevectors by the corresponding lattice lengths and . The energy spectrum of this Hamiltonian is given by

(3) |

where the sign in the subscript indicates conduction (valence) band. A plot of these bands as a function of (for ) and (for ) is shown in Fig. 1(b). The corresponding eigenvectors are given by

(4) |

where and denotes the signum function.

The Hamiltonian may be simplified in the low-energy low-momentum or continuum limit. Using the identities and , one can obtain the continuum version of Eq. (II): , to be

(5) |

where is the chemical potential, and , with eV. The other parameters in Eq. (5), can be expressed in terms of in a straightforward manner and are givenph-bs3-tb () by eV, eV, eV, and eV. Here we note that which allows one to have an approximately linear dispersion in for small . The eigenvalues and eigenfunctions of are given by

(6) | |||||

and . From Fig. 1(b), we find that the energy bands of the continuum model match with those of the lattice model. We note that in the continuum limit, since , the Hamiltonian is linear (Dirac like) in but parabolic (Schrödinger like) in ; consequently, the nature of the quasiparticle transport and the effect of a potential barrier depends crucially on the its orientation in the plane. We shall discuss this phenomenon in detail in the next section.

## Iii Transport across a single barrier: NBN junction

In this section, we present an analysis of transport across a single barrier. We first present the general formalism for both the lattice model (Eq. (II)) and its continuum approximation (Eq. (5)) in Sec. III.1. The numerical results of the analysis of the formalism developed in Sec. III.1 is presented in Sec. III.2.

### iii.1 Formalism

In this section, we study ballistic transport in phosphorene monolayer in the presence of a single barrier of strength and width as shown in Fig. 1(c). We will primarily be focussed on transport in the direction, which is likely to have signatures of Dirac quasiparticles.

Let us first discuss the formalism for calculating transmission and conductance along using the two band lattice Hamiltonian of Eq. (II). The wave function for the particles with transverse momenta and energy (which we choose to lie in the conduction band), where is the chemical potential and is the applied bias voltage, moving along direction can be obtained by diagonalization of Eq. (II). The corresponding eigenfunctions can be read off from Eq. (4); for an electron moving along (right (+) and left (-) moving electrons) with momenta , it is given by with

(7) |

where and is obtained from the solution of

(8) |

We first note that the sign here in the superscript indicates right(left) moving electrons in contrast to those in the subscript in Sec. II which indicated conduction and valence band energies. Since for transport we shall focus only on electrons in the conduction band in region I, we omit the subscript indicating conduction or valence band for brevity. We also note that Eq. (8) admits two solution for corresponding to any given energy and transverse momentum . One of these solutions is real and represents a propagating wave while the other is imaginary and represents an evanescent wave. While the expectation of the current operator and hence the conductance receives contribution from the propagating wave, the formalism that we use requires that the evanescent wave solutions are treated at an equal footing. In what follows, we denote the solutions, obtained by numerical solution of Eq. (8), as and with . In terms of these wavevectors, one can write the wavefunction in region I (which extends in the region as seen in Fig. 1(c)) as

(9) | |||||

where , , and is obtained using Eq. (4). In Eq. (9) we have kept only the decaying evanescent modes, and and denotes reflection coefficients corresponding to the propagating and the evanescent modes respectively.

The wavefunction in region II (as shown in Fig. 1(c)) consists of left and right propagating modes as well as evanescent waves which decays and grows within . Here depending on the magnitude of the applied barrier potential and , the longitudinal wavevector may be obtained from either the valence or the conduction band energy expressions (Eq. 3). This is in contrast to the situation in regions I and III where they are necessarily obtained using energy expressions for the conduction band. The longitudinal wavevector for these modes are obtained from the solution of

(10) |

and we denote these wavevectors as and . In terms of these, one can write the wavefunction in region II as

where , , , and . In Eq. (LABEL:wav2), and are coefficients corresponding to the left and the right propagating mode, and and are the coefficients corresponding to decaying and growing evanescent mode. Note that since region II extends within a finite span of , the solution with the growing evanescent mode is admissible in this region.

Finally, we consider the wavefunction in region III which extends for . In this region, one finds

where , and are the transmission amplitude corresponding to the propagating and the evanescent mode. Note that in this region, similar to region I, only the decaying evanescent mode is admissible.

To obtain the reflection and transmission coefficients, we use the standard procedure of imposing the continuity condition for the wavefunction as well as the current at and . This leads to the conditions

(13) |

where the velocity operator can be obtained using Eq. (II) and is given by

(14) | |||||

with the substitution . Substituting Eqs. (9)-(LABEL:wav3) in Eq. (13), we can then obtain the following eight equations,

(15) | |||

where and the quantities and are defined as

(16) |

Here can assume values depending on values of : or , takes values of or depending on or , and we have suppressed the and energy dependence of and for clarity. Note that the presence of the evanescent waves (complex solutions for ) is essential for unique solution of Eq. (15); without them, we would have an overdetermined set of equations.

Equation (15) is solved numerically to obtain the reflection and the transmission coefficients and . The transmission probability can be computed in terms of these as . Note that which corresponds to the transmission amplitude of the evanescent mode decays exponentially in region III away from the barrier and hence it does not contribute to the transmission for , where is the system size and the lead is placed at . To compute the conductance in the presence of a single barrier we first note that the current operator along is given by and thus the current flowing in regions I and III are given by

(17) |

The conductance of the system can then be computed as

(18) |

where , and denotes the maximum transverse momenta for which Eq. (6) admits a real solution for and the limits of integration are decided by the fact that the dispersion in Eq. (3) is symmetric in .

Having described the formalism to calculate the transmission and conductance for the two band lattice Hamiltonian, we now focus our attention on the low energy and small momentum limit, for which the system is described by . In what follows we shall neglect the term in the expression of ; this is justified by the fact that for low momenta and so thats the dependence of may be neglected. The motivation for solving the problem within this approximation is two-fold. First, as we shall see, the conductance and transmission coefficients computed within this approximation matches those from the exact lattice model at low applied voltages and second, this approximation yields an analytic expression for the transmission coefficient which allow further insight into the transport properties of the system.

For electrons described by Eq. (5) with , let us consider a wave incident on the barrier with energy and transverse wavevector . The wavefunction of the electron can be computed using Eqs. (5)-(6) as , where , and is given by

(19) |

The reflected wavefunction in this case is given by . Thus the wave function in region I (see Fig. 1(c)) can be written as

(20) |

We note that the continuum approximation does not support the evanescent modes found in the lattice formulation. This feature is consistent with the fact that for the effective Dirac-like Hamiltonian with linear dispersion that we found within this approximation, we only need continuity of the wavefunctions at the barrier edges ( and ); the continuity of the derivative of the wavefunction is no longer necessary since we are dealing with linear differential operators along . This situation is in contrast to Schrödinger electrons for which both the wavefunction and its derivative needs to be continuous across the barrier transref1 (); transref2 ().

In the barrier region, the wavefunction can be expressed as

(21) | |||||

where , and is given by

(22) |

Finally, in region III, the transmitted electron wavefunction is given by

(23) |

The reflection and the transmission coefficients can be determined using standard procedure of demanding wavefunction continuity at the interfaces at and : and . This yields

(24) | |||

A solution of Eq. (III.1) leads to the transmission and reflection amplitudes within the continuum approximation,

(25) |

where and . The transmission probability is then given by

(26) |

The conductance may then be obtained using

(27) |

Note that Eq. (27) is a simplified version of Eq. (18), for the case when the incoming and the outgoing regions across the barrier are identical. We shall analyze the results obtained from the formalism developed in this section for both the lattice model and its continuum approximation in Sec. III.2.

Before ending this section, we analyze Eq. (26) in the thin barrier limit where and with a fixed ratio , and . In this limit, Eq. (26) may be further simplified to obtain

(28) |

Equation (28) can be analyzed to obtain several characteristics of the transmission coefficient . We first note that as (which means ) and , which reproduces Klein paradox result for massless Dirac fermions seen in graphene. Second for , , and consequently as a function of the effective barrier strength oscillates with a frequency as or is varied; thus knowing and , one can estimate the mass of single-layer phosphorene by measuring the frequency of such oscillation. Third, for as is tuned such that , where

(29) |

becomes imaginary and hence changes from being oscillatory to a decaying function of . For , we find that

(30) |

Such a qualitative change in , which is generally not seen in other massless Dirac materials, may serve as an accurate measurement of the mass gap in phosphorene. It is to be noted that while indicates that the applied voltage is larger than the barrier height which is expected to lead to non-decaying behavior of (and hence ) with , the oscillatory behavior of for is a property of Dirac nature of the phosphorene electrons. The difference of such transport from that in graphene follows from the decaying behavior of for ; such a behavior is absent for transport in graphene and topological insulators where and are oscillatory functions of for any . We find that this property is not a consequence of the thin barrier limit can also be seen from Eqs. (26) and (22). Clearly the decaying behavior arises only when the wavevector in the potential region, is imaginary. Fourth, we note that the contribution to the transmission modes comes from transverse momenta modes for which away from the barrier. This requires the condition , where we have set the chemical potential so that corresponds to the bottom of the conduction band. Thus by tuning the applied voltage one may reach a regime where the conductance

(31) |

receives its contribution only from the quasi particles with near normal incidence. This leads to collimated transport; furthermore since the dependence of in this regime is negligible, the oscillation frequency of with the barrier width would mimic that of . This enables us to realize a setup where Klein paradox could be realized via measurement of ; we note that it is impossible to reach this regime in gapless Dirac systems such as graphene. We shall study the feasibility of this proposition in details in the next section. Finally, we note that in contrast to graphene, the presence of a finite mass gap allows us to tune transmission through the barrier. This can be seen by inspecting Eq. (22); we find that there are no propagating modes inside the barrier, for . The key point is that the value of can tuned to zero by choosing which allows one to tune, particularly for large , transmission through the barrier by tuning .

### iii.2 Results

In this subsection, we shall chart out the results corresponding to the theory of transport of monolayer phosphorene electrons through a single barrier along the direction developed in Sec. III.1. We shall first analyze the results for the continuum model (Eq. (II)) and then compare its prediction with those obtained from the lattice model (Eq. (5)).

To this end, we first plot the normal transmission as a function of the barrier width for several representative values of (barrier height) for fixed (incoming energy) as shown in Fig. 2. We note that for (or alternately ), displays oscillatory behavior with increasing which is similar to that of massless Dirac-like electrons as seen in graphene Klein_nature () and in contrast to that of conventional Schrödinger electrons with parabolic dispersion for which the transmission probability decreases monotonically with . However, in contrast to graphene, shows a decaying behavior as a function of for [or alternately for ; (Eq. (29))]. As is further decreased below (or ), becomes an oscillatory function of with frequency given by (green dashed plot in the left panel of Fig. 2). A similar behavior is seen for (Eq. 26) as shown in the right panel of Fig. 2. We note that this leads to tunability of transmission which is a result of both the band structure of phosphorene and the presence of the mass gap ; such behavior is therefore absent in graphene.

Next, we plot as a function of for a fixed , , and for several representative values of as shown in the left panel of Fig. 3. As we discussed in the last section, the transport becomes increasingly collimated as the applied voltage is tuned towards the bottom of the conduction band (which in our notation corresponds to ) since lesser number of modes satisfies Eq. (19) for real values of . This feature, which holds beyond thin barrier limit as shown from the plot of (Eq. 26) in the right panel of Fig. 3, can be shown to be related to the relative flatness of phosphorene bands near the band bottom along the -direction in contrast to that in the -direction; this naturally leads to collimated behavior. Thus, at sufficiently low , dominates the conductance as shown in the left panel of Fig. 4. At low , the oscillatory behavior of as a function of , which is a signature of Klein paradox for gapped Dirac systems Klein_nature (), may therefore be observed via measurement of tunneling conductance . We note that such a measurement would be impossible in graphene since one needs to be very close to the Dirac point to observe this phenomenon where the density of state is extremely small. Such an equivalence between and (Eq. 26) is lost for larger as shown in the right panel of Fig. 4.