Introduction

The unique electrical properties of graphene, such as high carrier mobility, µ > 104 cm2/Vs, at room temperature1,2,3, offer significant advantages for applications ranging from fast electronics to touch screens and ultrasensitive photon detection4,5,6. However, the emergence of graphene electronics on the market is limited by the absence of cost-effective large-scale production of high-quality graphene with reproducible electronic properties. The best results have been achieved in exfoliated suspended single-layer graphene (SLG) samples a few micrometres across, in which µ is limited only by the scattering of charge carriers (electrons and/or holes) by intrinsic phonons7,8. Epitaxial growth of graphene by chemical vapour deposition (CVD)9,10 or SiC-surface growth11,12 methods provides cost-effective growth of large (>10 mm) SLG layers. However, the mobility of suspended sheets of graphene is markedly different to that of graphene deposited on a substrate13,14. The presence of charged impurities near the graphene significantly reduces µ due to the associated long-range Coulomb-scattering centres7,15. If impurities are present, they often ionise and form charge-scattering centres, which deflect the trajectories of electrons and holes in the two-dimensional (2D) layer, thereby degrading the mobility. Often scattering by charged impurities has the dominant effect on the transport properties of graphene and effects due to ballistic transport are negligible16. Charged impurities are introduced in the substrate and/or in the SLG-capping layer during the device processing, e.g., created by the diffusion of metallic ions present in the solvents or etching solutions used17. Recent studies indicate that resonant impurities18 and neutral impurities and defects limit graphene mobility15,19,20 and can dominate at carrier densities away from neutrality point, resulting in a sub-linear dependence of the conductivity on carrier concentration.

Various theoretical models have been proposed to explain the effect of impurities on the carrier mobility in SLG. It is commonly accepted that µ is inversely proportional to the charged impurity density, nimp, but is independent of carrier density n21,22. Scattering by charge-neutral point defects can also affect µ, which in this case is inversely proportional to the carrier density23,24, making it the dominant scattering mechanism at large carrier densities. To date, qualitative and semi-quantitative models have been developed to describe the conductivity minimum; however, the physical mechanism for the minimum conductivity is not yet fully understood21,25. Owing to the 2D nature of SLG, the mobility is sensitive to the surrounding environment, in particular to the presence and position of the charged impurities. However, there is still limited understanding of the effect of the stand-off distance, d, of the impurities from the graphene plane on the carrier mobility and other transport parameters. Additional complications arise in graphene-based heterostructures where SLG is sandwiched between two other materials with the same or different dimensionality (three-dimensional (3D) bulk, 2D layers and/or 0D quantum dots (QDs)2,26,27,28).

Here we report a theoretical and experimental study of the effects of charged impurities and optical phonons on the charge transport properties of graphene. We develop a conductivity model for graphene based on the convolution of carrier density, which accounts for temporal and spatial fluctuations of the carrier density. This model accurately reproduces the experimental data and provides a robust, simple way to model the graphene conductivity as a function of the impurity density and carrier concentration. We show that the experimentally measured parameter δn, which is the full width at half maxima of the ρ(Vg), enables us to fit the whole conductivity curve σ(Vg) and to determine the exact shape of the conductivity minimum plateau for a wide range of graphene devices. Hence, several properties of electron transport in graphene can be determined using δn including the mobility and concentration dependence on Vg. This model is verified by numerical k-space simulations of carrier transport. We use the Discontinuous Galerkin (DG) technique to numerically solve the Boltzmann transport equation and investigate the effect of charged impurity scattering on the electron/hole mobility. We demonstrate that such processes give rise to universal mobility characteristics over a wide range of carrier densities and in the presence of multiple sources of scattering. The calculations are supported by experimental results obtained on both pristine and surface-decorated graphene devices, which have the following structures Si/SiO2/Graphene and Si/SiO2/Graphene/2D(0D), respectively. Our investigations show how these scattering processes give rise to mobility characteristics, which are universal over a wide range of graphene devices, and thus potential sources of scattering. Our results enable new understanding, based on first-principles calculations, of the link between different transport parameters, which are of fundamental and applied interest.

Results and discussion

Modelling the transport properties of graphene

In our work, we consider graphene sheets with charged impurities at a distance, dimp, from the plane of the graphene and optical phonons with energy ħω (Fig. 1a). We model the effect of impurities and optical phonon scattering on the following transport properties of graphene: the carrier concentration (n), mobility (µ), conductivity (σ) and resistivity (ρ) at the Dirac point (σmin and ρmax, respectively). The graphene conductivity in the vicinity of the Dirac point can be strongly affected by a number of different phenomena besides impurity scattering, including ballistic transport effects29, quantum capacitance7,30 and temperature31,32. As a result, the device conductivity and carrier density are non-zero even when the Fermi energy, εF, is at the Dirac (charge neutrality) point. Therefore, a simple model for the Drude conductivity:

$$\sigma _{\mathrm{D}} = e\mu n_{\mathrm{c}},$$
(1)

with a constant mobility, µ, is not applicable for small gate voltages, Vg, assuming the classical capacitance model for the graphene’s carrier number density

$$n_{\mathrm{c}} = \left| {CV_{\mathrm{g}}/e + n_0} \right|,$$
(2)

where n0 = nc(Vg = 0) is the sheet density of the graphene doping (Supplementary Note 1). Spatial fluctuations of the local electrostatic potential in the graphene layer and the presence of electron and hole ‘puddles’ (Fig. 1b) are thought to explain the non-zero conductivity and resistance (σmin and ρmax in Fig. 1c, d) observed at the Dirac point21,33. Electrons and holes play equal roles in determining the graphene conductivity with no scattering at the borders between the n- and p-type graphene areas due to the Klein paradox34.

Fig. 1: Convolution model.
figure 1

a Schematic diagram showing the impurity position with respect to the graphene sheet and a pictorial representation of spatio-temporal (angular frequency, ω) phonon oscillations over one unit cell. b Spatial inhomogeneity in the impurity potential and the resulting broadening of the energy distribution of electrons and holes near the Fermi level (Dirac cones). Conductivity (c) and resistivity (d) of graphene plotted vs. gate voltage. Field effect mobility, μFE, is estimated from the gradient of σ(Vg) dependence. Circles show experimental data, blue lines show linear fits, using Eq. (1), dashed pink curves are fits to the data found by including a short-range scattering resistivity, and red curves result from the convolution fit.

We now consider spatial fluctuations in the carrier number density and their dependence on the gate voltage. Combined with the Drude conductivity, this model accurately describes the shape of the measured ρ(Vg) curve (Fig. 1d).

To first approximation, the carrier number density at a given gate voltage, n(Vg), can be modelled as the moving average (convolution) of nc (Eq. (2)) over a window of width δn, which is the characteristic amplitude of the carrier density fluctuations in the graphene layer. This is equivalent to the convolution of the nc(Vg) function with a box function f(nc) of width δn (see Supplementary Video (convolution2.avi) and detailed description in Supplementary Note 1), which gives \(n = n_{\mathrm{c}} \ast\, f(n_{\mathrm{c}})\), where \(f(n_{\mathrm{c}}) = \left\{ \begin{array}{*{20}{c}} \frac{1}{{\delta n}} & {\mathrm{for}} - \frac{{\delta n}}{2}\,<\,n_{\mathrm{c}} < + \frac{{\delta n}}{2} \hfill\\ 0 & {\mathrm{for}} - \frac{{\delta n}}{2}\,> \, n_{\mathrm{c}}\;{\mathrm{or}}\;n_{\mathrm{c}} > + \frac{{\delta n}}{2} \end{array} \right.\). Using the linear nc(Vg) dependence in Eq. (2) gives

$$n\left( {V_{\mathrm{g}}} \right) = \left\{ \begin{array}{*{20}{c}} \frac{{\delta n}}{4} + \frac{{n_{\mathrm{c}}^2\left( {V_{\mathrm{g}}} \right)}}{{\delta n}} & {\mathrm{for}}\;n_{\mathrm{c}}(V_{\mathrm{g}})\,<\,\frac{{\delta n}}{2} \\ n_{\mathrm{c}}\left( {V_{\mathrm{g}}} \right) & {\mathrm{for}}\;n_{\mathrm{c}}(V_{\mathrm{g}})\,> \,\frac{{\delta n}}{2} \end{array} \right..$$
(3)

This expression for n(Vg) is equal to the constant-capacitance model for gate voltages where \(n_{\mathrm{c}}(V_{\mathrm{g}})\,> \,\frac{{\delta n}}{2}\) and has a parabolic form for gate voltages close to the Dirac point, where \(n_{\mathrm{c}}(V_{\mathrm{g}})\,<\,\frac{{\delta n}}{2}\). Using Eq. (3), we can determine the conductivity as

$$\sigma _{\mathrm{f}}(V_{\mathrm{g}}) = e\mu n(V_{\mathrm{g}}).$$
(4)

From this expression, it can be shown that δn equals the full width half maximum of the peak resistivity around the Dirac point. In addition, when nc = 0, δn = 4nNP, where nNP = σmin/, is the residual carrier density at the Dirac (neutrality) point. Our model thus provides a simple expression which enables us to extend the linear conductivity model8 to values of n close to the Dirac point, which was not possible with previous models21. Thereby, the full observed σ(Vg) and ρ(Vg) dependences can be accurately reproduced (Fig. 1c, d).

To model δn(Vg) and μ(Vg) curves in more detail, we consider semi-classical k-space simulations of the electron and hole dynamics. Graphene has a linear energy-wavevector dispersion relation, with electron energy in graphene ε = ±ħv|k|, where v = 106 ms−1 is the speed of electrons in graphene, for the electrons and holes at the two inequivalent valleys, K and K′, in reciprocal space. Charge carriers undergo diffusive scattering transport, which we describe using a semi-classical Boltzmann transport approach. The influence of perturbations, such as impurities and phonons on the scattering of electrons is calculated using the Fermi golden rule for transition rates between states. The electrons are initially assumed to obey a Fermi–Dirac distribution, f0(k). Inter-band transitions are neglected such that the valence band is assumed to be full throughout the time evolution when the gate voltage is positive, i.e., when the chemical potential lies within the conduction band. We assume full ionisation of all the impurities and their distribution to be independent on gate voltage. In the high-gate voltage regime, this assumption is confirmed by the linear dependence of n(Vg), where the value of n is determined using the equation for the field effect capacitance and verified using the Hall effect measurements.

The spatially homogeneous Boltzmann transport equation,

$$\frac{{\partial f(t,\,{{\mathbf{{k}}}})}}{{\partial t}} = - \frac{e}{\hbar }{{\mathbf{{E}}}} \cdot \nabla _{\mathbf{{k}}}f\left( {t,\,{{\mathbf{k}}}} \right) + \left( {\frac{{\partial f\left( {t,\,{{\mathbf{k}}}} \right)}}{{\partial t}}} \right)_{\mathrm{coll}},$$
(5)

describes the evolution of the occupancy, f(t, k), of state k at time t. The first term on the right-hand side of Eq. (5) describes the acceleration of electrons under an applied electric field, E, and the collision term is given by the detailed balance equation,

$$\left( {\frac{{\partial f\left( {t,\,{{\mathbf{k}}}} \right)}}{{\partial t}}} \right)_{\mathrm{coll}} = \, \mathop {\sum }\limits_{{{\mathbf{k}}}^{\prime} } \left[ S_{{{\mathbf{k}}}^{\prime} \to {{\mathbf{k}}}}f(t,\,{{\mathbf{k}}}^{\prime} )\left( {1 - f(t,\,{{\mathbf{k}}})} \right)\right.\\ \, - \left. S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }f(t,\,{{\mathbf{k}}})\left( {1 - f(t,\,{{\mathbf{k}}}^{\prime} )} \right) \right],$$
(6)

where Skk′ is the transition rate of carriers from a state of crystal momentum ħk to a new state with momentum ħk′. Equation (6) represents the collision integral. In the particular case of elastic scattering, Skk = Sk′→k and the products of the two distributions cancel. We solve Eq. (5) using the DG approach35 (for detailed solution, see Supplementary Note 2) for the steady-state distribution function, f(k). We then determine the mobility, μ, for an applied electric field, E, from the drift velocity:

$${\boldsymbol{v}}_{{\mathbf{d}}} = \frac{1}{{n\pi ^2\hbar }}{\int} {f\left( {{\mathbf{k}}} \right)\nabla _{\mathbf{{k}}}\varepsilon \left( {{\mathbf{k}}} \right){\mathrm{d}}{{\mathbf{k}}} = \mu {{\mathbf{E}}}} .$$
(7)

Alternatively, to find approximate analytical solutions to Eq. (5), we can assume a small shift in the initial distribution function, f0(k), proportional to the ensemble momentum relaxation time, τm. This results in the linearised Boltzmann (LB) approximation36 for the mobility, which, at zero-temperature, is related to the relaxation time at the chemical potential, εF, via

$$\mu = e\tau _{\mathrm{m}}\left( {\varepsilon _{\mathrm{F}}} \right)v^2/\varepsilon _{\mathrm{F}}.$$
(8)

We calculate the momentum relaxation time, τm(k), as the sum over all possible transition rates, Skk′, modified by the deflection angle, θk,k′, between the incoming and outgoing vectors:

$$\frac{1}{{\tau _m\left( {{\mathbf{{k}}}} \right)}} = \, \mathop {\sum }\limits_{{{\mathbf{k}}}^{\prime} } S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }( {1 - {\mathrm{cos}}\theta _{{{\mathbf{k}}},{{\mathbf{k}}}^{\prime} }} )\\ \approx \,\frac{A}{{\left( {2\pi } \right)^2}}\int S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }( {1 - {\mathrm{cos}}\theta _{{{\mathbf{k}}},{{\mathbf{k}}}^{\prime} }} ){\mathrm{d}}{{\mathbf{k}}}^{\prime} ,$$
(9)

where A is the area of graphene unit cell.

The effect of screening by the electron and hole gases is included by introducing a random phase approximation for the dielectric screening function37,38

$$\epsilon_{\mathrm{{sc}}} = 1 - \tilde v_{\mathrm{{2D}}}{\Pi}\left( q \right),$$
(10)

where 2D = e2/2ϵ0ϵrq is the unscreened Coulomb potential in Fourier space, Π(q) is the static polarisation function and the reciprocal space variable q = |k′ − k|. As the conduction band distribution function changes throughout the simulation, the screening function should be carefully considered. However, the valence band distribution is constant throughout as the band is assumed to be full. The polarisation function for screening by a full valence band is38

$${\Pi}_{\mathrm{val}} = - \frac{q}{{4\hbar v}}.$$
(11)

For the conduction band, at T = 0, maximum screening occurs in the Thomas–Fermi limit, q → 0, and is given by

$${\Pi}_{\mathrm{TF}} = {\int} {D\left( \varepsilon \right)\frac{{\partial f}}{{\partial \varepsilon }}{\mathrm{d}}\varepsilon = - \frac{2}{{\pi \left( {\hbar v} \right)^2}}{\int} {f{\mathrm{d}}\varepsilon = - \frac{2}{{\pi \left( {\hbar v} \right)^2}}\varepsilon _{\mathrm{F}},} }$$
(12)

where D(ε) = 2ε/π(ħv)2 is the density of states. Throughout the simulation, the integral in Eq. (12) does not change, due to conservation of charge. As both Eqs. (11) and (12) are independent of the evolution of the distribution function, we define a time-independent two-regime screening function, where Thomas–Fermi screening is assumed for low energy scattering and the valence electron screening is assumed for high-energy electrons21, i.e., we set

$$\epsilon _{\mathrm{sc}} = \left\{ {\begin{array}{*{20}{c}} {1 + \frac{{q_{\mathrm{s}}}}{q}\quad {\mathrm{for}}\;q \le \frac{8}{\pi }k_{\mathrm{F}}} \\ {1 + \frac{{\pi r_{\mathrm{s}}}}{2}\quad {\mathrm{for}}\;q\,> \,\frac{8}{\pi }k_{\mathrm{F}}} \end{array}} \right.,$$
(13)

where rs = e2/(4πϵ0ϵrħv) and qs = 4kFrs. For graphene on SiO2, we take ϵr ≈ 2.4522. A Coulombic scattering potential is assumed for charged impurities near the graphene plane:

$$U(r) = \frac{{e^2}}{{4\pi \epsilon _{\mathrm{0}} \epsilon _{\mathrm{r}}\sqrt {r^2 + d_{\mathrm{imp}}^2} }},$$
(14)

where dimp is the distance of the impurities from the graphene plane. It is noteworthy that in this model we consider randomly distributed impurities and disregard any possible spatial correlation of charges below and above the graphene plane39,40. Then, the transition rate is

$$S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }^{\mathrm{imp}} = n_{\mathrm{0}}\frac{\pi }{{A\hbar }}\left| {\frac{{2\pi e^2e^{ - qd_{\mathrm{imp}}}}}{{\kappa q \epsilon _{\mathrm{sc}}\left( q \right)}}} \right|^2\left( {1 + {\mathrm{cos}}\theta _{{{\mathbf{k}}},{{\mathbf{k}}}^{\prime} }} \right)\delta \left( {\varepsilon _{{{\mathbf{k}}}^{\prime} } - \varepsilon _{{\mathbf{k}}}} \right),$$
(15)

where q = 2k sin(θk,k′/2) for elastic scattering. Defects that perturb the band structure over a small spatial area are characterised by the short-range scattering potential U(r) = U0H(R − r), where H is the Heaviside step function and R gives the spatial extent of the perturbation. This potential represents any charge-neutral point defects within the lattice. The rate of carrier scattering transitions due to such defects is

$$S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }^{\mathrm{sr}} = n_{\mathrm{sr}}\frac{\pi }{{A\hbar }}\left| {\frac{{A_{\mathrm{sr}}U_0}}{{ \epsilon _{\mathrm{sc}}\left( q \right)}}} \right|^2\left( {1 + {\mathrm{cos}}\theta _{{{\mathbf{k}}},{{\mathbf{k}}}^{\prime} }} \right)\delta \left( {\varepsilon _{{{\mathbf{k}}}\prime } - \varepsilon _{{\mathbf{k}}}} \right),$$
(16)

where Asr = πR2 is the effective cross-section of defects with an areal density nsr.

In our calculations, we assume a low temperature and a phonon occupation of N ≈ 0 (kBTħω) and perform our numerical calculations using an initial Fermi-distribution of low finite temperature, T = 20 K, to avoid discontinuities over the discretised k-space. Therefore, one might not expect phonons to have a significant effect on the transport properties compared to that of scattering by impurities7. However, for low carrier densities, we find carriers can be accelerated to high energies (~100 meV) resulting in a ‘hot electron’ distribution (Joule heating), as was observed previously, e.g., in metals41. In this case, inelastic optical phonon scattering becomes important in relaxing the energy of the carriers. Hot electron phenomenon is a particularly important consideration for transport in graphene due to weak electron-acoustic phonon scattering and relatively high optical phonon energies. Our calculations show that the hot electron effect is significant even for electric fields as low as ~100 V/m, comparable to commonly used experimental values (for steady-state characteristics of the Boltzmann equation, see Supplementary Note 3). We use the optical phonon scattering rates calculated using density functional theory in refs. 42,43. Near the Γ− points of the reciprocal lattice, the energy and coupling strength of both transverse and longitudinal optical phonons are reported to be ħωO ≈ 165 meV and βO ≈ 10 eV/Å, respectively44. Therefore, we can combine the transition rates of the two modes to obtain a single overall optical scattering rate

$$S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }^{\mathrm{O}} = \frac{{2\pi \beta _{\mathrm{O}}^2}}{{A\rho _{\mathrm{m}}\omega _{\mathrm{O}}}}\delta \left( {\varepsilon _{{{\mathbf{k}}}^{\prime} } - \varepsilon _{{\mathbf{k}}} + \hbar \omega _{\mathrm{O}}} \right).$$
(17)

Phonons at the K-points cause inter-valley scattering at a rate

$$S_{{{\mathbf{k}}} \to {{\mathbf{k}}}^{\prime} }^{\mathrm{K}} = \frac{{2\pi \beta K^2}}{{A\rho _{\mathrm{m}}\omega _{\mathrm{K}}}}\left( {1 - {\mathrm{cos}}\theta _{{{\mathbf{k}}},{{\mathbf{k}}}^{\prime} }} \right)\delta \left( {\varepsilon _{{{\mathbf{k}}}^{\prime} } - \varepsilon _{{\mathbf{k}}} + \hbar \omega _{\mathrm{K}}} \right),$$
(18)

where ħωK ≈ 124 meV and βK ≈ 3.5 eV/Å44, and ρm = 7.6 × 10−7 kg/m2 is the mass density of graphene.

We consider a residual charge density of electron-hole puddles at the Dirac point, due to inhomogeneity in the impurity-induced potential (Fig. 1b), which limits the minimum conductivity. To calculate the residual charge, and thus the minimum chemical potential in our calculations, we use Eq. (8) in ref. 21, derived assuming a random distribution of impurities. Here we assume that the transition from the residual charge-dominated minimum carrier concentration to the linearly Vg-dependent concentration occurs when the gate-induced charged density, n(Vg − V0), is equal to the residual charge density, nNP, where V0 is the position of the Dirac point.

For all numerical simulations, we apply an electric field, E = 104 V/m (0.1 V drop across a 10 µm-long SLG), corresponding to a regime of low-field mobility, where μ is independent of the applied electric field strength (for details of numerical simulations, see Supplementary Note 3). For comparison, we also calculate the mobility using the LB formalism, Eq. (8), with the scattering time calculated using Eq. (9), in which the integrals are evaluated numerically. The LB method gives exact solutions when the electric field is sufficiently small and the hot electron effects are negligible. However, for small densities, carriers can be accelerated to high energies resulting in a hot carrier distribution which is far from thermal equilibrium. In this case, the LB approximation diverges from the accurate numerical solution provided by the DG approach (for details of DG approach, see Supplementary Note 3).

Figure 2a shows the calculated dependence of μ on n for n > nNP. The low carrier mobility μ(n ≈ nNP) corresponds to the regime where the chemical potential is near the Dirac point. With increasing n, we observe an initial increase of μ. This is followed by a peak and a monotonic decrease of μ at large n. This dependence arises from the competition between scattering by long-range Coulombic impurities and short-range defects. Short-range defect scattering is found to be dominant at large n, as expected from comparison of the momentum relaxation time for short-range defects, τsr ~ n−1/2, calculated using the Born approximation, and long-range impurities, τimp ~ n1/2. Beyond the Born approximation, for sufficiently strong defect scattering, the exponent of n in the momentum relaxation time, τsr, can increase towards that of long-range impurity scattering33. The dependence of mobility on carrier concentration, μ(n), is affected by the density of impurities, n0, and by their distance from the graphene, dimp. Hence, both δn and the low carrier mobility, μ(n ≈ nNP), depend on n0 and dimp. As shown in Fig. 2b, the mobility increases as dimp is increased. Furthermore, for low impurity densities, and thus small residual charge densities, the mobility given by the DG simulations differs from that obtained from the LB calculations, whereas the two methods give μ values that converge at higher impurity densities.

Fig. 2: Boltzmann equation results for the mobility characteristics of doped graphene.
figure 2

a Mobility, μ, calculated using our Discontinuous Galerkin (DG) model as a function of gate-induced carrier density, n, beyond the minimum carrier density, nNP (vertical dashed line), for varying short-range scattering strengths. The impurity density is n0 = 0.6 × 1016 m−2. b Mobility, μ, calculated as a function of impurity density, n0, for several stand-off distances, dimp. Dashed curves are calculated using the linearised Boltzmann (LB) approximation, using Eq. (8), the solid green curve shows the constant-capacitance approximation from ref. 21, where α = 20e/h. c Final distribution of electron momentum calculated for a small impurity density, n0 = 0.025 × 1016 m−2, at a distance dimp = 1 nm and carrier density, n = nNP = 0.016 × 1016 m−2 found using the DG simulation. Orange dashed circle represents the K-phonon level and the orange dashed-dotted circle represents the Γ-phonon level. d Cross-section of c for ky = 0 (solid blue line) compared to the result of Monte Carlo simulations. The Fermi wavenumber, kF, marks the bounding limit of the initial low temperature distribution.

Both the LB model and the DG simulations assume that initially, at t = 0, the electron gas is in thermal equilibrium and obeys the Fermi–Dirac distribution. Equation (8) assumes a linear shift in the momentum of the ensemble, proportional to the ensemble relaxation rate, τ(εF), whereas the DG simulations include the time evolution of the momentum distribution, described by the full Boltzmann transport equation, Eq. (5). Therefore, the discrepancy between the two methods seen at low impurity densities can be understood by consideration of the steady-state distribution functions. Figure 2c, d shows the final distribution of electrons obtained for a small impurity density, using the DG simulation. We obtain similar results by Monte Carlo simulations35,45 (Monte Carlo simulations are described in part (A) of Supplementary Note 2) as shown in Fig. 2d. In both the DG and Monte Carlo simulations, with increasing t we observe continuous spreading of the electron distribution in k-space, until the occupied k-values become limited by inelastic phonon scattering (see also Supplementary Note 3). Hence, a hot electron regime is realised, which is not captured within the LB approximation.

The effect of dimp on the electrical properties of graphene is summarised in Fig. 3. Our calculations demonstrate that the linewidth δn of the ρ(Vg) curve broadens with decreasing dimp (Fig. 3a). Combining the results of Figs. 2b and 3a, the mobility decreases with increasing δn (Fig. 3b), with the broadening of δn being larger for smaller dimp at a given value of mobility. At small values of δn, and hence small impurity densities n0, we observe discrepancy between the DG and LB calculations of μ (as shown in Figs. 2b and 3b). Despite the discrepancy in μ(n0) and μn), we find that both methods yield a similar δn(n0) profile. We now compare our calculations to experimental measurements.

Fig. 3: Mobility and conductivity characteristics of graphene devices.
figure 3

a Dependence of the FWHM of resistivity, δn, on impurity density, n0, for different stand-off distances dimp = 0.5 nm (blue curves), 1.0 nm (red curves), 2.0 nm (black curves). b Calculated mobility, μ, vs. the full width at half maximum (FWHM), δn, curves (dashed and solid curves) compared to data from pristine graphene samples grown by chemical vapour deposition (CVD) (filled circles) or exfoliated (filled triangles). The dashed lines are obtained from the linearised Boltzmann (LB) approximation, using Eq. (8), the solid lines are obtained using the Discontinuous Galerkin (DG) simulations.

Universal mobility characteristics

We apply our analysis to experimental results reported previously for >20 devices fabricated using exfoliated and CVD-grown graphene deposited on Si/SiO2 substrate. We use both pristine graphene devices and graphene heterostructures incorporating 2D layers (InSe, hBN) or 0D nanostructures (colloidal QDs, inorganic perovskites)46,47,48 (Fig. 3a). In these devices, impurities at a distance, dimp, from the 2D plane of graphene act as scattering centres. We fit the measured σ(Vg) dependencies and determine values of μ and δn (fit of σ(Vg) is described in Supplementary Note 1 and phenomenological fit of experimental data is in Supplementary Note 4). As shown in Fig. 3b, the mobility increases with decreasing δn. The experimental values measured in pristine graphene devices are in good agreement with the results of our DG simulations with dimp = 2 nm. Interestingly, our model provides good fit for high-mobility exfoliated graphene, where other scattering mechanisms could play a significant role. Since the convolution model is based on experimentally determined value of δn, it accounts for all different scattering mechanisms (for universality of analytical convolution model, see Supplementary Note 5).

We note that our fit (Fig. 3b) uses δn calculated from the full width at half maximum of the σ(Vg) curve rather than from the value of n0 extracted from the gate voltage at which σ(Vg) = σmax. By using Eq. (4) and assuming the universal minimum conductivity for pristine graphene as σmin ≈ 4e2/h29 and constant mobility (with respect to carrier density), we obtain a simple inverse power law for the dependence of μ on δn.

$$\mu = \frac{{4\sigma _{{\mathrm{min}}}}}{{e\delta n}} = \frac{{16e}}{{h\delta n}}.$$
(19)

Equation (19) includes one measured parameter δn, which simplifies the data analysis, as demonstrated on the experimental data from a wide range of devices (experimental results are included in Supplementary Note 4). Overall, our model, which considers the effect of impurity scattering to be dominant on mobility, describes well all examined types of graphene: high-mobility exfoliated graphene and low-mobility CVD-grown graphene. We stress that, remarkably, even in devices where other transport mechanisms are important, e.g., ballistic transport in high-purity exfoliated graphene, their electronic properties can be determined using the measured value of δn.

Recently, the decoration of graphene devices with other low-dimensional materials, such as 0D (colloidal PbS quantum dots46 or CsPbI3 perovskite47) and 2D (InSe flakes)48 materials has been used to functionalise these devices, e.g., for photon sensing5,47,48. The properties of the graphene heterostructures are greatly affected by both the unintentional presence of charge impurities in the vicinity of graphene (as described above by dimp) and those deliberately introduced by the top layer (dtop) in graphene heterostructures (Fig. 4a), which we model as a distribution of impurities at an effective distance, deff. We note that in surface-decorated graphene devices, the distance between the graphene plane and the top layer can be controlled, e.g., by introducing a dielectric layer such as hBN, thus providing a tool for tailoring the electrical properties. The relationship between mobility and the gate-voltage offset is μ 1/n0 for most pristine devices21. However, for devices with high densities of correlated unipolar charges39,40 or uncorrelated bipolar charges49, spatial correlation between charges must be considered. This is particularly important when the dopants are mobile and able to adopt low energy, correlated configurations. Such effects were recently demonstrated for quantum dot-decorated graphene and validated using Monte Carlo simulations40,49.

Fig. 4: Mobility characteristics of graphene devices and heterostructures.
figure 4

a Schematic diagram showing the position of impurities and surface-charges, due to 0D structures (e.g., perovskites) and 2D structures (e.g., InSe), with respect to the graphene sheet. b Relationship between mobility, μ, and the full width at half maximum (FWHM), δn, obtained using the Discontinuous Galerkin (DG) simulations, taking deff = 1.0 nm. In surface-decorated devices, the effective impurity distance, deff, describes combined effect of charges below (dimp) and above (dtop) graphene layer compared to data from multiple modified graphene samples (data points for each sample type are labelled as shown in the inset legend).

Despite the different μ(n0) characteristics of decorated and pristine graphene, remarkably, we find that both types of devices exhibit the universal scaling behaviour shown in Fig. 4b. Different surface-decorated devices follow a common trend observed in pristine graphene. In particular, the experimental results for the InSe, perovskite and PbS decorated SLG are best fitted by DG calculations when deff = 1 nm. Therefore, we find that the relationship between μ and δn is consistent throughout all of the devices, as can be expected from the analytical expression given in Eq. (19), with modifications to only the effective distance of the impurities. Flexibility to modify composition and/or geometry of a heterostructure offers opportunities to tune the distribution and stand-off distance of ionised impurities, hence changing deff and providing a tool to control transport properties of these devices. We note, that our model is valid for all devices where the position of ionised impurities is not affected by Vg. In rare cases, at high Vg regime, the ionisation of donor impurities can be affected by applied gate voltage (e.g., see ref. 50) and the corresponding change of dimp would need to be accounted for.

Our model links together three key transport parameters of SLG devices: μ, n0 and δn, where δn can also be used to calculate ρmax and σ(Vg) (for phenomenological equations for graphene transport parameters, see Supplementary Note 6). Remarkably, this model can be used to extrapolate the whole R(Vg) dependence from a single R(Vg) =Rmax measurement and for a wide range of different graphene devices (see Fig. 1c, d). Our approach is based on experimental value of δn, which accounts for presence of scattering centres, but does not distinguish their nature. We envisage that majority of ionised scattering centres present in our devices originates from substrate impurities and from impurities introduced from top layer (2D or 0D). The effects of other types of ionised impurities (substitutional doping, functional groups, etc.51) merits further studies.

Of particular interest is the applicability of our model to a wide variety of different graphene types and to different device structures and geometries. Consequently, the model has the potential to both predict and explain the observed behaviour of newly emerging device concepts and graphene types. Recently, the need for upscaling of graphene growth and device manufacturing has led to significant research focus on Molecular Beam Epitaxial growth52, liquid exfoliation of 2D materials53 and additive manufacturing (3D printing) of graphene devices54,55. Our preliminary results indicate that our model can be optimised and expanded to explain and predict the properties of 3D printed graphene devices, by accounting for flake-to-flake hopping of charges56.

Conclusions

We have developed a universal analytical convolution model of electron transport in graphene and graphene heterostructures, supported by numerical time-dependent analysis of the charge carrier distributions. Our model includes the effects of impurities and optical phonons on the observed charge transport properties of graphene devices. We find that the properties of a wide range of devices, from high-quality graphene with low carrier density to graphene-based heterostructures, exhibit universal behaviour that can be accurately described with this model. Our calculations combine multiple parameters that affect charge transport in graphene and facilitate the design, accurate ab initio prediction of key transport parameters and analysis of future electronic and optoelectronic devices based on 2D materials. Furthermore, our results may be generalised to predict and improve the electrical behaviour of 2D multilayers made by 3D printing or from reduced graphene oxide, which are promising candidates for the scalable high-yield manufacture of large-area optoelectronic devices that harness the unique properties of 2D materials.