{'fraction_non_alphanumeric': 0.027749563238938773, 'fraction_numerical': 0.0076585544224596755, 'mean_word_length': 4.516248319139399, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Eigenfrequencies of the randomly pinned drum and conductivity of graphene

Graphene is convenient material for nanomechanichal applications since high-frequency oscillations are easily accessible. In this Article, we consider graphene on a rough substrate attached to imperfections at random locations. We explore the statistics of low-lying phonon modes, which exert most influence on the conductivity of graphene. We find that the nearest neighbor spacings of low lying eigenfrequencies have the Wigner-Dyson probability distribution after averaging over the random configurations of disorder. Due to interaction of electrons with the oscillations of the membrane, an electron can be transfered to higher or lower energies, which is a manifestation of the phonon-assisted Tien-Gordon effect. The Tien-Gordon effect suppresses the conductivity of graphene. In the regime of low Fermi energies and small sizes of the sample an increase of conductivity is observed which we refer to Klein tunneling and electron-hole pair creation. Eventually, when the increase of the transmission becomes too prominent, the pair creation changes the ground state of the system, signalizing the limit of applicability of the single-particle Dirac equation used in this paper.


I. INTRODUCTION
Graphene is a monolayer of graphite with excellent elastic properties and low mass, which make it a prospective nanomechanical resonator 1 attractive for a wide range of applications from mass sensing in singlemolecule range 2 to quantum manipulation of the elementary mechanical vibrations 3 . They are combined with perfect electrical properties. The mobility of the carriers, both in suspended samples 4 and the samples on the substrate 5 , is one of the highest among all modern semiconductors. The combination of mechanical and electrical properties opens wide opportunities for manipulation and read-out of mechanical oscillations using electrical signals.
In most experiments, graphene is deposited on substrate. The substrate is often not ideal and can be considered as a random potential landscape. The graphene membrane interacts with the substrate, which may result to its attachment to the substrate at some areas or points, and free suspension between these areas. This arrangement was confirmed experimentally 6 . This means that the mechanical oscillation spectrum of such membrane is complex, with localized and extended phonons determined by the substrate profile. Moreover, since the profile is random, the frequencies of these excitations can be considered as random quantities. It is known that the frequencies do not uniquely characterize the shape of the oscillating membrane. Indeed, this question was studied in 1966 by Marc Kac in his famous paper "Can one hear the shape of the drum" 7 , with the conclusion that is possible to construct different in shape, but isospectral membranes. Therefore the statistical approach remains the only meaningful way to characterize the frequency spectrum of this system.
In this Article, we characterize statistical properties of these frequencies and show how they can be assessed via electric transport measurements.
The elastic energy of the oscillating membrane has two contributions: stretching and bending energies. We conclude that the bending energy is much smaller than the stretching one. Hence the dynamics of graphene membrane is described by the Helmholtz's equation with the fixed boundary conditions at the regions where the membrane is attached to the substrate. The same equations describes the oscillation of the classical membrane. In our work we show how to "hear" the disorder of attachments of the membrane and see the distribution of the low frequencies averaged over the disorder. The frequencies can be directly measured in the experiment, but for graphene it is more convenient to extract information on the oscillation of the membrane from the conductivity measurements.
The high-energy spectrum of the eigenvalues of the Laplace operator is well studied and is known to obey the Wigner-Dyson distribution 8 . This statement was proven by a variety of methods, analytical (random matrix theory and non-linear sigma-model) as well as numerical (various models of quantum billiards). Therefore in this Article we study the probability distribution of the low energy spectrum of the membrane on disordered substrate. This question was not previously addressed since most of the studies of eigenvalues concern electron systems, where the relevant energy levels lie around the Fermi surface and have very high energy. In contrast, the phonon frequencies have a clear cut-off at low energies. This cut-off frequency determines the scale for low-frequency modes. The second reason why low energy levels are beyond the scope of the usual methods is that no statistical analysis can be performed for only few levels. Thus, for a clean system the problem of the statistics of low-lying eigenvalues can not be formulated. One needs a disordered system to address this issue: Averaging over the disorder acts as averaging over the statistical ensemble. The problem can be addressed only numerically as the analytical methods mentioned above are not applicable for low energy spectrum. According to our numerical modeling the distribution of low frequency eigenvalues of the membrane randomly attached to the substrate also reveals Wigner-Dyson statistics.
Electrons in graphene are coupled to the deformation of the lattice. Due to the deformation of the lattice the overlap between the electron orbitals is changed, so that the band structure which describes the effective motion of the electrons is altered. Therefore the electron conductivity can be influenced by the oscillations of the membrane. In monolayer graphene, the deformations can be strong, which leads to strong electron-phonon coupling. The coupling of the out of plane phonons to the electrons results in electrostatic potential and pseudomagnetic field which depend on the deformation, acting on electrons 9 . Since electrons move in the oscillating potential created by the phonons, they can scatter with absoprtion or emission of energy from the membrane. If the current-voltage characteristic of the system has some structure at the voltage V , then this structure propagates also to the voltages V + n ω i , with ω i being the set of the membrane frequencies and n an integer number. This phenomenon is similar to the Tien-Gordon effect which arises if the one lead of the sample is irradiated by the photons 10 . In the situation we consider, however, the oscillating potential is created by the membrane itself, and the area where the electrons feel this potential is separated from the leads by potential barriers. Then the conductance exhibits Fabry-Perot resonances at the energies E l . Due to the oscillations with frequencies ω i each peak is enhanced at the energies E l + n ω i . If the conduction modes of graphene at different frequencies do not interact, the presence of the addition to the Fabry-Perot peaks provides us with an opportunity to measure the distribution of eigenfrequencies of the oscillating membrane. This way of detecting the oscillations of the membrane will work for narrow graphene strips. For wide strips, the mixture of Fabry-Perot peaks from different harmonics and their satellites will not give a clear picture.
A number of papers 11-16 studied phonons in graphene. They were focused on phonon spectrum either not influenced by disorder or influenced by weak disorder. We consider the special case, when the disorder changes macroscopic motion of the membrane, namely, the graphene sheet is attached to the substrate at random points.
This Article is structured as follows. In Section II we present the detailed results on the statistics of low-lying eigenfrequencies of the graphene membrane. The numerical method used to obtain the results is found in Appendix A. In Section III the conductivity modification due to the Tien-Gordon effect is considered first for one frequency, then for two frequencies, and finally the influence of the distribution of the phonon modes on the conductivity is addressed. In Section IV we discuss the results. In Appendix B we show the method used to solve the scattering problem for determining the conductivity of the oscillating membrane.

II. PHONONS IN DISORDERED GRAPHENE
Graphene on corrugated substrate can be modeled as a membrane which is randomly attached to some regions of the substrate. We describe this membrane using the boundary integral method and confirm the usual wisdom that the distribution of the spacings between the nearest neighbour phonon energy levels for high levels is the Wigner-Dyson distribution. Moreover, we investigate the statistics for the low energies and find that, after averaging over disorder, it is also described by the Wigner-Dyson type function. Scanning tunneling miscoscopy (STM) experiments show 6 that graphene on a rough substrate is attached to the substrate at some regions, at the hills of the substrate, and it is freely suspended over the valleys. To take this into account, we model a graphene sheet on a substrate as a membrane which is fixed at some regions. Due to external forcing or due to interaction with environment the membrane can vibrate. In this Section, we are interested in the statistics of eigenfrequencies of such membrane.
We assume that the membrane is attached to the substrate strongly enough, so that it can not be detached due to oscillations. The depinning of the membrane on the rough substrate was considered in Ref. 13. The authors concluded that for the large scale of the imperfections of the substrate the membrane indeed remains attached to the substrate, namely for h/h < h δ c /κ(l/h) 2 , where h is the height of the fluctuations, l is the characteristic length of the fluctuations, at which the attachment occurs, δ c is the coupling strength of the graphene to the substrate, estimated for graphene in SiO 2 as ∼ 2meV 17 , κ is the bending rigidity of the two-dimensional layer, κ ∼ 1eV, and the length scale h is the effective thickness of graphene layer h = κ/E 2D , h ≈ 1Å, E 2D is the two-dimensional Young modulus, E 2d = 340N/m. For example, if the height fluctuations are of the order of 5Å, the attachment region of 1nm is sufficient to keep the membrane attached to the substrate. We assume below that the condition is met, and depinning does not occur.
The Lagrangian of the oscillating membrane is with φ(t) being out-of-plane deviation. The first term in Eq. (1) stands for the kinetic energy, ρ 0 is the density of the two-dimensional graphene layer, ρ 0 = 7.6 × 10 −7 kg/m 2 , the second term is the bending energy. The third term is the energy of the deformed membrane with the tension δ, which can be estimated via the Hooke's law δ = E 2d ∆L/L with the relative elongation ∆L/L = δ 0 /E 2d + ξ/L. Here δ 0 is the initial pretension and ξ is the deformation of the layer (the amplitude of the oscillations). Hence we get the ratio of the bending and strectching energies E bend /E str ∼ (h/L) 2 (∆L/L) −1 . For a large membrane, this ratio is small, since, h/L ≪ 1.
In the following, in order to obtain the eigenfrequencies of the membrane, we only take into account the term determined by stretching.
Assuming that the oscillations are harmonic, φ ∝ exp(iωt), we obtain the equation for the spectrum of the membrane with the border Ω (at which the out-of-plane displacement vanishes), ∆φ(q) + λ 2 n φ(q) = 0, for q in Ω/∂Ω , with λ 2 n = ρ 0 ω 2 n /δ. The boundary Ω consists of the external border of the graphene sheet Ω 0 and of the areas {Ω 1 , Ω 2 , . . . , Ω N } where graphene is attached to the substrate. We note that formally λ 2 n correspond to the energy in the equivalent Schrödinger equation widely discussed in the context of quantum chaos. Therefore when we consider the statistics of the eigenfrequencies of the oscillations we always refer to λ 2 n , not to λ n . The eigenvalues of the system (2), (3) are found numerically by using the boundary integral method 18 , see Appendix A for more detail. This method is based on the re-writing the eigenvalue problem for the Laplace equation (2) and its boundary condition (3) in the form of the integral equation for the normal derivative of φ on the boundary. It leads to solving the eigenvalue problem for the system of linear equations in a numerical implementation of the method. We generalized the method for the case of several boundaries (see Appendix A). We tested the program on exactly solvable problems such as an eigenfrequencies of a circle and of a circle with a hole in the middle.
To model graphene on substrate we take a circle with circular regions of the attachments in the middle (to be referred below as impurities or holes). The model is general enough as its classical analogue is non-integrable already if two circles are randomly place inside the large circle. The impurities are placed at random positions. The eigenfrequencies are computed for each configuration. The averaging over disorder configurations is performed to determine the statistics.  (4), N = 444. We note that the distribution is quite noisy. There is no ensemble average for the integrable system, hence for the fixed range of energies there is a fixed number of levels. In the inset we show the distribution for larger range of energy, N = 7180, which is well fit by the Poissobian distribution. For the chaotic membrane (B) the distribution of the high energy levels reveals Wigner-Dyson probability distribution (5). The distribution is computed for the two randomly placed holes. The number of the levels is N = 401442.
Without any regions of the attachment in the middle the system (a circle) is integrable, and its eigenfrequencies obey Poissonian statistics, see Fig. 2B. A circle with two impurities in the middle becomes classically chaotic, namely a small initial difference between trajectories grows infinitely in time. It is well-known that the statistics of high-energy levels of classically chaotic systems (such as non-integrable billiards) is well fit by the same distribution function 8,19 .
The same statistics appears in the disordered solid state problem treated by non-linear sigma-model, see Ref. 19. Indeed, we find that for the high levels the statistics of the nearest level spacings is Wigner-Dyson, with s being the normalized spacing between two levels, see Fig. 2A.
To get higher eigenvalues, we need to discretize all boundaries very accurately, and it becomes computationally consuming to implement more holes. This is why we have chosen the minimal number of impurities (two) for which the system is chaotic.
We note that if one models the attachment regions as points, the system is classically integrable, hence the statistics of the levels for a correctly chosen window of the energies is Poissonian 20 , even though for a smaller window it is Wigner-Dyson. In this manuscript, we remain in the Wigner-Dyson regime.

C. Statistics of low eigenfrequencies of the membrane
The statistics of high-energy levels for our system is not different from any other system described by the Laplace equation. What is particular in a membrane over substrate is that we now have an access to statistics of lowlying levels as well, whereas for electron levels, which typically are the subject of investigation in the level statistics problems, low-lying levels are not accessible and are outside the scope of the problem. For high levels, the random matrix theory appears as collective dynamics of the levels, and the statistics is effectively described by the Wigner-Dyson distribution. For low levels we do not have such averaging over the level positions. To the best of our knowledge, statistics of low-lying levels has not been studied so far since the absence of large number of levels means that there is no ensemble over which one might perform averaging to reproduce the results of random matrix theory (which is where the Wigner-Dyson distribution comes from 21 ). It is not obvious what distribution the lowest lying levels obey.
We study this statistics in more detail than the one of high energy levels for different numbers of impurities. We treat the impurities as quenched disorder, hence we We averaged over configurations with different numbers of the impurities, N = 5, 6, 7. The oscillation in the density of levels are more pronounced for smaller size of the impurities (the number of impurities is the same). It is in agreement with the simple-minded argument that small impurities have less influence on the spectrum. We argue that that the oscillations in the density of low energy states do not disappear for the case of non-integrable boundary as the oscillations simply reveal the geometric shape of the membrane. Therefore for a non-integrable outer boundary of arbitrary shape the oscillations would be shifted according to the first few energy levels of the membrane without disorder. We obtain the Wigner-Dyson distribution for the level spacings of the first four-five low lying energy levels slightly depending on the number of impurities in the system, Fig. 4A. We average all obtained numbers for different number of impurities and get the distribution which is close to the Wigner-Dyson one for different sizes of the impurities, Fig. 4A,B. Let us note that we fit the numerical data with the dependence which is the same as in Wigner Dyson probability distribution. In the Wigner-Dyson probability distribution for high energy levels, the parameter s of the Wigner-Dyson probability function is normalized by the average level spacings δλ 2 , namely, s = (λ 2 n+1 − λ 2 n )/δλ 2 . In the case of low levels, the mean level spacing for the region of fit oscillates, see . We notice that it reproduces well the Wigner-Dyson distribution function (5). The normalization is done for the averaged to the constant density of states.
the fit in Fig. 4 is close to the Wigner-Dyson fit with the parameter s normalized by the mean value of density in the region of fitting. Furthermore, we explore the statistics of the position of the first energy level in the chaotic disordered cavity, Fig. 5.
First we compare the distribution of the energy levels for one impurity obtained by solving the corresponding equation numerically and computed from the Born approximation. We show the difference between conventional way of computing the shift of eigenenergies due to the disorder for the electron system with an impurity of finite potential V and the phonon system with fixed to zero boundary condition on the circumference of the impurity.
In the Born approximation, the shift of the energy of the first level for the impurity Ω i cdplaced at the random position at the distance ρ from the center is where V is the potential characterizing the degree of at-tachment of the membrane to the substrate, and J 0 is the Bessel function of the zeroth order, which represents the solution of the Laplace equation in a circle with the zero boundary conditions. The probability to find the impurity at the distance ρ from the center is We can get the dependence ρ(E 1 ) from Eq. (6) and consequently the probability distribution P 1 (E 1 ) in the Born approximation. Fig. 5A demonstrates that the Born approximation does not correctly describe the problem with the zero boundary conditions. The discrepancy can be attributed to the fact that we need an infinite potential to reproduce zero boundary conditions at the finite circumference of the attachement region, while the Born approximation is only applicable for V r 2 ≪ 1. For this reason, we proceed with numerical solution for a more complex system, with two and more impurities.
As we would expect from simple-minded considerations analogous to those in the quantum mechanics, the eigenfrequency can move only up in energy with increasing the amount of disorder (amount of impurities). Thus, the minimum eigenfrequency corresponds to the one of the ideal circle. On all our plots we display the distribution of the first energy level taking into account this finite energy offset. The distribution of the first level is smoothened with increasing the disorder. In Fig. 5B we show the evolution of the distribution with increasing number the impurities from N = 1 to N = 2 and N = 7. We clearly see how the distribution becomes more broad. In Fig. 5C it is seen that the distribution becomes more or less similar for larger number of impurities, N = 5, 6, 7 Microwave experiments 22 are an established way to measure the eigenfrequencies for the wave equation, system given by Eqs. (2), (3), in various geometries. Therefore we can propose studying microwave radiation in the cavity as a direct way of measuring the computed distributions. The boundary condition is imposed by placing highly reflective material. One can think of a simple model: a big cavity with randomly placed disks. While changing the positions of the disks and measuring the first few eigenfrequencies, the distribution for the average density and the spacing between nearest neighbors can be recovered. In this system, the lowest-lying eigenvalue statistics must be the same at the one we described above.

III. PHONON-ASSISTED TIEN-GORDON EFFECT
In this Section, we investigate how the low-lying eigenfrequencies of the oscillating graphene membrane can be experimentally assessed by means of conductance measurement.
The electrons traversing graphene are coupled to the phonon motion of the membrane as 11 where g 1 is the deformation potential, g 1 ∼ 20−30eV , affected by screening; g 2 is the coupling to pseudomagnetic fields, g 2 ≈ 1.5 eV , and f (u ij ) = 2u xy + i(u xx u yy ). For the second-order process involving out-of-plane phonons one obtains where the similarity sign shows that we do not take into account the angle of the scattering between different phonon modes. We can make the following estimations of the parameters in Eq. (9), q n ∼ n/L, g 1 ∼ 10eV (including some screening). The amplitude of the oscillations φ can be estimated in several different ways: as the amplitude of zero-point motion, φ = ( /(2M ω)) 1/2 , or as the amplitude of thermal fluctuations, M ω 2 t φ 2 = κ B T , or it can be varied by changing the amplitude of the oscillations to which one pumps the energy (for example, by placing the sample on the piezosubstrate). Let us consider electrons in graphene subject to oscillating potential. Electrons can emit or absorb quanta of energy corresponding to the oscillations due to the electron-phonon coupling. Hence an electron with the energy E which is coming to the region with a harmonically oscillating potential with frequency ω can be transmitted and reflected to the energies E ± n ω. A similar effect of the transmission of the electron from the energy E to the energy E + n ω due to microwave irradiation of the sample is called the Tien-Gordon effect 10 . Note however that there is an important difference, since in the Tien-Gordon effect the bias voltage oscillates. Our system is formally equivalent to the conductor with the oscillating barrier between the leads. Electrons in the leads are in thermal equilibrium described by the distribution functions f L (E) and f R (E + eV b ), with V b being the difference in the potentials between two leads.
The electron propagating from the left to the right can be transmitted from the energy E 1 to the energy E 2 or reflected to the energy E 2 . We denote the corresponding transmission and reflection coefficients from as T E1→E2 and R E1→E2 . For the electron propagating from the right to the left the similar coefficients are denoted as T E1→E2 and R E1→E2 . Transmission and reflection probabilities are determined from the solution of the scattering problem similarly to Ref. 23 (see Appendix B). The total transmission from the left to the right is determined by the following expression: where q is the transverse momentum. The first term represents electrons moving to the right. The second and the third terms represent reflected electrons moving back to the left and electrons transmitted from the right, respectively.
The change of integration variable gives: For the case of the transmission to different energies each wave function should be normalized by the flux in the direction of propagation, namely, where Ψ is wave function andĵ is current operator. The transmission and reflection probabilities defined in such basis obey generalized unitarity condition which is the conservation of the flux during the scattering process from the mode with energy E to the modes E + i ω (see also Appendix B, Eq. B8): Hence we can simplify the two first terms in the expression (10): We introduce the coefficients of total transmission from an energy E as E→E+i ω (14) and get the expression of the generalized Landauer formula 24,25 dE B. Transmission through the oscillating graphene membrane In further treatment, to get the conductivity of graphene with an oscillating barrier, we consider the rectangular graphene sheet. First, we get the general dependence of the conductivity on the strength of the oscillations which is determined by the influence of the substrate and the applied gate voltage (the Fermi energy in terms of the density of the electrons on the graphene flake is . This dependence is not crucially sensitive on the shape of the graphene sheet. Second, we connect the conductivity of oscillating graphene to the distribution of the phonon frequencies (which determine the frequency of the oscillations in our setup). For generic chaotic systems such as quantum billiards it is well known that the geometric shape should not influence the distribution of the eigenfrequencies after averaging over the ensemble (in our case, over disorder configuration).
As we have shown in Eq. (9), the oscillations of the graphene membrane enter the Dirac equation which describes electrons in graphene as oscillating electrostatic potential as well as pseudovector potential. Typically, the effect of the electrostatic potential is more significant 26 , and therefore we disregard below the effect of the pseudomagnetic fields. The electron Hamiltonian This problem was considered before in a number of studies: to investigate the direct analog of the Tien-Gordon effect in graphene with an oscillating lead 27 , to look at the properties of transmission amplitudes depending on the angle of the incoming wave 23 , and to implement quantum pumping via evenescent modes in graphene 24 . In our work we follow closely Ref. 23 solving numerically  the system of equations for transmission and reflection coefficients (as discussed in Appendix B), but we go beyond this work by summing up the transmission coefficients over the different transversal momenta q to obtain the conductance of the system. Our final goal is to relate the spectrum of phonon frequencies to the features of Tien-Gordon effect. The Tien-Gordon effect for graphene is characterized by the dimensionless parameter which is the ratio of the electron-phonon coupling V e−ph to the energy of the oscillations of the phonons ω n , v n = V e−ph ω n .
Based on the expression for electron-phonon cou- where we took into account the expression for the smallest frequency of the oscillating circular membrane λ 0 ≈ 2.4/L. For the case of the thermally driven oscillations with the amplitude of the oscillations with γ characterizing a tension of the membrane, γ = ∆L/L. For example v ∼ 1 for the realistic parameters of graphene L = 300 nm with the tension γ = 10 −4 with the coupling potential g 1 ∼ 10 eV at the temperature T = 300 K 28 . The parameter of the strength of the oscillations can be substantially modified by the tension and the size of the membrane, Fig. 7. The tension depends on the fabrication process. While increasing the tension the influence of the oscillations of the membrane on the conductivity decreases, v decreases. With increasing the size of the membrane v decreases as well. Decreasing temperature reduces the influence of the oscillations on the conductivity as the coupling is proportional to temperature. We would like to note that the higher oscillations of the membrane are less coupled to the electron motion, namely v n v 0 ∼ n −1/2 (20) We determine the influence of the oscillations of the potential on the conductivity using the generalized Landauer formula (15). For low temperatures, κT ≪ V bias , κT ≪ ω, we can approximate the distribution function in the leads by the constant, then the current through the system is determined by The transmission coefficient is already a sum over transverse momenta and over all final energies E 0 + i ω which scatter from the energy E 0 . Here we also took into account the inversion symmetry of the system setting T E LR = T E RL . For weak oscillations, v ≪ 1, we can make an analytic estimation of the effect (see also Refs. 23,24). The expression for the transmission probabilities bringing an electron from the energy here and below we use the simplified notations for the indices of the transmission coefficients, namely, denote the energy by numbers which are i in the expression E 0 + i ω). We can disregard all other probabilities as they are smaller in the parameter v. The current becomes then where I 0 is the current in the absence of the oscillations.
For moderate values of v ∼ 1 we solve the scattering problem numerically and then determine the conductivity by means of the generalized Landauer formula (15). In Fig. III B we present the dependence of the conductivity on the strength of the electron phonon coupling v for different energies on the graphene flake. We notice that if the strength of the oscillations grows, the condictivity decreases and becomes independent on the width of the graphene sheet for large values of v. The absolute value of the drop of the conductivity is approximately similar for different values of gate voltage on graphene.
For narrow graphene strips the conductance is quantized with increasing the gate voltage. Due to the oscillations the steps in conductance at gate voltages V step are reflected at the gate voltages V step + ω, see Fig. 10.  Fig. 9). Pair creation is important at low gate voltages where conductivity significantly increases. Here the oneparticle approximation with Dirac equation does not work and more precise consideration using quantum field theory should be applied.

C. Klein tunneling
The dependence of the conductivity on the gate voltage at high values of v (v > 1) and low gate voltages exhibits unexpected growth which reveals the existence of pair creation, see Fig. 11, as we explain below.
Graphene is a condensed matter material where Klein paradox can be measured directly 29,30 . Klein paradox arises in scattering of Dirac electrons on a high enough potential step: The tunneling probability is negative, and the reflection probability is greater then one. Another closely related counterintuitive result is that the transmission probability remains finite for arbitrary high and long barriers and does not depend on the length of the barrier. Both results come from the solution of the Dirac equation. It is a consequence of one-particle treatment of Dirac fermions. A fully consistent consideration of the problem requires the quantum field theory approach where many-particle states are taken into account. It was done for example in Refs. 31,32 and discussed in the context of graphene in Ref. 33.
Klein tunneling can be represented as tunneling with a creation of an electron-hole pair. When a particle-hole pair is created, the particle can be reflected backwards by the potential, which gives a positive contribution to the reflection probability and can make it larger than one. The hole moves in the same direction as the transmitted particle would move, thus contributing a negative term to the transmission probability 32 . In classical terms the motion of the hole is the same as motion of the particle in the opposite direction. The above picture stems from the single-particle (quantum-mechanical) understanding of the process. In fact, the particle-hole creation also changes the ground state of the system and introduces vacuum polarization. When these effects become significant, the calculation of the transport cannot be done in the simple way as above as now the potential itself in turn depends on the density of particles and holes.
The equilibration of the two continua (the electron and the hole one) happens very fast, so that it is difficult to observe the pair-creation in a p-n junction directly as a characteristic time is very short. In a time-dependent potential the pair creation occurs as well, as it is proven for the general case of slowly varying time-dependent Dirac equation 34 . An attempt to consider the time-dependent Dirac equation in the context of graphene was done in Ref. 35, however, the authors did not take into account possible transmission to different energies.
Obviously for low values of the gate voltage we encounter Klein paradox for the time-dependent potential. Indeed, an incident particle with the energy E scatters to the energies E + n ω. For high strength of the oscillations v the scattering to the states with large n can be strong. It imitates the scattering on the potential step, with n ω playing the role of the height of the potential barrier. If the energy |E − n ω| > E then the pair creation is allowed. The probability of the process becomes appreciable for sufficiently strong oscillations, v > 1. We identify the region of the parameters (V g , v) where the conductivity achieved from the scattering problem for the Dirac equation shows increased conductivity due to Klein paradox, Fig. 12. The increase of the transmission is a measurable effect, described still by one-particle Dirac equation for the moderate values of the tranmission increase. In the regime where transmission diverges, our results, based on single-particle Dirac equation, are not valid, and more precise quantum field theoretical treatment is needed in order to obtain correct value of the conductivity due to the pair-creation mechanism. Nevertheless, we can claim that the conductivity at small Fermi energies increases at low values of the gate voltage and high values of the amplitude of the oscillations.
Thus, the dependence of conductivity of the gate voltage presented at Fig. 11 is correct at large enough gate voltages (as it follows from Fig. 12A), where the single particle approach to determine the conductivity is valid, the transmission is not subject to many particle process which involve the pair creation. We also note that the plots representing the dependence of the conductivity on the strength of the oscillations v, Fig. III B, are made for the gate voltages large enough that the pair-creation process in not important.
The increase of the conductivity due to pair creation can be measured for small graphene samples where the coupling strength is large enough. Namely, we can approximately connect the region on the diagram Fig. 12A to the density of graphene n and the length of the membrane using the relation for the Fermi-energy E F = v F √ πn, frequency ω ∼ 2.4 δ ρ0 1 L and the coupling strength given by Eq. (19). After these substitutions the region where the Klein tunneling can be significant is: see Fig. 12B.

D. Tien-Gordon effect for several modes
Thus far, we only considered a monochromatic excitation of the membrane. In this Subsection, we look how the presence of several oscillation frequencies of the membrane influences the current. If several modes are essential, the largest contribution is given by the lowest modes, since for them v is the largest. The higher is the oscillation mode the less is the correction to the conductivity is contributes. As we show below there is no region for parameters in our problem where the interference between different frequencies is important for the conductivity, and the modes contribute independently: The influence of all modes is just the sum of individual contributions.
In this subsection we look first at the conductivity in the presence of two frequencies. We solve numerically the scattering problem for the Hamiltonian, H = −i v F σ · ∇ + V 0 + V e−ph (cos(ω 1 t) + cos(ω 2 t + ψ)) .
(24) using the same method which we set up for the case of one frequency.
For weak electron-phonon coupling, v 1 , v 2 ≪ 1 the influence of the two oscillations is additive in the lowest order approximation, see Eq. (22). Moreover, we see from the results of numerical treatment, displayed on Fig. 13, that this additivity property is approximately valid even for moderate values of the oscillation strength v i ∼ 1. Namely, the relative change of the conductivity for the oscillations with two frequencies (G(ω 1 , ω 2 ) − G 0 )/G 0 , G 0 = G(ω = 0), equals to sum of the relative changes of the conductivity for the oscillations with one frequency (G(ω) − G 0 )/G 0 , In the region of low energies the Klein tunneling is preserved and the conductivities have the additive property, Fig. 14. The slight deviation from the additivity occurs since we do not take into account evanescent modes for determining the conductivity. The conductivity for the case of two modes of the oscillating potential. The Fermi energy is fixed. We assume that the modes are characterized by the same coupling to the phonons V el−ph and different frequencies (the coupling strength depends on the phonon frequency, but for the illustration of the interaction of two frequencies this dependence is not important). Hence the dependence on v translates directly to the dependence on the frequency as 1/ω. A,B,C:The conductivity dependence on frequency ω1 at fixed frequency ω2 is shown. We show G(v1, v2 = const)/G0 and for comparison G(v1)/G0 and G(v2)/G0 as well as (G(v1) + G(v2))/G0 − 1. We see that the effect of two frequencies is almost additive, A,B. Fig. C presents G(v) for different sets of parameters: at smaller energies the correction to conductivity is more significant.
Taking into account this property, we can determine the disorder-average correction to the conductivity due to the presence of different phonon frequencies, The influence of the high-frequency modes can be easily taken into account since their density of states is constant, dN/d(ω 2 ) = C, which leads to P high (ω) = dN/dω = 2Cω, and the strength of the oscillations is small v = V el−ph / ω ≪ 1 for high frequencies, hence the relative change of the conductivity for every mode is proportional to v 2 , and for the range of high energy modes it is where the cut-off frequency ω 1 , which is the position of the first energy level of the oscillations, is given by the geometrical size of the membrane, and ω max is the maximum frequency of the oscillations which is determined by microscopic properties of the membrane at zero temperature or by κ B T / at finite temperature (we assume that the maximum frequency is given by comparission the amplitude of zero point fluctuations and the amplitude of the thermal fluctuations, namely M ω max φ 2 th ∼ κ B T ∼ ω max ). We also take into account the decrease of the coupling with increasing the frequency of the oscillation, Eq. (20). The constant density can be estimated from the general consideration: The density of states of the oscillating membrane according to the Weyl formula is dN/dλ 2 = (γ/ρ 0 )dN/dω 2 = (1 − r 2 n)L 2 /2 ≈ L 2 /2, where n and r are the number of holes and the radius. Therefore the oscillations at high frequencies can decrease the conductivity by We can express the change of the conductivity for the real sample: The larger is the membrane the smaller is the conductivity correction. For large values of the electron phonon coupling V el−ph > ω 1 one observes a noticeable decrease of the conductivity due to first several modes of the oscillations. For example, for the graphene membrane of the length L = 500 nm and the tension γ = 0.3 · 10 −3 the correction is δ high levels σ(ǫ) σ0 disorder ∼ 0.015. Let me note that the correction is very sensitive to the stretching. G/G 0 E F / − h ω G(ω 1 ,ω 2 )/G 0 G(ω 1 )/G 0 G(ω 2 )/G 0 (G(ω 1 )+G(ω 2 ))/G 0 -1 FIG. 14: Dependence of the conductivity on the Fermi energy in the case of a two-mode oscillating barrier. We clearly see that the additivity property is preserved even for low energies where Klein tunneling plays an important role. The deviation from the additivity occurs due to finite size of the system and lack of account of evanescent modes.

IV. DISCUSSION
In this Article we investigated the statistics of lowenergy spectrum of the graphene membrane suspended over disordered substrate and the effects of this statistics on the conductivity. We modeled the graphene sheet as a membrane attached to randomly located circular areas. The phonon frequencies in this system are given by eigenvalues of the Laplace equation with corresponding boundary conditions. It is known that the high-energy eigenvalues in this system obey the Wigner-Dyson statistics. This is valid for a generic chaotic system, where the statistics is found as the result of the averaging over the energy with the constant density of states. It was less obvious what the statistics of low-lying eigenfrequencies is, since conventional analytical methods are not applicable to this situation. For a single sample (with a given impurity configuration) one can not define the statistics, since one can not average over energy. However, in a disordered system one can instead average over impurity configurations (positions of the attachment areas). We performed this averaging numerically and found that the statistics of the nearest level spacing between the first five-six levels obeys the Wigner-Dyson probability distribution. The important difference between high-energy and low-energy case is that the density of eigenvalues is not constant for low energies -rather, it oscillates for the first few levels, which is a remnant of the existence of the positions of the eigenvalues of a clean circular membrane. At high energies, the density of eigenvalues saturates. Note also that the eigenvalues of a disordered membrane can never lie lower than the lowest eigenvalue of the clean mem-brane, which creates an offset energy in the density of states. While changing this cut-off, in physical terms changing the size of the membrane, the statistics of the levels is not changed much.
Our results can not be described by the non-linear σmodel, a standard tool to deal with disordered systems, because in our system the disorder is not strong enough. We note that whereas in the non-linear σ-model effective bosonic fields are introduced after averaging over the disorder, and these fields are confined to a sphere (compact sector) due to strong disorder, in our case of low-lying levels such fields would not be confined to the sphere. The indirect evidence for this is the oscillation in the density of levels for the first few levels.
The phonon levels could be accessed directly, for example, by coupling the graphene membrane to other degrees of freedom like microwave photons. However, these experiments presently look too involved, and therefore we investigated the effect of the phonons on electron transport through graphene. The coupling between electrons and phonons is provided by the deformation potential. The electron transport in this situation shows the phonon-assisted Tien-Gordon effect. Due to the oscillations of the membrane with the frequency ω the electron can be transmitted from the energy E to the energy E + n ω. At strong electron-phonon coupling and a small membrane we find that the influence of the phononassisted Tien-Gordon effect on the conductivity can be of order of ten percent. If the frequency increases, the influence of the phonon oscillations becomes less significant, and thus the lowest frequencies of the phonons have the most crucial influence on the change of the conductivity. The modes affect the current independently, and thus the influence of the superposition of the modes can be approximated by the sum of the influences of an individual mode.
The Tien-Gordon effect is the most pronounced in the case of the presence of structure in the conductivity. The structure, in our case the Fabry-Perot resonances, appears in the voltage dependence of the conductivity in the case when tunnel barriers are located at the ends of the sample. Then every resonance at the energy E has its sattelites at the energies E +n ω. This can be also a way to measure the phonon frequencies. This method can be applied only to the case of the narrow samples, as for real two-dimensional sample the Fabry-Perot resonances have different positions for different transverse momenta, and summing over the transverse momentum makes the Fabry-Perot peaks and their satellites indistinguishable.
At low gate voltages and strong electron-phonon coupling the oscillation potential can create the potential step large enough to detect Klein tunneling process, namely, increase of the conductivity due to electron-hole pair creation. There is a range of parameters where the single particle Dirac equation is applicable to describe this process. When the reflection and transmission probabilities start to diverge, many-particle quantum field theoretical treatment should be performed. Unexpected increase of the conductivity at low gate voltages can be a signal of strong electron-phonon coupling. These effects lie outside the scope of this Article.
Let us note that a suspended graphene membrane is not the only situation when low-lying eigenvalues of the Laplace equation can be assessed. For instance, our treatment also applies for a system of microwave cavities with randomly positioned reflective disks inside. The distribution of photon frequencies in these cavities will have the same features as described in this Article.
we are interested in a problem with multiple boundaries, namely, a region with holes inside it.
Let us first introduce the general formalism of the boundary integral method and then specify the problem for the case of several boundaries.
For numerical implementation the equation (A3) is discretized, the boundary is discretized, s (i) , and the corresponding values of u (i) = u(s (i) ) are taken. The condition that the solution exists is det Q n(i,j) − 1 = 0. (A7) We vary λ and monitor the left-hand side of the condition (A7), thus finding the eigenfrequencies of the Laplace equation. Let us note that numerically it is more efficient and more precise to look at the eigenvalues of Q n(i,j) − 1 then to look at the condition for the determinant.
The boundary integral method can be simply generalized for a region with holes. Let us represent the boundary as ∂Ω = ∂Ω 0 ∪ ∂Ω 1 ∪ · · · ∪ ∂Ω N where ∂Ω 0 is outer boundary of the area and {Ω i }, i = 1, . . . , N are holes inside of the area. The directions of the boundaries and normal vectors to the holes are consistent with the direction of the outer boundary. Now the function u(s) on the boundaries is represented as (u ∂Ω , u ∂Ω1 , . . . , u ∂ΩN ). For further convenience of the notations we use the indices Ω and Ω i instead of the precise notations ∂Ω 0 , ∂Ω i . The system of integral equations is: For example we present the kernels for our model of the circular outer membrane and the holes of radius r placed in the random positions inside of the membrane. We parametrize the hole position by the position of the center and the angle (r i , θ i ). The position on the outer membrane is given by angle φ R , on the small circle by φ r . We take the radius of the outer membrane R = 1 and for the holes r. All parameters are shown at Fig. A. By τ we denote the distance between the points for which the kernel element is given, τ (s, s ′ ) = |q − q ′ |.
• For q at the large circle and q ′ at small circle is (the negative sign originates from the circulation in another direction for the integral over the small circle) 1 (kτ (s, s ′ )), • for q at the small circle and q ′ at large circle is (the negative sign originates from the direction of normal) Q ΩiΩ k (s, s ′ ) = − iλ n 2 r + r i cos(φ r + θ i ) − cos(φ r + φ R ) τ H 1 (kτ (s, s ′ )) (A11) • for both q and q ′ at the small circle (one negative sign originates from normal, and another one from the direction) Q ΩiΩj k (s, s ′ ) = − iλ n 2 r i cos(θ i + φ ri ) − r j cos(φ rj + θ j ) − r cos(φ ri − φ rj ) τ H 1 (kτ (s, s ′ )) (A12) The same system of equations was considered in Ref. 23. We extend the analysis of that publication, since at the end we determine the full transmission probability, i.e. we sum over all propagating modes. We generalize the problem to several frequencies as well. For example for the case of two frequencies ω 1 and ω 2 , the scattering will happen to the energies E + l ω 1 + n ω 2 . If for the treatment of the system with only one frequency one needs the system of the equations of the size 4N , then for the treatment of the system with two frequencies it is 4N 2 . This growth of computational complexity prevented us from looking at the strong coupling region where both parameters V / ω 1,2 are large.