Electron Traversal Times in Disordered Graphene Nanoribbons

Using the partition-free time-dependent Landauer–Büttiker formalism for transient current correlations, we study the traversal times taken for electrons to cross graphene nanoribbon (GNR) molecular junctions. We demonstrate electron traversal signatures that vary with disorder and orientation of the GNR. These findings can be related to operational frequencies of GNR-based devices and their consequent rational design.


Introduction
A fundamental property limiting the operational frequency of a molecular device is the traversal time τ tr for electronic information to cross between the nanojunction terminals [1]. For instance, in graphene, the cutoff frequency f max is related to the traversal time as f max = 1/2πτ tr [2,3]. For the molecular electrician, this raises the key question of how long it takes for electronic information to propagate across a nanosized device, as this sets a fundamental limit on the speed of the device operation. In quantum mechanics, time does not have the same status as a dynamical variable such as the energy or particle position. In fact, much debate has centred around the correct definition of the traversal time through a generic potential barrier [4][5][6], as well as the relation of this quantity to the dwell time (time spent in the molecular region) [7], the Larmor clock time (the time taken to move between scattering channels) [8,9], the group delay time (the time delay in the nonlocal propagating wave packet caused by scattering off the potential barrier) [10], or to a generic description of probability distributions via path integrals [11,12]. Crucially, all the aforementioned times are defined in terms of the transmission probability, potential and incident energy of electrons moving in a static scattering theory picture [13], so that a theory which takes strong time-dependence into account is still needed. This is crucial for the understanding of laser-stimulated tunnelling processes and related to the problem of tunnelling times in strong field ionisation experiments [14][15][16].
Graphene nanoribbon (GNR)-based molecular junctions are excellent candidates for room-temperature transistors, i.e., graphene field-effect transistors (GFETs) [2], GHz-THz frequency modulators [17], and photodetectors [18], due to their high mobility and charge carrier saturation velocities. GNRs can be engineered with band gaps that are tuneable via the nanoribbon symmetry properties and widths [19], and currently, sub-10 nm nanoribbon widths are accessible from chemical fabrication techniques [20][21][22][23][24]. Recent experimental progress shows an inverse scaling of operational frequency with the nanoribbon length [3], or with the square of the nanoribbon length [20] for ribbons whose carrier drift velocity scales inversely with ribbon length. Typically, the maximum operational frequencies of radio frequency (RF) GFETs exceed those of Silicon-based transistors with the same dimensions [3], and can in principle be achieved in the 100 GHz range [25][26][27]. The cutoff frequencies of GFETs are strongly affected by the presence of defects and flexibility in the nanoribbon, but 100 GHz flexible nanoribbons have also recently been fabricated [28].
Disorder may have a profound impact on the operation of the graphene-based devices. For example, it has been investigated that edge disorder affects the armchair-oriented GNR (AGNR) more than the zigzag-oriented ones (ZGNR) [29,30]. This is because the edge states in ZGNRs are energetically protected against impurity perturbations. This is not the case for AGNRs, in which the edge states are less dominant so that disorder has a much larger effect on the conductance [31,32]. On the other hand, disorder-induced broken chiral symmetry in terms of random bond disorder was considered in Ref. [33].
In this paper, we expand upon a recent proposal [34] to investigate traversal times for electrons moving in disordered GNRs by looking directly at the dynamics of statistical correlations between electronic signals measured in different reservoirs connected to the GNR. We demonstrate that the traversal time has a clearer signature in AGNR than in ZGNR. This is because the charge densities in AGNR structures are more delocalised than in ZGNR, where the formation of standing-edge-state charge waves leads to wave fronts with a spatial orientation lying diagonal across the plane of the nanoribbon [35][36][37]. We also show that the effect of breaking chiral symmetry (on-site disorder) has less effect on the traversal times compared with the disorder that preserves chirality (hopping disorder). (A random on-site potential breaks the sublattice symmetry and can broaden possible Landau levels, i.e., the chiral symmetry is destroyed. This does not happen with nearest neighbour hopping disorder as hopping between different sublattices is by construction the same in both directions [33,[38][39][40]).

Model and Method
Our transport setup (cf. Figure 1) is described by the HamiltonianĤ =Ĥ lead +Ĥ lead-GNR +Ĥ GNR where the lead environment is given byĤ For times t ≥ 0, a voltage bias V(t) is applied in the leads and charge carriers start to flow through the graphene junction. We consider ribbons of varying lengths (2, 5, and 10 nm), and similarly also the zigzag orientation (not shown).
The operatorsĉ ( †) are electron annihilation (creation) operators. At the initial switch-on time t 0 , the energy dispersion is shifted by the bias voltage kα → kα + V α (t). The left-most atoms are connected to the left lead (L), whereas the right-most atoms are connected to the right lead (R). The coupling between the leads (α = L, R) and the GNR has the form where h.c. stands for the hermitian conjugate. We wish to emphasise that the coupling matrix elements T mkα are kept constant for all times (i.e., we work within the partition-free approach [41][42][43]). In practice, we choose the lead bandwidth, given by the energy dispersions kα , to be much larger than the coupling energies so that we may employ the wide-band limit approximation (WBLA), in which the level width matrix Γ α,mn (ω) = 2π ∑ k T mkα δ(ω − kα )T kαn is evaluated at the Fermi energy of lead α. The WBLA is justified because we are interested in a regime where the lead-GNR coupling is weaker than the internal hopping within the GNR [44][45][46][47], as this also enables us to focus on the effect on traversal time caused by the internal ribbon structure. The GNR is modelled by a single π-orbital tight-binding picturê where the intramolecular hopping parameters h mn are nonzero for the nearest neighbours only, and set by the typically used carbon-carbon hopping integral in graphene h mn = t C = 2.7 eV [48][49][50][51]. Longer-range hoppings could be included in the model similarly but here we wish to preserve the particle-hole symmetry of the undisordered GNR. While our system is completely described above, it would also be feasible to consider hydrogen passivation at the edges of the graphene nanoribbons, as is customary in typical experimental setups [52,53] and in their ab initio modelling [36,54,55], to remove dangling bonds at the edges. We employ the recently developed time-dependent Landauer-Büttiker (TD-LB) formalism [35,[56][57][58][59][60][61][62][63][64]] to compute the two-time current correlation function C αβ (t 1 , The current operator of lead α is related to the particle number there viaÎ α (t) ≡ qdN α /dt, where q is the charge of the particle. When there is no variation of current in one of the leads, i.e., ∆Î α (t 1 ) = 0 the correlation between this signal and the current variation in the other leads is trivially zero. The partition-free approach employed here also means that the whole system, described by Equations (1)-(3), is coupled during an initial equilibration period that occurs prior to the switch-on time t 0 . This leads to a global initial temperature T and lead-independent initial chemical potential µ.
In a two-terminal junction (cf. Figure 1), the labels α and β can refer to either the left (L) or right (R) terminals, whose energies are shifted symmetrically to create a voltage drop of V (t) across the system V (t) /2 = V L (t) = −V R (t). We note here that the driving bias voltage in our setup can be of dc type (sudden quench) or ac type (time-dependent modulation) [35,59,65]. In Ref. [34], it was shown that timescales associated with electron traversal times and internal reflection processes could be seen as resonances in the real part of symmetrised cross-lead correlations , as a function of the relative time τ. The traversal time τ tr is therefore defined by the following relation: The motivation behind this definition becomes more apparent when we look at the simulation data of current cross-correlations in Section 3. It is, indeed, observed that the currents in each lead I L (t) and I R (t) become most strongly correlated at the time taken for electrons to propagate between the leads following a voltage quench, which is then attributed to electron traversal events. Within the WBLA, the cross-correlations are evaluated in terms of Keldysh components of the Green's function G (t 1 , t 2 ) (projected onto the GNR subspace) [34] where the Λ ± α (t 1 , t 2 ) matrices are defined in terms of convolutions on the real and imaginary branches on the complex time contour [34]. We investigate the dynamics of steady state (the switch-on time is taken to t 0 → −∞) cross-correlations, which are accessible experimentally [66].
We note that there is some spreading in the individual resonant peaks associated with traversal times in the correlator, so that our proposal takes into account the probabilistic nature of electron propagation in accordance with realistic proposals for traversal time distributions [13,14,67]. This is particularly relevant for our approach which enables us to study arbitrary time-dependent biases, e.g., in which the drive is stochastic in time [61]. The Fourier transform of the real part of C (×) (t + τ, t), with respect to the relative time τ, is equivalent to the frequency-dependent power spectrum associated with cross-lead correlations: Here, P αβ (ω) is defined as the Fourier transform with respect to τ of the real part of C αβ (t + τ, t) [34].
In practice, the high-frequency component of the current fluctuations can be probed by studying the infrared-to-optical frequency range of light emitted by the junction [66]. It is worth noting that Equation (5) is valid in the transient regime for arbitrary time-dependent bias voltage profiles, whereas in Equation (6) the switch-on time t 0 is taken to minus infinity (t 0 → −∞) or, equivalently, the observation time t is taken to infinity (t → ∞). While all the quantities in Equation (5) depend only on the time difference τ = t 1 − t 2 at this "long-time limit" (and the Fourier transform in Equation (6) is taken with respect to this relative time), in the transient regime, it would be required to consider separate frequencies for each time t 1 and t 2 in Equation (5). The central idea of our work is that one should quantify the traversal time for electronic information to cross the system by looking directly at the correlations in the electronic signal itself, rather than trying to build an indirect definition of operational time from the calculation of transmission probabilities. The definition of traversal time here is closely related to the definition of Miller and Pollak, which makes use of flux-flux correlation functions [68]. However, the TD-LB formalism is valid for arbitrary lead temperatures, lead-GNR hybridisation strengths, and time-dependent biases.

Results and Discussion
As we consider the WBLA, the detailed electronic structure of the leads is not important for the description of the transport properties of the GNR. We then fix the coupling strength between the GNR and the leads by the frequency-independent resonance width Γ α = t C /10 corresponding to a weak-coupling regime where the WBLA is a good approximation [44][45][46][47]. This is further justified in typical transport setups where the bandwidth of the leads is sufficiently large (e.g., gold electrodes) compared to the applied bias voltage. As we wish to preserve the charge neutrality of the GNR in equilibrium, we set the chemical potential to µ = 0. The global equilibrium temperature is set by (k B T) −1 = 100t −1 C (T = 313 K).

Response to a dc Drive
It is instructive to first study the current correlations in GNRs without disorder. Figure 2 shows the current cross-correlations of undisordered AGNRs and ZGNRs of various lengths and time-independent bias voltages. We can make many general observations from the data:

•
The signal is more clear for the AGNR than ZGNR. In the AGNR case, the propagating wavefront is coherent [35], so that there is less spread in the resonant traversal time signal than in the corresponding ZGNR case. This relates to the shape of the propagating wavefront, since in AGNR it is flat, whereas in ZGNR it has a triangular shape [35]. The back-and-forth internal reflections of the wavepackets between the electrode interfaces have a fairly regular structure in AGNR which results in a clear signal in the current cross-correlation. This means devices based on ZGNR have a less well-defined operational frequency.

•
The current cross-correlations are mostly independent of the strength of the applied voltage. The voltage may affect the shape of the curves slightly, but not the location of the main resonance. This can be related to the group velocity of electrons crossing the GNR, v k = d k /dk, which should not depend on a k-independent shift in the energy dispersion k [34].

•
Evidently, there is a roughly linear increase of the time-difference between the first maxima with increasing L, due to the time taken for the propagating electron wavefront to cross the structure. The time-difference between the first maxima τ max is related to the traversal time of information through the GNRs via Equation (4).

•
Increasing the length in the AGNR does not increase the number of resonant peaks in the cross-correlations, but in ZGNR it leads to a broader range of resonances clustered about a mean traversal time. This dependence on the orientation of the GNR then affects the spread of operational device frequencies.

•
The low-frequency regions of the Fourier transforms show resonant frequencies at ω = nΩ L where n is a positive integer and Ω L is some intrinsic frequency depending on the length of the GNR. In particular, by increasing the length of the GNR more transport channels are opened in the bias window, and therefore more peaks appear in the Fourier spectra.

•
From the full frequency ranges of the Fourier transforms (insets), we observe a high-frequency operational cut-off which is smaller for the ZGNR case than for the AGNR case. This is itself an interesting effect, as it sets a limit on switches built with these kind of nanoribbons [2,3].

The Role of Disorder
As we have now established general principles for the traversal times, we concentrate our discussion on disordered GNRs of fixed length and fixed applied (time-independent) voltage. In Figure 3, we introduce disorder into the GNR without breaking chiral symmetry. This is done by drawing a random number from a uniform distribution of width w around the average hopping matrix value t C . We exclude second and third nearest neighbour hoppings, and we consider randomness only in the hoppings, for we wish to preserve particle-hole symmetry [49,69]. We see that the disorder appears to increase the traversal time and also has the effect of decreasing the quality of the signal, so that multiple side-peaks are visible. These are caused by internal reflections induced by the disorder. In addition, the intrinsic resonant frequencies for the disordered GNRs (shown by the Fourier spectra) are red-shifted due to the hopping disorder. This finding is consistent with the idea of reduced operational frequencies for the GNR devices due to disorder. In Figure 4, we break the chiral symmetry by adding a random term to the on-site energy levels of the GNR. This is also done by drawing a random number from a uniform distribution of width f around the zero on-site energy for the pristine GNR. In Figure 4a, we look at the AGNR case, and we observe in this case the deterioration of a clear traversal time signal as the Anderson localisation increases the average dwell time in the interior of the GNR. In Figure 4b, the average traversal time is once again seen to be larger in the ZGNR case. Interestingly, there is a crossover in both GNR configurations as the disorder destroys the coherence of the propagating wave packet around f = t C /2. In contrast to the case of hopping disorder, here the operational frequencies of the GNR device (shown by the Fourier spectra) remain roughly unchanged for the on-site disorder. This finding could be related to the character of the disorder: While both types of disorder may introduce an effective tunnel barrier around the disordered GNR that the propagating electrons must overcome, the hopping disorder case corresponds to deformation of the lattice geometry, whereas the on-site disorder corresponds to a change in charge neutrality or chemical potential.

Response to an ac Drive
We finally address the full two-time character of the cross-correlation in Equation (5). Specifically, we consider the case of ac driving by introducing a monoharmonic driving term to the voltage where the static part is set by V 0 = t C /2 and the amplitude of the ac driving is A = t C /2 with the driving frequency Ω = t C /10. To reduce the computational effort, we consider only the short nanoribbons in this case (L = 2 nm). In Figure 5, we show the propagation of the full two-time cross-correlation from the initial time t 0 . AGNR data are shown in Figure 5a-c and ZNGR data are shown in Figure 5d-f. In contrast to the previous steady state results, here we show the initial transient (up to 50 fs), which includes relaxation effects. We observe that the ac driving does not change the overall picture of traversal time, i.e., the time it takes for the information to traverse through the nanoribbons can be clearly read off from the separation of peaks along the anti-diagonal. Compared to the long-time limit in Figures 2-4, the initial transient only shows some additional oscillations but the main features seen in the steady state data are still visible. The two-time correlations also show the effect of disorder; as in the dc case, the signal gets considerably disturbed for hopping disorder (cf. Figure 3) and for on-site disorder (cf. Figure 4), but in the latter case the signal destruction is less severe. In the disorder energy scale considered in Figure 5 (hopping and on-site disorder strength: w, f = t C /2), it is also observed that the first traversal event peaks are more clearly visible for the case of on-site disorder (Figure 5c,f) compared to the hopping disorder (Figure 5c,e). However, the coherence of subsequent traversal events without disorder (shown in the side-peaks of Figure 5a,d), is strongly suppressed by either kind of disorder. In calculations not shown here, we have checked that other types of ac driving (biharmonic drive, faster/slower modulation, lower/higher intensity) have no effect on the qualitative behaviour of traversal times.

Conclusions
We have employed the recently developed TD-LB formalism to compute the two-time current correlation functions for disordered GNRs. This methodology is a fast and accurate way of addressing mesoscopic quantum transport phenomena out of equilibrium as it is well-supported by the underlying nonequilibrium Green's function theory [70]. By our analysis, we confirm that the current cross-correlation is a good measure of electron traversal time. We find that the traversal time scales roughly linearly with the length of the GNRs, and that the traversal time also depends strongly on the GNR orientation. In general, the results presented in our work indicate that the traversal time can be most sensibly connected to the group velocity of noninteracting electrons passing through the graphene nanoribbon. The distance between the first maxima in the cross-correlations corresponds to the time for an electron with velocity v ≈ 1 nm/fs to cross the ribbon. This velocity is consistent with the value of the Fermi velocity in graphene v F = 3t C a/2h where a = 1.42 Å is the carbon-carbon distance [49]. This can be related to the Büttiker-Landauer time for tunnelling [4] through a rectangular barrier as the time it would take a particle with velocity v F to traverse the barrier.
We found that disorder in GNRs increases the traversal times, in general, and ultimately destroys the whole picture of coherent information transfer over the GNR junction when the disorder-induced scattering is strong. The "rule-of-thumb" character of our findings is summarised in Figure 6. We considered two types of disorder, one that preserves (hopping disorder) and one that breaks (on-site disorder) the chiral symmetry of the GNR. In Figure 6b, we find that the intrinsic operational frequency of the GNR is redshifted for the hopping disorder while, in Figure 6c, we see that it remains roughly unchanged for the on-site disorder. However, in the latter case, the statistical spread of τ tr is significantly enhanced as the on-site disorder is increased. To measure the current cross-correlation and extract experimental values for τ tr , there exist a range of spectroscopic techniques which relate the field strength of photons emitted from each lead to the current. The zero-and finite-frequency current cross-correlations can then be extracted from these current measurements [66,71]. Our noninteracting approach is sufficient for graphene structures since monolayer graphene devices have been experimentally shown to have ballistic transfer lengths in the range from 100 nm at room temperatures to 1 µm at sub-Kelvin temperatures [2,72]. Even though our approach is limited to noninteracting electrons, we expect the current correlations and traversal times to be similarly related even when dealing with, e.g., electron-electron or electron-phonon interactions [47,[73][74][75]. If perturbation theory could be applied, i.e., when the interaction is weak, current correlation or noise simulations are still feasible to perform in terms of the one-particle Green's function [76][77][78][79][80][81]. Disordered interacting systems have also been studied using mean-field or density-functional theories [82][83][84]. Here, out-of-equilibrium dynamics only due to voltage bias was considered but also thermal gradients could be included similarly [46,[85][86][87][88]. At present, for the case of strong interaction, these approaches cannot yet be extended to realistic device structures since considerably more complicated and numerically expensive methods are required [89].
We also confirmed that the overall picture of electron traversal times is not qualitatively changed by introducing an ac driving voltage compared to the response to a dc drive. Possible quantitative differences in the response signals to an ac drive could be related to signatures of photon assisted tunnelling on traversal time [37,65] but, for now, will be left for future work. On the other hand, it would also be interesting to consider, e.g., a short laser pulse for exciting the system out of equilibrium [90][91][92] instead of the quench of the voltage bias employed in the present work. These topics will also be addressed more thoroughly in a forthcoming paper.