A comprehensive latticestability limit surface for graphene
Abstract
The limits of reversible deformation in graphene under various loadings are examined using latticedynamical stability analysis. This information is then used to construct a comprehensive latticestability limit surface for graphene, which provides an analytical description of incipient lattice instabilities of all kinds, for arbitrary deformations, parametrized in terms of symmetryinvariants of strain/stress. Symmetryinvariants allow obtaining an accurate parametrization with a minimal number of coefficients. Based on this limit surface, we deduce a general continuum criterion for the onset of all kinds of latticestabilities in graphene: an instability appears when the magnitude of the deviatoric strain reaches a critical value which depends upon the mean hydrostatic strain and the directionality of the principal deviatoric stretch with respect to reference lattice orientation. We also distinguish between the distinct regions of the limit surface that correspond to fundamentallydifferent mechanisms of lattice instabilities in graphene, such as structural versus material instabilities, and longwave (elastic) versus shortwave instabilities. Utility of this limit surface is demonstrated in assessment of incipient failures in defectfree graphene via its implementation in a continuum Finite Elements Analysis (FEA). The resulting scheme enables onthefly assessments of not only the macroscopic conditions (e.g., load; deflection) but also the microscopic conditions (e.g., local stress/strain; spatial location, temporal proximity, and nature of incipient lattice instability) at which an instability occurs in a defectfree graphene sheet subjected to an arbitrary loading condition.
keywords:
Graphene, ideal strength, latticestability limits, finite element analysis.1 Introduction
Limits to reversible deformation—the stress or strain at which elastictoinelastic transition takes place in a material — pose fundamental constraints on a material’s performance and determine its strength. The limiting stress or strain, in general, depends upon the loading path, and the dependence is often described by means of a phenomenological model termed a limiting (or failure) criterion. Mathematically, a limiting criterion is represented as a surface in stress or strain space, which separates the stable states of reversible deformation from the ‘failed’ or irreversiblydeformed states. For many materials, the limiting criterion is simple enough to be characterized by one or two material constants, which are readily determined from experiments. Examples include the Mises (or Tresca) yield criterion for metals, the MohrCoulomb criterion for cohesivefrictional solids (11), the DruckerPrager criterion for pressuredependent solids (15), and the HoekBrown criterion for rocks (22).
In crystalline materials that are free from lattice imperfections, the limit to elastic deformation sets an upper bound on the material strength (39). This upper bound, termed the ideal strength, depends on the intrinsic nature of bonding between atoms in the material. Because of the various types of lattice defects, such as dislocations, grain boundaries, interstitial impurities, and voids which normally exist in materials (40); (50), most conventional materials are irreversiblydeformed at stresslevels well below the ideal strength. However, in recent years, advancements in nanotechnology have enabled fabrication and growth of defectfree twodimensional crystals in which mechanical failure indeed occurs upon reaching stress levels near the ideal strength (8); (51).
Graphene, an atomic monolayer comprising a hexagonal network of covalentlybonded Catoms, is a representative example of such materials (see Fig. (2.5)). Experimental studies have shown that defectfree, singlecrystalline graphene can sustain nearidealstrength stresses while remaining within the reversible regime of deforequationmation (32); (49). Beyond the limit of elastic deformation, the fate of the material is determined by a strengthlimiting mechanism such as incipient plasticity or crackinitiation. Since, in the absence of latticeimperfections, a strengthlimiting mechanism can only be activated by a latticeinstability (33); (34), the incipient failure of a defectfree crystalline material is intrinsically related to the loss of internal latticestability. The point at which the loss of latticestability occurs is called the latticestability limit, and it varies from loading path to path. Further, the dependence of the latticestability limit on the loading condition in crystalline materials is generally too complex to be adequately characterized by one or two parameters.
Individually, latticeinstabilities of all kinds can be assessed via the latticedynamical stability analysis (see Born & Huang (7); Hill (21)), which asserts that a necessary and sufficient criterion for an ideal crystal under arbitrary uniform loading to be stable is that it exhibits stability with respect to bounded perturbations of all wavelengths. Integration of the latticedynamical stability analysis with a continuum analysis scheme such as FEA would be ideal for failure analysis of defectfree graphene crystals. This could enable assessment of incipient failures and ideal strength of graphene, under arbitrary loading conditions at realistic lengthscales and slow enough loading rates, directly in a continuumlevel simulation. However, stability analysis based on latticedynamics requires carrying out an elaborate sequence of computationallyexpensive steps that can not be treated within the confines of an analytical framework, making it difficult to integrate the latticedynamical stability analysis into a continuum scheme. Therefore, there remains a need for a general continuum criterion that could describe the onset of instability (of any kind) in graphene under an arbitrary state of deformation.
The aim of this work is to construct a comprehensive latticestability limit surface, which constitutes an analytical parametrization of incipient lattice instabilities of all kinds, over the space of all homogeneous deformations, in terms of stress/strain. There are two main difficulties in obtaining such a parametrization: First, crystalline materials are intrinsically anisotropic, so material response, including latticestability limit, varies with orientation (43). Secondly, two fundamentallydifferent types of latticeinstabilities govern strengthlimiting mechanisms under different loading conditions (36); (10): a longwave (or elastic) instability and a shortwave (or soft mode) instability. The condition for onset of an elastic instability can be parameterized in terms of strain via acoustic tensor analysis (see Kumar & Parks (28) for details), whereas the shortwave instabilities are much more complex since there is no continuum framework for parametrization of limiting conditions governing the onset of shortwave instabilities.
The proposed parametrization, in order to overcome the abovementioned difficulties, employs interpolation of the latticestability limits of graphene, corresponding to some representative homogeneous deformation modes, in the basis of symmetryinvariants of the strain/stress measure. Symmetryinvariants are special scalar functions of strain/stress that remain invariant under the point group operations. The use of symmetryinvariants ensures that the parametrization possesses the appropriate material anisotropy of graphene. It also introduces substantial functional simplification, reducing the requisite representative deformation modes needed for interpolation to a small set of biaxial deformations along the two special symmetrydirections in graphene: armchair and zigzag. The individual latticestability limits used in the interpolation are obtained by atomisticlevel latticedynamical stability analysis based on phonons to ensure that lattice instabilities of all kinds are captured in the analysis.
This paper is organized as follows. In Sec.(2), we briefly summarize the kinematics of graphene, outline the general framework of invariantbased representation theory of anisotropic materials, explicitly derive symmetryinvariants of the strain measure with respect to the symmetrygroup of graphene, and discuss the implementation of these ideas in the context of an (incipient) instability function of graphene. Employing a representation of the ideas outlined in Sec.(2), we systematically determine an analytical form for the limit surface for graphene in Sec.(3). Then, in Sec.(4), we map the instability function from strain space to stress space. The model is validated in Sec.(5): in Sec.(5.1), we assess the numerical accuracy of the representation of failure model based on interpolation; in Sec.(5.2), we compare the material instabilities predicted by the failure model with those from acoustic tensor analysis; and in Sec.(5.3), we compare the model predictions for buckling instabilities against the results obtained from buckling stability analysis based on the hyperelastic constitutive relation. FEA implementation of the failure function is discussed in Sec.(6): we discuss the procedure in Sec.(6.1), and illustrate its use in an example of limiting blistertype deformation in Sec.(6.2). Finally, we conclude in Sec.(7) by summarizing our results and their implications for analysis of strengthmeasuring nanoscale contact experiments.
2 Framework of representation
2.1 Kinematics
We consider graphene as a 2D deformable body denoted by unstressed reference configuration . Let and denote the coordinates of a material element of graphene in undeformed and deformed configuration, respectively. The convection of material points under deformation is described by a smooth, injective (onetoone) function called the motion. The nontranslational part of the motion can be equivalently defined by the positivedefinite secondorder deformation gradient tensor, . Then, the polar decomposition theorem provides the following factorizations of (19); (18); (16):
(1) 
where the orthogonal tensor characterizes rigidbody rotation, and (or ), termed the right (left) CauchyGreen tensor, characterizes shape and areachange. The deformation of a material point can be kinematically factored as the product of a purely dilatational (or shapepreserving, but areachanging) deformation , and a purely isochoric (or shapechanging, but areapreserving) deformation . Accordingly, the stretch tensor can be productdecomposed as
(2) 
where
(3) 
and
(4) 
here is the deviatoric stretch, is the 2D identity tensor, and and are the major and minor principal stretch directions, respectively.
Graphene is a composite of two Bravais lattices shifted with respect to each other by a shift vector . Upon deformation, the individual Bravais lattices follow the macroscopic deformation gradient , i.e., the lattice vectors in the deformed crystal are given as and . This is called CauchyBorn kinematics. However, the shift vector, separating the two lattices, does not obey the CauchyBorn kinematics, i.e., . The shift vector , in a deformed state, is determined by minimization of total energy of the deformed crystal. That is, under certain imposed deformations, graphene experiences sublattice shifts that lower the total energy of the crystal compared to the CauchyBorn unrelaxed crystal. These deformation modes are the ones that involve a deviatoric stretch, i.e., . Examples of this class of deformations include the uniaxial stress/stretch and simple shear. Conversely, because of symmetry, a purely volumetric deformation does not generate a sublattice shift.
The sublattice shift , associated with a deformation gradient is defined as the difference between the vectors connecting the two Bravais lattices in the relaxed and unrelaxed configurations. Mathematically,
(5) 
where denotes the shiftvector in the undeformed configuration.
The rigorous energetic and kinematic continuum description of a complex crystal, such as graphene, requires the definition of the strain energy density function as follows
(6) 
The above description implies the following variational form
(7) 
Localization of the principle of virtual work then requires the following energy balance
(8) 
where is the first PiolaKirchhoff stress (stress tensor workconjugate to ), and are the forces per unit reference volume acting on the atoms of the state defined by (, ). The conditions (8) imply the equilibrium equations
(9) 
and
(10) 
At equilibrium, the force on each atom in the lattice is zero, and the above equation becomes
(11) 
Thus equation (11) implicitly determines the equilibrium value of sublattice shift, , for which the crystal satisfies both external and internal equilibria. Inserting this form into the strain energy function allows representation of the equilibrium strain energy in terms of alone:
(12) 
2.2 Strain measure
In this work, the logarithmic strain measure is employed in the representation of the stressstrain constitutive response as well as of the limit surface. The spectral representation of is given by:
(13)  
where
(14) 
gives the areal logarithmic strain , and
(15) 
denotes the deviatoric part of . The orientation is measured from a fixed material axis aligned with the zigzag direction of the graphene lattice. In the rest of this article, we will refer to as the deviatoric stretch angle. Alternatively, we may also write the spectral representation of the logarithmic strain as
(16) 
where and are the major and minor principal logarithimic strains, respectively. In addition, we also define a mean hydrostatic strain , which characterizes the areal dilatation or contraction of the material.
2.3 Symmetryconstraints on the limit surface
The proposed approach is based on the notion of a limit function — a continuous scalarvalued function of the strain measure — such that all deformed configurations of the lattice that lie on the surface
(17) 
have reached a state of incipient lattice instability. All deformed states of the crystal lying in the interior of the limit surface satisfy , and are stable with respect to lattice perturbation of all wavelengths; conversely, for deformed states lying outside the surface, i.e., , the lattice is unstable with respect to an incremental perturbation of some wavelength.
For an anisotropic material, the limit function — in accordance with Neumann’s principle (43) — should include the material symmetry group of the underlying lattice, i.e.,
(18) 
where is an orthogonal tensor denoting the symmetry operations included in the material symmetry group ( in the case of graphene). A scalarvalued function of a tensor agency— such as — that remains invariant under a material symmetry group is called a invariant scalar function. The representation of a generic invariant scalar function involves using the isotropicization theorem and symmetryinvariants of the tensor agency with respect to the point symmetry group of the material.
2.4 Isotropicization theorem and symmetryinvariants
The isotropicization theorem — based on the notion of a materiallyembedded structure tensor — allows a invariant scalar function to be expressed in terms of a list of special scalar functions — , , …, — which are joint isotropic functions of and (4); (5); (31); i.e.,
(19) 
where
(20) 
Here denotes the transformation of the structure tensor under the orthogonal transformation . The functions are termed symmetryinvariants since they satisfy all the constraints belonging to the material symmetry group of the crystal. Smith (52); (53); (54) showed that the set of mutuallyindependent symmetryinvariants serves as a complete and irreducible basis for the representation of scalar constitutive functions of the anisotropic material. In the following section, we explicitly derive symmetryinvariants of for the structure tensor characterizing the material symmetry group of graphene.
2.5 Symmetryinvariants and functional basis for material symmetry group
The structure tensor characterizing the point group, the material symmetry group of graphene, is a sixthorder tensor given by (62); (63)
(21) 
where and are two traceless second order tensors given by
(22) 
and denote orthogonal material unit vectors fixed in the frame of the reference crystal such that at least one of them is aligned with an axis of reflection symmetry.
The complete and irreducible set of polynomial joint invariants of and constitute what is called an integrity basis, and are given by
(23) 
(24) 
where is the scalar tensor product,
and
(25)  
(26)  
(27) 
where indicates the orientation of maximum principal stretch.
The first two of these invariants, and , are simply two isotropic invariants of alone. Thus any material anisotropy in the constitutive response of graphene is captured solely by the third invariant . In a previous work (28), we showed that a constitutive function of graphene — such as its finite deformation hyperelastic response — is conveniently represented in terms of the integrity bases of the logarithmic strain with respect to the material symmetry group of graphene.
For the purpose of representation of , we employ an alternate set of symmetryinvariants of , comprising the mean hydrostatic strain , the magnitude of deviatoric strain ; and the symmetryreduced principal stretch direction defined as . Note that the elements of this set are not independent from the elements of the integrity basis, i.e., , , and ; and utilizing kinematical definitions of equations (1315), the two sets of symmetryinvariants can be shown to be related as follows:
(28) 
(29) 
(30) 
Thus, the failure function has the representation
(31) 
The set of symmetryinvariants, , and , which involve nonpolynomial combinations of strain components, is called a functional basis.
2.6 Features of the limit surface
The limit surface — as expressed in terms of the elements of the functional basis — prescribes the critical value of the magnitude of deviatoric strain as a function of the areal strain and the symmetryreduced principal stretch direction . Thus, the form of the limit function reads as
(32) 
such that if, at a given state of deformation, , then the material is stable with respect to incremental perturbations of all wavelengths; otherwise the instantaneous configuration of the lattice is unstable under some incremental perturbation. The limit surface, when graphically expressed in terms of , and , constitutes a 3D polar surface with on the radial axis, on the vertical axis and as the angular coordinate.
Boundedness of the limit surface — The proposed limit surface is bounded both on the vertical axis (i.e., axis) and the radial axis (i.e., axis), as explained in the following.
On the vertical axis:

By definition, the magnitude of deviatoric strain is a positive semidefinite entity, i.e., at all values of . Thus, the limit surface is bounded from below by the surface .

From above, the failure surface is bounded by the critical value of the magnitude of deviatoric strain beyond which the crystal is unstable, i.e., .
On the radial axis:

Because of graphene’s limiting thickness and extreme bending compliance, essentially any attempt to decrease the area of a sufficiently large, flat graphene sheet would result in outofplane buckling, graphene’s small, though nonzero bending stiffness notwithstanding; for this reason, we consider that the domain of the surface is bounded from below: , providing a graphical representation using as a ‘radial’ coordinate. Moreover, if is held at , imposition of a nonzero creates a negative principal stress, again leading to buckling in a sufficiently large, flat graphene sheet; hence we take for all . (These limits follow from the quadratic phonon dispersion relations of ZA phonons at in a periodicallyextended planar crystal.)

Upon a continued equibiaxial deformation, a state is reached, when graphene fails without any shear strain, i.e., in pure dilatation. Therefore, the radial domain of the surface is bounded from above as well: and for all .
3 Form of the limit surface in terms of the functional basis
Employing the representation ideas outlined in the previous section, we now systematically determine the functional form of the limit function for graphene in terms of the three symmetryinvariants. For this purpose, we first obtain the set of data containing the stabilitylimiting values over a range of values, calculated from latticedynamical stability analysis, for a set of representative homogeneous deformation paths.
3.1 Lattice dynamical stability analysis
The atoms in a crystalline material remain in a state of oscillatory motion about their equilibrium sites, constituting what is called the random lattice vibrations. Phonons are the normal modes of the lattice vibrations. Fundamentally, each phonon is a collective oscillation of atoms characterized by a wavevector and a oscillation frequency , where labels the branch index. The relation between and is called the dispersion relation of the crystal. The phonondispersion relation of a crystal contains a total of branches where is the number of atoms in the unit cell ( for graphene).
Three (3) of these branches are called acoustic and the remaining are called optical.
Phonons constitute a useful means of assessing the mechanical stability of a crystal lattice. Mechanical stability of a crystal requires that the underlying lattice remains stable against spatial perturbations of all wavelengths, and this is ensured if the following two conditions hold:

First, all the acoustic branches should have a positive slope at , i.e.,
(33) If at some point along the loading path, this criterion ceases to hold then it is an indication of an incipient instability, called a long wavelength instability, in the material. In monolayer graphene, a long wavelength instability may correspond to either an outofplane instability (buckling) or an elastic instability (34); (27). Buckling is a mode of structural instability wherein graphene explores the third spatial dimension via deformation under compression. Unlike a material instability, buckling does not involve material fracturing or undergoing plastic deformation. Further, since graphene is one atomic layer thin, a buckling instability manifests itself as a longwavelength instability in outofplane acoustic phonon (also called ZA branch) mode of graphene.
On the other hand, an elastic instability is a mode of material failure and manifests itself as a longwavelength instability in the inplane acoustic phonon (LA or TA branch) modes of graphene. An elastic instability in the material can also be captured by acoustic tensor analysis; therefore the onset of elastic instabilities in material can be parameterized with strain via the acoustic tensor route also. 
Second, all phonon frequencies for should be real and nonzero, i.e.,
(34) If at some point along the loading path, this criterion ceases to hold then it is an indication of an incipient soft mode instability in the material. Unlike an elastic instability, a soft mode instability does not have an equivalent continuum criterion, and therefore a mathematical framework parametrizing the onset of such instabilities in terms of stress or strain.
Latticestability analysis based on the above criteria is most effectively assessed by phonons, which, below the Debye temperature ( for graphene), constitute a complete normal basis for lattice vibrations. A deformationinduced instability of any kind — macroscopic or microscopic — is directly visible in the dispersion of phonons.
3.2 Basics of density functional perturbation theory
DFPT employs density functional theory (DFT) in conjunction with linear response perturbation theory to calculate the lattice dynamical properties of a crystal (see references (2); (13); (42)). In density functional theory (DFT) (23); (26), the ground state energy of a system of interacting electrons moving under the influence of an external potential is obtained by minimization of a universal functional, , of electronic density . The electronic density that minimizes the energy functional also happens to be the true ground state electronic density of the system.
Kohn and Sham (26) showed that the original system of interacting electrons can be replaced by an auxiliary system of noninteracting electrons moving under the influence of an effective potential such that the auxiliary system at the ground state possesses the same electronic density as the original system. The motion of electrons in such an auxiliary system are described by a set of oneelectron equations (also called KohnSham equations):
(35) 
where is a KohnSham orbital and is the corresponding eigenvalue. The electronic density is calculated from the KohnSham orbitals using the relation:
(36) 
where is the occupation function, and is the Fermi energy. The effective potential is a functional of electronic density, and is given by
(37) 
where the first term on the right, , is the potential due to interaction between electrons and ionic cores, the second term is the potential due to Coloumb interaction between electrons, and the last term, , is the functional derivative of the exchangecorrelation energy functional . Once an explicit form for has been determined, the set of equations (35  37) can be solved selfconsistently to determine the ground state energy and electronic density of the system.
Calculation of phonons requires knowledge of interatomic force constants, defined as the derivatives of the ground state total energy of the system with respect to ionic coordinates . Using FeynmanHellman theorem, we evaluate such derivatives as:
(38) 
where is the energy corresponding to Coloumbic interaction between ionic cores.
Thus, the evaluation of the interatomic force constants requires ground state electron density as well as the firstorder correction .
The DFPT method calculates the firstorder correction in electronic density from lattice’s response to a set of monochromatic lattice perturbations, each characterized by a vector. Provided the amplitude of such perturbations are small enough, the application of perturbation theory allows to recover a set of selfconsistent relations for the firstorder corrections in electronic density and wavefunctions. Let and be the unperturbed wave function and eigenenergy of an electron with wavevector and bandindex , respectively. Firstorder correction in the wave function then satisfies the following equation:
(39) 
where is the projector operator over the manifold of states of wavevector , and is the associated perturbation in the effective potential due to perturbation, given by:
(40) 
Perturbation in the electronic density is obtained as:
(41) 
being the volume of the system. Upon solving the set of equations (3941) selfconsistently, we obtain the firstorder corrections in wave function and electronic density of the system. Knowledge of such corrections allows direct access to the dynamical matrix of the system, defined as the Fourier transform of the interatomic force constant matrix, i.e.,
(42) 
where denotes the mass of atom. Phonon frequency and the associated eigenmode are then determined by solving the eigenvalue equation:
(43) 
3.3 Details of DFPT calculations
The exchange correlation energy of electrons is treated with the Local Density Approximation (LDA) of Perdew and Wang ((46)). The interaction between ionic cores and valence electrons is represented by an ultrasoft pseudopotential (57). KohnSham wave functions are represented using a planewave basis with an energy cutoff of 30 Ry and a charge density cutoff of 300 Ry. Integration over the irreducible Brillouin zone (IBZ) is performed with a uniform mesh of points, and occupation numbers are smeared using the MarzariVanderbilt cold smearing scheme (37) with broadening of 0.03 Ry. Errors in the Cauchy stresses and total energy due to basisset size, smearing parameter, and points are converged to less than 0.034 N/m and 0.01 Ry, respectively.
The phonon dispersion relations — of the undeformed and deformed graphene — are computed via linear response calculations as implemented in the density functional perturbation theory (DFPT) (2); (13). The dynamical matrix is calculated on an uniform grid of points in the IBZ — which is then fastFouriertransformed to calculate the interatomic force constants (IFC), corrected by the acoustic sum rule to ensure that for all the acoustic branches. The IFC’s are then used to interpolate phonon frequencies over a dense set of points in the IBZ. Both energies and phonon dispersion relations in this work are performed on a twoatom primitive unit cell of graphene shown in Fig. (2.5). All calculations are fully relaxed in terms of shift vector to ensure conditions of zero atomic force.
The LA and TA phonon branches, in the neighborhood of , obey linear dispersion: and ; where and and are longitudinal and transverse acoustic wavevelocities, respectively.
As a simple verification of the DFPT phonon calculations, Tab. (3.3) compares longitudinal and transverse acoustic wavevelocities obtained from the phonon dispersion with experimentallymeasured values; also, since the G band in the Raman spectra of graphene is associated with the doublydegenerate
LO  TO vibrational modes at , we also compare calculated LO/TO frequencies at with the Gband frequency measured in Raman spectroscopy.
(km/s)
(km/s)
Raman peak
Measured value (58)
1582
Continuum model (LDA) (28)
21.80
13.78
—
DFPT Phonons in this work (LDA)
21.15
13.50
1540
DFPT allows calculation of phonon dispersion relations with high accuracy at any arbitrarilydeformed state of the graphene lattice, whereas the errors in phonon frequencies based on empirical potentials can be forbiddingly large, rendering the predictions somewhat unreliable. For example, the sound velocities from an empirical manybody potential in the unstrained state of graphene differed from measured values (58) by nearly 50% (12); in contrast, the longitudinal and transverse elastic wavevelocities obtained from DFPT phonon calculations differ by less than 1% from measured values, as shown in Tab. (3.3).
3.4 Sampling scheme
For the purpose of sampling deformed states of the lattice, we use the deviatoric stretch angle to designate what we term a ‘deformation path’; i.e., increasing at a fixed is referred to as moving along the deformation path prescribed by . To determine the point of incipient lattice instability along a deformation path:

The lattice is first deformed via a pure equibiaxial strain, i.e., , followed by an isochoric shapechanging stretch of the form where and ; is kept fixed. Thus, total strain in the final configuration is , as schematically shown in Fig. (3.5a).

For each deformation pair along the sampled deformation path , we check the phonons to determine if the conditions for acoustic stability (given by equation (33)) and shortwave stability (given by equation (34)) hold. If at some critical magnitude of deviatoric strain, , either of the two conditions ceases to hold then the lattice is considered as incipiently unstable, and the corresponding triplet corresponds to a point on the instability surface.
Our sampling set includes values of ranging from the undeformed state to the critical equibiaxial deformation, , at which the lattice reaches a soft mode instability in the absence of deviatoric deformation. For each sampled value of , and along each sampled deviatoric stretch angle , the superimposed deviatoric strain is varied such that its magnitude varies over the range , the upper limit being the critical value of the stretch at the given value of and . Owing to the symmetry of graphene, the deviatoric stretch angle needs to be sampled only over the range . Further, we consider only deformation paths in the BZ, corresponding to the zigzag direction , and the armchair direction . Later, we will show that only two sampling directions suffice for an accurate angular interpolation of the failure function.
3.5 Representation and interpretation of the data
From the latticedynamical stability analysis, we obtain the limiting values, each in the form a triplet , representing the critical value of the deviatoric strain magnitude — — as a function of mean hydrostatic strain for both the sampled deformation paths: (shown in Fig. (3.5b)) and (shown in Fig. (3.5c)).
Shown in Fig. (3.5) are the phonon dispersion relations of graphene calculated with DFPT at certain critical deformed states, indicated as ‘1’, ‘2’ and ‘3’, along the sampled deformation path ; these particular critical states of deformation are those labeled on Fig. (3.5b). Depending upon the level of mean hydrostatic strain, a critical phonon state might be associated with buckling, as in ‘1’, for which the longwave ZA phonon slope at vanishes; with an elastic (material) instability, as in ‘2’, for which the longwave LA slope at vanishes; or with a shortwave material instability as in ‘3’, where the frequency of the TA phonon vanishes at . Thus, at different states of deformation, different kinds of mechanisms govern lattice instability in a graphene sheet. The proposed limit surface is a smoothed envelope of all the possible latticeinstabilities, as we will elaborate further in Sec.(3.7).
3.6 From discrete dataset to continuum latticestability limit function via interpolation
Construction of a continuumlevel limit function involves interpolating the discrete dataset of the preceding section in a manner that is consistent with the material symmetry () of graphene. This task is accomplished by interpolating the critical magnitude of deviatoric strain in terms of two of the symmetryinvariants of the strain tensor: and . The interpolation is performed in two steps: first, interpolation in terms of , called radial interpolation, and second, interpolation in terms of , called angular interpolation.
Radial interpolation (interpolation in terms of )
Along each sampled deformation path , we express by a continuous function — called a radial interpolation function. Each such function is obtained by fitting the sampled datapoints (plotted in Fig. (44a)), which are in the form of a doublet and denote the limiting strains along a sampled deformation path, to a polynomial function in of the form:
(44)  
where is the mean hydrostatic strain at which an incipient lattice instability emerges without any deviatoric strain.
The radial interpolation functions thus obtained, corresponding to the two sampled deformation paths, are shown in Fig. (44b).
Angular interpolation
The angular interpolation scheme involves approximating the failure surface by a weighted sum of the radial interpolation functions — with the weight functions being the angular shape functions along the sampled deformation paths. The following representation for the angular shape function conveniently satisfies the requirements that failure function is smooth, satisfies all the pointgroup symmetries of the lattice, and possesses periodicity with respect to :
(45) 
where the number of terms in the expansion is given by the number of deformation paths () being sampled, and the corresponding coefficients are determined using the following condition:
(46) 
Since the angular variation of the limit function is quite smooth, sampling along only two deformation paths, i.e., , suffices for accurate interpolation of the limit surface (see Fig. (5.1) in Sec.(5)). The corresponding two shapefunctions are given as:
(47)  
(48) 
Following the above interpolation scheme, is generically expressed as:
(49) 
where is the total number of deformation paths. Substituting from equation (44) and equation (45) into equation (49), we obtain the generic form of the lattice stability function as:
(50) 
The various coefficients appearing in the expression of the failure function have been tabulated in the Tab. (2).
1  0.21561   0.0403958   0.327806  7.51334   24.8089  

0.0251542   0.583471  4.55324  16.0798  26.6759   20.8193  6.22982 
The quadrant of the instability surface resulting from the interpolation scheme is shown in Fig. (2b).
3.7 Limit surface as a smoothed envelope of various possible instabilities
Different kinds of strengthlimiting mechanisms lead to instability of the graphene lattice under different modes of deformation. The limit surface is a smoothed representation of the envelope of all possible lattice instabilities: long wavelength as well as short wavelength, and structural as well as material failures.
Graphene, being a strictly twoD structure, has an exceedingly small bending rigidity. Therefore, when subjected to an inplane deformation, a graphene sheet can be brought to a state of incipient instability via one of two mechanisms:

a material instability in tension, which arises when the major principal strain (maximum strain over all material directions) reaches a limiting value, while the minor principal strain (minimum strain over all material directions) has not reached the critical value in compression so that both the Cauchy principal stress components are nonnegative, i.e., . This material instability can be either of the two types: longwave (also called elastic) or a shortwave (also called soft mode) or

a structural instability via outofplane buckling, which arises when the minor principal strain (minimum strain over all material directions) in the material reaches the critical value in compression such that one of the principal stress components or becomes incipiently compressive (), while the major principal strain has still not reached the limiting tensile strain.
The two principal strains, in terms of and , are given by and . Thus, at a fixed equibiaxial strain , with a superposed, increasing deviatoric strain , the major principal strain monotonically increases while the minor principal strain monotonically decreases. Now, whether a deviatoric strain superposed on an equibiaxial strain results in a material instability or a structural instability depends on the magnitude of , as schematically illustrated in Fig. (3.7). For example, when is small, the minor value reaches the critical value in compression, which is , before the major principal value reaches the limiting tensile strain. Thus, at small equibiaxial strain, the magnitude of the deviatoric strain is essentially limited by a buckling instability. This is also shown in the phonon dispersion relation of Fig. (3.51). Interestingly, when , the membrane is unstable in buckling as soon as any deviatoric strain is applied; thus the limiting value of .
On the other hand, when is large, reaches the limiting tensile strain before reaches the critical value in compression, and the failure is invariably a material failure (fracture). Further, as approaches closer to the critical value , the material starts to fail essentially in dilatation via a shortwave instability, as shown in Fig. (3.53). Thus, .
4 Latticestability limit surface in stress space
Often, it is useful to describe the latticestability limits in terms of stress. In the present case, although the limiting conditions are explicitly obtained in terms of the invariants of strain, the relationship can be numerically mapped to its counterpart in the stress space via a constitutive function appropriately describing the stressstrain response in graphene over the entire range of deformation up to the elastic stability limit. In a previous work (28), we derived a largedeformation hyperelastic constitutive response function for graphene based on ab initio calculations — which we briefly describe in the following.
The constitutive response function is derived within the framework of nonlinear hyperelasticity and invariantbased representation theory. Our formulation employs the symmetryinvariants of the log strain measure — — to express the dependence of the free energy per unit reference area on the state of strain , i.e., (see Kumar & Parks (28) for details). The strain energy density function accepts the additive decomposition: .
The term — the energetic response under pure dilation — is described by a function based on the universal binding energy relation (48); (47):
(51) 
The term — the energetic response under an isochoric deformation — is given by a simple additive form:
(52) 
where is a second order elastic constant, and is a thirdorder elastic constant. The coefficients are determined from fits to ab initio energies for a set of deformed configurations (see Tab.[3]).
(J/m)  (N/m)  (N/m)  (N/m)  (N/m)  

LDA  116.43  1.38  172.18  27.03  5.16  93.17  4408.76 
Differentiation of with respect to gives the workconjugate stress tensor, i.e., .
(53) 
and the Cauchy stress is given as (see Ogden (44))
(54) 
where denotes the right CauchyGreen stretch tensor. Employing equation (53) and equation (54), we numerically map the contours in the strain plane to stressplane. Once the Cauchy stress has been determined, we can obtain the principal stresses from the spectral representation of :
(55) 
where , and are the principal Cauchy stress components; and and are the corresponding principal directions. The orientation is measured from a fixed material axis aligned with the zigzag direction of the graphene lattice.
Employing the above transformation, we obtain the latticestability plot in stress space as a functional relation between the components of the functional bases comprising the maximum shear stress , the mean hydrostatic Cauchy stress , and the symmetryreduced deviatoric stress direction . The limit function— in terms of , , and — can be expressed in a form similar to the limit function in strain space, i.e.,
(56) 
such that if, at a given state of deformation, , then the material is stable with respect to incremental perturbations of all wavelengths; otherwise the instantaneous configuration of the lattice is unstable under some incremental perturbation.
The surfaces described by equation (56) constitute the latticestability limit surface in the space of symmetryinvariants of the Cauchy stress tensor, as shown in Fig. (56b).
5 Validation
The limit surface — given by equation (50) — contains information of all possible latticeinstabilities that can be triggered by a homogeneous deformation. Specifically, once this failure surface is determined, we can assess the stability of graphene subject to a biaxial state of strain — — for any biaxiality (the to ratio) and any deviatoric stretch angle . For such deformations, the principal strains are directly given by the strain components and as and . Substituting a particular value of in the limit surface, we can obtain the 2D radial section of the limit surface corresponding to that particular value of . For example, Fig. (5) shows the 2D radial sections of the limit surface corresponding to and , denoted by and , respectively. These contours characterize the stability of a graphene sheet subjected to a homogeneous biaxial state of strain referred to the set of axes defined by and , i.e., are aligned with (see Fig. (3.5a)).
On this contour, we have indicated points corresponding to certain special homogenous deformation cases such as uniaxial strain and uniaxial stress along the armchair and zigzag directions :

The point indicated as , given by the intersection of the line with the contour , denotes the uniaxial strain along the zigzag direction.

The point indicated as , given by the intersection of the line with the contour