Time-resolved impurity-invisibility in graphene nanoribbons

We investigate time-resolved charge transport through graphene nanoribbons supplemented with adsorbed impurity atoms. Depending on the location of the impurities with respect to the hexagonal carbon lattice, the transport properties of the system may become invisible to the impurity due to the symmetry properties of the binding mechanism. This motivates a chemical sensing device since dopants affecting the underlying sublattice symmetry of the pristine graphene nanoribbon introduce scattering. Using the time-dependent Landauer--B{\"u}ttiker formalism, we extend the stationary current-voltage picture to the transient regime, where we observe how the impurity invisibility takes place at sub-picosecond time scales further motivating ultrafast sensor technology. We further characterize time-dependent local charge and current profiles within the nanoribbons, and we identify rearrangements of the current pathways through the nanoribbons due to the impurities. We finally study the behavior of the transients with ac driving which provides another way of identifying the lattice-symmetry breaking caused by the impurities.


I. INTRODUCTION
Being under considerable research focus for the past two decades graphene [1] and carbon nanotubes [2] have been found to be extremely sensitive to external perturbations.For this reason, these nanomaterials have been proposed as ideal candidates for sensor technology [3][4][5].Based on these observations, carbon-based sensor devices have already been developed at the single-molecule resolution [6][7][8][9].Carbon-based transducers have been embedded in circuitries involving graphene nanopore platforms [10,11] and field-effect transistors [12], and these have been successfully applied to, e.g., disentangle biomolecules' rapid dynamics [13][14][15].This novel biosensor technology provides real-time information about the underlying physical and chemical mechanisms.These findings motivate ultrafast sensing devices as present transport measurements are able to resolve temporal information at sub-picosecond time scales [16][17][18][19][20].
These few examples include devices which are interacting with their persistently evolving environment for which a schematic in a quantum transport setup is shown in Fig. 1.The schematics depicts a quantum transport channel made of a graphene nanoribbon (GNR) of armchair configuration with a width characterized by the number of carbon-dimers arranged transversely (N ).The ribbon is subjected to a source-drain voltage (V SD ) that can vary on time and the whole device serves as a host for detecting the presence of impurities.The main purpose of the voltage bias is to excite the system away from its thermo-chemical equilibrium.We emphasize that the perturbation could be different in practice and still similar intrinsic dynamics would show up, irrespective of the specifics of the perturbation.It would also be feasible to use, e.g., short laser pulses which can perturb the system and then probe the consequent dynamics [21][22][23][24][25]. Nonetheless, the theoretical description of these processes is a challenge as these nanoscale devices are operating at high frequencies (THz) so the systems do not necessarily relax to a steady-state configuration instantly.In contrast, there are transient effects depending on, e.g., the system's geometry or topological character [26][27][28][29][30][31], its predisposition to external perturbations or thermal gradients [32][33][34][35][36][37], and the physical properties of the transported quanta and their mutual interactions [38][39][40][41][42][43][44].There is an increasing demand for theoretical and computational tools capable of addressing, in a general but computationally tractable level, the time-dependent responses of nanomaterials.
In this paper, we investigate how impurity invisibility in graphene [45,46] takes place in the transient regime.Depending on the dopant conformation with respect to the underlying graphene nanoribbon lattice symmetry, we identify whether the charge transport properties of the conducting device are modified due to impurity scattering.We observe that the time-resolved signals are highly sensitive to the impurity configurations: in the transient regime the charge-current pathways are reorganized due to the scattering introduced by the impurities.We argue that recognizing this mechanism is pivotal for efficient design of graphene-based ultrafast chemical sensing devices.

II. MODEL AND METHOD
We consider a system composed of metallic leads α connected to a central molecular structure, where we investigate the (charge) transport of noninteracting electrons.Even though, electron-electron and electron-phonon interaction should, in principle, influence the transport mechanisms, here we expect our noninteracting picture to be sufficient as recent studies on small monolayer graphene devices have revealed ballistic transfer lengths ranging from hundreds of nanometers to even micrometers at low temperatures [47,48].The transport setup (cf.Fig. 1) is partition-free [49][50][51] which means that the whole system is initially contacted in a global thermochemical equilibrium at unique chemical potential µ and at temperature (k B T ) −1 = β.The central molecular structure is modeled by a tight-binding Hamiltonian where T mn accounts for hoppings between the lattice sites m and n.The operator ĉn (ĉ † m ) annihilates (creates) an electron on site n (m) of the host lattice.In practice, we consider GNRs as the central molecular structure, and we set T mn = −γ for nearest neighbours m and n, and γ = 2.7 eV. Second and third nearest neighbour hoppings could be included similarly [52,53] but here we consider particle-hole symmetric cases and only take the first nearest neighbours into account.For impurities in the central conducting device we have similarly where the index n runs over the impurities, and a n is the host site in the (pristine) central region where the impurity is attached to.The operator ĉj (ĉ † j ) annihilates (creates) an electron on site j that can be on the host lattice or impurity site.Parameters for the on-site ( imp ) and hopping energies (γ imp ) for the impurities can be related to ab initio calculations [54,55] following Density Functional Theory (DFT).The impurity configurations considered here are shown in Fig. 2.An impurity can be positioned (i) right on the top of a carbon atom (T), (ii) positioned over a carbon-carbon bond characterizing a "bridge" configuration (B), (iii) placed on the center of an hexagonal ring (C), or (iv) substitute a carbon atom, yielding a "substitutional" configuration (S).The leads are described as (semi-infinite) reservoirs with kα corresponding to the energy dispersion for basis states k in lead α.For a one-dimensional tight-binding structure this is given by kα = 2t α cos k with t α the hopping energy between the lead's lattice sites.Here we assume that the density of states of the leads is smooth and wide enough allowing us to consider the wide-band approximation (WBA).In WBA the lead density of states is assumed independent of energy, which in practice means that we choose the energy scales in the leads much higher than other energy scales in the central region.As we concentrate on the effects between the graphene nanoribbon and the impurities this allows us to neglect the precise description of the electronic structure of the leads.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 [37,[56][57][58].The leads are connected to the central region by the coupling Hamiltonian with T m,kα the hopping energy coupling lead's states with the lattice sites of the molecular structure.The total Hamiltonian is then combined as Ĥtot = Ĥmol + Ĥimp + Ĥlead + Ĥcoupl .We consider a switch-on of a bias voltage V α in lead α at time t = 0 meaning that the lead energy dispersion in Eq. ( 3) becomes kα → kα + V α .Due to this nonequilibrium condition, charge carriers start to flow through the molecular conducting channel, in our case, a GNR.We stress that the coupling matrix elements between the central region and the leads, T m,kα , are constant at all times as the system is partition-free.Here we consider only voltage biases as mean of perturbation but recently it has also been shown that temperature gradients may be included in this consideration at an equal footing [37,59,60].
In the literature, a considerable amount of works uses the method of Landauer and Büttiker (LB) [61,62] to determine the transport properties of nanoscale devices as it provides a very simple and intuitive physical picture of the transport mechanism.The current I αδ in lead δ is calculated from the scattering states originating from lead α = δ.These scattering amplitudes are typically written as transmission probabilities for an electron to traverse from lead α to lead δ.The stationary current in lead δ is obtained from the difference α =δ [I αδ −I δα ].Here we also address time-resolved currents, and the nonequilibrium Green's function (NEGF) FIG. 2. Distinct impurity configurations on a GNR host (cyan atoms).Impurities (red atoms) can be positioned (i) right on the top of a carbon atom (T), (ii) positioned over a carboncarbon bond characterizing a 'bridge' configuration (B), (iii) placed on the center of an hexagonal ring (C), or (iv) substitute a carbon atom (S).
formalism [63] provides a natural framework for this as it is not limited to the stationary state.We describe the time-evolution (transient information) of the system by the NEGF formalism where the Hamiltonian, the Green's function, and correlation effects (self-energies) are coupled by the integro-differential Kadanoff-Baym equations of motion [63].For the model system described above, an analytic solution for the time-dependent one-particle density matrix of the central region and for the timedependent current between the central region and the leads can be found [64,65].This time-dependent extension to the LB formalism (TD-LB) shares the simple interpretation of the original LB formalism and does not increase the computational cost as it would be the case if one solves numerically the complete Kadanoff-Baym set of equations [66,67].In addition, an arbitrary timedependence may also be included in the bias voltage, e.g., ac driving [68,69].We emphasize that our method allows for studying the transient and stationary regimes at an equal footing since the stationary LB formula is recovered at the long-time limit t → ∞.

III. RESULTS
We consider GNRs of varying widths with armchair edges in the transport direction (see Figs. 1 and 2).The widths of the GNRs studied here are N = {11, 12} with N indicating the number of carbon-dimers across the ribbon width representing, respectively, the metallic and semiconducting families of armchair GNRs: N = 3p − 1 and N = 3p with p an integer number [70,71].In addition to the pristine GNR, we consider impurities being adsorbed or substitutionally placed over the GNR host as shown in Fig. 2. The impurities can connect to the pristine GNR in four different configurations: 'Center' (C), 'Bridge' (B), 'Top' (T), and 'Substitutional' (S) [45].For the impurities we set imp = 0.66γ and γ imp = −2.2γ in Eq. ( 2) [54].The left-most atoms of the graphene nanoribbon are connected to the left lead and the rightmost atoms to the right lead, see Fig. GNR in the transport direction is 10 hexagons (≈ 4 nm) which is large compared to the impurity section.This choice is justified also because the overall transient features have been shown to scale with the length of the GNRs [64].For the sake of simplicity, we study transient responses when the GNR is subjected to adsorption of 4 impurities; these correspond to the four centermost hexagons in the GNR (Fig. 2).More complex chemical perturbations such as increase in the number of impurities, asymmetric bonds, random distribution of impurities etc. [72][73][74] could also be addressed with the same theoretical toolbox.However, our understanding can benefit from such simplified picture that enables us to address this problem with mathematical transparency.In this way, we can identify with clarity on typical transport patterns and universal responses in a class of impurities that attach to the graphene lattice following these particular bonding symmetries.
Current-voltage characteristics.-Westart by setting a source-drain voltage, V SD ≡ V L − V R , over a GNR section sandwiched by left (L) and right (R) leads and evaluate the stationary current with I = (I L + I R )/2 where I L and I R are the currents through the left and right interface, respectively.This gives I-V curves as seen in Fig. 3. From the current-voltage characteristics, we see that placing the impurities on the center configuration, in general, corresponds to impurity-invisibility [45] as the curves are essentially on top of the pristine ones.We have tested with different on-site and hopping energy parameters for the impurity that this symmetric configuration gives rise to vanishingly small scattering regardless of the specifics of the impurity.The only deviation between the pristine and center configurations is seen in a very small bias voltage window in the N = 12 ribbon (see inset in Fig. 3(b)).This effect could be related to Anderson localization, due to disorder induced by the impurities, leading to a metal-insulator transition [75,76].In general, from Fig. 3 we confirm that the P N = 11 ribbon is metallic (nonzero slope at zero bias), and doping with the B configuration makes the ribbon more metallic-like while doping with the T configuration makes the ribbon more semiconducting-like.Also, the P N = 12 ribbon is semiconducting and the other configurations show that the absolute values of the stationary currents are smaller, i.e., the impurities introduce scattering to some extent.
Current transients.-Nowwe investigate, how the stationary state in Fig. 3 is reached from the transient regime.As we observed some albeit small different behavior at small and large voltages, we performed transient calculations also in these two regimes by fixing the bias voltage to V SD /2 = 0.05γ and V SD /2 = 1.2γ, respec- tively.The time-dependent current signals are shown in Fig. 4 where we depict the initial transient behaviour (up to 30 fs) and also the long-time limit (up to 500 fs) at which we observed saturation of the currents.At small voltages there is considerable difference between the P and C configurations during the transient.As expected from Fig. 3, they saturate to the same value for the N = 11 ribbon and to a different value for the N = 12 ribbon.Larger voltages (Fig. 5) bring more transient features but we see that in both N = 11 and N = 12 ribbons the P and C configurations saturate essentially to the same value, although during the transient they can be very different.At higher voltages some configurations (T and B) take very long times to saturate (hundreds of femtoseconds); these impurity states introduce a considerable amount of back-and-forth scattering, which is not coupled to the leads, and these states have a long lifetime resulting in very slow damping.
Local charge densities and bond currents.-Aswe have seen above, even if the time-dependent currents through the nanoribbons saturate to expected values from stationary calculations, during the transient the currents oscillate significantly.In addition to the interface currents between the nanoribbons and the leads, we now investigate the local charge fluctuations and bond current patterns within the entire samples.In order to access this information, we need to evaluate the full one-particle density matrix, where the diagonal elements correspond to the local site densities and the off-diagonal elements correspond to the bond currents between the sites [64].We concentrate on the initial transient to understand how the role of impurities affects the formation of the stationary state, and we focus our discussion on two representative impurity configurations, C and B. Complete results are shown in the Supplementary Information [77].
From the time-dependent currents in Figs. 4 and 5 and from the snapshots in Figs. 6 and 7 we see that, even though the stationary current through the C configuration is mostly unaffected by the impurity sites, in the transient regime the impurities provide a "shock absorber" for the initial density wavefront.In the C configuration the density wavefronts undergo a symmetrydriven destructive interference, and the opposing bond currents cancel each other.This transparency is lifted once the lattice symmetry of the system is broken in other impurity configurations.This effect is observed by the decreased initial current peak at small voltages in Fig. 4, and as modified transient oscillations at high voltages in Fig. 5.The snapshots in Figs. 6 and 7 show the density variation (with respect to the ground state) and bondcurrent profiles before the first collision of the density wavefronts at the middle of the ribbons (t = 2 fs), and later when the wavepackets are reflecting from the lead interfaces back to the middle of the ribbons (t = 22 fs).We see how the B configuration introduces a peculiarly locked current pattern around the impurity atoms.This effect results in a remarkably partitioned charge distribution compared to the C configuration.
The full dynamics is better visualized by animations in the Supplementary Information [77].From the animations we may also see how the overall current pathways through the ribbon are modified by the impurities in real time.We note in passing that we observe, similar to Ref. [45], that the top-bonded impurities are strong scatterers compared to other impurity configurations, and the local bond-current profiles significantly rearrange and focus due to the impurities [30,77].
Transient charge pumping.-Wehave seen above how breaking the lattice symmetry of the unbiased graphene sample by introducing impurities leads to different signals both in the transient and stationary regimes.The symmetry of the transport setup can also be broken by the driving mechanism, and thus, we consider charge pumping through the graphene samples [69,[78][79][80][81][82][83] in the transient regime.In contrast to the previous sections, we now introduce a harmonic driving V L = −V R = V (t) with the voltage profile where V 0 is the source-drain dc voltage and A the amplitude of the ac driving.We set the driving frequency to be Ω = γ/10 which corresponds to a period of 2π/Ω ≈ 15 fs.We consider two types of ac driving: (1) V 0 = 0, A = γ, i.e., ac driving around zero dc voltage with odd inversion symmetry of the voltage profile, V (t + π/Ω) = −V (t); and (2) V 0 = γ/2, A = γ/2, i.e., breaking the oddinversion symmetry of the applied bias with a constant shift term.In Fig. 8 we show, for the N = 11 ribbons, the current responses to the two different drives described above, and also the corresponding Fourier spectra.For better frequency resolution the Fourier transform is calculated from an extended temporal window up to 500 fs, and Blackman-window filtering is used.We see from Fig. 8(a) and Fig. 8(c) that even if the time-dependent signals are essentially on top of each other for all the impurity configurations (within this temporal window), the frequency content is still different.The pristine ribbon expectedly excludes the even harmonics, showing pronounced peaks at ω = (2n + 1)Ω only, due to the odd-inversionsymmetric drive.However, introducing any impurities (even in the C configuration) breaks the corresponding symmetry of the time-independent Hamiltonian, and peaks at even multiples (ω = 2nΩ) of the driving frequency also appear.In Fig. 8(b) the odd-inversion symmetry of the drive is already broken, so all the ribbons including the P show pronounced peaks in the Fourier spectra also at the even multiples of the driving frequency.In addition, we see that in both cases peaks up to very high harmonic order are visible, indicating operation far beyond the linear-response regime.

IV. CONCLUSIONS
We presented a time-resolved characterization of impurity invisibility in graphene nanoribbons.Our transport setup of graphene nanoribbons supplemented with impurities was described by a single π orbital tight-binding framework where the impurity atoms were modeled by modified tight-binding parameters compared to the pristine graphene nanoribbons.We accessed the transport properties both at stationary and transient regimes by the TD-LB formalism, allowing for a fast and accurate simulation based on the NEGF method [63,84].
Our stationary results showed that the center-bonded impurities in graphene are invisible to conduction electrons being unable to scatter them [45,46,76].We then compared the time-dependent build-up of a steadystate current after a sudden quench of the bias voltage for different impurity configurations, and we discovered that the dynamics for different configurations look significantly different.Further, our spatio-temporal-resolved results showed that the impurities induce rearrangement and focusing of the current pathways along the graphene nanoribbons.In addition to the stationary picture, we further argue that graphene nanoribbons could serve as excellent probes or chemical sensors via ultrafast transport measurements [16][17][18][19][20].
Driving the graphene samples with strong ac bias voltage was shown to lead in highly nonlinear behavior.The resulting high-harmonic responses were shown to contain selective even-odd signals implying a generation of distinct on-off signals from an analog source which could be further realized as ac-dc conversion or rectification.These findings further highlight the great potential sensing devices probing ultrafast modifications in the sample, the general design of efficient circuitries, and engineering of, e.g., plasmonic and optical nanoscale devices [85][86][87][88][89].

FIG. 1 .
FIG. 1. Transport setup of an N = 11 armchair graphene nanoribbon with adsorbed impurity atoms (red spheres) and N indicating the number of carbon-dimers across the ribbon width.Contacts to the metallic leads are from the terminal sites of the nanoribbon.The leads are further connected to a source-drain voltage (VSD).

2 . 3 IFIG. 3 .
FIG. 3. Current-voltage characteristics of a (a) N = 11 and (b) N = 12 GNR in pristine and adsorbed/doped forms.Adatoms are placed on the center, bridge, top, and substitutional configurations as specified in Fig.2.The insets show a zoom-in at the low-voltage regime.

FIG. 4 .
FIG. 4. Time-dependent currents driven by a small bias voltage through the (a) N = 11 ribbons and (b) N = 12 ribbons.The long-time limit of the currents is shown as a cutout on the right-hand side.

FIG. 6 .
FIG. 6. Snapshots of the local charge densities and bond currents along the nanoribbons during the initial transient due to a small bias voltage, VSD/2 = 0.05γ.The charge densities are calculated as the difference from the ground-state density (color map).The strength of the bond current is indicated by the width of the black arrows.

FIG. 8 .
FIG. 8. Time-dependent currents through the N = 11 ribbons driven by ac bias voltage.Current responses to (a) an odd-inversion-symmetric drive, and (b) a broken-inversionsymmetric drive.(c) The corresponding Fourier transforms.