{'fraction_non_alphanumeric': 0.029093027772382988, 'fraction_numerical': 0.011730433093804622, 'mean_word_length': 4.429354702657107, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective medium theory for disordered two-dimensional graphene

We develop an Effective Medium Theory to study the electrical transport properties of disordered graphene. The theory includes non-linear screening and exchange-correlation effects allowing us to consider experimentally relevant strengths of the Coulomb interaction. Assuming random Coulomb impurities, we calculate the electrical conductivity as a function of gate voltage describing quantitatively the full cross-over from the fluctuations dominated regime around the Dirac point to the large doping regime at high gate voltages. We find that the conductivity at the Dirac point is strongly affected by exchange correlation effects.


I. INTRODUCTION
Since the experimental realization of graphene 1,2,3 the problem of understanding the electronic conductivity of graphene has generated a large experimental and theoretical effort, given its fundamental importance and its technological relevance. The unusual transport properties of graphene 1 arise mostly from its "Dirac spectrum" -a linear, zero-gap dispersion relation of chiral Fermion carriers (electrons and holes), with the charge neutral "Dirac point" defined as the singular point (with vanishing density of states) where the electron and hole bands touch. In this work we present a theory that is able to describe at a qualitative and semi-quantitative level most of the current bulk transport experimental results on exfoliated graphene. The theory is designed specifically to handle the disorder induced density inhomogeneities present in graphene and dominating the physics at the Dirac point where the average carrier density is zero. The theory is general: it can be applied to systems other than graphene (as formulated, it applies readily to graphene bilayers); it does not depend directly on the source of the spatial inhomogeneities; and it can be used to calculate other graphene transport properties such as thermo-transport 4 and magnetotransport. The theory takes into account non-linear screening effects and exchange-correlation terms via a Hartree-Fock Local Density Approximation (LDA); in this sense the theory is non-perturbative and therefore, within the limits of this approximation, is valid for large values of the graphene's fine structure constant r s ≡ e 2 / v F ǫ, where v F is the Fermi velocity and ǫ the dielectric constant. For graphene charge transport, the theory is powerful being able to: 1) Predict values of conductivity at the Dirac point, σ min , in very good agreement with experiments; 2) Describe quantitatively the dependence of σ on the doping ∝ V g including the crossover or plateau region, i.e. the region where V g is finite but smaller than the range of values for which σ(V g ) exhibits a linear behavior; 3) Show that many-body terms are essential to understand the transport properties of graphene close to the Dirac point; 4) Explain the observed dependence of σ min on the sample quality; 5) Explain the dependence of σ min on r s observed experimentally 5 and predict its behavior over a wide range of r s , a prediction that could be tested experimentally; 6) Show that at large doping the spatial density fluctuations do not modify the linear dependence of σ on V g ; and 7) Predict the transport properties at very low but finite long-range disorder; a regime not accessible in previous theories.
The success of the transport theory presented in this work relies on the proper treatment of the effects of the spatial carrier inhomogeneity. In a metal or semiconductor, defects play a double role 6 : they act as scattering centers and locally modify the conduction-band carrier density. A defect effectively shifts the local Fermi level away from its average value, which has dramatic consequences on the properties of graphene close to the Dirac point where the average density, n , is zero and therefore the spatial density fluctuations completely dominate the physics. The source of the disorder 7 can be of different nature: short-range scatterers, such as atomic defects in the graphene lattice structure, ripples 8 or long range scatterers such as charged impurities, but all will induce spatial inhomogeneities of the carrier density that close to the Dirac point result in the formation of electronhole, e-h, puddles 9,10,11,12 . This situation has emerged as the one realized in current experiments on exfoliated graphene. The ideal case of zero average density and zero amplitude of the density fluctuations has been considered in previous works that predicted σ min to be either 0 or ∞, depending on the nature of the disorder, or, in the limit of vanishing disorder, finite and universal 13,14,15,16,17,18,19 .
The formalism that we present does not depend on the source of the inhomogeneities, however recent experimental results 5,20,21 provide convincing evidence that random charged impurities -located in the graphene environment -are the dominant source of disorder in graphene. Therefore for concreteness we assume the carriers in graphene to be subject to a disorder potential due to a 2D random distribution of uncorrelated Coulomb impurities with surface density n imp , placed at a distance d from the graphene layer. The paper is organized as follows: in section II we develop and justify the effec-tive medium theory for graphene in presence of a 2D random distribution of charged impurities, in section III we present the results for the conductivity of disordered graphene and finally in section IV we summarize the main results and discuss the range of validity of the approach presented.

II. THEORETICAL MODEL
As said in the introduction, in the remainder of the paper we only consider the disorder induced by random charged impurities in the graphene environment, motivated by the strong experimental evidence 20,21 that this type of disorder dominates in current experiments on exfoliated graphene on SiO 2 substrate. In order to be able to correctly describe the properties of graphene close to the Dirac point an accurate characterization of the strong, disorder induced, density fluctuations is necessary. This requires the calculation of the ground state for interacting massless Dirac electron in presence of disorder, a very challenging problem. In absence of disorder one can use a Kohn-Sham-Dirac Density-Functional-Theory, DFT, that appropriately takes into account the fact that in graphene the low energy quasiparticles behave as chiral massless Dirac fermions 22 . In this approach many-body effects enter via a local-density approximation exchange-correlation potential. In presence of disorder the DFT-LDA approach is computationally very expensive and as a consequence can only be used to calculate the ground state for few disorder realizations 22 and does not allow for a statistical characterization of the carrier density fluctuations that requires the calculation of the disorder averaged quantities. A computationally cheaper approach is the microscopic Thomas-Fermi-Dirac, TFD, theory that recently some of us have developed and used to study the electronic structure of disordered graphene 12 . The TFD approach is very similar to the DFT-LDA one, the difference being that in the TFD theory also the kinetic energy is approximated by a local-density functional, whereas in the DFT-LDA the kinetic term is a differential operator that operates on the Kohn-Sham orbitals. For single disorder realizations the TFD results can be directly compared to the DFT-LDA results and we found them to be in very good agreement for impurity density values, n imp , relevant for current experiments on exfoliated graphene: n imp > 10 9 cm −2 . The Thomas-Fermi theory in general is accurate when the ground-state density varies on length scales bigger than the the Fermi wavelength, λ F , i.e.: Because k F = √ πn inequality (1) can be rewritten as: Inequality (2) is equivalent to say that in any region with homogeneous carrier density the number of carriers has to be bigger than 1. As shown below the non-linear screening properties of graphene guarantee that in presence of charged impurities this condition is satisfied for graphene also close to the Dirac point.
To summarize the main properties of the graphene carrier distribution close to the Dirac point in the lower color plot of Fig. 1 we show a typical result obtained using the TFD theory for a given disorder realization with n imp = 5 × 10 11 cm −2 , d = 1 nm and r s = 0.8, values that are typical for current exfoliated graphene samples and accurately estimated by comparing the Boltzmann theory and the experimental results at high doping 20 . This plot exemplifies the main qualitative properties of the carrier distribution close to the Dirac point: I) the distribution is characterized by e-h puddles in quantitative agreement with surface probe experiments 11,23 ; II) the typical size of an e-h puddle, defined as a region with same-sign charges, is of the order of the sample size, L, as expected for a semimetal close to the neutrality point; III) n(r) is characterized by two types of inhomogeneities, namely, wide regions (i.e. big puddles spanning the system size) of low density (an example is shown in Fig 1 with the blue contour that contains ≈ 10 electrons); and few narrow regions of high density (an example is shown in Fig 1 with the white contour containing 2 electrons); IV) the narrow puddle regions (peaks, dips) have a correlation length ξ that depends on n imp , r s and d. 12 For the parameter values used in Fig. 1, typical for current experiments, ξ ≈ 10nm. The combination of the high density in the peaks/dips and the fact that in the low density regions n(r) varies over scales much bigger than ξ guarantees that the inequality (2) is satisfied over the majority of the graphene sample and therefore justifies the TFD theory. The properties I)-IV) are essential for the development and justification of the effective medium theory, EMT, for graphene presented in this work. In presence of disorder, the problem of graphene transport then becomes one of how to properly take into account the strong density inhomogeneities shown in Fig. 1. At high doping the graphene carrier density can be assumed to be homogeneous. In this regime, σ has been shown experimentally to scale linearly with n . The theory that has been able to describe most successfully, away from the Dirac point, the large number of experimental results for graphene transport is scattering from random charge impurities 9,24,25,26 , where, within the RPA-Boltzmann approximation, we have: where F (r s , d) was given in Ref. 10. An important conceptual step to derive a correct transport theory valid close to the Dirac point is to introduce a local spatially varying "puddle" conductivity σ(r). The local conductivity is a well defined quantity as long as σ(r) varies on length scales that are bigger than the mean free path l i.e.: Because, as shown in Fig. 1, close to the Dirac point the density profile is characterized by large puddles with n = 0 the RPA-Boltzmann theory is expected to be valid locally inside the e-h puddles. To estimate σ in a puddle we therefore use Eq. (3) replacing n by its local value n(r) so that in Eq. (4) we can replace ∇σ(r)/σ(r) with ∇n(r)/n(r). Considering that l = hσ/(2e 2 k F ), using Eq. (3), for disordered graphene inequality (4) takes the form: For graphene on SiO 2 (r s = 0.8) 1/F (r s , d) ≈ 10. Close to the Dirac point there are few narrow regions of high density and size ξ(d, n imp , r s ). For r s = 0.8, d = 1 nm, n imp ∼ 10 11 cm −2 (parameter values relevant for current experiments on exfoliated graphene) is ξ ∼ 10 nm l as also shown in the lower panel of Fig. 1. However, the regions where the carrier density changes over length scales of the order of ξ are very sparse and the density landscape is characterized by wide regions where the density is approximately constant on much bigger length scales, of the order of the system size. Because the small puddles of size ≈ ξ are isolated (i.e. do not form a path spanning the whole sample) and occupy a small area fraction (we estimated it to be less than 20% 12 ) we can neglect their contribution to the graphene conductivity. In addition, given their high density, steep gradients at the boundaries, and the fact that ξ ≪ l, the small puddles should provide very small resistance to the carriers' motion. In the large regions where n(r) is almost uniform n is approximately equal to the root mean square of the density distribution n rms ∼ n imp 10,12 , and as a consequence in these regions the inequality (4) is satisfied for experiments on bulk exfoliated graphene for which we have L ∼ 1µm and n imp is in the range [10 11 − 10 12 ] cm −2 . We conclude that the non-linear screening properties of graphene that justify the use of the TFD theory, and the linear dependence of σ on n (Eq. (3)), allow for the definition of a local conductivity for graphene at the Dirac point. The conductivity spatial profile corresponding to the density profile of the lower color plot of Fig. 1, is shown in the upper color plot of the same figure.
Close to the Dirac point graphene can then be thought of as a highly heterogeneous material for which the EMT should be appropriate to calculate its transport properties. Even though the EMT was proposed long ago for heterogeneous materials 6,27 , the development and justification for graphene is not obvious and is the central result of this work. To derive the EMT one considers a single homogeneous region (e.g. area encircled by the white perimeter in the color plot for σ in Fig. 1) with conductivity σ embedded in a uniform medium with effective conductivity, σ EMT and then requires the spatially integrated inhomogeneities of the electric field E, due to the inclusion of the region with σ = σ EMT , to be equal to zero. For a 2D system this requirement translates into the equation: The EMT is valid if the conductance across the regions where σ is homogeneous is higher than the conductance inside these regions. For graphene at the Dirac point this requires that the conductance across the disorder induced p-n junctions, PNJ, (shown as black lines in the lower color-plot of Fig. 1) is higher than the conductance inside the puddles. Close to a PNJ Eq. (3) is not valid.
If D is the length scale over which n changes across a PNJ and k F = √ πn is the Fermi wave vector at the sides of the PNJ, then using TFD results for n(r) we find k F D ∼ 1; for single disorder realizations this result can also be inferred from the DFT-LDA calculations 22 . In this limit, the conductance per unit length of a PNJ,g PNJ is given by 28,29 : Let p be the perimeter of a typical puddle and Γ its form factor then the condition that the conductance across the PNJ must be bigger than the conductance of a puddle takes the form: . i.e.
where n is the density at the sides of the PNJ: n ∼ n rms ∼ n imp 10,12 . Given the macroscopic size of the majority of the e-h puddles the length of the PNJ, p, is of the order of the sample size L ∼ 1µm, and then the inequality (9) is satisfied. We can therefore conclude that the EMT developed here properly describes the electronic transport properties of bulk disordered graphene samples and that close to the Dirac point, transport in bulk graphene is dominated by scattering processes inside the puddles and not across their boundaries as was assumed in previous works 30 . Disorder averaging Eq. (6) we find: where P (σ) is the probability for the local value of σ. Using Eq. (3) we have: To obtain an accurate value of σ EMT it is crucial to select the correct form of P [n]. Choosing P [n] to be a Gaussian distribution of width n rms , for σ min we find σ G EMT ≈ (e 2 /h)[1.1460/F (r s , d)][n rms /n imp ]. Using a Lorentzian approximation for P [n] with width n L , we find σ L EMT = (e 2 /h)[2/F (r s , d)](n L )/(n imp ). These analytical results allow for a direct comparison with other results for σ min in the literature 10 . In our approach P [n] is accurately calculated using the TFD theory 12 that takes into account nonlinear screening, exchange and correlation contributions 22,31 , and can handle the crossover plateau region where P [n] has a bimodal character that is not captured by the Gaussian and Lorentzian approximations. Only close to the Dirac point and for small values of r s do the Gaussian and Lorentzian approximations give results that are comparable to the ones obtained using P [n] calculated with the TFD theory, P [n] T F D . The TFD approach allows us to go to higher values of r s and lower n imp thereby also providing a range of validity for analytical theories.

III. RESULTS
In this section we present the results for the conductivity of disordered graphene obtained using the effective medium theory developed in the previous section and the probability distribution, P [n] T F D , obtained using the TFD theory. Shown in Fig. 2 is σ(V g ). We see that our theory correctly predicts a finite value of σ at the Dirac point that is of the same order as that measured experimentally. At high gate voltages the linear scaling of σ as a function of V g is recovered. One important achievement of our theory is the ability to correctly describe the crossover of σ from its minimum, σ min , at V g = 0, to its linear behavior at high gate voltages. Moreover our work shows the importance at low gate voltages of the exchange term. In Fig. 3 we present the result for σ min as function of n imp . We see that our approach explains the non-universality of σ min from sample to sample and that exchange reduces the dependence of σ min on n imp for impurity densities in the range [10 10 −10 12 ] cm −2 . For n imp 10 10 cm −2 σ min grows quite rapidly as n imp is decreased. The same qualitative dependence of σ min on n imp has been recently observed experimentally as shown in Fig. 5 (a) of Ref. 21. In Fig. 4, the dependence of σ min on d is shown for a given value of n imp . We see that σ min depends weakly on d and this dependence is even weaker if exchange-correlation effects are taken into account. In Fig. 5 we show the results for σ min as a function of r s . The solid (dashed) line shows the values of σ min obtained including (neglecting) exchange. We see that σ min has a non-monotonous behavior due to the fact that r s strongly affects both the carrier density spatial distribution (by controlling the strength of the disorder potential, screening and exchange); and the scattering time τ . From the RPA-Boltzmann theory in presence of charged impurities we have that τ = √ n/(πF (r s , d)v F n imp ), so that τ depends on r s only through F (r s ). In Fig. 6 we show 1/F (r s ) for d = 1, to show how τ decreases as r s increases. On the other hand as r s increases the width of P  the density fluctuations is plotted as a function of r s . We see that in presence of exchange n rms grows much more slowly than when exchange is neglected. This is due to the fact that in graphene, contrary to parabolic band materials, close to the Dirac point, exchange penalizes density fluctuations. As a consequence for r s 0.1 the behavior of σ min is qualitatively different when exchange is taken into account. In Ref. 5, r s was varied from 0.56 to 0.8 and σ min was found to remain constant in agreement with the results shown in Fig. 5 in presence of exchange. We can then draw the important conclusion that to understand, even at a qualitative level, the transport properties of current graphene samples exchange contributions must be taken into account.

IV. DISCUSSION AND CONCLUSIONS
In this paper we have derived and justified the effective medium theory for disordered bulk graphene. In particular we have considered the case when the disorder is due to random charged impurities in the graphene environment. Charge impurities have been shown to be the most important source of disorder for current experiments on exfoliated graphene. There are however other sources of disorder like ripples 32,33,34,35 and atomic defects 36 that might play an important role in determining the transport properties close to the Dirac point of graphene systems other than bulk exfoliated graphene. For example there is evidence that atomic defects might be the dominant source of disorder in epitaxial graphene. We should emphasize that our motivation was to calculate the intrinsic conductivity of graphene close to the Dirac point. For this reason we have neglected the effects of contacts that have been shown to be important 37,38,39,40 but that do not modify the value of the graphene intrinsic con- ductivity and in experiments can be minimized by performing four-terminal measurements. Our theory does not apply to graphene systems with a band gap 41 and is only valid in the diffusive limit. To model the disorder we have assumed the impurities to be randomly distributed in a 2D plane placed at a distance d from the graphene. The only two input parameters that enter the theory are the charge impurity density, n imp , and d. These two parameters have been accurately extracted by comparing the RPA-Boltzmann theory and the experimental results at high doping 20,21 . We have used the EMT developed here to study in details the transport properties of graphene close to the Dirac point, where the physics is dominated by the strong density fluctuations induced by the disorder. The parameters n imp and d that enter the EMT should not be regarded as fitting parameters: to test our results for a given sample, n imp and d should be extracted by comparing the value of the mobility at high doping with the one provided by the RPA-Boltzmann theory, and then used to obtain the EMT results at the Dirac point. A key ingredient for the EMT is the accurate calculation of the density probability distribution, P (n), close to the Dirac point taking into account non-linear screening and exchange-correlation effects. We have used P (n) provided by the TFD theory 12 . The TFD theory is justified for most of the current experiment on exfoliated graphene and gives results that, for single disorder realizations, are in agreement with the DFT-LDA approach developed in Ref. 22. On the other hand due to its computational cost the DFT-LDA does not allow the calculation of disorder averaged quantities and therefore a statistical characterization of the density profile. The TFD-EMT approach that we have developed is valid when the inequalities (1), (5), (9) are satisfied. In these inequalities n is the density inside the electronhole puddles. In the small regions, of size ξ ≈ 10 nm, because of the high carrier density, inequality (1) is always satisfied. As we have shown, most of the graphene is covered by large puddles of size L p of the order of the system size, L, and with density of the order of n rms . As shown in Ref. 10,12 n rms ≈ ηn imp , where η is a number of order 1. Using these results we find that (5) and (9) are equivalent apart from a factor of order 1 and that the inequalities (1), (9), take the form: respectively. Given that for the cases of interest 1/F (r s , d) ≈ 10, see Fig. 6, the second inequality in (12) is the most stringent: the condition that the conductance across the PNJ must be bigger than the conductance inside the puddles, or equivalently the condition that the mean free path must be smaller than the puddles, sets the minimal size that the puddles must have for the EMT to be valid. Considering that L p ∼ L > 1 µm and n imp ∈ [10 10 − 10 12 ] cm −2 , we see that for bulk samples of exfoliated graphene the inequalities (12) are satisfied. We conclude that for sufficiently disordered and large graphene samples the TFD-EMT approach is accurate to estimate the graphene transport properties, this conclusion is consistent with the results of Ref. 42 in which the graphene conductivity in presence of Gaussian disorder obtained using a full quantum mechanical calculation was found to be in agreement with the semiclassical Boltzmann one even at zero doping provided the disorder is strong enough.
Using the TFD-EMT approach we find that at the Dirac point the theory gives a finite minimum conductivity in very good agreement with that observed experimentally 5,20,21 . In addition the theory is able to describe the crossover of the conductivity from its minimum value at the Dirac point to its linear behavior at high doping. We have also calculated the dependence of σ min on n imp and found it to be in good agreement with recent experimental results 21 . We obtain the dependence of σ min on r s , finding that for r s 0.1, σ min decreases very rapidly with r s but for r s 0.2, exchangecorrelation terms are qualitatively important and must be included to correctly describe the minimum conductivity. This result is in agreement with the experimental results presented in Ref. 5 in which the minimum of conductivity is found to be the same for r s = 0.5 and r s = 0.8. By comparing our results for σ min as a function of n imp and r s to the experimental ones we find that exchangecorrelation effects must be included when evaluating the transport properties of graphene at the Dirac point. The reason is that exchange-correlation effects, close to the Dirac point, at low n imp and high r s , strongly suppress the amplitude of the density fluctuations and therefore strongly affect the transport properties. The success of the EMT, coupled with the microscopic disorder-induced graphene electronic structure calculation, in obtaining the graphene conductivity near the Dirac point indicates that this technique should be useful in calculating many other properties of graphene in the theoretically difficult inhomogeneity-dominated regime near the charge neu-