Clustered impurities and carrier transport in supported graphene
CClustered impurities and carrier transport in supported graphene
N. Sule, S. C. Hagness, and I. Knezevic ∗ Department of Electrical and Computer Engineering, University of Wisconsin–Madison, Madison, WI 53706-1691, USA (Dated: October 9, 2018)We investigate the effects of charged impurity distributions and carrier-carrier interactions on electronic trans-port in graphene on SiO by employing a self-consistent coupled simulation of carrier transport and electrody-namics. We show that impurity clusters of characteristic width 40–50 nm generate electron–hole puddles ofexperimentally observed sizes. In the conductivity versus carrier density dependence, the residual conductivityand the linear-region slope are determined by the impurity distribution, and the measured slope can be usedto estimate the impurity density in experiment. Furthermore, we show that the high-density sublinearity in theconductivity stems from carrier-carrier interactions. PACS numbers: 72.80.Vp, 81.05.ue, 72.10.-d
I. INTRODUCTION
Graphene, a single sheet of carbon with a honeycomb lat-tice, is a two-dimensional (2D) material whose high carriermobility and carrier density tunable by a back gate [1–5] makeit attractive for electronic device application [6–13]. Large-area, good-quality graphene is commonly fabricated by chem-ical vapor deposition (CVD) on metal substrates [7, 14, 15],followed by transfer onto insulating substrates using poly-mers, such as poly-dimethyl siloxane (PDMS) or poly-methylmethacrylate (PMMA). An important concern with these pro-cessing methods is the contamination of graphene with or-ganic molecules [16], residues of the transfer polymer andmetal ions [17], or charged impurities trapped in the support-ing substrate [18].Impurities near graphene are believed to be responsiblefor several observed transport properties. Spatial inhomo-geneities in the carrier density, known as electron–hole pud-dles, are formed due to the presence of charged impurities inthe substrate [19–21]. The charged impurities and the result-ing electron–hole puddles have been linked to the observednon-universal minimum conductivity (also known as residualconductivity) of graphene close to the Dirac point [22]. How-ever, high-resolution scanning tunneling microscopy (STM)studies [23] have shown that electron–hole puddles near theDirac point are typically
20 nm in diameter, while theoret-ical calculations using a random charged impurity distribu-tion near graphene result in electron–hole puddle sizes of onlyabout [21]. This evidence suggests that the underlyingcharged impurities may be clustered. It has also been shownthat PMMA and metal ion residue can persist on graphenesamples even post-annealing [17] and transmission electronmicroscopy (TEM) images [17] show that the residue is notuniformly distributed, but forms clusters. Furthermore, theformation of gold clusters has been shown to affect the elec-tron mobility in graphene [24].The linear dependence of conductivity, σ , on carrier den-sity, n , has been attributed to carrier scattering with chargedimpurities [25, 26]. However, experimental measurements ∗ Email: [email protected] distinctly display a sublinear σ ( n ) dependence away fromthe charge-neutrality point [27–29]. The origin of the sub-linear σ ( n ) behavior is still under debate: it has been ascribedto different physical mechanisms, such as electron scatteringwith residual organic molecules [30] or the effect of spatialcorrelations in the distribution of the charged impurities neargraphene [31, 32].In this paper, we employ numerical simulation of coupledcarrier transport and electrodynamics to investigate the roleof carrier-carrier and carrier-ion Coulomb interactions on theroom-temperature, low-field transport in graphene on SiO ,with focus on the effect of impurity clustering. We solve theBoltzmann equation for carrier transport by using the ensem-ble Monte Carlo (EMC) method, coupled with the electro-dynamics solver that incorporates the finite-difference time-domain (FDTD) solution to Maxwell’s curl equations andmolecular dynamics (MD) for short-range carrier-carrier andcarrier-ion interaction. We show that clustered distributions ofimpurities with an average cluster size of 40–50 nm result inthe formation of
20 nm -wide electron–hole puddles, the sizeobserved in several experiments [20, 23, 33]. We demonstratethat the sublinear behavior of conductivity at high carrier den-sities, which becomes more pronounced with decreasing im-purity density [27, 28], stems from short-range carrier-carrierinteractions. Also, we show that the linear portion of the con-ductivity versus carrier density curve is governed by carrier-ion interactions, with the slope and the residual conductivitydependent on both the sheet impurity density and the impuritydistribution. We characterize the dependence of the conduc-tivity slope on the impurity density for uniform random andclustered distributions, which can be used to estimate the im-purity density in experiment.This paper is organized as follows: In Sec. II, we overviewthe EMC, FDTD, and MD techniques and their coupling (Sec.II A), and describe the generation of a clustered impurity dis-tribution (Sec. II B). In Sec. III, we discuss electron–holepuddle formation (Sec. III A), the role of impurity cluster-ing in low-carrier-density transport (Sec. III B, sublinearityin conductivity and its connection to the short-range carrier-carrier interaction (Sec. III C, and how to estimate impuritydensity from the linear-region conductivity slope (Sec. III D).We conclude with Sec. IV. a r X i v : . [ c ond - m a t . m e s - h a ll ] M a y II. THE SIMULATION FRAMEWORK
Our goal is to accurately simulate room-temperature elec-tron and hole transport in supported graphene with chargedimpurities in the substrate, with focus on impurity clusteringand Coulomb interactions (carrier–ion and carrier–carrier).Experiments have shown that charged impurities are the dom-inant source of disorder in supported graphene [27, 29, 34].As shown by Kohn and Luttinger [35], the Boltzmann trans-port equation can be derived quite generally from the density-matrix formalism for electrons in the presence of dilute uncor-related charged impurities. Indeed, at moderate carrier den-sities in graphene, transport is diffusive and well-describedby the Boltzmann transport equation, with the conductivitybeing linear in the carrier density owing to carrier–ion inter-actions [2, 25]. In the vicinity of the Dirac point, the aver-age carrier density can be considerably lower than the impu-rity density and charge inhomogeneities referred to as pud-dles govern transport. However, the effective medium theory[2, 21, 36, 37] argues that, while the average carrier densityfor the entire sample may be low, carrier density within anindividual puddle is fairly uniform and on the order of theimpurity sheet density, and the Boltzmann transport pictureremains applicable [22].Therefore, we will assume the diffusive transport regime,captured through the Boltzmann transport equation, through-out the range of carrier and impurity densities and distribu-tions considered here. In fact, we find that clustered impuri-ties result in sizeable puddles with the carrier density that isnearly uniform and is of order the impurity density, in agree-ment with the effective medium theory. The assumption ofdiffusive transport is further strengthened by the fact that weare at room temperature and working with macroscopic sam-ples with size greater than the mean free path [22].
A. EMC/FDTD/MD for Graphene on SiO In order to simulate diffusive carrier transport and elec-trodynamics in supported graphene, we employ a cou-pled EMC/FDTD/MD technique [38, 39]. In a nutshell,EMC solves the Boltzmann transport equation, FDTD solvesMaxwell’s curl equations, while MD accounts for the interac-tion of charges when very close to one another. The coupledEMC/FDTD/MD technique was successfully used to calcu-late the high-frequency conductivity of bulk silicon, with verygood agreement to experimental data [38, 40]. Below, webriefly describe the key elements of the constituent techniquesand refer the interested reader to references [38] and [39] forextensive computational detail.EMC is a stochastic numerical technique widely used forsolving the Boltzmann transport equation [41]. In EMC, alarge ensemble (typically of order ) of carriers is trackedover time as they experience periods of free flight interruptedby scattering events. Free-flight duration, the choice of the re-laxation mechanisms, and carrier momentum direction post-scattering are sampled stochastically according to appropri-ate distributions. During free flight, carriers interact with GrapheneSiO Air Impurity ioncluster
FIG. 1. A schematic of the simulated structure, depicting a mono-layer of graphene on an SiO substrate, with air on top. Clusters ofsubstrate impurities near the graphene sheet are also shown. the local electromagnetic fields via the Lorentz force, (cid:126)F = q ( (cid:126)E + (cid:126)v × (cid:126)B ) , where q , (cid:126)v , (cid:126)E , and (cid:126)B are the carrier charge,carrier velocity, electric field, and magnetic flux density, re-spectively. The fields are calculated using the electrodynam-ics solver that includes the FDTD and MD components. Theevolution of physical properties of interest, such as the car-rier average drift velocity or kinetic energy, are calculated byaveraging over the ensemble.The FDTD method [42] is a popular and highly accurategrid-based technique for solving Maxwell’s curl equations. InFDTD, Maxwell’s equations are discretized in both time andspace by centered differences using the fully explicit Yee al-gorithm [43]: the components of electric and magnetic fields, (cid:126)E and (cid:126)H , are spatially staggered and solved for in time usinga leapfrog integration method, where the (cid:126)E and (cid:126)H updatesare offset by half a time step, yielding second order accuracyof the algorithm. The spatial grid cell size and the time stepin FDTD must be chosen such that they satisfy the Courantstability criterion [42].Carrier motion in EMC gives rise to a current density, (cid:126)J ,which acts as a field source in FDTD; in turn, fields calcu-lated by FDTD affect the motion of carriers in EMC. How-ever, grid-based methods such as FDTD do not account forfields on the length scales shorter than a grid-cell size [44],so we use the MD technique [45, 46] to calculate the short-range, sub-grid-cell fields stemming from pair-wise Coulombinteractions among electrons, holes, and ions. Carrier–ion,direct carrier–carrier, and exchange carrier–carrier (electron–electron and hole–hole) interactions are included [38, 39].The simulated structure, shown in Figure 1, consists ofa monolayer of graphene placed on a silicon-dioxide sub-strate that contains charged impurities. On the four verti-cal planes that bound the simulation domain perpendicular tothe graphene layer, we apply periodic boundary conditions tothe fields and carrier momenta. The top and bottom planesthat bound the simulation domain parallel to the graphenelayer are terminated using convolutional perfectly matchedlayer (CPML) absorbing boundary conditions [42]. In theFDTD/MD electrodynamic solver, the monolayer of grapheneis defined by one plane of grid points with a dielectric constantof . , while the grid points above and below this plane aregiven dielectric constants of (air) and . (SiO ), respec-tively.We assume that the Fermi level and carrier density ingraphene can be modulated by a back gate, located at thebottom of the SiO substrate. For a given Fermi level andtemperature, the electron and hole densities are given by n = n i F ( η ) / F (0) and p = n i F ( − η ) / F (0) , respectively.[47]Here, n i = π (cid:16) kT ¯ hv F (cid:17) , η = E F /kT , and F j ( η ) is theFermi integral of order j . E F is the tunable Fermi level and v F = 10 cm / s is the Fermi velocity in graphene on SiO [48]. The carrier ensemble in the 2D plane of graphene, com-prising electrons and holes, is initialized by using randomnumbers to assign a position, momentum, charge, and spinto each carrier, taking into account the appropriate statisticalprobabilities. For the calculation of the grid-based charge den-sity, carriers localized throughout the simulation domain areassigned to the grid using the cloud-in-cell method [49]. Theinitial electric field distribution is calculated by solving Pois-son’s equation using the successive-over-relaxation method[50]. We use the tight-binding Bloch wave functions [51]to calculate the electron-phonon scattering rates in graphene,accurately reproducing the rates from first-principles calcula-tions [52], and to compute the electron-surface-optical (SO)phonon scattering rates [53]. The scattering rates for holes areassumed to be the same as those for electrons. These initial-ization steps are followed by a time-stepping loop in whichEMC and FDTD/MD source each other and which terminatesonce a steady state is achieved, as identified by the saturationof the ensemble-averaged carrier velocity and energy. B. Generating Clustered Impurity Distributions in theSimulation
In order to capture the influence of charged impurities onelectron and hole transport in graphene, we generate differ-ent impurity distributions throughout the SiO substrate. Thetype and charge of relevant impurities vary with the process-ing details [54]; for simplicity, we use generic impurity ionswith unit positive charge. The impurity ions in the simulationare distributed in three dimensions (3D); however, impuritiesin the graphene literature are typically described via a cumu-lative sheet density, N I , in units of cm − . For a generated 3Ddistribution of ions, the sheet density is obtained by integrat-ing over a depth equal to r d , where r d represents the effectivesize of an impurity ion in the MD calculation (see Willis et al. [38, 40] for more details), followed by averaging over the to-tal depth of the 3D distribution. r d is typically between . and . . We have observed that charged impurities placeddeeper than
10 nm do not significantly affect carrier transportfor reasonable impurity sheet densities ( N I < cm − ).The problem of positioning individual impurities in 3D toachieve a predetermined cluster size distribution is related to3D Voronoi tessellation [55–57]. Here, we have developeda relatively simple algorithm that enables us to generate anapproximately Gaussian distribution of individual impurities y [ n m ] y [ n m ] (a) (b)(c) (d) L c [nm] 0 10 20 30 40 50 60 [ n m ] λ (cid:1) c y [ n m ]
50 0 5050050 y [ n m ] fi tSACFL c
10 nm 50 nm0 50 100 150 200 250x [nm]00.20.40.60.81 A r b . un i t s (e) FIG. 2. Examples of clustered impurity distributions generated forclustering parameters (a) L c = 10 nm and (c) L c = 50 nm . Thecorresponding normalized SACF are shown in (b) and (d), respec-tively. The average impurity cluster size, λ c , is estimated from theFWHM (yellow ring) of the SACF. (e) Gaussian fits (orange and bluecurves) to the SACFs from (b) and (d) (red and purple dots, respec-tively). (f) λ c versus L c . Each data point corresponds to the averageof simulation runs for a given L c , while the error bars denote thestandard deviations. The dashed line is a quadratic fit to guide theeye ( λ c = 0 . L + 0 . L c + 22 . ). starting from a single numerical parameter, L c , which we re-fer to as the clustering parameter . For L c = 0 , we distributeall the impurity ions stochastically according to a uniform ran-dom distribution. For a non-zero L c , we generate N c = A/L impurity clusters, where A is the 2D area of the graphenelayer in the simulation. To initialize the positions of individ-ual impurities, we first distribute the centers of the N c clus-ters stochastically. Secondly, we pick the characteristic sizeof each individual cluster from a uniform random distributionbetween L c / and L c / , the average being L c / . Next, wedistribute individual impurity ions around each cluster centerfollowing a Gaussian distribution whose standard deviationequals half of the cluster size. Examples of clustered impuritydistributions are shown in Figure 2a ( L c = 10 nm ) and Figure2c ( L c = 50 nm ), with the corresponding spatial autocorrela-tion functions (SACFs) depicted in Figures 2b and 2d, respec-tively. As shown in Figure 2e, normalized Gaussians (orangeand blue curves) fit the SACFs (red and purple dots) well.Moreover, the full-width at half-maximum (FWHM) of theSACF agrees well with the correlation length extracted fromthe Gaussian fits. Henceforth, the FWHM of the impurity-distribution SACF will be referred to as the average impuritycluster size and denoted by λ c . Figure 2f presents λ c ver-sus L c . Each data point in Figure 2f represents the averageof fourteen slightly different impurity ion configurations ob-tained stochastically for a given value of L c (ranging from 0to
60 nm in the increments of ) and the error bars on thedata points denote the standard deviations.It is important to note that λ c is conceptually different fromthe correlation length r used by Li et al. [31] r representsthe extent to which impurity ions can interact with one an-other and diffuse; as a result, a larger r results in an impuritydistribution that is more spread out than clustered. In con-trast, a larger λ c (stemming from a larger L c , see Figure 2f)represents a more clustered distribution. III. RESULTS AND DISCUSSIONA. Formation of electron–hole puddles
Figure 3 shows the formation of electron–hole puddles inthe presence of clustered impurity distributions. We simu-late carrier transport at room temperature, for the Fermi levelat the Dirac point ( E F = 0 ), and without external fields.The initial positions of the charge carriers in the simulationare generated randomly based on a uniform spatial distri-bution and the calculated electron and hole sheet densities n = p = 8 × cm − . As the simulation progresses,carriers move and scatter until a steady state is reached. Themotion of carriers under the influence of the other charges inthe domain (the clustered ions as well as other carriers) resultsin a charge redistribution and the formation of electron–holepuddles. The average electron–hole puddle size is estimatedfrom the FWHM [21] of the SACF of the carrier density dis-tribution. In Figures 3a and 3b, we contrast the carrier densitydistributions that stem from the underlying uniform random( L c = 0 , λ c = 22 nm ) and clustered impurity distributions( L c = 50 nm , λ c = 46 nm ). The corresponding SACFs ofthe carrier density are shown in Figs. 3c and 3d; the corre-sponding average electron–hole puddle sizes, estimated fromthe FWHM of these SACFs, are λ p = 5 nm and
20 nm , re-spectively. These examples show a very significant differencein the sizes of electron–hole puddles that result from randomand clustered impurity ion distributions. Figure 3e shows theaverage electron–hole puddle size, λ p , as a function of the av-erage impurity cluster size λ c . Different simulation runs forthe same n , p , and N I produce slightly different puddle andimpurity cluster sizes owing to the stochastic nature of theimpurity position initialization and the EMC routine. There-fore, each data point in Figure 3e represents the average offourteen simulations for a given value of L c (ranging from 0to
60 nm in the increments of ) and the error bars onthe data points denote the standard deviations. A uniform ran-dom impurity distribution results in an average puddle size of − −
1x 10 −
50 0 50x [nm] − − −
1x 10 y [ n m ] − −
50 0 50x [nm] 0.80.60.40.201 −
50 0 50x [nm] − −
50 0 50x [nm] y [ n m ] − λ , e - h pudd l e s i z e [ n m ] λ , impurity cluster size [nm] c (a) (b)(c) (d)(e) p L c = 10 nm L c = 50 nm20 30 40 50 60510152025 FIG. 3. Carrier density distribution (blue: electrons, red: holes) de-picting the electron–holes puddles formed in graphene at the Diracpoint for (a) uniform random ( L c = 0 , λ c = 22 nm ) and (b) clus-tered ( L c = 50 nm , λ c = 46 nm ) impurity distributions, both withimpurity sheet density equal to × cm − . The average sizeof the electron–hole puddles, λ p , is estimated from the FWHM (yel-low ring) of the normalized carrier-density SACF, shown in (c) and(d), corresponding to the random and clustered impurity distributionsfrom (a) and (b), respectively. The estimated puddle size from (c) is λ p = 6 nm and that from (d) is λ p = 20 nm . (e) Characteristicelectron–hole puddle size λ p as a function of the average impuritycluster size λ c . Each data point corresponds to a single value of L c (swept from 0 to
60 nm in the increments of ) and is theaverage of simulation runs; the error bars denote the standard de-viations. The insets show illustrative impurity distributions, nearlyrandom on the left ( L c = 10 nm ) and highly clustered on the right( L c = 50 nm ). only , while impurity clusters with an average size of40–50 nm give rise to electron–hole puddles with an averagesize of
20 nm , in agreement with experimental observations[20, 23, 33].
B. Role of impurity distribution in carrier transport. Residualconductivity
Next, we examine the effect of random and clustered impu-rity distributions on carrier transport in supported graphene.We calculate the conductivity, σ , as a function of electron den- σ [ e / h ] cm –2 × cm –2 cm –2 uniform random x 10 –2 ] no impurities10 cm –2 × cm –2 cm –2 clustered –2 ](a) (b) FIG. 4. Conductivity of graphene on SiO for (a) uniform ran-dom ( L c = 0 , λ c = 22 nm ) and (b) clustered ( L c = 50 nm , λ c = 46 nm ) impurity distributions, at impurity sheet densities of cm − (triangles), × cm − (squares), cm − (dia-monds), and without impurities (circles). sity n for various spatial formations and total sheet densitiesof impurity ions. The electron density is varied by varyingthe Fermi level, mimicking the effect of a back gate. An ex-ternal dc electric field is applied in the plane of the graphenesheet. The field is introduced using a total-field scattered-fieldincident-wave source condition for a uniform plane wave witha half-Gaussian temporal variation [42]; the magnitude of thesource remains constant once the peak value is achieved. Theconductivity is calculated from σ = (cid:126)J · (cid:126)E/ | (cid:126)E | , where (cid:126)E isthe local electric field and (cid:126)J is the current density. As (cid:126)E and (cid:126)J are noisy, we find σ in the steady state, upon averagingover position and time. In the following simulation results,we have used L c = 0 ( λ c = 22 nm ) for a uniform randomand L c = 50 nm ( λ c = 46 nm ) for a clustered impuritydistribution.In Figure 4, we present σ ( n ) for graphene on SiO atseveral impurity sheet densities, ranging from impurity-freeto N I = 10 cm − , with uniform random and clusteredimpurity distributions. At low impurity densities ( N I < cm − ), the carrier-density dependence of conductivityis nearly the same for the random and clustered impurity dis-tributions, which is not surprising and agrees with the workof Li et al. [31]: with few impurities present, their effecton transport is minor, while carrier interactions with phononsand other carriers dominate. In contrast, at impurity densitieshigher than cm − , uniform random and clustered im-purity distributions result in significantly different σ ( n ) vari-ations. The most significant difference is seen at low carrierdensities, where the conductivity for randomly distributed im-purities increases nearly linearly with increasing carrier den-sity, while that for clustered impurities remains flat. Theslow increase in the conductivity near the charge neutralitypoint has also been observed in experimental measurements[27, 28], notably for samples with considerable impurity con-tamination.In Figures 5a and 5b, we zoom in on the low-density behav-ior of σ ( n ) from Figure 4. The low-density limit of conductiv-ity, σ , known as the residual conductivity, has been observedin experiment [27] and attributed to charged impurity scatter- (cid:1) [ e / h ] n [cm –2 ]no impurities10 cm –2 cm –2 cm –2 clustered no impurities10 cm –2 cm –2 cm –2 uniform random (cid:1) [ e / h ] n [cm –2 ] 6(c) (d) FIG. 5. Conductivity at low carrier densities ( < cm − ) for (a)uniform random ( L c = 0 , λ c = 22 nm ) and (b) clustered ( L c =50 nm , λ c = 46 nm ) impurity distributions with sheet densities of cm − (triangles), × cm − (squares) and cm − (diamonds), and no impurities (circles). Paths of sample carriers ingraphene for (c) uniform random and (d) clustered impurity distri-butions, for N I = 5 × cm − and n = 8 . × cm − ( E F = 0 . ). ing [22]. Here, we see that the value of σ depends on theimpurity sheet density and distribution, with higher impuritydensity and more clustered distributions resulting in a lower σ . We attribute the low- n flattening of conductivity and thelower value of σ for clustered distributions to carrier trap-ping. Figures 5c and 5d depict the paths of sample carriersin graphene with underlying random and clustered substrateimpurity distributions, respectively. A large impurity clustereffectively traps an electron, localizing the electron’s trajec-tory to the cluster vicinity and preventing it from participatingin the current flow. C. Sublinearity in σ ( n ) and carrier-carrier interactions In Figure 6, we examine the role of short-range Coulombinteractions (carrier-carrier and carrier-ion) on dc transport ingraphene on SiO . We account for these effects via the MDpart of the simulation and can selectively turn them on or off tobetter elucidate their importance. Figure 6a presents σ ( n ) forimpurity-free graphene, with MD (circles) and without MD(diamonds); without impurities, MD accounts only for theshort-range, direct and exchange carrier-carrier interactions.We deduce that the sublinearity in σ ( n ) at high carrier den-sities occurs largely due to carrier-carrier interactions: whenwe exclude their short-range component by turning off MD, σ ( n ) becomes nearly linear. Any remaining sublinearity inthe “no MD” results can be attributed to the long-range, di-rect carrier-carrier Coulomb interaction that is captured by theFDTD solver. Carrier-carrier Coulomb interactions do not di-rectly affect conduction (the total momentum of an interacting with MDno e-eno MD uniform random –2 ] x 10 with MDno MD01020304050 σ [ e / h ] –2 ] x 10 no impurities with MDno e-eno MD010203040 clustered –2 ] x 10 a r b . un i t s FIG. 6. Effect of short-range Coulomb interactions, accounted for via MD in the simulation, on transport in supported graphene (a) withoutsubstrate impurities, as well as with (b) uniform random ( L c = 0 , λ c = 22 nm ) and (c) clustered ( L c = 50 nm , λ c = 46 nm ) impuritydistributions. In all three panels, “no MD” indicates simulation results without any short-range interactions. In panel (a), “with MD” denotessimulations with the short-range direct and exchange carrier-carrier interactions included via MD. Inset to (a) shows a representative kineticenergy distribution of electrons with and without carrier-carrier interaction ( n = 7 × cm − , E F = 0 . ). In (b) and (c), “no e-e”indicates results of simulations with short-range carrier-ion interaction but without short-range carrier-carrier interactions, while “with MD”indicates simulations with the full account of all short-range interactions through the coupled EMC/FDTD/MD simulation. Impurity sheetdensity in panels (b) and (c) is × cm − . pair is conserved, as is the pair’s total energy), but redistributethe momentum and energy among the pair and therefore affectthe single-particle distribution function, pushing it towards ashifted Fermi-Dirac distribution [58–61]. The inset to Fig-ure 6a presents the computed distribution of electrons overkinetic energy with and without carrier-carrier interaction forthe electron density of × cm − ( E F = 0 . ). Thiscurve corresponds to g ( E ) f ( E ) , where g ( E ) is the electrondensity of states and f ( E ) is the distribution function, andcarrier-carrier interaction clearly leads to a greater abundanceof higher-energy carriers. Since electron and hole scatteringrates with phonons increase with increasing energy, the redis-tribution of carriers over energy effectively raises the averagecarrier-phonon scattering rate and leads to a reduction in con-ductivity that we observe as the slopeover in σ ( n ) .In Figures 6b and 6c, we plot σ ( n ) for uniform random andclustered impurity distributions with all short-range interac-tions accounted for through MD (circles, “with MD”), withshort-range carrier-ion but without carrier-carrier interactions(triangles, “no e-e”), and without any short-range interactions(diamonds, “no MD”). We have already discussed the low- n region (see Figure 5) and will focus here on the medium-to-high electron density range. In both Figures 6b and 6c, thesheet density of impurities is appreciable ( × cm − ),so carrier-ion interactions govern transport in the medium andthe σ ( n ) dependence is largely linear [25]. Turning off short-range carrier-carrier interactions causes insignificant changeto the slope in either panel, while turning off short-rangecarrier-ion interactions significantly affects the slope. D. Estimating impurity density from the inverse slope of σ ( n ) The slope of the σ ( n ) curves in the linear region is gov-erned by the short-range carrier-ion interactions, and is de-pendent on both the impurity density (Figure 4) and distri- bution (Figures 6b and 6c). As the slope can be accuratelymeasured in experiment, we can use it to indirectly extractthe impurity density and cluster size. In Figure 7, we presentthe EMC/FDTD/MD simulation results for the inverse slopeof σ ( n ) in the linear region as a function of the sheet impu-rity density, with the cluster size as a parameter. The solidmarkers represent the simulation results for uniform random( L c = 0 , λ c = 22 nm , denoted by squares) and clusteredimpurity distributions ( L c = 50 nm , λ c = 46 nm , denotedby triangles). The curves are polynomial fits to guide the eyeand indicate the range of results for different impurity dis-tributions. As we discussed earlier, impurity cluster sizes of40–50 nm correspond to electron–hole puddle sizes obtainedin experiment(see Fig. 3), so it is likely that a reasonable sheetimpurity density estimate can be obtained from the clusteredimpurity curve in Figure 7. As examples, we present the in-verse slopes extracted from several room-temperature mea-surements on supported graphene (a – Ref. [62], b – Ref.[12], c – Ref. [63], and d – Ref. [64]). The interceptsof each inverse-slope horizontal line with the clustered andrandom distribution curves in Fig. 7 indicate an estimate ofthe impurity-density range, with the clustered-curve interceptlikely yielding a good approximate value for N i . Note that themore recent experiments (a from 2012 [62] and b from 2007[12]), which arguably had samples with fewer impurities thanthe early ones owing to the advances in processing, indeedcorrespond to lower sheet impurity densities than the earliermeasurements (c in 2005 [63] and d in 2004 [64]). IV. CONCLUSION
In summary, we have employed EMC/FDTD/MD coupledsimulation of carrier transport and electrodynamics to inves-tigate the effects of carrier-carrier and carrier-ion Coulombinteractions on the transport properties of graphene on SiO , clustereduniform random23456 x 10 I n v e r s e s l ope [ h / ( e · c m ) ] N , Impurity density [cm –2 ] i FIG. 7. Inverse slope of σ ( n ) as a function of the sheet impurity den-sity for graphene on SiO at room temperature. Squares denote theuniform random impurity distribution ( L c = 0 , λ c = 22 nm ), whiletriangles correspond to clustered impurity distributions that wouldgive realistic e-h puddle sizes ( L c = 50 nm , λ c = 46 nm ). The hor-izontal lines (a–d) correspond to the inverse slope values obtained inseveral experiments: a – Ref. [62], b – Ref. [12], c – Ref. [63], and d– Ref. [64]. The N i -range between the intercepts of an inverse-slopehorizontal line with the clustered and random distribution curves (i.e.the range within the lightly shaded area) yields an estimate of the im-purity density range. with focus on the role of substrate impurity clustering. Whilecorrections due to many-particle correlations [19] and coher-ent transport features [65, 66] may play an important role inextremely clean suspended graphene at low temperatures, oursimulations accurately capture the physics of diffusive, room-temperature carrier transport in supported graphene, which isrelevant for device applications. We have shown that clusteredimpurity distributions with an average cluster size of 40–50nm result in the formation of electron–hole puddles with a typ-ical size of
20 nm , comparable to observed values. We havealso demonstrated that high-density clustered impurities leadto carrier trapping and a flattening of the low- n σ ( n ) depen-dence. By selectively controlling the short-range Coulomb in-teractions of the carriers in the coupled EMC/FDTD/MD sim-ulation, we have shown that the sublinear σ ( n ) dependence athigh carrier densities can be attributed to carrier-carrier inter-actions [27, 28]. The slope of the linear-region σ ( n ) relatesto the strength of the carrier-ion Coulomb interactions, andwe have characterized its dependence on the impurity den-sity and distribution. The computed dependence of the linear-region slope of σ ( n ) on the impurity density might be used asa noninvasive technique for estimating the impurity density inexperiment. ACKNOWLEDGMENT
This work has been primarily supported by AFOSR, awardNo. FA9550-11-1-0299. I.K. acknowledges partial support byNSF, award No. 1201311. [1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Rev. Mod. Phys. , 109 (2009).[2] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Rev. Mod.Phys. , 407 (2011).[3] P. Avouris, Nano Lett. , 4285 (2010).[4] A. F. Young and P. Kim, Annu. Rev. Cond. Mat. Phys. , 101(2011).[5] F. Schwierz, Nat. Nano , 487 (2010).[6] F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, Nat. Pho-ton. , 611 (2010).[7] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim,J.-H. Ahn, P. Kim, J.-Y. Choi, and B. H. Hong, Nature , 706(2009).[8] B. J. Kim, H. Jang, S.-K. Lee, B. H. Hong, J.-H. Ahn, and J. H.Cho, Nano Lett. , 3464 (2010).[9] Y.-M. Lin, C. Dimitrakopoulos, K. A. Jenkins, D. B. Farmer,H.-Y. Chiu, A. Grill, and P. Avouris, Science , 662 (2010).[10] B. Sensale-Rodriguez, R. Yan, M. M. Kelly, T. Fang, K. Tahy,W. S. Hwang, D. Jena, L. Liu, and H. G. Xing, Nat. Commun. , 780 (2012).[11] T. Otsuji, S. A. B. Tombet, A. Satou, H. Fukidome,M. Suemitsu, E. Sano, V. Popov, M. Ryzhii, and V. Ryzhii,J. Phys. D: Appl. Phys. , 303001 (2012).[12] F. Schedin, A. K. Geim, S. V. Morozov, E. W. Hill, P. Blake,M. I. Katsnelson, and K. S. Novoselov, Nat. Mater. , 652(2007). [13] M. Pumera, Mater. Today , 308 (2011).[14] Q. Yu, J. Lian, S. Siriponglert, H. Li, Y. P. Chen, and S.-S. Pei,Appl. Phys. Lett. , 113103 (2008).[15] X. Li, W. Cai, J. An, S. Kim, J. Nah, D. Yang, R. Piner, A. Ve-lamakanni, I. Jung, E. Tutuc, S. K. Banerjee, L. Colombo, andR. S. Ruoff, Science , 1312 (2009).[16] J. C. Meyer, A. Geim, M. Katsnelson, K. Novoselov, T. Booth,and S. Roth, Nature , 60 (2007).[17] Y.-C. Lin, C.-C. Lu, C.-H. Yeh, C. Jin, K. Suenaga, and P.-W.Chiu, Nano Lett. , 414 (2012).[18] C. Casiraghi, S. Pisana, K. Novoselov, A. Geim, and A. Ferrari,Appl. Phys. Lett. , 233108 (2007).[19] J. Martin, N. Akerman, G. Ulbricht, T. Lohmann, J. H. Smet,K. von Klitzing, and A. Yacoby, Nat. Phys. , 144 (2008).[20] Y. Zhang, V. W. Brar, C. Girit, A. Zettl, and M. F. Crommie,Nat. Phys. , 722 (2009).[21] E. Rossi and S. Das Sarma, Phys. Rev. Lett. , 166803(2008).[22] S. Adam, E. H. Hwang, V. M. Galitski, and S. Das Sarma, Proc.Natl. Acad. Sci USA , 18392 (2007).[23] A. Deshpande, W. Bao, Z. Zhao, C. N. Lau, and B. J. LeRoy,Phys. Rev. B , 155409 (2011).[24] K. McCreary, K. Pi, A. Swartz, W. Han, W. Bao, C. Lau,F. Guinea, M. Katsnelson, and R. Kawakami, Phys. Rev. B , 115453 (2010).[25] T. Ando, J. Phys. Soc. Jpn. , 074716 (2006). [26] E. H. Hwang, S. Adam, and S. Das Sarma, Phys. Rev. Lett. ,186806 (2007).[27] J. H. Chen, C. Jang, S. Adam, M. S. Fuhrer, E. D. Williams,and M. Ishigami, Nat. Phys. , 377 (2008).[28] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgen-frei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, andJ. Hone, Nature Nanotechnology , 722 (2010).[29] Y.-W. Tan, Y. Zhang, K. Bolotin, Y. Zhao, S. Adam, E. Hwang,S. D. Sarma, H. Stormer, and P. Kim, Phys. Rev. Lett. ,246803 (2007).[30] T. Wehling, S. Yuan, A. Lichtenstein, A. Geim, and M. Kat-snelson, Phys. Rev. Lett. , 056802 (2010).[31] Q. Li, E. H. Hwang, E. Rossi, and S. Das Sarma, Phys. Rev.Lett. , 156601 (2011).[32] T. Radchenko, A. Shylau, and I. Zozoulenko, Phys. Rev. B ,035418 (2012).[33] J. Xue, J. Sanchez-Yamagishi, D. Bulmash, P. Jacquod,A. Deshpande, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero,and B. J. LeRoy, Nat. Mater. , 282 (2011).[34] C. Jang, S. Adam, J.-H. Chen, E. D. Williams, S. Das Sarma,and M. S. Fuhrer, Phys. Rev. Lett. , 146805 (2008).[35] W. Kohn and J. M. Luttinger, Phys. Rev. , 590 (1957).[36] E. Rossi, S. Adam, and S. Das Sarma, Phys. Rev. B , 245423(2009).[37] S. Adam, E. Hwang, E. Rossi, and S. D. Sarma, SolidState Communications , 1072 (2009), recent Progress inGraphene Studies.[38] K. J. Willis, S. C. Hagness, and I. Knezevic, J. Appl. Phys. ,063714 (2011).[39] N. Sule, K. J. Willis, S. C. Hagness, and I. Knezevic, J. Com-put. Electron. (2013), published online, DOI: 10.1007/s10825-013-0508-1.[40] K. J. Willis, S. C. Hagness, and I. Knezevic, Appl. Phys. Lett. , 122113 (2013).[41] C. Jacoboni and L. Reggiani, Rev. Mod. Phys. , 645 (1983).[42] A. Taflove and S. Hagness, Computational Electrodynamics:The Finite-Difference Time-Domain Method (Artech House,2005).[43] K. Yee, IEEE T. Antenn. Propag. , 302 (1966).[44] M. V. Fischetti and S. E. Laux, J. Appl. Phys. , 1205 (2001).[45] R. P. Joshi and D. K. Ferry, Phys. Rev. B , 9734 (1991).[46] D. Rapaport, The Art of Molecular Dynamics Simulation (Cam-bridge University Press, 2004).[47] T. Fang, A. Konar, H. Xing, and D. Jena, Appl. Phys. Lett. ,092109 (2007). [48] K. R. Knox, S. Wang, A. Morgante, D. Cvetko, A. Locatelli,T. O. Mentes, M. A. Ni˜no, P. Kim, and R. M. Osgood, Phys.Rev. B , 201408 (2008).[49] S. Laux, IEEE Trans. Comput-Aided Des. Integr. Circuits Syst. , 1266 (1996).[50] W. H. Press, Numerical recipes : the art of scientific computing (Cambridge University Press, New York, N.Y., 1989).[51] N. Sule and I. Knezevic, J. Appl. Phys. , 053702 (2012).[52] K. M. Borysenko, J. T. Mullen, E. A. Barry, S. Paul, Y. G. Se-menov, J. M. Zavada, M. B. Nardelli, and K. W. Kim, Phys.Rev. B , 121412 (2010).[53] A. Konar, T. Fang, and D. Jena, Phys. Rev. B , 115452(2010).[54] X. Liang, B. A. Sperling, I. Calizo, G. Cheng, C. A. Hacker,Q. Zhang, Y. Obeng, K. Yan, H. Peng, Q. Li, X. Zhu, H. Yuan,A. R. Hight Walker, Z. Liu, L.-m. Peng, and C. A. Richter,ACS Nano , 9144 (2011).[55] A. Priolo, H. M. Jaeger, A. J. Dammers, and S. Radelaar, Phys.Rev. B , 14889 (1992).[56] Z. Fan, Y. Wu, X. Zhao, and Y. Lu, Comp. Mat. Sci. , 301(2004).[57] J.-S. Ferenc and Z. N´eda, Physica A , 518 (2007).[58] A. M. Kriman, M. J. Kann, D. K. Ferry, and R. Joshi, Phys.Rev. Lett. , 1619 (1990).[59] P. Lugli and D. K. Ferry, Phys. Rev. Lett. , 1295 (1986).[60] P. Lugli and D. K. Ferry, IEEE Trans. Electron Devices , 2431(1985).[61] M. Lundstrom, Fundamentals of Carrier Transport (CambridgeUniversity Press, Cambridge, UK, 2000).[62] L. Vicarelli, M. S. Vitiello, D. Coquillat, A. Lombardo, A. C.Ferrari, W. Knap, M. Polini, V. Pellegrini, and A. Tredicucci,Nat. Mater. , 865 (2012).[63] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotke-vich, S. V. Morozov, and A. K. Geim, Proc. Natl. Acad. Sci.U.S.A. , 10451 (2005).[64] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov,Science , 666 (2004).[65] A. F. Young and P. Kim, Nature Phys. , 222 (2009).[66] E. R. Mucciolo and C. H. Lewenkopf, J. Phys.: Condens. Mat-ter22