EMC/FDTD/MD simulation of carrier transport and electrodynamics in two-dimensional electron systems
aa r X i v : . [ c ond - m a t . m e s - h a ll ] M a y EMC/FDTD/MD simulation of carrier transport and electrodynamics in two-dimensional electronsystems
N. Sule, ∗ K. J. Willis, S. C. Hagness, and I. Knezevic † Department of Electrical and Computer Engineering, University of Wisconsin-Madison, Madison, WI 53706, USA
We present the implementation and application of a multiphysics simulation technique to carrier dynamicsunder electromagnetic excitation in supported two-dimensional electronic systems. The technique combines en-semble Monte Carlo (EMC) for carrier transport with finite-difference time-domain (FDTD) for electrodynamicsand molecular dynamics (MD) for short-range Coulomb interactions among particles. We demonstrate the useof this EMC/FDTD/MD technique by calculating the room-temperature dc and ac conductivity of graphenesupported on SiO . I. INTRODUCTION
Electronic properties of supported two-dimensional (2D)materials, such as the single layer graphene [1] or
MoS [2],and quasi-2D materials, such as semiconductor membranes[3], have generated a lot of interest in recent years. These2D electronic systems (2DESs) have potential applications aselectronic [4, 5] and optoelectronic devices [6, 7], THz detec-tors [8], as well as chemical and biologicial sensors [9]. Re-alizing these applications requires an understanding of carriertransport in 2D materials in the presence of electromagneticfields, while accounting for the strong influence of the sup-porting substrate [10, 11] and the impurities found near the2DES/substrate interface [12, 13].A multiphysics numerical solver that combines ensembleMonte Carlo (EMC) simulation of carrier transport with thefinite-difference time-domain (FDTD) technique for solvingMaxwell’s curl equations and molecular dynamics (MD) forshort-range particle interactions [14] can provide insight intothe carrier transport and electrodynamics of 2DESs. Un-like most device simulation tools that implement EMC cou-pled with a quasi-electrostatic solver of Poisson’s equation[15, 16], EMC/FDTD/MD couples EMC with a fully elec-trodynamic solver [17–20], which enables simulation of car-rier dynamics under electromagnetic excitation, from low fre-quencies (including dc ) to the THz frequency range. (At fre-quencies above THz, interband transitions in semiconductorsbecome important and the classical view of carrier–field inter-action is no longer sufficient.)Grid-based solvers, such as Poisson solvers or FDTD, can-not accurately capture the forces among charges on spa-tial scales smaller than the size of a grid cell. In or-der to include short-range (sub-grid cell) interactions, EMC-FDTD has been extended through coupling with MD [21–24].The EMC/FDTD/MD technique includes accurate pair-wise,short-range, real-space Coulomb forces among carriers andbetween carriers and charged impurities, together with thefull electrodynamics solution for long-range Coulomb fields[14]. EMC/FDTD/MD has been used to accurately calculatethe conductivity of bulk Si in the THz frequency range, where ∗ [email protected] † [email protected] FIG. 1. Schematic of the simulated structure: single-layer graphenerests on an
SiO substrate, with charged impurities present near theinterface between the two. Carrier transport is simulated by includingboth electrons and holes in the graphene layer, while the positivelycharged ions near the interface and within the SiO substrate remainstationary. the usual Drude model fails, with good agreement to experi-mental data [14, 25].In this paper, we describe the EMC/FDTD/MD techniqueas applied to simulating carrier transport in 2DESs. We sim-ulate a structure with a single graphene layer resting on topof an SiO substrate, with impurity ions near the interface(Fig. 1). We describe the constituent techniques (Sec. II)and the procedure for coupling them (Sec. III). As the dis-tribution of impurity ions is important for the overall carriertransport properties [26], we describe the generation of a de-sired impurity distribution, from uniformly random to clus-tered (spatially correlated). In Sec. IV, we give examples of dc and ac conductivity of supported graphene in the presenceof charged impurities, as calculated using EMC/FDTD/MD.We conclude with Sec. V. II. CONSTITUENT TECHNIQUES
This section provides a brief overview of the constituenttechniques of the combined EMC/FDTD/MD solver, with fo-cus on the implementation details relevant for the simulationof carrier transport in graphene.
A. Ensemble Monte Carlo (EMC)
Ensemble Monte Carlo simulates carrier dynamics in thediffusive transport regime [16]. This method yields a solutionto the Boltzmann transport equation by using statistically ap-propriate stochastic sampling of the relevant relaxation mech-anisms, free-flight times, and angular distributions of mo-menta [27]. In an EMC simulation, the evolution of a largeensemble of carriers [typically O (10 ) ] is tracked over time.The evolution of physical properties of interest, such as thecarrier average drift velocity or kinetic energy, are calculatedby averaging over the ensemble.During the simulation, each carrier undergoes periods of”free flight”, or drift, under the influence of the Lorentz force, ~F = q ( ~E + ~v × ~B ) , (1)interrupted by instantaneous scattering events. In Eq. (1), ~E and ~B are the electric and magnetic field vectors, respec-tively, q is the carrier charge, and ~v is the carrier velocity.In our EMC simulation, we include both electrons and holesin graphene. The carrier velocity in graphene is given by ~v = v F ~k | ~k | , where v F is the Fermi velocity and ~k is the car-rier momentum. For a free flight of duration t d (obtainedstochastically [27]), the momentum and energy of a carrierare updated based on the Lorentz force and E - k dispersion, as ~r new = ~r old + Z t d ~v ( t ) dt , (2a) ~k new = ~k old + ¯ h − Z t d ~F ( ~r ( t )) dt , (2b) E = ¯ hv F | ~k new | . (2c)The electron-phonon scattering rates in graphene are calcu-lated based on the third-nearest-neighbors tight-binding Blochwave functions (3NN TBBW) [28] and are shown in Fig. 2.(Near the Dirac point, electron and hole dispersions, as wellas their rates for scattering with phonons, are considered to beidentical.) The deformation potential constants ( D ac = 12 eVand D op = 5 × eV/m) were determined based on fittingthe longitudinal acoustic (LA) and optical (LO) phonon scat-tering rates to the rates calculated from first principles [29].The surface optical (SO1 and SO2) phonon scattering ratesare calculated from the interaction Hamiltonian following thedielectric continuum model of surface phonons [30]. B. Finite-Difference Time-Domain (FDTD) method
In FDTD [31], the time-dependent Maxwell’s curl equa-tions µ ∂ ~H∂t = −∇ × ~E − ~M (3a) ǫ ∂ ~E∂t = ∇ × ~H − ~J (3b) LALOSO1SO210 sc a tt e r i ng r a t e [ s − ] FIG. 2. Electron–phonon scattering rates in graphene, calculated us-ing the third-nearest-neighbor tight-binding Bloch wave functions(3NN TBBW) [28]. Scattering due to longitudinal acoustic (LA)and optical (LO) phonons intrinsic to graphene, as well as the sur-face optical phonons (SO1 and SO2) between
SiO and graphene, isincluded. are discretized using a centered-difference scheme for the par-tial derivatives in both space and time [32]. The field compo-nents ( ~E and ~H = µ − ~B ) are spatially staggered. Equations(3) are solved by leapfrog time integration: the ~E -field and ~H -field updates are offset by half a time step, yielding a fullyexplicit scheme with second-order accuracy in time.Fig. 3 shows a schematic of two FDTD grid cells above andbelow the plane of graphene. The field components ( ~E and ~H = µ − ~B ) are spatially staggered (Fig. 3). The E x and E y field components in the ( k +1) -th plane, as well as the E z fieldcomponents in the ( k + 1 / -th plane, are updated assumingmaterial properties corresponding to air ( ǫ a = 1) ; the E x and E y field components in the ( k ) -th plane are updated assuminggraphene properties ( ǫ g = 2 . ; and the E x and E y fieldcomponents in the ( k − -th plane, as well as the E z fieldcomponents in the ( k − / -th plane, are updated assuming SiO properties ( ǫ s = 3 . .The convolutional perfectly matched layer (CPML) absorb-ing boundary condition [33], with a thickness of – gridcells, is applied at the top and bottom horizontal boundariesof the domain. We use periodic boundary conditions for thefour vertical boundary planes perpendicular to the graphenesheet. An incident plane wave is introduced via the total-fieldscattered field (TFSF) framework [31]: electric and magneticcurrents ( ~J and ~M ) are calculated using the surface equiva-lence and applied at the boundary between the total-field andscattered-field regions in order to source a propagating planewave. For dc excitation, the electric field component is forcedto remain constant once the peak magnitude of the plane waveis attained. (i,j,k) (i+1,j,k) (i+1,j+1,k)(i,j,k+1) (i,j+1,k+1)(i,j,k-1) (i+1,j,k-1) (i+1,j+1,k-1)(i+1,j,k+1) (i+1,j+1,k+1) ρ AirSiO2GrapheneEx, Jx Ey, JyEz, Jz Hy, My Hx, MxHz, Mz ε xg ε xs ε zs ε yg ε ys ε za ε xa ε ya FIG. 3. Schematic of two FDTD grid cells at the air/graphene/
SiO interface. The top cell is assumed to be filled with air and the bottomcell is assumed to be filled with SiO . The central plane betweenthese two cells represents the graphene layer. The FDTD field andcurrent vectors ( ~E , ~J , ~H , ~M ), staggered in space, are shown witharrows on the grid edges and grid faces. The charge density ( ρ ) isdefined on the grid cell corners shown with open circles. C. Molecular Dynamics (MD)
Molecular dynamics simulates short-range interactions inclassical many-particle systems [34]. For a collection of elec-trons, holes, and charged ions, the particle-particle short-rangeinteractions we include are the direct and exchange Coulombforces among carriers (electrons and holes), and the directCoulomb forces between carriers and ions [14]. We only cal-culate the pairwise interactions among the particles presentwithin a × × -cell volume of one another, in order to min-imize the computational burden in MD, which scales as N , N being the number of interacting particles [35].The carriers in MD are described by Gaussian wave packets[36, 37] with a finite size r c , corresponding to the effectiveradius of the Hartree-Fock exchange-correlation hole [14, 25,38] φ ˜k i ( ~r i ) = 1(2 πr c ) / exp (cid:20) − ~r r + ~k i · ~r i (cid:21) (4)where ~k i and ~r i are the wave vector and position of the i -thcarrier, respectively. The charged impurity ions are also de-scribed by a Gaussian profile with a characteristic radius of r d = 3 . ˚A. Considering that the type and charge of impu-rity ions appear to be strongly dependent on processing [39],we considered positive ions with a unit charge. We swept dc conductivity as a function of r d and picked an r d value froma range over which the dc conductivity does not significantly vary with r d . These Gaussian profiles for charges in MD avoidlarge unphysical forces between pointlike particles that canlead to instability and errors [38]. The Coulomb forces be-tween particles with such Gaussian profiles are given below[14]: ~F D , dij = − q i Q j πǫ g ∇ ~r i (cid:20) r ij erf (cid:18) r ij r d (cid:19)(cid:21) , (5a) ~F D , cij = − q i q j πǫ g ∇ ~r i (cid:20) r ij erf (cid:18) r ij r c (cid:19)(cid:21) , (5b) ~F Exij = − q π / ǫ g r ~r ij k ij (5c) × exp − r r − k r ! Z k ij r c d t e t ,~G ij = − q π / ǫ g r ¯ h (5d) × ∇ ~k ij " k ij exp − r r − k r ! Z k ij r c d t e t . In the above equations, ~F D , dij is the direct Coulomb force be-tween the i -th carrier and the j -th ion, while ~F D , cij is the di-rect force between the i -th and j -th carriers. ~F Exij is the “ex-change force” between the i -th and j -th carriers having thesame charge and spin, and ~G ij is a small correction to the ve-locity of the i -th carrier due to the j -th carrier, stemming fromthe exchange interaction. Also, ~r ij = ~r i − ~r j , ~k ij = ~k i − ~k j ,erf( x ) denotes the error function, ǫ g is the relative permittiv-ity of graphene, while q and Q are the carrier and impuritycharge, respectively. These forces, given in Eq. (5), are cal-culated numerically at the beginning of the simulation for asmall fixed volume in the real and momentum spaces, andstored in lookup tables. III. COUPLED EMC/FDTD/MD
At the beginning of the coupled simulation, the carrierensemble is initialized based on the equilibrium Maxwell-Boltzmann distribution. Poisson’s equation is solved to calcu-late the initial microscopic field distribution stemming fromall the charges in the domain (electrons, holes, and chargedimpurities). As time-stepping commences, carriers in theEMC module drift under the action of the fields and scatter ac-cording to the appropriate scattering mechanisms. Carrier mo-tion results in a current density that can be calculated from car-rier velocities. The positions of the carriers also change, lead-ing to different short-range electrostatic interactions. Thus,the current densities and the new carrier positions can now beused to adjust the fields acting on the carriers (in the FDTDand MD modules). We ensure that the fields from FDTD andMD are not double counted in the vicinity of the charges [14],by subtracting the grid-based (FDTD) contribution of fields inthe vicinity of the charges from the total pair-wise MD con-tribution of the fields. Moreover, to correctly represent fieldsarising from non-uniform and time-varying charge densities inFDTD, the initial field distribution must satisfy Gauss’s lawand the continuity equation must be enforced at each time-step [14]. Thus, accurate and stable coupling of the EMC,FDTD, and MD methods requires proper initialization and as-signment of charges to the grid, initialization of the fields in-cluding the grid-based fields of the impurity ion distribution,calculation of the current density, and calculation of the MDfields for updated positions. In the following subsections, wedescribe the four requirements for coupling in further detail.
A. Charge initialization and assignment
1. Initialization
We assume that the Fermi level and charge density ingraphene can be modulated by a back gate, located at thebottom of the SiO substrate. The simulation domain is notcharge neutral due to the assumption of a back gate, whichis unlike the previous applications of the EMC/FDTD/MDmethod [14, 25]. For a given Fermi level and temperature, thedensity of carriers (electrons and holes) in graphene is givenby [40] n = π (cid:18) kT ¯ hv F (cid:19) R d u u [1 + exp ( u ∓ η )] − R d u u [1 + exp( u )] − , (6)where u = Ek B T and η = E F k B T . Here, the minus (plus) signcorresponds to electron (hole) density. The size of the simula-tion domain is chosen such that the total number of carriers is O (10 ) ; molecular dynamics calculation necessitates that onenumerical particle correspond to one physical particle [14].The momentum and energy of the carriers are defined accord-ing to the Maxwell-Boltzmann distribution. The carriers areinitially positioned according to a uniform random distribu-tion throughout the graphene plane.The impurity ions are distributed below the graphene plane,down to a depth of
10 nm . In our tests, we have ob-served that charged impurities, for reasonable sheet densities( < cm − ), do not appreciably affect transport in thegraphene layer beyond a depth of about
10 nm . The typeand charge of the relevant impurities vary with the process-ing details [39]; for simplicity, here we use positive impurityions with unit charge. In the literature, density of impuri-ties is typically described via a cumulative sheet density, inunits of cm − , however, these ions are distributed throughoutthe three-dimensional substrate. For a generated 3D distribu-tion of ions, the sheet density is obtained by integrating overa depth equal to r d (i.e. twice the typical ion radius) and av-eraging over the
10 nm depth. The positions of impurity ionscan be generated based on a variety of volumetric distribu-tions, from a uniform random to more clustered ones, basedon a correlation length parameter λ . For a non-zero λ , thenumber of impurity clusters, N c , is calculated by N c = A/λ ,where A is the two-dimensional area of the graphene layer.The size (or diameter) of each individual cluster is pickedfrom a uniform random distribution between λ/ and λ/ , so the average cluster size is λ/ . To initialize the impuri-ties, we first distribute the position of the centers of the N c clusters stochastically and then distribute individual impurityions around these centers with a Gaussian distribution, whosemean is the cluster center and the standard deviation equalshalf of the individual cluster size. This procedure results inan overall distribution that has a spatial autocorrelation func-tion (SACF) very close a Gaussian, exp( − r /λ ) , as shownin Fig. 4. λ extracted from the Gaussian fit [Fig. 4(b)] andthe impurity cluster size estimated directly from the full widthat half maximum (FWHM) of the SACF of the impurity dis-tribution are in good agreement. For λ = 0 , we distribute allthe impurity ions stochastically, obtaining a uniform randomdistribution.
2. Assignment
In order to calculate the electric field that results from thecharge distribution in the domain, the charges first have to beassigned to the grid. This is done using the cloud-in-cell (CIC)charge assignment scheme, in which the charges are repre-sented by finite-volume charge clouds [41]. The CIC schemeresults in a smoother field distribution than the commonlyused nearest-grid-point scheme, where, as the name indicates,the total charge of each particle is assigned to the nearest gridpoint. In the CIC scheme, for each particle in a given grid cell,a portion of the particle’s charge is assigned to each one of thecell’s grid points. The fraction of the charge, or weight, of aparticle located at ( x, y, z ) on the n -th grid point with position ( x n , y n , z n ) is given by w n = (cid:18) − | x n − x | ∆ x (cid:19) (cid:18) − | y n − y | ∆ y (cid:19) (cid:18) − | z n − z | ∆ z (cid:19) , (7)where ∆ x , ∆ y , and ∆ z are the grid cell dimensions along x , y , and z . For carriers in graphene, the weights are non-zeroonly in the plane of the sheet, since their motion is confinedto that plane. B. Field initialization
An initial electric field distribution satisfying Gauss’s lawcan be calculated as the gradient of the electrostatic poten-tial Φ , which is obtained by solving Poisson’s equation forthe initial charge distribution. The initial charge distribution ρ ( i, j, k ) at grid point ( i, j, k ) is given by the sum of theweights of all the charges in the cells surrounding that point, ρ ( i, j, k ) = N X c=1 q c w c ( i, j, k )∆ x ∆ y ∆ z , (8)where N represents the total number of charges in the gridcells surrounding ( i, j, k ) , q c is the charge of particle c , and w c ( i, j, k ) is the weight of particle c at point ( i, j, k ) , givenby Eq. (7). Using ρ ( i, j, k ) , we solve Poisson’s equation withthe successive-over-relaxation (SOR) method [42] to get the y [ n m ] Gaussian correlationfunction fitSpatial autocorrelationfunction00.20.40.60.81 a r b . un i t s FIG. 4. (a) Example of a numerically generated clustered impuritydistribution with a sheet density of × cm − . (b) The spatialautocorrelation function (SACF) of the numerically generated dis-tribution (orange circles) is fitted with a Gaussian correlation func-tion (blue line) of the form exp( − r /λ ) , where λ is the correlationlength. The average cluster size/correlation length estimated fromthe full width at half maximum (FWHM) of the SACF is . ,while that estimated from a Gaussian fit is . . electrostatic potential Φ( i, j, k ) . We use periodic boundaryconditions at the four bounding planes perpendicular to thegraphene layer and the Dirichlet boundary condition with avanishing potential on the top and bottom planes. The initial electric field distribution is then calculated from E x (cid:18) i + 12 , j, k (cid:19) = − [Φ( i + 1 , j, k ) − Φ( i, j, k )] / ∆ x, (9a) E y (cid:18) i, j + 12 , k (cid:19) = − [Φ( i, j + 1 , k ) − Φ( i, j, k )] / ∆ y, (9b) E z (cid:18) i, j, k + 12 (cid:19) = − [Φ( i, j, k + 1) − Φ( i, j, k )] / ∆ z. (9c)In addition to the initial electric field distribution, we also needthe grid-based electric field contribution of the impurity ions,which is subtracted from the MD fields in the vicinity of theions to avoid double counting [14]. Since the ions are fixedin space, the grid-based contribution from the ions does notchange in time; therefore we only have to calculate it once be-fore commencing with the time-stepping. We rely on the lin-earity of Poisson’s equation and simply use the solution for asingle ion instead of solving the equation for each ion. More-over, we use the MD contribution to the force on a carrieronly for charges within a × × -cell volume surround-ing the given carrier and therefore only ions near the interfacethat are present within this cell volume are treated with MD.This approach is illustrated for a 2D grid in Fig. 5(a)–(e) anddescribed in further detail below.We start by placing a single ion at grid point a near thecenter of the domain, solve Poisson’s equation and calculatethe electric field, shown in Fig. 5(a) with orange arrows. Fromthe complete field solution, we store the field values in thevicinity of cell abcd , marked in Fig. 5(b)–(d) by the brownarrows. With these stored fields, we can determine the short-range grid-based fields for an ion at points a , b , c and d simplyby correctly shifting the stored fields to different points on thegrid. For example, as shown in Fig. 5(c), marked by dark bluearrows are local fields for a single ion located at grid point a .For the purpose of illustration, the grey dashed box representsa secondary cell and marks the relative position of the ion.Now, by moving the grey dashed box to the cell above abcd ,the position of the ion (at point a ) relative to the grey boxis equivalent to that of grid point c relative to the cell abcd .Thus, using the original position of the ion (at grid point a ),the local fields due to an ion at grid point c (marked by darkblue arrows) can be found simply by shifting the grey dashedbox, as shown in Fig. 5(d). Once, the local fields due to asingle ion at all the corners of the grid cell abcd are known,the field due to an ion at any arbitrary position within cell abcd [Fig. 5(e)] can be calculated as a weighted sum of theseshifted fields, where the weights are given by Eq. (7) for thatimpurity ion.This procedure of storing and shifting short-range gird-based fields [14] is applicable only within a single uniformmedium. However, in order to have correct continuity in thefields near the interface of air, graphene and SiO , the poten-tials in the respective mediums are required. A discontinu-ity in the fields at the interface results in residual fields thatproduce unphysical dc current components that persist evenwithout any externally applied field. Therefore we solve Pois-son’s equation three times – in air, graphene, and SiO – forthe impurity ions near the interface of graphene and the sub-strate. The grid-based electric field contributions of the ionsnear the interface are then calculated by taking a combinationof the appropriate fields from the three mediums, as shown inFig. 5(f).Calculating the grid-based field contributions from carriersusing the above method is more complicated than for ions be-cause of carrier motion. An implementation similar to the ionswould be computationally burdensome as it would require re-calculating the weights and shifting of the stored local fieldsat each time step. Therefore, we use the corrected-Coulombscheme [14] for determining the local grid-based field contri-butions of the carriers. In this scheme, a carrier is placed at afixed location at a grid point. The grid-based field contributionis found from the solution to the Poisson’s equation. The samePoisson’s solution is used to determine the fields when the car-rier is displaced by a small amount within the grid cell usingthe new weights of the carriers. The corrected-Coulomb fieldsare then calculated by subtracting these grid-based fields fromthe MD fields at the carrier locations. Although this method isnot exact, the reduction in computational burden is significantand errors introduced have very little impact on the results[14]. C. Current density calculation
At the initialization stage, Poisson’s equation is solvedto obtain the electric field profile consistent with the initialcharge distribution. Thereafter, the continuity equation needsto enforced at each time step because FDTD and the continu-ity equation together ensure that Gauss’s law remain satisfiedthroughout the simulation [14, 31]. We achieve this here byusing the Villasenor-Buneman method [43] to calculate thecurrent density, ~J , from the change in the carrier position overa time step and assign it to the grid using the same CIC schemeas before. The current densities are defined at the same loca-tions on the grid as the corresponding fields (see Fig. 3) andare given by [19] J x (cid:18) i + 12 , j, k (cid:19) = N X c=1 q c ∆ x ∆ yt g x f − x i ∆ t (cid:18) − y f + y i y + j (cid:19) , (10a) J y (cid:18) i, j + 12 , k (cid:19) = N X c=1 q c ∆ x ∆ yt g y f − y i ∆ t (cid:18) − x f + x i x + i (cid:19) , (10b)where subscripts f and i represent the final and initial posi-tions, respectively, and t g is the thickness of graphene. Al-though the carriers move in a 2D plane, the current densitymust have the units of Am − for the sourcing of FDTD fields.Therefore, to calculate the current density in the correct units,we divide by t g . We use t g ≈ ˚A to represent the approxi-mate thickness of the graphene layer [44, 45]. Motion of car-riers into the neighboring grid cell is treated by dividing the path into sections, such that the motion in each cell is treatedindividually. D. Lorentz force calculation
We require the fields at the location of the carriers to calcu-late the force acting on the carriers drifting in EMC, Eq. (1).To determine the fields at the location of a carrier that is foundinside a given grid cell, for each of the 8 grid points of thatcell we first average the fields on the grid faces and grid lines(Fig. 3) surrounding it. From the fields at each grid point, wecan use the same CIC weights of each carrier to interpolatethe fields and obtain the values at the carrier’s actual location.The total electric field E cT that accelerates carrier c in EMCthus consists of the following contributions: E cT = P E nFDTD w cn + P N c ′ =1c ′ =c (cid:16) E c - c ′ MD − E c ′ grid (cid:17) + P M i=1 (cid:16) E c - iMD − E igrid (cid:17) , (11)where n enumerates the corners of a grid cell containingcarrier c , and w n are the CIC weights. N and M are the totalnumbers of carriers and ions, respectively, within the gridcells surrounding c . E nFDTD is the FDTD field contribution, E c - c ′ MD and E c - iMD are the MD field contributions due to carrier-carrier and carrier-ion interactions respectively, while E c ′ grid , E igrid are the local grid-based field contributions due to carri-ers and ions, respectively. The MD fields are pre-calculatedbefore starting the time-stepping loop for a dense mesh in thereal space and k -space within a volume of × × grid cellsand stored in look-up tables. At any given time step, we thenlook up the MD fields based on the pairwise differences be-tween carrier positions. IV. EXAMPLE: CONDUCTIVITY OF GRAPHENE
In this section, we use the coupled EMC/FDTD/MD solverto calculate the dc and ac conductivity of graphene. The com-plex conductivity is computed from the spatially averaged val-ues of the current density ˆ J ( ω ) and electric field ˆ E ( ω ) phasorsas σ ( ω ) = ˆ E ( ω ) · ˆ J ∗ ( ω ) | ˆ E ( ω ) | . (12)The phasor quantities are calculated at each grid point in thegraphene plane by using on-the-fly discrete Fourier transformof the time-dependent vector components after a steady statehas been reached, then spatially averaging the components foruse in (12). For example, the current density phasor is givenby ˆ J ( ω ) = T s X n = t s ~J [cos (2 πf n ∆ t ) − i sin (2 πf n ∆ t )] , (13) (a) a bc d a bc d(b) a bc d(c)a bc d(d) a bc d(e) A i r g r aphene S i O (f) FIG. 5. (a)–(e) Illustration of the calculation of the grid-based electric field from an impurity ion. (a) A single ion (dark blue circle) placed atone of the corners of the cell abcd and the corresponding electric field, marked by orange arrows, as calculated from the solution to Poisson’sequation. (b) A smaller array of field values in the grid cells close to point a , marked by brown arrows, used for subsequent calculations. (c)The local fields (marked by dark blue arrows) due to a single ion at point a found by properly shifting the stored fields around appropriate gridpoints. (d) The local fields due to a single ion at point c using the solution for the ion at point a . (e) Local fields for a single charge locatedat an arbitrary position within the grid cell abcd , calculated using a weighted sum of the shifted potentials that correspond to a single ion ateach of the corners. The weights are given by Eq. 7. (f) Side view of the local 3D grid-based fields due to an impurity ion in the substrate nearthe interface of graphene, calculated as the weighted sum of the shifted fields in air (blue), graphene (green), and the SiO substrate (orange).Poisson’s equation for a single ion is solved three times, for calculating fields due to an ion close to the interface, in air, graphene, and SiO .The weighted sum is calculated based on the combination of fields in the appropriate medium for each impurity ion. where t s is the time to reach a steady state, T s is the total sim-ulation time, f is the frequency of external excitation ( f = 0 for dc excitation), and ~J = J x ~x + J y ~y is the current densitycalculated from Eq. (10). Being that graphene is a 2D ma-terial, its conductivity is typically measured and presented inthe units of e /h , the quantum of conductance. In order toconvert the conductivity calculated in Eq. (13) to those units,we multiply by a factor of t g h/e , where, as before, t g ≈ ˚A is the effective thickness of the graphene electron system[44, 45]. A. dc Conductivity
We calculate the dc conductivity of graphene as a func-tion of the carrier sheet density, shown in Fig. 6, for an im-purity sheet density of × cm − with a uniform ran-dom distribution (blue squares) and a clustered distribution(red diamonds; correlation length of
40 nm ) and compareit with the conductivity of impurity-free graphene (black cir-cles). These results reproduce important features of the con-ductivity vs. carrier density curve observed in experiment[46]. The curve displays a sublinear increase at high carrierdensities ( > × cm − ) for “clean” graphene. More- over, for a given sheet density of charged impurities, the clus-tered impurity distribution results in lower conductivity thanthe uniform random one. This behavior has also been ob-served in experiment [47] and predicted in the calculationsof carrier-impurity scattering rates with a structure factor de-scribing correlations [26]. Our EMC/FDTD/MD simulationmakes no assumptions about the screening length or struc-ture factor, and uses real-space impurity positions and the cor-responding carrier-impurity interactions to calculate the con-ductivity. Our results also show a flattening of the conductiv-ity curve near the Dirac point for clustered impurity distribu-tions, similar to that observed in conductivity measurementsinvolving intentional potassium doping [46]. B. ac Conductivity
The frequency-dependent ac conductivity, shown in Fig. 7,is calculated for the same impurity density and distributionsas the dc case (Fig. 6). Here we use a carrier density of × cm − . The frequency of the external excitation isvaried from
500 GHz to
13 THz . In this range, carrier trans-port is dominated by intraband processes [48] and is capturedvery well in our simulation. These results are in line with n [cm −2 ]0510152025303540 σ [ e / h ] Impurity freeClustered impurity distributionRandom impurity distribution
FIG. 6. dc conductivity of supported graphene as a function of thecarrier density. Black circles denote the results for impurity-freecase, blue diamonds for a uniform random distribution, and red di-amonds for a clustered distribution (
40 nm average cluster size) ofcharged impurities with a sheet density of × cm − . The blackline is a linear fit to the low-density ( < × cm − ) part of theimpurity-free curve. R e { σ } [ e / h ] FIG. 7. Frequency-dependent ac conductivity of supported graphenefor the same charged impurity distribution as in Fig. 6. Carrier den-sity is assumed to be × cm − . experimental measurements of frequency-dependent conduc-tivity [48, 49]. For frequencies greater than , our re-sults show that the total impurity density and distribution donot affect the conductivity of graphene. However, for lowerfrequencies ( < ), there is a significant dependence ofconductivity on the impurity density and distribution. (As ex-pected, the low-frequency conductivity limit obtained from ac calculations is very close to the values calculated in the dc simulations.) V. SUMMARY
We have presented the implementation and application ofthe coupled EMC/FDTD/MD simulation technique to carriertransport in supported graphene in the presence of charged im-purities. We have described the constituent techniques, as wellas the important steps for their self-consistent coupling, suchas charge initialization and assignment to the grid, field ini-tialization, current density calculation based on particle mo-tion in the gird, and avoiding double counting of the fieldsfrom FDTD and MD. The general implementation can alsobe applied to transport simulations of other 2D or quasi-2Dmaterials.We have demonstrated the use of the EMC/FDTD/MDmethod by calculating the dc and ac conductivity of supportedgraphene. The calculated dc conductivity as a function of thecarrier density reproduces the important features observed inexperiments, such as the sublinear increase at high carrier den-sity in clean samples [46] and flattening of the curve near theDirac point for clustered impurity distribution [47]. The cal-culated ac conductivity agrees with experimental observations[48, 49]. [1] S. Das Sarma, S. Adam, E.H. Hwang, E. Rossi, Rev. Mod. Phys. , 407 (2011) [2] M. Fontana, T. Deppe, A.K. Boyd, M. Rinzan, A.Y. Liu,M. Paranjape, P. Barbara, Sci. Rep. (1634) (2013) [3] P. Zhang, E. Tevaarwerk, B.N. Park, D.E. Savage, G.K. Celler,I. Knezevic, P.G. Evans, M.A. Eriksson, M.G. Lagally, Nature , 703 (2006)[4] Y.M. Lin, C. Dimitrakopoulos, K.A. Jenkins, D.B. Farmer, H.Y.Chiu, A. Grill, P. Avouris, Science (5966), 662 (2010)[5] R. B., R. A., B. J., G. V., K. A., Nat. Nano. , 147 (2011)[6] F. Bonaccorso, Z. Sun, T. Hasan, A.C. Ferrari, Nat. Photon. ,611 (2010)[7] Z. Yin, H. Li, H. Li, L. Jiang, Y. Shi, Y. Sun, G. Lu, Q. Zhang,X. Chen, H. Zhang, ACS Nano (1), 74 (2012)[8] L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H.A.Bechtel, X. Liang, A. Zettl, Y.R. Shen, F. Wang, Nat. Nano. ,630 (2011)[9] Y. Liu, X. Dong, P. Chen, Chem. Soc. Rev. , 2283 (2012)[10] S. Sonde, F. Giannazzo, C. Vecchio, R. Yakimova, E. Rimini,V. Raineri, Appl. Phys. Lett. (13), 132101 (2010)[11] S. Fratini, F. Guinea, Phys. Rev. B , 195415 (2008)[12] S. Adam, E.H. Hwang, V.M. Galitski, S. Das Sarma, P. Natl.Acad. Sci. USA (47), 18392 (2007)[13] E.H. Hwang, S. Adam, S. Das Sarma, Phys. Rev. Lett. ,186806 (2007)[14] K.J. Willis, S.C. Hagness, I. Knezevic, J. Appl. Phys. (6),063714 (2011)[15] K. Tomizawa, Numerical Simulation of Submicron Semicon-ductor Devices . Electronic Materials and Devices Library(Artech House, 1993)[16] C. Jacoboni, P. Lugli,
The Monte Carlo Method for Semicon-ductor Device Simulation . Computational Microelectronics(Springer, 1989)[17] J. Ayubi-Moak, S. Goodnick, S. Aboud, M. Saraniti, S. El-Ghazaly, J. Comput. Elec. (2-4), 183 (2003)[18] J. Ayubi-Moak, S. Goodnick, M. Saraniti, J. Comput. Elec. (4), 415 (2006)[19] K. Willis, J. Ayubi-Moak, S. Hagness, I. Knezevic, J. Comput.Elec. (2), 153 (2009)[20] K.J. Willis, S.C. Hagness, I. Knezevic, Appl. Phys. Lett. (6),062106 (2010)[21] P. Lugli, D.K. Ferry, Phys. Rev. Lett. , 1295 (1986)[22] D.K. Ferry, A.M. Kriman, M.J. Kann, R.P. Joshi, Comput.Phys. Commun. (1), 119 (1991)[23] C. Wordelman, U. Ravaioli, Electron Devices, IEEE Transac-tions on (2), 410 (2000)[24] D. Vasileska, H. Khan, S. Ahmed, J. Comput. Theor. Nanosci. (9), 1793 (2008)[25] K.J. Willis, S.C. Hagness, I. Knezevic, Appl. Phys. Lett. (12), 122113 (2013)[26] Q. Li, E.H. Hwang, E. Rossi, S. Das Sarma, Phys. Rev. Lett. , 156601 (2011) [27] C. Jacoboni, L. Reggiani, Rev. Mod. Phys. , 645 (1983)[28] N. Sule, I. Knezevic, J. Appl. Phys. (5), 053702 (2012)[29] K.M. Borysenko, J.T. Mullen, E.A. Barry, S. Paul, Y.G. Se-menov, J.M. Zavada, M.B. Nardelli, K.W. Kim, Phys. Rev. B , 121412 (2010)[30] A. Konar, T. Fang, D. Jena, Phys. Rev. B , 115452 (2010)[31] A. Taflove, S. Hagness, Computational Electrodynamics: TheFinite-Difference Time-Domain Method . The Artech Houseantenna and propagation library (Artech House, Incorporated,2005)[32] K. Yee, IEEE T. Antenn. Propag. (3), 302 (1966)[33] J.A. Roden, S.D. Gedney, Microw. Opt. Techn. Lett. (5), 334(2000)[34] D. Rapaport, The Art of Molecular Dynamics Simulation (Cam-bridge University Press, 2004)[35] B.J. Alder, T.E. Wainwright, J. Chem. Phys. (2), 459 (1959).doi:10.1063/1.1730376[36] A.M. Kriman, M.J. Kann, D.K. Ferry, R. Joshi, Phys. Rev. Lett. , 1619 (1990)[37] R.P. Joshi, A.M. Kriman, M.J. Kann, D.K. Ferry, Appl. Phys.Lett. (21), 2369 (1991)[38] P. Gori-Giorgi, F. Sacchetti, G.B. Bachelet, Phys. Rev. B ,7353 (2000)[39] 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, C.A. Richter, ACS Nano (11), 9144 (2011)[40] T. Fang, A. Konar, H. Xing, D. Jena, Appl. Phys. Lett. (9),092109 (2007)[41] S. Laux, IEEE Trans. Comput-Aided Des. Integr. Circuits Syst. (10), 1266 (1996)[42] W.H. Press, Numerical Recipes : The Art of Scientific Comput-ing (Cambridge University Press, 1989)[43] J. Villasenor, O. Buneman, Comput. Phys. Commun. (2–3),306 (1992)[44] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science (5887), 385(2008)[45] M.C. Lemme, T. Echtermeyer, M. Baus, H. Kurz, IEEE Electr.Device L. (4), 282 (2007)[46] J.H. Chen, C. Jang, S. Adam, M.S. Fuhrer, E.D. Williams,M. Ishigami, Nat. Phys. (5), 377 (2008)[47] J. Yan, M.S. Fuhrer, Phys. Rev. Lett. , 206601 (2011)[48] H. Choi, F. Borondics, D.A. Siegel, S.Y. Zhou, M.C. Mar-tin, A. Lanzara, R.A. Kaindl, Applied Physics Letters (17),172102 (2009)[49] J. Horng, C.F. Chen, B. Geng, C. Girit, Y. Zhang, Z. Hao,H.A. Bechtel, M. Martin, A. Zettl, M.F. Crommie, Y.R. Shen,F. Wang, Phys. Rev. B83