Atomistic deconstruction of current flow in graphene based hetero-junctions
Abstract
We describe the numerical modeling of current flow in graphene heterojunctions, within the Keldysh Landauer Non-equilibrium Green’s function (NEGF) formalism. By implementing a -space approach along the transverse modes, coupled with partial matrix inversion using the Recursive Green’s function Algorithm (RGFA), we can simulate on an atomistic scale current flow across devices approaching experimental dimensions. We use the numerical platform to deconstruct current flow in graphene, compare with experimental results on conductance, conductivity and quantum Hall, and deconstruct the physics of electron ‘optics’ and pseudospintronics in graphene junctions. We also demonstrate how to impose exact open boundary conditions along the edges to minimize spurious edge reflections.
Keywords:
NEGF, RGFA, graphene, electron ’optics’pacs:
85.80.Fi 81.07.Nb 73.63.Rt 72.80.Vp∎
1 Introduction
Ever since its inception (1), graphene as a material has occupied a unique place between one dimensional pi-conjugated organic molecules and three dimensional bulk crystalline solids with well defined bandstructures. The quasi-ballistic scattering lengths (2); (3), photon like dispersion (4) and chiral electron flow (5) together promote nontrivial transport physics in graphene, opening up various device possibilities that necessitate detailed numerical modeling. While the photonic bandstructures argue for a continuum model of modest sized device segments, scattering at heterojunction interfaces (6) and edges requires a detailed atomistic treatment. The aim of this paper is to describe practical techniques for simulating quantum transport through graphene based heterojunctions, compare various transport regimes with experiments associated with the chiral nature of electron transmission at graphene heterojunctions. We also point out some of the unique electronic properties of bilayer graphene, but unless otherwise stated, the results are for single layer graphene.
2 Simulation platform
The Non-equilibrium Green’s function (NEGF) formalism (7) provides a unified, ‘bottom-up’ platform for modeling quantum flow of electrons. Indeed, it has been successfully applied to understand transport physics in materials and systems as diverse as organic molecules (8), carbon nanotubes (9); (10), graphene (11); (12); (13); (14); (15); (16), silicon nanowires (17); (18); (19); (20); (21); (22), spintronics, nanomagnets (23); (24), nanoscale phonon transport (25); (26); (27). A challenge however is the sheer problem size associated with atomistic deconstruction of experimentally relevant dimensions, typically hundreds of nanometers, as well as the scattering physics at the atomistically rough edges. We employ a recursive technique to simplify the problem size in order to make such a deconstruction tractable.
2.1 Recursive Green’s function Algorithm (RGFA)
The central entity for quantum kinetics in a weakly interacting system is the retarded Green’s function, defined as,
(1) |
is the Hamiltonian matrix of graphene, described in this paper using a minimal one orbital basis per carbon atom with eV being the hopping parameter. are the self energy matrices from the source and drain contacts, and are the corresponding anti-Hermitian parts representing the energy level broadening associated with charge injection and removal. The corresponding spectral function and quantum mechanical total transmission are given by
(2) |
The local density of states is obtained from the diagonal elements of divided by .
The contact self energy matrices are calculated from, , where is the surface Green’s function of the contact and is the coupling matrix between the contact and graphene. The surface Green’s function is found by solving the following equation recursively
(3) |
where for the layered contact unit cell. For speedy convergence of , we use a decimation technique outlined in (28).
Inverting the matrix in Eq. 1 is computationally expensive and goes out of hand even for a small structure (width5nm). The Recursive Green’s Function Algorithm (RGFA) is a fast practical method to avoid inverting the entire matrix by brute force, but instead extract only select, relevant blocks of the matrix, as described in Ref. (29), (30). For instance the transmission across a layered device involves only a corner block of the matrix, , while the density of states involves only a diagonal block , enabling such partial inversion.
Evaluating the Green’s functions piece by piece
We consider the device area in layers of N as shown in Fig. 1. The right connected Green’s function (30) is calculated from
(4) |
where is the Hamiltonian of the th layer, is the coupling matrix between adjacent layers and is the electrostatic potential of the th layer. The calculation of is done by stepping through the above equation from right to left starting at layer , with the right connected Green’s function for the th layer
(5) |
Once we reach layer 1 and extract , the device Green’s function component for that layer is then calculated from
(6) |
where . The remaining layer block diagonal matrices are calculated by stepping through from left to right using,
(7) |
with .
From Green’s function to conductance and non-equilibrium carrier density
From the diagonal blocks of the layer Green’s functions, we can calculate the spectral function
(8) |
The left connected spectral function for layer is given by,
(9) |
which requires evaluating a corner block of
(10) |
once again by stepping through from left to right, starting with in layer L=2. All the matrices from Eqs. 4-10 are , where is the number of atoms in one layer. Since matrix inversion is an procedure, inverting substantially smaller matrices at a time leads to considerable computational savings overall.
The total conductance is calculated from the transmission over all modes
(11) |
where takes care of spin degeneracy, while valley degeneracy is captured with the graphene Hamiltonian
The non-equilibrium carrier density at each layer is calculated from the above quantities as
(12) |
where and are Fermi functions at the source and drain respectively.
Extracting current density
The current from th atom to h atom is calculated from (7),
(13) |
where the electron correlation function, and scattering function, . The source and drain Fermi levels are at and . To see the current distribution in the device, we apply a small drain bias . is nonzero only if the th atom and h atom are neighbors. The total current at one atomic site can be found by adding all the components vectorially, , where is the angle of the individual current vectors with respect to a reference direction. As expected, the NEGF expression for satisfies the steady-state Kirchhoff’s law, . The terminal current is the sum of all currents in one atomic layer, and should be equal to the current from the Landauer formalism using the total transmission,
(14) |
While the total current (Eq. 14) can be calculated efficiently using the partial blocks in the RGFA, the spatially resolved current density (Eq. 13), useful primarily for visualization purposes, employs the full matrix and is thus computationally expensive. The recursive algorithm to calculate blocks of is described in (29). The current density from layer to is given by,
(15) |
The matrix has the size of carrying the atom to atom current components between the two layers.
2.2 space formalism (KSF)
We can further cut down the problem size when the system is periodic along the direction perpendicular to electron transport, by employing Bloch’s theorem. The periodicity allows us to decompose the system so that the transverse k-points act as decoupled 1-D chains or modes, whose contribution to the total current can be calculated independently. We can thereafter calculate the conductance through bulk graphene or graphene nanoribbons (GNRs) with mode by mode resolution. As we shall shortly see, it also allows us to visualize the transition from the smooth transmission profile of bulk graphene to the stepwise transmission of GNRs driven by width quantization.
Consider a graphene sheet with an applied potential varying only in the direction (Fig. 2a). After dividing the sheet into blocks, we can split the hamiltonian into block matrices () describing the interaction between the wave functions and belonging to blocks and , with subscripts denoting the rows labeled along the transverse direction, and superscripts denoting layer columns labeled along the transport direction (Fig. 2b). Using this pictorial representation of the Hamiltonian, the Schrödinger equation for the block is given by
(16) |
To take advantage of the manifest periodicity in , we assume Bloch type solutions , with the vector describing the position of the block . Substituting this solution into Eq. 16, we get
(17) |
with
(18) |
(19) |
(20) |
For each , Eq. 17 with a varying index can be interpreted as a one dimensional chain (Fig. 2c) decoupled from chains with other s. In other words, we have split the graphene sheet into a set of independent chains or modes. Now, we can use the usual NEGF formalism to find the transmission of each chain or mode ().
(21) |
The quantities (, , and ) in real space are found by inverse Fourier transform, since our Bloch assumption is equivalent to Fourier transforming from the discrete -space to -space. For example, the Green’s function is given by
(22) |
with , the number of points in space, and the dependent Green’s function given by
(23) |
The sigma matrices are calculated as before, but the surface Green’s function is now dependent.
(24) |
where . takes care of all the couplings in the transverse direction.
Thus the combination of semi-infinite contacts in the direction and Fourier transformation along the direction amounts to a change in basis set for the completely periodic bulk 2D graphene. Since trace is preserved under the basis transformation, we can calculate transmission directly in space
(25) |
where
(26) |
The distinction between bulk graphene and GNR arises primarily from the different sets of ’s (in practice, one needs to worry about edge reconstruction, next nearest neighbor interactions, strain and roughness at the edges of GNRs. For an ideal GNR of width the hard wall boundary conditions restrict to the discrete set
In the case of bulk graphene we can assume periodic boundary conditions every , which restrict the ’s to
and let . Note that as increases, the transmission evolves from a stepwise curve showing width quantization to a smooth curve.
2.3 Applying open boundary conditions
Since simulations deal with finite sized matrices, it is important to make sure the results do not get influenced by edge effects. While realistic physical devices do have edges, their impact on transport is complicated, dominated by strain and reconstruction, roughness, possible charge trapping, edge dipoles, localized vibrons and other subtleties. In fact, measurements of tilt-dependent junction resistance in graphene show a surprisingly weak influence of edge scattering (16); (31). For modeling purposes therefore, it is important to ensure that additional edge effects do not creep in simply because of the finite-sized domain of our simulation grid.
Section 2.1 deals with hard-wall boundary conditions along the transverse direction. For such structures, eliminating edge effects requires making the widths very long (often longer than device length), which may have added repercussions on material properties (see section 3.1). Section 2.2 implemented periodic boundary conditions (PBC) instead with a k-space formalism, allowing us to work with just a couple of nearest neighbor unit cells along the transverse direction, but requiring strict periodicity along that direction. In this section, we will discuss how to implement open boundary conditions (OBC) that would allow us to irreversibly remove an electron impingent at the edges. The trick is to treat the edges as virtual ‘contacts’ with self-energy matrices built out of complex numbers, acting like energy-dependent lossy potentials.
In this section, we will show how to generate exact open boundary conditions (OBC) at the
edges of a device. It is worth emphasizing at this point that these are quantum boundary
conditions, i.e., acting on the retarded Green’s function . What is much harder to solve is a thermal open boundary condition acting on , because this
depends on details of the scattering processes outside the simulation regime. In metallurgical
contacts, we impose by fiat Fermi-Dirac distributions, but in the extended geometry just
below our simulation regime, the distribution is non-equilibrium and the corresponding
boundary conditions need to be estimated self-consistently with the inside the device in order to prevent the outside regime from sucking out the electrons too
aggressively or else inadequately.
The ‘scooped’ open boundary self-energy
The ‘obvious’ way to implement some kind of open boundary is to attach virtual leads along the transverse direction and extract their self-energies using recursion. This is schematically described by Fig. 3(h). While this method will take the electrons away, the leads do not completely span the outside regime relative to the central device represented by the Hamiltonian [H], and the missing chunks at the corners are expected to create spurious reflections. One could extend the contacts or laterally using space formalism to span those corner regions, but the problem is that the four virtual leads are ultimately decoupled from each other with no bonds running in between, and thus do not quite represent a straightforward extension of the device domain beyond the simulation regime.
What a proper OBC will need is a geometry like Fig. 3(g), where the virtual
‘contact’ should be a single monolithic structure with the
central device region ‘scooped’ out. The single self-energy matrix for this
outside regime must be reverse engineered so that
the Green’s function of the central region matches that of the infinite system periodically
extended along all axes. Once we find this self energy (section 2.3.2), we can use it for
any modified condition (e.g. a gate voltage, a molecular adsorbate or an injecting contact)
that only influences the central region but keeps the outside unaltered.
Calculating self energy
To find the boundary free real-space Green’s function in the central region, we first extract the space version corresponding to the infinite periodic system.
(27) |
where the k-dependent quantities are Fourier transformed sums over nearest neighbors. From this we can calculate the real-space Green’s function projected onto the device region by inverse transforming
(28) |
for the atoms ‘m’ and ‘n’ spanning the central device segment. For a system with atoms, is .
We can now reverse engineer the open boundary self-energy, armed with this real space Green’s function and the device Hamiltonian , giving us
(29) |
This self energy can be used to calculate transmission through the structure with open boundary conditions, additional injecting/removing leads and any other perturbations acting only on the central region through a localized potential that does not influence the boundaries
(30) |
Verifying the open boundary condition
In Fig. 3, we verify the efficiency of the open boundary condition by launching a Gaussian wavevector with a spread and initial quasi-momentum set by the wave-vector ,
(31) |
about an injection point . We calculate the evolution of the gaussian wave-packet with time by numerically solving the time-dependent Schrödinger equation within the Crank Nicholson algorithm (described below), acting on the time-dependent vector with row entries corresponding to the spatial coordinates
In the absence of any open boundaries, the wavepacket stays confined within the simulation domain and bounces around (not shown). With the scooped self-energy (Fig. 3g), the open boundary condition takes the electron out at the boundary (Figs. 3a-e). Note that the self-energy we used works only at a specific energy set by the pair , although we can extend it to all energies by Fourier transforming.
3 NEGF simulation of electron transport in graphene
In the following sections, we will use a combination of numerical techniques outlined so far to simulate current flow through graphene segments comparable in size to experimental dimensions. We show how these simulations accurately capture the nuances of conductivity, conductance, magnetotransport, pseudospintronics and Klein tunneling within a common, unified simulation platform.
3.1 Minimum conductivity vs conductance
One of the unique properties of graphene is that its minimum conductivity (rather than conductance ) is quantized in units of , when the width to length aspect ratio is large. In other words, is proportiontal to the aspect ratio in wide graphene samples. The width-dependence is expected in a large ballistic device as the number of modes is roughly proportional to the number of Fermi half-wavelengths one can fit into . However, the average transmission per mode, and thus the conductance, is usually length-independent. This is because propagating modes have transmissions of unity in a ballistic channel, while evanescent tunneling modes have transmissions that are nearly zero. However since graphene is semi-metallic, a wide sheet supports a nearly continuum set of sub-bands with ultralow band-gaps and tunneling probabilities that are not ignorable. In fact the tunnel transmission terms () from these closely spaced modes all add up to an overall dependence of the total transmission (more precisely, a factor ) and therefore the conductance , leading to the conductivity quantization. Experiments on dirty samples show such a quantized minimum conductivity, albeit off by a factor of , usually attributed to Coulomb scattering from charge puddles (32); (33); (34).
We will use the NEGF technique to extract the conductivity of graphene at the Dirac point. To recap, for large scale graphene devices (), the conductance per mode is expected to be (the mode count is given by the integer part of ), while for smaller lengths () and larger widths, the conductivity is expected to approach . We will simulate the structure shown in Fig. 4(a) to calculate the minimum conductivity of a clean graphene sheet. The graphene device is connected to two heavily doped graphene contacts. As long as the contacts provide a large number of modes and thus minimal series resistance, the details of the contact metallicity should not matter. The energy band diagram for the problem is shown in Fig. 4(b). We model the system atomistically with RGFA as described in section 2.1.
We use Eq. 11 to calculate the total conductance . Minimum conductivity is shown in Fig. 4(c) along with the experimental data from Ref. (33). We take two different widths = 100, 200 nm to calculate and for various lengths (). is found to be very close to . As expected, the minimum conductivity stays close to this number as long as is large and a continuum evaluation of the tunnel contribution is a valid approach. In the opposite limit the conductivity increases with length and we reach the conventional regime of conductance quantization (Fig. 4(d)). Note that for a graphene ribbon of finite width, the valley degeneracy of the fundamental mode is broken (35), so that the minimum conductance is . At higher gate voltages as more modes come into play, the conductance per mode reaches , with the degeneracy of 4 coming from the presence of 2 spins and 2 valleys at the same energy.
3.2 RGFA for finite sized graphene devices
We will now apply the RGFA technique discussed in section II to simulate three types of finite sized graphene based devices - uniformly doped, step junction and barrier junction.
Uniformly gated graphene
We apply RGFA to calculate the density of states (DOS) and conductance of graphene sheets of various widths (). Fig. 5(a) shows the DOS (and conductance in inset) for = 10 and 100nm. Narrower nanoribbons show a gap as well as Van Hove singularities due to quantization for select chiralities while wider sheets show a quasi-linear bulk graphene density of states. In practice, nanoribbons also need attention to edge state dynamics, particularly the presence of strain and roughness (12).
A striking property of graphene is its anomalous integer quantum Hall effect. When the Fermi energy lies between two Landau levels (LL) in presence of a magnetic field, conventional two-dimensional electron gases show a vanishing longitudinal resistance and a quantized Hall conductance where is a non-negative integer representing the number of filled Landau levels. In contrast, chiral quasiparticles in graphene exhibiting Berry phase show (36); (37) a Hall conductance , with LL defined as , where is the cyclotron radius. This non-conventional sequence can be explained with the presence of a LL at E = 0 resulting in a four fold degeneracy at zero carrier density (from spin and valley degeneracy). On the contrary, Bilayer graphene LLs are defined as, , where is the cyclotron frequency. Now both and are at zero energy producing thereby an overall eightfold degeneracy at zero carrier density and QHE plateaus of (38).
Numerical modeling magnetotransport in graphene requires a minor modification to the transport scheme outlined earlier. We replace the kinematic momentum with the quasi-momentum , so that the plane wave terms in the Bloch representation pick up an additional phase of where is the magnetic vector potential such that . Thus the hopping parameters between atoms ‘m’ and ‘n’ are now given by
(33) |
For directed magnetic field, our guage is, . We can now turn on a magnetic field perpendicular to the sheet, modify the hopping terms as above, and extract the transvese (Hall) conductance as a function of varying Fermi energy location. The Hall conductance , where is the current, is the number of filled Landau levels and the Hall voltage the difference between the electrochemical potentials of the and edge states. These states represent skipping orbits along the edges created by the cyclotron orbits in the bulk, and are separately in equilibrium with the left and right contact Fermi energies.
Fig. 6b shows the calculated Hall conductance yielding a series of plateaus for both single layer and bilayer graphene. Fig. 6c shows the local carrier distribution functions obtained from the ratio of the local carrier density and the local density of states, in other words, the ratio
(34) |
where is the index for an atom belonging to either top or bottom edge. The Hall voltage turns out to be (Fig. 6c) giving conductance plateaus for . The electrochemical potentials () don’t change along transport direction so that the longitudinal resistance (given by the drop in electrochemical potential along the channel) is zero.
Single p-n junction
junction heterostructures in graphene act qualitatively different from conventional junctions. In normal semiconductors, the presence of a band-gap blocks electron flow, so that reducing (increasing) the built-in potential across a junction with a drain bias leads to an exponentially increased (reduced) current, creating a rectifying current-voltage (I-V) characteristic. The lack of a bandgap removes such rectification in graphene junctions (GPNJs). However, the chiral nature of the electron flow provides fascinating electron ‘optics’ behavior, reminiscent of Snell’s law, albeit with some notable twists. Across a junction for instance, the conservation of transverse quasi-momentum () leads to anomalous Snell’s law
(35) |
where the refractive index is set by the Fermi wavevector and thus the gate voltage applied. The opposite sign of the voltage across a junction, arising when the quasimomentum component flips sign in going from conduction to valence band, means that the system acts like a negative index metamaterial.
While the electron trajectories are set by the above Snell’s law, i.e., the conservation of transverse quasi-momentum, the actual transmission probabilities are set by the conservation of pseudospins. The pseudospins arise from the eigenvectors corresponding to the photon-like graphene E-k. The eigenfunctions are given by
(36) |
where is the 2-D angle setting the quasi-momentum and denotes the conduction or valence band (An additional sign flip with respect to occurs near the other valley). Up and down pseudospin states are given by and , corresponding to the bonding and antibonding combinations of the dimer basis sets in graphene. Analogous to optics, we can derive the equivalent ’Fresnel equation’ for transmission probabilities, by matching the pseudospinor components across the junction. Assuming a split distance ‘d’ between the backgates defining the junction, and a voltage difference between the gates, the transmission works out to be
(37) |
below a critical angle defined from the Snell’s law (Eq. 35), () representing the incident (refracted) angle for electrons. The first term shows that at normal incidence () the transmission function is unity regardless of the size of the barrier . The absence of normal reflection is a manifestation of Klein tunneling, and arises because backscattering requires the flipping of pseudospins, in other words, a fast spatially varying potential. The second, exponential factor (drops out for unipolar or junctions) in the transmission equation comes from tunneling across the slowly varying junction, where the transverse modes face a momentum-dependent band-gap, reminiscent of cut-off frequencies in a wave-guide.
The overall conductance of a GPNJ can be calculated by summing over all transmissions, , which can be thought of as the total number of modes times an average transmission. Thus the average transmission over all modes at a particular energy is given by , where number of modes, . From RGFA calculation, we get the total conductance in units of .
Figs. 7 (a,b) show the doping-dependent resistance of a GPNJ for a 100nm wide graphene sheet, for abrupt and slowly varying junctions respectively. In fig. 7, we plot the variables against the shifts in the Dirac point in both regions, i.e. and . The resistances at the upper left and lower right corners of the plot are higher than the other two ( vs uniformly doped), resulting in an asymmetric resistance vs doping in fig. 7(d) (plotted for specific doping values as indicated with horizontal lines in 7(b)).
The WKB term in Eq. 37 is only present in the junction regime, and that is why only the junction resistance is affected while going from abrupt to smooth junctions (a-b). For a or junction the Fermi energy does not cross the smoothly varying Dirac point anywhere in the device and the transmission expression only includes the wave-function mismatch term. Fig. 7(c) shows the resistance variation for a fixed built-in potential (along black and red lines in fig. 7(a) and (b)). In the conductance variation (inset), we see effectively two Dirac points, which has a simple physical explanation. Recall that the normalized conductance can be decoupled into the mode count from one end times the average transmission over to the other side
(38) |
The left conductance minimum at is where becomes zero, while the right one at is where the average transmission for all modes, , becomes zero. To make this clear, we calculate numerically (Fig. 7(e)), clearly showing a vanishing transmission at the second Dirac point. To calculate , we first simulate a graphene device with uniform doping and extract the overall mode count from the ballistic conductance (). We then simulate the device with different dopings (finite built-in potential, , and conductance ). The ratio of and at each energy gives us . For example, for a symmetric GPNJ (at ), we have
(39) | |||||
leading to , which is close to what we get in the numerical calculation. Thus the conductance of a symmetric GPNJ is 2/3 of the uniformly doped graphene sheet with the same Fermi energy (). In Fig. 7(e) we also show the analytical calculation of (circles) in the same manner as done above but this time using the transmission in Eq. 37. The analytical calculation matches very closely with the numerical calculation (solid lines).
A more direct evidence of chiral tunneling is seen from the resistance across a tilted junction. We use the following equation to extract the junction resistance (7), that eliminates the quantized resistance from the contacts.
(40) |
We calculate and numerically from NEGF as discribed earlier. Fig. 7f shows the variation of plotted against (along orange line) for a 100nm wide graphene sheet. For low carrier density in the middle, the background doping () makes the device an or junction and is low. For , i.e., in the wings of the plot, the device reaches the junction regime and becomes high. This characteristic doping dependence of junction resistance has been seen in experiments in the past (39); (40), usually captured by the odd resistance, . The chiral tunneling can be seen more directly by varying the tilt angle, as seen in experiments recently (16). The purple squares show the trend for a tilted GPNJ, where a pronounced peak is seen in . With increasing tilt between the junction and the transport axis normal to the contacts, the incident angless of the incoming modes are effectively increased, leading to higher resistance. Such an increase in junction resistance is seen in a recent experiment (31) and a detailed theoretical treatement can be found in (16). Interestingly, the results suggest the absence of specular reflection from the edges of the graphene sample.
We will next apply the current density formalism in RGFA to show the electron trajectories inside the GPNJ device. The total current at an atomic site is equal to vector sum of all current components to nearest neighbor atoms (from Eq. 15). The classical ray tracing analysis from Snell’s law is shown in Fig. 8(a). For equal dopings on both sides (), the angle of refraction is exactly equal (with a negative sign) to the incident angle. As a result, a group of electrons originating from a point contact focus back to one point on the refracted side (6). Fig. 8(b) shows the trajectories when the doping at the refracted side is smaller than that at the incident side () making the critical angle smaller than . The geometrical ‘optics’ trajectories corresponding to electron focusing and total internal reflection are reproduced by the calculated current density withing an atomistic RGFA simulation, shown in Fig. 8(c), (d). A small source contact (10nm wide) is placed nm to the left of the pn junction, and electrons are injected with a small drain bias ( = 0.08 V) around a Fermi energy. When the gates are biased symmetrically around the junction, although the electrons see a voltage bias along the drain that spans the entire device width, the junction Hamiltonian and associated pseudospin conservation causes the electrons to focus to a small point at the drain. For an np junction, the electrons incident above the critical angle are unable to preserve their transverse quasi-momentum and reflect back, while those within critical angle tunnel to the opposite band on the other side.
We now investigate the quantum Hall plateaus for a single junction. According to recent experiment (41), the conductance in the junction regime follows the Ohmic conductance rule
(41) |
and are the filling factors in the two segments of the device. This was explained by Abanin et. al. (42) in terms of mode mixing at the interface for large systems with diffusive transport. As graphene filling factors are , it produces junction plateaus as for filling factor combinations of , etc. matching closely with the experiment (41).
This result changes considerably for coherent ballistic transport for smaller structures. Tworzydlo et. al. (43) showed that for armchair nanoribbons the conductance plateaus become
(42) |
independent of individual filling factors , . is the lowest Hall plateau for graphene. is the angle between valley isospins at the top and bottom edges and depends on the chirality of the nanoribbon. Depending on the number of hexagons along transverse direction, for and otherwise. This leads to
(43) |
The ballistic to ohmic cross over for the plateaus can be recovered by using interface and edge disorder (44).
The magnetotransport for zigzag ribbons cannot be explained with valley isospin arguments. Akhmerov et. al. (45) proposed the theory of valley valve effect with intervalley scattering which leads to complete suppression of conductance in the zigzag junction, if the number of atoms across the ribbon is even.
(44) |
Impact of depletion width (spacing between two gates) on GPNJ ballistic magneto-conductance for both armchair and zigzag ribbons are shown in Fig. 9. Magnetic field B = 10T and dopings are eV producing . Conductance values approach analytically predicted numbers only when the depletion width is sufficiently large ( 30nm), filtering out higher Landau Levels (LL) (44). It requires larger depletion width as individual filling factors get higher. In most cases, the lower plateaus ( and 0 for armchair and zigzag respectively) are more difficult to achieve than the higher one (2).
n-p-n junction
By applying a step like potential, we can realize an junction. The dopings can be done by either a combination of global back gate and top gate (40) or by selective chemical doping. In fig. 10(a), we show the energy band diagram for the device. We consider both junctions as abrupt, , while the extent of the type barrier region is . In fig. 10(c), we show an NEGF calculation of the npn conductance for various values. For larger , two Dirac points are evident from the pinched off conductance, but for smaller s we see a considerable tunneling around the second Dirac point ( = 0.5eV), where the critical angle for the incident electrons is supposed to be very small. The first dip happens at the Dirac point of the incident n region where as before, while the second corresponds to the electrons in the n region aligned with the Dirac point of the barrier p region and the transmission is small. In a single junction, the modes with higher angles than the critical angle reflect back as they do not have any propagating states to tunnel into while preserving their transverse quasimomentum. But in the junction case, the length of the forbidden region is finite, and the electrons even while aligned with the Dirac point of the central p region can still preserve their quasi-momenta by tunneling to the other side. This increased conductance has been seen experimentally (46) in the past (Fig. 10). In the limit when is very large, we approach the single junction case we worked out in Fig. 7(c) inset. And when is very small, we approach the uniform doping ( doped here) case through significant tunneling. Note that we also see oscillations for junctions, which originate from the Fabry Perot cavity formed by the electrostatic barrier (46); (47). Both the tunneling and the resonant oscillations can be captured analytically by matching eigenvector components across the pnp junction, as worked out in Ref. (5).
(45) |
giving us a transmission, . is the wave vector inside the barrier, , and and . We then sum over different modes, set by the width of the graphene sheet, to get the total conductance in units of (Fig. 10 (c), ,circles), in excellent agreement with the NEGF atomistic calculation and qualitatively with experiments (Fig. 10)).
3.3 KSF to simulate bulk graphene and GPNJ
The inset in Fig. 5(b) shows the calculated conductance of bulk graphene from KSF. Note that the peaks due to finite width quantization in Fig. 5 are removed by the application of periodic boundary conditions. The black circles show the trend that is expected from a linear approximation of the number of modes at low energy. The main figure in Fig. 5(b) shows the bulk graphene density of states calculated from three methods. The red line shows the result from an atomistic NEGF calculation. The black circles are from an analytical formula (48), while the blue line is from a numerical extraction of the density of states (DOS) from the graphene relation. All three calculations match very well with one another.
KSF approach also allows us to show the angle (mode) dependent chiral tunneling in graphene. We use Eq. 21 to get which is converted into a polar plot (Fig. 11) using , for transmission per spin per valley. For single layer graphene the transmission becomes unity for normal incidence (known as Klein tunneling). On the contrary, for bilayer graphene, the reflection becomes unity (Klein reflection). The tranmission ultimately depends on how well the wave-functions match on both sides of the junction (branches with same pseudospin components have same color in Fig. 11).
3.4 Combining KSF and RGFA
We now combine the two transport formalisms KSF and RGFA to a graphene junction and show how this method allows us to simulate very long, micron-sized devices with periodic boundary conditions along the transverse direction. We apply this technique to numerically reproduce the tunnel transmission across a GPNJ with split gates, which requires simulating a long device to extract the proper decay constant. We assume that the Dirac point varies linearly across the junction for simplicity (Fig. 7(a)) In presence of a slowly varying potential across the GPNJ, the transverse modes with finite wave-vector act like waveguides with a cut-off frequency and thus need to overcome a real tunnel barrier. A smooth GPNJ thus selectively transmits the low angle electrons (49), with suppression of the higher angle modes with increasing length of the junction transition region. The fundamental normal mode still does not see a barrier and thus transmits perfectly, reflection once again eliminated by a pseudospin dominated selection rule, retaining the Klein tunneling properties of GPNJ.
The conductance of a symmetric (equal and opposite doping on both sides) smooth GPNJ, using the expression from Eq. 37, becomes
(46) |
where is the effective critical angle imposed by the smooth junction (49) and is the angular separation between adjacent modes. Numerically capturing this exponential decrease in conductance with gate split distance requires simulating a long device. The definitions of the unit cell Hamiltonian are now changed by the following.
(47) |
where is the layer to layer coupling matrix in the transverse direction. This will describe a layer Hamiltonian with periodic boundary condition. Similarly we use to calculate as described in section II.
The vs is shown in inset of Fig. 12 with extended in the micron regime.
4 Conclusions - the role of numerical simulation
Graphene and its heterojunctions constitute a fascinating system that demonstrate rich physics such as anomalous quantum Hall, chiral tunneling, metamaterial behavior and pseudospintronics. Extended to other sub-systems such as bilayer graphene, this richness proliferates with Lifschitz transitions, anti-Klein tunneling, more anomalous Hall signatures and so on. The simplicity of the low-energy 2-D bandstructures of graphene and its progeny allows pen and paper demonstration of these physical concepts. The NEGF approach provides a common computational framework to verify these concepts atomistically, whereby a small modification to the Hamiltonian (e.g. a barrier or a field-induced phase in the coupling parameter) can manifest these physical concepts with very little effort, even in tougher geometries such as multiple junctions and tilted heterostructures. When comparing with experiments, however, one needs to worry about the sensitivity of these effects to scattering, edge states, quantization, smooth junctions, strain, roughness and so on. This is where computation plays a definitive role. It was hard to anticipate a-priori what the effect of specular edge scattering would be on the resistance of a tilted junction, but we now know through numerical simulations (16) that it creates a sweet-spot in the scattering profile that experiments do not show. Such ‘smoking-gun’ demonstrations allow us to judge the feasibility of devices that utilize unconventional electronic switching in graphene, such as the creation of transmission gaps with subthermal switching (15). We can now atomistically simulate graphene systems whose sizes are comparable to experimental dimensions, and directly demonstrate many of the physical effects that such devices are based on in presence of non-idealities. Further extensions will involve the role of self-consistent screening (Poisson), charge puddles, defects and incoherent scattering from remote optical phonons.
Acknowledgements.
The authors acknowledge financial support from INDEX-NRI. Redwan Sajjad would like to thank Khairul Alam and Golam Rabbani for the assistance in implementing RGFA. Authors also thank Eugene Kolomeisky, Branislav Nikolic, Farzad Mahfouzi for useful discussions.Footnotes
- email: redwan@virginia.edu
- email: redwan@virginia.edu
- email: redwan@virginia.edu
- journal: Journal of Computational Electronics
References
- K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric field effect in atomically thin carbon films,” Science, vol. 306, no. 5696, pp. 666–669, 2004.
- S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, D. C. Elias, J. A. Jaszczak, and A. K. Geim, “Giant intrinsic carrier mobilities in graphene and its bilayer,” Phys. Rev. Lett., vol. 100, no. 1, p. 016602, Jan 2008.
- K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. Stormer, “Ultrahigh electron mobility in suspended graphene,” Solid State Communications, vol. 146, no. 9-10, pp. 351 – 355, 2008.
- J. C. Slonczewski and P. R. Weiss, “Band structure of graphite,” Phys. Rev., vol. 109, pp. 272–279, 1958.
- M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, “Chiral tunnelling and the klein paradox in graphene,” Nat Phys, vol. 2, no. 9, pp. 620–625, 2006.
- V. V. Cheianov, V. Fal’ko, and B. L. Altshuler, “The Focusing of Electron Flow and a Veselago Lens in Graphene p-n Junctions,” Science, vol. 315, no. 5816, pp. 1252–1255, 2007.
- S. Datta, Electronic Transport in Mesoscopic Systems. Cambridge University Press, 1997.
- P. Damle, A. Ghosh, and S. Datta, “Unified description of molecular conduction: From molecules to metallic wires,” Physical Review B, vol. 64, no. 20, p. 201403, 2001.
- J. Guo, S. Datta, M. Anantram, and M. Lundstrom, “Atomistic simulation of carbon nanotube field-effect transistors using non-equilibrium greenâs function formalism,” Journal of Computational Electronics, vol. 3, pp. 373–377, 2004.
- S. Koswatta, S. Hasan, M. Lundstrom, M. Anantram, and D. Nikonov, “Nonequilibrium green’s function treatment of phonon scattering in carbon-nanotube transistors,” Electron Devices, IEEE Transactions on, vol. 54, no. 9, pp. 2339 –2351, sept. 2007.
- Y. Ouyang, Y. Yoon, and J. Guo, “Scaling behaviors of graphene nanoribbon fets: A three-dimensional quantum simulation study,” Electron Devices, IEEE Transactions on, vol. 54, no. 9, pp. 2223 –2231, sept. 2007.
- F. Tseng, D. Unluer, K. Holcomb, M. R. Stan, and A. W. Ghosh, “Diluted chirality dependence in edge rough graphene nanoribbon field-effect transistors,” Applied Physics Letters, vol. 94, no. 22, p. 223112, 2009.
- T. Low, S. Hong, J. Appenzeller, S. Datta, and M. Lundstrom, “Conductance asymmetry of graphene p-n junction,” Electron Devices, IEEE Transactions on, vol. 56, no. 6, pp. 1292 –1299, 2009.
- J. R. Williams, T. Low, M. Lundstrom, and C. Marcus, “Gate-controlled guiding of electrons in graphene,” Nat Nano, vol. 6, no. 4, pp. 222–225, 2011.
- R. N. Sajjad and A. W. Ghosh, “High efficiency switching using graphene based electron ”optics”,” Appl. Phys. Lett., vol. 99, no. 12, p. 123101, 2011.
- R. Sajjad, S. Sutar, J. U. Lee, and A. Ghosh, “Manifestation of chiral tunneling at a tilted graphene pn junction,” Phys. Rev. B, vol. 86, p. 155412, 2012.
- J. Wang, E. Polizzi, and M. Lundstrom, “A three-dimensional quantum simulation of silicon nanowire transistors with the effective-mass approximation,” Journal of Applied Physics, vol. 96, no. 4, pp. 2192–2203, 2004.
- M. Bescond, K. Nehari, J. Autran, N. Cavassilas, D. Munteanu, and M. Lannoo, “3d quantum modeling and simulation of multiple-gate nanowire mosfets,” in Electron Devices Meeting, 2004. IEDM Technical Digest. IEEE International, dec. 2004, pp. 617 – 620.
- M. Luisier, A. Schenk, W. Fichtner, and G. Klimeck, “Atomistic simulation of nanowires in the tight-binding formalism: From boundary conditions to strain calculations,” Phys. Rev. B, vol. 74, p. 205323, 2006.
- K. Nehari, N. Cavassilas, F. Michelini, M. Bescond, J. L. Autran, and M. Lannoo, “Full-band study of current across silicon nanowire transistors,” Applied Physics Letters, vol. 90, no. 13, p. 132112, 2007.
- M. Shin, “Quantum simulation of device characteristics of silicon nanowire fets,” Nanotechnology, IEEE Transactions on, vol. 6, no. 2, pp. 230 –237, march 2007.
- R. N. Sajjad, K. Alam, and Q. Khosru, “Parametrization of a silicon nanowire effective mass model from sp3d5s* orbital basis calculations,” Semiconductor Science and Technology, vol. 24, p. 045023, 2009.
- A. R. Rocha, V. M. Garcia-suarez, S. W. Bailey, C. J. Lambert, J. Ferrer, and S. Sanvito, “Towards molecular spintronics,” Nat Mater, vol. 4, no. 4, pp. 335–339, 2005.
- S. Salahuddin and S. Datta, “Self-consistent simulation of hybrid spintronic devices,” in Electron Devices Meeting, 2006. IEDM ’06. International, dec. 2006, pp. 1 –4.
- N. Mingo and L. Yang, “Phonon transport in nanowires coated with an amorphous material: An atomistic green’s function approach,” Phys. Rev. B, vol. 68, p. 245406, 2003.
- N. Mingo, “Anharmonic phonon flow through molecular-sized junctions,” Phys. Rev. B, vol. 74, p. 125402, 2006.
- X. Wang, Y. Ouyang, X. Li, H. Wang, J. Guo, and H. Dai, “Room-temperature all-semiconducting sub-10-nm graphene nanoribbon field-effect transistors,” Phys. Rev. Lett., vol. 100, no. 20, p. 206803, May 2008.
- M. Galperin, S. Toledo, and A. Nitzan, “Numerical computation of tunneling fluxes,” The Journal of Chemical Physics, vol. 117, no. 23, pp. 10 817–10 826, 2002.
- A. Svizhenko, M. P. Anantram, T. R. Govindan, B. Biegel, and R. Venugopal, “Two-dimensional quantum mechanical modeling of nanotransistors,” Journal of Applied Physics, vol. 91, no. 4, pp. 2343–2354, 2002.
- R. Lake, G. Klimeck, R. C. Bowen, and D. Jovanovic, “Single and multiband modeling of quantum electron transport through layered semiconductor devices,” Journal of Applied Physics, vol. 81, no. 12, pp. 7845–7869, 1997.
- S. Sutar, E. S. Comfort, J. Liu, T. Taniguchi, K. Watanabe, and J. U. Lee, “Angle-dependent carrier transmission in graphene pân junctions,” Nano Letters, vol. 12, no. 9, pp. 4460–4464, 2012.
- Y.-W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. H. Hwang, S. Das Sarma, H. L. Stormer, and P. Kim, “Measurement of scattering rate and minimum conductivity in graphene,” Phys. Rev. Lett., vol. 99, p. 246803, 2007.
- F. Miao, S. Wijeratne, Y. Zhang, U. C. Coskun, W. Bao, and C. N. Lau, “Phase-coherent transport in graphene quantum billiards,” Science, vol. 317, no. 5844, pp. 1530–1533, 2007.
- J.-H. Chen, C. Jang, S. Adam, M. S. Fuhrer, E. D. Williams, and M. Ishigami, “Charged impurity scattering in graphene,” Nat Phys, vol. 4, pp. 377–381, 2008.
- J. Tworzydło, B. Trauzettel, M. Titov, A. Rycerz, and C. W. J. Beenakker, “Sub-poissonian shot noise in graphene,” Phys. Rev. Lett., vol. 96, p. 246802, 2006.
- 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, vol. 438, no. 7065, pp. 197–200, 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, vol. 438, no. 7065, pp. 201–204, 2005.
- E. Mccann and V. I. Fal’ko, “Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer,” Phys. Rev. Lett., vol. 96, p. 086805, Mar 2006.
- B. Huard et al., “Transport measurements across a tunable potential barrier in graphene,” Phys. Rev. Lett., vol. 98, no. 23, p. 236803, Jun 2007.
- N. Stander, B. Huard, and D. Goldhaber-Gordon, “Evidence for klein tunneling in graphene p-n junctions,” Phys. Rev. Lett., vol. 102, p. 026807, Jan 2009.
- J. R. Williams, L. DiCarlo, and C. M. Marcus, “Quantum hall effect in a gate-controlled p-n junction of graphene,” Science, vol. 317, no. 5838, pp. 638–641, 2007.
- D. A. Abanin and L. S. Levitov, “Quantized transport in graphene p-n junctions in a magnetic field,” Science, vol. 317, no. 5838, p. 641, 2007.
- J. Tworzydło, I. Snyman, A. R. Akhmerov, and C. W. J. Beenakker, “Valley-isospin dependence of the quantum hall effect in a graphene junction,” Phys. Rev. B, vol. 76, p. 035411, 2007.
- T. Low, “Ballistic-ohmic quantum hall plateau transition in a graphene junction,” Phys. Rev. B, vol. 80, p. 205423, 2009.
- A. R. Akhmerov, J. H. Bardarson, A. Rycerz, and C. W. J. Beenakker, “Theory of the valley-valve effect in graphene nanoribbons,” Phys. Rev. B, vol. 77, p. 205416, 2008.
- B. Özyilmaz, P. Jarillo-Herrero, D. Efetov, D. A. Abanin, L. S. Levitov, and P. Kim, “Electronic transport and quantum hall effect in bipolar graphene junctions,” Phys. Rev. Lett., vol. 99, p. 166804, 2007.
- A. F. Young and P. Kim, “Quantum interference and klein tunnelling in graphene heterojunctions,” Nat Phys, vol. 5, no. 3, pp. 222–226, 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., vol. 81, no. 1, pp. 109–162, Jan 2009.
- V. V. Cheianov and V. I. Fal’ko, “Selective transmission of dirac electrons and ballistic magnetoresistance of junctions in graphene,” Phys. Rev. B, vol. 74, p. 041403, Jul 2006.