Multiscale charge transport in van der Waals thin films: reduced graphene oxide as case study

Large area van der Waals (vdW) thin films are assembled materials consisting of a network of randomly stacked nanosheets. The multi-scale structure and the two-dimensional nature of the building block mean that interfaces naturally play a crucial role in the charge transport of such thin films. While single or few stacked nanosheets (i.e. vdW heterostructures) have been the subject of intensive works, little is known about how charges travel through multilayered, more disordered networks. Here we report a comprehensive study of a prototypical system given by networks of randomly stacked reduced graphene oxide 2D nanosheets, whose chemical and geometrical properties can be controlled independently, permitting to explore percolated networks ranging from a single nanosheet to some billions with room temperature resistivity spanning from 10-5 to 10-1 ohm m. We systematically observe a clear transition between two different regimes at a critical temperature T*: Efros-Shklovskii variable range hopping (ESVRH) below T* and power law (PL) behavior above. Firstly, we demonstrate that the two regimes are strongly correlated with each other, both depending on the charge localization length xi, calculated by ES-VRH model, which corresponds to the characteristic size of overlapping sp2 domains belonging to different nanosheets. Thus, we propose a microscopic model describing the charge transport as a geometrical phase transition, given by the metal-insulator transition associated with the percolation of quasi-1D nanofillers with length xi, showing that the charge transport behavior of the networks is valid for all geometries and defects of the nanosheets, ultimately suggesting a generalized description on vdW and disordered thin films.

domains belonging to different nanosheets. Thus, we propose a microscopic model describing the charge transport as a geometrical phase transition, given by the metal-insulator transition associated with the percolation of quasi-1D nanofillers with length ξ, showing that the charge transport behavior of the networks is valid for all geometries and defects of the nanosheets, ultimately suggesting a generalized description on vdW and disordered thin films.
The development of cheap techniques to produce and to process large quantities of monoatomic thick materials such as graphene 1 and related 2D-materials allowed to create synthetic structures, i.e. van der Waals (vdW) materials, with tuned properties. [2][3] Using solutionprocessing approaches, 4-6 the produced nanosheets can be arranged forming macroscopic vdW thin films composed of billions of 2D sheets randomly stacked, in which the structural and electrical continuity of the 2D/3D architectures is provided by the contact regions between the sheets. The versatility and the processability of such assembles allow to control with great flexibility both the crystal quality of the sheets and their stacking geometry (e.g. thin films, membranes, foams, etc.), holding a significant potential for the emerging applications based on large-area flexible and wearable devices, 7 coatings and advanced composites. 8 Despite of the great technological developments, a general framework describing the charge transport (CT) in vdW thin films is still lacking. Most of the studies carried out during the last 15 years provided an in-depth insight into CT in individual almost defect-free nanosheets or stacks obtained by superimposing only a finite number of them fin the so-called vdW heterostructures. 9 The CT in a network differs intrinsically from that of the single sheets because of the major role played by inter-sheet processes. A model taking into account the superimposition of individual sheets represents an over demanding computational task since it would require a detailed description of the mechanisms involved at interfaces on different length and time scales.
Moreover, due to the need of finding the ideal trade-off between costs of production on large scale and competitive properties, the individual sheets forming the film necessarily possess intrinsically lower quality than the pristine (greater amount of defects and chemical functionalization, higher size polydispersity), thus further increasing the resources required for their simulation.
The experimental approaches used so far suffered primarily from the poor processability of graphene. In particular, macroscopic films made by exfoliated graphene dispersions contain significant amounts of multilayers. The presence of residual contaminants or dopants is also a factor that limits the ability to develop comprehensive models. Instead, the use of pure 2dimensional, monoatomic nanosheets as (semi)conductors would enable to unravel the correlation between system dimensionality, nanostructuring of the interfaces and charge transport in complex networks.
It is worth taking into account that networks of 2D materials obtained by overlapping sheets are characterized by 2D interfaces which can be as large as the area of the individual sheets.
Hence, each of the atom composing the material will lie at the interface, in contact with the sheets above or beneath it (figure 1a), being in sharp contrast with conventional 3D assemblies or composite materials. Networks of stacked 2D nanosheets can thus be defined as a multi-junction material where the CT mechanisms acting on the single sheet are strongly involved with those at the interface (sheet-to-sheet).
As prototypical material we used single monolayer sheets of reduced graphene oxide (RGO), which consists of a (semi)conductive graphene lattice comprising oxygen-containing defects. [10][11] Differently from solution-processed exfoliated graphene, such sheets can be obtained in the monodisperse form of single layers with their lateral size tunable from >100 nm to ≈ 100 µm by sonication treatments. 12 These nanosheets can thus be dispersed in various liquids, processed varying their size, film thickness hence their conductivity, 13 and easily arranged on a substrate forming networks with partially ordered structure, where all the single sheets are randomly distributed in the plane orientation and stacked perfectly parallel at a fixed distance (i.e. RGO interlayer distance) forming thus RGO thin films. In addition to the study of the CT mechanisms of vdW thin films, the capability to tune the chemico/physical properties of the building block and as well as to fabricate large-area networks in controlled way allowed to use RGO networks as test-bed systems to investigate the CT mechanisms in disordered and amorphous materials and as well as to verify the wide range of theories developed; see 14  Here we present a scale-independent model for the charge transport in highly disordered networks of defective 2D materials based on a robust data analysis of the electrical resistivity vs temperature ρ(T). This resistivity is a reliable physical observable to achieve in depth insight into single sheet to large-area aggregates. Moreover, we show that CT properties of such network can be described in a more general framework of the geometrical phase transition. By taking full advantage of the controlled structural and morphological composition of the films we carried out a multiscale quantitative analysis to elucidate the different transport regimes taking place in such materials. In particular, we investigated the CT mechanisms occurring at the sheet-to-sheet interfaces, typically considered as bottlenecks, as well as the role of the geometrical complexity of the network on the overall electrical conductivity of the nanosheets assemblies. The use of RGO enabled us to tune independently i) the conductivity of each nanosheet, ii) the lateral size of the nanosheets, and iii) the thickness of the material, i.e. the number of nanosheets stacked over each other. 15

RESULTS AND DISCUSSIONS
The starting material were water solutions of single-layer GO nanosheets. These solutions were spun onto SiO2 substrates at different concentrations, followed by a sample's reduction via thermal annealing at increasing temperature (T ann = from 200 °C up to 940 °C) for 60 minutes at a vacuum of 10 -6 mbar to transform them in conductive RGO. Finally, gold electrodes were deposited by thermal evaporation through a shadow mask. Upon changing systematically such experimental parameters, a set of 28 different devices was produced. At lower nanosheet density, micrometer sized highly-reduced RGO sheets (sp 2 content = 96±2% as determined by X-ray Photoelectron Spectroscopy) were connected by rare and clearly visible contact points with a size to form networks on micrometer scale (figure 1b). Differently, at higher deposition density, continuous layers of stacked sheets (figure 1c), were produced with tunable oxidation degrees, sheet size tuned from <0.3 µm to >17 µm, and film thickness from 1 up to 35 layers. All the fabricated samples (a.k.a. devices) were studied by microscopic and spectroscopic techniques (see Methods) at each step of preparation.
Preliminary Field-Effect Transistor (FET) measurements performed on single and sparse network of micrometer sized highly-reduced RGO sheets (figure S5) showed hole transport with linear region mobility ( , ) of ca 1 -10 cm 2 ⋅V -1 ⋅s -1 (figure 1e). Despite of the analysis of the , vs temperature is the most common approach to study CT mechanism, such physical observable could not be the most useful to compare systems with different geometries and scalelengths. Moreover, thin films with channel length and width of ca. 1 cm are strongly affected by the edge effects. For this reason, we focussed our attention to the study of the electrical resistivity vs temperature for each device from ca 10 K to 300 K, to detect the existence of different transport regimes. In all the 28 different devices we found that decreases with the increasing of temperature, clearly indicating that all the networks possess a semiconducting behavior, with room temperature resistivity ρ RT lying in the range between 2 and 2⋅10 -5 Ω⋅m. All the acquired curves are reported in the Supporting Information.

Observation of CT regimes at different temperatures. Typically, CT in disordered and
amorphous organic semiconductors is modeled as hopping transport between localized states.
Although the literature related to the transport phenomena is enormously rich and several models have been developed to study this physical framework, we can identify two "common" functional dependences that were used and combined to describe the measured ρ(T) trends: i) a stretched exponential and ii) a power law (PL). The first form is used to describe Variable Range Hopping (VRH) models [16][17] and Fluctuation-Induced Tunneling (FIT) 18 transports, typically observed at low temperatures. The latter form is observed in a wide range of systems and temperatures and it is related to high-density conductivity of disordered semiconductors (see Ref. From the experimental and data analysis point-of-view, a critical issue is related to the fitting analysis since the use of stretched exponential or PL functional forms describing ρ(T) can lead to systematic errors. For example, the application of the commonly used least-squares procedures for fitting a Poisson-distributed data is known to lead to biases. 19 where ∆= ln18 • # . 3.
Both VRH and PL regimes showed similar functional forms being defined by three parameters: a prefactor (ρ 0,VRH or ρ 0,PL ), a characteristic temperature (T 0 or T 1 ), and a characteristic exponent (β or m).
In the following paragraphs, we describe the role and the physical details of such parameters (i.e. the six ones defined in eqn.1 and T*) using a generalized description of the experimental features and correlating the data obtained in all the RGO networks. The data analysis procedure was performed as follows: first, we fitted the W(T) curves by calculating β, T 0 and m as free parameters, then we set these parameters and calculated the resistivity prefactors by fitting the ρ(T) curves. Finally, we fitted both W(T) and ρ(T) using the parameters previously calculated as initial values, verifying their convergence on successive reiterations.
VRH regime. In VRH the stretching exponent p is strongly dependent on the shape of density of states at Fermi energy, g(E F ). As reported in figure 3b, all the devices showed the same β values = 0.52±0.06, being in excellent agreement with results obtained on single RGO sheets [28][29] and with the values measured on hydrogenated graphene (β = 0.47 -0.58). 30 All the β values were always close to 1/2, providing an unambiguous answer to the long- conductors (e.g., quantum wires, nanotubes or polymers), 32 but is apparently insignificant in the case of the carbon network structures. 33 Conversely, ES-VRH dominates at all measurable temperatures in case of high-disorder systems and it has been already observed in RGO sheets. 28 The characteristic temperature of the hopping mechanism # is inversely proportional to the localization length (ξ) defined as the average spatial extension of the charge carrier wavefunction. According to the ES-VRH regime, the relation is given by # = Differently from the VRH regime, in PL regime only m values could be evaluated using the W(T) analysis (figure S7-S14). The m values observed ranged between 0.2 and 4.0 (table S3) increasing with the measured ρ RT of each network ( Figure 3c) and the PL characteristic temperature T 1 could not be directly decoupled from the prefactor ρ 0,PL . The resistivity measured in the PL regime is often described, in literature using a general function = L • 45 , where B is the scale factor. We could combine such general function with eqn. 1 to obtain: It is noteworthy to underline that eqn.2 is a mathematical equivalence between the parameters of each individual i-th RGO network. Therefore, we should consider a sequence of independent equivalences, one for each measured sample, where m and ln L parameters are directly calculated by W(T) and ρ(T) analysis, respectively. Moreover, each equivalence showed a solution given by a linear combination of ln and ln #,/0 . Another approach is based on the use of a universal scaling curve, 41 being such curve commonly observed in low-defects semiconductive polymers. The CT along the polymeric backbone is influenced by the coupling strength of the charge and the nuclear vibrations (i.e. a phonon bath). This mechanism is described by the quantum nuclear tunneling (NT) theory, in which the carriers tunnel through the potential barrier formed by the coupling of the electronic charge to its nuclear environment. Such coupling constant is called Kondo parameter, defined as O P = 5 + 1. 41 The higher the value, the stronger the coupling. Typical values reported in literature in π-conjugated polymers for bath phonon energies are ∼100 meV 42 and α K within the range 1.6 -6.75, 43 corresponding to m values between ca. 1.2 and 11.5. An alternative description of PL transport, assuming multiphonon tunnelling of localized electrons with a weak electron-lattice interaction, 44 could be discarded because in our case the phonon bath energy was >> k B T. We found that the power exponent m in RGO networks spans between the ranges defined by the values calculated by MIT and measured using NT models ( Figure 3c). Such experimental evidence i) highlight the analogy between PL transport in RGO networks and conjugated polymers thin films, and ii) suggest that both classes of materials can be studied with similar approaches. Thus, we could describe the parameter m similarly to the case of polymers, in terms of reorganization energy (R) of phonon bath upon electron transfer. 45 We generalize the NT theory defining the phonon bath energy in RGO networks = M N . Thus, the reorganization energy term describes the strength of the electron-phonon bath and can be reliably estimated as twice the relaxation energy of a polaron localized over the region of two overlapped RGO sheets, i.e. the twice the polaron binding energy: R = 1 5 + 13 • M N . 46

Correlation between ES-VRH and PL regimes.
A natural question to be addressed is whether the two CT regimes that we observed at low and high T (ES-VRH and PL) are correlated. We observed that: i) the PL exponent m (or the correlated reorganization energy λ) decreases with the increase of the localization length calculated in ES-VRH (figure S10), and ii) the measured ρ(T) data and the corresponding W(T) analysis do not show discontinuities at the transition temperature T* between ES-VRH and PL regimes.
The first experimental evidence confirms a strong analogy with polymeric semiconductors, where the reorganization energy decreases as the size of the π conjugated system is increased. 47 Moreover, the measured trend ( Figure S10) roughly agrees with the relation K ∝ 1/< computationally obtained for highly doped 1D semiconductors by Rodin et al. 48 The second experimental result implies that, at the transition temperature T*, the resistivity calculated in the VRH regime corresponds to that in PL one: $%& * = /0 * . As shown in the study of the functional forms of the ρ(T) curves, also in this case we used the W(T) method focusing the analysis of the equivalence achieved for W curves: $%& * = /0 * . The use of linearized curves made it possible to derive a general formula comparing the parameters of the two CT regimes: T*, β, m and ξ, as calculated using eqn. 1 for each i-th RGO network. Calculating all the terms as sum of logarithms (for more details see SI, Section S6.1), we obtained: Similarly to the correlation analysis discussed before in eqn. 2 and figure 3d, we should consider a set of independent equations, one for each measured i-th RGO network. Figure 3e shows the correlation plot ln K vs 8 • ln * + ln < where each point corresponds to the data of each RGO network and the red line corresponds to a linear fit from eqn. 3. The plot shows an excellent linear behavior, with slope = −8 5@T = -0.55±0.03, in good agreement with β = 1/2 (ES-VRH). The good correlation achieved is an additional proof that ES-VRH and PL regimes should be strongly correlated, even if the two regimes are based on different mechanisms and act at different temperature ranges. A similar approach related to equivalence of the resistivity values is reported in SI (figure S12 and Section S6.2).
A clarification is needed on the physical meaning of the transition temperature T*.
Experimentally, T* values range between 200 K and 20 K with the increase of the localization length ( Figure S11). This is in qualitative agreement with the model developed for disordered wires by Gornyi et al., 49 where T* is related to the disorder and it is inversely proportional to the localization length. Noteworthy, the model previously proposed by Rodin, 48 where the PL exponent m is proportional to the number of hopping events, can be seen as a special case of eqn.

Dependence of ξ on the structural and geometrical properties of the RGO nanosheets
network. In the case of single RGO sheet the transport is purely 2D, the localization length ξ corresponds approximately to the half of the size of the aromatic domains (φ): U ≈ 2K, 28 and charges "hop" from one domain to the next one. Differently, when two or more sheets are overlapped (even partially) we can distinguish two components in the transport: in-plane and out-of-plane. The case in which two adjacent sheets (i.e. belonging to the same plane) only touch the edges is statistically negligible. Thus, the sheet-to-sheet transport is mainly out-of-plane and consequently the in-plane component is related to the single sheet.
Experimentally, we varied such two components tuning independently different parameters which define the material: the oxidation degree (i.e. the corresponding intrinsic conductivity), the lateral size of the nanosheets composing the network, and the thickness of the thin film (i.e. the number of stacked nanosheets layers N layer ). In general, we observed that ξ increases with the decrease of the room temperature resistivity ( Figure 4a). First, we varied the intrinsic conductivity of each nanosheet composing the network. This could be achieved by changing the content of the sp 2 carbons on the RGO basal plane by thermal annealing, as evidenced by XPS measurements (Figure 4b). The average film thickness was kept constant with the average number of stacked sheets being N layer = 8±1. The average lateral size of the single RGO sheets was <s RGO > = 428±14 nm. Second, we varied the lateral size of the RGO sheets ( Figure 3b). We tested nanosheets having 3 different lateral sizes: S1 = 17.2±0.6 µm, S2= 428±14nm and S3= 380±7nm. For each size, three different sp 2 contents: 77% (■), 86% (■) and 96% (■) were tested for a total of 9 different size-conductivity combinations. The localization length was found to increase significantly with the aromatic content and to a lesser extent, with the nanosheet lateral size.
The results of the electrical characterizations of these films are provided in figure 4b: for low sp 2 content (≲80%), we obtain 1.3 < K < 5 nm, in agreement with the values typically observed in a single RGO sheet. 11,50 For higher aromatic contents (≳80%) the parameter K reaches a size up to 300±20 nm evidencing that charges can be delocalized over the single aromatic region on the RGO sheet. We want to underline that such result is peculiar to the macroscopic network because in the case of an assembly of a few sheets partially in contact, the measured K amounts to a few nm (see For each sample, we measured the dependence of room temperature resistivity (ρ RT ) on N layer (figure S13) showing, as expected, an opposite trend with respect to that observed on ξ in the case of RGO thin films. Such behavior is similar to that observed in thin metal films (<100 nm) where the electrical resistivity become larger as the film thickness decreases in size. In such systems this change occurs because the mean free path of charge carriers is reduced due to increased scattering effects. 51 The inter-sheet CT (out of plane-plane) is favored respect to intrasheet CT (in-plane). This result agrees with the observation of coherent commensurate electronic states at the interface between sp 2 regions, recently observed. 52 In general, all the layered graphitic materials have similar electrical properties along the out-of-plane direction; for example turbostratic and Kish, high oriented pyrolytic and natural graphites usually show metallic-or semi-metallic-like behavior. [53][54][55][56] RGO networks as "composite materials". The CT of RGO networks is governed by πconjugated regions given by the overlapping sp 2 domains connected by a network of random paths with ξ as a characteristic length. Thus, the longer the π-conjugated domains due to increased amount of aromatic content or the number of RGO layers, the greater the localization length, i.e. the lower are both the parameter m and the corresponding reorganization energy, as observed in semiconducting π-conjugated systems. 57 We developed a purely geometrical approach to describe ξ as a characteristic size and to reproduce as well the behavior of such parameter with the number of layers. RGO   Since the localization length becomes orders of magnitude larger than the average size of a sp 2 region in a single sheet: K ≫ U, we can neglect the sp 2 domain size and consider only the distance between two overlapping aromatic regions (d) corresponding to the layer-layer distance.
Thus, we can write a general form to describe the number-of-layer -dependence of the localization length as: where the first term corresponds to the contribution of the single layer (i = 1) and the latter to the following layers. n(i) is the number of steps (i.e. the number of conductive sp 2 regions involved in the path, being the stack) along the i-th layer, such parameter follows roughly a Poisson distribution because of the external bias is parallel to the film (k = l ∥ ).
Given a random path in (3D) thick film, the decrease of layers corresponds to a reduction of the number of steps. When the film thickness decreases down to a critical thickness (in our case, Heeger defined ^ as an "effective order parameter" for metal-insulator transitions (MIT) in the case of semiconducting π-conjugated polymers. 61 However, such parameter can give several issues. It is clear that this approach could not be quantitative since the arbitrary choice of the chosen temperatures. Moreover, in the case of RGO networks, ρ r roughly indicates the ratio between resistivity corresponding to VRH and PL regimes; for each RGO network the two regimes were sampled in different ways, since for each sample T* lies in different positions within the 50 K-room temperature range (see Figure 3a). For this reason, we propose to use the prefactors ρ 0,VRH and ρ 0,PL as the characteristic resistivity values of the two regimes and to define a "general effective order parameter" as #,^= #,$%& #,/0 ⁄ . Experimentally, we achieved that i) ρ 0,VRH tends to ρ 0,PL with the increase of ξ ( Figure S14), corresponding to #,^→ 1, and ii) ρ 0,PL is the same for all the RGO networks (Figure 3d).
We could thus describe the CT of RGO networks in terms of geometrical phase transition (i.e. percolation) 8 of random networks of conductive wires with mean length equal to K. The complete disorder-driven transition to metal corresponds to the case in which the critical regime is achieved for all the temperatures and the electrical resistivity no longer depends on temperature.
In such case, no VRH regime could be measured and the electrical resistivity ρ(T) trend could be We modelled the disorder-driven MIT in terms of the critical exponent formalism in continuous phase transition theory, assuming that: i) ξ corresponds to the mean cluster size and is correlated to the order parameter; ii) ρ 0,PL corresponds to the critical threshold of the transition.
Mathematically, a percolating system is described by a parameter p which controls the occupancy of bonds (or sites) in the system. At a critical value p c , the mean cluster size goes to infinity and the percolation transition takes place. As one approaches p c , the mean cluster size diverges following a power law as |* − * q | 4r , where γ is the critical exponent depending on the dimensionality. γ= 43/18 for 2D-and γ=1.8 for 3D-systems, as reported in Table 2 in ref. 58 . Such mathematical approach is general, the parameter p is a probability and shows the same role as T in thermal phase transitions -a short list of examples is reported in ref. 62 , Table 1.1 -or the filler loading in a composite in MIT, for instance.
We applied the percolation formula to follow the transition between the -: ρ 0,VRH and ρ 0,PL . In our approach, p corresponds to ρ 0,VRH and the percolation threshold p c is ρ 0,PL , defining the variable Π as: where at the percolation threshold ρ 0,VRH = ρ 0,PL (ρ 0,r = 1), the difference is zero (Π = 0) and the mean wire length is infinite (K → ∞). Conversely, ρ 0,VRH >> ρ 0,PL (ρ 0,r → ∞) and Π → ∞ when K → 0. The approach to the threshold is governed by power laws, which hold asymptotically close to critical value ρ 0,PL : K~Π 4r , as a continuous transition. This expression is analogue to the magnetic fluctuations / susceptibility x~| − q | 4r , where T c is the threshold temperature. In the case of random networks of conductive wires (i.e. spaghetti-like structure) the cluster size can be easily approximated to the wire length and thus, γ is the critical exponent describing the mean wire length. Figure 6c depicts the relation between ξ and Π in log-log scale over two orders of magnitude. The observed negative linear trend confirms the power law behavior with γ = 1.7±0.2. Such result agrees with the critical exponent of 3D systems and allows to describe RGO networks as analogous to entangled disordered systems of conductive 1D fillers with fractal dimension amounting to 2.53, close to the transition point (see Table 2 in ref. 62 ).
It is noteworthy to underline that the percolation approach model does not depend on EXPERIMENTAL SECTION Device preparation. Graphene Oxide (GO) was produced by modified Hummers methods, via the intercalation of strong oxidizing agents into Graphite, sonication in to ultrasonic bath and purified by dialysis. 65 The average size of the GO flakes was tuned by using ultrasonic bath of the solution for different amounts of time, following a previously reported procedure. 12 Thin films were produced by spin-coating the GO water suspension on a clean SiO 2 /Si substrate (2,000 rpm for 60s), the concentration of GO in water was 2.0 g/L for the high nanosheet density and 0.3 g/L for the low nanosheet density. To tune the thickness, the procedure could be repeated multiple times (see below). Reduced Graphene Oxide (RGO) thin films were produced by thermal annealing in high vacuum (10 -7 mbar) at different temperatures the GO films. In order to prevent the potential occurrence of a short circuit between the RGO film and the doped silicon, the portion of film close to the edges was mechanically removed for a width of about 1 mm ( Figure S3a). All the 28 devices showed very low leakage currents i leak (10 3 times lower than source-drain current).  Figure S1). The XPS spectra corresponding to different reduction degrees, and the details of the fitting procedure are given in our previous publications. 66 In particular, the shape of aromatic carbon peak is chosen as asymmetric following our recent work, 67 Table S4. Contact resistance was measured using the Transfer Line Method (TLM) and can be neglected ( Figure S4b).

ρ(T) equations, a brief overview
At low temperature, charge transport in disordered materials is typically occurring via charge hopping in a disorder-broadened density of states near the Fermi level, g(E F ), 2 as demonstrated in previous works (a comprehensive list is reported in the SI). In the Ohmic regime, the resistivity can be modelled by a stretched exponential behavior: where ρ 0,VRH is a prefactor and β a characteristic exponent. T 0 represents a characteristic temperature correlated to the localization length ("), the higher the first one, the lower the latter. The analytic expression depends on the model, as reported in Table S1. " is defined as the average spatial extension of the charge carrier wave function: the lower the T 0 , the larger the ". The stretching exponent β is strongly dependent on the shape of g(E F ) and can be used to infer the correct CT between any two states i and j for the inequality to be true. As a consequence of this, in the ES variable range hopping regime the characteristic exponent β of equation 1 is always β =1/2 and does not depend on the system dimensionality. 6 When CT proceeds through thermally activated hopping between nearest-neighbours (NNH model), 7 it can be modelled with an Arrhenius expression, giving β = 1.
A further model describing CT is the so-called Fluctuation-Induced Tunneling (FIT) transport. 8 The model was developed to describe the conduction mechanism in disordered systems, i.e. conducting polymers and nanocomposite materials, featuring metal pathways separated by small insulating barriers, and recently it has been considered for disordered materials with metallic and semiconducting regions. Also, FIT can be described by equation1 introducing an effective temperature 1 ∆, where ∆ is a characteristic temperature, and β = 1. The characteristic parameters of all CT models described above are summarized in Table S1. where ρ 0,PL is a prefactor and T 1 is a term that mathematically allows us to define a correct dimensional analysis and that physically corresponds to the characteristic temperature associated with the CT mechanism.
Historically, most of the studies on CT with PL model focused on studying the exponent m only, which is the parameter, describing the functional shape of the measured . Conversely, only fewer works investigated the characteristic temperature T 1 , although the latter can provide direct information on the transport mechanism. In the specific case materials of conjugated materials, different models are required to describe different ranges of m. In the case of electronic transport, for 1/3 < m < 1 the model describes disordered systems with large localization length, and T 1 is equal to the Fermi temperature. [11][12] When m ≥ 1, (integer values) the transport is characterized by a single-/multi-phonon assisted process and m corresponds to the number of phonons involved in the process with characteristic energy = k B T 1 . For instance, in the case of amorphous carbon films m ranges between 15 and 17. 13 In the case of strong electron-lattice coupling, multiphonon processes are described in terms of polaronic transport. In such case, m can range between 0.3 and 11, 14 as observed in disordered semiconducting polymers and described by different models (i.e. Nuclear tunnelling, Miller-Abrahams and Marcus theory).
A further aspect related to the wide range of cases described by the equation 2 is that also the stretched exponential form of equation S2.1 can be approximated to a power law in a wide range of parameters. For instance, by numerical simulation Rodin et al. 15 showed that equation S2.1 collapsed to equation S2.2 for high doping while Gornyi et al. 16 modelled the transition between the two regimes at increasing disorder.

Charge transport properties of RGO thin films reported in literature
The published results lie in a wide range, as summarized in Table S2. This issue can be ascribed to the complexity of the chemical structure of GO, whose properties strongly depend on the mechanisms of reduction, and the lack of a standard procedure to perform quantitative analysis of the R(T) measurements. Moreover, the use of samples with different sizes and geometries, and as well as different temperature ranges, makes the comparison of results complicated.
Regarding the quantitative analysis, different approaches are used to manage of the acquired data.
Only few works report a systematic W(

FIT transport. Comments on W(T) analysis
Taking into account the expression reported in Equation S2.2, the corresponding reduced activation energy of FIT model assumes the form: The use of T eff allows to use a clear and compact analytic form, however the experimental parameter is the temperature T and in particular an analytical expansion in term of lnT has to be provided. The logarithmic expressions can be written as: In the limit T>>∆, we can use the Maclaurin series and equation S3.6 becomes gh 1 ≅ gh ∆ .
Thus, the Equation S3.3 assumes the form: It is noteworthy to remind that the curve described by the Equation S3.5 is a concave function and describes a line in the case ∆ = 0 (i.e. Arrhenius-like case), as reported in figure S6. Figure S6. Plot of W(T eff ) in lnW vs lnT space. FIT regime is described by a concave curve (bold curve) collapsing to a line (dot line) when ∆ = 0, corresponding to a VRH-like case with β = 1. Figure S7. Device #1-#26. RGO thin-films.     Table S3. Comprehensive results of the analysis shown in Figure 1 for all the reported devices. The two regimes are separated at the characteristic temperature T*, as directly obtained by W(T) plots and using equation 3. The carbon sp 2 fraction was estimated by XPS analysis, the number of effective layers (N eff ) by AFM. †The relative error amounts to 5% in log-space (log T*). ‡ W(T) of device #25 shows a high-noise level because of the low reduction degree and the corresponding higher resistivity. For this reason, the parameter β was fixed to 0.5 in order to evaluate the compatibility of the data with the ES-VRH.  It is a phenomenological evidence that there is abrupt variation in ρ(T) and W(T) at the transition temperature. For this reason, we assume the equivalence: i * d i ef * 6 (S6.1a)

Summary of experimental parameters of all the prepared devices
The equivalence is conserved using the logarithms: ln i * d ln i ef * 6 (S6.1b) Thus, by taking into account the analytic form of W of the two regimes we obtain: −# ln * ln # # ln ln q (S6.1c) Using the formula u • ", as reported in Table 1, we obtain the form shown in equation 3.
As performed in Section 6.1, we assume the equivalence of the resistivity functions: Using the formula u • " and calculating the logarithms of the two terms we obtain: −# • ln * − # • ln " # • ln u ln −q • ln * q • ln C − l (S6.2c) Using the eq. S6.1d, the first term of eq. S6.2c assumes the form: ln q − ln #. Thus, we calculate the exponential forms obtaining the equivalence: q • @ln H * − C B l (S6.2d) In the case q 0 we obtain a trivial solution as: l 0, corresponding to , ,ef

RGO films as an amorphous polymer network material
Using the approach to describe the RGO thin films in term of composite materials, the localization length is defined as the product of the number of distinct overlapped sp2 regions (n tot ), called nodes, where the cumulative function of the Poisson distribution is the incomplete gamma function Γ , " . 28 In the case of 3D system (thick RGO film) the localization length amounts in the order of magnitude of microns. Taking into account that the interlayer distance is lower than 1 nm, the number of nodes is ca 1,000. In the limit n ≈ 1,000, the Poisson distribution can be approximated to the Gaussian one. In this case, taking into account that where the comulative distribution of the  Taking into account that " depends exponentially on disorder, the equation S6.4d assumes the form: