Spatial Current Patterns, Dephasing and Current Imaging in Graphene Nanoribbons

Using the non-equilibrium Keldysh Green's function formalism, we investigate the local, non-equilibrium charge transport in graphene nanoribbons (GNRs). In particular, we demonstrate that the spatial current patterns associated with discrete transmission resonances sensitively depend on the GNRs' geometry, size, and aspect ratio, the location and number of leads, and the presence of dephasing. We identify a relation between the spatial form of the current patterns, and the number of degenerate energy states participating in the charge transport. Furthermore, we demonstrate a principle of superposition for the conductance and spatial current patterns in multiple-lead configurations. We demonstrate that scanning tunneling microscopy (STM) can be employed to image spatial current paths in GNR with atomic resolution, providing important insight into the form of local charge transport. Finally, we investigate the effects of dephasing on the spatial current patterns, and show that with decreasing dephasing time, the current patterns evolve smoothly from those of a ballistic quantum network to those of classical resistor network.

Previous theoretical studies investigating the transport properties of GNRs have predominantly focused on the bias dependent conductance [1,2,3,4,5,6,7,8,10]. Spatial current patterns were investigated in the vicinity of the Fermi energy of half-filled GNRs in the wide-lead limit both for zero-magnetic field [3] and for non-zero magnetic fields and disorder [28]. Recently, however, it was argued [29] that in nanoscopic networks with narrow constrictions and/or narrow leads, spatial current patterns emerge which are qualitatively different from those in the wide-lead limit. In particular, these spatial current patterns exhibit clear signatures of quantum behavior: they (a) possess coherent "current riverbeds", i.e., spatial regions of large current density, whose widths are of the order of the Fermi wave-length, and (b) are strongly dependent on boundary conditions, such as the position of the leads, the geometry and size of the network, as well as the gate voltage. These results clearly suggest that understanding the quantum nature of local charge transport in GNRs, as reflected in the form of spatial current patterns, and its dependence on a GNR's geometry and size, or the presence of dephasing, and developing a method to visualize it, is of utmost importance for the further development of graphene based nanoelectronics and DNA sequencing.
In this article, we address this open question by demonstrating how the spatial current patterns in GNRs are determined by the interplay between the GNRs' geometry, aspect ratio and size, by the location and number of leads, and the presence of dephasing. By using the non-equilibrium Keldysh Green's function formalism [30,31] we find that the GNR's unconventional electronic structure [11] is reflected not only in the bias dependence of the conductance, but also in the form of the spatial current paths. In particular, we find that each of a GNR's discrete transmission resonances possesses a characteristic spatial current pattern (hence we refer to these resonances also as current eigenmodes in analogy to the equilibrium eigenmodes in the local density of states of nanoscopic systems). We demonstrate that the spatial form of these current patterns is determined by the number of degenerate states participating in the charge transport, which in turn allows us to predict how these current patterns evolve with increasing size or aspect ratio of a GNR. We show that current patterns can include closed loops of circulating currents, as well as exhibit backflow, i.e., the flow of charge through certain links opposite to the direction of the net charge flow [32]. We demonstrate that the spatial current paths are qualitatively different between leads attached to the zig-zag or armchair edges of the GNR and explore the form of charge transport parallel and perpendicular to edge states. Moreover, we show that while current patterns in certain 4-lead configurations can arise from the superposition of two 2-lead configurations, the conductance does not necessarily obey the superposition principle. Furthermore, we demonstrate that the quantum behavior of local charge transport can be visualized using scanning tunneling microscopy (STM) [33]: it represents an essentially non-intrusive method to image spatial current paths in GNR with atomic resolution, providing important insight into the form of local charge transport. Finally, we investigate more realistic models including the effects of dephasing and next-nearest-neighbor hopping, and show that with decreasing dephasing time, the current patterns evolve smoothly from those of a ballistic quantum network to those of a classical resistor network. These results provide important insight into the fundamental aspects of charge transport in nanoscopic graphene lattices.

Theoretical Model
To study electron transport in a graphene nanostructure, we consider two onedimensional leads coupled to a finite (N a × N z ) honeycomb lattice with N a hexagonal cells in the armchair direction, and N z cells in the zig-zag direction, shown schematically in figure 1. The graphene nanoribbon is described by the Hamiltonian [35] where −t r,r is the electronic hopping matrix element between sites r, r in the GNR, and c † r,σ creates an electron with spin σ at site r. In what follows, we take only the nearest neighbor hopping element, t, to be non-zero; the effects of a next-nearest neighbor hopping t will be discussed in Sec. 3.8. In order to account for the effects of dephasing, we assume that the electrons interact locally (with coupling constant g) with a phonon mode of energy ω 0 , as described by the last two terms in equation (1), where a † r creates a phonon at site r. The coupling between the GNR and the leads is described by the Hamiltonian where the primed sum runs over all sites r and l in the GNR and leads, respectively, that are coupled by a hopping element t h , and d † l,σ creates an electron with spin σ at site l in the leads. Below, we assume that each of the two leads is coupled to a single GNR site only, labeled L and R (see figure 1). Therefore, the only relevant property of the leads entering our calculations is the local Green's function at the lead sites that are coupled to the GNR. Finally, for the purpose of imaging spatial current patterns in the GNR [33], we consider the tunneling of electrons from an STM tip to a single site T in the GNR, a process which is described by the Hamiltonian where f σ annihilates an electron with spin σ in the STM tip. The spatial current patterns in a GNR are obtained by computing the current, I rr , between adjacent sites r,r in the GNR. This current is induced by different chemical potentials, µ L,R in the left and right leads, and given by [31] where G K rr is the full Keldysh Green's function between sites r and r which accounts for the electronic hopping within the GNR and between the GNR and the leads, as well as the electron-phonon interaction. We next introduce a matrix notation such that the Keldysh Green's function matrix is given byĜ K , and its (ij) element,Ĝ K ij , is the Keldysh Green's function between sites in the GNR denoted by i and j. In the absence of an electron-phonon interaction, one findsĜ Here,ĝ r,a,K are the diagonal retarded, advanced and Keldysh Green's function matrices, respectively, containing the Green's functions of the decoupled (t, t h = 0) GNR sites and leads: for the GNR sites, we have g r = 1/(ω + iδ) with δ = 0 + while for the leads, we assume the wide-band limit and take g r = −iπt −1 , yielding a constant density of states N L = t −1 in the leads. Finally,n F is a diagonal matrix containing the Fermi-distribution functions andt is the symmetric hopping matrix. Moreover, the local density of states at site r in the GNR is obtained via N (r, ω) = −ImG r rr (ω)/π. For the results shown below, we have used for numerical convenience δ = k B T = 10 −5 t, ∆µ = µ L − µ R = 2 × 10 −5 t and t h = 0.1t, unless otherwise specified. All spatial currents, I rr , presented below are normalized to the largest current in the GNR.
To account for the effects of the electron-phonon interaction for general g, ω 0 and T is computationally demanding and beyond the scope of this article. However, since we are mainly interested in the effects of dephasing on the spatial current patterns, we can consider the high-temperature approximation introduced in Ref. [34]. In this limit with k B T ω 0 , implying a large thermal population of the phonon mode, the calculation of the electronic self-energy correction is greatly simplified and computationally possible even for larger GNR sizes. Retaining in the Dyson equation only those terms that contain a factor of n B (ω 0 ) 1, the electronic self-energy in the self-consistent Born approximation is given by where i denotes a site in the GNR and α = K, r, a. Moreover, is the Keldysh phonon Green's function, which we assume to remain unchanged by the electron-phonon interaction, and n ph B (ω) is the phonon Bose distribution function (we assume that the phonons remain in thermal equilibrium). Note that due to the coupling to local phonon modes, the electronic self-energy is entirely local. A further simplification is achieved by considering the limit ω 0 → 0 such that the self-energy, to leading order in k B T /ω 0 , is given by As shown before [29,34], the solution of the Dyson equation for the retarded Green's function is then given bŷ whereD is a superoperator introduced in Ref. [34] which, when operating on a Green's function matrix, returns the same matrix with all elements set to zero except for the diagonal elements in the matrix that correspond to sites in the GNR, e.g., After self-consistently solving equation (8) forĜ r ,Ĝ K can be obtained in a closed expression viaĜ K =Ĝ rΣĜa where the diagonal matrixΣ is defined viã Here, the vector λ has components λ m =Λ mm withΛ = (ĝ r ) −1ĝ K (ĝ a ) −1 being a diagonal matrix whose only non-zero elementsΛ ii are those where i denotes one of the two leads. The matrixQ contains the elementŝ

Current Patterns in the (11 × 5) GNR
We begin by discussing charge transport, and in particular the relation between the conductance and spatial current patterns, in an (11 × 5) GNR in the absence of the electron-phonon interaction. In figure 2(a), we present the GNR's total conductance G(V c ) = I(V c )/∆V in the limit ∆V = (µ L − µ R )/e → 0 which is given by where V c = µ c /e is the bias midpoint with µ c = (µ L + µ R )/2, N 0 is the leads' local density of states, and G r L,R is the non-local Green's function between sites L and R where the current enters and exits the GNR (see figure 1), respectively [31]. Due to its finite size, the GNR possesses discrete energy levels [35] which for t h = 0 are given by As a result, the conductance exhibits discrete transmission resonances whenever µ c equals the energy a state whose wavefunction does not vanish at the site that the leads are coupled to. For each of these resonances, the current flowing through the GNR exhibits a different spatial pattern (the resonances are therefore also referred to as current eigenmodes [29], in analogy to eigenmodes in the density of states). The spatial current patterns [as obtained from equation (4)] for a number of current eigenmodes are shown in figures 2(b) -(h) [due to the particle-hole symmetry of the GNR's electronic structure, the conductance is symmetric around V c = 0 and the eigenmodes at ±V c exhibit identical current patterns; we therefore restrict our discussion below to the case V c ≥ 0]. Note that these eigenmodes can in general be accessed experimentally via gating of the GNR. Similar to prior results [29], we find that there exists current eigemodes which exhibit circulating current loops [as indicated by a yellow arrow in figure 2(f)] or which possess a net current through the GNR which is significantly smaller than some of the currents within the GNR [see figure 2(b)]. Of particular interest is the wellordered current pattern that occurs for V c = t/e shown in figure 2(d). To understand the spatial form of this current pattern, we note that charge transport for V c = t/e with leads attached to the zig-zag edge involves a 5−fold degenerate E = t state (not including N a = 5 states with vanishing amplitude along the zig-zag edge [35]), whose wave-vectors are shown in figure 3(a) (solid and open circles); for comparison, we also present the E = t equal energy contour for an infinitely large graphene sheet (solid line). Only three of these 5 states [green circles in figure 3(a)] possess a non-zero wavefunction at the GNR sites L and R that the leads couple to in figure 2 and thus contribute to the charge transport at V c = t/e. The Fermi velocity at E = t is given by and thus parallel to the primitive lattice vectors of one of the triangular sublattices of graphene. As follows from figure 2(d), the current propagates generally along the same direction. Thus, the spatial current pattern is reminiscent of the motion of a ballistic particle that is specularly reflected off the walls of the GNR. We find that this relation between the spatial current pattern and v F holds whenever there are degenerate states with the same v F participating in the charge transport. The relation between the spatial current pattern for V c = t/e, the aspect ratio of the GNR, and the degeneracy of the E = t state will be discussed in more detail in section 3.  Finally, we note that for finite-sized GNRs, there exist two states near zero energy at ±E l , that are delocalized along the zig-zag edge, but localized in the armchair direction [35]. For V c = ±E l /e, the GNR therefore exhibits strongly anisotropic transport properties in the armchair and zig-zag directions, as discussed in more detail in section 3.3.
3.2. Spatial current patterns, aspect ratios, and the degeneracy of the E = t state In the ballistic limit, spatial current patterns in nanoscale systems sensitively depend on boundary conditions such as the location of the leads or the geometry or aspect ratio of the system [29]. To investigate this issue in GNRs, we next study the dependence of spatial current patterns on the aspect ratio, N a /N z of the GNR, as well as the location of the leads. To exemplify this dependence, we consider the current eigenmode at V c = t/e, since it does not only exhibit a spatially highly ordered structure, but also occurs in all GNRs considered below. In figure 3, we present the evolution of the spatial current pattern with increasing N a in (N a × 5) GNRs for V c = t/e. We note for the discussion below that in all of these GNRs, there exist N a states at E = t whose wavefunction vanishes at the sites that the leads are attached to [35]; those states are therefore irrelevant for the purpose of the charge transport considered here and will not be discussed below.
The spatial current patterns in the GNRs exhibit two types of characteristic forms. In particular, in GNRs that satisfy N a +1 = p(N z +1) with integer p, the current follows the direction of the Fermi velocity, v F , being specularly reflected off the walls of the GNR, as shown in figures 3(b) and (e). In these GNRs, there are 3 states at E = t whose wavefunctions do not vanish at sites L and R; their wave-vectors are shown in figure 3(a) (green circles), together with the E = t equal-energy contour of an infinite graphene sheet. With increasing p, certain elements of the spatial current pattern are repeated [for example, the current pattern in figure 3(e) with p = 2 consists of two copies of the current pattern of figure 3(b) plus an additional current loop connecting these two elements]. In contrast, the GNRs shown in figure 3(c),(d), and (f) do not satisfy the above requirement, possessing only a single E = t state [whose wave-vector is shown as the open circle in figure 3(a)]. Consequently, these structures exhibit current patterns that are significantly more complex and possess a series of traits characteristic of quantum mechanical charge transport in nanostructures [29]. In particular, these current patterns exhibit circulating current loops that are isolated [see red arrow in figure 3(c)] and thus do not contribute to the net current through the GNR. Such current loops give rise to magnetic fields that could potentially be detected experimentally: we estimate that the magnetic field at the center of the current loop indicated by an arrow in figure 3(c) is approximately 5.5 × 10 −6 T. In addition, these current patterns exhibit backflow: while a net current flows from the left to the right (µ L > µ R ), there exist certain bonds in the GNR where the current flows from right to left (and thus opposite to the applied bias), such as the ones indicated by a yellow arrow in figure 3(d). We find that circulating current loops and backflow occurs in all GNRs in which the current pattern does not simply reflect the direction of the Fermi velocity, v F , as is the case in figures 3(b) and (e). Finally, in figure 3(g) we present the local density of states for the (11 × 5) GNR (with leads attached) at E = t, which has no resemblance to the form of the spatial current pattern shown in figure 3(e), a conclusion that was previously also drawn for nanoscale systems with a square lattice geometry [29]. This result stands in contrast to earlier findings relating either the spatial current pattern at V c directly to that of the modulus of the wave-function, |Ψ(r, E = eV c )| [36] (which possesses the same spatial structure as the density of states, N (r, E), since N (r, E) ∼ |Ψ(r, E = eV c )| 2 ), or the spatial flow of electrons to the spatial form of |Ψ(r, E = eV c )| 2 [38]. More generally, our results demonstrate that conductance variations measured in SPM experiments cannot simultaneously reflect the spatial flow of electrons [39] and the spatial form of |Ψ(r, E = eV c )| 2 [42].
Different types of current patterns emerge when the size of the GNR is increased along the zig-zag direction, as shown in figures 4(a) -(d). Here, we present the evolution of the current patterns with increasing N z in (7 × N z ) GNRs at V c = t/e. The current patterns for the (7 × 7) GNR (which possesses a 7-fold degenerate state at E = t, with 4 of these states possessing a non-zero wavefunction at GNR sites L and R) is similar to that of the (5 We can identify an interesting relation between the number of degenerate states at E = t, and the form of the current pattern: while the (7 × 9) and (7 × 13) GNRs possess only a single state at E = t, and their spatial current patterns exhibit some similar subpatterns, the (7 × 11) GNR possesses three degenerate states (two of which possess a non-zero wavefunction at sites L and R), resulting in a current pattern that is qualitatively different from that of the (7 × 9) and (7 × 13) GNRs. All of these current eigenmodes exhibit a different response to (symmetric) changes in the location of the leads, as shown in figures 4(e)-(h). For the (7 × 7) GNR, the current pattern simply deforms, but the current still propagates along the direction of the Fermi velocity. In contrast, for the (7 × 9), (7 × 11) and (7 × 13) GNRs, the current pattern undergoes a qualitative change when the leads are symmetrically displaced. This reiterates the fact that charge transport through confined systems in the quantum regime is highly sensitive to boundary conditions, and the spatial form of the resulting coherent current eigenmode need not necessarily conform to expectations based on v F . The question naturally arises whether the same type of spatial current patterns emerge when the leads are attached to the armchair edges. To investigate this question, we consider a (7 × 7) GNR [the same GNR as shown in figure 4(a)] and present in figure 5 the current patterns at V c = t/e for different locations of 4 leads attached to the armchair edges. When the leads are attached to the center sites, as shown in figure 5(a), the dominant contribution to the current flows along a line directly connecting the leads, rather than following the direction of the Fermi velocity. However, when the leads are displaced symmetrically from these high symmetry sites as shown in figures 5(b) and (c), the current flows again primarily along the direction of the Fermi velocity. It is interesting to note that in all three cases [figures 5(a)-(c)], the current exhibits two disjoint current paths, connecting either the left (leads 1 and 3) or right (leads 2 and 4) pair of leads. When one of these pairs of leads is removed (and the spatial symmetry of the system is thus broken), as shown in figures 5(d) where only leads 2 and 4 remain, the current pattern becomes rather diffuse (a very similar current pattern is found when only leads 2 and 3 remain). This not only demonstrates that the symmetry of the armchair edges requires two leads to be attached to either armchair edge for a well-defined current  pattern to emerge, but also that quantum interference effects between all four leads are crucial for creating the well-ordered current patterns shown in figures 5(b) and (c).

Current Through Localized Edge States
GNRs possess two low energy states near the middle of the band at ±E l , i.e., in close proximity to E = 0, that are delocalized along the zig-zag edge, and localized along the direction of the armchair edge [35]. The form of charge transport through these localized states in the wide-lead limit has recently attracted some attention [1,3,9], in particular in view of its possible application for DNA sequencing [23]. In figure 6(a), we present the local density of states for a (15 × 7) GNR at E = E l = 4.6 × 10 −13 t, that clearly demonstrates the localized nature of these low energy states. As a result, we expect that the charge transport at V c = E l /e is highly anisotropic. To investigate this anisotropy, we consider the (15 × 7) GNR with four leads attached to sites along their zig-zag edges, as shown in figures 6(b) and (c). When a voltage bias is applied between the armchair edges with V c = E l /e, (with chemical potential µ L for leads 3 and 4, and µ R for leads 1 and 2), the current flow is highly localized along the zig-zag edges, as shown in figure 6(b), where the E = E l state is delocalized. In this case, the conductance of the E = E l state (in the four-lead configuration) is equal to the quantum of conductance (similarly, the state at E = −E l possesses a conductance given by the quantum of conductance). This result holds for all values of N a or N z that we have considered. Moreover, for an infinitely large GNR in the armchair direction, i.e., N a → ∞, one has a doubly degenerate state at E l = 0 (one localized state at each of the zig-zag edges), and the conductance of the GNR at E = 0 is as expected equal to twice the quantum of conductance. Note that the spatial current pattern in figure 6(b) is qualitatively different from that obtained when wide leads are attached to the armchair edges: in this case, the largest current density occurs in the center of the GNR [1,3,9] and not along the zig-zag edges. On the other hand, when a voltage bias is applied between the zig-zag edges (with chemical potential µ L for leads 1 and 3, and µ R for leads 2 and 4), the current is strongly suppressed [see figure 6(c)], and scales as ∼ E 2 l . Since E l decreases exponentially with increasing length N a [35], we find that the current along the armchair direction also decreases exponentially with increasing length of the GNR, as expected from the localized nature of the ±E l -states.

Current Patterns in Four Leads Configurations
Do quantum mechanical currents obey the principle of superposition? In particular, can spatial current patterns in GNRs that are attached to four leads, simply be obtained by superposing two-lead current patterns such as the ones discussed above? To address this question we compare the spatial current patterns at V c = t/e in 4-lead and 2-lead configurations for two GNRs with different aspect ratios.
In figure 7, we present the current patterns for a (9×9) GNR in 4-lead configurations with the leads attached symmetrically around the center of the GNR, as shown in figures 7(a) -(c), or with two leads attached, as shown in figures 7(d) -(f). Since the GNR possesses 9 degenerate states at E = t, we again find that the current follows predominantly the direction of the Fermi velocity, thus propagating along the direction of the primitive lattice vector. Moreover, changing the separation between the leads simply deforms the current pattern, but does not introduce any qualitatively new spatial elements. A comparison of the 4-lead current patterns in figures 7(a) and (b) with the current patterns in a symmetric 2-lead configuration shown in figures 7(d) and (e), respectively, demonstrates that the 4-lead current pattern indeed represents a superposition (i.e., is the sum) of two 2-lead current patterns. Moreover, the current flowing in the 4-lead configuration is twice that of the 2-lead configurations. This implies that in the 4-lead configuration, there are no interference effects between the two superposing 2-lead current patterns, and that, in particular, no current flow occurs between the upper (leads 1 and 3) or lower (leads 2 and 4) pair of leads. We also note that the 4-lead current pattern in figure 7(b) is not the superposition of the 2-lead current pattern with just the upper two leads (leads 1 and 3) in figure 7(f), with its respective counterpart. In this case, the spatial current pattern is very diffuse, since the two leads cannot be connected by a current path that propagates along v F , and the total current is only approximately 58.5% of that in the 2-lead configuration of figure 7(e). In contrast, a (7 × 13) GNR possesses only a single state at E = t, and the spatial form and complexity of the current pattern at V c = t/e changes qualitatively when the separation between the leads on each side of the GNR is varied, as shown in figures 8(a)-(c). When the leads are attached to the corners [see figure 8(a)], the current flows along the edges and through the center of the GNR. At the same time, there are two large, circulating current loops in the upper and lower parts of the GNR, that do not contribute to the net current through the GNR. In contrast, when the separation between the leads is reduced [see figure 8(b)], the current flows along the direction of v F , which directly connects the leads, while most of the GNR does not exhibit any charge flow at all. Upon reducing the separation between leads even further, the current pattern again becomes more complex [see figure 8(c)]: the current now flows primarily along the outer perimeter, while connected to a larger circulating current loop in the center of the GNR. In figure 8(d), we present the corresponding 2-lead current pattern; the superposition of this pattern with its counterpart yields a spatial current pattern that is identical to that of figure 8(c). However, the total current through the 4-lead and 2-lead configurations is the same since [in contrast to the (9 × 9) GNR] there is only a single state that contributes to the charge transport at V c = t/e, thus limiting the conductance for both configurations to a single quantum of conductance. In other words, the spatial current in the 4-lead configuration is half the sum of the currents of the two 2-lead configurations. The degeneracy of the E = t state therefore does not only affect the spatial form of the current patterns in 2-lead configurations, but also whether the current patterns and total conductance in 4-lead configurations arises as a superposition of current patterns and conductances of 2-lead configurations.

Diamond-shaped GNRs
The observation that localized states exist near zig-zag edges in GNRs raises the interesting question of the form of localized state in other GNR geometries [1]. To investigate this question, we consider the diamond-shaped GNR shown in figure 9 whose edges are only of the zig-zag type. A plot of the local density of states at E = 0 [see figure 9(a)] reveals that a localized state exists, but that it is not confined to the zigzag edges of the GNR [as was the case for the GNR shown in figure 6(a)] but to its two corners with acute angles. Moreover, while the GNR possesses three degenerate states at E = t, we find that the current does not follow the direction of the Fermi velocity (which would require it to flow along the edges of the GNR), but rather flows through the GNR's center, as shown in figure 9(b). Finally, a plot of the local density of states at E = t in figure 9(c), which is the largest in the GNR's center, again reveals no resemblance to the spatial current pattern, as previously discussed in the context of figure 3.

Imaging of Spatial Current Patterns
The ability to experimentally visualize spatial current patterns is crucial for understanding and manipulating transport at the nanoscopic or atomic scale. In mesoscopic systems, such as quantum point contacts [37,38,39,40], quantum rings [41,42] and DNA [43], imaging of spatial current paths was successfully achieved using scanning probe microscopy (SPM) [44,45]. The spatial resolution of SPM, however, is insufficient to image spatial current patterns predicted to exist in nanoscopic systems [29,36,46] which vary on the atomic scale. Two of us therefore recently proposed a novel method for the imaging of spatial current paths in nanoscopic systems based on scanning tunneling microscopy (STM) [33], a method that can resolve spatial current patterns on the scale of a lattice constant. In addition, in the experimentally realized weak tunneling limit [47], this STM method only probes but essentially does not perturb the GNR's electronic structure, in contrast to SPM.
To demonstrate that this method for the spatial imaging of currents can also successfully be applied to GNRs, we consider as an example the current patterns in figures 4(a) and (e) for the (7 × 7) GNR, and in figure 3(e) for the (11 × 5) GNR. To visualize these spatial current patterns, we plot the spatial dependence of the current, I L,R (T), that flows from an STM tip (held above the GNR at potential µ T ) through the GNR into the left (L) or right (R) lead as a function of tip position T (for a more detailed description, see Ref. [33]). The tunneling of electrons from the STM tip into the GNR is described by the Hamiltonian in equation (3). Here, we set the bias between the leads to zero, i.e., µ L,R = µ 0 , and require that the bias difference between the tip and the system, ∆V T = (µ T − µ 0 )/e be equal to ∆V , and that the bias midpoint V T c = (µ 0 + µ T )/2e be equal to V c (here, V c and ∆V were used to generate the actual current patterns). In figures 10(a) -(c), we present the resulting STM image, i.e., the spatial plots of I L (T), which correspond to the current patterns in figures 4(a) and (e) and in figure 3(e), respectively. A comparison of the actual current patterns with the STM images shows that the latter provide an accurate and atomically resolved image of the spatial current patterns. This agreement holds even when the leads are not located at high symmetry positions, as shown in figures 4(e) and 10(b), or when the size of the system increases, and the current pattern becomes repetitive, as in figures 3(e) and 10(c). Moreover, we demonstrated in Ref. [33], that the form of the total conductance, G(V c ), is an independent, experimentally verifiable, criterion for the success of the STM method in imaging spatial current patterns. In particular, we found that when G(V c ) is close to the maximal allowed conductance, and varies weakly around a given V c , the STM method successfully images the spatial current patterns at V c in the entire nanoscopic system. On the other hand, if G(V c ) is sharply peaked, the agreement breaks down in at least part of the network. The (7 × 7) GNR considered here satisfies this criterion, as follows from figure 10(d), where we present its conductance G(V c ): it is close to the quantum of conductance at V c = t/e (see arrow) where we image the current pattern and its width (as a function of V c ) is larger than the distance between two peaks. Note that the zeros in the conductance G(V c ) correspond to zeroes of the real part of the retarded Green's function for t h = 0 which occur near, but never exactly on, conductance resonances. More generally, we find that the above criterion for G(V c ) is satisfied in nanoscopic networks which (a) possess a continuum of energy eigenstates through which transport can take place or (b) are strongly coupled to leads.
Finally, we note in passing that understanding the magnitude and spatial paths of currents injected from the STM tip into the GNR, will likely also have relevance for DNA sequencing where currents are induced locally by pulling a DNA through a hole in a GNR [20,21,22,23].

Dephasing and the Classical Limit
To study the effects of dephasing on the form of spatial current patterns, we consider the coupling of electrons to a local phonon mode [48,49], as described by the Hamiltonian in equation (1). As discussed in section 2, we employ the high-temperature approximation, where the dephasing is controlled by a single parameter, γ [see equation (7)]. In figure 11 we present the evolution of the spatial current pattern with increasing γ for an (11 × 5) GNR that contains a constriction shown in gray. In the absence of dephasing, i.e., for γ = 0, [see figure 11(a)], the current pattern is similar to that of the (11×5) GNR without a constriction, as shown in figure 3(e). We note that while the constriction is located in a region of the GNR through which no current flows in figure 3(e), it nevertheless leads to a change in the actual current pattern. As γ increases, the current pattern becomes more diffuse, and evolves smoothly from the ballistic limit [see figure 11(a)], to that of a classical resistor network [see figure11(d)]. Moreover, with increasing γ, the current pattern remains well defined in the vicinity of the leads, but becomes more diffuse as one moves further away from the leads [see regions indicated by yellow arrows in figure 11(c)]. This is expected since the diffuse current pattern arises from multiple scattering of the electrons off phonons, and the resulting dephasing. While required by symmetry, it is nevertheless interesting to note that the effects of dephasing are reversed (with the Figure 11. Evolution of spatial current patterns with increasing dephasing parameter, γ for an (11 × 5) GNR with a constriction (shown in gray): (a) γ = 0 with mean free path l = ∞ (b) γ = 0.005t 2 , l = 11.5a 0 , (c) γ = 0.02t 2 , l = 5.8a 0 , (d) γ = 0.5t 2 , l = 1.2a 0 . spatial current pattern again becoming more coherent) as the current approaches the site where it exits the GNR. This phenomenon of current patterns exhibiting coherent (well-defined) and incoherent (diffusive) spatial regions occurs when the mean-free path is considerably shorter than the size of the GNR [in figure 11(c), one has l = 5.8a 0 ], but larger than a few lattice spacings. As the mean-free path becomes shorter, the region around the leads where the current forms a well defined, coherent pattern shrinks and eventually vanishes [see figure 11(d)].

Effects of Next-Nearest Neighbor Hoppings on Current Patterns
In the preceding sections, we neglected the effects of a possible next-nearest neighbor hopping, t , on the form of the electronic excitation spectrum and the resulting conductance and spatial currents patterns in GNRs. The question naturally arises to what extent a non-zero t , with estimates ranging from t = 0.02t to t = 0.2t [51], qualitatively affects the results presented above.
In macroscopic graphene lattices, the existence of a non-zero t is often neglected, not only because of the large uncertainty in its value, but also because its effects on the electronic excitation spectrum are rather trivial. In particular, the energy dispersion for a non-zero t is given by [50] where f k = 2 cos( √ 3k y a 0 ) + 4 cos As a result, all momentum states located at an energy E 0 for t = 0 are simply shifted to a new energy for t = 0. This, in particular, implies that the degeneracy of states is not lifted by a non-zero t . For finite-size GNRs, however, the situation is different since a non-zero t lifts the degeneracy of the states which for t = 0 are located at E = t states. The extent of the energy splitting between these states is non-universal and depends not only on the aspect ratio of the GNR, but also on its overall size, since, according to Eq.(15), the splitting vanishes in the limit of infinitely large graphene layers where degeneracy is restored. In figure 12(a), we present a conductance scan for a (5 × 5) GNR with t = 0.1t over an energy range that contains all three states (and only those) that are located at E = t for t = 0, and possess a non-zero wavefunction at L and R. Choosing a left and right chemical potential that brackets these three states [see vertical red lines in figure 12(a)], we find that the resulting spatial current pattern shown in figure 12(b) still agrees very well with that obtained for the case t = 0, as follows from a comparison of figures 12(b) and 3(b). The same conclusion also holds for larger GNRs and different values of t , as follows from a comparison of the current pattern for a (11 × 5) GNR with t = 0.06t shown in figure 12(c) with that for t = 0 shown in figure 3(e). This implies that for both cases, t = 0 and t = 0, the same (well-ordered) spatial current pattern occurs if the left and right chemical potentials are chosen such that they bracket all states which for t = 0 occur at E = t.
Finally, we found that a non-zero t shifts the energy of the state localized along the zig-zag edge near E ≈ 0 for t = 0, but does not change its localized nature or spatial structure. In figure 12(d) we present the local density of states for t = 0.06t at E l ≈ 0.128t, which demonstrates that this state remains localized along the zig-zag edge. Figure 12. (a) Conductance trace near the states which occur at E = t when t = 0. If µ L,R are chosen to bracket all three conductance resonances, well-ordered current patterns are produced for the (b) (5 × 5) GNR at t = 0.1t with δ = 10 −11 t, ∆µ = 0.1t, and V c = 1.17t/e, and (c) (11 × 5) GNR at t = 0.06t with δ = 10 −13 t, ∆µ = 0.055t, and V c = 1.1075t/e. (d) The LDOS for the edge states at energy E l ≈ 0.128t.

Summary
In this article, we investigated the quantum nature of local charge transport in graphene nanoribbons, as reflected in the form of spatial current patterns. By using the nonequilibrium Keldysh Green's function formalism [30,31], we demonstrated that the form of spatial current patterns is determined by the interplay between the GNRs' geometry, size and aspect ratio, by the location and number of leads, and the presence of dephasing. In particular, we identified a crucial relation between the spatial form of current patterns, and the number of degenerate states participating in the transport. This insight, in principle, provides us with the opportunity to predict and customdesign spatial current patterns in GNRs. Furthermore, we showed that the number of degenerate states participating in charge transport also determines whether current patterns and conductances in GNRs with 4-lead configurations can be considered as arising from the superposition of current patterns and conductances in 2-lead configurations. In addition, we showed that spatial current patterns in GNR can be spatially imaged with atomic resolution using scanning tunneling microscopy, allowing us to gain unique insight into the nature of local charge transport. We also demonstrated how spatial current patterns evolve with increasing dephasing from the ballistic limit, where they are coherent and spatially well defined, to the classical limit, where they are incoherent and spatially diffuse. Finally, we showed that our conclusions remain valid in the presence of a realistic next-nearest neighbor hooping element. These results represent an important first step towards understanding and manipulating charge transport at the atomic level in GNRs, which are a necessary requirement for the further development of graphene based nanoelectronics and DNA sequencing.