Range and correlation effects in edge disordered graphene nanoribbons

In this paper, we investigate the impact of edge disorder on the transport properties of graphene nanoribbons with zigzag and armchair symmetries. The diffusive and localization conduction regimes are analysed by performing a mesoscopic study on long disordered ribbons and by extracting elastic mean free paths and localization lengths. At fixed defect density and depending on specific edge disorder profile and ribbon symmetry, we observe strong transport fluctuations resulting in large mobility gaps or robust quasi-ballistic transport. These features are shown to be connected with the topology of edge irregularities as well as their correlation degree. Zigzag nanoribbons are also shown to be more robust than armchairs for similar disorder parameters.


Introduction
In the last five years, graphene has attracted much attention from the scientific community [1] for its interesting fundamental properties, most of them resulting from the massless Dirac fermion nature of the electrons in such a material, and for its potential applications, including the development of highly integrated microelectronics.
Several different techniques [2]- [7] have enabled the fabrication of graphene nanoribbons (GNRs) with width varying from several tens of nanometres down to 10 nm or less. On the one hand, GNRs show some similarities to carbon nanotubes due to their one-dimensional geometry [8], but on the other hand they also exhibit unique properties due to the presence of edges. Two main edge terminations for GNRs have been predicted theoretically [9]- [11]: a zigzag one and an armchair one. The corresponding ribbons are indicated by zGNRs and aGNRs, respectively. Experimental observations confirmed the existence of these highly symmetric GNRs and estimated a conduction gap that scales as E g ∼ α/W , with α = 0.2-1.5 eV nm −1 , where W stands for the ribbon width. There is, however, a considerable discrepancy between this experimental estimation and the theoretical band gap deduced from quantum confinement effects, which is much smaller (especially for wide ribbons). Additionally, a large number of resonances are found inside the gap, thus suggesting that the conduction gap and the band gap have two different origins and, in particular, the first one must be a consequence of disorder, which results in a strongly fluctuating tunnelling transport. In addition to usual disorder sources, especially charged ions in the substrate above which graphene is deposited, GNRs are subject to edge disorder due to the high reactivity of the edges, which are exposed to chemical passivation, roughness and topological reconstruction. While passivation can be reduced or controlled by annealing processes and exposure to chosen molecules, the geometrical quality of the edges is much more difficult to bring under control, even in the case of the ultrasmooth GNRs obtained by the most recent and innovative techniques [7]. Experimental characterizations [12]- [19] have provided direct evidence of a multitude of edge defects with different topologies. Among them, the Klein defect [14,15] consists of an extra atom linked to a single carbon of the edge. Some indirect lines of evidence of edge disorder come from the measurement of the transport gap [20]- [24], whose origin and features have been justified in terms of edge-disorder-induced localization effects and whose large width has been theoretically predicted to some extent [25]- [32]. From both a theoretical and an experimental viewpoint, the presence of a completely random disorder might conceal the relation between transport properties and the different edge disorder topologies [33]. In this perspective, we here consider different defect geometries separately and investigate their specific impact on the transport properties of aGNRs and zGNRs. It turns out that, for a chosen defect density and termination, even a slight difference in the disorder profile can drive the system from a quasi-ballistic to a localized regime.
The paper is organized as follows. In section 2, we present the adopted tight-binding model and briefly recall Green's function method for the calculation of transport properties. This is followed by some explanations about the performed scaling analysis that is employed to calculate mean free paths and localization lengths. Section 3 reports the transport properties of aGNRs and zGNRs with completely random disorder (with two degrees of randomness) on the edges. The different behaviours of aGNRs and zGNRs are evidenced and analysed in relation to the existing literature, which presents some controversial aspects [34]. In section 4, we show the results for a zGNR with edge defects with different topologies and varying density.

3
The scaling properties of the results are then analysed and connected with the disorder features. The conclusion is given in section 5.

Methodology and scaling analysis criteria
In order to investigate the transport properties of disordered GNRs, we consider the standard configuration that consists of the disordered system connected to two semi-infinite pristine GNR-based leads that serve as electron reservoirs. In this paper, we consider quasi-equilibrium transport, i.e. assuming low bias voltage between the two external reservoirs.
Pristine GNRs are described by a tight-binding Hamiltonian with one orbital per site, corresponding to p z electrons where c † i and c i are the creation and annihilation operators for an electron on the ith site of the ribbon honeycomb lattice, the sum is performed over first-nearest-neighbour sites i, j and we have chosen zero on-site energy and t = −2.7 eV as the π-π hopping parameter [1]. In the case of aGNRs, the hopping parameter between orbitals on the edge has been scaled by a factor of 1.12 in order to account for the shrinkage of the bonding length at the edges [35]. This correction entails a slight modification of the band structure and, in particular, the presence of a band gap for aGNRs of any width. In the topological model of edge disorder, we account for removed carbons by setting to zero the hopping parameters between the corresponding orbitals and their neighbours. The specific disorder-building procedure varies according to the different defect topologies and it will be explained in the following sections.
The conductance G of the system is evaluated by means of the commonly used Green's function approach adapted to the underlying honeycomb geometry, by using the Hamiltonian defined in (1) [36]. The result is the standard Landauer formula where T is the transmission coefficient, G R and G A are the retarded and advanced Green's functions evaluated on the disordered region and (left) and (right) are the rate operators corresponding to the left and right pristine leads. Thanks to the decimation-renormalization technique, the rank of all the matrices involved in the calculation can be limited to the number N of zigzag chains or dimer lines that compose the ribbon. In our calculations, we consider a large number of different disorder configurations and obtain an average conductance G = (2e 2 / h) × T , which can show the typical features of different transport regimes. In particular, the diffusive regime is characterized by where N ⊥ is the number of active conduction channels, L is the length of the disordered region, e is the mean free path and R is the resistance of the system without contacts. According to [37]- [39], we consider that the system is in the diffusive regime when the statistical properties of T are well defined, i.e. where T is the standard deviation of the transmission T . In the localized regime, transport takes place by tunnelling through localized states and T fluctuates considerably between very small values and values close to 1. In this case, a well-defined statistical quantity is provided by ln T with ln T ∝ −L/ξ, where ξ is the localization length. The average conductance is expected to decrease as T ∝ exp[−L/(2ξ )] for sufficiently long system length. The localized regime is well established when ln T ln T < 1 and T T > 1.
By means of these discrimination criteria and by analysing the length scaling of the conductance or its logarithm, e and ξ are numerically extracted. Another parameter useful for our analysis is the reference conductance G * = 10 −2 × 2e 2 / h below which the conductance is not considered experimentally significant [30]. This allows us to identify the mobility edge, to discriminate between regions inside and outside it and, in certain cases, to estimate the width of the transport gap.

Completely random disorder
Let us consider GNRs with a width of about 3.2 nm, which correspond to aGNRs made up of 27 dimer lines (denoted by 27-aGNRs) and zGNRs made up of 16 zigzag chains (denoted by . In the completely random disorder model, each carbon atom on the edges has a certain probability (which we fix as 7.5%) of being removed without any constraints, as shown in the top regions of figure 1. This type of disorder is classified as A1 and Z1 for aGNRs and zGNRs, respectively. The presence of single pending atoms, a kind of Klein defect, makes these structures highly irregular and exposes the edges to topological reconstruction. We can regularize the geometry of the edges by forbidding the presence of pending atoms, as reported in the lower region of figure 1. We indicate this new disorder profiles as A2 and Z2.
To understand conductance scaling in such systems, an analysis of transport length scales is performed, as explained in the previous section. Figure 2 shows the conductance for the four types of edge disorder profiles, as a function of the disordered region length L and the energy of the injected electrons E. The conductance is averaged over 1000 different random configurations. The case of A1 (figure 2(a)) and A2 (figure 2(b)) clearly shows the different impacts of the two topologically different disorder profiles on the conductance, as is evident by looking at the black lines that identify the (L,E) regions where G = G * . For A2, the first and second plateaus remain fairly conductive up to a length L ∼ 200 nm. In contrast, the case A1 shows a low conductance G < G * for very short lengths, although in both cases the density of defects is the same. A similar behaviour is observed in the case of the 16-zGNRs with disorder type Z1 and Z2, reported in figures 2(c) and (d), respectively. This is particularly marked within the first conductance plateau, where the onset of the low-conductance region takes place at lengths much longer for the case Z2 than for the case Z1. From such a comparison, it is clear that the specific local defect topology is crucial to tune the conduction ability of a long disordered ribbon.
Moreover, even for exactly the same types of defect and defect density, zGNRs and aGNRs show large differences. Specifically, the region of the first plateau is highly robust to randomness in the case of zigzag symmetry, as observed earlier for other types of model disorders [25]- [27], [40]. Note that other results in the literature [29] report an equivalent behaviour of aGNRs and zGNRs when the disorder erodes more than one dimer line or zigzag chain at the edges.
The procedure used for discriminating between different transport regimes and to extract the transport length scales is exemplified in figures 3(a)-(d) in the case of Z2 for an energy  E ≈ 250 meV. Figure 3(a) shows the typical decay of T versus disordered ribbon length L. The ratios T / T and ln T / ln T are reported in figure 3(b), where the vertical line marks the border between the diffusive regime to the left and the localized regime to the right, according to the criteria of (4) and (6). Figure 3(c) reports the resistance R and the corresponding fitting used to evaluate e . In figure 3(d), we show both ln T and ln T . From the fit of the first quantity, we obtain the localization length ξ , whereas, according to scaling theory, the slope of ln T should give twice ξ . As is evident from the two curves, this relation is more or less satisfied even though some deviations (up to 20%) are observed at different energies.
With the same procedure as above, we analyse the scaling properties in the cases A1, A2, Z1 and Z2 throughout the whole energy spectrum. The estimated e and ξ are reported in figure 3(e) for A1 and A2 and in figure 3(f) for Z1 and Z2.
Let us start by analysing the result for Z1 and Z2 reported in figure 3(f). As already observed from figure 2, the typical lengths are quite large in the region corresponding to the first plateau with respect to the other energy regions. This can be explained by considering that, in this range of energy, the current mainly flows in the bulk of the ribbons and thus it is not much affected by the defects located at the edges. This is also true in the regions close to E = 0, where the states are located at the edges but the current still flows through the bulk [41] due to the absence of intersublattice hopping in the first nearest-neighbour tightbinding model adopted here. A more sophisticated second (or third) nearest-neighbour tightbinding model might lead to a different result in this region around zero energy [31]. Moreover, we note that e ≈ ξ in the first plateau. This is a consequence of the more general behaviour that is clearly visible in figure 4(a), where we reported the ratio ξ/ e for the whole energy spectrum. According to the Thouless relation, the ratio between the localization length and the mean free path is proportional to the number of conducting channels (N ⊥ ). The different scaling lengths in the region of the first conductance plateau highlight the strikingly different impacts of Z1 and Z2 disorder types on the transport property in this energy range. In the higher subbands, Z 1 e ≈ Z 2 e and ξ Z 1 ≈ ξ Z 2 , thus suggesting a substantial equivalence of Z1 and Z2 at higher energies.
The situation for A1 and A2 disorder topologies is reported in figures 3(f) and 4(b). The different scales of figures 3(e) and (f) confirm that, at least for the considered disorder profiles, aGNRs are less robust than zGNRs. Again, the features of the evaluated scaling lengths evidence 8 a different impact of the A1 and A2 on the conductance. In contrast to zGNRs, these differences go beyond the low-energy region and give another indication of the higher robustness of zigzag symmetries that are less sensitive to small changes in the allowed defect zoology. Figure 4(b) shows how the ratio ξ/ e scales with the number of available conduction channels N ⊥ . As for zGNRs, a linear scaling is found in good agreement with prior results [25]- [27].

Scaling analysis of a specific disorder topology
Driven by the results of the previous section, we now focus on the case of 16-zGNRs where specific disorder types are investigated separately. The most irregular of the studied topologies corresponds to the Klein defect (K). Its geometry is shown in figure 5(a). In addition, there are other types of defects denoted by H1, H2 and H3 disorder types and they consist of 1, 2 and 3 consecutive missing hexagons with a separation of at least 1, 2 or 3 hexagons between adjacent defects, as shown in figures 5(b)-(d). The spatial regularity of the considered edge defects decreases from case K to case H3. Defects are randomly distributed along the edges of the ribbons with a density P with respect to the number of carbon atoms on each edge. For fixed P and length L, the number of defects is the same for the different topologies, while the number of atoms removed is clearly larger for more regular defects. Figure 5 provides the average conductance as a function of L and E for K, H1, H2 and H3 with P = 2.5, 5 and 7.5%. For each case, a blue (respectively yellow) line is superimposed on the 3D plot to pinpoint the frontier between the diffusive and the localization regime (respectively the border between regions where G G * and G > G * ). The difference between cases K and H3 is particularly striking. While the K profile induces a localization regime even for short lengths L, and consequently large transport gaps, the more regular H3 profile has a much weaker effect on the conductance with a robustness of the quasi-ballistic regime even for high defect densities. The associated mean free paths for K, H1, H2 and H3 profiles at fixed P are given in figure 6. The logarithmic scale emphasizes how e varies by orders of magnitude when changing the spatial regularity of the defects. As already pointed out, this effect is more significant in the first plateau and vanishes, or even reverses its tendency, when many conduction channels are active. The only exception to this behaviour is given by the profiles K and H1, which have very similar mean free paths that even coincide in a narrow region around E = 0. Another peculiarity is observed for the very regular profile H3 at high density. In fact, while a higher P generally corresponds to a shorter mean free path, from figure 6(c) it is obvious that e shows an opposite trend in the case of H3, at least in some regions of the spectrum.
We can justify all these properties in terms of the topological features of the considered disorder. The most important factor is the spatial regularity of the defects. The larger the defect width, the more the electronic structure of the disordered zone resembles that of a pristine ribbon, with, for example, preserved zigzag edge states. This phenomenon reduces backscattering with a consequent enhancement of the mean free path from the K to the H3 defect profile. Another point to be taken into account is the spatial correlation between defects. When increasing the density P of regular defects (such as H3 or more missing hexagon-defects), the geometrical constraints entail the formation of highly correlated structures, similar to a superlattice geometry. This gives rise to the peculiar behaviour observed in figure 6(c) and, again, contributes to reduced backscattering. The K and H1 defect types present some analogies that explain their similar behaviour, especially close to E = 0. First of all, both of them involve single carbons on the edges. In fact, single atoms have been added in K and single atoms have been removed in H1. However, when a few conductive channels are active, e is usually slightly longer in the case of H1; this means that the impact of K on the conductance is stronger. To understand this feature, one can consider that, in the calculation of Green's functions, the extra atoms in K act on the connected atoms P=7.5% P=5% P =2.5% Figure 6. Comparison of the mean free path e for a 16-zGNR with different disorder topologies (K, H1, H2 and H3) at chosen defect density P = 2.5% (a), P = 5% (b) and P = 7.5% (c). The vertical lines indicate the activation of new conduction channels. by means of the self-energy = t 2 /(E + η), where η → 0 + . As a consequence, Klein defects are seen by the system as energy-dependent on-site disorder randomly located along the edges. This resembles the Anderson disorder [40], which has a strong impact on the conductance of zGNRs, especially when only one conduction channel is active, due to the edge nature of the states. When E is very close to zero, diverges, thus blocking the access of the electrons to the site where the Klein defect is attached. As a consequence, in this energy region, K and H1 are completely equivalent and, as observed, e coincides in the two cases.
A scaling analysis of e as a function of P is reported in figure 7. In this figure, we show the ratio between the mean free paths 2.5%,5% e evaluated for P = 2.5 and 5%, and the mean free path 7.5% e evaluated for P = 7.5%, for each disorder profile. When the number of active channels is high, e ∼ P −1 as expected for a uniform disorder distribution. This is shown by the dashed horizontal lines in figure 7, which correspond to 7.5/5 = 1.5 and 7.5/2.5 = 3. The strongest deviations from this scaling law are observed in the region where only one channel is active, especially in the case of K and H1. In this energy window, in fact, the spatial distribution of the eigenstates is highly inhomogeneous and varies considerably with the energy. As a consequence, the strength of the impact of a 'short range' disorder (such as K or H1) on the conductance is different at different energies. The more regular H2 disorder profile better fulfils the scaling law, especially for higher P, when the disorder is distributed more homogeneously. The H3 profile shows some anomalies and deviations mainly due to the presence of highly correlated structures.

Conclusions
In conclusion, we have presented a mesoscopic transport study in GNRs with different symmetries and edge disorder profiles. Large fluctuations in transport regimes have been found depending on the specific local nature of edge defects, as well as on the degree of correlation between their distributions. Spatial correlations occur when constraints are introduced in the defect topology, which further reduces the backscattering efficiency. In contrast to the case of armchair ribbons, the robustness of quasi-ballistic transport has been found to be particularly strong for zigzag ribbons and for energy within the first conductance plateau. These results give evidence of the richness of transport behaviours in such fascinating low dimensional structures.