Critical radius and temperature for buckling in graphene
Abstract
In this work, we find an analytical flatmembrane solution to the saddle point equations, derived by Guinea et al. [Phys. Rev. B 89, 125428 (2014)], for the case of a suspended graphene membrane of circular shape. We also find how different buckled membrane solutions bifurcate from the flat membrane at critical temperatures and membrane radii. The saddle point equations take into account electronphonon coupling and this coupling provides a residual stress even for a flat graphene layer. Below a critical temperature (which is exceedingly high for an infinite layer) or above a critical size that depend on boundary conditions, different buckling modes that may be the germ of rippling appear. Our results provide the opportunity to develop new feasible experiments dealing with buckling in small suspended graphene membranes that could verify them. These experiments may also be used to fit the phononelectron coupling constant or the bending energy.
pacs:
63.22.Rc, 68.60.Dv, 65.80.Ck, 62.20.mqI Introduction
The discovery of graphene has spawned one of the most fertile scientific fields of today. Its unique electrical nov05 () and mechanical properties lee08 () promise to revolutionize current technology, making extremely important to characterize the properties of graphene in order to design and optimize new devices. Among mechanical properties, the presence of ripples in suspended graphene mey07 () has aroused a great theoretical effort focused on their explanation from first principles fas07 (); abe07 (); kim08 (); gazit09 (); sanjose (); gui14 (); gonzalez (), and by using simple statistical mechanics models jsm12bon (); pre12bon (); RuizGarcia (); prb12bon (); schoelz15 (). Thermal and external stresses have been used to create large ripples (wrinkles) in a layer of graphene which is suspended over trenches bao09 (). These ripples have been explained using classical elasticity bao09 () and seem to be qualitatively different from the nanometer sized ripples observed in mey07 (); kiri (); locatelli (). The latter tend to be considered inherent to suspended graphene. Even when ripples in graphene are not fully understood, one of the most promising approaches to ripple formation is the quantum mechanical study of the membrane and the phononelectron interaction as a mechanism for rippling gazit09 (); sanjose (). Since the graphene Debye temperature is K, the quantum treatment of graphene seems justified at room temperature. Perturbation treatments point to the vanishing of the renormalized bending energy of the membrane as a possible mechanism for ripple formation sanjose (); gonzalez (). The destabilizing effects of quantum fluctuations on the bending rigidity of crystalline membranes (without electronphonon interaction) have been considered in kats (). Within equilibrium theory and without phononelectron interaction, the effect of (Matsubara) frequencydependent renormalization of the anharmonic coupling after elimination of inplane phonons has been considered in amo14 ().
Recent experimental studies emphasize buckling effects in suspended graphene, as buckling is more easily characterized than ripple formation and can strongly affect the design of new graphenebased devices. Longrange buckling of the graphene sheet is pronounced near defects and dislocations SN88 (); CN93 (); leh13 (); kot14 (); rob14 (); BCGW15 (). In lin12 (), buckled bilayer and monolayer graphene are fully clamped to circular holes in an electrode. Buckling is induced by the fabrication process so that the sheets are buckled and strainfree at zero applied electric field. In the experiment, an external electric field pressures the membrane and induces a sudden change in the buckling direction (a “snapthrough” effect). A relation between the critical pressure and the bending rigidity of bilayer graphene provides a measure of the latter. Even when the graphene layer is not buckled initially, a strong enough electrostatic force can produce irreversible buckling svensson11 (). A spontaneous buckling effect (“mirror buckling”) is also found by Xu et al. in xu14 (). This effect is systematically studied in schoelz15 (), in a complementary approach to that in lin12 (). In schoelz15 (), the authors study a suspended graphene monolayer using scanning tunnelling microscopy (STM). The microscope tip is located on the centre of an initially nonbuckled layer. The STM keeps a constant current and a variable potential between the tip and the sample. Once the current is fixed at a sufficiently high value, the sample buckles when the potential increases, similarly to the experiments in lin12 (); svensson11 ().
These experimental efforts have their theoretical counterpart through different approaches to the buckling of two dimensional (2D) layers such as graphene. A broad overview of the effects of strain in graphene and its relation with its electrical properties can be found in amorim15 (). The study of spontaneous buckling in graphene due to doping is carried out in gazit09 (). In polymerized membranes, buckling occurs below a critical temperature Guitter88 (). Phenomenological models of crystalline membranes coupled to spins undergoing Glauber dynamics also indicate that there are critical temperatures below which the membranes buckle jsm12bon (); RuizGarcia ().
In gui14 (), Guinea et al. have developed a model for suspended graphene based in the coupling between flexural phonons and electrons. They consider a 2D membrane in equilibrium that is embedded in a dimension space, where there are dimensional outofmembrane displacements ( in the physical situation). Then they eliminate the inplane phonons, ignore the frequency dependence of the resulting couplings, use the selfconsistent screening approximation and perform a saddlepoint analysis of the free energy in the limit as . The resulting saddlepoint equations (SPEs) are timeindependent and will be analyzed in the present paper. They consist of one von Kármán type plate equation for the physical outofplane displacement coupled to two equations for two auxiliary fields: one “scalar” stress associated to membrane curvature and a field associated to charge fluctuations. We expect these equations to describe buckling and ripples of the graphene membrane in a stable stationary (equilibrium) configuration after all possible transients have decayed. We find flat membrane solutions with constant nonzero auxiliary fields. The linearization of the stationary SPEs about these solutions accompanied by appropriate boundary conditions provide an eigenvalue problem that yields critical values of temperature and membrane size corresponding to the bifurcation of buckling states from the flat membrane. Why?
The stationary SPEs provide stationary solutions to not yet derived dynamic SPEs that should describe the graphene membrane out of equilibrium. Among them, we have found the stationary flat membrane solution. Other solutions may bifurcate from it at appropriate parameter values. To find them, we have to solve the eigenvalue problem that governs the linear stability of the stationary flat membrane solution to dynamic SPEs. The corresponding eigenvalues give the growth of disturbances about the stationary solution. We do not know the precise shape of the dynamic SPEs but we may surmise that stationary solutions bifurcating from the flat membrane appear when these eigenvalues are zero. But linearized dynamic SPEs with zero eigenvalues are the same as linearized stationary SPEs, which we know from gui14 (). These linearized stationary SPEs have nonzero solutions only for particular values of the constant auxiliary fields that correspond to critical temperature and membrane sizes. The eigenmodes solving these equations give the shape of the buckling states near bifurcation points. A combination of eigenmodes may characterize ripples in graphene.
The rest of the paper is as follows. In section II we briefly revise the derivation of the stationary saddlepoint equations gui14 () and write them in real (not Fourier) space. In section III, we find the flat membrane solution of the SPEs for constant auxiliary fields. In Section IV, we linearize the plate equation that forms part of the stationary SPE about the flat solution and add appropriate boundary conditions for a finite circular graphene sheet. Thus we obtain an eigenvalue problem for critical values of the constant auxiliary fields. From the critical auxiliary fields, we get critical values of temperature and membrane size at which nonflat buckled solutions may bifurcate from flat ones. The linearized plate equations are solved for a circular monolayer of graphene with two different boundary conditions corresponding to a free graphene layer and to a clamped sample. We obtain different eigenmodes (corresponding to vertical deformations of the layer) together with their respective critical temperatures and radii. The critical radii allow us to predict the minimal size that allows a graphene layer to buckle. Moreover the combination of the different bifurcating modes could be used to characterize ripples in graphene. The last section contains our conclusions.
Ii Model and saddlepoint equations
In gui14 (), Guinea, Le Doussal and Wiese have considered the graphene sheet to be a 2D membrane embedded in a larger space of dimension and interacting with copies of a free Dirac fermion. The real physical system has (four flavors, two valleys and two spins) and but Guinea et al have derived saddlepoint equations involving the outofplane displacement in the limit as . Here we analyze buckling of the graphene membrane using the saddlepoint equations. Let us briefly recall the model, the saddlepoint equations and their meaning. The model Hamiltonian consists of:

The deformation energy of the graphene sheet is the sum of kinetic energy, , and of curvature and elastic energy, , with
(1) Here , , are the inplane phonon displacements, , are the outofplane flexural phonon modes, , and are the Lamé constants, and is the bending energy.

The energy of the free Dirac fermions (in units such that ),
(2) Here are the Pauli matrices, is the Fermi velocity of Dirac electrons in graphene such that, in these units, eV ( Å is the side of a graphene hexagon).

The electronphonon coupling
(3) (4) where is the the charge distribution on the graphene layer, is the equilibrium carrier density, and is in the range between 4 and 50 eV.

Finally the Coulomb interaction has the form
(5) where is the Fourier transform of the electrostatic potential, , is the dielectric constant of the environment and is the charge of the electron.
Once all interactions are defined, the resultant Hamiltonian is integrated over the inplane phonons, leading to a coupled theory of flexural phonon modes and electrons. The frequency dependence of the resulting couplings is ignored, as it is important only for temperatures below 90 K amo14 (), that we do not consider here. In this process, the Coulomb interaction becomes
(6) 
The effective Hamiltonian, depending only on the electrons and flexural phonon modes, is used to build a Matsubara equilibrium partition function. For later convenience, this partition function is transformed using two fluctuating auxiliary fields, and , defined through the relation,
(7) 
Here , , , and table 1 gives the dimensions of parameters and variables. In (7), is
(8) 
where is the transversal projector. With this definition, is related to the Gaussian curvature () through the relation . In (7), is a local linear combination of and the charge disturbance , whereas is a nonlocal linear combination, as it involves the potential energy of one electron in with respect to the charge disturbance distribution in , . The action corresponding to the Matsubara partition function is
(9) 
where is the chemical potential. We now decompose the vertical displacements into an average and a fluctuating part, , and assume symmetry breaking in the direction , . Integrating over the fermions and the fluctuating part of the vertical displacements, the action is gui14 ()
(10) 
Dimensions of variables  
To obtain the saddlepoint equations, we vary (10) with respect to , and , thereby obtaining
(11) 
(12) 
(13) 
(14) 
where we have replaced and not yet made as in gui14 (). Eq. (11) provides the average fields and in terms of . The average auxiliary fields are related to the averages of the charge distribution (4) and of the field (8),
(15) 
by the linear equation (7), which also holds for the averages of the corresponding quantities. (12), (13) and (14) are von Karman plate equations (they appear in Fourier transform form in gui14 ()). The stress tensor is generated by the average auxiliary field and it is automatically in equilibrium as . is the transversal projector, and and are the bosonic and fermionic Matsubara frequencies. The matrix Green function appearing in (11) satisfies
(16) 
and can be computed from (16) as
(17) 
It is clear from (17) that . Using this translation invariant Green’s function, we can write (11) as
(18) 
Henceforth from this equation we suppress the subscript 0 in all average quantities.
Let us now set as in gui14 (). It is important to note that the flat membrane is always a solution of (12). However, the right side of (18) is not zero for and therefore the trivial solution () is not a valid solution of the SPEs. As we will see in the next section, and nonzero and is a valid solution. Nonzero auxiliary fields and provide a residual stress that induces plate buckling and rippling below a certain critical temperature and over a critical plate size. In a mathematically related phenomenon, a growing bacterial biofilm can be modeled as a plate with a growth tensor that modifies the elastic part of strain, acts as residual stress and may trigger ripples (called wrinkles in that application) esp15 (). The critical growth term may then be sought by solving an appropriate eigenvalue problem for coming from the linearized plate equations der08 ().
Iii Flat membrane with constant auxiliary fields
We now seek simple solutions of equations (12)(14) and (18) with , (a constant), (a constant) and . Using that
(19) 
and (17), Equation (18) now becomes
(20) 
where ( is an infrared cutoff resulting from the graphene sheet radius) is its regularized version and the constant appearing in (6) has been neglected. The right side of (20) is computed separately in Appendix A with the result
(21) 
where is an ultraviolet momentum cutoff ( is the side of a graphene hexagon). The term in (21) is proportional to the area of the Brillouin zone, and it will play a key role in the computation of the critical temperature and radius for the graphene layer. We discuss in Appendix B other possible interpretations of (20) that do not involve setting the chemical potential and turn out to be unphysical. From (20) we obtain
(22)  
(23) 
The left side of (23) is much smaller than the right side and it can be ignored, thereby producing the solution
(24) 
which gives by insertion in (22).
Iv Linearized equation and eigenvalue problem
We want to ascertain the linear stability of the solution found in the previous section, and constant auxiliary fields. For this we need to know the dynamics and to linearize the corresponding SPEs about this solution. Let us assume that the dynamic SPEs are
(25) 
so that
(26) 
are the stationary SPEs (12)(14) and (18). Let be the flat solution and constant and of (22) and (24). Here is a bifurcation parameter, later to be identified as the temperature or the radius of a circular graphene membrane. Linear stability about follows from substituting in (25) and keeping terms of order in the result. We obtain the eigenvalue problem
(27) 
where , and are functionals acting upon . For appropriate values of , the real parts of all eigenvalues are negative and the flat constant solution is linearly stable. Then zero eigenvalues or eigenvalues with zero real part give the critical values at which nonflat solutions bifurcate from the flat constant solution. Stationary solutions bifurcate from at those critical values of for which . At such , (27) becomes
(28) 
While the dynamic SPEs are not known, (28) are the known linearized stationary SPEs about the flat constant solution. We consider (28) as an eigenvalue problem for the critical values and proceed to determine “eigenvalues” . The corresponding eigenvectors are buckling modes of the membrane issuing forth from the flat constant solution. When many of these modes become active we may have generated a variety of rippling states of the graphene sheet.
To solve the linearized stationary SPEs, we need appropriate boundary conditions. To solve the linearized plate equation (12), we need the value of in (13) for the constant and of (22) and (24). Now solves the equation with appropriate boundary conditions. A solution is
For a circular plate with zero data Dirichlet boundary conditions at the solution is . Then, by (13) and (14), , and (12) becomes
(29) 
There are different eigenvalue problems associated to solving Equation (29) with different boundary conditions for a finite graphene membrane. Each of these problems yield its critical temperature and size and allows us to know the shape of different eigenmodes that appear as buckled solutions of the suspended layer of graphene. We consider a finite circular layer of graphene with a free border (natural boundary conditions) or a graphene sheet clamped to a circular hole. Suspended graphene sheets are clamped at their boundaries but the case of a free border is easier to treat mathematically, so we start by analyzing it.
iv.1 Membrane with free border
The natural boundary conditions to solve (29) for a free circular graphene layer of radius are
(30) 
Let us first solve the eigenvalue problem for : with at . We find the solution
(31) 
where , , . The lowest possible value of in (31) corresponds to the first zero of the Bessel function . Note that . Thus there are two eigenvalues corresponding to azimuthal eigenfunctions () between the first two eigenvalues corresponding to radially symmetric eigenfunctions with .
The solution of that satisfies the other boundary condition is
(32) 
for (normalized so that ), and
(33) 
for .
We are interested in computing the critical temperature and radii. To this end, we use (31) with equations (22) and (24), and find that ripples appear below the critical temperature
(34) 
Note that the critical temperature decreases as the size decreases. As , the critical temperature tends to eV (about 186,000 K). For finite membrane size, (34) with instead of produces the following numerical estimate,
(35) 
Here is the temperature (measured in Kelvin) at which the mode appears and is measured in nanometers. The values considered for the parameters are eV (the lowest value considered in gui14 ()) and eV. At room temperature, equation (35) gives us different critical radii above which the different modes (characterised by ) appear, see table 2 and figure 1.
:  

(nm): 
iv.2 Clamped Membrane
In this more realistic case, the boundary conditions are and at instead of (30). We use instead of to distinguish the outofplane displacement in the clamped case. Following the same procedure as for equations (31) and (32), we get
(36) 
instead of (31), and in which is not yet determined. Now the radial part satisfies
(37) 
For , this equation is
in which we have used the boundary condition . Integrating once and using , we obtain
The vertical displacement is unbounded at unless
(38) 
Then
(39) 
Similarly, for the solution of is
(40) 
Clamped boundary conditions yield
(41)  
(42) 
The condition that be bounded at produces . Thus , i.e., and therefore
(43) 
Equation (40) becomes
(44)  
which agrees with (39) for .
According to (43), the critical radii are given by (35) with instead of . Then the critical radii of table 3 are greater than those in table 2, and buckling of a clamped graphene membrane should be observable only for radii over nm. Thus a clamped membrane remains flat for larger radii than in the case of natural free boundary conditions. This is according to our intuition that it is harder to buckle a clamped membrane than a membrane with a free border. The lowest possible value of in (31) corresponds to . Note that . Again there are two eigenvalues corresponding to azimuthal eigenfunctions () between the first two eigenvalues corresponding to radially symmetric eigenfunctions with . These first four modes given by (44) are depicted in figure 2.
:  

(nm):  43  66  80 
V Conclusions
In conclusion, the stationary saddlepoint equations give the vertical displacement of a graphene membrane and auxiliary fields associated to curvature and to charge fluctuations, provided elasticity and phononelectron interaction are considered gui14 (). We have solved these equations for a flat circular graphene sheet under constant auxiliary fields. We have then solved the SPEs linearized about the flat membrane solution and subject to either free or clamped boundary conditions. The latter problem has nonzero solutions only for discrete values of the auxiliary fields. These discrete values correspond to critical values of temperature and membrane radius at which buckling membrane solutions issue forth from the flat membrane. The flat membrane is stable only for temperatures larger than critical and radii smaller than critical. For an infinite membrane, the predicted critical temperature is extremely large and unphysical (186,000 K). However, critical radii at room temperature are in the nanometer range. Different modes of membrane buckling appear below specific critical temperatures and over specific critical radii, see tables 2 and 3. The stability of these modes is to be analytically determined, but experiments describing rippling mey07 () and buckling schoelz15 () of graphene show that the flat solution is not stable in most cases: Flat graphene has only been observed when it stands over a substrate lui09 ().
The calculated critical radii (tables 2 and 3) suggest that the different modes appearing in figures 1 and 2 could be observed for a specific temperature and membrane size. This brings the opportunity to develop new experiments, as no buckling should appear below the minimum critical radius . Similar experimental procedures to the ones in lin12 (); svensson11 (); xu14 (), where suspended graphene layers were forced to buckle using an external electrical potential or a STM, could be used to determine the value of and the validity of our results. For this purpose, suspended graphene layers may be clamped to holes of different dimensions, under different temperatures. Checking when the layer may or not buckle, it should be possible to get a better estimation of and by fitting the results of experiments with equation (34).
A different approach to study the effects due to the existence of nonflat modes would be experiments with graphene resonators bunch07 (); eichler11 (). In this setting, when the external driving force is sufficiently strong to deform the graphene layer close to one of the stable buckling modes (away from the flat configuration), the response of the resonator should be extremely nonlinear. Instead of oscillating around a single potential well (more or less parabolic) corresponding to the rippled configuration, the system may oscillate under influence of two different potential wells (the buckling modes). This effect could be measured by observing oscillation damping once the driving force is turned off.
Characterization of ripples in suspended graphene is not yet accurate enough to compare with our results. Once the experimental techniques have been improved, experiments with suspended graphene layers clamped to holes of different radii should show ripples that are combination of the possible bifurcating modes for these membranes (table 3).
The study of different solutions of the SPEs with positiondependent auxiliary fields and is left as future work. From a theoretical point of view, it is also interesting to take into account the influence of external forces, such as an electrostatic force, and to study the resulting dynamical effects. Then we may be able to reproduce the transition from the flat membrane with superimposed ripples to the buckled membrane and compare with existing STM experimental results schoelz15 ().