Coherent transport through graphene nanoribbons in the presence of edge disorder

We simulate electron transport through graphene nanoribbons of experimentally realizable size (length L up to 2 micrometer, width W approximately 40 nm) in the presence of scattering at rough edges. Our numerical approach is based on a modular recursive Green's function technique that features sub-linear scaling with L of the computational effort. We identify the influence of the broken A-B sublattice (or chiral) symmetry and of K-K' scattering by Fourier spectroscopy of individual scattering states. For long ribbons we find Anderson-localized scattering states with a well-defined exponential decay over 10 orders of magnitude in amplitude.


I. INTRODUCTION
The experimental realization of graphene, i.e., of a monolayer of carbon atoms 1-3 has opened up a rapidly developing field of fundamental and applied physics. The topology of the planar honeycomb lattice [figure 1(a)] with the resulting peculiar band structure near the K and K points [ figure 1(b), for a review, see 4,5 ] gives rise to many novel and intriguing physical properties, including the room temperature quantum Hall effect, minimum conductivity at the Fermi energy as well as possible applications for spintronics. Recent advances in fabricating width-modulated graphene nanoribbons helped to overcome intrinsic difficulties in creating tunnelling barriers and confining electrons in graphene, where transport is dominated by Klein tunneling-related phenomena 6,7 . Graphene quantum dots have been fabricated and Coulomb blockade [8][9][10] , quantum confinement 11 and charge detection 12 have been demonstrated.
The electronic properties of the perfect honeycomb lattice are meanwhile theoretically well understood 4 . However, in realistic graphene devices finite-size effects and imperfections play an essential role, especially for transport through confined structures. The importance of such effects results from the gapless band structure of graphene which does not allow straightforward confinement by electrostatic potentials. Devices thus have to be cut or etched resulting in rough edges. In turn, properties of the ideal graphene band structure cannot be invoked when simulating quantum transport through realistic devices in the presence of randomly shaped boundaries. Moreover, recent results underline that an incoherent Boltzmann transport approach, unlike a full quantum-mechanical calculation, fails to reproduce experimentally observed conductance signatures for impurity scattering 13 . However, the application of numerical methods for a full quantum mechanical simulation of graphene ribbons of realistic size constitutes a considerable challenge. A method of choice is the widely-used recursive Green's function technique 14 which is well suited to treat scattering structures such as wires or ribbons extended in one of their dimensions and usually implemented within the Landauer-Büttiker framework 15 for calculating transport coefficients. This technique has meanwhile been incorporated successfully into both an equilibrium [16][17][18][19][20] and a non-equilibrium [21][22][23][24] (Keldysh) description. Several variants of this method have been put forward employing specific symmetries of the system [25][26][27] if applicable, or using recursive algorithms [28][29][30][31][32][33] . In this article we present an extension of the modular recursive Green's function method (MRGM) 27 designed to treat graphene nano-ribbons with edge or bulk disorder. We address disorder scattering in graphene nanoribbons both on the microscopic level of specific lattice defects as well as on the macroscopic level of current measurements for which, as we will demonstrate, the details of the underlying disorder scattering play a crucial role. This paper is organized as follows: We briefly review key properties of the band structure of "ideal" infinitely extended graphene in the absence of disorder in section II. In section III, we introduce the application of the MRGM to finite-size graphene structures which allows us to treat extended structures efficiently due to the favorable scaling of the numerical effort with the linear dimensions of the ribbon. Applications to transport through rough-edged graphene nanoribbons will be presented in section IV followed by a short summary (section V).

II. TIGHT-BINDING SIMULATION OF GRAPHENE BAND STRUCTURE
The ideal, infinitely extended graphene sheet features a honeycomb lattice made up of two (A and B) interleaved triangular sublattices. It can be described in tight-binding (TB) approximation by the Hamiltonian 34 where the sum (i, j) extends over pairs of lattice sites, |φ j,s is the tight-binding orbital with spin s at lattice site j, V i is a locally varying potential (onsite energy), and γ i,j is the hopping matrix element between lattice sites i and j. Within our TB approximation, we include third-nearest-neighbour coupling [see figure 1(a)] using orthogonal tight-binding orbitals. This allows for four free parameters, namely the site-energy ε 0 and the overlap integrals γ i , i = 1, 2, 3, representing the interaction with the first, second and third nearest neighbour, respectively. We choose the γ i by fitting to ab-initio calculations, taken from Reich et al. [35][36][37] , arriving at γ 1 = −3.145, γ 2 = −0.042, and The dispersion near the non-equivalent K and K points [ figure 1(b)] resulting from the diagonalization of equation (1) features [for large distances k from K(K )] deviations from a perfect cone reflecting the influence of the hexagonal lattice. The cone becomes squeezed along the K − K directions, an effect known as trigonal warping 4,38 . Near the K point and for small k the band structure of equation (1) can be approximated (assuming that V i t) by a conical dispersion relation around the K point 39 , where we have set E(k K ) = 0. Note that the above expansion ignores both the length scale of the graphene lattice constant a = 1.4Å and the preferred directions of the lattice due to the discrete lattice symmetry. In this low-k limit, the Hamiltonian, equation (1), can be approximated by the Dirac Hamiltonian, with σ and τ being the Pauli spin matrices acting on the pseudo-spin and valley degrees of freedom. Analytic solutions for an infinitely extended graphene sheet described by equation (3) yields plane waves |k where the angle of the k vector θ k , connects relative amplitudes on the A and B sublattice 4 , Consequently, the pseudo-spin projection along the direction of propagation, the "helicity", is conserved reflecting the chiral symmetry of the ideal graphene sheet in the low-k limit. The additional degeneracy of two non-equivalent cones ("valleys") at the K and K points in the reciprocal lattice allows to formally represent the low-energy band structure near E = 0 in terms of Dirac-like four-spinors |ψ = (ψ K A , ψ K B , ψ K A , ψ K B ) with amplitudes for the A − B sublattice in real space and K −K in reciprocal space. The sign of θ k [equation (4)] is reversed upon transition from K to K . (Note that physical spin is not included in the present analysis.) This Dirac-like picture may serve as valuable starting point for the analysis of finite-size and edge effects on transport incorporated within the tight-binding Hamiltonian [equation (1)]. Chirality [equation (6)] is preserved in the presence of slowly (on the scale of the C-C bondlength) varying perturbations, suppressing backscattering 4 Conversely, rough edges may break chirality resulting in non-vanishing backscattering on the same cone and, at the same time, in coupling of the K and K cones. In the following, we will investigate in detail the effects of short-range defects 40 (e.g. edges) and deviations from the continuum Dirac picture on the transport properties of rough-edged graphene nanoribbons.

III. NUMERICAL METHOD
For the numerical treatment of finite-size graphene flakes and ribbons a number of simulation algorithms have meanwhile been proposed 14,[16][17][18][19][20][21][22][23][24] . We use in the following an extension of the Modular Recursive Green's function Method (MRGM) 26,27,41 applied to the third-order TB Hamiltonian [equation (1)]. The key idea of the MRGM is to break down a large device into independent smaller modules, each of which can be computed efficiently [see figure 2]. The Green's functions G 2 of the different modules with width W and length L are then combined to the desired device geometry using a small number of Dyson equations. In this Article, we introduce an efficient method to calculate the  (15)], a rectangular region can be separated from the half-infinite ribbon. (d) A rough-edged ribbon can now be assembled by randomly combining rectangles of variable length L2 and width W2. (e) For an arbitrarily shaped graphene scattering structure, the modular approach can still be used, although the Green's function of each module has to be calculated by direct inversion. G 2 : the associated numerical effort becomes independent of module length L. The algorithm involves the somewhat counter-intuitive steps to first calculate infinitely and semi-infinitely extended ribbons, i.e., modules of width W 2 but L = ∞, from which rectangular modules of finite length L 2 are "cut out" as needed by applying the Dyson equation "in reverse". The obvious advantage of this approach is that the computational effort becomes independent of L 2 . This approach is particularly advantageous for the simulation of weakly disordered graphene nano-ribbons: for weak disorder the spacing between individual defects both in the bulk and at the edges of the ribbon is large as compared to the lattice spacing. We can thus simulate the region between neighbouring defects by a single graphene module with perfect boundaries and place adjacent defects at the module boundaries [see red dotted lines in figure 2(d)]. This efficient calculation of the extended rectangular modules in between two defects is key to simulating large devices with length L total up to several micrometers (or ≈ 10 4 hexagons in one direction).
As a prototypical example, we build up an infinitely long nanoribbon with ideal zigzag boundaries (along thê x direction) by periodic repetition of a chain of carbon atoms of width W in transverse (ŷ) direction [see figure 2(a)]. Other geometries and boundaries can be treated analogously. The Hamiltonian H of the ribbon can thus be decomposed into a matrix H 0 describing the Hamiltonian of the vertical chain, and the coupling matrix H I describing the connection between two adjacent chains 42 , The solution of the Schrödinger equation for the infinite ribbon can be written in terms of an ansatz for a Bloch wave with |χ n the transverse eigenfunction. For expanding the Green's function in terms of χ n we need a complete set of transverse eigenfunctions including all evanescent modes in the sum [equation (9)] and the subsequent equations. The resulting generalized eigenvalue problem for e ik∆x and |χ n gives n left (right)-moving states |χ j (|χ  ), with corresponding momentum k j (k  ) in x direction. In the following, we introduce the shorthand notation D j (x) = |χ j e ikj x χ j | (D  (x) = |χ  e ikx χ  |) for the projections onto the right (left) moving Bloch states. From the Bloch states the Green's function of the infinite ribbon follows as 42 with the hopping matrix The Green's function of the half-infinite ribbon G L (G R ) extending from x 0 to −∞ (or +∞) can be written as with satisfying the boundary conditions G R,L (x, x ) = 0 for all x or x located at the end of the half-infinite ribbon (x, x = x 0 ). For an intuitive interpretation of equation (12) consider a disturbance from a point source at x . It reaches x by two paths: the direct propagation from x to x, given by the Green's function of the infinite ribbon, and the propagation from (14)], where the wave is reflected at the end of the ribbon and then propagates from (14)]. Note that the numerical effort to calculate G ∞ and G R,L is controlled by the transverse width of the ribbon W 2 and the number of transverse modes |χ n to be included while the x-dependence is given analytically. This scaling behaviour is key to calculate G 2 for the rectangular ribbon of arbitrary length L 2 by solving the Dyson equation in reverse for G 2 instead of for G L,R . Consequently, the numerical effort to calculate G 2 is independent of L 2 . In the final step, a rough-edged nanoribbon can now be assembled by successively joining rectangular ribbons G (i) 2 of varying length L 2 and width W 2 (average width W = 60nm) using the Dyson equation in forward direction, In our simulations, we treat ribbon lengths of several micrometers, and average over 100 random realizations ξ of edge roughness to eliminate non-generic features of particular ribbon configurations. To assemble such very long disordered ribbons, we start with a set of N B different modules M 1 , . . . , M N B and combine them to obtain a larger module M N B +1 . We connect the calculated modules in a random permutation P, formally equivalent to the composition rule of a Fibonacci sequence], creating an exponentially growing, pseudorandom sequence of modules. The interfaces between modules that include the disorder are randomly determined at each iteration step to avoid periodic repetition. If a more general shape of the scattering geometry is desired [e.g., for a non-separable disorder potential or curved boundaries, see figure 2(e)], the partitioning into modules and the subsequent efficient buildup of long structures is still readily possible. Only the first step of our algorithm has to be modified. The Green's function of individual modules is directly calculated by inversion, i.e., G = (E − H) −1 , using, e.g., a parallelized sparse-matrix solver 43 . Subsequent application of Dyson equations allows to assemble complex scattering geometries.
We find that the computing time τ of our approach scales as ∼ W 3 due to the cubic dependence of the eigenvalue problem and of the matrix multiplications. τ scales linearly with the number of building blocks N B used to set up the geometry, and logarithmically with L/L 2 . Numerically, we find for τ , with W given in nm. We have determined a prefactor a ≈ 5 from calculating the full scattering problem averaged over 100 configurations, when computing on 3 AMD Opteron processors (24 cores) at 2.2 GHz. Clearly, the prefactor strongly depends on the details of the employed hardware (i.e., network speed, cache size, compilation flags, etc.), while the scaling (17) does not. We note that the application of the algorithm presented here is not restricted to graphene nanostructures: any modular scattering system, which is build from modules along the lines of figure 2 (a)-(c) can be treated analogously. Possible applications include acoustic cavities, conventional semiconductors, topological insulators, or neutron scattering devices.

A. Transport coefficients
For conventional semiconductor heterostructures (e.g., quantum dots made of GaAs-AlGaAs), confinement is usually achieved by electrostatic gates resulting in smooth dot boundaries. Such confinement is not realizable for graphene due to its gapless band structure. While several theoretical concepts for opening a band gap have been proposed, the majority of experiments have achieved confinement by patterning of graphene nanodevices with oxygen plasma etching, chemical vapor deposition, specially prepared SiC substrates 44 , or chemical etching. These techniques, however, do not result in well-defined armchair or zigzag edges but in a rough-edge pattern featuring armchair and zigzag elements as well as adsorbates at the dangling carbon bonds [8][9][10]45 leading to an irregular edge structure. Edge effects can thus be expected to strongly influence the properties of graphene nanodevices.
We simulate the influence of edge scattering on transport through graphene nanoribbons by randomly varying the widths of the rectangular modules which build up the ribbon in the range W = 40±1 nm [see figure 3(a)]. Numerical tests show that a random sequence of N B = 5 different module widths represents a good compromise between a high degree of randomness and limited computational effort. In order to suppress correlations in the x-dependence of the roughness, the length of each rectangular module is chosen at random in the range of 0.24nm (one unit cell) to 10nm. We then use the above Fibonacci-like procedure to assemble a scattering geometry [see figure 3(b)] of up to several µm in length. Finally, all modules are connected to two ideal half-infinite graphene waveguides. We average over 100 realizations ξ of nanoribbons to eliminate non-generic features of particular ribbon configurations.
In addition to the quantization steps due to transverse confinement [dashed black line in figure 3(c)], a graphene nanoribbon with a perfect zigzag boundary of fixed width W features edge states with finite dispersion [figure 3(d)] since the coupling between the outermost carbon atoms is non-zero 5 . Consequently, the edge states of an ideal nanoribbon give rise to a peak in conductance G just below the Dirac point [shaded area in figure 3(c,d)]. In contrast to first-nearest-neighbor tight binding, and in line with the full ab-initio bandstructure and experiment 5 , our third nearest neighbor approach accounts for the breaking of electron-hole symmetry. Conductance is thus only approximately symmetric relative to E = 0.
In the presence of edge disorder, G undergoes several pronounced changes [figure 3(c)]: overall, G decreases with increasing distance in energy from the Dirac point relative to the ideal ribbon. The quantization steps due to the transverse confinement are strongly suppressed. Moreover, the edge disorder completely removes the sharp conductance peak attributed to edge states, as they are no longer conducting but become localized parallel to the ribbon 37 , i.e. along the direction of transport. Consequently, corresponding signatures are difficult to observe in transport measurements of realistic samples. Scanning tunneling spectroscopy provides an alternative approach: peaks in the local density of states at energies slightly below the Dirac point have been recently, indeed, observed in STS experiments 46 .
In the limit where quantization steps due to the transverse confinement are strongly suppressed [figure 3(c)] pronounced broad dips in transmission (see arrows in figure 3) replace the original steps in conductance. This counterintuitive reduction of transmission with increasing energy in the vicinity of steps can be qualitatively understood by considering Fermi's golden rule for the scattering of mode |nk into mode |n k 47,51 , Two trends contribute to this effect: firstly, strong fluctuations of the ribbon width broaden the DOS ρ n (E) and smoothen the steps. Secondly, as the ribbon locally narrows, backscattering via scattering into evanescent modes is enhanced. This occurs preferentially for energies close to the opening of a new mode and results in a reduction in transmission causing the dips. It is instructive to compare the present results to calculations for edge-and bulk-disordered semiconductor nanowires featuring a parabolic dispersion relation in the long-wavelength (continuum) limit. While a reduction of quantization steps by disorder is observed for edge-disordered semiconductor ribbons 49 resembling the present results, the prominent transmission dips observed for graphene (arrows in figure 3) appear not to be present in such a system (compare with figure 2 in 49 ). However, dips have been found in other disordered semiconductor nanostructures that are associated with resonances supported by attractive (bulk) disorder potentials 50 . The distance in energy between these resonances and the quantization step (i.e. the subband minimum) corresponds to the binding energy of the (quasi) bound state 50 . In the present case of graphene with rough edges, the enhancement of the local density of states near the Dirac point resulting from localized states at the edges is well known 19 . Their statistical weight has been found to be much higher for graphene than for conventional nanostructures 37 . In the absence of a local attractive potential, these localized states could take on the role of resonances: for edge structures similar to the ones we investigate, the resonance energies E i are statistically distributed in the range E i ∈ [−80, 0] meV 37 , and could, possibly, give rise to the observed broad dips.

B. Localized scattering states at edges
To gain a deeper understanding of the transport characteristics of edge-disordered graphene ribbons, we now analyze also individual scattering states. We find that states with energies where the conductance is only weakly perturbed by edge disorder [e.g. open triangle in figure 3(c)] feature a low amplitude at the edges [see figure 4(a)]. The overall probability density of these scattering states remains concentrated near the center of the ribbon [see figure 4(a)] and is therefore only slightly affected by edge disorder. When only a single mode is open in the leads (at energies close to the Dirac point), modes located at K and K in momentum space are not coupled in a zigzag graphene nanoribbon 4 . Only one cone contributes to transport in each direction [see figure 3(d)]. This imbalance in the number of leftand right-moving channels on each cone is a special property of zigzag graphene nanoribbons 52 , similar to the band structure of topological insulators. Backscattering is only possible in this energy window by inter-valley scattering at the rough edges. Since we observe a nearly perfectly conducting channel [ figure 5(b),(c)], inter-valley scattering that requires momentum transfers of the order |K − K | is obviously suppressed at low energies 53 .
By contrast, scattering states at energies where transmission (and conductance) is considerably reduced [solid triangle in figure 3(c) and figure 4(b)] feature a strong enhancement of their wavefunction near corners of the edges originating from one sub-lattice only [see arrows in figure 4(e, f)]. Projections onto the A and B sublattices [ figure  4 (e, f)] show pronounced differences reflecting the violation of pseudo-spin conservation [equation (7)]. We find enhancements of the A (B) sub-lattice scattering wave function at the upper (lower) edges of the ribbon, i.e., at those edges where the outermost carbon atom is of type A (B) [see zoom-ins in figure 4(d)-(f)], in line with a strong enhancement of the local DOS near rough edges 4,5,19 . However, the pronounced differences in the wavefunction patterns near the center of the ribbon are not accounted for only by localized edge states since their decay length into the ribbon interior is much smaller than the ribbon width. We therefore attribute the dramatic drop in conductance to pronounced intra-valley and inter-valley backscattering at the edge corners, since the suppression of backscattering associated with the conservation of pseudo-spin [equation (7)] no longer holds.
As reported in earlier work on edge disorder in rough-edged graphene nanoribbons, transmission is strongly suppressed close to the Dirac point, leading to the formation of a transport gap 17,19 . Atomic-scale defects on the edges of wide ribbons may lead to exponential (i.e., Anderson) localization due to destructive interference 19,49,53 . We use our modular approach to calculate scattering states on mesoscopic length scales [ribbon length L = 2µm, see figure 5(a, b)]. By averaging over many realizations of edge disorder, we can thus explicitly probe for exponential localization and determine the localization length. Looking at the longitudinal dependence of the scattering state, we observe an exponential decay over up to 10 orders of magnitude [ figure 5(c)]. Fitting to the functional form ψ (x) 2 ∝ exp(−x/l A ) we can numerically extract the localization length l A . We find l A to scale as l A ≈ αW/∆W , i.e., l A increases linearly with ribbon width and is inversely proportional to the disorder amplitude ∆W [ figure 5(d)]. The localization length l A is found to increase with increasing distance (in energy or in k) from the Dirac point (not shown), as suggested by the disorder-induced formation of a transport gap [17][18][19] . Superimposed on the exponential decay are oscillations on two shorter length scales: (i) a short beating period of λ = 0.7nm due to interference between the K and K cones 37 [λ in figure 5(a)] and (ii) a much slower variation with the length scale Λ ≈ 30nm [Λ in figure 5(a)] which corresponds to the wavelength Λ = 2π/k associated with the linear dispersion relation E = v F k, i.e., the distance in k space from the K point.
For comparison we also plot for the nearly perfectly conducting channel its wavefunction and its projection according to equation (19) [ figure 5(b),(c)]. If the incoming scattering wave couples to the near-perfectly conducting channel, this contribution will be dominant after a certain ribbon length as all other contributions quickly die out. The oscillations due to K − K interferences (i) are also present for this conducting state, though at reduced amplitude. While we observe Anderson localization for incoming scattering states at energies where more than one mode is open per cone, near-perfect conduction 52,53 appears to be confined to the topologically insulating part of the band structure. We expect these states to have a localization length that exceeds the dimension of our structure, if it is, at all, finite.

C. Variations of edge roughness
To investigate to what extent the above results depend on our particular choice of rectangular edge roughness we generalize our approach to include randomly jagged edges. We combine graphene segments featuring horizontal zigzag edges with segments featuring a boundary profile tilted by an angle β with respect to the horizontal zigzag direction [see insets in figure 6]. As outlined in section III, we calculate a set of modules (we use modules with length L ∈ [20,40]nm) by direct inversion of a finite-sized Hamiltonian, and combine these modules to efficiently generate very long structures (total length > 1µm).
Qualitatively, we find the same Anderson localization behaviour [see figure 6(a)] as a function of ribbon width W and roughness amplitude ∆W (for fixed β) as in the case with rectangular modules. However, unlike the case of free-particle dispersion 49 , graphene nanostructures feature an interesting interplay between lattice orientation and surface roughness. As this interplay is determined by the alignment angle β between the lattice orientation and the roughness, we can explicitly study its influence on transmission through the ribbon. We observe, indeed, that the value of the localization length strongly depends on the shape of the boundary with respect to the discrete lattice: edges consisting of randomly concatenated zigzag-edges only (i.e., with β = 60 • ) show substantially longer localization lengths than edges formed by an even mixture of zigzag and armchair edges [i.e., with β = 75 • ]. The dependence of the localization length on β can be understood in terms of the length of undisturbed zigzag (or armchair) edges: cutting close to a symmetry plane of the lattice (i.e. 60 • or 30 • ) results in comparatively longer segments of zigzag (or armchair) edges. By contrast, a cut at 75 • yields an irregular sequence of very short segments of armchair and zigzag boundaries and thus strongly breaks the translation symmetry of a clean zigzag (or armchair) edge. As an aside we note that the data presented in the previous subsection includes a variation in ribbon direction, i.e., the ribbons are not perfectly straight, notably in figure 5(a, b). This amounts to an effective increase of edge roughness. As a result, the localization length is further decreased in that case [compare figure 5(d) with figure 6 The observation of the relative change in resistance as a function of β, i.e., of the angle between the graphene lattice and the atomic-scale edge, might have implications for experiments. Measuring the atomic-scale roughness is difficult requiring an STM setup. Measuring localization length for different ribbon widths might provide an alternative probe for the atomic-scale edge roughness. Conversely, our results could be tested by comparing transport measurements for nanoribbons fabricated with different methods (i.e. etching, growth on Si-C substrates 44 , unzipping of graphene nanotubes 54 ) resulting in (known) different edge characteristics.

D. Fourier analysis of channel states
We explore now the interplay between short-range defects in real space and the absence (presence) of K − K inter-cone scattering in k-space. For this purpose we analyze the Fourier transforms of the asymptotic scattering state in the semi-infinite entrance (exit) waveguides. The Fourier transform is calculated as where we extend the integral over a finite area A in the asymptotic region of the waveguide, i.e., far away from the scattering region. Three different classes of asymptotic scattering states need to be considered: (i) the incoming Bloch states propagating in x-direction with wavenumber k n , ψ n (r) = e ikxnx χ n (y), where χ n (y) represents the transverse eigenfunction of mode n of the semi-infinite nanoribbon, (ii) the waves transmitted through the disordered region with transmission amplitude t mn , and (iii) the reflected waves with reflection amplitude r nm . The corresponding Fourier components are given by For a hexagonal lattice, the real and the reciprocal lattice are rotated by 90 degrees with respect to each other 56 . The first Brillouin zone for the ideal zigzag ribbon is thus given by a hexagon resting on a side rather than on a tip [see white hexagon in figure 7(a)]. To better visualize the enhancement of ψ(k) near the K or K points, we integrate ψ(k) over the transverse direction, In perfect zigzag ribbons, incoming modes feature non-vanishing amplitudes either near the K or the K point, i.e., there is no coupling (scattering) between K and K . For the incoming Bloch state the projected Fourier transform equation (24) features peaks at the k x values corresponding to K [see figure 7(d)]. The close-up of the peak in ψ(k) near the K point [inset of figure 7(a)] is structureless for the incoming Bloch wave with fixed transverse quantum number n. The horizontal and vertical lines are finite-size effects of the Fourier-transformed sample. Likewise, the origin of the additional bright spots inside the Brillouin zone is zone folding (the Brillouin zone of the ribbon is smaller than the graphene Brillouin zone). The interesting physics, on the other hand, is contained in finite amplitudes at both K and K points of the scattered wave [see figure 7(b,c)] which are induced by K−K scattering at rough edges. The relative strength of the integrated K and K peaks [see figure 7(e,f)] is a direct measure for the amount of inter-valley scattering. Furthermore, we observe a pronounced fine structure near the K and K points: enhancement along a half-circle forms around the Dirac points [see inset in figure 7(b,c)]. The surface of section of the double-cone band structure of constant energy is approximately a circle, the diameter of which is proportional to the energy. In the reflected (transmitted) part of the wavefunction, we only see the left half (right half) of this circle being populated, corresponding to negative (positive) group velocities. Enhancement along the full semicircle is due to inter-mode scattering n → m between transverse modes within the same valley (intra-valley scattering). We can thus conclude that pronounced intra-valley scattering at the rough edges distributes the reflected (or transmitted) wave almost uniformly over the energetically accessible half-circle of the band structure compatible with their propagation direction. The Fourier transform thus allows us to assess the amount of both inter-valley K −K scattering (by the relative amplitude around the K and K points in the reciprocal lattice) and inter-mode scattering by the angular distribution on the half-circle of a single cone for the incoming and reflected (or transmitted) states.

V. CONCLUSIONS AND OUTLOOK
We have presented a novel numerical approach to efficiently calculate the Green's function for extended nanoribbons. Key is the build-up of the ribbon by a random assembly of modules. Connecting these modules by way of a Dyson equation allows us to calculate the transport properties of long graphene ribbons. We find the conductance to be suppressed by rough edges relative to that of the perfect ribbon. For low energies, we observe near-perfectly conducting channels due to the band structure of zigzag graphene nanoribbons. Quantization steps are washed out and, in part, replaced by dips due to scattering into evanescent modes 48 , in contrast to edge-disordered semiconductor nanoribbons with free-particle dispersion 49 . An analysis of individual scattering states in both real space and Fourier space reveals pronounced A − B sublattice asymmetries and K − K scattering. We determined specific signatures of inter-and intra-valley scattering by Fourier transform spectroscopy of scattering states. We also identified Anderson localized states for different disorder configurations, extending over several micrometers with an exponential decay spanning 10 orders of magnitude. The corresponding localization length was calculated as a function of both the magnitude of edge roughness, and its alignment with the graphene lattice. We find that the latter plays a significant role in determining the localization length hinting at the importance of correctly modeling microscopic details of edge disorder beyond its amplitude and correlation length. We conclude by pointing to possible future applications. While early transport measurements were strongly affected by substrate interactions resulting in puddles of electron and hole conductivity due to bulk disorder, recent advances in the manufacturing of much cleaner graphene nanostructures by growing on Si-C substrates 44 , unzipping nanotubes to arrive at smooth-edged ribbons 54 as well as suspended graphene 57 have shifted the focus to edge disorder investigated in the present work. Indeed, the measurement of size quantization plateaus has been surprisingly elusive in graphene nanoribbons [58][59][60][61] , in qualitative agreement with our present findings. Only recently, by prolonged annealing and suspending graphene nanoconstrictions, first signatures of size quantization could be found 57 . Our findings regarding the sublattice sensitivity of the wavefunction (figure 3) could be tested by STM scans of bound states in graphene nano-islands 55 : in these measurements strong enhancements of wavefunction amplitudes were found on one sublattice, resulting in trigonal patterns close to edges that are quite similar to our numerical findings. Our simulations predict similar STM patterns for scattering states in extended nanostructures. In particular, such measurements could elucidate the precise nature of the scattering mechanisms encountered at edges prepared wth different techniques. Indeed, we expect defects affecting both sublattices, as recently investigated 40 , to exhibit different signatures than e.g. single vacancies. Finally, magnetic field effects allow for an additional external parameter more easily tunable in experiments than lattice geometries. States associated with different K points react differently to magnetic fields. This dependence might help to disentangle contributions from different K-points in scattering accessible by our Fourier analysis. We note that our algorithm can be easily adapted to accomodate magnetic fields while retaining its favorable scaling properties. Investigations in this direction are currently under way.