A comprehensive lattice-stability limit surface for graphene

# A comprehensive lattice-stability limit surface for graphene

## Abstract

The limits of reversible deformation in graphene under various loadings are examined using lattice-dynamical stability analysis. This information is then used to construct a comprehensive lattice-stability limit surface for graphene, which provides an analytical description of incipient lattice instabilities of all kinds, for arbitrary deformations, parametrized in terms of symmetry-invariants of strain/stress. Symmetry-invariants 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 lattice-stabilities 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 fundamentally-different mechanisms of lattice instabilities in graphene, such as structural versus material instabilities, and long-wave (elastic) versus short-wave instabilities. Utility of this limit surface is demonstrated in assessment of incipient failures in defect-free graphene via its implementation in a continuum Finite Elements Analysis (FEA). The resulting scheme enables on-the-fly 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 defect-free graphene sheet subjected to an arbitrary loading condition.

###### keywords:
Graphene, ideal strength, lattice-stability limits, finite element analysis.
1

## 1 Introduction

Limits to reversible deformation—the stress or strain at which elastic-to-inelastic 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 irreversibly-deformed 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 Mohr-Coulomb criterion for cohesive-frictional solids (11), the Drucker-Prager criterion for pressure-dependent solids (15), and the Hoek-Brown 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 irreversibly-deformed at stress-levels well below the ideal strength. However, in recent years, advancements in nanotechnology have enabled fabrication and growth of defect-free two-dimensional 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 covalently-bonded C-atoms, is a representative example of such materials (see Fig. (2.5)). Experimental studies have shown that defect-free, single-crystalline graphene can sustain near-ideal-strength 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 strength-limiting mechanism such as incipient plasticity or crack-initiation. Since, in the absence of lattice-imperfections, a strength-limiting mechanism can only be activated by a lattice-instability (33); (34), the incipient failure of a defect-free crystalline material is intrinsically related to the loss of internal lattice-stability. The point at which the loss of lattice-stability occurs is called the lattice-stability limit, and it varies from loading path to path. Further, the dependence of the lattice-stability limit on the loading condition in crystalline materials is generally too complex to be adequately characterized by one or two parameters.

Individually, lattice-instabilities of all kinds can be assessed via the lattice-dynamical 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 lattice-dynamical stability analysis with a continuum analysis scheme such as FEA would be ideal for failure analysis of defect-free graphene crystals. This could enable assessment of incipient failures and ideal strength of graphene, under arbitrary loading conditions at realistic length-scales and slow enough loading rates, directly in a continuum-level simulation. However, stability analysis based on lattice-dynamics requires carrying out an elaborate sequence of computationally-expensive steps that can not be treated within the confines of an analytical framework, making it difficult to integrate the lattice-dynamical 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 lattice-stability 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 lattice-stability limit, varies with orientation (43). Secondly, two fundamentally-different types of lattice-instabilities govern strength-limiting mechanisms under different loading conditions (36); (10): a long-wave (or elastic) instability and a short-wave (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 short-wave instabilities are much more complex since there is no continuum framework for parametrization of limiting conditions governing the onset of short-wave instabilities.

The proposed parametrization, in order to overcome the above-mentioned difficulties, employs interpolation of the lattice-stability limits of graphene, corresponding to some representative homogeneous deformation modes, in the basis of symmetry-invariants of the strain/stress measure. Symmetry-invariants are special scalar functions of strain/stress that remain invariant under the point group operations. The use of symmetry-invariants 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 symmetry-directions in graphene: armchair and zigzag. The individual lattice-stability limits used in the interpolation are obtained by atomistic-level lattice-dynamical 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 invariant-based representation theory of anisotropic materials, explicitly derive symmetry-invariants of the strain measure with respect to the symmetry-group 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 blister-type deformation in Sec.(6.2). Finally, we conclude in Sec.(7) by summarizing our results and their implications for analysis of strength-measuring nano-scale 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 (one-to-one) function called the motion. The non-translational part of the motion can be equivalently defined by the positive-definite second-order deformation gradient tensor, . Then, the polar decomposition theorem provides the following factorizations of (19); (18); (16):

 F=RU=VR, (1)

where the orthogonal tensor characterizes rigid-body rotation, and (or ), termed the right (left) Cauchy-Green tensor, characterizes shape- and area-change. The deformation of a material point can be kinematically factored as the product of a purely dilatational (or shape-preserving, but area-changing) deformation , and a purely isochoric (or shape-changing, but area-preserving) deformation . Accordingly, the stretch tensor can be product-decomposed as

 U=Ua~U=~UUa, (2)

where

 Ua≡J1/2I (3)

and

 ~U≡λr1⊗r1+λ−1r2⊗r2; (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 Cauchy-Born kinematics. However, the shift vector, separating the two lattices, does not obey the Cauchy-Born 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 sub-lattice shifts that lower the total energy of the crystal compared to the Cauchy-Born 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 sub-lattice shift.
The sub-lattice 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,

 s=d−dCB;where dCB=Fd0, (5)

where denotes the shift-vector 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

 ψ=ˇψ(F,s). (6)

The above description implies the following variational form

 δψ=∂ψ∂F:δF+∂ψ∂s:δs. (7)

Localization of the principle of virtual work then requires the following energy balance

 0=⎧⎪⎩∂ψ∂F−T(1)⎫⎪⎭:δF+⎧⎪⎩∂ψ∂s−f⎫⎪⎭:δs, (8)

where is the first Piola-Kirchhoff stress (stress tensor work-conjugate 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

 ∂ψ∂F=T(1), (9)

and

 ∂ψ∂s=f. (10)

At equilibrium, the force on each atom in the lattice is zero, and the above equation becomes

 ∂ψ∂s=0. (11)

Thus equation (11) implicitly determines the equilibrium value of sub-lattice 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:

 ψ=ˇψ(F,seq(F))≡^ψ(F). (12)

### 2.2 Strain measure

In this work, the logarithmic strain measure is employed in the representation of the stress-strain constitutive response as well as of the limit surface. The spectral representation of is given by:

 E(0) = 12lnJIlnUa+lnλ(r1⊗r1−r2⊗r2)ln~U (13) ≡ 12ϵaI+E(0)0,

where

 ϵa=trE(0)=lnJ=ln(detU), (14)

gives the areal logarithmic strain , and

 E(0)0=ln~U=lnλ(r1⊗r1−r2⊗r2), (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

 E(0)=Emaxr1⊗r1+Eminr2⊗r2. (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 Symmetry-constraints on the limit surface

The proposed approach is based on the notion of a limit function — a continuous scalar-valued function of the strain measure  — such that all deformed configurations of the lattice that lie on the surface

 F(E(0))=0, (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.,

 F(QTE(0)Q)=F(E(0)) ∀ Q∈ G, (18)

where is an orthogonal tensor denoting the symmetry operations included in the material symmetry group ( in the case of graphene). A scalar-valued 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 symmetry-invariants of the tensor agency with respect to the point symmetry group of the material.

### 2.4 Isotropicization theorem and symmetry-invariants

The isotropicization theorem — based on the notion of a materially-embedded 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.,

 F(E(0))=^F(J1,J2,...,Jn), (19)

where

 Ji(E(0);H)=Ji(QTE(0)Q;PQ(H))∀ Q∈SO2. (20)

Here denotes the transformation of the structure tensor under the orthogonal transformation . The functions are termed symmetry-invariants since they satisfy all the constraints belonging to the material symmetry group of the crystal. Smith (52); (53); (54) showed that the set of mutually-independent symmetry-invariants 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 symmetry-invariants of for the structure tensor characterizing the material symmetry group of graphene.

### 2.5 Symmetry-invariants and functional basis for C6v material symmetry group

Figure 1: (a) Graphene lattice with orientations of the material unit vectors — and — and the Cartesian unit vectors — and — indicated. The dashed blue lines denote the unit cell used in the ab initio calculations. The GGA (LDA) value of the lattice parameter is also indicated. The armchair and zigzag directions are along the and axes, respectively. (b) Brillouin zone of graphene with high symmetry points and irreducible wedge indicated.

The structure tensor characterizing the point group, the material symmetry group of graphene, is a sixth-order tensor given by (62); (63)

 H=M⊗M⊗M−(M⊗N⊗N+N⊗M⊗N+N⊗N⊗M), (21)

where and are two traceless second order tensors given by

 M=^x⊗^x−^y⊗^y;N=^x⊗^y+^y⊗^x. (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

 J1≡ϵa=trE(0)=lnJ, (23)
 J2≡(γi/2)2=12 E(0)0:E(0)0=(lnλ)2, (24)

where is the scalar tensor product, and

 J3≡(γθ/2)3 =18H[E(0)0,E(0)0,E(0)0] (25) =18[(M:E(0)0)3−3(M:E(0)0)(N:E(0)0)2] (26) =(lnλ)3cos6θ, (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 symmetry-invariants of , comprising the mean hydrostatic strain , the magnitude of deviatoric strain ; and the symmetry-reduced 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 (13-15), the two sets of symmetry-invariants can be shown to be related as follows:

 1. Mean hydrostatic strain ¯E =ϵa2=lnJ1/2. (28)
 2. Magnitude of deviatoric strain γ =√E(0):E(0)2=γi2=|lnλ|≥0. (29)
 3. Directionality Θ(θ) =arccos[γ3θγ3i]=arccos(cos6θ). (30)

Thus, the failure function has the representation

 F(E(0))=^F(γ,¯E,Θ(θ)). (31)

The set of symmetry-invariants, , and , which involve non-polynomial 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 symmetry-reduced principal stretch direction . Thus, the form of the limit function reads as

 F(E(0))≡γc(¯E,Θ(θ))−γ; γ≥0, (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 semi-definite 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., .

• 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 out-of-plane buckling, graphene’s small, though non-zero 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 non-zero 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 periodically-extended planar crystal.)

• Upon a continued equi-biaxial 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 symmetry-invariants. For this purpose, we first obtain the set of data containing the stability-limiting values over a range of values, calculated from lattice-dynamical 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 wave-vector and a oscillation frequency , where labels the branch index. The relation between and is called the dispersion relation of the crystal. The phonon-dispersion 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.,

 dωnd|q|∣∣∣q=0>0 ∀ n∈1,...,3. (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 out-of-plane 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 long-wavelength instability in out-of-plane 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 long-wavelength instability in the in-plane 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 non-zero, i.e.,

 ω2n(q≠0)>0 ∀ n∈1,...,3N. (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.

Lattice-stability 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 deformation-induced 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 non-interacting 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 one-electron equations (also called Kohn-Sham equations):

 [−ℏ22m∇2+VEff(r)]ψn(r)=ϵnψn(r), (35)

where is a Kohn-Sham orbital and is the corresponding eigenvalue. The electronic density is calculated from the Kohn-Sham orbitals using the relation:

 n(r)=∑n|ψn(r)|2g(ϵF−ϵi), (36)

where is the occupation function, and is the Fermi energy. The effective potential is a functional of electronic density, and is given by

 VEff(r)=VIE(r)+e2∫n(r′)|r−r′|dr′+vxc(r), (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 exchange-correlation energy functional . Once an explicit form for has been determined, the set of equations (35 - 37) can be solved self-consistently 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 Feynman-Hellman theorem, we evaluate such derivatives as:

 ∂2E(R)∂RI∂RJ=∫[∂n(r)∂RJ]n0(r)∂VIE(r)∂RIdr+∫n0(r)∂2VIE(r)∂RI∂RJdr+∂2VII(R)∂RI∂RJ, (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 first-order correction .
The DFPT method calculates the first-order 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 self-consistent relations for the first-order corrections in electronic density and wavefunctions. Let and be the unperturbed wave function and eigenenergy of an electron with wave-vector and band-index , respectively. First-order correction in the wave function then satisfies the following equation:

 [−ℏ22m∇2+VEff]Δψk+qv=ϵkvΔψk+qv−Λk+qΔVEff(r)ψkv, (39)

where is the projector operator over the manifold of states of wave-vector , and is the associated perturbation in the effective potential due to perturbation, given by:

 ΔVEff(r)=∫Δn(r′)|r−r′|d3r′+Δn(r)[dvxc(r)dn]n0(r)+ΔVIE(r). (40)

Perturbation in the electronic density is obtained as:

 Δn=4V∑kψkve−iq.rΛk+qΔψk+qv, (41)

being the volume of the system. Upon solving the set of equations (39-41) self-consistently, we obtain the first-order 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.,

 DIaJb(q)=1√MIMJ∑RCIaJb(R)e−iq.R; CIaJb=∂2E(R)∂RIa∂RJb, (42)

where denotes the mass of atom. Phonon frequency and the associated eigenmode are then determined by solving the eigenvalue equation:

 ω2(q)uq=D(q)uq. (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). Kohn-Sham wave functions are represented using a plane-wave 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 Marzari-Vanderbilt cold smearing scheme (37) with broadening of 0.03 Ry. Errors in the Cauchy stresses and total energy due to basis-set 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 fast-Fourier-transformed 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 two-atom 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 wave-velocities, respectively. As a simple verification of the DFPT phonon calculations, Tab. (3.3) compares longitudinal and transverse acoustic wave-velocities obtained from the phonon dispersion with experimentally-measured values; also, since the G band in the Raman spectra of graphene is associated with the doubly-degenerate LO - TO vibrational modes at , we also compare calculated LO/TO frequencies at with the G-band 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 Table 1: In-plane longitudinal and transverse wave velocities calculated from phonon dispersion closely agree with measured values as well as with values calculated from the continuum model (28). The LO/TO frequencies at are also in good agreement with the G-peak of the Raman spectra of graphene. Note: is the mass per reference area of graphene which is calculated from the atomic mass of carbon and the undeformed lattice constant of graphene; denotes the elastic modulus in uniaxial strain; and denotes the shear modulus.

DFPT allows calculation of phonon dispersion relations with high accuracy at any arbitrarily-deformed 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 many-body potential in the unstrained state of graphene differed from measured values (58) by nearly 50% (12); in contrast, the longitudinal and transverse elastic wave-velocities 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 equi-biaxial strain, i.e., , followed by an isochoric shape-changing stretch of the form where and ; is kept fixed. Thus, total strain in the final configuration is , as schematically shown in Fig. (3.5-a).

• 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 short-wave 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 equi-biaxial 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 lattice-dynamical 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.5-b)) and (shown in Fig. (3.5-c)).

Figure 2: (a) Notation scheme illustrating deformation of the graphene lattice: the undeformed lattice (faint black) is first subjected to a uniform dilatation (resulting intermediate configuration shown in faint blue), and then to an isochoric shape-changing deformation (resulting final configuration shown in dark red). (b) Cartesian plot of as a function of for the sampled orientation , which corresponds to the zigzag direction. (c) Cartesian plot of as a function of for the sampled orientation , which corresponds to the armchair direction.

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.5-b). Depending upon the level of mean hydrostatic strain, a critical phonon state might be associated with buckling, as in ‘1’, for which the long-wave ZA phonon slope at vanishes; with an elastic (material) instability, as in ‘2’, for which the long-wave LA slope at vanishes; or with a short-wave 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 lattice-instabilities, as we will elaborate further in Sec.(3.7).

Figure 3: Phonon dispersions of graphene at certain critical states of deformation, marked as ‘1’, ‘2’ and ‘3’ in Fig. (3.5b), along the armchair deformation path (). Also shown are the deformed Brillouin zone (green polygon), and the irreducible wedge (red polygon). The sampled direction in each case is shown by an arrow. The labels ‘ZA’, ‘LA’ and ‘TA’ show the branch that goes unstable in the three cases.

### 3.6 From discrete dataset to continuum lattice-stability limit function via interpolation

Construction of a continuum-level 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 symmetry-invariants 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 ¯E)

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. (44-a)), 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:

 γc(¯E,Θ=Θn) ≡Rn(¯E/¯Ec) (44) =β1n(¯E/¯Ec)+β2n(¯E/¯Ec)2+β3n(¯E/¯Ec)3+...+βMn(¯E/¯Ec)M,

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. (44-b).
Figure 4: (a) Sampled valued of along two deformation paths: and (straight-line segments connect data points). (b) -order polynomial interpolation functions for along the armchair and zigzag directions — obtained by fitting functions described by equation 44.

#### 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 point-group symmetries of the lattice, and possesses periodicity with respect to :

 Kn(Θ)=α0n+α1ncos(Θ)+α2ncos(2Θ)+...; (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:

 Km(Θn)={0;Θn≠Θm1;Θn=Θm. (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 shape-functions are given as:

 K1(Θ)=12(1+cos(Θ)), (47) K2(Θ)=12(1−cos(Θ)). (48)

Following the above interpolation scheme, is generically expressed as:

 γc(¯E,Θ(θ))=n=N∑n=1Rn(¯E/¯Ec)Kn(Θ(θ)), (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:

 γc(¯E,Θ(θ))=n=N∑n=1M∑m=1γmn(¯E/¯Ec)mcos(n−1)Θ(θ); (50)

The various coefficients appearing in the expression of the failure function have been tabulated in the Tab. (2).

The quadrant of the instability surface resulting from the interpolation scheme is shown in Fig. (2-b).
Figure 5: (a) The radial interpolation of the lattice-stability limits along the two sampled directions, and . (b) A quadrant cut out from the limit surface. From below, the limit surface is bounded by the zero surface shown in green; while from above the failure surface is bounded by the critical surface , shown in red. The directions and correspond to the zigzag and armchair directions, respectively.

### 3.7 Limit surface as a smoothed envelope of various possible instabilities

Different kinds of strength-limiting 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 two-D structure, has an exceedingly small bending rigidity. Therefore, when subjected to an in-plane 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 non-negative, i.e., . This material instability can be either of the two types: long-wave (also called elastic) or a short-wave (also called soft mode) or

• a structural instability via out-of-plane 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 equi-biaxial 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 equi-biaxial 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 equi-biaxial 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.5-1). 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 short-wave instability, as shown in Fig. (3.5-3). Thus, .
Figure 6: Schematic illustration showing the lattice-stability plot as a smoothed representation of the envelope of the various possible lattice instabilities.

## 4 Lattice-stability limit surface in stress space

Often, it is useful to describe the lattice-stability 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 stress-strain response in graphene over the entire range of deformation up to the elastic stability limit. In a previous work (28), we derived a large-deformation 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 invariant-based representation theory. Our formulation employs the symmetry-invariants 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):

 ^ψDil(ϵa)=ψo[1−(1+αϵa)exp(−αϵa)]. (51)

The term  — the energetic response under an isochoric deformation — is given by a simple additive form:

 ^ψDev(ϵa,γ2i,γ3θ)=12μ(ϵa)γ2i+18η(ϵa)γ3θ, (52)

where is a second order elastic constant, and is a third-order elastic constant. The coefficients are determined from fits to ab initio energies for a set of deformed configurations (see Tab.[3]).

Differentiation of with respect to gives the work-conjugate stress tensor, i.e., .

 T(0)(ϵa,γ2i,γ3θ)=∂ψ∂E(0)=⎧⎪⎩∂^ψDil(ϵa)∂ϵa+∂^ψDev(ϵa,γ2i,γ3θ)∂ϵa⎫⎪⎭I+4∂^ψDev(ϵa,γ2i,γ3θ)∂(γ2i)E(0)0+∂^ψDev(ϵa,γ2i,γ3θ)∂(γ3θ)SE(0)0=⎧⎪⎩ψoα2ϵaexp(−αϵa)+12μ′(ϵa)γ2i+18η′(ϵa)γ3θ⎫⎪⎭I+2μ(ϵa)E(0)0+18η(ϵa)SE(0)0, (53)

and the Cauchy stress is given as (see Ogden (44))

 σ=1JFT(2)FT=21JF[∂E(0)∂C   ]T(0)FT, (54)

where denotes the right Cauchy-Green stretch tensor. Employing equation (53) and equation (54), we numerically map the contours in the strain plane to stress-plane. Once the Cauchy stress has been determined, we can obtain the principal stresses from the spectral representation of :

 σ=Smaxs1⊗s1+Smins2⊗s2, (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 lattice-stability 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 symmetry-reduced 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.,

 T(σ)≡τcmax(¯σ,Φ(ϕ))−τmax=0, τmax≥0; (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 lattice-stability limit surface in the space of symmetry-invariants of the Cauchy stress tensor, as shown in Fig. (56-b).
Figure 7: Numerical mapping of the instability surface from strain space (shown in (a)) to stress space (shown in (b)) via the hyperelastic constitutive relation of equation(53). The axes in stress space are in units of N/m. The directions and correspond to the zigzag and armchair directions, respectively.

## 5 Validation

The limit surface — given by equation (50) — contains information of all possible lattice-instabilities 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.5-a)).

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 :

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

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