## Abstract

Elastic material instabilities are precursors to failure in defect-free graphene single crystals. Elastic instabilities originate from softening in the material response (decay of tangent moduli) induced by dilatant mechanical deformation. Here, we characterize the softening in the constitutive response of graphene within the framework of hyperelasticity based on symmetry-invariants of the two-dimensional logarithmic strain tensor **E**^{(0)}. The use of symmetry-invariants provides significant functional simplification in representation of the strain energy function of graphene; *ab initio* calculations of deformation energy are well-fit using half the number of elastic constants of previous formulations. For a set of large homogeneous deformations comprising uniaxial stretch/stress along the armchair and the zigzag directions, and equi-biaxial tension, stress values predicted by the model compare well with those directly calculated from *ab initio* solutions. Using acoustic tensor analysis, we show that the constitutive model accurately captures elastic stability limits for a number of biaxial deformation modes, providing results that compare well with independent phonon calculations carried out using linear response density functional perturbation theory. In the case of equi-biaxial deformation, an elastic shearing instability is identified which occurs prior to the configuration of maximum true stress. Potential implications of the present results on the interpretation of limiting deformations achieved in nanoindentation experiments are briefly noted.

## 1. Introduction

Graphene is an atomic monolayer comprising a hexagonal network of *sp*^{2}-bonded C-atoms. Nanoindentation experiments on suspended graphene sheets suggest that defect-free, single-crystalline graphene can sustain near-ideal-strength stresses while remaining within the reversible regime of deformation [1–3]. These studies also indicate that graphene is intrinsically brittle and exhibits hyperelastic softening. The experimentation involves instrumented indentation of a suspended graphene sheet by a nanometre-sized diamond indenter. The nanoscale size of the indenter ensures that only a submicrometre region of the graphene sheet is probed to the largest deformation levels, and this region often comprises a single defect-free grain. Unfortunately, the stress and strain beneath the indenter are not directly measurable in such nanoindentation experiments; instead, the indentation load (*F*) as a function of indentation depth (*d*) during the course of indentation is recorded. Continuum-level finite-element analysis (FEA) is then used to estimate the stress and strain beneath the indenter, especially values corresponding to the measured fracture load and indentation depth [1,3,4]. For this purpose, a constitutive function for graphene appropriately describing the stress–strain response over the entire range of deformation up to failure is desired.

One approach to crystal hyperelasticity is based on Cauchy–Born lattice kinematics and tailored interatomic potentials. In applications near dislocations or crack tips, the use of interatomic potentials is generally required, as first principles calculations prove prohibitively costly. But continuum crystal hyperelasticity can be directly calibrated on the energetics of deformed but fully relaxed first principles models of small unit cells, bypassing approximations inherent in the definition of interatomic potentials.

Previously, numerous constitutive models of graphene have been proposed. Within the framework of hyperelasticity, constitutive response is characterized by the dependence of the strain energy density, or free energy per unit reference area, *ψ*, on the measure of imposed elastic strain, with all symmetries belonging to the material symmetry group duly incorporated [5]. Cadelano *et al.* [6] developed a nonlinear hyperelastic constitutive model of monolayer graphene by expressing *ψ* as a cubic polynomial in components of the Green–Lagrange strain measure, and then imposing the symmetry to reduce the number of independent fitting coefficients [7]. Wei *et al.* [1,8] noted that the model of Cadelano *et al.* failed to reproduce with sufficient fidelity the response at both infinitesimal and finite strains. They proposed a fifth-order series expansion for *ψ*, again based on components of the Green–Lagrange strain measure, to model the in-plane elastic properties of graphene. Their symmetry-reduced expression contained fourteen independent elastic constants. Xu *et al.* [9] modelled the constitutive response in terms of polynomials based on isotopic invariants of Green–Lagrange strain and accounted for the anisotropy evidenced at larger deformations by means of cumbersome switch-functions. Certain of these constitutive models offer better numerical accuracy than the others; however, owing to the polynomial nature of the representations, the underlying material physics has been—to some degree—obscured in all of these models. This obscurity is a by-product of the common modelling paradigm based on *reduction-by-symmetry*. The reduction-by-symmetry approach incorporates material symmetry into the constitutive modelling by extrinsically imposing symmetry constraints to reduce the number of independent fitting coefficients [7]. A major limitation, in this context, arises from the fact that the reduction-by-symmetry method is *a priori* restricted to the use of polynomials in the components of strain as the basis of representation—with no scope for including arbitrary non-polynomial functions in the representation.

One objective of this work is to present, and to illustrate with the example of graphene, the idea that an adequate elucidation of the underlying physics requires an appropriate choice of representation of the material response. To this end, we present a constitutive modelling scheme for crystal hyperelasticity based on symmetry-invariants of the logarithmic measure of strain **E**^{(0)}. The *symmetry-invariants* are scalar-valued functions of the tensor agency (here, strain) that satisfy all the symmetry constraints of the crystal’s material symmetry group [10,11]. Compared to reduction-by-symmetry, the symmetry-invariants-based approach offers substantial advantages. First, the resulting description of material response automatically satisfies the material symmetry requirements without further constraints. Secondly, in terms of the symmetry-invariants, *ψ* is expressed in a fashion independent of any particular choice of coordinate frame. Lastly, and most importantly, the approach does not restrict the representation to be only in terms of polynomials but allows for usage of more general functions. Using symmetry-invariants of **E**^{(0)} for the *C*_{6v} symmetry of graphene, we obtain a constitutive model for deformations involving both large area-change and shape-change. The proposed representation, involving both polynomial and more general functions, not only results in a functional simplification with fewer coefficients in the representation, but also introduces a certain transparency that affords clearer insights into the deformation physics of the material.

The other major objective of this work is to develop a predictive capability—using the obtained hyperelastic constitutive response—which allows the identification of elastic material instabilities in graphene. Elastic instabilities are the major mechanisms governing the elastic-to-inelastic transition in defect-free single crystals [12–14]. Such instabilities are precursors to incipient plasticity or the initiation of fracture in these crystals and originate from hyperelastic softening in the material response (decay of tangent moduli) induced by mechanical deformation. Provided this softening is well captured in the hyperelastic model, acoustic-tensor analysis enables identification of such instabilities. However, construction of the acoustic tensor can be a challenging task unless a convenient representation is employed. The present approach enables straightforward evaluation of the acoustic tensor in a frame-independent form. By evaluating elastic-stability limits for a number of deformation modes, results of which compare well with phonon-calculations carried out using density functional perturbation theory (DFPT), we demonstrate that the acoustic-tensor analysis constitutes a rigorous and effective means of predicting deformation-induced elastic instabilities in graphene.

In §2, we briefly discuss the kinematics of a two-dimensional deformable body and introduce relevant field variables and notations used in the formulation. In §3, we outline the general framework of hyperelasticity, Boehler’s principle of isotropy of space based on the structure tensor [15,16], and invariant-based representation theory [10,11]. The structure tensor and symmetry-invariants of **E**^{(0)} are explicitly derived for graphene in §4. Form-invariant tensors are obtained as derivatives of the symmetry-invariants with respect to **E**^{(0)}. In §4*a*, we propose a representation for *ψ* in the basis of the symmetry-invariants. The expression for work-conjugate stress in terms of the form-invariants and the conversion from the work-conjugate stress measure to Cauchy stress is presented in §4*b*. We obtain expressions for the work-conjugate tangent modulus tensor in §4*c*, and for the acoustic tensor in §4*d*. The model is validated in §5 by computing the stress for a set of finite deformations comprising uniaxial stretch/stress along the armchair and the zigzag directions; and equi-biaxial tension. The predictive capability of the model is demonstrated by obtaining elastic stability limits—based on acoustic tensor analysis—for a number of biaxial deformation modes in §6. We conclude in §7 by summarizing our results and considering certain of its implications.

## 2. Kinematics

We consider graphene as a two-dimensional deformable body denoted by unstressed reference configuration **X** an arbitrary material point of **X** moves to another point in the two-dimensional space, characterized by its deformed coordinate **x** at current time *t*. The convection of material points under deformation is described by a smooth, injective (one-to-one) function *χ*(**X**,*t*) called the motion. The non-translational part of the motion can be equivalently defined by the positive-definite second-order deformation gradient tensor, **F**=∇*χ*(**X**,*t*). Notationally suppressing this functional dependence for convenience, the polar decomposition theorem provides the following factorizations of **F** [17,18].
**R**∈*SO*_{2} characterizes rigid-body rotation, whereas **U** (or **V**=**RUR**^{T}), termed the right (left) Cauchy–Green tensor, characterizes shape- and area-change. Physically, deformation in the neighbourhood of a material point in the body can be kinematically considered as stretching followed by a superimposed rigid-body rotation, or vice-versa. **U** is a symmetric tensor having two real positive eigenvalues, the principal stretches *λ*_{1} and *λ*_{2}. Using spectral decomposition, **U** can be expressed as
**r**_{1} and **r**_{2} are orthogonal principal unit vectors in the plane. The deformation of a material point can be kinematically factored as the product of a purely dilatational (or shape-preserving, but area-changing) deformation **U**^{a}, with det **U**^{a}=det **U**, and a purely isochoric (or shape-changing, but area-preserving) deformation *J*=det **U**=*λ*_{1}*λ*_{2}, **I** is the two-dimensional identity tensor.

The spectral representation of **E**^{(0)}.

## 3. Invariant-theoretic approach to constitutive modelling

In a hyperelastic modelling framework, the strain energy density *ψ* is formally expressed as a scalar-valued function of deformation gradient **F**. Material frame indifference requires that this function should remain invariant under superimposed rigid body motion; i.e.
**Q**∈SO_{2} denotes a rigid-body rotation. Such objectivity is automatically satisfied if *ψ* is functionally dependent on one of the Seth–Hill strain measures — **E**^{(m)}≡((**F**^{T}**F**)^{m/2}−**I**)/*m*, for real *m*. In particular, for *m*=0, the corresponding strain measure is the logarithmic strain

The use of any particular strain measure is, in principle, arbitrary. However, different strain measures can lead to differing levels of complexity or simplification in accurately describing material response. For example, Anand [21] explored the extension of the classical quadratic strain energy function of isotropic linear elasticity based on two small-strain Lamé constants, replacing invariants of the infinitesimal strain tensor with corresponding invariants of various of the Seth–Hill strain measures. The formulation based on invariants of **E**^{(0)} (Hencky’s strain energy function) most accurately captured the initial constitutive nonlinearities at moderate deformations, with results clearly superior to those obtained by using invariants of the Green–Lagrange strain measure **E**^{(2)}.

Here, for reasons further described in §4*a*, we adopt the logarithmic strain tensor **E**^{(0)} to construct a hyperelastic constitutive response of graphene, i.e.
*ψ* due to material symmetry are expressed as
*ψ*, that remains invariant under a material symmetry group

The isotropicization theorem—based on the notion of a materially embedded structure tensor **E**^{(0)} and **Q**. The functions **E**^{(0)} for the structure tensor characterizing the material symmetry group of graphene.

## 4. Hyperelastic constitutive response of graphene

Following the approach outlined in §3, we now systematically construct a hyperelastic strain energy function for arbitrary in-plane deformation of graphene. First, we explicitly obtain the structure tensor characterizing the material symmetry group of graphene. The structure tensor ^{⊗n}=(…)⊗(…)⊗…⊗(…)(*n* times), and 2*n* denotes the order of the principal rotation axis. **M** and **N** are dimensionless symmetric traceless tensors given by

From equation (4.1), we obtain the sixth-order structure tensor characterizing

The complete and irreducible set of polynomial joint invariants of **E**^{(0)} and **E**^{(0)}. Following the procedure of Zheng & Betten [24], we obtain three independent scalar joint invariants of this two-dimensional system as
**A**:**B**=tr(**A**^{T}**B**) is the scalar tensor product, and
*ϵ*_{a} and *γ*_{i}≥0, are simply two isotropic invariants of **E**^{(0)} alone. Thus, any material anisotropy in the constitutive response of graphene, evidenced only at large deformations, is captured solely by the third invariant *γ*_{θ}.

In order to represent the second-order tensorial quantities, we need the following list of form-invariants:
**E**^{(0)} are obtained as follows:

### (a) Strain energy per unit reference area *ψ*

The proposed hyperelastic model is based on representation of the strain energy per unit reference area *ψ* in terms of the symmetry invariants of **E**^{(0)}; i.e. **E**^{(0)} offers substantial simplifications in terms of formulation. First, as is evident from its spectral representation (equation (2.6)), **E**^{(0)} additively decomposes areal (**U**^{a}) and isochoric (*ψ* in a hyperelastic material enables calculating *ψ* by integrating *dψ*=(∂*ψ*/∂**E**^{(0}):*d***E**^{(0)} along any convenient strain path, where it is understood that *T*^{(0)}≡∂*ψ*/∂**E**^{(0)} is the work-conjugate stress measure. Let a first, purely isotropic strain path (Path 1) correspond to areal deformation **U**^{a} while holding **U**^{a}=*J*^{1/2}**I**. The isotropic/deviatoric decomposition of the work-conjugate stress is **E**^{(0)}—can be additively decomposed into a term *ψ*^{Dil}(*ϵ*_{a})—corresponding to a pure areal deformation **U**^{a} while *ψ*^{Dev}(*γ*_{i},*γ*_{θ};*ϵ*_{a})—corresponding to a superimposed isochoric deformation **E**^{(0)} remains fixed at *ϵ*_{a}. Thus, we write
*ψ*=0 in the undeformed configuration, the areal contribution *ψ*^{Dil} equals the isotropic stress working along isotropic strain Path 1, along which *ψ*^{Dev} is numerically equal to the stress working along the subsequent deviatoric strain Path 2, along which *dϵ*^{a}=0; thus *ψ*^{Dev} depends on invariants *γ*_{i} and *γ*_{θ} of the imposed deviatoric strain *ϵ*_{a} along Path 2, with ‘initial’ condition *ψ*^{Dev}(0,0;*ϵ*^{a})=0.

From elasticity theory, we recall that the flexural stiffness *D* of a thin structural element scales with its thickness *h* as *D*∼*h*^{3}, whereas its stretching stiffness *C* scales linearly as *C*∼*h*. Notably, graphene is just one atomic layer thick, so in a superthin structure like graphene, the ratio *D*/*C* is exceedingly small and, accordingly, we assume that the contribution of bending to the strain energy is negligible compared to that of in-plane strain. Furthermore, and for the same reason, a suspended graphene sheet under a compressive in-plane loading, i.e. a state of Cauchy stress ** σ** with

**n**.

**.**

*σ***n**<0 for some in-plane direction

**n**, is structurally unstable in the limiting case of vanishing flexural stiffness, and will immediately buckle out-of-plane. Therefore, the scope of the modelling effort is limited to only those in-plane deformation states for which

**n**.

**.**

*σ***n**≥0 ∀

**n**.

#### (i) Energetic response under pure dilation, *ψ*^{Dil}

The energetic response under ‘Path 1’ deformations characterized by **U**=**U**^{a}=*J*^{1/2}**I** is well described by a function based on the universal binding energy relation (UBER) proposed by Rose *et al.* [25,26]. The UBER relation is
*α* as determined by fitting the *ab initio* energies, calculated using both LDA and GGA formulations, for the graphene lattice subjected to pure dilatory deformations. As shown in figure 2*a*, the UBER-based model for the fitted values of *α* accurately describes the *ab initio* energies for pure areal change.

#### (ii) Energetic response under shape-changing deformations, *ψ*^{Dev}

For reasons outlined above, the energies associated with various ‘Path 2’ deviatoric deformations are expected to depend on the value of *ϵ*^{a} at the end of ‘Path 1’; accordingly, we assume that

A simple linear combination of monomials in *ϵ*_{a}, i.e.
*ab initio* calculations well using simple functional forms for *μ*(*ϵ*_{a}) and *η*(*ϵ*_{a}). The shear modulus *μ* is well-fit by an exponentially decreasing function of the areal strain, *ϵ*_{a} (figure 2*b*),
*η*(*ϵ*_{a}) is fit by an even quadratic function of *ϵ*_{a} (figure 2*c*),
*ψ* is
*μ* and *η* are determined by fitting the set of *ab initio* energies calculated for a number of deformed states described in the electronic supplementary material; fitted values are summarized in table 2.

Using a total of only seven scalar fitting parameters, half the number used in the formulation of Wei *et al.*, the proposed functional form fits the entire DFT dataset very well.

### (b) Work-conjugate stress tensor T^{(0)}

Using the proposed functional forms for *ψ*, the stress measure *T*^{(0)} work-conjugate to **E**^{(0)} is calculated as
^{′} denotes differentiation with respect to *ϵ*_{a}. Following the hyperelastic work-conjugacy relation [27,28], the conversion from work-conjugate stress *T*^{(0)} to Cauchy stress ** σ** is obtained as follows. Letting a superposed dot denote the material time derivative, the power balance of isothermal hyperelasticity identifies measures of stress that are power-conjugate to differing measures of strain-rate by

*T*^{(2)}is its power-conjugate stress tensor, often denoted as the second Piola–Kirchhoff stress tensor. Using the chain rule, and viewing

**E**

^{(0)}as a function of

**E**

^{(2)}, we obtain

**can be similarly obtained from power balance relations [28] as**

*σ*As a particular illustration, we can use equation (4.19) to evaluate **T**^{(0)} for equi-biaxial stretch. In the absence of shape change, *γ*_{i}=*γ*_{θ}=0, so that
** σ** under equi-biaxial stretch becomes

*κ*is

### (c) Work-conjugate tangent moduli tensor L ( 0 )

The fourth-order tensor of tangent moduli connecting

### (d) Acoustic tensor A

The acoustic tensor is defined as the second derivative of the *ψ* with respect to the deformation gradient tensor **F** [7,28,29]:
*ψ* with respect to **F** gives the generally non-symmetric first Piola–Kirchhoff stress tensor, *T*^{PK1}, which can be evaluated using the chain rule as:

In the present applications, the acoustic tensor **e**_{i}⊗**e**_{j}⊗**e**_{k}⊗**e**_{l}}_{1≤i,j,k,l≤2} is given by
*i*,*j*,*k*,*l*≤2. By definition, *f*(1,1)=1, *f*(2,2)=2, *f*(1,2)=3, and *f*(2,1)=4. In this alternate basis,

## 5. Validation of the constitutive model

As detailed below, we validate our continuum model for a number of homogeneous deformations by comparison with small-strain elastic constants inferred from experiments, comparison of stress–strain curves with independent *ab initio* calculations, and comparisons of predictions of elastic stability limits to independent phonon calculations.

### (a) Comparison of small-strain elastic constants and wave-velocities with measured values

The in-plane elastic constants—recovered from the constitutive model in the limit of infinitesimal strain—are compared with the measured values of Lee *et al.* [1]. The predicted values for the in-plane Young’s modulus *Y* ^{(0)}, shear modulus *μ*^{(0)}, bulk modulus *κ*^{(0)} and Poisson’s ratio *ν*^{(0)} are all in good accord with the experimentally reported values (table 3). The acoustic wave-velocities *c*_{LA} and *c*_{TA} obtained from the continuum model also compare very well with measured values, as well as with values obtained from phonon dispersion, as shown in table 4.

### (b) Comparison with directly calculated *ab initio* stress response curves

The constants of the continuum model are determined by a least-squares fit to the strain energies alone (stress values were not used in the fitting), so the agreement of predicted stress values with those determined directly from *ab initio* calculations remains to be examined. Next we compare the continuum model’s predictions of stress response with corresponding *ab initio*-calculated values in a few important modes of large homogeneous deformation, including equi-biaxial tension, and uniaxial stretching and stressing in the armchair and the zigzag directions. In obtaining the stress response, the elastic and soft-mode instabilities—discussed in previous sections—have been suppressed in both *ab initio* and continuum calculations.

#### (i) Pure biaxial stretch

Cauchy stress as a function of strain for pure equi-biaxial deformation as obtained from the continuum model (equation (4.28)) is shown in figure 3*a*. The Cauchy biaxial stress as a function of areal strain *ϵ*_{a}, as predicted by the UBER-based model, compares well with the *ab initio* values. Our model predicts that the Cauchy stress reaches its maximum value at nearly 42% areal strain. In figure 3*b*, we show the softening of the dilatant tangent modulus, defined in equation (4.29), with areal strain *ϵ*_{a}. The areal tangent modulus vanishes when the Cauchy stress approaches its maximum value.

#### (ii) Uniaxial stretching along the armchair and the zigzag directions

Cauchy stress components as a function of stretch, obtained from the continuum model for the case of uniaxial stretching along the armchair and the zigzag directions, are shown in figure 4*a*,*b*, respectively. Orientation-dependent differences in response at large deformation evidence a limited anisotropy that is fully captured by the strain energy contribution ‘

#### (iii) Uniaxial tension along the armchair and the zigzag directions

Cauchy stress components as functions of stretch obtained from the continuum model for the case of uniaxial tensile stress along the armchair and the zigzag directions are shown in figure 5*a*,*b*, respectively. The ratios of transverse strain (due to Poisson contraction) to longitudinal strain obtained from the continuum model for uniaxial stress along the armchair and zigzag directions are also in good agreement with the *ab initio* values (see the electronic supplementary material, Fig. [S2]).

## 6. Predictive capability: elastic stability limits for some biaxial deformation modes

Material elastic stability requires that the speed of propagation of all acceleration waves in a solid should be non-negative [32]. This condition is satisfied when the acoustic tensor is positive-definite, which is equivalent to the condition that the continuum elasto-dynamic equations are hyperbolic.

The acoustic tensor-based prediction of elastic instability involves detecting a deformed state at which

Employing the continuum expression for *ab initio* linear response phonon calculations.

### (a) Pure biaxial stretch

For this deformation state, Cartesian components of **m** and **n**. If, for any two unit vectors, equation (6.2) does not hold, then the deformed crystal is said to have lost strong ellipticity. Substituting for **m** and **n** be represented by
*ϵ*_{a}<1/(1+*α*)=*ϵ*_{a}|_{κ=0}, the minimum value of the term containing the factor ** m**.

**n**=0; therefore, in this range of deformation

*ϵ*

_{a}given by

*ϵ*

_{a}|

_{μ=0}<

*ϵ*

_{a}|

_{κ=0}, as assumed. Since the associated directions

**m**and

**n**are perpendicular, the predicted dynamic instability is of transverse acoustic nature. The model’s capture of this instability is confirmed in independent linear-response-based phonon calculations for equi-biaxially strained graphene. The LDA phonon dispersion shows a vanishing long-wavelength transverse acoustic wave speed emerges at

*μ*(figure 6). Based on our constitutive modelling and on independent phonon calculations, the loss of elastic stability under equi-biaxial loading occurs prior to the stress peak associated with the zero tangent modulus condition.

### (b) Uniaxial strain in armchair and zigzag directions

Uniaxial stretching along both armchair and zigzag directions preserves reflection symmetries w.r.t the *λ*^{(c)} at which *Λ*_{1}, *Λ*_{2}, *Λ*_{3} and *Λ*_{4} with increasing uniaxial stretch.

First we consider uniaxial stretching along the zigzag direction, oriented along the Cartesian **e**_{1}-axis, with **U**=*λ*_{s}**e**_{1}⊗**e**_{1}+1**e**_{2}⊗**e**_{2}. For this case, the loss of strong ellipticity occurs at a critical stretch, *Λ*_{2}-eigenvalue of *a*. The associated eigenvector corresponds to ** m**=

**n**=

**e**

_{1}(see equation (6.2)), so the unstable longitudinal mode coincides with the maximum principal eigenvector of

**U**. Independent phonon dispersion calculations shown in figure 7

*b*, confirm the occurrence of a longitudinal acoustic instability in this direction, at

*λ*

^{(c)}

_{s}≈1.188. Now consider uniaxial stretching along the armchair direction, which is taken to parallel Cartesian axis

**e**

_{2}, with

**U**=1

**e**

_{1}⊗

**e**

_{1}+

*λ*

_{a}

**e**

_{2}⊗

**e**

_{2}. In this case,

*λ*

^{(c)}

_{a}≈1.23 (GGA) and ≈1.24 (LDA), when

**=**

*m***n**=

**e**

_{2}, again being a longitudinal elastic instability in the direction of maximum principal stretch. This is confirmed by independent LDA phonon calculations shown in the electronic supplementary material, Fig. (S3-b), which exhibits a vanishing slope of the LA branch at

*Γ*when

*λ*

^{(c)}

_{a}=1.252.

### (c) Uniaxial stress in armchair and zigzag directions

In cases of uniaxial stress along the armchair and zigzag directions, **U**=*λ*_{s}**e**_{1}⊗**e**_{1}+*f*(*λ*_{s})**e**_{2}⊗**e**_{2}, where the transverse stretch *f*(*λ*_{s})≤1, and its value, for a given *λ*_{s}, is determined by setting *σ*_{22}=**e**_{2}.** σ**.

**e**

_{2}=0. For uniaxial stress along the zigzag direction, the loss of strong ellipticity occurs at

*λ*

^{(c)}

_{s}=1.19, when

**=**

*m***n**=

**e**

_{1}, an unstable longitudinal mode parallels to the principal stretch direction. Independent phonon calculations confirm the occurrence of an acoustic instability in the longitudinal branch of the phonon dispersion at

*λ*

^{(c)}

_{s}≈1.20, as was also indicated by phonon calculations of Liu

*et al.*[14].

For uniaxial tensile stress along the armchair direction, **U**=*g*(*λ*_{a})**e**_{1}⊗**e**_{1}+*λ*_{a}**e**_{2}⊗**e**_{2}, where the transverse stretch *g*(*λ*_{a})≤1, its value, for a given *λ*_{a}, being determined by setting *σ*_{11}=0. In this case, acoustic tensor analysis gives an elastic instability at *λ*^{(c)}_{a}=1.24, as shown in the electronic supplementary material, Fig. (S4-a). This instability is also longitudinal, in the **e**_{2}-direction, and is also in good agreement with phonon calculations shown in the electronic supplementary material, Fig. (S4-b). The predicted elastic stability limits for various deformation modes considered in this work are summarized in table 5. Respective stretch limits for elastic stability under both uniaxial stress and stretch, and in both zigzag and armchair directions, are quite close; this is presumably related to graphene’s relatively small Poisson ratio, and the diminishing ratio of transverse to axial strain magnitude as the lattice is deformed under uniaxial stress (see the electronic supplementary material, Fig. (S2)).

## 7. Discussion and conclusion

Using as a basis scalar-valued symmetry invariants of the logarithmic strain tensor, we derived a nonlinear hyperelastic constitutive response for graphene. Because the model employs symmetry invariants, the material symmetry group of the underlying lattice is built-into the model, and the need for externally imposing the symmetry restrictions is eliminated. The derivation is coordinate-free, facilitating computational implementation. Derivation of higher order tensor variables such as work-conjugate tangent moduli tensors and the acoustic tensor is likewise coordinate-free. The model elucidates the respective contributions to the strain energy density due to purely equi-biaxial area change, and due to purely isochoric shape-changing deformations. The model correctly follows the stress–strain response of graphene in both uniaxial stretching and tension, along both the zigzag and armchair directions, and in equi-biaxial tension. Values of the isotropic small-strain elastic constants deduced from the model are also in good agreement with experimentally reported values. Acoustic-tensor-based stability analysis of the constitutive model provides predictions of the deformation at the limits of strong ellipticity, and the corresponding critical modes, that are in good agreement with the results of independent lattice dynamics calculations based on linear response perturbation theory.

Our model predicts that the initial loss of elastic stability in graphene under pure biaxial stretch occurs at a logarithmic areal strain of *et al.* [8], but it could occur only at a larger deformation than that of the presently identified elastic shearing instability. Moreover, as noted previously by Yevick & Marianetti [34], phonon calculations for equi-biaxial deformation of graphene also show that, prior to the onset of the long-wavelength shear instability, a short-wavelength instability (or ‘soft mode’) at *K* emerges at *ϵ*_{a}=0.28−0.30, as seen also in figure 6*c*.

The experimental results of Lee *et al.* [1] have been interpreted, based in part on continuum finite-element modelling [4], as suggesting that local equi-biaxial Cauchy stress levels in the suspended graphene beneath the tip of a frictionless diamond nanoindenter reach the ideal peak value occurring at

- Received July 25, 2014.
- Accepted October 29, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.