Multi-scale electronics transport properties in non-ideal CVD graphene sheet

In this work, we benchmark non-idealities and variations in the two-dimensional graphene sheet. We have simulated more than two hundred graphene-based devices structure. We have simulated distorted graphene sheets and have included random, inhomogeneous, asymmetric out-of-plane surface corrugation and in-plane deformation corrugation in the sheet through autocorrelation function in the non-equilibrium Green’s function (NEGF) framework to introduce random distortion in flat graphene. These corrugation effects inevitably appear in the graphene sheet due to background substrate roughness or the passivation encapsulation material morphology in the transfer step. We have examined the variation in density of state, propagating density of transmission modes, electronic band structure, electronic density, and hole density in those device structures. We have observed that the surface corrugation increases the electronic and hole density distribution variation across the device and creates electron-hole charge puddles in the sheet. This redistribution of microscopic charge in the sheet is due to the lattice fields’ quantum fluctuation and symmetry breaking. Furthermore, to understand the impact of scattered charge distribution on the sheet, we simulated various impurity effects within the NEGF framework. The study’s objective is to numerically simulate and benchmark numerous device design morphology with different background materials compositions to elucidate the electrical property of the sheet device.

Since the isolation of single-layer graphene in 2004, graphene has shown promising research results in laboratories across the globe 1 , Furthermore, tremendous progress was made in fabricating and characterization techniques to fabricate large-area graphene sheet devices 2,3 . However, large-scale commercial production and market adoption are still very far from the horizon. The possible utilization of graphene sheets is in various applications from the transparent electrode, electronics, telecommunication equipment, and internet of thing (IoT) sensors. To realize this commercial adoption in the electronics and optoelectronics industry, the main challenge for graphene is in the controllability, homogeneity, and uniformity in the electronic property in the mass roll-to-roll production process. The growth and transfer process, such as temperature variation and gas flow fluctuations in the chemical vapor deposition reactor, give rise to long-range spatial variation in the electronics performance parameter of the large-area graphene sheet. It prevents the seamless stitching of graphene single grain, which gives rise to imperfect lattices and grain boundary interfaces, causing a decrease in conductance [4][5][6][7][8] . Moreover, the geometrical circumstance of the graphene sheet transfer scheme gives rise to spatial variation in the doping density. Therefore, for homogeneous and uniform graphene sheets, defect density is an important measure. However, we must consider the graphene sheet's average defect density and spatial defect distribution profile to quantify the effect 9 . Moreover, all those issues, as mentioned above, are a challenge to scale up the production. There are various methods proposed to fabricate the graphene sheet on an industrial scale, such as exfoliation, the epitaxial method 10 , chemical vapor deposition (CVD) 11 , Roll-to-Roll plasma-enhanced chemical vapor deposition (PECVD) method [12][13][14][15] . However, the synthesis of graphene on metal catalysts by CVD has been the most scalable process until now. Besides, all these process steps have their underline chemistry and influence on the quality of the graphene sheet 16 . Therefore, the electrical quality of CVD/PECVD grown graphene sheets varied broadly from process to process and batch to batch. Besides, the transfer and isolation process are bottlenecks in the production scale and incorporate defects in the sheet and degrade the sheet's electrical quality [17][18][19][20] . These nanoscale atomistic effects during the growth and transfer steps influence the vital electrical properties of the sheet through variation in background carrier density, carrier mobility, and sheet resistivity. Electrical property mapping is critical for the process optimization and quality control of large-area graphene sheets in Roll-to-Roll

Method
At the nano-scale atomic composition, crystal symmetry and the spatial disorder affect the material's bulk properties as quantum mechanics effects come into the picture and modify the electronic-phononic structure. Therefore, atomistic simulations are more appropriate to model the quantum device's electronic properties in the entire Brillouin zone 29 . The electronic band structure derives from the tight-binding method. The method is similar to the linear combination of atomic orbitals (LCAO) used to construct molecular orbitals 30,31 . However, in the tight-binding approximation, electron-electron interactions of the orbital are neglected, but it gives an excellent approximation to electronic band structure. Furthermore, for more rigorous treatment, the Hubbard model was employed. The graphene is a two-dimensional sheet of carbon atoms arranged in a hexagonal lattice. 2S, 2P x , and 2P y orbitals of carbon atom in graphene are SP 2 hybridize, resulting in strong σ bonds where each carbon atom is bonded to all of its three neighbors carbon atoms. The P z or π orbital of the carbon atoms defines the low-energy electronic bandstructure in graphene.
The primitive lattice vectors in graphene are, The graphene primitive unit cell comprises two carbon atoms, and each atom has one valence 2P z orbital. At the C1 atoms site valence φ 2p z1 orbital is centered and at the C2 atoms site valence φ 2p z2 orbital is centered, The tight binding wavefunction for graphene is, After multiplying from the left by both of P z orbitals to time-independent Schrödinger equation and integrating overall space, the dispersion relation is, In tight-binding approximation, on-site and nearest-neighbor matrix elements are retained, and all the other terms are assumed too small enough to ignore the equation, (2) ψ � k (� r) = 1 √ N h,j e i h � k·� a 1 +j � k·� a 2 c 1 φ 2p z1 � r − h� a 1 − j� a 2 + c 2 φ 2p z2 � r − h� a 1 − j� a 2 .
(3) �φ 2p z1 |Ĥ|ψ k � = E�φ 2p z1 |ψ k �; �φ 2p z2 |Ĥ|ψ k � = E�φ 2p z2 |ψ k �. www.nature.com/scientificreports/ where Graphene electronic band structure is defined in the tight-binding approximation while considering interactions up to third-nearest neighbors as 32,33 , where Ĥ 1 is the first nearest neighbor Hamiltonian and Ĥ 3 third nearest neighbor Hamiltonian. In the second quantization language, Hamiltonians are in real space representation is expressed as creation and annihilation operator acting on the π state on each carbon atom site as follow, where c † creation and c annihilation operators, and summation runs over entire n lattice point, and l is first and m is third nearest neighbor site of n lattice point. First nearest-neighbor hopping parameter t n,l = −2.74 eV and third nearest neighbor t n,m = −0.3 eV . We have also passivated the dangling bond at the device edge by hydrogen passivation treatment in the device simulation 34 . In the tight-binding framework dangling bonds at the surface or edge is passivated by primarily two numerical methods. In the first strategy, passivation atoms are implicitly incorporated without distinguishing passivation atom types and add a passivation potential to the dangling bonds' orbital energies. This method works well with a relatively large system, arbitrary crystal structures, and hybridization symmetries. Furthermore, with appropriate parameters, it is applied to any passivation scenario 35 . The second approach is explicit, including the passivation atoms and their coupling to the surface or edge atoms Hamiltonian matrix and limited to small molecules and systems 36 . In this explicit treatment, ab-initio results for different passivation atoms used to fitting targets. The unsaturated dangling bonds at the edge or surface will result in edge or surface states at the electronic band structure. These unwanted states iron out by coupling the hydrogen passivation atoms to the surface's or edge's unsaturated dangling bonds in the device. We have used the P − D orbital tight-binding model, which represents the edge effects by explicitly including the passivated hydrogen in the Hamiltonian matrix. The carbon atom is represented by three P z , D yz , and D zx orbitals. The simple single orbital P z tight-binding model works well for the two-dimensional graphene sheet. More further details about the tight-binding framework address in the article's supplementary information. In 1960 Keldysh, Kadanoff, and Baym first developed the non-equilibrium Green's function formalism (NEGF) 37,38 . The adaption of NEGF formalism to semiconductor devices first demonstrated by Lake and Datta [39][40][41][42] , in 1992 and later in 2002 by Wacker 43 . Through NEGF formalism, the time evolution of many-body quantum fields in thermodynamic equilibrium and non-equilibrium state can be investigated [44][45][46][47][48][49] . These quantum fields constitute carriers such as electrons, phonons, spin, and electric-field in semiconductor devices [50][51][52][53][54][55][56][57][58][59][60] . The theoretical framework of Non-equilibrium Green's function formalism with detailed reference is given in the supplementary information. Here we briefly summarize the NEGF formalism on a discrete device grid as, where the bold symbol represent the vector space quantities, G R retarded, G A advanced, G < lesser, G > greater Green's functions, R retarded, A advanced, � < lesser, � > greater self-energies of interaction between various www.nature.com/scientificreports/ quantum field, A is spectral function related to the spectral peak of Green's functions in energy spectrum, k t = (k x , k y ) is parallel momentum in vector space, E is energy space, δ is Kronecker delta in discrete energymomentum space where Green's functions realized it's value, system Hamiltonian h nm = drφ * n (r)H 0 (r)φ m (r) with hermitian h * ml = h lm property, φ is tight-binding orbital describing the system from Eq. (3), r is position space, and n, m, l, v are discrete lattice point of device where Green's propagator travel in coupled system of Eq. (7). From Eq. (7) the stationary state solution of non-equilibrium Green's functions, the Density of state D(E) is defined by taking the trace of lesser G < Green's functions as, The carrier density N(r, t) in the non-equilibrium Green's functions formalism is defined as, In the stationary solution regime from Eq. (7) of NEGF in the (z) directional transport, eigenfunction expansion as orthonormal eigenstate φ k t n (z) , φ * k t m (z) , the carrier density as, For the contacted device with an external source and continuous supply of electrons in the non-equilibrium configuration G < nm (k t ; E) is defined at different points in the device. Where l belongs to the active part of the device, m belongs to any point in the contacts with defined carriers Fermi distribution, and n belongs to the interface between both regions. Using the corresponding boundary conditions for h ml,lm and with equilibrium contacts assumption as in bulk contacts, the Fermi distribution is in equilibrium. By using the fluctuation-dissipation theorem, the charge exchange rate Ŵ Contact where A mm 1 is spectral function correspond to lattice indices point at mm 1 and h lm , h m 1 l 1 is corresponding points at surface Hamiltonian in the boundary condition. From Eq. (10), and using the argument of Eq. (11), The Transmission probability T at (k t ; E) is defined as, where in the non-interacting active part of the device indices l, l 1 , l 2 , l 3 run covering all the points. And, Ŵ L ll 1 is connected towards the left contact and Ŵ R l 2 l 3 connected through the right contact. From the argument of Eqs. (11) and (12), the total number of the transmissible propagating modes of the wave-function in the device is defined as Density of Modes or Modes Density M at energy (E) is, Though NEGF formalism is computationally intensive to solve compared to classical drift-diffusion and semiclassical Boltzmann transport approach and has simulation overhead, it is a complete quantum treatment. Therefore, it gives excellent detail of insight into device operation in spatial resolution. We have employed Nano-Electronics Modeling Tools (NEMO5) , a multi-scale quantum transport kernel for nano-electronic device modeling 61,62 . Its modular architecture parallelizes in a five-level message passing interface (MPI) in the position, momentum, energy, bias, and random seeding space. Due to the modular architecture of the transport kernel, various physical models add up and extend into it, the detailed discussion provided in the article's supplementary information.

Result and discussion
We have employed multi-scale, multi-physics-based bottom-up, non-equilibrium Green's function mechanismbased quantum transport simulation techniques to investigate various device aspects and deduce the observable and measurable at the final stage in the computationally simulated devices. Moreover, most reported publications in the graphene and two-dimensional material field calculate electrical parameters of coherent transport limit and leave scattering interaction behind as a residual. The scattering nature and origin depend on graphene's different fabrication process steps, the device orientation and substrate, and surface encapsulation in the device preparation step. The complete stepwise description of the theoretical framework provides in the article's method section and supplementary information. This industrial research work was envisioned to show the opportunity of improving the electrical characterization of the roll-to-roll CVD graphene production process. Therefore, our results in this work are directed to corrugations and impurities effect in the reasonably clean state-of-art industrial production and transfer process. Here we will discuss the main results from numerical calculation.

Surface corrugation effect.
We have investigated the out-of-plane random, inhomogeneous, asymmetric structural corrugation effect in the electrically open infinite reservoirs contacted geometry of graphene (9) N(r, t) = −i G < (r, t; r, t). www.nature.com/scientificreports/ sheet to examine the large variability across the electric properties of CVD graphene sheet. Graphene is a single atom thin sheet of the carbon atom and needs a substrate to work as a device. The substrate's effect is manifold, such as corrugation in the graphene sheet by the substrate surface profile and scattering through the substrate's impurities 63,64 . Due to thermodynamic instability, the corrugation and ripple have been observed in the free-standing and substrate-born graphene sheet devices. There is modest work reported on deformed graphene. Most of the previous work concentrated on a uniform, uniaxial or biaxial strain and centrosymmetric deformations [65][66][67][68][69] . Katsnelson et al. derived the mathematical formulation for microscopic defects and generic random potential by ripples for the Dirac fermions from the semi-classical Boltzmann transport view 70 . The fluctuating interatomic distances and angles between carbon bonds are described by in-plane and out-of-plane atomic displacements deformation tensor in the corrugated graphene. This work tried to elongate their original idea in the non-equilibrium Green's functions formalism (NEGF) based quantum transport domain. Based on the substrate material type and fabrication process, the corrugation profile has considerable variation from 0.30 nm to 0.005 nm when placed on the substrate, graphene follows the underlying substrate's surface morphology, and these ripples induce stress in the sheet. Therefore, it will influence the electronic property of the graphene sheet. Variation in bonding lengths and angle between the carbon atoms in the corrugated graphene sheet depends on substrate morphology. In tight biding calculation, it modulates the hopping parameters in the hamiltonian. From Harrison's model, the hopping parameter is inversely proportional to the square of the bond length 71 . Due to corrugation, p z -orbitals are also bent, but hopping parameter variation due to bond length variation is more prominent compared to orbital bending 72 . The surface corrugation effect in the graphene sheet modeled through the gaussian auto-correlation function in the transverse x and y-direction as, where δh is the root mean square amplitude of fluctuations height. The correlation length l x and l y are in spatial domain morphology repetition length scale. The surface corrugation effect in the lattice is created in simulation by taking the Fourier transform of auto-correlation function from Eq. (14), applying a random phase, and taking an inverse Fourier transformation. In experiment, depending upon the fabrication and transfer process, and underline substrate morphology, the correlation length varies from 2 to 30 nm and corrugation amplitude 0.30 nm to 0.005 nm from the SiO 2 to h-BN substrate. To numerically simulate these substate effects and variation, we have calculated the density of state from Eq. In the simulated device, the primitive unit cell has two atoms per cell, and a finite element mesh simulates a total of 200 atoms with a domain size of 800 finite element mesh grid points. The P-D tight-binding model contains three orbitals, namely carbon P z , and carbon-hydrogen passivated D yz , D xz orbitals. Therefore total degree of freedom in hamiltonian is 600 variable-sized orbitals. The K ′ and M ′ are high symmetric point corresponds to the folded reduce BZ-zone of graphene supercell. In the article's supplementary information, from Fig. 1a-ad is also rendered in the full-text page resolution in the corresponding Supplementary Figs. 1-30 for the reader's reference. We have observed from Fig. 1f-j, the density of states fluctuates from the flat to deformed structure progressively gives rise to additional states to scattered and decreases conductance. Also, in the deformed sample, sharp peaks appear in the density of states at lower energies compared to undeformed flat samples. We have also observed the disorder-induced Dirac point in the density of state shifting from zero energy point towards the right due to a shift in Fermi-energy also visible in from band structure plots Fig. 1z-ad as beforehand announced and reported in Raman spectra of the experimental graphene sheet 73 . Moreover, the sharply peaked Van-Hove singularity point slowly shifts from low energy towards higher energy and slowly smoothing at these higher energies with the increased deformation in the sample. Even a tiny randomized inhomogeneous in-plane or out-ofplane deformation strain will break the local sublattice and inversion symmetry in the one atom thin graphene. It will give rise to an inhomogeneous pseudo magnetic field that produces the additional density of the state's low-energy peaks and increases the scattering probability with a deteriorating conductance. Moreover, inhomogeneous pseudo field incipient a non-equivalent charge density localization and distribution across the asymmetric sublattices of deformed regions. Furthermore, creating the electron-hole puddles in the graphene sheet can be seen from the numerical calculation from the electronic density Fig. 1p-t and hole density Fig. 1u-y. This will be more visible next for a larger size 100x100 graphene supercell sample simulation for electronic density Figs. 2g-3h and hole density Figs. 2i-3i. In the deformed sheet where sublattice and inversion symmetry is locally broken, energy bandgap opening is not visible from Fig. 1z-ad. However, as discussed in other studies, a sizeable large gaussian deformation can introduce an energy bandgap in the corrugated sheet 74 . From density of modes M(E) Fig. 1k-o, we have observed the first transmission plateau weakly affected for the small deformations sample. Figure 1. Device structure 10x10 graphene supercell (a-e) corrugate-5pm to −20pm, Density of state (f) flat to (j) corrugate-5pm to −20pm, Electronic density of modes M(E) in x-y Direction (k) flat to (o) corrugate-5pm to -20pm, Electronic density (p) flat to (t) corrugate-5pm to −20pm, Hole density (u) flat to (y) corrugate-5pm to −20pm, and Electronic band structure (z) flat to (ad) corrugate-5pm to −20pm. www.nature.com/scientificreports/ However, progressively growing deformations induce a random fluctuation in the higher energy transmission plateau by reducing the plateau's energy steps size following the density of state trajectories. It will mix up the various transmitting modes and wash out the clearly defined energy step size of the flat ballistic transmission modes. Therefore, it will lead to injected propagating electrons in graphene sheet with clearly defined energy that will not propagate in their own modes but have the further opportunity to interact with another mode through newly available deformed state and gain or release energy and transit between different conducting modes. It will provide them access to the new density of state to scattered further before reaching to drain contact from source contact. This verticle energy coupling will give rise to a phase incoherence in the ballistic regime. Overall electrical conductance in an energy window will degrade with increasing deformed samples from the idealistic flat scenario. First, the extra density of the state gives rise to an additional scattering probability, and second these energies mix up of different transmitting modes in the corresponding window. Moreover, degradation from ideal characteristics and fluctuation in density of state and modes density is more severe in the high-energy windows than in the low-energy region, and this trend is more severe for larger corrugated devices. Phenomenologically corrugated deform graphene sheet can be interpreted as a flat graphene sheet with local defect and scattering  www.nature.com/scientificreports/ center at each site and delocalized spread across the sheet. Deformation will introduce variation in density of state and transmission modes density in the entire scattering wavelength spectrum from the idealized flat graphene characteristic. Therefore, high-energy electronic states with small scattering wavelengths have more variation than low energy states. In near-zero energy fields, there is less distortion for more minor deformed structures. As the deformation grows at lattice size, It also starts distorting the low energy idealize characteristic of flat graphene through a large scattering wavelength. Next, we have simulated the large 100x 100 graphene supercell to observe the above-discussed effect in greater spatial detail. In the simulated device, the primitive unit cell has two atoms per cell, and a total of 20000 atoms are simulated by a domain size of 80000 finite element mesh grid points. The total degree of freedom in hamiltonian is 60000 variable-sized orbitals. The density of state from Figs. 2a-3b plotted for the larger device, and as previously discussed, we have observed density of state fluctuation; however, because of the relatively larger size in comparison to the correlation length of deformation density of state variation relatively smoothen in compare to small size cell, Also, The Van-Hove singularity in the DOS curve is not spiked to very high in the curve due to the numerical resolution of computational calculation. Moreover, the slight shift in the Dirac point in the density of state with the corrugation follows the Fermi level calculation. The effective fermi level value follows the band structure calculation.  , d) in-plane deformation-5pm to −10pm, Electronic density (e, f) in-plane deformation-5pm to -10pm, Hole density (g, h) in-plane deformation-5pm to −10pm, Density of state (i, j) in-plane deformation-5pm to −10pm, and Electronic band structure (k, l) in-plane deformation-5pm to −10pm. www.nature.com/scientificreports/ Also, the density of state shifts towards the right is due to a shift in Fermi-energy, as seen in the band structure plots and reported in Raman spectra of the experimental graphene sheet 75 . However, this shift is minuscule in comparison to a small geometry corrugated structure from Fig. 1. From electronic density Figs. 2g-3h and hole density Figs. 2i-3i we represented a real-space description of confined states concentrating in the region where the pseudo magnetic field acquires its minima or maxima. Furthermore, this will provide inhomogeneous charge redistribution in the deformed lattice. Moreover, these features are independent of the lattice orientation of the graphene sheet and within reach of detailed STM characterization of the CVD graphene sheet for industrial-scale qualitative benchmarking. The emergence of a localized and delocalized state in graphene is due to symmetry breaking in the crystal. We have observed electronic bandstructure, from Figs. 2e -3f and the electronic state gets delocalized with more roughness. However, there is no bandgap opening at the K-point for these deformations values, but it can be possible, as discussed previously. Nonetheless, we have observed that energy splitting started at www.nature.com/scientificreports/ calculation from comparing Figs. 2g-j and 3g-i on 25nm X 35nm graphene sheet, we have observed comparable ripples latent size with fluctuating electron-hole charge density puddle in our NEGF calculation. We have also simulated a two-dimensional orthorhombic device structure to investigate whether a four atom per unit supercell will influence similarly or have different characteristics. This geometrical structure can be further tailored to make an armchair and zigzag nanoribbon configuration based upon confinement direction, boundary condition, and edge topology. In the article's supplementary information, simulation results are rendered in the full-text page resolution from Supplementary Figs. 56-80 for the reader's reference. We have observed a similar characteristic in the orthorhombic device structure as in the primitive device. The surface roughness in the corrugated device will shift the Fermi energy level from zero and a shift in the high symmetric points of the Brillouin zone of the ideal graphene sheet as reported in the experimental study of the twisted graphene sheet 75 . The In-plane deformation profile also influences the electronic properties of the narrow-width graphene device and graphene nanoribbon. However, for the wide graphene sheet device on the substrates, the surface corrugation will always impact and fluctuate the device performance characteristics. Furthermore, in the corrugated device structure, the electron density distribution is nonuniform in the sheet implies that scattering rate and carrier conductivity are also spatially varying phenomena in the device. A highly corrugated device will complicate the band diagram with the mixing and split of various branches. In the corrugated device structure, we have observed the density fluctuation with the increase of corrugation compared to the flat configuration. www.nature.com/scientificreports/ In addition, the local density of state variation suggests that the scattering probability should be in the spatial domain, momentum-dependent, and energy-dependent. Moreover, the current path will be highly disordered and nonlinear in such a device and follow the minimum resistance path. In the corrugated structure, density variations from the average density valuation in the flat configuration arise due to random quantum fluctuation due to variable atomic bond length and bond angle in the microscopic description of the material at the atomistic scale, and inversion symmetry and lattice symmetry breaking. These effects influence the macroscopic electronic and charge transport properties of the sample.
In-plane deformation roughness effect. We have investigated the in-plane deformation effect on the electronics properties of the graphene sheet. In-plane, inhomogeneous deformation acts as a scattering source and affects the graphene sheet's transport property 78 . Auto-correlation function is used to statistically model the in-plane deformation effects 79 . In the transport direction of z, the graphene sheet has a uniform ideal width of www.nature.com/scientificreports/ W. At position z 1 and z 2 , the width deviates from the ideal width by δW(z 1 ) and δW(z 2 ) respectively. By autocorrelation function the correlation between δW(z 1 ) and δW(z 2 ) is defined as, By using the argument of correlation function R(z 1 , z 2 ) = R(z 1 − z 2 ) and taking the Fourier transform of autocorrelation function the power spectral density is, The graphene sheet in-plane roughness profile is approximated as Gaussian or exponential nature 80 . By using the exponential nature of the auto-correlation function, where W rms is fluctuation amplitude root mean square. Smoothness of the edge is measure by roughness correlation length l c and z = nδz is total sampling length divide into δz = a cc /2 sampling interval, and a cc is carboncarbon bond length. By using Eq. (17) in the (16) and taking the Fourier transformation, in real space in-plane deformation corrugation defined as the auto-correlation power spectrum is,  www.nature.com/scientificreports/ Finally, the real space roughness is calculated by applying a random phase to Eq. (18) and taking inverse Fourier transformation 81 . Using NEGF formalism, we have created various numerical roughness profiles and evaluated the graphene sheet's electronic transport properties. We will discuss the results of 10x10 graphene supercell Fig. 4a-b, correspond to the in-plane deformation roughness device structure with roughness varying from 5 pm to 10 pm, mimicking the in-plane deformation profile of h-BN and SiO 2 substrate. Similarly, Fig. 4i-j represents the corresponding density of state, Fig. 4c From Fig. 4, we have observed that an increase in the in-plane deformation from Fig. 4i 5 pm to Fig. 4j 10 pm, the density of state fluctuate more from the ideal flat structure of Fig. 1f, however, in comparison to corresponding  www.nature.com/scientificreports/ surface roughness structures distortion in the density of the state is more diminutive. Also, from Fig. 4c-d, we have observed that an increase in the in-plane deformation will mix up more energy modes progressively by reducing the energy plateau and reduce the average transmission probability by providing an extra state to scatter. Furthermore, from Fig. 4e-l we have observed that increasing in-plane deformation will fluctuate electron-hole density in the sample. The correlation length of deformation will also change the localization length scale to mean free path for transport. Furthermore, depending upon the correlation length, more modes can mix up in a specific energy window of transport 82 . Also, the modest increase in deformation amplitude will introduce these localized states, and transmission current through these new states will progress 83 . However, a radical increase in the in-plane deformation will effectively decrease the sheet width available for propagating modes, increase the transport bandgap, and decrease the current value. One of the main contributions of roughness scattering through non-equilibrium Green's functions formalism (NEGF) transport is that we have shown the corrugations scatterers is a long-range scattering potential, as depicted in Figs. 2g-j and 3g-i on 25nm X 35nm sheet with varying corrugations from 5 pm to 20 pm. Furthermore, with the increasing strength of corrugations, the scatterer's potential heavily disturbs graphene's ideal electronic transport properties. In out-of-plane roughness due to the fluctuating scattering potential distribution, it creates electron-hole puddles in the sheet; further inplane roughness is not so dominant to perturb the ideal graphene's transport as discussed above. Furthermore, to overcome these effects, graphene's conductivity increases by intentional doping, but it will lower mobility. Also, conductivity is enhanced by using multiple graphene layers, but it will trade-off with the lower optical transmission in the device, restricting its application as the transparent electrode. In the in-plane deformation roughness effect, we have observed that only the carbon-carbon bond length is modified. However, in the surface roughness effect, both carbon-carbon bond length and bond angle are simultaneously randomized due to corrugation in the substrate. Therefore, the electronic band structure and other transport properties are more sensitive to these surface fluctuations and deviate further from the idealized flat graphene sheet.

Impurities effect.
We have also investigated the effect of unwanted impurities on the electronic properties of the graphene sheet to elucidate the non-idealities that arrive during the fabrication and transfer step in a natural CVD-grown graphene sheet. Nitrogen and Phosphorus impurity atoms numerically substitute on the graphene sheet with corresponding bond radius and excess charge of the impurity atom. In the realistic scenario, numerous methods may arise in spatial charge distribution fluctuation in the graphene sheet. For example, an inhomogeneous electrical field applies at the gate or source/drain contact in the hall device through bias voltages is one such feasibility. Also, an inhomogeneous spatial distribution of impurity atoms in the sheet or below in the substrate or an inhomogeneous strain field due to inherent ripple of flake and topological morphology of underline substrate or overlapping encapsulation may introduce these effects. Ziegler et al. mathematically formulated the impurity scattering using linear response theory Kubo formalism. Furthermore, gauged robust minimal conductivity independent of disorder of sample in the Boltzmann conductivity limit by applying meanfield approximation to the Kubo conductivity tensor in a weak scattering regime exploiting the liner Dirac electronic bandstructure of graphene 84 . However, our simulation work is based upon the Landauer formula in the non-equilibrium Green's functions formalism, which is more generic to handle linear as well as non-equilibrium transport regimes and derived through Keldysh's generalized perturbation theory for the non-equilibrium system. In the non-equilibrium Green's functions formalism, the generic interacting current ansatz has been given by Caroli et al. [85][86][87] from the fluctuating electronic density between the interacting contact, and it corresponds to equation (5) of the Landauer formula for the current through an interacting electron region 88 . We have numerically incorporated the spatial distribution of impurities atoms in graphene by solving the Poisson equation of spatially fluctuating impurities potential with NEGF formalism. Henceforward, we calculate the spatial charge distribution in terms of density of state, electronic and hole density, and variation in transmitting density of modes. Furthermore, we have investigated the change in electronic properties with the interpretation of impurity density variation. We have varied the number of carbon atoms involving a substitutional impurity atom to consider these effects, effectively changing the device's doping density. Furthermore, We have changed the impurity atom type, e.g., Nitrogen and Phosphorous atoms, with varying atomic radius. Besides, in the Phosphorous impurity study, we intentionally change the electronic valency of the ion with ten electronic charges. This fictitious super ion configuration is simulated to mimic a larger diameter atomic or molecular defect with the high electronic charge concentrated at a particular lattice site on graphene sheet in the substitutional or adsorbent configuration pattern during the fabrication or transfer step. To circumvent the high computational cost of non-local scattering self-energy term in the NEGF framework. Adsorbed or substrate charged impurities trap are modeled in the graphene sheet as charge impurity center and truncated at the U 0 screened Coulomb potential energy superimposed upon the potential energy profile of impurity-free graphene sheet. Furthermore, the potential profile outside the screening length of impurity charges does not vary significantly from the impurityfree device configuration as, where r is the position of the current atom, �(r) is the screened Coulomb potential, Q is element charge, and r 0 is the position of the superimposed impurity atom, ε is the dielectric constant, U 0 onsite core potential modified to adjust impurity atom state energy, and D Debye screening length, www.nature.com/scientificreports/ where T is temperature, k is the Boltzmann constant, and N d is the impurity doping concentration. The charge placement can be done randomly and explicitly at a particular atomic location with a random dopant distribution, sampling, and determining a stochastic model. With the added computation overhead, charge impurities can also be added as a scattering self-energy similar to an explicit way without requiring a sampling step in the NEGF coupled Eq. (7) of motion. However, the convergence achieves fast in the self-consistent born calculation with Poisson's equation in a stochastic scheme. As the electrons and holes occupy low energy bands, long-range impurity potential interacts more critically with the device's long-range externally applied electrostatic potential. Therefore from a theoretical perspective in the NEGF framework, a more rigorous model with multiple impurity scattering, impurity-impurity interaction, and interaction with the device's electrostatic potential should be incorporated for better accuracy and added computational cost. Previously, Chen et al. demonstrated in-situ controlled potassium doping on pristine graphene devices in an ultrahigh vacuum chamber at a low T = 20 K temperature and measured gate-voltage-dependent conductivity, and investigated carrier density-dependent conductivity near Dirac-cone 89 . In our impurities, NEGF calculation similarly minimum conductivity value is observed. However, our calculation scenario is more complex. We calculate the transmission at room temperature and mix up all the propagating modes at thermal energy as it happens in the room temperature measurement in the CVD graphene sheet for any practical application. We have calculated the Density of state from Eqs. (8) and (19), electronic density and hole density from Eqs. (10) and (19), density of propagating modes M(E) in x-, y-direction from Eqs. (13) and (19), Electronic band structure from Eqs. (7) and (19), Spatially resolved electronic state and self-consistent Poisson potential from Eqs. (19) and (21) on a graphene sheet with Nitrogen and Phosphorus atom substitute impurity in the Figs. 5, 6 and 8, 9 respectively for the various supercell size to investigate the effect of impurity density by varying the supercell size. As observed in the CVD-grown graphene samples, the doping density varied batch to batch based on the recipe. To investigate the correlation between the level of impurity effect on the relative size of the graphene sheet, we have used a single impurity with varying sizes of graphene structure. Next in the figure Fig. 5, Fig. 5a correspond to bigger 20x20 graphene supercell flat structure with substitutional Nitrogen impurity, fig. 5b corresponding density of state, Fig. 5c electronic band structure, Fig. 5d density of modes, fig. 5e self-consistent Poisson potential due to Nitrogen impurity, Fig. 5f electronic density, Fig. 5g hole density. In the article's supplementary information, from Fig. 5a-g is also rendered in the full-text page resolution in the corresponding figures of Figs. 93-99 for the reader's reference. Also, spatially resolved electronic orbital wave-function amplitude |ψ| 2 for first seven eigenvalues of stationary solution of Schrödinger equation from |ψ 0 | 2 to |ψ 6 | 2 in the corresponding graphene supercell is represented in the article's supplementary information from Supplementary Figs. 100-106. In the simulated device, the primitive unit cell has two atoms per cell, and a total of 800 atoms are simulated with a domain size of 3200 finite element mesh grid points. The total degree of freedom in hamiltonian is 2400 variable-sized orbitals. In Fig. 5, from Fig. 5b, we have observed that Nitrogen impurity will induce additional state and density of state fluctuation and hence provide more scattering probability to propagating electron in a graphene sheet. Also, from Fig. 5c there is a small energy bandgap opening in the band structure and change in dispersion curvature of Dirac fermion and a dip in Fermi level due to additional Nitrogen electronic charge in sheet. Furthermore, from Fig. 5d there is a more density of modes mixing up in comparison to idealistic characteristics as distinct energy level plateaus are vanishing. Moreover, from Fig. 5f-g electronic and hole density significantly fluctuate from the constant neutral density due to impurity induced self-consistent Poisson potential of 0.3 eV at the impurity site and spread in the vicinity of next to the next nearest neighbor in the lattice as observed in Fig. 5e. Next in the figure Fig. 6, Fig. 6a correspond to biggest 50x50 graphene supercell flat structure with substitutional Nitrogen impurity, Fig. 6b corresponding density of state, Fig. 6c electronic band structure, Fig. 6d density of modes, Fig. 6e self-consistent Poisson potential due to Nitrogen impurity, fig. 6f electronic density, Fig. 6g hole density. In the article's supplementary information, from Fig. 6a to Fg. 6g is also rendered in the full-text page resolution in the corresponding figures of Figs. 107-113 for the reader's reference. Also, spatially resolved electronic orbital wave-function amplitude |ψ| 2 for first seven eigenvalues of stationary solution of Schrödinger equation from |ψ 0 | 2 to |ψ 6 | 2 in the corresponding graphene supercell is represented in the article's supplementary information from Supplementary Fig. 114-120. In the simulated device, the primitive unit cell has two atoms per cell, and a total of 5000 atoms are simulated with a domain size of 20000 finite element mesh grid points. The total degree of freedom in hamiltonian is 15000 variable-sized orbitals.
In Fig. 6, from Fig. 6b, we have also observed that Nitrogen impurity will induce additional state and density of state fluctuation and hence provide more scattering probability to propagating electron in a graphene sheet; however, the severity of distortion is less in comparison to the smaller geometry case as discussed previously. Also, from Fig. 6c there is a restoration of the Dirac fermion band structure, and from Fig. 6d However, due to additional availability of state, there is a still density of modes mixing up in comparison to idealistic characteristics as distinct energy level plateau are fading. However, distortion is less in comparison to the smaller geometry case of Fig. 5. Moreover, from Fig. 6f-g electronic and hole density significantly fluctuate from the constant neutral density due to impurity induced self-consistent Poisson potential of 0.3 eV at the impurity site and spread in the vicinity of next to the next nearest neighbor in the lattice as observed in surface projected Poisson potential landscape plot of Fig. 6e. Also, in comparison to Fig. 5 from Fig. 6, we have observed that a single impurity state in the larger device does not affect the energy band edge significantly on the energy scale. Also, it induces insignificant energy broadening compared to contact induces energy broadening as the system is in a non-equilibrium open quantum state. The weak perturbation effect of a single impurity state validates these estimates. Therefore, the device's electrostatic potential landscape did not further calibrate for this new impurity www.nature.com/scientificreports/ state induces a potential loop, and overall propagating states are unimpaired by this variation. However, such a numerical truncation situation rapidly changes for a smaller device geometry with a single impurity state or a larger device with multiple impurity states. Therefore, a more detailed self-consistent calculation incorporating an impurity potential loop is required to reach higher accuracy. Next, we have doped the corrugated graphene sheet of 100x100 supercell with substitutional Nitrogen impurities with varying corrugation from 5 pm to 20 pm in the Fig. 7 Fig. 7, we observed that the transmission profile and atomistic electronic density were disturbed significantly by the dual effect from the flat idealistic and uniform charge distribution graphene sheet and the formation of electronic-hole puddle formation in the graphene flake. In the article's supplementary information, from Fig. 7b-i is also rendered in the full-text page resolution in the corresponding figures of Figs. 121-132 for the reader's reference. In the Fig. 7 with the progressively increasing corrugation profile with impurities atoms doping in the large 100x100 graphene supercell from corrugate-5pm Fig. 7b to 20pm Fig. 7e, We have observed the electronic density of transmission modes M(E) in x-y direction heavily distorted from the idealistic condition. Moreover, the electronic density for corrugate-5pm Fig. 7f to 20pm Fig. 7i, also fluctuated, as also seen in Fig. 7j view. The impurity atoms, inhomogeneous impinging applied electrical field, and various strain fields through ripples, corrugation in the graphene sheet increase the energy mixing of propagating density of modes. Furthermore, the extra state density assists in propagating electrons to scatter more, therefore reducing overall conductance. Nevertheless, at the same time, even the slight variation in band structure dispersion gives rise to variation in the effective mass of Dirac fermion from idealistic values, and overall observed mobility and conductivity deteriorate furthermore from the flat characteristic. We have observed the slight kink in the density of state around zero-energy Dirac point due to substitutional impurity and the start of band branch splitting around the high symmetric points of the Brillouin zone, which broadens the original graphene band at these points. These additional band branches will influence the energy landscape in the next few nearest-neighbor carbon atoms. Hence, more accuracy can be achieved in band structure calculation by adopting a higher-order tight-binding model for an impurity-driven sheet in an open boundary device or using an electrically close boundary sheet in the ab-initio calculation with few ten of thousand atoms device with additional computation overhead. Next, we have a doped graphene sheet with substitutional Phosphorus impurity. Although intentionally replace the electronic valency of the ion with ten electronic charges to create a fictitious super ion to mimic a larger diameter atomic or molecular defect with the high electric charge concentrated at a specific lattice site on the graphene sheet. In the figure Fig. 8, Fig. 8a correspond to 20 x 20 graphene supercell flat structure with substitutional Phosphorus impurity, Fig. 8b corresponding density of state, Fig. 8c electronic band structure, Fig. 8d  In Fig. 8, from Fig. 8b, we have observed that Phosphorus impurity with excessive charge will induce extra state density and hence provide more scattering probability to propagating electrons in a 20x20 graphene sheet. The density of state shifted significantly towards the left due to shifting in Fermi value due to excessive electronic impurity charge doped in this sample; a similar trend is visible in the band structure plot from Fig. 8c. There is a tiny band opening in the band structure, and the dispersion of graphene has heavily deteriorated from Dirac fermion. From Fig. 8d there is a more density of modes mixing up in comparison to idealistic characteristics as distinct energy level plateau dissolved. Moreover, from Fig. 8f electronic density significantly fluctuates from the constant neutral density due to impurity induced self-consistent Poisson potential of 0.85 eV at the impurity site and spread in the vicinity of next to the next nearest neighbor in the lattice as observed in surface projected Poisson potential landscape plot of Fig. 8e. The increase in the Poisson potential in this case as compared to Figs. 5e-6e is due to the excess charge we intentionally doped at a lattice site to investigate the effect of a heavily contaminated site in a CVD Graphene sheet. Next, in the Fig. 9, Fig. 9a correspond to 50 x 50 graphene supercell flat structure with substitutional Phosphorus impurity, Fig. 9b corresponding density of state, Fig. 9c electronic band structure, Fig. 9d density of modes, fig. 9e self-consistent Poisson potential due to Phosphorus impurity, Fig. 9f electronic density. In the article's supplementary information, from Fig. 9a-f is also rendered in the full-text page resolution in the corresponding figures of Figs. 146-151 for the reader's reference. Also, spatially resolved electronic orbital wave-function amplitude |ψ| 2 for first seven eigenvalues of stationary solution of Schrödinger equation from |ψ 0 | 2 to |ψ 6 | 2 in the corresponding graphene supercell is represented in the article's supplementary information Supplementary Fig. 152-158.
In Fig. 9, from Fig. 9b, we have observed that Phosphorus impurity with excessive charge will induce extra state density kink and hence provide more scattering probability to propagating electron in a 50x50 graphene sheet. However, again in comparison to Fig. 8b fluctuation in density of state of becoming smooth. The density of state shifted towards the left due to shifting in Fermi value because of excessive electronic impurity charge in the sample. In the band structure plot from Fig. 9c. There is a mixing of states, and the dispersion of graphene has heavily deteriorated. From Fig. 9d there is a more density of modes mixing up in comparison to idealistic characteristics as distinct energy level plateau dissolved. However, the situation improves from the smaller sample calculation with the same number of impurity charges. Moreover, from Fig. 9f electronic density significantly fluctuates from the constant neutral density due to impurity induced self-consistent Poisson potential as observed in Fig. 9e,  www.nature.com/scientificreports/ in Fig. 8e. The phosphorus atom, due to its relative atomic size, bond radius, excess orbitals, and additional states compared to the carbon atom, the impurity effect spread in the graphene lattice. Also, these excessive impurity state influences the energy landscape in a couple of next nearest-neighbor carbon atoms. These excess impurity charges are included in numerical simulation by replacing the carbon atom with Nitrogen and Phosphorus atom with an adequate charge. From the above observation, we conclude that the enhancement of conductivity of graphene sheet through the substitutional and adsorbent impurity doping route will come with its tradeoff cost of the enhancement of additional density of states, deterioration of Dirac band dispersion, and fluctuation of transmission and modes density. Therefore reduce mobility values in the sheet. On the other hand, enhancement of conductivity of graphene by depositing multiple sheets in the device will reduce the transparency of the sheet and its possible application in display devices. Moreover, the enhancement of conductivity of graphene sheet through indirect induce doping by gate field will trade off with Fermi level shift and shift in conduction from valance band branch to conduction band branch. Furthermore, giving rise to ambipolar characteristics, and turning off such device by gate filed, especially in transistor application, will be difficult. The combined impact of corrugations and impurities may arise in practical fabrication and transfer processes. The practical impact of various nanoscale atomistic scattering mechanisms in the semiconductor consider by Matthiessen's rule and total scattering length and time is the summation of all the underline individual effects 53,90-92 , Furthermore, the most prominent mechanisms of these various physical effects determine the final device characteristic. In the NEGF loop, these various mechanisms are adopted as Dyson's self-energy loop to solve self-consistently. The macroscopic sample's characterization properties are defined by the microscopic averaging out of competing scattering statistics rather than individual scattering events of the order of the femtosecond process. Finally, we conclude the article with a summary of the results and observations.

Summary
We have numerically investigated various random, non-idealities, e.g., inhomogeneous out-of-plane surface corrugation, in-plane deformation, and excess atomistic charge impurities in the 2-dimensional graphene sheet. These non-idealities primarily arise in the Roll-to-Roll chemical vapor deposition (CVD) and plasma-enhanced chemical vapor deposition (PECVD) manufacture process. Mitigating these variances is necessary for the practical application success of the graphene industry, and the high fidelity and uniformity require in the graphene production, transfer, benchmark process is the utmost requirement for developing the graphene-based industry.
We have employed a multiscale, multi-physics-based non-equilibrium Green's function framework in a third nearest-neighbor tight-binding configuration. The Tight-binding model is computationally efficient to incorporate atomistic effects on mesoscopic electrical properties. Our detailed study provides an essential understanding for evaluating and benchmarking the electronic properties of graphene sheets. It also oversees the nanoelectronics device design aspect on graphene and 2-D material for practical application in industry. Our calculation has not explicitly treated surface roughness scattering as a self-energy term in the NEGF framework, which is numerically cumbersome, affecting the device's performance. In future investigations, we will attempt to incorporate that in an equal mathematical footing in the framework. Also, a more sophisticated local density of state-dependent scattering model should be employed for a more detailed investigation of the highly corrugated device. Also, surface passivation treatment can influence the device's electrical characteristics. Theoretical modeling of capping material effect, doping efficiency, long-range electrostatic screening, interface effect, interlayer phenomena, contact quality, and heat transfer in the low-dimensional material-based device should be examined further for quantitative evaluation. However, scattering is omitted and treated as residue in the most published articles in the 2-D materials domain and ignored without sufficient justification. In addition, the atomistic resolutions are critical for the charge in the self-consist computations. The theoretical treatment of defect situations is even challenging at the multiscale simulation. However, these are challenging transport situations at the nanoscale's theoretical modeling in the graphene and 2-D material domain. www.nature.com/scientificreports/ Correspondence and requests for materials should be addressed to B.B.
Reprints and permissions information is available at www.nature.com/reprints.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.