Effects of the structure of charged impurities and dielectric environment on conductivity of graphene

# Effects of the structure of charged impurities and dielectric environment on conductivity of graphene

R. Aničić Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1    Z. L. Mišković Department of Applied Mathematics, and Waterloo Institute for Nanotechnology, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1
July 13, 2019
###### Abstract

We investigate the conductivity of doped single-layer graphene in the semiclassical Boltzmann limit, as well as the conductivity minimum in neutral graphene within the self-consistent transport theory, pointing up the effects due to both the structure of charged impurities near graphene and the structure of the surrounding dielectrics. Using the hard-disk model for a two-dimensional (2D) distribution of impurities allows us to investigate structures with large packing fractions, which are shown to give rise to both strong increase in the slope of conductivity at low charge carrier densities in graphene and a strongly sub-linear behavior of the conductivity at high charge carrier densities when the correlation distance between the impurities is large. On the other hand, we find that a super-linear dependence of the conductivity on charge carrier density in heavily doped graphene may arise from increasing the distance of impurities from graphene or allowing their clustering into disk-like islands, whereas the existence of an electric dipole polarizability of impurities may give rise to an electron-hole asymmetry in the conductivity. Using the electrostatic Green’s function for a three-layer structure of dielectrics, we show that finite thickness of a dielectric layer in the top gating configuration, as well as the existence of non-zero air gap(s) between graphene and the nearby dielectric(s) exert strong influences on the conductivity and its minimum. While a decrease in the dielectric thickness is shown to increase the conductivity in doped graphene and even gives rise to finite conductivity in neutral graphene for a 2D distribution of impurities, we find that an increase in the dielectric thickness gives rise to a super-linear behavior of the conductivity when impurities are homogeneously distributed throughout the dielectric. Moreover, the dependence of graphene’s mobility on its charge carrier density is surprisingly strongly affected, quantitatively and qualitatively, by the graphene-dielectric gap(s) when combined with the precise position of a 2D distribution of charged impurities. Finally, we show that the conductivity minimum in neutral graphene is increased by increasing the correlation distance between the impurities, reduced by increasing the graphene-dielectric gap, and increased by decreasing the dielectric thickness in a top-gated configuration, even though the corresponding residual charge carrier density is reduced by decreasing the dielectric thickness.

graphene, conductivity, charged impurities, dielectric screening
###### pacs:
73.22.Pr, 72.80.Vp, 81.05.ue

## I Introduction

Graphene is a realization of a two-dimensional (2D) material made of carbon atoms strongly bonded in a honeycomb–like lattice, exhibiting a Dirac-like spectrum for low-energy excitations of its electrons, which has been under intense scrutiny for possible applications in electronics, photonics,Avouris_2012 () and biochemical sensing.Allen_2010 () Being an all-surface material renders graphene extremely sensitive to the incident electromagnetic fields and to the dielectric properties of the surrounding matter,Newaz_2012 () which is both a blessing and a curse from a technological point of view. While the use of external gates and/or controlled adsorption of atomic and molecular species present an efficient means for inducing precise concentrations of charge carriers in graphene,Chen_2008 () the presence of indeterminate amounts of charged impurities, which may be trapped in a substrate or directly adsorbed on graphene, render quantitative details of many measurements of graphene’s electronic and optical properties ”sample dependent”.Tan_2007 () In addition, integrating graphene in layered structures with different material properties may bring additional issues due to uncertainties in the geometric structure and the chemical composition of such structures.Fallahazad_2012 (); Hollander_2011 ()

Possibly the most intriguing manifestation of the presence of charged impurities is the famed minimum in the DC conductivity of single-layer graphene in the limit of vanishing doping, i.e., when the average density of induced charge carriers in graphene approaches zero.Tan_2007 (); Sarma_2011 () It was shown that the minimum conductivity may be explained by the manifestation of a system of electron-hole puddles in graphene due to corrugation of the electrostatic potential that arises from a spatial distribution of the charged impurities in a substrate. PNAS_2007 () On the other hand, the conductivity in heavily doped graphene layers often exhibits sub-linear behavior, or saturation with increasing charge carrier density, which is often explained by the presence of short-range scatterers in graphene, presumably arising from atomic-size defects in the carbon lattice. Yan_2011 () However, it turned out that spatial correlation among the nearby charged impurities may provide an alternative and more plausible explanation of the conductivity saturation in single-layer graphene. Li_2011 (); Yan_2011 () Moreover, the atoms adsorbed on graphene often show tendency of clustering and forming islands, which may additionally affect the mobility of charge carriers in graphene. McCreary_2010 ()

As far as the structure and composition of the surrounding material is concerned, preference is usually given to insulators and metals that only engage in weak interactions with graphene of the van der Waals type, leaving the structure of its electron bands largely intact in the vicinity of the Dirac point.Wehling_2009 () Those interactions are characterized with relatively large spatial gaps between graphene and the nearby material, on the order of several Ångströms, which reduce the dielectric screening by that material and often exhibit significant fluctuations in their size due to the surface roughness of the material.Ishigami_2007 () Furthermore, when graphene is top gated with a layer of high- dielectric material, the mobility of its charge carriers may be affected by a strong image interaction with the metallic top gate.Fallahazad_2012 (); Hollander_2011 (); Ong_2012 () Finally, for electrolytically top-gated graphene, the presence of mobile ions in the nearby electrolyte may provide additional screening of the charged impurities in a solid substrate.Chen_2009 (); Miskovic_2012 ()

All of the above examples of the effects of charged impurities near graphene and the structure of the surrounding dielectrics play important roles in its charge carrier transport, plasmon dispersion in doped graphene, and graphene’s capacitance, which are of interest in electronics, photonics, and sensing, respectively. It was recently shown that those effects may be conveniently modeled by using Green’s function (GF) for the Poisson equation for a layered structure, Ong_2012 (); Miskovic_2012 () which is easily combined in a self-consistent manner with the polarization function of graphene within the random phase approximation (RPA) when graphene is modeled as a zero-thickness material. Castro_2009 () In this work, we illustrate such approach to modeling the conductivity of single-layer graphene with large area by considering a three-layer structure of the surrounding dielectrics and using an expression for the conductivity that results from the semiclassical Boltzmann transport (SBT) theory for doped graphene. Sarma_2011 () However, that expression is derived here via the Energy loss method (ELM),Gerlach_1986 () which explicitly evaluates the friction force on a system of external charges with the spatial distribution that moves rigidly parallel to graphene.Allison_2009 (); Allison_2010 (); Radovic_2012 () Hence, the ELM has an added utility as it may be used in studying other processes, such as sliding friction of molecular layers physisorbed on graphene,Krim_2012 () or probing the streaming potential in a flowing electrolyte by a graphene based sensor,Newaz_2012_b () which will be tackled in future work.

In this work we focus on several effects in the DC conductivity of graphene. First, we explore the effects of long correlation distances among impurities that give rise to large packing fractions, which cannot be described by a simple step-like pair correlation function.Yan_2011 (); Li_2011 () For that purpose we use an analytically parameterized model of hard disks (HD) due to Rosenfeld,Rosenfeld_1990 () which gives reliable results for packing fractions up to the freezing point of a 2D fluid. Next, whereas all the previous studies assumed that charged impurities reside in a plane parallel to graphene, our statistical formulation of the theory allows for a fully three-dimensional (3D) spatial distribution of impurities that may reside at a range of distances from graphene. In addition, we allow that individual impurities may be characterized by atomic-like form factors, which include a finite dipole moment and a spatial spread that accounts for the existence of disk-like clusters near graphene. Furthermore, by taking advantage of the electrostatic GF for a three-layer structure, we also study the effects that arise in conductivity of graphene due to finite thickness of a nearby dielectric and a finite gap between graphene and the nearby dielectrics. Finally, the above effects are also studied in the context of the conductivity minimum within the Self-consistent transport (SCT) theory.PNAS_2007 ()

Specifically, in this paper we show via the HD model that large correlation distances between charged impurities may give rise to significantly larger initial slopes of the conductivity (or larger mobility) at lower charge carrier densities, as well as to a more pronounced saturation, or sub-linear behavior of conductivity at higher densities than in the case of small correlation distances. Next, the effects of clustering of charge impurities, as well as the increasing distance from graphene are confirmed to give rise to super-linear dependence of conductivity on charge carrier density in heavily doped graphene, in agreement with observationsMcCreary_2010 () and modeling,PNAS_2007 () respectively. Impurities with finite dipolar polarizability are shown to give rise to electron-hole asymmetry in the conductivity as the sign of charge carrier density changes, which may be related to experimental observations in some graphene samples.Tan_2007 () Regarding the geometrical factors of a nearby dielectric layer, we find an increase in both the conductivity and mobility of graphene when the layer thickness decreases in the case of a 2D distribution of impurities, whereas a homogeneous 3D distribution of impurities gives rise to a super-linear behavior of the conductivity with increasing layer thickness. Most intriguingly, we find a strong effect on the mobility of graphene due to the presence of a finite gap between graphene and the nearby dielectrics in conjunction with the varying position of impurities, which was not previously considered in the modeling of the transport properties of graphene, but was observed in studying the polarization forces on external charges. Allison_2009 () Finally, a similarly strong effect of the finite gap between graphene and a nearby dielectric is also demonstrated in the minimum conductivity within the SCT theory.PNAS_2007 ()

After outlining the theoretical model in the next section, we discuss our numerical results, and give concluding remarks. In the Appendices we outline a derivation of the electrostatic GF and provide details for several models of the charged impurity structure. Note that, unless otherwise explicitly stated, we use gaussian electrostatic units where , with being the dielectric permittivity of vacuum.

## Ii Theory

We assume that a single-layer graphene sheet of large area is embedded into a stratified structure so that it lies parallel to layers of various dielectrics with abrupt interfaces among them, as shown in Fig. 1. Using a 3D Cartesian coordinate system with coordinates , the entire structure may be then considered translationally invariant (and is assumed to be isotropic) in the directions of a 2D position vector . Furthermore, assume that a system of charged particles is distributed throughout the structure and is moving rigidly at a constant velocity parallel to graphene. If the stationary volume density of charges in the moving frame of reference is given by , then the corresponding volume density in the rest frame of graphene (the laboratory frame of reference) is given by .

This notion of a rigidly moving distribution of external charges may be related to several realistic physical situations where the relative motion of particles with respect to each other may be treated as adiabatic at the time scale of the charge carrier dynamics in graphene. Examples include sliding of a film of adsorbed molecular layers across graphene,Krim_2012 () flow of a molecular fluid that contains dissolved ions in thermal equilibrium,Newaz_2012_b () or propagation of ionized fragments that result from planar Coulomb explosion of a cluster grazingly scattered from graphene. Song_2005 () In each of those examples, the movement of external charged particles gives rise to energy dissipation due to excitations of charge carriers in graphene.

Conversely, one my reverse the frames of reference and consider the regime of steady-state electric conduction in graphene where its charge carriers move with a constant drift velocity . In this case the distribution of external particles is static in the laboratory frame and hence may be used to model fixed charged impurities near graphene. If the speed is sufficiently low, then the electrical resistivity of graphene may be related to energy dissipation due to scattering of its charge carriers on external charged impurities, giving rise to Ohmic heating of graphene. This idea of reversing the frames of reference is a basis of the ELM that was developed for studying the transport properties of semiconductor heterostructures by means of the dielectric response formalism for their conducting electrons.Gerlach_1986 () This method was used successfully in studying the scattering of conduction electrons on interface roughnessKaser_1995 () and polarizable scattering centers,Kaser_1997 () as well as in discussing vibrational damping in adsorbed layers due to surface resistivity,Persson_1991 () and in studying optical properties of thin films for solar energy materials.Jin_1988 () Moreover, this same idea of the equivalence of a drag force on a uniformly moving system of impurities and the total force on the electron fluid in doped graphene was recently applied to evaluate the conductivity of graphene within the semiclassical hydrodynamic model for its charge carriers.Mendoza_2013 ()

We note that the ELM gives an expression for the conductivity of doped graphene, which is identical to that obtained by the SBT theory,Sarma_2011 () but we chose ELM because it yields the drag force on externally moving charges as a side result that may be more directly used in modeling other processes, such as sliding friction of molecular layers physisorbed on graphene Krim_2012 () or probing the streaming potential in a flowing electrolyte by a graphene based sensor,Newaz_2012_b () to mention a few.

### ii.1 Energy loss method

To be specific, we assume that the system of external charges consists of particles, each carrying a total charge of (where is the proton charge) that is distributed around the center of the particle according to some function , such that with . If the th particle is centered at the position in the moving frame of reference, we may write for the total density of charges in that frame

 ρ0(r,z)=eN∑j=1Δj(r−rj,z−zj). (1)

Given that the positions of external particles, as well as their individual charge densities are statistically distributed, we denote their joint ensemble average by . Assuming that this distribution is translationally invariant in the directions of , we note that can only be a function of the perpendicular coordinate . Therefore, assuming that the equilibrium areal number density of charge carriers is uniform across graphene, its value will be determined by both the function and the potential applied through the external gates. We assume that has a sufficiently large value allowing us to neglect the effects of fluctuations in the charge carrier density in graphene on its screening properties. On the other hand, we assume to be small enough to allow the use of a 2D response function for graphene’s electrons in the approximation of Dirac fermions. Wunsch_2006 (); Hwang_2007 () Those requirements practically limit our considerations of graphene’s DC conductivity within the SBT theory to an approximate range of doping densities 10 cm 10 cm (we assume unless stated otherwise).

We further define the fluctuation in the charge density of external particles by and use it in the Poisson equation, allowing us to express the resulting fluctuation of the electrostatic potential, , in terms of the electrostatic GF for the entire system, , as

 δΦ(R,t)=∫d3R′∞∫−∞dt′G(R,R′;t−t′)δρ(R′,t′). (2)

Using a tilde to denote the Fourier transform (FT) of various quantities with respect to the 2D position () and time (), the above expression is recast in the form

 δ˜Φ(q,z,ω)=∞∫−∞dz′˜G(q;z,z′;ω)δ˜ρ(q,z′,ω), (3)

where

 δ˜ρ(q,z,ω) = ∫d2r∞∫−∞dte−iq⋅r+iωtδρ0(r−vt,z) (4) = 2πδ(ω−q⋅v)δ˜ρ0(q,z)

defines the relation between the FTs of the fluctuations of the external charge densities in the two reference frames. Here, is defined via the FT of the external charge density in the moving frame of reference,

 ˜ρ0(q,z)=eN∑j=1˜Δj(q,z−zj)e−iq⋅rj. (5)

It may be shown that the ensemble average of the energy loss rate is given by Mowbray_2010 ()

 ⟨dWdt⟩ = −∫d3R⟨δρ(R,t)∂∂tδΦ(r,z,t)⟩ (6) = i∫d2q(2π)2(q⋅v)∫dz∫dz′˜G(q;z,z′;q⋅v) ×⟨δ˜ρ0(−q,z)δ˜ρ0(q,z′)⟩.

On using the symmetry properties of the FT of the full GF, and , where denotes complex conjugation, one notices that only the imaginary part of the factor in Eq. (6) contributes to the energy loss. Furthermore, assuming that graphene has a zero thickness and is placed in the plane , we may express in terms of the (real valued) 2D FT of the GF for the dielectric environment without graphene, as given in Eq. (35). Thus, Eq. (6) may be rewritten as

 ⟨dWdt⟩ = ∫d2q(2π)2VC(q)(q⋅v)I[−1ϵ(q,q⋅v)] (7) ×⟨δ˜N(−q)δ˜N(q)⟩,

where we have defined a dielectric function that describes the dynamic screening of external electrostatic fields in the plane due to the polarization of the entire system as

 ϵ(q,ω)=ϵbg(q)+VC(q)χ(q,ω), (8)

with being an effective background dielectric function due to the polarization of the system without graphene, the in-plane FT of the Coulomb potential, and a 2D polarization function of noninteracting electrons in graphene. Wunsch_2006 (); Hwang_2007 () Moreover, in Eq. (7) we have introduced the fluctuation in an effective areal (or surface-projected) number density of external particles, , which is defined via its 2D FT as , with

 ˜N(q)≡1e∫dzψ(q,z)˜ρ0(q,z)=N∑j=1Fj(q)e−iq⋅rj, (9)

where

 ψ(q,z)=˜G(0)(q;zg,z)˜G(0)(q;zg,zg) (10)

is a profile function that takes into account the decay of the Coulomb interaction throughout the system with increasing distance from graphene, and

 Fj(q)=∫dzψ(q,z)˜Δj(q,z−zj) (11)

may be considered to be a weighted form factor of the th particle.

### ii.2 Friction regime and conductivity of graphene

In order to use the ELM to obtain the DC conductivity of graphene, we require an ensemble average of the energy loss rate to the lowest order in speed , which corresponds to the friction regime for slowly moving external charges. This is easily accomplished by expanding the loss function in Eq. (7) to the leading order in frequency by using the truncated expansion for the polarization function of doped graphene,Allison_2010 ()

 χ(q,ω)=χs(q)+iωπℏv2FU(2kF−q)√(2kFq)2−1, (12)

where is the static polarization function and is an average value of the Fermi wavenumber for Dirac electrons in graphene. We further define the auto-correlation function of charged impurities in Eq. (7) by

 S(q)≡1N⟨δ˜N(−q)δ˜N(q)⟩, (13)

and note that it only depends on the magnitude when the distribution of impurities is isotropic in the directions parallel to graphene. This allows us to finally obtain from Eq. (7)

 ⟨dWdt⟩=2ℏkFNv2r2s2kF∫0dq√1−(q/2kF)2[ϵbg(q)+4kFrs/q]2S(q), (14)

where with being the Fermi speed of Dirac electrons. Note that the quantity is an average total stopping, or drag force that acts on the moving system of external charges, Radovic_2012 () which may be used to, e.g., evaluate the friction coefficient for an adsorbed layer on graphene from the expression in the limit .Allison_2010 (); Krim_2012 ()

Within the ELM, by reversing the frames of reference one may express the energy loss rate in graphene by the standard expression of classical electrodynamics,

 (15)

where is the current density of charge carriers in graphene, induced by a constant electric field applied across graphene, and is its DC conductivity. Assuming a uniform charge carrier density across graphene, we may write in a steady-state regime, which gives , where is the macroscopic area of graphene. We discard possible contribution to the conductivity of graphene coming from charge carrier scattering on short-ranged impurities, and we limit our considerations to sufficiently low temperatures to be able to neglect the contribution from scattering on phonons. Thus, the final expression for the DC conductivity takes a form that is familiar from the SBT for doped graphene,Sarma_2011 (); Li_2011 ()

 σ=e2h¯nnimp21∫0duu2√1−u2[2+ursϵbg(2kFu)]2S(2kFu), (16)

where is the mean areal number density of external charged particles.

### ii.3 Variance of the potential in graphene and minimum conductivity

Equation (16) implies that the conductivity obtained within the SBT theory as a function of the average equilibrium charge carrier density in graphene, , should vanish in a linear manner close to the neutrality point, i.e., when , as long as and remain finite. However, experiments show that the conductivity reaches a minimum value at the neutrality point due to electron-hole puddles in the charge carrier density across graphene, which are caused by fluctuations of the electrostatic potential in the plane of graphene due to spatial inhomogeneity of the external charged impurities. Chen_2008 (); Tan_2007 () An estimate of may be found according to the SCT theory as , where is referred to as a residual charge carrier density that gives a measure of the width of the plateau near the neutrality point where the conductivity minimum is reached.PNAS_2007 () It was shown that may be found as a solution of an equation involving the square of graphene’s Fermi energy, , and the variance of the fluctuating electrostatic potential in graphene, , that arises from a distribution of immobile external charges,

 (ℏvF)2π¯n=C0(¯n), (17)

where . We note that the SCT theory extends the applicability of the SBT result for the conductivity of graphene down to lower charge carrier densities with typically cm.PNAS_2007 ()

Working in the time-independent regime, we use the 2D spatial FT to express the fluctuating potential in graphene in terms of the 2D FT of the fluctuating charge density as

 δ˜ϕg(q) = ∞∫−∞˜G(0)(q;zg,z)1+e2χs(q)˜G(0)(q;zg,zg)δ˜ρ0(q,z)dz (18) = 2πeqδ˜N(q)ϵs(q), (19)

where is the total dielectric function of the entire system in the static limit. By invoking the translational invariance of the distribution of external charges in the directions of , we may use a general relation,

 ⟨δ˜N(q′)δ˜N(q)⟩=nimpδ(q′+q)S(q), (20)

that allows us to write

 (21)

### ii.4 Statistical description of external charges

It is important to make distinction between the geometric structure of the external particle system and the statistical distribution of the charge density functions for individual particles. Assuming that those two characteristics of the system are statistically independent, the geometric structure may be modeled by using the one- and two–particle distribution functions for their positions

 F1(r,z)=NAf1(z), (22)

and

 F2(r1,r2;z1,z2) = N(N−1)A2f1(z1)f1(z2) (23) ×g(r2−r1;z1,z2),

where describes the distribution of particle positions along the axis and is normalized to one, whereas is the usual pair correlation function. A significant further simplification may be achieved by assuming that the charge densities of individual particles are identically distributed, so that for all . Still, Eqs. (9) and (11) show that the corresponding individual particle form factors generally remain entangled with the dependence of the geometric arrangement of particle positions, unless all the particles reside in the same plane, say .

Accordingly, we first consider a 2D geometric model with , which is commonly used in all theoretical modelings of the effects of correlated charged impurities on the conductivity of graphene.Sarma_2011 (); PNAS_2007 (); Yan_2011 (); Li_2011 () In that case, we find that the auto-correlation function from Eq. (13) may be written as

 (24)

where each particle is characterized by an ”atomic” form factor

 F0(q)=∫dzψ(q,z)˜Δ(q,z−z0), (25)

and

 S2D(q)=1+nimp∫d2reiq⋅r[g2D(r)−1] (26)

is a ”geometric” structure factor that describes the arrangement of external particles in the plane . As regards the corresponding pair correlation (or radial distribution) function , in addition to uncorrelated particles with , we consider two models that contain a single parameter characterizing the inter-particle correlation distance: a step-correlation (SC) model with , where is a Heaviside unit step function, which was often used in the previous studies of charged impurities in graphene,Yan_2011 (); Li_2011 () and the HD model, in which particles interact with each other as hard disks of the diameter . Rosenfeld_1990 ()

There are several advantages to using the HD model over the SC model. First, the former model is based on a Hamiltonian equation for the thermodynamic state of a 2D fluid with a well-defined pair potential between impurities, whereas the latter model is an ad hoc description of the impurity distribution, made-up for simple, analytic results. That is not to say that the SC model is poor at capturing the interesting results in the conductivity of graphene with correlated impurities.Yan_2011 (); Li_2011 () However, from Eq. (16) it is obvious that, with , the initial slope of is strongly influenced by the limiting value of the structure factor as , that is, by the value of via Eq. (24). It is well known that is related to the isothermal compressibility of a 2D fluid,Hansen_1986 () which may be expressed as a function of the packing fraction defined by . Thus, is a key measure of performance of the two models. It was recently shown by Li et al.Li_2011 () that the SC model gives reliable results for packing fractions by comparing the analytical result for the 2D structure factor in that model, , with a numerically calculated structure factor of a hexagonal lattice of impurities. However, the analytical limit shows that the SC model already breaks down for because the corresponding compressibility becomes negative at higher packing fractions. On the other hand, it was recently shown that the interaction potential between two point ions near doped graphene is heavily screened and, moreover, exhibits Friedel oscillations with inter-particle distance, giving rise to a strongly repulsive core region of distances on the order of that resembles the interaction among hard disks with diameter .Radovic_2012 () Therefore, we may estimate that the packing factor could reach values on the order that may not always be negligibly small, necessitating the use of a model that goes well beyond the SC model, at least for systems of adsorbed alkali-atom submonolayers on graphene.Yan_2011 () In that respect, we note that various parameterizations of the HD model extend its applicability to include phase transitions in a 2D fluid as a function of the packing fraction, Mak_2006 () even going up about , corresponding to a crystalline closest packing where hard disks form a hexagonal structure in 2D.Guoa_2006 () In this work, we use a simple analytical parametrization for the 2D structure factor in the HD model, , provided by RosenfeldRosenfeld_1990 () (see Appendix B) which works reasonably well for packing fractions up to about , just near the freezing point of a 2D fluid.

Regarding the structure of individual charged particles within the 2D geometric model, we study a few specific examples. First we consider a point particle of charge that carries a dipole moment with the density function

 Δp(R)=(Z−D⋅∇R)δ(R), (27)

where is an effective dipole length and is a 3D delta function, which gives a form factor from Eq. (25) as

 Fp(q)=(Z+iq⋅D∥)ψ(q,z0)+D⊥∂ψ(q,z)∂z∣∣∣z=z0, (28)

where and are the effective dipole lengths in the directions parallel and perpendicular to graphene, respectively. We note that, having in mind that the first two terms in the right-hand side of Eq. (24) represent the variance of the form factor , all of the three parameters of the point particle model, namely, , and may exhibit fluctuations about their respective means (with the mean due to the presumed isotropy), as well as mutual cross-correlations. In addition, assuming to be small enough, the perpendicular dipole moment component may depend on the local electrostatic field according to , where is an effective dipole polarizability near graphene.

We also consider a cluster of uniformly distributed charge within a disk of radius parallel to graphene with

 Δcl(R)=ZπR2clU(Rcl−r)δ(z), (29)

giving

 Fcl(q)=2ZqRclJ1(qRcl)ψ(q,z0), (30)

where is a Bessel function of order one. We limit our considerations to cases with , validating the perturbative treatment of charge carrier scattering on such clusters,Katsnelson_2010 () and we also assume to avoid the interference in scattering patterns from neighboring clusters.

On the other hand, it is of interest to explore the effects a fully -dependent geometric structure of particle positions in 3D, with arbitrary distribution function and the pair correlation function that depends on the coordinates, . In this case, we only consider point charges with and obtain the auto-correlation function from Eq. (13) as

 S(q) = ∫dzf1(z)ψ2(q,z)+∫dzf1(z)ψ(q,z) (31) ×∫dz′f1(z′)ψ(q,z′)[S3D(q;z,z′)−1],

where partial structure factor in the 3D case is defined by

 S3D(q;z,z′)=1+nimp∫d2reiq⋅r[g3D(r;z,z′)−1]. (32)

Any realistic modeling of the 3D pair correlation function in the presence of charged graphene is beyond the scope of the present study, so we only consider uncorrelated point charges with , and focus on the effect of their distribution over the depth . In a first study of this type, we only consider the case for a homogeneous distribution of point charges throughout a dielectric slab of finite thickness . In Appendix B we also provide a result for semi-infinite region () based on a pair correlation function for a bulk one-component plasma (OCP) with the volume density of charged particles , which may be of interest in future work.

## Iii Results and discussion

In this section, we study several special configurations of graphene with the surrounding dielectrics by using the electrostatic GF, which is derived in Appendix A for a three-layer structure of Fig. 1, defined on the intervals , and along the axis that are characterized by the relative bulk dielectric constants with , respectively.

In Fig. 2 we consider a two-layer structure that consists of a semi-infinite SiO substrate ( with ) and a semi-infinite layer of air ( with , or with ) with graphene placed right on their boundary at . We show the dependence of graphene’s conductivity on its average charge carrier density for a planar distribution of charged impurities with fixed and no dipole moment, having the areal number density cm, which are all placed a distance away from graphene. We show the results for several values of the correlation length among the impurities, which are obtained by using the SC and the HD models for their 2D structure factor, and note that the SC model only yields physical results for nm for the given value of . In addition to the case of point-like impurities being placed directly on graphene ( and ), we also show in Fig. 2 the effects of point-like impurities embedded at nm inside the SiO substrate, as well as disk-like impurities with fixed radius nm placed on graphene ().

As regards the effects of finite and , one notices in Fig. 2 that they both contribute to an increase in the slope of conductivity at higher values, as expected, where they even give rise to a super-linear dependence of conductivity on for smaller values of the correlation length . (Note that the case of uncorrelated disks with is somewhat unphysical as the disks are allowed to overlap.) However, the effects of finite and are relatively weak and only affect quantitative details of conductivity at higher , whereas comparison among the insets in Fig. 2 shows that their effects are barely noticeable at cm.

The most prominent effect in Fig. 2 is a strong increase of the initial slope of conductivity as a function of (and hence an increase in mobility of graphene, ) at low values of as the correlation length increases. One notices from the insets in Fig. 2 that the initial slopes from the SC model are higher than those from the HD model for the same value of because , but the latter model permits the use of much larger values of than the former model, hence giving rise to rather large initial slopes of the conductivity at the largest packing fractions shown. (Notice that the case with a maximum packing fraction of that is shown in Fig. 2 is still well within the interval of confidence for the HD model used here.Rosenfeld_1990 ()) As the charge carrier density increases, the conductivity shows a sub-linear dependence on that becomes more pronounced as the correlation length increases. In the case of and the sub-linear dependence occurs for all , whereas in the cases of finite or values the sub-linear dependence may even overcome the opposite effect of super-linear dependence for sufficiently large s. For the largest value shown in Fig. 2, the sub-linear behavior even gives rise to a pronounced saturation effect in the conductivity of graphene with increasing , which is sometimes observed in experiments.Tan_2007 (); Yan_2011 () Thus, high packing fractions that result from long correlation distances among the charged impurities can give rise to both higher initial slope of conductivity at lower and a more pronounced sub-linear dependence of conductivity at higher with the HD model than those that can be achieved with the SC model. We pause to discuss those two effects in some detail.

Various models that attempt to reproduce the experimental dependence of graphene’s conductivity on its charge carrier density use the areal density of charged impurities as free parameter to fit the slope of conductivity in the range of values where that dependence is found to be predominantly linear. Ignoring the relatively narrow region of values around zero where the conductivity of a nominally neutral graphene reaches a minimum, one sees that Eq. (16) implies a linear dependence of conductivity in the form when , where is constant when the dielectric media are semi-infinite. For a system of uncorrelated impurities that may be described as a 2D gas, one simply finds because . However, when impurities are strongly correlated, one should consider their number to be a random variable because different samples of graphene flakes with fixed area may cover different regions of a much larger area of the substrate plagued by varying concentrations of impurities. Then, the impurity density should be defined in terms of the average number of impurities covered by the graphene flake, . On the other hand, the long wavelength limit of the structure factor may be expressed as the ratio , where the numerator is the variance in ,Hansen_1986 () with being the fluctuation in the number of impurities that are covered by the graphene flake. Therefore, from the statistical point of view, the limit of the SBT conductivity should be reinterpreted as , where we define to be an effective density of impurities rather than the average density. In general, unless is Poisson distributed, i.e., the impurities behave as an ideal 2D gas. Clearly, the distinction between and should be borne in mind when attempting to use as a fitting parameter in modeling the slope of graphene’s conductivity in the presence of a liquid-like distribution of charged impurities.

On the other hand, the sub-linear dependence of graphene’s conductivity on at large doping densities is often modeled by combining the scattering processes of its charge carriers on both charged impurities and short-ranged impurities via the Matthiessen’s rule.Yan_2011 () However, the density of atom-size defects in graphene that could give rise to short-range scattering is extremely low due to the structural and compositional resilience of graphene’s atomic lattice, so that ”the source of the proposed weak short-range scattering is mysterious.”Yan_2011 () Another contender for the explanation of the sub-linear conductivity is the resonant scattering model that invokes the existence of bound-state resonances in the electron bands due to chemisorbed molecules on graphene.Ferreira_2011 () However, the fact that graphene is chemically inert also makes this mechanism unlikely in most situations. On the other hand, it was recently shown that the charge carrier scattering on charged impurities in a substrate may also give rise to the sub-linear behavior of conductivity in highly doped graphene in the presence of a strong spatial correlation among the impurities.Yan_2011 (); Li_2011 () Noting that the sub-linear behavior was demonstrated in simulations based on the SC model with small packing fractions,Yan_2011 (); Li_2011 () we follow the same idea and suggest that, by being able to go to much larger packing fractions in the HD model than in the SC model, one may include large enough values of in simulations that could even give rise to saturation of graphene’s conductivity at high enough charge carrier densities, thus eliminating the need to invoke the existence of resonance scatterers or atom-size defects in graphene. Namely, one may verify that, with increasing packing fraction the structure factor develops a very pronounced peak at the wavenumber corresponding to the first coordination shell due to the nearest neighbors.Rosenfeld_1990 (); Guoa_2006 () So, from Eq. (16) it follows that a relatively sudden increase in the value of the integral over may be expected with the HD model when surpasses the value , causing a slowdown in the increase of when that is reminiscent of the saturation in conductivity. For example, in the case of the largest correlation distance shown in Fig. 2, nm, one finds that a strong saturation of the conductivity indeed occurs at about cm.

In Fig. 3 we consider the same configuration of single-layer graphene atop a semi-infinite SiO substrate with a semi-infinite layer of air above it as in Fig. 2, and attempt to model the experimental data for conductivity versus charge carrier density from Ref. Tan_2007 () by using the HD model for a 2D distribution of point charges with . We select two graphene samples from Ref. Tan_2007 () labeled K17 and K12, which both exhibit sub-linear behavior with increasing , with K17 being symmetric and K12 showing an electron-hole asymmetry (i.e., asymmetry with respect to the sign of ). The physical mechanism(s) that occasionally give rise to this kind of asymmetry in graphene are still unclear, so we explore here the possibility that the presence of the perpendicular component of dipole moment in each impurity, , may give rise to a sizeable asymmetry, as that seen in Fig. 3 for the sample K12. We assume , where is the effective polarizability and is the total perpendicular electric field near graphene. Assuming to be small enough, we may neglect mutual depolarization among the impurities and simply write , with being positive (negative) for electron (hole) doping of graphene.Maschhoff_1994 () The two samples were fitted in Ref. Tan_2007 () by assuming that the impurities reside in graphene () and are uncorrelated, and the optimal linear symmetric fits were found with cm for K17 and with cm for K12. We also assume the impurities to lie in graphene (), and we use , and as fitting parameters. In the case of the symmetric K17, the best fit is found for cm with nm () and , whereas for the asymmetric case of K12 the best fit is found for cm with nm () and Å. Both fits obtained with the HD model in Fig. 3 are quite satisfactory as far as the sub-linear behavior of conductivity is concerned, and the relatively large values of packing fractions used in both cases suggest the necessity of using the HD rather than the SC model. On the other hand, a good fit in the asymmetric case can only be achieved with a rather large value of , which indicates that the dipole mechanism may not be the primary cause of the electron-hole asymmetry in conductivity, at least for the experimental setting of Ref. Tan_2007 () However, we note that the effective polarizability of a single impurity may be significantly increased by the presence of a nearby conducting surface.Maschhoff_1994 ()

In Fig. 4 we consider a structure that consists of a dielectric material of finite thickness (we choose HfO with ) and a semi-infinite layer of SiO (either with , or with ) with graphene placed right on their boundary at . This configuration may represent the physical situation where single-layer graphene sits on a thick SiO substrate (with typically nm) and is top-gated through a thin layer of HfO (with nm). We show the dependence of the conductivity on charge carrier density for several model distributions of point charge impurities in the HfO layer with fixed and no dipole moment, having the areal number density cm. We consider a homogeneous 3D distribution of uncorrelated charges throughout the HfO, which extends up to a distance from graphene, as well as a 2D planar distribution placed in HfO a distance away from graphene, with both uncorrelated () and correlated ( nm, ) charges that are described with the HD model.

One notices in Fig. 4 that finite thickness exhibits strong effects on conductivity, both in quantitative and qualitative aspects, which are dependent on the underlying structure of charged impurities. First noted is that the overall conductivity is generally increased compared to that seen in Figs. 1 and 2, which is expected due to the more efficient screening of charged impurities by a high- material such as HfO. Moreover, the conductivity is seen to increase with decreasing for all in the 2D cases and only for lower in the 3D case, which may be explained by the more efficient screening of impurities due to the proximity of a metal gate. Furthermore, the conductivity is larger in the 3D case than in the corresponding uncorrelated 2D case because the same number of impurities is spread over larger distances from graphene so that the resulting scattering potential in graphene is weaker. As regards the distance , one notices similar trends as in Fig. 2, namely, a finite increases both the value of conductivity and its slope (i.e., mobility) in both 3D and 2D models. However, as regards the effects of finite correlation length in the 2D models with finite , one sees little evidence to the increase in the initial slope of conductivity at lower , in contrast to the trends seen in Fig. 2, whereas saturation of conductivity at higher seems to get stronger than in Fig. 2 as decreases. In fact, for the shortest thickness of nm for both and nm, this saturation turns into a broad maximum of conductivity around cm, followed by a still broader minimum at higher values.

One remarkable feature seen in Fig. 4 is that the conductivity generally does not vanish in the SBT limit when for finite , but rather reaches a minimum value . This minimum may be easily estimated for by using the limiting form of the background dielectric constant when in Eq. (16), which then gives

 σ(0)=(ϵ1πrsL)2e2/(2h)nimpS(0)=4vFπrsC2LnimpS(0), (33)

where for the 3D case, for the uncorrelated 2D case, and for the correlated 2D case in the HD model. In the second expression for in Eq. (33) we emphasize that the minimum conductivity in the SBT limit for neutral graphene is governed by the geometric capacitance per unit area, , of the dielectric with finite thickness used in top-gating the graphene.

Finally, one notices in Fig. 4 that, as the thickness increases in the 3D case, the conductivity gains quite strong super-linear dependence with increasing . This dependence may be estimated by considering Eq. (16) in the limit of large but finite , such that . In that case, the background dielectric constant becomes , whereas the 3D structure factor, which is determined by the first term in Eq. (31), goes as , so that Eq. (16) gives , where is the volume density of charge impurities. We note that this behavior of conductivity in graphene at large is a consequence of the 3D nature of a distribution of uncorrelated charges that gives rise to the special form of structure factor, . The lack of experimental observations of such super-linear dependence of conductivity in graphene should not be taken as evidence to rule out the role of 3D distributions of impurities, because both the correlation among impurities, as that described in the Appendix B for a OCP, as well as their clustering close to graphene seem to be capable of eliminating the super-linear dependence.

In Fig. 5 we consider a three-layer structure that consists of a HfO layer () with finite thickness , a layer of air () of thickness nm, and a semi-infinite layer of SiO (), with graphene placed in the air at nm, midway between the two dielectrics. This configuration is similar to that in Fig. 4 with graphene sandwiched between the HfO and SiO dielectrics, but we introduce in Fig. 5 gaps of air of equal thickness 0.3 nm on both sides of graphene. We investigate the effects of finite thickness on the mobility of graphene, , as a function of charge carrier density for a 2D planar distribution of uncorrelated point charges with and no dipole moment, having the areal density cm. We consider three configurations, with the impurities placed either (A) on graphene () or (B) on the surface of the HfO layer a distance nm away from graphene, both in the presence of the 0.3 nm gaps, as well as the case (C) from Fig. 4(c) having zero gaps between graphene and the HfO and SiO dielectrics with the 2D distribution of uncorrelated charges placed on graphene (). One may see in Fig. 5 that the mobility generally increases with decreasing within each of the three configurations, (A), (B) and (C), but that there are remarkable differences between them in the magnitude of the mobility and its dependence on . In the configurations (A) and (C) with charge impurities placed on graphene, the mobility generally decreases with increasing , whereas in the configuration (B) with the impurities placed on the surface of the HfO layer with a finite gap relative to the graphene, the mobilities with higher values pass through a minimum at a low value and further increase as increases. Moreover, the magnitudes of the mobility with equal values are seen in Fig. 5 to increase in the order of configurations (A)(C)(B), which is also the order of increasing spread of the curves with different values within each configuration. Finally, it is interesting to notice that differences between the magnitudes of the mobility in the three different configurations with become diminished as decreases.

One may conclude from Fig. 5 that the existence of a finite gap between graphene and the nearby dielectric, as well as the precise location of impurities within that gap (with the extreme positions being on graphene and on the surface of the dielectric) both have decisive influences on the mobility. Noting that the configuration (A) with impurities on graphene in the presence of finite gaps was considered in Ref.Ong_2012 (), it is remarkable how closing the gaps increases the magnitude of the mobility and increases the spread of its values for different values, whereas moving the impurities to the surface of a HfO layer in the presence of finite gaps further accentuates those two effects, and even gives rise to a non-monotonous dependence of the mobility on for thicker HfO layers. While the role of the distance of impurities from graphene was discussed in detail for the case of zero gaps,PNAS_2007 () one may conclude from our analysis that the size of the gap(s) between graphene and the nearby dielectric(s) plays equally important role in modeling the conductivity of graphene in a broad range of charge carrier densities.

We next turn to studying the conductivity minimum as due to the presence of electron-hole puddles by using Eqs. (17) and (21) based on the SCT theory.PNAS_2007 () We only consider a 2D planar distribution of point charges with having no dipole moment and note that, unlike the integral in Eq. (16) for conductivity, in order to render the integral in Eq. (21) convergent one must assume that charged impurities are placed a finite distance away from graphene.

In Fig. 6 we consider a configuration similar to that in Fig. 2, with a semi-infinite SiO substrate ( with ) and a semi-infinite layer of air ( with ), with graphene placed in the air at a distance above SiO. We show in the main panel of Fig. 6 the dependence of the variance of the potential in the plane of graphene from Eq. (21) for a 2D distribution of charged impurities with density cm that are placed in/on SiO at a fixed distance nm below graphene. Specifically, we explore the effects of the size of the gap between graphene and the SiO substrate by considering both the zero gap case with (impurities embedded at the depth of 0.3 nm inside SiO) and the finite gap case with nm (impurities placed on the surface of SiO). In addition to considering uncorrelated impurities, we use a finite correlation length of nm () allowing us to compare in the main panel the effects of the SC and the HD models on . In the insets of Fig. 6, we show the dependence of the residual charge carrier density and the corresponding minimum conductivity on for both the HD and the CS models, in the presence of both zero and finite gaps.