pubs.acs.org/NanoLett Letter Toward Optimized Charge Transport in Multilayer Reduced Graphene Oxides Mustafa Neşet Ç ınar, Aleandro Antidormi, Viet-Hung Nguyen, Alessandro Kovtun, Samuel Lara-Avila, Andrea Liscio, Jean-Christophe Charlier, Stephan Roche,* and Hal̂ dun Sevinçli* Cite This: Nano Lett. 2022, 22, 2202−2208 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: In the context of graphene-based composite applications, a complete understanding of charge conduction in multilayer reduced graphene oxides (rGO) is highly desirable. However, these rGO compounds are characterized by multiple and different sources of disorder depending on the chemical method used for their synthesis. Most importantly, the precise role of interlayer interaction in promoting or jeopardizing electronic flow remains unclear. Here, thanks to the development of a multiscale computational approach combining first-principles calculations with large- scale transport simulations, the transport scaling laws in multilayer rGO are unraveled, explaining why diffusion worsens with increasing film thickness. In contrast, contacted films are found to exhibit an opposite trend when the mean free path becomes shorter than the channel length, since conduction becomes predominantly driven by interlayer hopping. These predictions are favorably compared with experimental data and open a road toward the optimization of graphene-based composites with improved electrical conduction. KEYWORDS: disordered van der Waals thin films, reduced graphene oxides, charge transport, quantum transport, interlayer transport, multilayer transport scaling law U nderstanding charge transport in multilayered van der Waals materials has become an attractive and challenging problem, in the perspective of both fundamental and applied rGO nanosheets randomly stacked. Besides, the resulting localization length is also found to increase (by 3 orders of magnitude) with both the aromatic content and the thickness research on graphene-based composites.1,2,34 Indeed, gra- of the thin films as well, while being roughly independent of phene-related materials (including chemically disordered the lateral size nanosheet. Accordingly, the main contribution graphene like rGO) have shown remarkable capability to to transport properties in multilayered rGO with random improve charge and thermal conductivities of many insulating stacking likely stems from a bulk contribution, with marginal flexible materials such as organic polymers, suggesting them as nanosheet edge effects.5 the privileged filler material to reinforce, diversify, and improve Here, we investigate charge transport in rGO thin films the properties and performances of traditional materials used using state-of-the-art modeling techniques and compared the in wearables, flexible electronics, conducting textiles, and results with experimental measurements. The multilayered thermoplastics. However, accurate understanding of the rGO models are constructed using a rectangular ribbon-like microscopic mechanisms leading to transport in these complex geometry with a width of 20 nm and periodic boundary systems is still a matter of debate. Indeed, recent experimental conditions in the transverse direction. Disorder is introduced studies3−5 suggested the key role played by various types of by incorporating random distributions of chemical defects defects, as well as interlayer interaction in the transport (Figure 1a) such as divacancies (0.29%), Stone−Wales defects mechanisms taking place inside the graphene-based compo- (0.01%), and epoxides (4.7%). Such chemical nature and sites. densities of defects have been extracted from atomistic samples More specifically, transport characteristics in multilayered of rGO obtained by classical molecular simulations of the rGO5 were observed to change from a conventional Efros- Shkloskii variable range hopping (ES-VRH)35,36 to a temper- Received: October 11, 2021 ature-dependent power-law regime, with increasing number of Revised: February 22, 2022 stacked layers. Such findings do not apparently depend on any Published: March 1, 2022 length scale of the system, since they are observed in both micrometric networks of few nanosheets partially overlapping, as well as in centimeter-scale thin films built from billions of © 2022 American Chemical Society https://doi.org/10.1021/acs.nanolett.1c03883 2202 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter Figure 1. Reduced graphene oxide hosts different types of defects that are randomly distributed (a), where carbon atoms are colored with gray, and oxygen with red. Pentagon, heptagon, and octagons are indicated with red, blue, and green shades. Charge transport calculations in typical multilayered rGO models are investigated in rGO using (b) the Kubo formalism to estimate the bulk electronic conductivity and (c) the Landauer- Büttiker approach in a conventional device geometry. (d) Atomic structures of various incorporated defect types (pristine divacancy (585), reconstructed divacancy (555−777), Stone−Wales (55−77), and epoxy defects) and their corresponding tight-binding models (e) extracted from ab initio band-structure calculations. thermal reduction process of graphene oxide (GO) sheets.6 nearest-neighbor hopping energy γ0 = −2.6 eV, nearest- More specifically, to produce the models, we employ the neighbor distance r0 = 1.42 Å, and the decay parameter β = thermal annealing protocol described in ref 6, since it 3.37.15 The local doping due to the localized states induced by reproduces the main structural features observed in measured defects and impurities is also included by modulating onsite rGO samples. In these MD simulations, large-scale GO energies as a distance-decay function centered at the defect samples with an initial equivalent concentration of epoxide positions. Since the atomic positions are altered around and hydroxyl groups, amounting to a O/C ratio of 35%, were defects, the parametrization for interlayer coupling also annealed at 900 °C. At the end of the annealing process, the requires one to account for changes in interatomic distances final concentration of oxygen atoms was found to be compared to Bernal graphite (AB-stacking). Consequently, to approximately 5%, in excellent agreement with the exper- determine TB couplings in the presence of varying bond imental results reported in ref 5. lengths around structural defects, an exponential decay-based These chemical defects are known to induce lattice model formula is used following γ(r) = γ1 exp(βz(1 − r/z)) distortion and charge redistribution locally around their spatial with interlayer coupling energy γ1 = 0.36 eV, the corresponding location.7−12 Since these detrimental effects are local, rGO decay parameter βz = 24.99, and interlayer distance z = 3.34 models with defect concentrations similar to those investigated Å.16,17 Further details and comparison against density in this study can be modeled using a conventional tight- functional theory and other methods are given in the binding (TB) Hamiltonian (as for graphene) but with specific Supporting Information. The different types of considered adjustments made locally around the defect position. Within defects used in the model are illustrated in Figure 1d. The such a framework, the Hamiltonian of the rGO system presents electronic bands obtained from TB parametrization are in a simple form allowing one to further perform large-scale excellent agreement with first-principles density functional transport calculations in realistic samples containing more than theory results (Figure 1e). Dependence of electronic structure 106 atoms and with varying defect concentrations. The and transport properties of twisted graphene layers has been parametrized pz TB Hamiltonian reads studied intensely during the last years, and it was shown that even small twist angles give rise to substantial changes.18−20,32 H0 = ∑ γpqcp†cq A similar effect is not expected in randomly stacked rGO, ⟨p , q⟩ (1) because defects are strong scattering centers compared to the where in-plane couplings are limited to nearest-neighboring weak interlayer coupling. Hence, scattering from impurities is interactions. In contrast with other models in the literature the dominant mechanism for the randomization of carriers’ (see refs 13, 14), the effects of chemical defects are described momenta causing more than 100 scatterings within the by properly adjusting the TB parameters to recover ab initio interlayer diffusion length. Indeed, conductivity was shown results. The optimization of these TB parameters for each to be practically independent of the twist angle in twisted specific single defect are obtained by fitting first-principles bilayer graphene in the presence of lower concentrations of electronic band structures of the graphene supercell containing random defects.21 a single defect (see Supporting Information). Kubo-Greenwood (KG) and Landauer-Büttiker (LB) trans- Importantly, the change in C−C bond length due to the port formalisms are used to study the electronic conductivity lattice distortion around the defects is accounted for by and conductance in multilayered rGO models and for varying computing the hopping energies as a function of C−C bond transport geometries (see Figure 1b,c, respectively). The LB length (rpq) and defining γpq = γ0 exp[−β(rpq/r0 − 1)] with method22−24 allows one to include the charge injection from 2203 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter Figure 2. (a) Electronic conductance of mono/bi/trilayer rGO devices (blue, red, green curves, respectively) as a function of energy. The conductance is found to increase with the number of rGO layers, in contrast with stack of pristine graphene (inset). (b) Electronic conductance of mono- and trilayered rGO systems for E = 6.7 meV (lower region) and E = −73.8 meV (upper region) estimated by LB (markers) and eq 2 (solid lines). Mean free paths used in eq 2 are determined by using KG method. Diffusive versus localized regimes in multilayered rGO devices are distinguished by referring to the β-function (c). Transport regimes in pristine and defective stacks are illustrated in (d,e).(f) Local density of states (LDOS) around the Fermi level in multilayer rGO is shown for a relatively smaller sample. LDOS for individual layers and the multilayer are shown in the below panels. Red, green, and blue correspond to individual layers from bottom to top. contact electrodes (Figure 1c), whereas the KG method gives interference processes taking place across the layers and access to bulk properties25 (Figure 1b). These two techniques throughout the central region. enable contrasting bulk properties with “device” related In sharp contrast, in strongly disordered rGO systems, the transport (see Supporting Information for details). Note that electronic behavior is opposite, in the sense that the higher the neither electron−phonon coupling nor many-body effects are number of added layers on the stack, the larger the included in the present study. conductance (Figure 2a, main frame). This counterintuitive Figure 2a presents the energy-dependent conductance for result suggests that the interlayer hopping is opening more the mono-, bi-, and trilayer rGO for identical defect conductive channels once the localization length of the bottom concentration, calculated within the LB approach. When layers are free from defects (pristine graphene), adding new layers layer is short enough due to strong in-plane disorder. The two will act as a supplemental scattering source, hence reducing the opposite roles played by interlayer coupling in defect-free and total conductance (see Figure 2a-inset). Indeed, in defect-free defective multilayers are pictured in Figure 2d,e. monolayer graphene, transport is ballistic and a maximum V- To better understand the electronic transport mechanisms in shaped conductance curve is obtained. When more clean layers mono- and multilayered rGO systems, we perform a scaling are stacked on top of this graphene layer, new scattering study of the conductance to better differentiate between channels are opened and an interface resistance is formed. diffusive and localized regimes. Figure 2b shows results for the Consequently, the conductance of bi- and trilayers (multi- mono- to multilayered rGO systems presenting the same layered stack) is reduced below the single layer of pristine defect density. Clear evidence for different scaling behaviors is monolayer graphene. Actually, such conductance decay is not observed (note that conductance values from LB simulations specific to disorder-free systems but is also observed in low- are marked with circles (monolayer) and crosses (trilayer)). In defect concentrations as well (see Figure S6). Importantly, this the diffusion regime, the conductance dependence with the reduction of the electronic transmission is found when the electrodes are attached to the bottom layer, and it decreases in length is theoretically expected to scale as a function of the a non-monotonic way with the number of added layers. This mean-free-path (Smfp) and the defect-free conductance behavior can be understood as driven by complicated (GGraphene) of the system (at a given energy) following 2204 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter ij y jj1 + 2 L zzz jj z π S mfp z{ 1 1 GGraphene k = G (2) where L is device length.23 For short enough channel length, the computed con- ductance scaling is indeed well described by the diffusive formula (eq 2), allowing the extraction of Smfp from simulations. For the chosen disorder features, Smfp in the range of 1−2 nm is obtained (Figure 2b, inset), regardless of the chosen transport formalism. As the number of layers is increased, Smfp decreases slightly (Figure S5). We note that we report here the scaling behavior of G values for two selected energies (6.7 meV and −73.8 meV), but the trends are similar Figure 3. Localization lengths for varying numbers of layers. Circles at all energies around the charge neutrality point (CNP). represent experimental values, crosses are simulation results. Nlayer ≤ 6 We further observe that for channel length L > 10 nm, the stands for randomly stacked rGO sheets, and the last data point is for conductance values for the monolayer are below the diffusion partially oxidized graphite. Experimental Sloc values are calculated curve (which are plotted by solid lines). This pinpoints the using different dielectric constants. In the case of a single or a few onset of localization effects. However, for the trilayer case, a stacked layers, we estimate ϵr = 2.5 from the average of the dielectric puzzling sudden change of the conductance behavior is seen constant of the substrate (SiO2 = 3.9) and vacuum/air (= 1),26 while for a channel length of 10 nm. To further substantiate this for RGO thin film we assume ϵr = 3.5,27 and ϵr = 15 for partially striking difference of the conductance scaling between oxidized graphite. monolayer and multilayered stack, we investigate the local- ization regime and capture quantitative information about the for resistivity, ϵr is the relative dielectric constant of the localization lengths. LB simulations are performed on rGO material, and A = 2.8e2/4πϵ0kB = 0.021 μm−1 K−1 for a 2D systems with long channel lengths and using the scaling system (see Supporting Information for more details). rGO function defined as β = ∂ ln G/∂ ln L.28 Substituting the devices have Sloc values varying from ca. 4 to 30 nm, showing corresponding expressions for diffusion and localization, β can the same trend and the same order of magnitude when be written as βdif = −(1 + π Smfp/2L)−1 in the diffusive regime, compared with the simulations. It is noteworthy to emphasize whereas βloc = −L/Sloc in the localization regime with Sloc that in the case of partially oxidized graphite the corresponding being the localization length. Both scaling functions (βdif and Sloc value is 1 order of magnitude larger amounting to 250 nm, βloc) are presented in Figure 2c versus the device length. The clearly evidencing the combined role of the crystalline solid curves represent predictions of the aforementioned structure and the dielectric properties. Such aspects are out analytical formulas, whereas the markers are directly obtained of the scope of this work, and some details related to the from the simulation (for an energy of 6 meV) without referring experimental setup and the device characterizations are to Smfp or Sloc , since the scaling functions can be expressed as βdif reported in the Supporting Information. = G/GGraphene − 1 and βloc ∝ ln G. The substantial increase in the localization lengths is thus a On one hand, for the short-channel rGO devices, the key transport feature to distinguish between various multi- diffusive regime is confirmed since the 1/L behavior is layered rGO devices and which can be rationalized withing the observed when using Smfp (see Figure 2c, inset). On the variable range hopping (VRH) framework.29 Indeed, at zero temperature, the energy separation between two localized other hand, localization regime is obtained for longer-channel states should be very small to enable significant hopping, with rGO devices, and the evaluated localization lengths varies with probability proportional to e−2αR (with 1/α being the the number of layers: 7.8 nm for monolayer, 13.7 nm for attenuation length for the localized states and R being the bilayer, and 22.9 nm for trilayer. Thus, when increasing the spatial separation of states). The hopping probability between number of layers, hopping transport gives rise to larger two states ψA and ψB, which have similar energies and are transmission amplitudes compared to diffusion because of the localized at the bottom layer well-separated from each other, enhancement of Sloc , in very good agreement with the should thus be strongly enhanced if a third state ψC is localized experimental findings as shown in Figure 3. at the upper layer and between ψA and ψB. Therefore, layers Such behavior agrees well with experimental findings with localized states have enhanced transmission when stacked. achieved on rGO devices with different numbers of layers In Figure 2f, the local density of states (LDOS) is plotted for a (Nlayer) and similar chemical structure (i.e., sp2 content 96 ± trilayer rGO sample.33 It is clearly visible that states localized at 1%).5 Figure 3 collects the results obtained comparing seven different layers tend to fill the spatial gaps when they are different systems ranging from the single nanosheet (Nlayer = 1) superimposed. This can qualitatively explain the enhancement to a flake of partially oxidized graphite (Nlayer = 6). Differently, of Sloc with number of layers. all the other devices (2 ≤ Nlayer ≤ 6) are assemblies of rGO The interlayer coupling thus plays a crucial role in multilayer sheets randomly stacked. According with the results reported transport. In low defect concentrations, it is a source of in ref 5, all six devices reveal ES-VRH transport mechanisms at scatterings which impedes transport, whereas in highly low temperatures (10 K < T < 100 K). Thus, the defective rGO multilayers, it tends to promote longer state corresponding localization length (Sloc ) is calculated using the delocalization. The contact geometry is also a major factor in temperature-dependent electrical resistivity curves transport across the layers. To gain further information about ρ(T ) = ρ0,VRH exp{A ϵr SlocT}−1/2 where ρ0,VRH is a prefactor the combined roles of interlayer coupling and contact 2205 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter geometry, we simulate a situation where the device is made of the monolayer as ;multi = ; 1/ Nlayer mono . Correspondingly, the from overlapping bilayer rGO in which the injection and resistivities (ρ) satisfy the following relation 1 jij ρmulti zyz jij ρmono zyz collection takes place at different layers (inset of Figure 4). We jj zz = jj zz Nlayer jk ρ0 z{ jk ρ0 z{ 1/ N layer (3) where ρ0 = hAmono/2e L is a system-wide constant Amono being 2 the cross section area for a monolayer. Specifically, we predict a linear dependence of ln(ρmulti/Nlayer) ∼ 1/Nlayer for film thickness larger than the mean free path. For the sake of comparison, we analyze a data set of 11 rGO devices with similar chemical structure (i.e., sp2 content = 96 ± 1%) and film thickness ranging from 2 to 13 nm, i.e., 5 < Nlayer < 35. Note that we do not need to assume any values for Amono and ρ0 or measure them, but they are used in order to relate resistivity and conductance so as to compare experimentally measured values with simulated ones through the scaling relation (eq 3). Figure 5 shows the correlation plot ln(ρmulti/Nlayer) vs 1/ Nlayer for each device displaying the resistivity values acquired Figure 4. Conductance as a function of overlap length in bilayer reduced graphene oxide when the electrodes are connected to different layers. Dotted curves represent exponential decay. compute the conductance with different overlap distances, L, which is also the length of the central region. As shown in Figure 4, the conductance rapidly increases with L at short distances because the overlap area enhances the probability for a carrier to diffuse from one layer to the other. Differently, for overlap distances longer than 10 nm, the conductance decays exponentially with L, indicating that the localization behavior prevails over the interlayer diffusion. The maximum con- ductance is achieved at around 7 nm, which marks the interlayer diffusion length (Sinter ). We have checked the dependence of Sinter in clean structures, and found very similar Figure 5. Scaling of multilayer resistivity with the number of layers is values, which proves that Sinter is dictated by the strength of plotted as ln(ρmulti/Nlayer) versus 1/Nlayer. Experimental data of interlayer coupling and not by the disorder content. We can multilayer rGO (circles for devices reported in ref 5 and square for actually estimate a length scale for interlayer diffusion as hvF/γ1 partially oxidized graphite) acquired at different temperatures show = 11.5 nm, where the Fermi velocity is vF = 106 m/s and the linear dependence, in good agreement with the theoretical prediction interlayer coupling strength is γ1 = 0.36 eV. The comparison for scaling (cf., eq 3). All the linear fitting curves calculated at different temperatures are included between the two curves acquired at 5 and Sinter ≫ Smfp in the simulated structures is the reason behind the 300 K (dashed area). fact that Smfp is weakly affected from the increase in the number of layers. On the other hand, since Sinter is comparable with the monolayer localization length, there is room for localized at different temperatures: 5, 100, and 300 K; the lowest, an carriers to further spread over the neighboring layers. To clarify intermediate, and the highest measured, respectively. A total of the effect of contact geometry on transport, we have data sets corresponding to 43 different temperatures were considered the electrodes with the same number of layers as analyzed, and the remaining 40 curvesnot depicted in the in the central region. In this geometry, as expected, the figureare included between the two curves acquired at 5 and conductance increases with the number of layers in both clean 300 K (dashed area). In all the cases we observe a linear trend and defective samples, but the behavior of neither Smfp nor Sloc y = m·x + q in excellent agreement with eq 3. Moreover, we is considerably affected. determined that the slope m = ln(ρmono/ρ0) decreases with Finally, based on these findings, the scaling of the resistivity increasing temperature, while the y-intercept (q) is a constant with the number of layers can be expressed quantitatively. In value corresponding to ρ0 = exp(q) = (2.1 ± 0.2) × 10−7 Ω·m. highly defective rGO multilayers Smfp changes only slightly with Similarly, a linear behavior is achieved in the case of rGO devices with lower amount of the aromatic content (77% and the number of layers, but Sloc increases linearly with Nlayer. 86%, see Figure S9) where the resistivity values increase with Using the fact that Sloc is proportional to the number of the oxidation degree, as expected. To summarize the transmission channels and Smfp,30,31 one can approximate the experimental findings, ρmono is the single-layer resistivity transmission probability across multilayer rGO in terms of that (temperature-dependent), while ρ0 does not depend on the 2206 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter temperature, being therefore a kind of resistivity scaling factor Authors only depending on the aromatic content of the device. Mustafa Neşet Ç ınar − Department of Materials Science and ■ CONCLUSION We have reported quantum simulations on realistic models of Engineering, Izmir Institute of Technology, 35430 Urla, Izmir, Turkey Aleandro Antidormi − Catalan Institute of Nanoscience and multilayered rGO which reveal the complex interplay between Nanotechnology, CSIC and The Barcelona Institute of Science disorder and interlayer interactions in dictating the dominant and Technology, 08193 Bellaterra, Barcelona, Spain; transport mechanism. Depending on the concentration of orcid.org/0000-0002-5266-8147 defects, multilayer interactions can enhance or suppress the Viet-Hung Nguyen − Institute of Condensed Matter and system conductance, which results from the competition Nanosciences, Université catholique de Louvain between the mean free path Smfp and the interlayer diffusion (UCLouvain), B-1348 Louvain-la-Neuve, Belgium; orcid.org/0000-0001-6729-3520 length Sinter . When about 5% of the carbon atoms are involved Alessandro Kovtun − Consiglio Nazionale delle Ricerche, in defected regions, Sinter becomes much longer than Smfp. In Istituto per la Sintesi Organica e la Fotoreattivitá, (CNR- that case, intralayer scattering largely dominates over interlayer ISOF), 40129 Bologna, Italy; orcid.org/0000-0002- diffusion, leading to a weak dependence of Smfp on Nlayer. 7614-7100 On the other hand, Sloc ∼ Sinter , so that a localized state in Samuel Lara-Avila − Department of Microtechnology and one of the layers has enough extension for tunneling to an Nanoscience, Chalmers University of Technology, 41296 adjacent one. If Sinter was much larger than Sloc , the tunneling Gothenburg, Sweden; orcid.org/0000-0002-8331-718X rates would be much smaller. Once Sinter ∼ Sloc , tunneling rates Andrea Liscio − Consiglio Nazionale delle Ricerche, Istituto are appreciable and charge delocalization is promoted. While per la Microelettronica e Microsistemi, Roma Unit (CNR- Smfp is weakly dependent on Nlayer, Sloc increases with Nlayer as IMM), 00133 Rome, Italy; orcid.org/0000-0003-2986- 3398 expected from a generalization of the Thouless relationship for Jean-Christophe Charlier − Institute of Condensed Matter one-dimensional conductors. This unprecedented interplay and Nanosciences, Université catholique de Louvain between transport length scales is a specific result of 2D (UCLouvain), B-1348 Louvain-la-Neuve, Belgium layered nature of the multilayer rGO systems. Such a mechanism enables hopping transport to overcome the Complete contact information is available at: diffusion limit, which is usually the upper bound in bulk https://pubs.acs.org/10.1021/acs.nanolett.1c03883 systems. Our theoretical analysis enables us to derive a novel scaling rule, which is in perfect agreement with experimental Notes data at various temperatures and consistent with the Thouless The authors declare no competing financial interest. ■ relationship. The fundamental findings of this study are not limited to multilayered reduced graphene oxide but could find ACKNOWLEDGMENTS applications in other two-dimensional stacks as well. ■ The authors acknowledge support from the Flag-Era JTC 2017 ASSOCIATED CONTENT project ‘ModElling Charge and Heat trANsport in 2D- materIals based Composites  MECHANIC’. M.N.Ç . and * sı Supporting Information H.S. acknowledge support from TÜ BITAK ̇ (117F480). The The Supporting Information is available free of charge at numerical calculations reported in this paper were partially https://pubs.acs.org/doi/10.1021/acs.nanolett.1c03883. performed at TÜ BITAK ̇ ULAKBIM, ̇ High Performance and Kubo-Greenwood and Landauer-Bü ttiker methods; Grid Computing Center (TRUBA resources). A.A. and S.R. simulated structures; system setups; tight-binding are supported by MECHANIC reference number: PCI2018- parameters; supplementary mean-free-path plots; results 093120 funded by Ministerio de Ciencia, Innovacion y with lower defect concentrations; supplementary LDOS Universidades and the European Union Horizon 2020 research plots; temperature dependence of electrical resistivity and innovation programme under Grant Agreement No. ρ(T) ; Efros−Shklovskii variable range hopping model; 881603 (Graphene Flagship). ICN2 is funded by the details of electrical resistivity measurements; table of CERCA Programme/Generalitat de Catalunya, and is experimental parameters of all the studied devices supported by the Severo Ochoa program from Spanish (PDF) MINECO (Grant No. SEV-2017-0706). V.-H.N. and J.-C.C. acknowledge financial support from the Fédération Wallonie- ■ AUTHOR INFORMATION Corresponding Authors Bruxelles through the ARC Grants (No. 16/21-077 and No. 21/26-116), from the European Union’s Horizon 2020 Research Project and Innovation Program  Graphene Stephan Roche − Catalan Institute of Nanoscience and Flagship Core3 (No. 881603), from the Flag-ERA project Nanotechnology, CSIC and The Barcelona Institute of Science (No. R.8010.19), and from the Belgium FNRS through the and Technology, 08193 Bellaterra, Barcelona, Spain; research projects (No. T.0051.18 and No. T.029.22F). ICREA−Institució Catalana de Recerca i Estudis Avançats, Computational resources have been provided by the CISM 08010 Barcelona, Spain; orcid.org/0000-0003-0323- supercomputing facilities of UCLouvain and the CECI 4665; Email: stephan.roche@icn2.cat consortium funded by F.R.S. -FNRS of Belgium (No. Hâldun Sevinçli − Department of Materials Science and 2.5020.11). Computational resources have been provided by Engineering, Izmir Institute of Technology, 35430 Urla, the CISM supercomputing facilities of UCLouvain and the Izmir, Turkey; orcid.org/0000-0002-1896-2588; CECI consortium funded by F.R.S. -FNRS of Belgium (No. Email: haldunsevincli@iyte.edu.tr 2.5020.11). Authors are particularly grateful to Prof. Paolo 2207 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208 Nano Letters pubs.acs.org/NanoLett Letter Samorí and Marco Gobbi for the scientific assistance to (19) Carr, S.; Massatt, D.; Fang, S.; Cazeaux, P.; Luskin, M.; Kaxiras, prepare the devices, Valentina Mussi for some supporting E. Twistronics: Manipulating the electronic properties of two- measurements and for enlightening discussions. dimensional layered structures through their twist angle. Phys. Rev. ■ B 2017, 95, 075420. (20) Cao, Y.; Fatemi, V.; Fang, S.; Watanabe, K.; Taniguchi, T.; REFERENCES Kaxiras, E.; Jarillo-Herrero, P. Unconventional superconductivity in (1) Ferrari, A. C.; et al. Science and technology roadmap for magic-angle graphene superlattices. Nature 2018, 556, 43−50. graphene, related two-dimensional crystals, and hybrid systems. (21) Namarvar, O. F.; Missaoui, A.; Magaud, L.; Mayou, D.; Nanoscale 2015, 7, 4598−4810. Trambly de Laissardière, G. Electronic structure and quantum (2) Mohan, V. B.; Lau, K.-t.; Hui, D.; Bhattacharyya, D. Graphene- transport in twisted bilayer graphene with resonant scatterers. Phys. based materials and their composites: A review on production, Rev. B 2020, 101, 245407. applications and product limitations. Composites Part B: Engineering (22) Landauer, R. Spatial Variation of Currents and Fields Due to Localized Scatterers in Metallic Conduction. IBM J. Res. Dev. 1957, 1, 2018, 142, 200−220. 223−231. (3) Turchanin, A.; Weber, D.; Büenfeld, M.; Kisielowski, C.; Fistul, (23) Datta, S. Electronic Transport in Mesoscopic Systems; Cambridge M. V.; Efetov, K. B.; Weimann, T.; Stosch, R.; Mayer, J.; Gölzhäuser, Studies in Semiconductor Physics and Microelectronic Engineering; A. Conversion of Self-Assembled Monolayers into Nanocrystalline Cambridge University Press, 1995. Graphene: Structure and Electric Transport. ACS Nano 2011, 5, (24) Ryndyk, D. Theory of Quantum Transport at Nanoscale: An 3896−3904. Introduction; Springer, 2015. (4) Silverstein, K. W.; Halbig, C. E.; Mehta, J. S.; Sharma, A.; Eigler, (25) Fan, Z.; Garcia, J. H.; Cummings, A. W.; Barrios-Vargas, J. E.; S.; Mativetsky, J. M. Voltage-reduced low-defect graphene oxide: a Panhans, M.; Harju, A.; Ortmann, F.; Roche, S. Linear scaling high conductivity, near-zero temperature coefficient of resistance quantum transport methodologies. Phys. Rep. 2021, 903, 1−69. material. Nanoscale 2019, 11, 3112−3116. (26) Jang, C.; Adam, S.; Chen, J.-H.; Williams, E. D.; Das Sarma, S.; (5) Kovtun, A.; Candini, A.; Vianelli, A.; Boschi, A.; Dell’Elce, S.; Fuhrer, M. S. Tuning the Effective Fine Structure Constant in Gobbi, M.; Kim, K. H.; Lara Avila, S.; Samorí, P.; Affronte, M.; Liscio, Graphene: Opposing Effects of Dielectric Screening on Short- and A.; Palermo, V. Multiscale Charge Transport in van der Waals Thin Long-Range Potential Scattering. Phys. Rev. Lett. 2008, 101, 146805. Films: Reduced Graphene Oxide as a Case Study. ACS Nano 2021, (27) Jung, I.; Vaupel, M.; Pelton, M.; Piner, R.; Dikin, D. A.; 15, 2654−2667. Stankovich, S.; An, J.; Ruoff, R. S. Characterization of Thermally (6) Antidormi, A.; Roche, S.; Colombo, L. Impact of Oxidation Reduced Graphene Oxide by Imaging Ellipsometry. J. Phys. Chem. C Morphology on Reduced Graphene Oxides upon Thermal Annealing. 2008, 112, 8499−8506. J. Phys.: Mater. 2020, 3, 015011. (28) Abrahams, E.; Anderson, P. W.; Licciardello, D. C.; (7) Cresti, A.; Nemec, N.; Biel, B.; Niebler, G.; Triozon, F.; Ramakrishnan, T. V. Scaling Theory of Localization: Absence of Cuniberti, G.; Roche, S. Charge transport in disordered graphene- Quantum Diffusion in Two Dimensions. Phys. Rev. Lett. 1979, 42, based low dimensional materials. Nano Research 2008, 1, 361−394. 673−676. (8) Zhang, Y.; Brar, V. W.; Girit, C.; Zettl, A.; Crommie, M. F. (29) Mott, N. F. Conduction in non-crystalline materials: III. Origin of spatial charge inhomogeneity in graphene. Nat. Phys. 2009, Localized states in a pseudogap and near extremities of conduction 5, 722−726. and valence bands. Philos. Mag. 1969, 19, 835−852. (9) Wang, S.; Wang, R.; Wang, X.; Zhang, D.; Qiu, X. Nanoscale (30) Thouless, D. J. Localization distance and mean free path in one- charge distribution and energy band modification in defect-patterned dimensional disordered systems. Journal of Physics C: Solid State graphene. Nanoscale 2012, 4, 2651−2657. Physics 1973, 6, L49−L51. (10) Ma, C.; Sun, H.; Zhao, Y.; Li, B.; Li, Q.; Zhao, A.; Wang, X.; (31) Leconte, N.; Lherbier, A.; Varchon, F.; Ordejon, P.; Roche, S.; Luo, Y.; Yang, J.; Wang, B.; Hou, J. G. Evidence of van Hove Charlier, J.-C. Quantum transport in chemically modified two- Singularities in Ordered Grain Boundaries of Graphene. Phys. Rev. dimensional graphene: From minimal conductivity to Anderson Lett. 2014, 112, 226802. localization. Phys. Rev. B 2011, 84, 235420. (11) Gao, Y.; Qin, C.; Qiao, Z.; Wang, B.; Li, W.; Zhang, G.; Chen, (32) Trambly de Laissardiere, G.; Mayou, D.; Magaud, L. R.; Xiao, L.; Jia, S. Observing and tuning the density distribution of Localization of Dirac Electrons in Rotated Graphene Bilayers. Nano localized states of monolayer graphene oxide by using external electric Lett. 2010, 10, 804−808. field. Appl. Phys. Lett. 2015, 106, 131103. (33) Datta, S. Nanoscale device modeling: the Green’s function (12) Zhang, W.; Lu, W.-C.; Zhang, H.-X.; Ho, K. M.; Wang, C. Z. method. Superlattices Microstruct. 2000, 28, 253−278. Lattice distortion and electron charge redistribution induced by (34) Cheah, C. Y.; Gómez-Navarro, C.; Jaurigue, L. C.; Kaiser, A. B. defects in graphene. Carbon 2016, 110, 330−335. Conductance of partially disordered graphene: crossover from (13) Leconte, N.; Moser, J.; Ordejón, P.; Tao, H.; Lherbier, A.; temperature-dependent to field-dependent variable-range hopping. J. Bachtold, A.; Alsina, F.; Sotomayor Torres, C. M.; Charlier, J.-C.; Phys.: Condens. Matter 2013, 25, 465303. Roche, S. Damaging Graphene with Ozone Treatment: A Chemically (35) Mott, N. F. On the transition to metallic conduction in Tunable Metal-Insulator Transition. ACS Nano 2010, 4, 4033−4038. semiconductors. Can. J. Phys. 1956, 34, 1356−1368. (14) Lherbier, A.; Dubois, S. M.-M.; Declerck, X.; Niquet, Y.-M.; (36) Efros, A. L.; Shklovskii, B. I. Coulomb gap and low temperature Roche, S.; Charlier, J.-C. Transport properties of graphene containing conductivity of disordered systems. Journal of Physics C: Solid State structural defects. Phys. Rev. B 2012, 86, 075402. Physics 1975, 8, L49−L51. (15) Pereira, V. M.; Castro Neto, A. H.; Peres, N. M. R. Tight- binding approach to uniaxial strain in graphene. Phys. Rev. B 2009, 80, 045401. (16) Suárez Morell, E.; Correa, J. D.; Vargas, P.; Pacheco, M.; Barticevic, Z. Flat bands in slightly twisted bilayer graphene: Tight- binding calculations. Phys. Rev. B 2010, 82, 121407. (17) Reich, S.; Maultzsch, J.; Thomsen, C.; Ordejón, P. Tight- binding description of graphene. Phys. Rev. B 2002, 66, 035412. (18) Suárez Morell, E.; Correa, J. D.; Vargas, P.; Pacheco, M.; Barticevic, Z. Flat bands in slightly twisted bilayer graphene: Tight- binding calculations. Phys. Rev. B 2010, 82, 121407. 2208 https://doi.org/10.1021/acs.nanolett.1c03883 Nano Lett. 2022, 22, 2202−2208