Coherent transport through graphene nanoribbons in the presence of edge disorder

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 have helped to overcome intrinsic difficulties in creating tunnelling barriers and confining electrons in graphene, where transport is dominated by the Klein tunnelling-related phenomena [6, 7]. Graphene quantum dots have been fabricated and Coulomb blockade [8–10], quantum confinement [11] and charge detection [12] have been demonstrated.

Zoom In Zoom Out Reset image size

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 (GF) technique [14], which is well suited to treating 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–20] and a non-equilibrium [21–24] (Keldysh) description. Several variants of this method have been put forward employing specific symmetries of the system [25–27] if applicable, or using recursive algorithms [28–33]. In this paper, we present an extension of the modular recursive Green's function method (MRGM) [27] designed to treat graphene nanoribbons with edge or bulk disorder. We address disorder scattering in graphene nanoribbons both on the microscopic level of specific lattice defects and 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 2. In section 3, we introduce the application of the MRGM to finite-size graphene structures which allows us to treat extended structures efficiently due to the favourable 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 4 followed by a short summary (section 5).

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 the tight-binding (TB) approximation by the Hamiltonian [34]

$$\begin{equation} H = \sum_{i,s}\left| \phi_{i,s} \right>V_i\left< \phi_{i,s} \right|- \sum_{(i,j),s}\gamma_{i,j}\left| \phi_{i,s} \right>\left< \phi_{j,s} \right| + \mathrm{h.c.}\,, \end{equation} \tag{ 1 }$$

where the sum (i,j) extends over pairs of lattice sites, $\left | \phi _{j,s} \right >$ is the TB orbital with spin s at lattice site j, V_i is a locally varying potential (on-site 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 TB orbitals. This allows for four free parameters, namely the site-energy ε₀ and the overlap integrals γ_i, i = 1,2,3, representing the interaction with the first-, second- and third-nearest neighbours, respectively. We choose the γ_i by fitting to ab initio calculations, taken from Reich et al, arriving at γ₁ = −3.145, γ₂ = −0.042 and γ₃ = −0.35 [35–37].

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],

$$\begin{equation} E(k + k_K) = E(k_K) + k \partial_k E(k_K) + \mathcal{O}(k_K^2)\approx v_{\mathrm{F}}\hbar|k|, \end{equation} \tag{ 2 }$$

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

$$\begin{equation} H=v_F\left(p_x\sigma_x\otimes\tau_0+p_y\sigma_y\otimes\tau_z\right) \end{equation} \tag{ 3 }$$

with $\vec {\sigma }$ and $\vec {\tau }$ 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) yield plane waves |k〉 where the angle of the k vector θ_k,

$$\begin{equation} \theta_k = \tan^{-1}(k_y / k_x), \end{equation} \tag{ 4 }$$

connects relative amplitudes on the A and B sublattice [4],

$$\begin{equation} \left| \mathbf{k} \right>= \mathrm{e}^{\mathrm{i} \mathbf{k\cdot r}}\left(\left| A \right> + \mathrm{e}^{\mathrm{i} \theta_k}\left| B \right>\right)/\sqrt 2. \end{equation} \tag{ 5 }$$

Consequently, the pseudo-spin projection along the direction of propagation, the 'helicity',

$$\begin{equation} \hat{h} = (\mathbf{\sigma}\;\cdot\;\mathbf k)/\left|k\right|, \end{equation} \tag{ 6 }$$

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 us 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 a valuable starting point for the analysis of finite-size and edge effects on transport incorporated within the TB Hamiltonian (equation (1)). Chirality (equation (6)) is preserved in the presence of slowly (on the scale of the C–C bond length) varying perturbations, suppressing backscattering [4]

$$\begin{equation} P(\mathbf{k}\rightarrow\mathbf{k}') =\left|\left< \mathbf k' \right|V\left| \mathbf k \right>\right|^2 \propto \cos^2[(\theta_{\mathbf k} - \theta'_{\mathbf k'})/2]. \end{equation} \tag{ 7 }$$

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.

For the numerical treatment of finite-size graphene flakes and ribbons, a number of simulation algorithms have, meanwhile, been proposed [14, 16–24]. We use in the following, an extension of the modular recursive GF 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 GFs $G_\Box$ of the different modules with width W and length L are then combined with the desired device geometry using a small number of Dyson equations. In this paper, we introduce an efficient method to calculate the $G_\Box$: 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_{\Box }$ but L = ∞, from which rectangular modules of finite length $L_\Box$ 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_{\Box }$. This approach is particularly advantageous for the simulation of weakly disordered graphene nanoribbons: for weak disorder the spacing between individual defects both in the bulk and at the edges of the ribbon is large as compared with 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 the 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 micrometres (or ≈10⁴ hexagons in one direction).

Zoom In Zoom Out Reset image size

As a prototypical example, we build up an infinitely long nanoribbon with ideal zigzag boundaries (along the $\hat x$-direction) by periodic repetition of a chain of carbon atoms of width W in the transverse ($\hat y$) 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₀ describing the Hamiltonian of the vertical chain and the coupling matrix H_I describing the connection between two adjacent chains [42],

$$\begin{equation} H = \left(\begin{array}{ccccc} \ddots & & & & \\ & H_0 & H_I & 0 & \\ & H_I^\dagger & H_0 & H_I & \\ & 0 & H_I^\dagger & H_0 & \\ & & & & \ddots\\ \end{array}\right). \end{equation} \tag{ 8 }$$

The solution of the Schrödinger equation for the infinite ribbon can be written in terms of an ansatz for a Bloch wave

$$\begin{equation} \left| \psi \right> = \sum_n \mathrm{e}^{\mathrm{i}k_nx}\left| \chi_n \right> \end{equation} \tag{ 9 }$$

with $\left | \chi _n \right >$ the transverse eigenfunction. For expanding the GF 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 $\left | \chi _n \right >$ gives n left (right)-moving states $| \chi _j \rangle$ ($| \chi _{{\overline \jmath }} \rangle$), with the corresponding momentum k_j ($k_{{\overline \jmath }}$) in the x-direction. In the following, we introduce the shorthand notation $D_j(x) = \left | \chi _j \right >\mathrm {e}^{\mathrm {i} k_j x}\left < \chi _j \right |$ ($D_{\overline \jmath }(x) = \left | \chi _{\overline \jmath } \right >\mathrm {e}^{\mathrm {i}k_{\overline \jmath } x}\left < \chi _{\overline \jmath } \right |$) for the projections onto the right (left) moving Bloch states. From the Bloch states the GF of the infinite ribbon follows as [42]

$$\begin{equation} G^{\infty}(x,x') = \left\{ \begin{array}{rr} \sum_{j=1}^N D_j(x-x') V^{-1}, & x\geqslant x',\\ \\ \sum_{j=1}^N D_{\overline\jmath}(x-x') V^{-1}, & x\geqslant x',\\ \end{array} \right. \end{equation} \tag{ 10 }$$

with the hopping matrix

$$\begin{equation} V=\sum_{j=1}^N H_I^\dagger \left[ D_j(-\Delta x) - D_{\overline\jmath}(-\Delta x)\right]. \end{equation} \tag{ 11 }$$

The GF of the half-infinite ribbon G_L(G_R) extending from x₀ to −∞ (or +∞) can be written as

$$\begin{equation} G_{\mathrm{R,L}}(x,x') = G^{\infty}(x,x') + G^0_{\mathrm{R,L}}(x,x') \end{equation} \tag{ 12a }$$

with

$$\begin{eqnarray} G_{\mathrm{L}}^0(x,x') &=& \sum_j D_{\overline\jmath}(x-x_0) D_j(x_0-x') V^{-1}, \end{eqnarray} \tag{ 12b }$$

$$\begin{eqnarray} G_{\mathrm{R}}^0(x,x') &=& \sum_j D_j(x-x_0) D_{\overline\jmath}(x_0-x') V^{-1}, \end{eqnarray} \tag{ 12c }$$

$D_{\overline \jmath }(x-x_0)$

$W_\Box$

$\left | \chi _n \right >$

$G_{\Box }$

$L_{\Box }$

$$\begin{equation} G_{\mathrm{L,R}} = G_\Box + G_\Box H_I G_{\mathrm{L,R}} \end{equation} \tag{ 13 }$$

in reverse for $G_\Box$ instead of for G_L,R. Consequently, the numerical effort to calculate $G_\Box$ is independent of $L_\Box$. In the final step, a rough-edged nanoribbon can now be assembled by successively joining rectangular ribbons $G_\Box ^{(i)}$ of varying length $L_\Box$ and width $W_\Box$ (average width $\overline W = 60\,\mathrm {nm}$) using the Dyson equation in the forward direction,

$$\begin{equation} G = G_\Box^{(i)}+G_\Box^{(i)}H_IG. \end{equation} \tag{ 14 }$$

In our simulations, we treat ribbon lengths of several micrometres, 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₁,...,M_NB and combine them to obtain a larger module M_NB + 1. We connect the calculated modules in a random permutation $\mathcal {P}$, $M_{N_{\mathrm {B}}+1}=\mathcal {P}(M_1+\cdots +M_{N_{\mathrm {B}}})$ (e.g. for N_B = 5, M₆ = M₃ + M₁ + M₅ + M₄ + M₂). This procedure is repeated iteratively (i.e. $M_7 = \mathcal {P}(M_2 + \cdots + M_6)$, formally equivalent to the composition rule of a Fibonacci sequence), creating an exponentially growing, pseudo-random 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 build-up of long structures is still readily possible. Only the first step of our algorithm has to be modified. The GF of individual modules is directly calculated by inversion, i.e. G = (E − H)⁻¹, using, e.g., a parallelized sparse-matrix solver [43]. Subsequent application of the Dyson equations allows us to assemble complex scattering geometries.

We find that the computing time τ of our approach scales as ∼W³ 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_\Box$. Numerically, we find, for τ,

$$\begin{equation} \tau\, (\mathrm{ms}) \sim a \; N_{\mathrm{B}}\,\ln\left(\frac{L}{L_\Box}\right) W^3, \end{equation} \tag{ 15 }$$

with W given in nm. We have determined a prefactor a ≈ 5 by calculating the full scattering problem averaged over 100 configurations, when computing on three 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 (15) does not.

We note that the application of the algorithm presented here is not restricted to graphene nanostructures: any modular scattering system which is built from modules along the lines of figures 2(a)–(c) can be treated analogously. Possible applications include acoustic cavities, conventional semiconductors, topological insulators or neutron scattering devices.

4.1. 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 vapour 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–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.24 nm (one unit cell) to 10 nm. 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.

Zoom In Zoom Out Reset image size

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 figures 3(c) and (d)). In contrast to first-nearest-neighbour TB, and in line with the full ab initio band structure and experiment [5], our third-nearest-neighbour 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 tunnelling spectroscopy provides an alternative approach: peaks in the local density of states at energies slightly below the Dirac point have indeed been recently 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 counter-intuitive 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 $\left | n k \right >$ into mode $\left | n' k' \right >$ [47, 51],

$$\begin{equation} \Gamma(E) \propto \sum_{n'} \left|\left<nk|H|n'k'\right>\right|^2 \rho_{n'}(E). \end{equation} \tag{ 16 }$$

Two trends contribute to this effect: firstly, strong fluctuations of the ribbon width broaden the density of states (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 with 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.

4.2. Localized scattering states at edges

To gain a deeper understanding of the transport characteristics of edge-disordered graphene ribbons, we now analyse 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 centre 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 left- and 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 (figures 5(b) and (c)), inter-valley scattering that requires momentum transfers of the order $\left |K-K'\right |$ is obviously suppressed at low energies [53].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By contrast, scattering states at energies where transmission (and conductance) is considerably reduced (solid triangle in figure 3(c) and 4(b)) feature a strong enhancement of their wavefunction near corners of the edges originating from one sub-lattice only (see arrows in figures 4(e) and (f)). Projections onto the A and B sublattices (figures 4(e) and (f)) show pronounced differences reflecting the violation of pseudo-spin conservation (equation (7)). We find enhancements of the A (B) sub-lattice scattering wavefunction 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 figures 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 centre 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 instead 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 figures 5(a) and (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,

$$\begin{equation} |\bar{\psi}(x)|^2 = \int_0^W |\psi(x,y)|^2\, \mathrm{d}y, \end{equation} \tag{ 17 }$$

we observe an exponential decay over up to ten orders of magnitude (figure 5(c)). Fitting to the functional form $\left |\bar {\psi }(x)\right |^2\propto \exp (-x/l_{\mathrm {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–19].

Superimposed on the exponential decay are oscillations on two shorter length scales: (i) a short beating period of λ = 0.7 nm due to interference between the K and K' cones [37] (λ in figure 5(a)) and (ii) a much slower variation with length scale Λ ≈ 30 nm (Λ 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 (17) (figures 5(b) and (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, although 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.

4.3. 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 3, 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).

Zoom In Zoom Out Reset image size

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 include a variation in ribbon direction, i.e. the ribbons are not perfectly straight, notably in figures 5(a) and (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(b) for β = 90°).

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 a scanning tunnelling microscope. Measuring localization lengths 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] and unzipping of graphene nanotubes [54]) resulting in (known) different edge characteristics.

4.4. 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 analyse the Fourier transforms of the asymptotic scattering state in the semi-infinite entrance (exit) waveguides. The Fourier transform is calculated as

$$\begin{equation} \widetilde\psi(\mathbf k)=\mathcal F[\psi(\mathbf r)] = \int_{\mathcal{A}}\mathrm{d}\mathbf r \;\psi(\mathbf r)\, \mathrm{e}^{\mathrm{i}\,\mathbf{k\cdot r}}, \end{equation} \tag{ 18 }$$

where we extend the integral over a finite area $\mathcal {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 the x-direction with wavenumber k_n,

$$\begin{equation} \psi_n(\mathbf r) = \mathrm{e}^{\mathrm{i}k_{xn}x}\chi_n(y), \end{equation} \tag{ 19 }$$

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

$$\begin{eqnarray} \widetilde\psi_n^{\mathrm{T}}(\mathbf k)&=&\sum_m\mathcal F[t_{mn}\cdot\psi_m(\mathbf r)], \end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray} \widetilde\psi_n^{\mathrm{R}}(\mathbf k)&=&\sum_m\mathcal F[r_{mn}\cdot\psi_m(\mathbf r)]. \end{eqnarray} \tag{ 20b }$$

For a hexagonal lattice, the real and the reciprocal lattice are rotated by 90° 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 the white hexagon in figure 7(a)). To better visualize the enhancement of $\widetilde \psi (\mathbf k)$ near the K or K' points, we integrate $\widetilde \psi (\mathbf k)$ over the transverse direction

$$\begin{equation} \left|\widetilde{\overline{\psi}}(k_x)\right|^2 = \int \mathrm{d}k_y |\widetilde\psi(\mathbf k)|^2. \end{equation} \tag{ 21 }$$

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 (21) features peaks at the k_x values corresponding to K' (see figure 7(d)). The close-up of the peak in $\widetilde \psi (\mathbf k)$ near the K' point (inset of figure 7(a)) is structureless for the incoming Bloch wave with a 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 figures 7(b) and (c)) which are induced by K–K' scattering at rough edges. The relative strength of the integrated K and K' peaks (see figures 7(e) and (f)) is a direct measure of 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 the inset in figures 7(b) and (c)). The surface of the 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.

Zoom In Zoom Out Reset image size

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 micrometres with an exponential decay spanning ten 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 modelling 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 this work. Indeed, the measurement of size quantization plateaus has been surprisingly elusive in graphene nanoribbons [58–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 scanning tunnelling microscopy (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 with 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 accommodate magnetic fields while retaining its favourable scaling properties. Investigations in this direction are currently under way.

We thank K Anson, L Chizhova, J Güttinger, and C Stampfer for valuable discussions. Support from the Austrian Science Foundation (grant no. FWF-P17359), the Max Kade Foundation and the SFB 041-ViCoM (FWF) is gratefully acknowledged. SR acknowledges support from the Vienna Science and Technology Fund (WWTF) through project no. MA09-030 and the Austrian Science Fund (FWF) through project no. P14 in the SFB IR-ON. Numerical calculations were performed on the Vienna Scientific Cluster.