**2. Model and Method**

Our transport setup (cf. Figure 1) is described by the Hamiltonian *H*ˆ = *H*ˆlead + *H*ˆlead-GNR + *H*ˆ GNR where the lead environment is given by

**Figure 1.** Transport setup for an armchair graphene nanoribbon where the left-most carbon atoms are connected to the *α* = *L* lead and the right-most carbon atoms are connected to the *α* = *R* lead. 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 *c*ˆ (†) 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

$$
\hat{H}\_{\text{lead-GNR}} = \sum\_{mka} (T\_{mka}\hat{c}\_m^\dagger \hat{c}\_{ka} + \text{h.c.}),
\tag{2}
$$

where h.c. stands for the hermitian conjugate. We wish to emphasise that the coupling matrix elements *Tmkα* are kept constant for all times (i.e., we work within the partition-free approach [41–43]). In practice, we choose the lead bandwidth, given by the energy dispersions  *k<sup>α</sup>*, 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 Tmkα<sup>δ</sup>*(*<sup>ω</sup>* −  *<sup>k</sup>α*)*Tkα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–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 picture

$$
\hat{H}\_{\text{GNNR}} = \sum\_{mn} h\_{mn} \mathfrak{E}\_{m}^{\dagger} \mathfrak{E}\_{n\prime} \tag{3}
$$

where the intramolecular hopping parameters *hmn* are nonzero for the nearest neighbours only, and set by the typically used carbon–carbon hopping integral in graphene *hmn* = *tC* = 2.7 eV [48–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–64] to compute the two-time current correlation function *<sup>C</sup>αβ* (*<sup>t</sup>*1, *<sup>t</sup>*2) ≡ <sup>Δ</sup><sup>ˆ</sup>*Iα* (*<sup>t</sup>*1) <sup>Δ</sup><sup>ˆ</sup>*Iβ* (*<sup>t</sup>*2) between the current deviation operators <sup>Δ</sup><sup>ˆ</sup>*Iα* (*<sup>t</sup>*1) ≡ ˆ*Iα* (*<sup>t</sup>*1) − ˆ*Iα* (*<sup>t</sup>*1) between different leads labelled by *α* and *β*. The current operator of lead *α* is related to the particle number there via ˆ*Iα* (*t*) ≡ *qd N*ˆ *<sup>α</sup>*/*dt*, where *q* is the charge of the particle. When there is no variation of current in one of the leads, i.e., <sup>Δ</sup><sup>ˆ</sup>*Iα* (*<sup>t</sup>*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 = *VL* (*t*) = − *VR* (*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 *C*(×) (*t* + *τ*, *t*) ≡ [*CLR* (*t* + *τ*, *t*) + *CRL* (*t* + *τ*, *t*)] /2, as a function of the relative time *τ*. The traversal time *τ*tr is therefore defined by the following relation:

$$\max \left| \text{Re} \left[ \mathbf{C}^{(\times)} \left( t + \tau, t \right) \right] \right| \equiv \left| \text{Re} \left[ \mathbf{C}^{(\times)} \left( t \pm \tau\_{\text{tr}}, t \right) \right] \right| \,. \tag{4}$$

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 *IL*(*t*) and *IR*(*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* (*<sup>t</sup>*1, *<sup>t</sup>*2) (projected onto the GNR subspace) [34]

$$\begin{aligned} \mathbb{C}^{(\times)}\left(t\_1, t\_2\right) &= 2q^2 \sum\_{\mathfrak{a}' \mathfrak{a}'} \text{Tr}\_{\text{GNR}} \left\{ \Gamma\_{\mathfrak{a}} G^> \left(t\_1, t\_2\right) \Gamma\_{\mathfrak{a}'} G^< \left(t\_2, t\_1\right) \\ &+ \quad \text{i} \mathbf{G}^> \left(t\_1, t\_2\right) \left[\Lambda\_{\mathfrak{a}'}^+ \left(t\_2, t\_1\right) \Gamma\_{\mathfrak{a}'} + \Gamma\_{\mathfrak{a}} \left(\Lambda\_{\mathfrak{a}'}^+\right)^\dagger \left(t\_1, t\_2\right)\right] \\ &+ \quad \text{i} \left[\Lambda\_{\mathfrak{a}}^- \left(t\_1, t\_2\right) \Gamma\_{\mathfrak{a}'} + \Gamma\_{\mathfrak{a}} \left(\Lambda\_{\mathfrak{a}'}^-\right)^\dagger \left(t\_2, t\_1\right)\right] \mathbf{G}^< \left(t\_2, t\_1\right) \\ &- \quad \text{A}\_{\mathfrak{a}}^+ \left(t\_2, t\_1\right) \Lambda\_{\mathfrak{a}'}^- \left(t\_1, t\_2\right) - \left(\Lambda\_{\mathfrak{a}'}^+\right)^\dagger \left(t\_1, t\_2\right) \left(\Lambda\_{\mathfrak{a}}^-\right)^\dagger \left(t\_2, t\_1\right) \right), \end{aligned} \tag{5}$$

where the *Λ*<sup>±</sup> *α* (*<sup>t</sup>*1, *<sup>t</sup>*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:

$$\lim\_{t \to -\infty} \mathcal{F} \left\{ \text{Re} \left[ \mathbf{C}^{(\times)} \left( t + \tau, t \right) \right]; \omega \right\} = \frac{P\_{\rm LR} \left( \omega \right) + P\_{\rm RL} \left( \omega \right)}{2}. \tag{6}$$

Here, *<sup>P</sup>αβ* (*ω*) is defined as the Fourier transform with respect to *τ* of the real part of *<sup>C</sup>αβ* (*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 (*<sup>t</sup>*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.
