Localisation and finite-size effects in graphene flakes

We show that electron states in disordered graphene, with an onsite potential that induces inter-valley scattering, are localised for all energies at disorder as small as 1/6 of the band width of clean graphene. We clarify that, in order for this Anderson-type localisation to be manifested, graphene flakes of size of approximately 200 x 200 nm^2 or larger are needed. For smaller samples, due to the surprisingly large extent of the electronic wave functions, a regime of apparently extended (or even critical) states is identified. Our results complement earlier studies of macroscopically large samples and can explain the divergence of results for finite-size graphene flakes.


I. INTRODUCTION
In two-dimensional (2D) quantum systems, uncorrelated potential disorder has been shown to lead to complete (Anderson) localisation of single particle states [1][2][3][4] . This statement has been supported by a wealth of experimental, numerical and theoretical results, including the celebrated scaling hypothesis 5 and seminal works based on the nonlinear σ model 4,6 . States in a 2D system are marginally localised even for small disorder and d = 2 is the lower critical dimension of the Anderson transition in time-reversal invariant systems. However, while this statement is true in general, it has also been shown that the situation is more complex when correlations in the disorder 7,8 or many-body interactions 6,9 have to be taken into account. Even without these additional factors, the 2D situation remains challenging since the extent of the localised states for weak disorder can become much larger than the available system sizes, which might lead to results of a feigned extended behaviour.
In graphene, as prototypical 2D material 10,11 , one naturally expects disorder to lead to localisation as well. However, due to its linear dispersion relation around the Dirac point at energy E = 0 and non-zero momentum, the resulting absence of backscattering in clean samples 12 , might lead to a somewhat unusual behaviour. The localisation properties of graphene in the vicinity of the Dirac point have been studied intensively. It was found that strong disorder leads to localisation at E = 0 13,14 , while disorder that does not lead to intervalley mixing does not 15,16 . The direction of transport along graphene 17 and graphene nanoribbons 18,19 was shown to modify the quantitative strengths of the localisation effects. On the other hand, many, mainly numerical, results have indicated the existence of localisation that is unusually weak at E = 0 20-24 or close to E = 0 25 . Some results supporting mobility edges 20,22 , critical states 21,23,24 and a metallic-like to insulating transition 25 have been put forward. Recent discussions of results at E = 0 26,27 or for strong disorder at E 0 28 indicate complete localisation for disorder with inter-valley mixing, in agreement with the earlier studies 14-16 and a true metal-insulator transition has only been observed in hydrogenated graphene 29,30 .
Nevertheless, these studies still leave the regime of small energies that are close to but away from E = 0, for weak but inter-valley mixing onsite disorder unresolved, where Ref. 25 (see their fig. 2) found evidence for a transition-like behaviour. In fig. 1 we show this behaviour for 2D graphene flakes with 700 2 lattice sites. Clearly, increasing the size M 2 of the graphene samples leads to increasing localisation lengths around E ≈ 0.25, with energy in units of the hopping energy between carbon atoms, while around E = 0.9 the trend seems to have reversed. In this paper, we will show that fig. 1 does not indicate the existence of a transition to delocalised states. Rather, we find that the finite-size trend reverses towards localised behaviour upon further increasing the system size. However, we will need to go to very large system sizes of the order of 2.25 × 10 6 lattice sites to show this. For smaller system sizes from about 360, 000 to about 10 6 20-24 , scaling results indicate roughly a system size independence of Λ M . Hence our results explain why there is such a diversity of results for the localisation properties of graphene at and close to E = 0, i.e. we find that one needs very large system sizes, larger than 2 × 10 6 lattice sites, to reach the asymptotic regime.

II. NUMERICAL APPROACH
Our calculation is based on the standard 2D singleparticle Hamiltonian The off-diagonal elements model the hopping in transverse direction while t l ≡ t C l is the hopping along the l direction with C l denoting the connectivity matrix between layers l and l + 1 13,28,32 . All energies are measured in units of the hopping energy, t. The electronic problem defined by the Schrödinger equation Hψ = Eψ for the Hamiltonian (1) can be studied conveniently by the transfer-matrix method (TMM) 3,28 . However, since we are not interested in the quasi-1D problem of graphene nanoribbons with L M 18,19 , we need to modify the TMM to allow the treatment of 2D M ×M graphene samples. 33 This has implications for the convergence of standard TMM calculations since we can no longer use the self-averaging property normally used for L → ∞. Our modification involves the definition of forward and backward transfer matrix multiplications 34,35 . The method also yields the inverse localisation length 1/λ M (E, W ), but only for a single M × M graphene sample. Afterwards, the 1/λ values need to be averaged for many M × M disorder configurations with the same parameters M , E and W .
The TMM must be adapted to handle the hexagonal structure of the graphene lattice 13,32 by suitably chosen l and C l matrices. We distinguish between transport directions parallel to armchair (AC) and zig-zag (ZZ) edges. Our approach is similar to Ref. 28 and for more details see Ref. 36. A pictorial representation is shown in the inset of fig. 1 for AC and ZZ graphene. 37 We chose hard wall boundary conditions for all results presented here. In order to have the same number of atoms for both ZZ and AC edges, the width of the AC sample should be chosen as M AC = L ZZ /2 and the length as L AC = 2M ZZ . In this way we ensure that we are studying the same sample but in both directions of transport.
The scaling hypothesis for finite-sized systems implies The λ M data can be rescaled numerically by a least-squares fitting procedure to obtain the scaling function f 3,38 . In the case of the 2D Anderson model on a square lattice, this function has a single finitesize scaling (FSS) branch with decreasing Λ M for increasing M -indicating the localised regime. For the 3D Anderson model, the same procedure leads to two branches, the first one denoting the localised regime and the second one indicating the extended regime with increasing Λ M values as M increases. This two-branch behaviour is the signature of the transition from localised to extended states 3 . Alternatively, we can try to assume an analytical form for f and test whether this form fits the data with the required accuracy 39,40 . Assuming e.g. the power-law behaviour f ∝ |1 − E/E c |L 1/v of the 3D Anderson transition, then this approach allows not only to construct f , but also determines the critical exponent ν and the energy E c (or disorder W c ) at which the transition occurs 39,41 . We will use both FSS approaches below.

III. RESULTS
In fig. 2 we show the variation of the disorder-averaged localisation length λ M (E) for different values of disorder W . The lattices correspond to square lattice, AC and ZZ graphene. In each case, the system sizes where chosen such that M × L = 10 4 lattice sites, corresponding to M = 100 and L = 100 for the square lattice and the ZZ graphene, but M = 50 and L = 200 for the AC graphene lattice. We first note that at weak disorder (W ∼ 1) the half of the bandwidth reflects the number of nearest neighbours and hence tends to 4 for the square lattice and tends to 3 for AC and ZZ graphene 42 . Furthermore, there is the usual approximate symmetry between positive and negative energies. When the strength of the disorder increases the λ values decrease for all lattices as the wavefunctions become more localised. For very strong disorder, the localisation lengths are much smaller than the system sizes M and L and the states are exponentially localised with λ representing the decay length. On the other hand, for weaker disorder, the localisation lengths are comparable or larger than the system sizes, and we can no longer assume that the exponential decay implicit in the use of λ is still justified. Then λ is simply a convenient measure of the spatial extent of the wave functions, but not necessarily linearly related to a localisation length. Still, a larger such extend will imply larger λ values. With this in mind, we see in fig. 2 that, for W 4, the localisation lengths increase rapidly as we decrease W for the square lattice. However, for the case of AC and ZZ graphene lattices, we observe that in the vicinity of E = 0, the λ values again decrease, leading to values of λ M (E ≈ 0) which seem very similar for W = 1 and 2. Clearly, the drop in λ M in the graphene lattices at E = 0 is a signature of the Dirac point with reduced density of states 43,44 .
In standard quasi-1D TMM, an increasing value of Λ M for weak disorder as M → ∞ signals the start of the extended regime. Even with Λ M > 1, λ M can still be interpreted as a localisation length since we have L M and the localisation in the l direction is well-defined. As discussed before, the situation might be different for our modified TMM. Nevertheless, we already see from fig. 2 that for energies |E| 1, the λ values for the square lattice and AC/ZZ graphene behave similarly. If any new, graphene-specific, finite-size behaviour can be expected, it should be around E ≈ 0. Therefore we have studied in fig. 1 the finite size behaviour of Λ M in ZZ graphene for energies 0 ≤ E ≤ 1 at weak disorder W = 1.5 when Λ M ≥ 1. As one can see from this figure, for energies larger then E = 0.9, increasing M (and L) leads to a decrease of Λ M , the traditional signature of localisation. However, for energies E 0.6, increasing M gives increasing Λ M values. Such a behaviour for M → ∞ would indicate extended states. Quite similar findings have been reported previously in the same energy range for smaller systems up to M = 252 25 .
Clearly, the existence of extended states in the vicinity of the Dirac point in weakly disordered (but with inter-valley mixing) graphene would be surprising. However, let us already note several suspicious observations, namely (i) there is no clear crossing point, rather a series of not well-defined crossing points in the region E ∈ [0.7, 0.9]. Furthermore, (ii) increasing the system size does not lead to a clearer crossing, and we can also not identify a simple, monotonic in M (irrelevant) shift of such a crossing point. Let us also emphasize that system widths of M = 700 as used in fig. 2 are already reasonably large for TMM 39 . If there truly was a metal-insulator transition in the indicated energy range, then we would expect to see good quality FSS. On the other hand, if the behaviour of fig. 1 was simply due to finitesize effects, then we should see the increase in Λ M vanish for large enough M . Since the increase seems largest at energy E = 0.25, we shall study this energy in detail for a square lattice as well as AC/ZZ graphene.
In fig. 3 we show FSS results for Λ M in square lattices and ZZ graphene with M and L values chosen such that the number of sites M × L ranges from from 100 2 to 700 2 . For strong disorder, we have Λ M ∝ 1/M as expected since states are highly localised and λ M is constant for M λ M as indicated. Decreasing the disorder -or, equivalently, decreasing M -leads to deviations from the simple 1/M behaviour and indicates that ξ(W ) starts to increase. In the standard quasi-1D square lattice TMM, this leads to an evermore flat behaviour for Λ M (W ) as W → 0. We indeed observe this behavior for E = 0 in square lattices, AC and ZZ graphene (not shown) 36 . For smaller disorder, W 2, we find the reconstruction of a well-defined FSS curve becomes numerically difficult. Nevertheless, the estimated scaling parameter ξ(W ) agrees very well with a previous highprecision FSS from a quasi-1D TMM 45 . Furthermore, the ξ(W ) behavior for squares and ZZ graphene shows a single branch only, consistent with complete localisation.
The situation is rather different for E = 0.25 as shown in fig. 3. We see that FSS gives rise to localised branches as well as the beginnings of what look like extended branches. Here it is intriguing to see that even for a square lattice, for the range of available system sizes and disorder -determined by the longest TMM runs available to us -we find an apparent transition-like behavior. Obviously, this would be in disagreement with the scaling theory and of course also to the body of numerical results based, among others, on quasi-1D TMM 3,4 . Similarly, we observe transition-like behavior also for ZZ graphene at E = 0.25. As in the square lattice case, the onset of the extended branch is around W 2. We have found similar results also for AC graphene. 46 We have also tried to apply FSS assuming the expansions of the power-law behaviour 39,41 . However, we never     find an acceptable fit to the data, although we vary not only the expansion coefficients, but also the initial values used in the non-linear fits for W c , ν, etc. Upon closer inspection, we find that most such attempts to fit the data lead to W c ∼ 0 and large values of ν > 5. But even with these large ν values, the Λ M values rise much faster for small disorder. This suggests that the true behaviour is not a power-law but rather an exponential as in the well-known square lattice 48 .
The FSS results of fig. 3 for E = 0.25 and W 2 do not show a very clear formation of extended branches, particularly for the square case. In order to test the stability of these branches in FSS, we would need even larger system sizes for all disorder W 2. This is, however, numerically prohibitive. 49 Thus we have chosen to restrict ourselves to two disorder strengths, W = 1 and 1.25 for E = 0.25. Even with this restriction, a considerable number of runs for M > 900 do not finish within our chosen maximum time limit of about one week. Such λ M values have therefore a relative error n , with n denoting the sample, larger than the target of 0 ≡ 5 × 10 −5 . Hence we weigh such results less when computing an average. With (i) w n = 1/ 2 n or (ii) w n = max(1, 0 / n ), we define the averaged Lyapunov exponents as γ M = n w n γ n / n w n with weighted standard-deviations n w n (γ n − γ M ) 2 / n w n . In case (ii), samples, which have converged better than the target, are given less weight in order to test the robustness of our results.
We show  A qualitative argument can be put forward to motivate our results. Without disorder the density of states (DOS) at E = 0 for a square lattice diverges whereas it is zero for graphene (ZZ or AC). Upon increasing the disorder, the DOS for the square lattice decreases as well as the localisation length. For graphene, the same happens for the van-Hove singularities at E = ±1 44 . On the other hand, at E = 0 the DOS increases 44 , which is well-known to correlate with large localisation lengths. The crossover between these two regimes should be expected around E ≈ 0.5 which is similar to what we observe. For larger disorder W 2 or, equivalently, larger system sizes M 1400, we recover the expected localised regime.
Once the modified TMM has reached convergence, the wave functions (ψ l , ψ l−1 ) are true eigenfunctions of the global 2M × 2M forward-backward transfer matrix T † L T L for a given sample. Hence ψ l,m , l, m = 1, . . . , M , is the eigenfunction of H. In fig. 5 we show the results for ZZ graphene at four different values of disorder at E = 0. For weak disorder W = 0.5 and 1, one can clearly see the enduring presence of edge states previously predicted for clean ZZ samples 50 . For stronger values of disorder, the spatial disorder distribution itself becomes dominant. At E = 0.25 there is no evidence of edge states. Results for AC graphene are similarly consistent with the literature, i.e. we find an absence of edge states for the chosen AC graphene lattice sizes consistent with semiconducting behaviour on finite width samples 50 . As expected for square lattices, we do not observe those strong edge states.

IV. CONCLUSIONS
Our results show that up to lengths scales of 1500 times the C-C distance in graphene, i.e. up to 213 nm, onsite disordered graphene, even with inter-valley scattering, exhibits surprisingly delocalised states in the vicinity of the Dirac point. This corroborates the trend towards similar such delocalisation-like behaviour found previously [20][21][22][23][24][25] , while also reaffirming that the true infinite system limit obeys the localisation predictions [14][15][16][17][18][19] . In fact, the tendency for large localisation lengths is so strong that even FSS can mislead to construct seemingly extended branches, although a very large system size analysis shows that only the localised behaviour corresponds to the true thermodynamic behaviour 15,16 . We emphasise that our results also explain graphene's robustness against defects in similarly sized ribbons 51,52 , billiards 53 and quantum dots 54 . Our approach is based on a modified TMM which allows to study "square" flakes of graphene. This TMM can convincingly reproduce the infinite-size estimates of localisation lengths obtained from standard TMM and we expect the method to be useful in other contexts as well.