Electronic transport in graphene with out-of-plane disorder

Real-world samples of graphene often exhibit various types of out-of-plane disorder -- ripples, wrinkles and folds -- introduced at the stage of growth and transfer processes. These complex out-of-plane defects resulting from the interplay between self-adhesion of graphene and its bending rigidity inevitably lead to the scattering of charge carriers thus affecting the electronic transport properties of graphene. We address the ballistic charge-carrier transmission across the models of out-of-plane defects using tight-binding and density functional calculations while fully taking into account lattice relaxation effects. The observed transmission oscillations in commensurate graphene wrinkles are attributed to the interference between intra- and interlayer transport channels, while the incommensurate wrinkles show vanishing backscattering and retain the transport properties of flat graphene. The suppression of backscattering reveals the crucial role of lattice commensuration in the electronic transmission. Our results provide guidelines to controlling the transport properties of graphene in presence of this ubiquitous type of disorder.

Being the first and the most investigated twodimensional (2D) material, graphene continues attracting attention as a platform for exploring novel physics and realizing prospective technological applications [1].The 2D nature of graphene gives rise to soft flexural modes that result in low-energy out-of-plane disorder otherwise absent in bulk, three-dimensional materials [2][3][4][5].The interplay between bending upon in-plane compression and the interlayer adhesion results in several distinct types of out-of-plane disorder: ripples, wrinkles and folds (see Refs. [2,6] and Figs.1(a,b)).The out-of-plane disorder has a prominent effect on the electronic structure and transport properties of graphene [7][8][9][10].Finite curvature of the deformed region results in pseudo-gauge fields [11,12], while the collapsed regions in wrinkles and folds provide a pathway for electronic tunnelling between layers [6,13].In addition, out-of-plane disorder locally accumulates charges and act as scattering centers [6,[14][15][16], subsequently having an impact on the operation of graphene-based nanoscale electronic devices [13,17,18] as well as electrical characteristics of large-scale graphene samples.
Out-of-plane disorder in graphene may occur for several reasons.For instance, graphene grown using the chemical vapour decomposition (CVD) process develops wrinkles and folds as a result of the thermal contraction of substrate during the cooling stage [19][20][21].The out-of-plane disorder may also be introduced during the transfer procedure [22,23].Significant efforts have then be devoted to eliminating wrinkles [20,24], e.g. using the substrates with matching thermal expansion coefficients [22], strain engineering [25] and tailored temperature control protocols [20].Experimental studies of the electronic transport in graphene with out-of-plane disorder have also been published [6,26].It was proposed that controlled folding of graphene can be used for engineering charge-carrier dynamics [27][28][29][30].No question, future applications of graphene in electronics call for a detailed understanding of the effect of this ubiquitous type of disorder on the electronic transport.
In this work, we systematically investigate the electronic transport across wrinkles and folds in graphene using first-principle computations.For commensurate graphene wrinkles, in which the interlayer stacking corresponds to the energetically favorable Bernal stacking configuration, we find that the electronic transmission oscillates over wide energy ranges.The observed oscillation patterns are attributed to quantum interference between the inter-and intralayer transport channels.In incommensurate wrinkles and folds, the mismatch between the layers is found to suppresses the interlayer tunneling resulting in transmission probabilities close to the limit of flat, pristine graphene.

A. Construction of models
The atomistic models of graphene with out-of-plane disorder considered in our work are defined by a compressive displacement of length ∆ W (see Fig. 1(a)) forming a wrinkle or a fold along crystallographic vector v = (a, b).The considered configurations are thus assumed to be periodic along v.The interplay between the bending energy and attractive interlayer interactions of graphene layers define the evolution across the three types of outof-plane disorder realized upon increasing ∆ W as shown in Fig. 1(b).While ripples are formed at small ∆ W , interlayer attraction collapses such structures to wrinkles for larger values of ∆ W , and further increase of ∆ W leads to folds, in which the contact area between graphene layers is further increased.Extremities of wrinkles and folds have loop-like structures free of interlayer coupling [6].All atomistic models of wrinkles and folds considered in our work have been constructed with the help of classical force-field relaxation (see the Methods section for details).

B. Electronic transport across commensurate wrinkles
We first consider the special case of wrinkles defined by v = (1, 0) and v = (1, 1), referring to them as zigzag and armchair, respectively.The collapsed regions of such wrinkles are compatible with the energetically favorable Bernal interlayer stacking configuration [31][32][33][34], and hence referred to as commensurate in the rest of our paper.For these relaxed models, we calculated ballistic charge-carrier transmission from first principles, using the combination of density functional theory (DFT) and the non-equilibrium Green's function formalism implemented in the TranSIESTA package [35,36] (see Methods).The results of DFT calculations are discussed in comparison with the tight-binding (TB) approximation calculations employing the Slater-Koster formalism [6,37] (see the Supplementary Materials).Figures 2(a  E and momentum parallel to the wrinkle k // .Furthermore, each panel shows transmission T (E) plotted at a specific k // =±2π/(3a 0 ) (a 0 = 2.46 Å is the lattice constant of graphene), which corresponds to the momentum of projections of the Dirac cone band degeneracies.
There are two striking observations in the presented transmission plots.Firstly, both in DFT and TB results, we observe a pronounced electron-hole asymmetry in the charge-carrier transmission.The electron-hole asymmetry has an origin in the interlayer stacking of zigzag wrinkles.The collapsed region assumes Bernal stacking configurations AB 1 or AB 2 [38], as illustrated in Fig. 1(c), in which one of the graphene sublattices couples to itself upon folding since the two layers are mirror-symmetric with respect to each other.Such a coupling breaks the sublattice symmetry and hence the electron-hole symmetry [39,40].
Secondly, ballistic transmission T (E, k // ) shows pronounced oscillations over broad energy ranges.Apart from making transmission highly energy-dependent, such oscillations also affect average conductance at a finite bias.These oscillations are clearly visible in the side panels of Figs.2(a-d) that show transmission at a fixed momentum k // = 2π/(3a 0 ) that corresponds to the projections of the Dirac points.Further analysis shows that the energy separation ∆ E between the peaks has an approximately linear dependence on ∆ W (Fig. 2(e)).Such a dependence is the signature of the interference between the interlayer and intralayer transport channels, as found by some of us previously in the case of electromechanical response of bilayer graphene [13].This transport phenomenon is further addressed in Section I C.
The second family of investigated commensurate configurations is defined by v = (1, 1), that is wrinkles are oriented along the armchair direction.Atomic relaxation effects are more complex in such wrinkles.Unlike in the zigzag case, realizing the lowest-energy Bernal stacking is possible only at a cost of introducing shear deformation as shown in Fig. 3(a).Consequently, the Bernal stacking is not achieved at small values of ∆ W , and the collapsed region assumes the saddle-point (SP) stacking configuration [41] that does not break sublattice symmetry.Figure 3 c-d) present the transmission maps for the armchair wrinkles with ∆ W = 40 Å and ∆ W = 120 Å.In the case of v = (1, 1), the Dirac points are projected onto k // = 0. Similar to the case of zigzag wrinkles, oscillations with the ∆ E ∝ 1/∆ W period are observed in the transmission maps.The oscillation pattern is more regular than in the case of ∆ W = 40 Å armchair wrinkle, which assumes the SP stacking and hence preserves electron-hole symmetry.In contrast, the ∆ W = 120 Å wrinkle is significantly closer to the Bernal stacking (see Fig. 3(b)) and the electron-hole symmetry appears to be well visible in this case.

C. Conductance oscillations in the atomic chain model
In order to further address the physical mechanism underlying the conductance oscillations observed in both the zigzag and armchair wrinkles, we introduce a simple one-dimensional model treated using the tight-binding approximation.The presence of interlayer conductance channels is defined by ∆ W , and also l that represents the absence of interlayer hopping in the loop-like region as shown in Fig. 4(a).At the same time, we observe that k // does not have any significant effect on the oscillation period, hence we introduce a one-dimensional chain described using the nearest-neighbor tight-binding model with an extra hopping t that models interlayer coupling in graphene wrinkles.Schematic diagram of this model with hopping t represented by a rainbowlike graph is shown in Fig. 4(b).The ratio of the newly introduced hopping t to the nearest-neighbor hopping t is chosen to resemble that of graphene wrinkles t /t = 0.48 eV/−2.7 eV [6,33].mission T as a function of energy E at a fixed ∆ W =12 in units of intersite distance, while parameter l is varied.We observe that oscillation peaks have the same positions, which indicates that l is of little effect on the oscillation period.Combined with the results of DFT calculations we conclude that the oscillations are defined by the largest path difference ∆ W .We further analyze the transmission oscillations in the atomic chain model using the non-equilibrium Green's functions (NEGF) approach, in which hoppings t are treated as a perturbative correction to the transmission.
First, we define an infinite atomic chain with the Hamiltonian where c i (c † i ) is the annihilation (creation) operator on the ith site.This Hamiltonian commutes with the translation operator, thus the energy eigenstates are also momentum eigenstates.
In the NEGF formalism [42,43], the transmission is calculated as where G is the Green's function The coupling matrices Γ i are given by Γ i =i(Σ i − Σ † i ), with Σ i being the self-energies of the two semi-infinite leads.
Green's function G 0 describes the chain in absence of t , while adding coupling t that models interlayer coupling in wrinkles adds an additional term ∆h ∆h = t The Green's function is then Keeping only the first order of correction G 0 ∆hG 0 , the transmission becomes The Green's function can be written as an expansion involving eigenstates |ψ n of the chain with no hoppings t and the correction term G 0 ∆hG 0 becomes As the simplest case, we analyze the E n = E m correction G 0 ∆hG 0 = ψ n |∆h|ψ m G 0 that gives an E i -dependent prefactor to the Green's function.We write the factor as a function δ(E) as The leading order of transmission correction is Γ hence the correction to transmission contains δ 2 + 4δ + 1.
We then evaluate the correction δ(E), keeping in mind that the eigenstates of the pristine chain are also momentum eigenstates.The correction factor δ represents the phase difference between wavefunctions: connected by the additional hoppings t .It can then be approximated by a sum of sinusoidal functions The results of the summation shown in Fig. 4(d) suggests that the highest-frequency component in Eq. (11), which corresponds to the interference path ∆ W , defines the oscillation peaks.Our first-principles results are consistent with the conclusions of this simple model.

D. Transport across incommensurate wrinkles
We will now discuss graphene wrinkles formed along general crystallographic directions v = (a, b) other than high-symmetry zigzag and armchair orientations.In these cases, the collapsed region locally forms twisted bilayer graphene with matching vectors (a, b) and (b, a).The resulting twist angle is while the translational vector along the wrinkle has a length of d = √ a 2 + b 2 + ab.
We discuss the effect of wrinkle direction (a, b) on the transmission T (E, k // ).Translational vector (a, b) defines a one-dimensional mini Brillouin zone (mBZ) obtained by projecting the 2D Brillouin zone of graphene onto the k // direction in momentum space.The Dirac cones of graphene are projected onto either k // = 0 (class Ia) or k // = 2π/(3|v|) (class Ib) of the mBZ according to the classification introduced in Ref. [44].Class Ia is defined by |a − b| mod 3 = 0, class Ib otherwise.The projections of the Dirac cones define the regions in the T (E, k // ) maps where transmission is allowed and limited by n conductance channels in case of n-fold degeneracy of bands at given E and k // in the ballistic regime.
The periodic structure of wrinkles results in consequences deeper than just the conservation of momentum k // upon ballistic transmission.We stress that semi-infinite graphene sheets on both sides of wrinkles of constant width have the same crystallographic orientation.The momentum conservation implies suppressed backscattering at the Dirac point, which can be observed by evaluating contribution to the transmission from the first-order correction G 0 ∆hG 0 .Starting with the pristine graphene and a simple interlayer containing only hopping between aligned atoms the effective ∆G writes which becomes most significant at E m = E n = z.Recalling the fact that |ψ m and |ψ n are eigenstates of pristine graphene, ψ m |∆h|ψ n gives an exp(2πi Integrating over r ij , vanishes if k m = k n , while the wrinkle enforces a transformation k m = M x k n due to its mirror-symmetric stacking configuration of the two layers as shown in Figs.5(a,b).Here, M x denotes the mirror-reflection with respect to transport direction x: M x (k x , k y ) = (−k x , k y ).From the above rules of momentum conservation, we conclude that the transmission is only affected in the overlapping region of the Dirac cones.In the non-overlapping region, the correction G 0 ∆hG 0 is vanishing, and the transmission retains the value of ideal, defect-free graphene.These results are verified by the explicit DFT transport calculations as shown in Fig. 5(d-e) for class Ia and class Ib wrinkles, respectively.The transmission maps T (E, k // ) have overall shape of the Dirac cone projections.Transmission values near the charge neutrality are T ≈ 2 and T ≈ 1 for class Ia and Ib configurations, respectively, indicating that interlayer tunnelling plays a minor role.At higher energies where the Dirac cones overlap, e.g.near E ≈ 2 eV in Fig. 5(e), backscattering becomes significant leading to a series of transmission dips.We also point out that class Ia presents larger backscattering from the interlayer coupling since the projected Dirac cones overlap with each other.

E. Transport across graphene folds
We will now discuss folds as the ultimate regime of out-of-plane disorder in graphene.Folds realize triplelayer graphene configurations in their collapsed regions (Fig. 6(a-c)).Importantly, adjacent layers (pairs 1-2 and 2-3) in incommensurate folds are twisted with respect to each other, while the outside layers 1 and 3 are aligned.This configuration is equivalent to mirrorsymmetric twisted trilayer graphene.While we still expect the effect of interlayer coupling to be weakened by the incommensuration, our DFT calculations predict a larger degree of backscattering in folds than in wrinkles (compare Figs. 5(e) and 6(d) for the the (1,2) direction).For the folded region of width l f = 40 Å, the average transmission in the energy interval (−0.15 eV, 0.15 eV) is 0.727, while in the wrinkle of equivalent ∆ W = 80 Å it is 0.908.The observed transport behaviour raises the question of whether the enhanced backscattering in incommensurate folds as compared to wrinkles originates from the direct coupling of the outmost layers 1 and 3.The corresponding matrix elements of the Hamiltonian in localized-basis-set first-principles calculations [35,45], are found to be negligible.The estimated Slater-Koster coupling also has a negligible magnitude of 10 −4 eV.Therefore, we attribute the enhanced scattering to the fact that the number of interlayer tunneling channels is doubled in the folds.As expected, for a commensurate zigzag fold (Fig. 6(d)) we observe strong backscattering with transmission magnitudes lower than in the equivalent zigzag wrinkles (Fig. 2).

II. DISCUSSION
We investigated the effect of our-of-plane disorder on the electronic transmission in graphene.Different forms of the our-of-plane disorder exist in graphene, depending on the compressive displacement and the orientation of the deformation.Our work studied ballistic transmission through the wrinkles and folds using first-principle calculations, taking into account their width and interlayer commensuration.
The interlayer coupling was found to cause substantial oscillations the electronic transmission across commensurate wrinkles.Such oscillations were found to originate from the quantum interference involving the interlayer tunneling channels.Based on DFT calculations, we propose a simple one-dimensional model that fully captures the observed oscillations.On the other hand, in incommensurate, "twisted" wrinkles the interlayer coupling is effectively weaker, and the transmission near the Fermi level preserves that of pristine, flat graphene.We have also found enhanced backscattering in folds that was attributed to the doubled contact region in this type of the out-of-plane disorder.
Our results offer an approach toward understanding the transport in mesoscopic graphene samples containing out-of-plane disorder of different type and arbitrary orientation.The theory of transmission across graphene wrinkles and folds is thus useful for designing graphenebased devices as well as fold-engineering of graphene.As a generalization, the principles presented in our work are expected to apply also to other types of 2D materials.Formation of locally twisted bilayers in the wrinkles and folds provides an interesting outlook for further studies, e.g. the "twisted" wrinkles in the smaller-angle regime.

A. Structure relaxation with classical force fields
The atomic structures of models of the out-of-plane disorder in graphene were obtained by means of the classical force field simulations using LAMMPS [46,47].The classical force field includes the bond-order potential for describing covalent bonding [48] as well as the modified version of the Kolmogorov-Crespi registry-dependent potential [49] for describing the interlayer van der Waals interactions.The energy minimization was performed using the conjugate-gradient and fire algorithms.

B. The tight-binding model calculations
In order to describe both the interlayer coupling and the effect of curvature in the tight-binding calculations of graphene with out-of-plane disorder, we employ the Slater-Koster model [6,37].The p z atomic orbitals of carbon atoms form the intralayer π bonds and the interlayer σ bonds.The general form of the Hamiltonian including both contributions is Explicit expressions for the hoppings t ij π and t ij σ are [6] Following the previous [6] Slater-Koster parametrization, we set V 0 π = −2.7 eV, V 0 σ = 0.48 eV, characteristic distances a 0 = 1.42 Å, d 0 = 3.35 Å and the decay length r 0 = 0.184a (a = √ 3a 0 = 2.46 Å as defined in the main text).In the orientation-dependent terms, angles θ i and θ j are defined as the angle between r ij and the local normal vector at atomic positions r i and r j , that is θ i = ∠(r ij , n i ).These terms accounts for the effect that the local curvature of graphene sheet on the overlap between p z orbitals.

C. Recursive Green's function methods
The ballistic transmission was calculated using the non-equilibrium Green's function methods in both the TB model and DFT calculations.The transmission probability is expressed as where G is the Green's function as used in Eq. ( 2): G = [G 0 −Σ] −1 .The Γ matrices contains the self-energy terms of the two leads The self-energy from the ith lead is calculated as Σ i = h i G i h † i , where h i is the coupling matrix between the lead and the scattering region.For each of the semi-infinite leads, Green's function G i is obtained through the recursive Green's function methods.In each step one layer is added to the lead, and the Green's function iterates as g j =[E − h − T g j−1 T † ] −1 .G i is taken as the converged value of g, that is G i = g j→∞ j .

D. First-principles electronic transport calculations
First-principles transport calculations were performed with TranSIESTA package [35,36].We used the double-ζ plus polarization basis set combined with the local density approximation exchange-correlation functional [50].The energy shift for constructing the localized basis was set to 275 meV, and the real-space cutoff to 250 Ry.The estimation of the direct coupling between the top and bottom layers in graphene folds was extracted from the localized basis set Hamiltonian using the sisl package [45].

FIG. 1 .
FIG. 1.The structure of out-of-plane disorder in graphene.(a) Definition of compressive displacement ∆W relative to the flat, unstrained graphene.(b) Formation of the three distinct types of out-of-plane disorder upon increasing ∆W .The curves show a schematic illustration of the dependence of energy E on ∆W for the three deformation regimes.Yellow color exposes the collapsed regions where the interlayer coupling is enabled.(c) Illustration of the interlayer coupling between the atoms belonging to the same sublattice in commensurate zigzag wrinkles and folds.
)-(d) present the ballistic transmission T (E, k // ) for the models of zigzag wrinkles defined by ∆ W = 40, 60, 120 and 240 Å as a function of energy
(b) presents the evolution of shear deformation ∆ y upon the change of ∆ W with ∆ y = a 0 /(2 √ 3) representing the pure Bernal stacking configuration.Figures 3( FIG. 3. (a) Schematic illustration of the shear deformation in armchair wrinkles.The shear is characterized by displacement ∆y.(b) Evolution of shear deformation ∆y versus compressive displacement ∆W .At small values of ∆W , shear deformation ∆y is small, which corresponds to to the SP stacking configuration (∆y = a0/(2 √ 3) corresponds to pure Bernal stacking.(c-d) Ballistic transmissions T (E, k // ) across armchair wrinkle models defined by (c) ∆W = 40 Å and (d) ∆W = 120 Å.The T (E) cross sections are taken at k // =0 that corresponding to the projected Dirac points.

5 FIG. 4 .
FIG. 4. Transmission oscillations in atomic chain model.(a) Cross-section drawing of the trivialized graphene wrinkle and (b) its unfolded representation equivalent to atomic chain with additional hoppings.(c) Transmission T as a function energy E in units of t calculated using the TB model Hamiltonian.In this plot ∆W = 12 in units of intersite distance is fixed, while different curves correspond to difference values of l.(d) First-order correction to the Green's function δ(E) = G0∆hG0/G0 plotted for different l and constant ∆W = 20 reveals that the period of oscillations is governed by ∆W .

FIG. 5 .
FIG. 5. Atomic structure of incommensurate wrinkle defined by the (1,2) direction: (a) local structure of the collapsed region equivalent to twisted bilayer graphene (unit cell is shown with the shaded region), (b) side-view with the sketch of the Brillouin zones and the Dirac cones of adjacent layers, and (c) top-view of the wrinkle illustrating the conservation of crystallographic orientation of graphene leads.Transmission maps T (E, k // ) for wrinkle models defined by (d) v=(1,4) and ∆W = 80 Å(class Ia), (e) v=(1,2) and ∆W = 80 Å(class Ib).

FIG. 6 .
FIG. 6.(a) Atomic structure of an incommensurate fold defined by v = (1,2) as an example.(b) Side-view of the fold with layers numbered and Brillouin zone orientations indicated.Transmission maps T (E, k // ) for (c) the incommensurate fold shown in the above panels and (d) zigzag fold characterized by ∆W = 80 Å.