In fe w -la y e r g r a p h e ne ( FLG) systems on a die le ct ric substrate such as SiO 2 , t he addit io n of each extra layer ofgraphene can drastica lly a lt e r t he ir e le ct ro n ic and structural properties . He re , we map t he chargedistributio n a mo ng the ind ivid ua l la yers of fin ite -size FLG systems using a no ve l spatial discrete model t hat describes bot helectrost at ic int e rla ye r s cre e ning and fringefield effects. Our results reve a l that the charge density in the re gio n ve ry closeto the edges is screened out anorderof magnit ud e mo re we a kly than that acrossthe central region of the laye rs. Our discrete model suggests that the interla ye r charge screening le ngt h in 1-8 la ye r t hick graphene syst ems depends most ly on the ov era ll gat e/ mo le cu la r doping level rat he r t ha n on temperature, in pa rt icula r at an induced charge density > 5 X 1 0 12 cm -2 and can relia b ly be determined t o be larger t ha n half the int e rla ye r spacing but shorter than t he bila ye r t h ickne ss . Our model can be used for designing FLG-based devices, and offersa simple rule regarding the charge distribution in FLG: a pprox imat e ly 70%, 20%, 6% and 3%(99% overall) of the total induce d cha rge density reside within the four in ne rm ost laye rs, im plying t hat the gate-induced ele ct ric fie ld is not defi nit e ly felt by > 4th laye r.
Since its discovery in 2004, single-layer graphene (SLG) has become the most studied nanomaterial due to its exceptional mechanical 1 electrical2 an d opticaP properties. Although several physical properties are shared between SLG and few-layer grap hene (FLG), increasing layer thickness can give rise to a unique range of electronic and structural properties that has not yet been sufficiently understood, in particular for FLG systems with more than 3 layers. More specifically,electrical noise, charge transport and nonlinear optical properties of FLG on substrates (usuallySiOifSi) exhibit strongdependence on the number oflayers, gate-induced chargedensities and underlying oxide substrates. It is therefore crucial in the design of FLG-based high-speed transistors 4 , terahertz plasmonics 5 , photonicsand optoelectronicdevices6 to quant itatively unde rstand the role of the number of layers in the charge distribution and the electric field screening of the FLG/SiOifSi systems and also to explore the unclear relationship between the excess gate-ind uced chargedensities and the layer-by-layer Fermi level and charge density profiles in the FLG systems.
Owing to the importance of the subject, the question of interlayer charge screeni ng length >.. in th e FLG systems has been addressed by several exper imental methods, including angle-resolved photoemission spectroscopy (>..= 0. 1 4 -0 . 19 nm )7, no ndegenerate ultrafast mid-infrared pump-probe spectroscopy (>..= 0.34nm) 8 , Kelvin probe force microscopy9• 13 ( >..= 1.36 -l. 7 0 nm 11 , >..= 0 . 42 nm 1 2 an d >..= 2.4 nm 13 ) , single -gated fie ld effect trans istor (>..= 0 . 6 n m) 1 4 , d oub le -ga ted fie ld e ffect transistor (>..= 1.2 nm ) 15 and dark -fie l d sc atte r ing spec t roscop y (>..= 1 .2 ± 0 .2 n m ) 16 However, a relatively wide range of experimental values for>.. (from less than a single layer to seven layers) is observed, which is not yet fully understood . Nevertheless, we believe that a part of this data scattering maybe attributed to the dependenceof the screening length on the device quality and experimental conditions, such as sample preparation processes, the presence of defects and impurities in graphene, the intrinsic charge density in each graphene layer and the actual doping level of the system. This diversity in the reported values of>.. is also seen in theoretical approaches 1 • 7 2 • 0 Depending on whether the inter-layer electron tunneling is takeninto account or not, >..between 0.54 nm 17 and 0.7 nm 18 is obtained using a random phase approximation.
Kuroda and coworkers theoretically reported that both the gate charge and temperature could highly influence >.. , whose value may range from ~0.2nm to 3.1nm1 9 . We will later show in this paper that the presence of the effective m ass, a key missing param eter in Kuroda's model' 9 not onlyleadsto a much narrower range of>..values (= 0.2 -0.7 nm), but also rules out the possible effect of temperature on the reported values of>... We also note De partme nt of Me cha nical Eng inee ring, Unive rsit y of Michigan, Ann Arbor, Mich igan 48109, United Sta t es.
Correspondenceand requestsfor materials shou ld be add ressed to W.L.(em a il: weilu@ umich.edu) Finite-size FLG flakes and graphene nanoribbons in actual devices exhibit an intriguing dependence of the electrostatic and electrical conductivity response on their geometrical parameters (e.g., lateral sizes, thicknesses, shapes and edge types) 21 .2 2 _ Both experimental and theoretical studies have demonstrated that a strong charge accumulation takes place at the edges of the finite-size graphene flake due to the electrostatic fringe field effects 2 >-3 • 1 Scanning gate microscope measurements of a monolayer graphenedeviceon a SiOifSi substrate reveal significant conductance enhancement at the edge of the graphene devicedue to the strong chargeaccumulation 23 • Similar observations of inhomogeneous chargedensityand capacitance profilesnear the edgesof bothsuspended and hBN-supported mono/bilayergraphene devices have been reported usingquantum Hall edge channels 24 .2 5 Among different theoretical models on the charge distribution of the finite-sized graphene, we particularly note a strong charge accumulation at the edges and the corners of a positively charged rectangular graphene sheet using the charge/dipole molecular dynamics modeJ2 6 .2 7 and a long the edges of a graphene nanoribbon using the tight-binding modeJ2 8 Despite recent progress, a detailed understanding of the electrostatic charge distribution in connection with the actual electronic structure of FLG is still lacking. Here, we exploit the layered nature of FLG to develop a novel spatial discrete model that successfulyl accounts for both electrostatic screening and fringe field effects on the charge distribution of the finite-size FLG system . To this end, an effective bilayer model based on two tight-binding parameters is utilized to accurately describe electronic band structures and thus density of states (DOS) of one to eight Bern al-stacked graphene layers. We then explore the unclear relationship between the gate-induced charge densities and layer-by-layer Fermi level and charge density profiles in FLG systems usin g a global energy minimization, where its total energy is calculated based on electrostatic interaction between graphene layers and band-filling energy in each layer. Our discrete model offers a unique capability to quantify the nonlinear charge density profile, interlayer capacitance, quantum capacitance, and local surface electrostatic potential of FLG by showing a verygoodqualitative and quantitative agreement between the previously measured work functions in FLG and our theoretical results.
We first examine the charge distribution of an FLG/SiOifSi system containin g N (up to 8) layers of finite-size graphene sheet with desired shapes (i.e., square, rectangle, circle or ribbon), as schematically illustrated in Fig. 1(a). Each graphene layer is labeled by an integer number starting from i = 1 for the layer closet to the substrate (hereafter referred to as the innermost layer) to i = N for the top layer (as the outermost layer). Applying a biasvoltage V 0 between the highly-doped Si substrate and N-layer graphene (N-LG) inducesa total excesscharge density of Q 0 in N-LG, whose layer i can carry a charge density of Q 1 such that the following constraint holds Qo = I: f: 1 Q-;
The electronic bands of N-LGcan be modeled by two tight-binding parameters, namely, the nearest neighbor hopping parameter 10 (which defines the Fermi velocity vf = (3l 2 )'Y 0 a /h , wh ere a= 0.142 nm is the C-C bond length) and the nearest neighbor interlayer coupling constant 1 1 . We take 10 = 3.14 eV an d 11 = 0.4eV as typical values of bulk graphite. The energy dispersion in Bernal-stacked N-LG, obtained from 2D cuts in the electronic dispersionof graphite, perpendicular to the graphene planes at specific values of(}=j7rl2 (N + 1) , can be given by kJ = E 2 f,y 2 ± 2mj E/h 2 , where, = v_l,(Ii being the reduced Planck constant ), mj = {, / v}) sin(} is the effective mass, j (= 1, 3, 5, ... , N -1 for even layers and 0, 2, 4, .. ., N-1 forodd layers) istheindex offue energy band with which are produced by the band extrema at the K-point, followed by a linear increase with kinetic energy lE. Of particular importance for the electronic structures of N-LGat low energies is the excitation energy from the ground state (Diracpoint) to the first excited state (denoted byJE ), as explicitly shown in Fig. 20).
We next determine the charge distribution profile in a finite-size N-LG stack with a circular shape of radius R, based on the method of images, followed bysolving the Love equation (Section Sl.1ofSupplemental Material (SM)). Thechargedensity profile in the circular layer i can then be expressed by (2) where
is the charge distribution profile, normalized to its average value (j} for generality purposes; the index notation i varies from l to N; If" (= r !R) is a dimension less parame ter; r de notes the radial coordinate of atom and g (lf") is a polynomial function ofR which onlydepends on the ratio of the graphene size to the dielectric thickness (Fig. Sl of SM). A new parameter a 1 (>0) isintroduced by Eq. ( 3) in order to determine the amount of charge densityat the edgeof the layer i (11" = 1 ) . Although the focusof the present work is on graphene flakes with a circular shape, we note that the charge distribution of circular graphene flakes and graphene nanoribbons is ofa similar form as given by Eq. ( 3) and, therefore, doesnot qualitatively and pretty much quantitativelyalter the main resultsof this paper(Section Sl.2ofSM). We also refer the interested reader to Section Sl.3ofSMfor the corresponding charge distribution profile of rectangular/square graphene flakes.
As we already discussed, in practice, the charge distribution in electrostaticallydoped graphene devices is inhomogeneous, yielding a non-uniform Fermi level profile. For instance, scanning gate microscope measurements ofa monolayer graphene deviceon a SiOz/Si substrate reveala strongshift of the local Dirac point from the Fermilevel at the grapheneedge due to the contribution of both localized edgestates (i.e., zigzag or armchair) and accumulated charge along the edge 2 • 3 The Fermi energy profile eF; across the layer i can be expressed in terms of theconstant Fermi energy eF, as follows (Section S2 of SM)
where(lEFi) is the average value oflEFi in terms oflE F ; and a •1 The avera ge charge density Q,can be obtained by minimizing the total energy of the system with respect to eF; and a 1 as the varia tional parameters under the constraint thatQ 0 = I;Z,. 1 Q,. In the N-LG/SiO;JSi system, the total energy can be split as, U, = U,+ U, + Ub,where -----------------.
C o m p ari son Stud ies. In order to verify the accuracy of the results presented in this paper, wefirst compare our local work functions (cl>, = -eV;) with those measured by angle-resolved photoemission spectroscopy (ARPS)7and Kelvin probeforce microscopy (KPFM) 1 We note that sincethe accurate work function of the tip under the ambient conditions and also the accurate valueof the dielectric constant for the N-LG/SiO2 interface are unknown, the difference of the work function is used to achieve more accurate comparison purposes. We begin bycomparing <1> 1 in a 4 -LG system with that measured byARPS 7 , as shown in Fig. 3(a).The resultsare given relative to the work function of the outermost layer <1> 4 as the zero -reference level and Q0 is set to be 2.2 x 10 13 cm -2 . It is evident from Fig. 3(a) that a very good agreement exists between the proposed discrete model and those measured byOhta et al. 7.Another comparison study is conductedin Fig. 3(b) between the present discrete model and KPFM resultsof Ziegler et a /. 11 , who measured cl>, in t he 1-6 -LG systems relative to that of bulkgraphitecl> 00 . F i gure 3(b) clearly demonstrates that the measured work functions are generally in much better agreement with our results than those obtained by ab initio OFT calculations 11 when assuming a total induced charge density of 4.85 x 10 12 cm -2 We further perform a similar comparison in Fig. 3(c) between the present work functionsat the uppermost layer of N-LG (clN > ) relative to those of (N-1)-LG (clN > _ 1 ) with KPFM results measured for N-LGwith layer number rangingfrom 1 to 8 12 . It is indicated that the present work functionsclosely match with the experimental observations for Q0 = 1.7 x 10 13 cm -2 Further comparison study is performed in Fig. 4 to investigate the influen ce of the effective mass m• on the charge distribution of an 8-LG system. It is seen from Fig. 4 that the model based on the monolayer-uke band structurefailsto accurately predict the chargedistribution of the 8-LG system,in particular at the smaller induced charge densities. This figure also shows a significant deviation in the charge densities of layers i > 5 for Q0=10 13 cm - 2 Also, our energy evaluations of N-LG systems under a given Q 0 for three possible charge distribution scenarios-(a) optimum distribution given in Eq. ( 3), (b) non-uniformdistribution with the chargesingularity at the very edge (i.e., a 1 =0), and (c) fully uniform distribution (i.e.,q,= Q,)reveal that the minimum energyis only achieved by the present optimum charge distribution model, further indicating its merit in predicting the charge distribution of other families of atomically thin layeredmaterials. 7 and Wang et al. 12 , the charge density is drastically reduced as one move away from the innermost toward the outermost layer. However, the charge dens ity in the region very close to the edges is screened out an order of magnitude more weakly than that across the central region of the layer, as shown in Fig. S(b), which can be explained by the presence of the strong fringe field along the edges, as schematically shown in Fig. 1. Our results in Fig. S(a) also suggest that the innermost layer plays the most important role in the electrostatic charge distribution of the N-LG systems by hosting~70% of the gate charge density Q 0 • Hence, it is worth looking into its Fermi level profile more in detail, as illustrated in Fig. S(c). By following the evolution of the Fermi level along the innermost layer, it is observed that a strongcharge accumulation and thus sufficiently large shift in the Fermienergy at the edge can give rise to a jump in the electronic band structures of 5-LG toward the first excited state, 0.4eV (as shown in green solid curve in Fig. S(c) and in green dashed curve in the inset, which shows the energy band structure of the 5-LGsystem). However, our Fermi level analyses in the inn ermost layer of 6-and 8-LG systems exhibit few jumps in the Fermi level of the regions both close to and away from the edges when Q 0 =10 13 cm -2 (see Fig. S2(b) of SM for detailed discussions). This can be attr ibuted to the fact that the lowest energy of the first excitation band decreases for the N-LG system with a larger number of graphene layers, as shown in Fig. S(b).
To quantitativelyelucidate the correlation between the magnitude of the gate charge density Q 0 and the averagecharge distribution Q; through the 5-LG thickness, Fig. S(d) shows Q;IQ 0 ratio as a function of the layer positions for three different values of Q 0 (= 10 12 , 10 13 and 10 14 cm -2 ) . It is seen th at a larger value of Q 0 leads to a strongercharge screening normal to the layers,however, this effect dimin ishes when Q0 <10 12 cm-2 This figure also demonstrates that almost 90% of the excess charge density resides in the first two layers, implying that the interlayerscreening length can reliably be determined to be less than ~0.7nm. Having Q; data for each layer enables us to calculate the "local" (interlayer) charge screening >.. ;,+ ; 1 as Q; + 1 fQ ; = exp(-df\ ,;+i ) based on Thomas-Fermi charge screening theory (seeSection S4 of SM for the calculation of the interlayer screening). It is deduced from Fig. 2(d) that the charge screening length between the first and second layers >..1.2 may redu ce from ~Id at Q 0 = 10 12 cm -2 to ~0.Sd at Q 0 = 10 14 cm -2 , while a smaller variation in >..; ;+ 1 is observed for the layers farther from the substrate due to the reduction in their DOS at the Fermilevel.
Layer-Dependent Charge Screen ing in N-LG Syst ems. We now turn to a discussion of the laye-rdepen dentcharge distribution/chargescreening in 1-8-LGsystems for a given gate-indu ced charge density of 10 13 cm -2 . Figure 6(a) presents a plot of Q;IQ0 versusthe layer posi tions in 1-8-LG systems, indica ting that approximately 70%, 20%, 6% and 3% (99% overall ) of Q 0 sit in layers i = 1 to 4, respectively, and thus the gate-inducedelectric field is not definitely felt by i > 4 layers. Interestingly, weobserved that thechargedensity of the layers located in the same position in N-LG systems decreases in a sawtooth-like fashion, as shown in the insets of Fig. 6(a) for the normalized charge density of the innermost Q/ Q 0 and second innermost Qif Q 0 layers. This saw-tooth pattern which is associated with the presence of the lin ear energy dispersion in N-LG with odd layer numbe r has been experimentallyconfirmed through the measurement of the electric double-layer capacitance between an ionic liquid and l -6-LG 33 . The results in Fig. 6(a) provide an important piece of information about the charge screening effect of the innermos t layer on different layers of 2-8-LG. Hence, we first define a "global" (effective) charge screening >-. as Q ;IQ 1 = exp[-d(i-1)/>-.] . This new definition of the "global" charge screening length allows us to explore how the innermost layer impacts the surface potential drop across the FLG thickness and also provides a single value of the screening length to predict the charge distribution of all layers relative to that of the innermost layer. Keeping both global and local screening definitions in mind, we observe from Fig. 6(b) that our global charge screening can be well fitted by the simple exponential decay function (in particular for Q 0 ::; l0 13 cm -2 , see Fig. S3ofSM) when>-. i::, d. We next address the problem of the charge accumulation along the graphene edge, focusing first on very limited publications that have quantitatively studied the charge density at the edge of graphene thus far. From prior experimental work, a nearly th ree-fold increase in capacitance and thus the charge density near the edge of a suspended bilayer flake (0.4µm wide and 2.6µm long) was observed using quantum Hall edge channels 2 4.
Fro m t heoretical pointsofview, the charge/dipole molecular dynamics model predicts a seven-fold (fifteen-fold) enhancement of the charge density at the edge (corner) over that at the center ofa charged 8.5 nm x 4.8 nm rectangular graphene sheet 26 and a similar eight-fold enhancement of the chargedensity in a 20-nm-wide graphene nanoribbon 27 . This model also suggests that the charge enhancement is more significant in multi-layered graphene in such a way that thechargedensityat the edge relativeto that at thecenter canvary from 9 in the inner layer to > 14 in the outer layer of a 4-LG nanoribbon system 2 • 7 Also, using the tight-binding Hartree model, the charge density along the edge of a 20-nm-wide graphene nanoribbon enhances up to five times over that at the center 28 . Having this quantitative description of the charge accumulation at the graphene edge in mind, we present in Fig. 6(c) the chargedensity at the edge relative to that at thecenter,q.•I q.C,as a fu n ction of the layer position in the 1-8-LG systems for Q 0 = l0 13 cm-• 2 As is evident from the figure, o discrete model predicts the edge-to-center charge density ratio for monolayer graphene to be ~7.5 which is consistent with the theoretical results 2 6-28 . Surprisingly,the addition ofeach extra layer reduces the charge accumulation at the edge of the innermost layer from 7.5 in 1-LGdown to ~5 in 8-LG, whereas an inverse trend is observed for the charge accumulation at the • 5 -------------------------------N-=-1--N ----8-0.
Qo (cm -2) Temperature-Dependent Charge Screening Model. While the present study has focused on the charge distribution of N-LG at absolute zero temperature, wenote that a variation in temperature from zero to room temperature has no appreciable effect on the charge screening length, more specifically at the higher gate electric field. Following a tempera ture-dependent model of the charge distribution detailed in Section S5 of SM, the local chargescreening between the first and second layers of an 8-LG system is plotted in Fig. 7 as a function of Q0 at T = 0 and 300K. For comparison purposes, the results of Kuroda et al. 19 based on the linear energy dispersion are reproduced bysettingmj = 0,as indicated bydashed curves with opensymbols in Fig. 7. It isevident from Fig. 7 that the interlayerchargescreening is insensitive to the temperature variation when Q 0 ?_ 5 x 10 1 2 cm -2 and onlya slight changein >.,.2 is observed at smaller gate charge densities (see lower inset) and ultimately saturates to >.1,2 d. Consistent with our temperature-independent charge screening length, Yang and Liu reported using the first-principles calculations that the interlayer screening, static perpendicular dielectric function and density of states of bi-and tri-layer graphene slightly changes as temperature increases from OK to 300 K to 600K 34 It is alsoobserved from Fig. 7 that thelinear dispersion model fails to predict the interlayer chargescreening between the two innermost layers for Q 0 :;;: 10 12 c-m 2 such that >. 1 ,2 goes to infinity (i.e. Q 2 Q 1 ) at T = 0 as Q 0 -+ 0. Interestingly,a layer-by-layer inspection of the charge density in a similar 8-LGsystem for different values of Q 0 reveals that the linear dispersion model not only yields inconsistent charge density profiles in almost all 12 2 layers for Q 0 ::; 10 cm -but also shows a significant deviation in the charge densities of outer layers for Q 0 >10 12 cm -2 , as shown earlier in Fig. 4. Thisdeviation from our model can be understood in terms of the effective mass in N-LG with N?_ 2: an essential ingredient that is not captured in Kuroda's model where an N-LGsystem is considered as N parallel single layers with a massless linear energy dispersion (upper inset for N = l ), rather than a single 2D system with the actualenergy dispersion (upper inset for N = 8).
We developed a novelspatial discrete model to unravel the relationship between the macroscopic induced charge density and microscopic (layer-by-layer) charge distribution in finite-size FLG through considering the effects of both electrostatic interlayer screening and fringe field. We showed that adding each extra layer reduces the charge accumulation at the edge relative to that at the center of the innermost layer up to 20% (from ~7.5 in 1-LG down to ~5 in 8-LG). Our model offers a sim ple rule of thumb regarding the charge distribution in FLG: approximately 70%, 20%, 6% and 3% (99% overall) of the total induced charge density reside within the four innermost layers (layers i= l to 4, respectively),implying that the gate-induced electric field is not definitely felt by layers i> 4. We finally found that a variation in temperature from zero to 300K has no appreciable effect on the interlayer charge screening when the gatecharge density is larger than ~5 x 10 12 cm -2 • Although our study is concerned with FLG systems, the generalityof our spatial discrete model suggests that the chargedensity profile, interlayer screening, quantum capacitance, and local surface potential of other atomically thin layered materials (ATLMs),such as semiconductingtransition metal dichalcogenides(e.g., MoS 2 , WSe 2 and WSi) and heterostructures (e.g., graphene/MoS 2 and MoSifWSe 2 ) , can be characterized by feeding relevant electronic band structures of ATLMs into our model. In addition, the effect of structura l defects (e.g., vacancies, adatoms, dislocations and grain boundaries) and stacking faults on the chargedistribution ofdefective FLG systems can be studied by modifying DOS of pristine FLG.
© The Author(s) 2017