Structural deformations of two-dimensional planar structures under uniaxial strain: The case of graphene

TL;DR

This paper presents a molecular mechanics method to predict structural deformations of two-dimensional materials under uniaxial strain, validated on graphene with improved accuracy using second nearest neighbor interactions.

Contribution

The study introduces a modified stick and spiral model including second nearest neighbor interactions for better deformation predictions in 2D structures.

Findings

Original model insufficient for graphene deformation prediction.

Second nearest neighbor interactions significantly improve accuracy.

Method can be extended to other strain conditions.

Abstract

In the present work, a method for the study of the structural deformations of two dimensional planar structures under uniaxial strain is presented. The method is based on molecular mechanics using the original stick and spiral model and a modified one which includes second nearest neighbor interactions for bond stretching. As we show, the method allows an accurate prediction of the structural deformations of any two dimensional planar structure as a function of strain, along any strain direction in the elastic regime, if structural deformations are known along specific strain directions, which are used to calculate the stick and spiral model parameters. Our method can be generalized including other strain conditions and not only uniaxial strain. We apply this method to graphene and we test its validity, using results obtained from {\it ab initio} Density Functional Theory calculations. What we find is that the original stick and spiral model is not appropriate to describe accurately the structural deformations of graphene in the elastic regime. However, the introduction of second nearest neighbor interactions provides a very accurate description.

Tables

Table 1. Table 1: Values of λ i subscript 𝜆 𝑖 \lambda_{i} , μ i ​ j subscript 𝜇 𝑖 𝑗 \mu_{ij} and A 𝐴 A obtained from the fittings for the four strain directions.

$n$	$m$	${\theta }_0$ (o)	$i$	${\theta }_{0i}$ (o)	${\cos }^2{\theta }_{0i}$	${\lambda }_i$	$\cos 2{\theta }_{0i^{\prime }}$	${\mu }_{j^{\prime }{k}^{\prime }}$
1	-1	90.000000	3	270.000000	0.000000	-0.001556	-1.000000	1.315279
4	-5	100.893395	3	259.106605	0.035714	0.008633	-0.928571	1.221761
2	1	10.893395	1	109.106605	0.107143	0.027796	-0.785714	1.032964
1	1	0.000000	1, 2	120.000000	0.250000	0.066506	-0.500000	0.654704
2	1	10.893395	2	229.106605	0.428571	0.116258	-0.142857	0.185116
4	-5	100.893395	2	139.106605	0.571429	0.156141	0.142857	-0.190542
1	-1	90.000000	1, 2	30.000000	0.750000	0.206426	0.500000	-0.657640
4	-5	100.893395	1	19.106605	0.892857	0.246905	0.785714	-1.031171
2	1	10.893395	3	349.106605	0.964286	0.267113	0.928571	-1.218509
1	1	0.000000	3	360.000000	1.000000	0.277621	1.000000	-1.309408

Equations

$$U=U_s+U_b+U_{\omega }+U_{\tau }+U_{vdw}+U_e,$$

$$U=U_s+U_b.$$

$$U=\frac{1}{2}\sum _i\left(k_{s,i}\delta l_i^2+\frac{1}{2}\sum _j{k}_{b,ij}\delta {\varphi }_{ij}^2\right),$$

$${\varphi }_{i{j}_1}+{\varphi }_{i{j}_2}+{\varphi }_{j_1{j}_2}={\varphi }_{i{j}_3}+{\varphi }_{i{j}_4}+{\varphi }_{j_3{j}_4}=2\pi ,$$

$$\delta {\varphi }_{i{j}_1}+\delta {\varphi }_{i{j}_2}+\delta {\varphi }_{j_1{j}_2}=\delta {\varphi }_{i{j}_3}+\delta {\varphi }_{i{j}_4}+\delta {\varphi }_{j_3{j}_4}=0.$$

$$U=\frac{1}{2}\left(\sum _i{n}_i{k}_{s,i}\delta l_i^2+\frac{1}{2}\sum _i\sum _j{m}_{ij}{k}_{b,ij}\delta {\varphi }_{ij}^2\right),$$

$$\Lambda =U+\lambda (\epsilon -\delta L/L_0),$$

$$\delta L$$

$$\delta x$$

$$\delta y$$

$$\delta x$$

$$\delta L$$

$$\delta L=\sum _i{q}_i\left(\delta l_i\cos {\theta }_{0i}-l_{0i}\sin {\theta }_{0i}\delta {\theta }_i\right)$$

$$\epsilon -\sum _i{q}_i\left(\delta l_i\cos {\theta }_{0i}-l_{0i}\sin {\theta }_{0i}\delta {\theta }_i\right)/L_0=0.$$

$$\forall {\varphi }_{ij}\in (0,\pi ],\delta {\varphi }_{ij}^2=(\delta {\theta }_j-\delta {\theta }_i{)}^2,$$

$$\Lambda$$

$$\sum _i{q}_i(\delta l_i\sin {\theta }_{0i}+l_{0i}\cos {\theta }_{0i}\delta {\theta }_i)=0.$$

$$\partial \Lambda /\partial \delta l_i=0\text{and}\partial \Lambda /\partial \delta {\theta }_i=0.$$

$$\delta l_i=\frac{\lambda }{L_0}\frac{q_i}{n_i}\frac{\cos {\theta }_{0i}}{k_{s,i}}.$$

$$\frac{1}{2}\sum _{k=1}^4{m}_{i{j}_k}{k}_{b,i{j}_k}(\delta {\theta }_i-\delta {\theta }_{j_k})=-\lambda q_i{l}_{0i}\sin {\theta }_{0i}/L_0.$$

$$U_{min}=\lambda \epsilon /2,$$

$$U=U_1+U_2=U_1+(1/2)\sum _i\sum _j{p}_{ij}(1/2)k_{s,ij}\delta r_{ij}^2,$$

$$n_i{k}_{s,i}\delta l_i+\frac{1}{2}\sum _{k=1}^4{m}_{i{j}_k}{k}_{s,i{j}_k}\delta r_{i{j}_k}\frac{\partial \delta r_{i{j}_k}}{\partial \delta l_i}=\lambda q_i\cos {\theta }_{0i}/L_0$$

$$\frac{1}{2}\sum _{k=1}^4{m}_{i{j}_k}\left[k_{b,i{j}_k}(\delta {\theta }_i-\delta {\theta }_{j_k})+k_{s,i{j}_k}\delta r_{i{j}_k}\frac{\partial \delta r_{i{j}_k}}{\partial \delta {\theta }_i}\right]$$

$$=-\lambda q_i{l}_{0i}\sin {\theta }_{0i}/L_0,$$

$$r_{0ij}\delta r_{ij}$$

$$\partial \delta r_{ij}/\partial \delta l_i=\left[l_{0i}\mp l_{0j}\cos ({\theta }_{0i}-{\theta }_{0j})\right]/r_{0ij},$$

$$\partial \delta r_{ij}/\partial \delta {\theta }_i=\pm l_{0i}{l}_{0j}\sin ({\theta }_{0i}-{\theta }_{0j})/r_{0ij}.$$

$$U=U_1$$

$$U_2=k_{s2}(\delta r_{1,2}^2+\delta r_{2,3}^2+\delta r_{3,1}^2),$$

$$\mathbf{a}={\mathbf{r}}_3-{\mathbf{r}}_2\text{and}\mathbf{b}={\mathbf{r}}_3-{\mathbf{r}}_1.$$

Full text

I Introduction

Undoubtedly, graphene is one of the most studied materials in recent years. This is due to its exotic properties, like for instance its high carrier mobility graphene_electron_mobility and high thermal conductivity DT217 ; ghosh ; susp_graphene_exper_TC at room temperature, its high strength PCCP_Fthenakis ; strength , etc, which makes graphene one of the most interesting materials for future nanoelectronic and nanomechanic applications. Following graphene, several two dimensional (2D) materials have also gained interest, exhibiting interesting mechanical mech_properties ; mech_BN ; mech_C2F_BN ; mech_MoS2 ; mech_silicene_1 ; mech_silicene_2 ; mech_phosphorene ; Fthenakis_Si2BN and electronic properties. The world of 2D materials that have been brought to the center of attention recently review_2D_1 ; review_2D_2 includes several transition metal dichalcogenides TMDs ; mech_MoS2 , (like for instance MoS2 or WS2), hexagonal BN (h-BN) mech_BN ; mech_C2F_BN ; hex_BN ; many , Si2BN Fthenakis_Si2BN ; Madhu_Si2BN , SinBm SiB ; Si_nB_m ; BSi3 , SiX and XSi3 (X=B, C, N, Al, P) SiX , CdS CdS_planar , AlN hex_AlN ; hex_AlN_theory ; many ; hex_AlN_mechanical , SiC, InN and GaN many , C2F mech_C2F_BN , Silicene mech_silicene_1 ; mech_silicene_2 ; silicene ; silicene_germanene , Germanene silicene_germanene , Siligene (SiGe) SiGe_2D , Phosphorene mech_phosphorene ; phosphorene , as well as several graphene allotropes, like pentaheptites and octagraphene PCCP_Fthenakis ; enyashin , or other Carbon 2D allotropes, like pentagraphene pentagraphene , graphyne, graphydine enyashin ; allotrop , or graphene-based derivatives, like graphane and graphone graphane ; allotrop etc.

A special class of these materials are those which are entirely planar, like for instance several graphene allotropes (pentaheptites, octagraphene, etc) PCCP_Fthenakis ; enyashin ; haekel , as well as h-BN hex_BN , Si2BN Fthenakis_Si2BN ; Madhu_Si2BN , AlN, SiC, SinBm SiB ; Si_nB_m ; BSi3 , CdS CdS_planar , XSi3 with X=B,C,Al SiX etc. In this work, we present a method for the study of the mechanical response, of these materials, e.g. bond stretching and angle bending deformations, in the presence of uniaxial tensile strain, providing analytic expressions for these deformations along any strain direction. Our method can be generalized including any other strain condition (i.e. not only uniaxial strain) and is based on molecular mechanics assuming two different versions of the so called stick and spiral model stick_spiral , which has been employed previously for the study of the mechanical properties of Carbon nanotubes Geng ; Chang ; Zhao1 ; Zhao2 ; Zhao3 ; Jiang .

As an example, we apply our method to graphene, providing analytic expressions for bond length and bond angle deformations under tensile strain. We test the accuracy of these expressions using results we obtain from ab-initio density functional theory (DFT) calculations. In particular, we calculate the structural deformations of graphene under tensile strain along the high symmetry arm chair and zig-zag directions, as well as two other randomly selected directions, which are perpendicular to each other. According to our findings, the original stick and spiral model is not sufficient to provide an accurate description of the mechanical deformations of graphene under tensile strain in the elastic regime, since the DFT results can not be reproduced accurately by the analytic expressions provided by that model. However, due to the coupling between the bond stretching and angle bending terms, which is inherently included in the modified stick and spiral model, this modified model provides a quite accurate description. Moreover, fitting these analytic expressions to the DFT results we calculate the force constants for bond stretching and angle bond bending for graphene, thus allowing the prediction of the mechanical response of graphene in the elastic regime for strain on any direction.

II The deformation energy

In molecular mechanics approach the deformation energy $U$ is a sum of energy contributions from different deformation modes stick_spiral . In particular, $U$ is written as

[TABLE]

where $U_{s}$, $U_{b}$, $U_{\omega}$, $U_{\tau}$, $U_{vdw}$ and $U_{e}$ correspond to the energy contributions from bond stretching, bond angle bending, bond inversion, bond angle torsion, Van der Walls interactions and electrostatic interactions, respectively. Since tensile strain in a 2D planar structure is in-plane strain, the terms $U_{\omega}$ and $U_{\tau}$ vanish. Moreover, since there are no interactions between different sheets of those 2D structures, the terms $U_{vdw}$ and $U_{e}$ also vanish. Thus, the deformation energy becomes

[TABLE]

$U_{s}$ and $U_{b}$ may be expressed in several different ways (see for instance Refs. Davydov2010, ; kalosakas, ; lobo, ). However, the simplest way is to be expressed as a sum of harmonic terms constituting the so-called stick and spiral model.

According to the stick and spiral model, the deformation energy per unit cell is written as a sum of energy contributions from each bond length and bond angle deformation. Each of these contributions has a quadratic dependence on the corresponding deformation, i.e. it is either of the form $(1/2)k_{s}\delta l^{2}$ (for bond stretching), or $(1/2)k_{b}\delta\phi_{ij}^{2}$ (for bond-angle bending), where $k_{s}$ and $k_{b}$ are the corresponding force constants, and $\delta l$ and $\delta\phi$ the bond length and bond-angle deformations for each specific bond and bond angle, respectively. Thus, the deformation energy per unit cell is

[TABLE]

where $i$ counts all the bonds inside the unit cell and $j$ counts the bonds which form bond angles with bond $i$. The $1/2$ factor of the second sum is to avoid double counting of the bonds.

In the description provided by the stick and spiral model, bond stretching and bond angle bending are not coupled. The energy provided by Eq. (3) does not have any terms mixing these deformations. In addition, as we will see later, in the minimization of the deformation energy under constant strain these deformations remain decoupled. More specifically, one arrives at two independent systems of analytic equations one for stretching and one for bending. A more accurate description would include a coupling term between these deformations. This can be achieved by introducing extra terms describing the stretching of second nearest neighbor interatomic distances. In the present work, we study both cases.

For a planar structure with three-fold coordinated atoms, there are three bonds and three bond angles per atom (see Fig. 1(a)). If we label $i$, $j_{1}$ and $j_{2}$ the bonds of atom A and $i$, $j_{3}$ and $j_{4}$ those of atom B, (the two atoms share the bond $i$), then the index $j$ of Eq. (3) takes the values $j_{1}$, $j_{2}$, $j_{3}$ and $j_{4}$. Moreover, since the structure is planar, and all atoms remain in the plane under tensile strain

[TABLE]

where $\phi_{ij_{1}}$, $\phi_{ij_{2}}$, $\phi_{j_{1}j_{2}}$ are the bond angles of atom A and $\phi_{ij_{3}}$, $\phi_{ij_{4}}$, $\phi_{j_{3}j_{4}}$ the bond angles of atom B. Consequently,

[TABLE]

In the present work we study structures with only 3-fold coordinated atoms, since this is the most common case. However, the generalization of our method to structures with $n$-fold coordinated atoms, with $n\neq 3$, is obvious.

Due to symmetry reasons (if any), several bonds length deformations (as well as bond angle deformations) may be equivalent with each other under specific strain conditions. In that case, $U$ can be written as a function of only the independent bond length and bond angle deformations per unit cell, and Eq. (3) can be rewritten as

[TABLE]

where $n_{i}$ is the number of equivalent bond length deformations of type $i$ and $m_{ij}$ the number of equivalent bond angle deformations formed by the bonds which have independent bond length deformations of type $i$ and $j$. $i$ runs over the independent bond deformations only.

Under uniaxial strain, the deformation energy and the corresponding deformations $\delta l_{i}$ and $\delta\phi_{ij}$ at the strained equilibrium can be found from the minimization of the deformation energy subject to constrains describing the strain condition. These constraints can be incorporated using the Lagrange multipliers technique. For constant uniaxial tensile strain $\varepsilon$ there is only one constraint described by $\varepsilon=\delta L/L_{0}$, where $L_{0}$ is a length along the strain direction and $\delta L$ the elongation of $L_{0}$ upon that strain, which should be expressed as a function of the independent variables $\delta l_{i}$ and $\delta\phi_{ij}$. Thus, the function which should be minimized becomes

[TABLE]

with $\lambda$ the corresponding Lagrange multiplier. Obviously, for different strain conditions, different constrains will apply, which can be incorporated in Eq. (7) using the corresponding Lagrange multipliers. Thus, our method can be easily generalized to describe the structural deformations of a 2D planar structure, not only under uniaxial strain, but under any strain condition.

In order to minimize $\Lambda$ in Eq. (7), with respect to the bond stretching and angle bending deformations, one needs to express $\delta L$ in terms of these deformations.

II.1 $\delta L$ as a function of bond deformations

Without loss of generality, we may assume that the structure is periodic. A non-periodic (i.e. amorphous) structure could be considered as periodic with infinite periodicity. For convenience, let us assume that the unit cell vectors for $\varepsilon=0$ are $\mathbf{a}_{0}=a_{0}\hat{\mathbf{i}}$ and $\mathbf{b}_{0}=b_{x0}\hat{\mathbf{i}}+b_{y0}\hat{\mathbf{j}}$, as shown in Fig. 2. Let us apply tensile strain by stretching the structure along the line connecting two equivalent atoms in different unit cells. The vector connecting those two atoms, (which determine the strain direction), is $\mathbf{L}_{0}=n\mathbf{a}_{0}+m\mathbf{b}_{0}$, where $n$ and $m$ are integers. Under the applied strain the vector $\mathbf{L}_{0}$ will be deformed to $\mathbf{L}$, so that the vectors $\mathbf{L}$ and $\mathbf{L}_{0}$ are parallel, i.e. $\mathbf{L}_{0}$ will be just elongated. The unit cell vectors $\mathbf{a}_{0}$ and $\mathbf{b}_{0}$ will be also deformed to $\mathbf{a}$ and $\mathbf{b}$, respectively, so that $\mathbf{L}_{0}=n\mathbf{a}_{0}+m\mathbf{b}_{0}\parallel\mathbf{L}=n\mathbf{a}+m\mathbf{b}$.

If $L_{0}$ and $L=L_{0}+\delta L$ are the lengths of the vectors $\mathbf{L}_{0}$ and $\mathbf{L}$, respectively, and $\hat{\boldsymbol{\varepsilon}}$ is the unit vector directed along the strain direction (i.e. $\hat{\boldsymbol{\varepsilon}}=(n\mathbf{a}_{0}+m\mathbf{b}_{0})/(n^{2}a_{0}^{2}+m^{2}b_{0}^{2}+2nm\mathbf{a}_{0}\mathbf{b}_{0})^{1/2}$, where $b_{0}=(b_{x0}^{2}+b_{y0}^{2})^{1/2}$), then $L=\hat{\boldsymbol{\varepsilon}}(n\mathbf{a}+m\mathbf{b})=n(\hat{\boldsymbol{\varepsilon}}\mathbf{a})+m(\hat{\boldsymbol{\varepsilon}}\mathbf{b})$ and $L_{0}=\hat{\boldsymbol{\varepsilon}}(n\mathbf{a_{0}}+m\mathbf{b_{0}})=n(\hat{\boldsymbol{\varepsilon}}\mathbf{a_{0}})+m(\hat{\boldsymbol{\varepsilon}}\mathbf{b_{0}})$, i.e. $L$ ($L_{0}$) depend on the projections of $\mathbf{a}$, and $\mathbf{b}$ ($\mathbf{a}_{0}$ and $\mathbf{b}_{0}$) on the strain direction.

The vectors $\mathbf{a}_{0}$ and $\mathbf{b}_{0}$ can be expressed as a sum of bond vectors $\mathbf{r}_{0ai}$ and $\mathbf{r}_{0bi}$, respectively, ($i=1,2,3,...$), which correspond to specific bonds of the undeformed structure, constituting a crooked line connecting the tails of $\mathbf{a}_{0}$ and $\mathbf{b}_{0}$ with their heads, i.e. $\mathbf{a}_{0}=\sum_{i}\mathbf{r}_{0ai}$ and $\mathbf{b}_{0}=\sum_{i}\mathbf{r}_{0bi}$. Thus, if the bond vectors $\mathbf{r}_{0ai}$ and $\mathbf{r}_{0bi}$ are deformed under strain into $\mathbf{r}_{ai}$ and $\mathbf{r}_{bi}$, respectively, then $\mathbf{a}=\sum_{i}\mathbf{r}_{ai}$ and $\mathbf{b}=\sum_{i}\mathbf{r}_{bi}$. This is shown schematically in Fig. 2, where the sum of the red colored vectors, (denoted as $\mathbf{r}_{ai}$, $i=1,2,3,...$), constitute $\mathbf{a}$, while the sum of the green colored vectors, (denoted as $\mathbf{r}_{bi}$, $i=1,2,3,...$), constitute $\mathbf{b}$. Obviously, the corresponding sums of the projections of $\mathbf{r}_{ai}$ and $\mathbf{r}_{bi}$ along the strain direction equals the projection of $\mathbf{a}$ and $\mathbf{b}$, respectively, along the same direction. These projections of $\mathbf{r}_{ai}$ and $\mathbf{r}_{bi}$ are shown as black arrows in Fig. 2, and should be considered as positive or negative. Thus,

[TABLE]

i.e. $\delta L$ can be expressed as a function of the differences of the projections of the $\mathbf{r}_{0ai}$, $\mathbf{r}_{ai}$ and the $\mathbf{r}_{0bi}$, $\mathbf{r}_{bi}$ vectors, along the strain direction. We should note that, although the vectors, $\mathbf{a}$, $\mathbf{b}$, $\mathbf{a}_{0}$, $\mathbf{b}_{0}$ are not uniquely expressed in terms of bond vectors, the sums of the projections are unique and one could always choose optimal paths (e.g. of minimal length) of bond vectors. Let us now see how the differences of those projections depend on the bond deformations.

II.2 The strain constrain

Let us assume that strain along a specific direction is applied to a bond, as shown in Fig. 1(b). For convenience we have assumed that the strain direction coincides with the x-axis direction. Let us further assume that at equilibrium for $\varepsilon=0$, the bond length and the angle between the bond and the strain direction are $l_{0}$ and $\theta_{0}$, and under strain they become $\theta_{0}+\delta\theta$ and $l_{0}+\delta l$, respectively. If the projections of the bond along and normal to the strain direction for $\varepsilon=0$ are $x_{0}$ and $y_{0}$, respectively, and under strain they are $x_{0}+\delta x$ and $y_{0}+\delta y$, respectively, then $x_{0}=l_{0}\cos\theta_{0}$, $y_{0}=l_{0}\sin\theta_{0}$, $x_{0}+\delta x=(l_{0}+\delta l)\cos(\theta_{0}+\delta\theta)$ and $y_{0}+\delta y=(l_{0}+\delta l)\sin(\theta_{0}+\delta\theta)$.

Thus the projection of the bond deformation along the strain direction is

[TABLE]

and the projection normal to the strain direction is

[TABLE]

According to Eq. (9), the projection $\delta x$ of the deformation of $\mathbf{r}_{0ai}$ along the strain direction $\hat{\boldsymbol{\varepsilon}}$ is

[TABLE]

where $l_{0ai}=|\mathbf{r}_{0ai}|$, $\theta_{0ai}$ is the angle between $\mathbf{r}_{a0i}$ and the strain direction (i.e. $\cos\theta_{0ai}=\hat{\boldsymbol{\varepsilon}}\mathbf{r}_{0ai}/l_{0ai}$), and $\delta l_{ai}$ and $\delta\theta_{ai}$ are the deformations of $l_{0ai}$ and $\theta_{0ai}$, respectively. Changing the index ”$a$” with ”$b$”, we get the corresponding relation for $\mathbf{r}_{0bi}$. Consequently,

[TABLE]

As a function of the projections of independently deformed bonds, this equation is written as

[TABLE]

where here index $i$ is the same as in Eq. (6), (i.e. it runs over the bond vectors of the independently deformed bonds) and $q_{i}$ is the number of the bond vectors $\mathbf{r}_{0a}$ and $\mathbf{r}_{0b}$ with equivalent deformations, which contribute to the sums in Eq. (II.1). Obviously, if $\mathbf{r}_{i}$ does not contribute to the sums in Eq. (II.1), then $q_{i}=0$, and if $-\mathbf{r}_{i}$ contributes to the sums in (II.1) instead of $\mathbf{r}_{i}$, then the angle $\theta_{0i}$ of the above equation should be replaced by $\theta_{0i}+\pi$, which changes the sign of both $\cos\theta_{0i}$ and $\sin\theta_{0i}$. This sign change can be absorbed in $q_{i}$, and therefore, the constrain of our case has the form

[TABLE]

As one can see, the deformation energy in Eq. (6) is expressed as a function of the deformations $\delta l_{i}$ and $\delta\phi_{ij}$, while the constrain in Eq. (14) is expressed as a function of $\delta l_{i}$ and $\delta\theta_{i}$. As we show in the Sec. A,

[TABLE]

and therefore, the function $\Lambda$ in Eq. (7), which has to be minimized, can be rewritten as

[TABLE]

where by $\{\delta l_{i}\}$ and $\{\delta\theta_{i}\}$ we denote all the $\delta l_{i}$ and $\delta\theta_{i}$ independent variables, respectively, (i.e. $\{\delta l_{i}\}=\delta l_{1},\delta l_{2},\ldots$ and $\{\delta\theta_{i}\}=\delta\theta_{1},\delta\theta_{2},\ldots$), and therefore $\Lambda$ becomes a function of only $\delta l_{i}$, $\delta\theta_{i}$ and $\lambda$.

It is worth noting that the projection of $\delta\mathbf{L}=\mathbf{L}-\mathbf{L}_{0}$ normal to the strain direction should be zero, i.e. (according to Eq. (10))

[TABLE]

As we will see, minimizing $\Lambda$ in (II.2) we will be able to calculate the differences of $\delta\theta_{i}$ for the same atom, (i.e. the bond angle deformations $\delta\phi_{ij}$), but not the deformations $\delta\theta_{i}$ themselves, which give the direction of the bonds with respect to the strain direction. However, using (17) and the results of the minimization in (II.2), the deformations $\delta\theta_{i}$ can be also determined and we can have a complete figure for the deformations of the structure.

III Minimization of $\Lambda(\{\delta l_{i}\},\{\delta\theta_{ij}\},\lambda)$

The steady state of $\Lambda$ occurs at the specific $\delta l_{i}$ and $\delta\theta_{i}$ values for which

[TABLE]

$\delta l_{i}$ appears only in one term of $U$, namely in $(1/2)k_{s,i}\delta l_{i}^{2}$. Consequently, from $\partial\Lambda/\partial\delta l_{i}=0$ we obtain

[TABLE]

On the other hand, $\delta\theta_{i}$ appears in 4 terms of $U$ (see Fig. 1(a)), namely in $m_{ij_{1}}k_{b,ij_{1}}(\delta\theta_{i}-\delta\theta_{j_{1}})^{2}$ and $m_{ij_{2}}k_{b,ij_{2}}(\delta\theta_{i}-\delta\theta_{j_{2}})^{2}$ for the angles $\delta\phi_{ij_{1}}$ and $\delta\phi_{ij_{2}}$ of atom A, and $m_{ij_{3}}k_{b,ij_{3}}(\delta\theta_{i}-\delta\theta_{j_{3}})^{2}$ and $m_{ij_{4}}k_{b,ij_{4}}(\delta\theta_{i}-\delta\theta_{j_{4}})^{2}$ for the angles $\delta\phi_{ij_{3}}$ and $\delta\phi_{ij_{4}}$ of atom B. From $\partial\Lambda/\partial\delta\theta_{i}=0$ we obtain the linear system

[TABLE]

Substituting the expressions for $\delta\theta_{i}$ obtained from Eq. (20) and the expressions for $\delta l_{i}$ shown in Eq. (19) into (14), we obtain an equation for $\lambda$. Solving this equation with respect to $\lambda$, we obtain $\lambda$ as a function of the strain $\varepsilon$ and the strain angle $\theta_{0}$.

As we show in the Section B,

[TABLE]

where $U_{min}$ is the minimum of $U$ subject to the constrain $\varepsilon=\delta L/L_{0}$. Thus, if $\lambda$ is determined, then $U_{min}$ can also be determined. Eq. (21) gives a physical meaning in the Lagrange multiplier $\lambda$ and minimizes the effort to find a convenient expression for $U_{min}$ as a function of $k_{s,i}$ and $k_{b,ij}$ for strain $\varepsilon$.

IV Including second nearest neighbor stretching terms

As we can see from Eqs. (19) and (20), the original stick and spiral model, expressed utilizing (6), does not provide any coupling between $\delta l_{i}$ and $\delta\phi_{ij}$. However, as already mentioned, including energy terms which describe stretching from second nearest neighbor interactions, we obtain a more accurate model, since it provides coupling between $\delta l_{i}$ and $\delta\phi_{ij}$.

Let us assume that atoms B and C are second nearest neighbors, forming bonds $i$ and $j$, respectively, with atom A. If $\mathbf{r}_{0i}$ and $\mathbf{r}_{0j}$ are the bond vectors of bonds $i$ and $j$, at equilibrium for $\varepsilon=0$, then, depending on the orientation of $\mathbf{r}_{0i}$ and $\mathbf{r}_{0j}$, the interatomic distance $r_{0ij}$ between atoms B and C is either the magnitude of the vector $\mathbf{r}_{0j}-\mathbf{r}_{0i}$ (if both heads or tails of $\mathbf{r}_{0i}$ and $\mathbf{r}_{0j}$ are at the position of atom A), or the vector $\mathbf{r}_{0j}+\mathbf{r}_{0i}$ (if the tail of the one and the head of the other are at the position of atom A).

If the interatomic distance $r_{0ij}$ is deformed upon strain by $\delta r_{ij}$, then the deformation energy per unit cell $U$ is

[TABLE]

where $U_{1}$ is the deformation energy of the original stick and spiral model in Eq. (6) and $U_{2}$ describes the contribution due to stretching deformations of second nearest neighbor interatomic distances. The factor $1/2$ in the second term of Eq. (22) is inserted to avoid double counting, the notation $i$ and $j$ is the same as in (6) and $p_{ij}$ is the number of the equivalent second nearest neighbor interatomic distances in the unit cell with a $\delta r_{ij}$ deformation. Obviously, $p_{ij}=m_{ij}$, because each specific bond angle $\phi_{ij}$ corresponds to a specific second nearest neighbor interatomic distance $r_{ij}$.

Consequently, for the atomic arrangement shown in Fig. 1(a), Eqs. (19) and (20) should be replaced by

[TABLE]

and

[TABLE]

which have to be solved.

As we show in the Sec. C,

[TABLE]

and

[TABLE]

[TABLE]

The upper signs, (wherever $\pm$ and $\mp$ appear), occur when $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ have their tails (or their heads) at the position of the same atom and the lower signs, when the tail of the one and the head of the other are at the position of the same atom, as explained in Sec. C.

Obviously, if $k_{s,ij}=0$, then $U_{2}=0$ and the modified stick and spiral model reduces to the original one. Thus, we can treat both models by solving the system of Eqs. (23) and (IV) of the modified model. Then, by setting $k_{s,ij}=0$ in these solutions, we directly get the solutions of (19) and (20) of the original model. This is the subject of the next section specified for graphene.

V Application to graphene

Bellow, as well as in the appendices, whenever the indices $i^{\prime}$, $j^{\prime}$ and $k^{\prime}$ are used, $(i^{\prime},j^{\prime},k^{\prime})=(1,2,3)$, or $(2,3,1)$, or $(3,1,2)$.

V.1 The energy

Fig. 3 shows the unit cell of graphene, which is defined by the lattice vectors $\mathbf{a}=(\sqrt{3}/2)(\sqrt{3}\hat{\mathbf{i}}+\hat{\mathbf{j}})a_{0}$ and $\mathbf{b}=(\sqrt{3}/2)(\sqrt{3}\hat{\mathbf{i}}-\hat{\mathbf{j}})a_{0}$, where $a_{0}$ is the bond length of graphene. In this figure, A and B are the 2 atoms of the lattice base. As one can see, there are 3 bonds per unit cell, which can be deformed independently, corresponding to the bond vectors $\mathbf{r}_{1}$, $\mathbf{r}_{2}$ and $\mathbf{r}_{3}$ of atom A, or the bond vectors $\mathbf{r}_{4}=\mathbf{r}_{1}$, $\mathbf{r}_{5}=\mathbf{r}_{2}$ and $\mathbf{r}_{6}=\mathbf{r}_{3}$ of atom B. Consequently, in Eqs. (6) and (22), $n_{i}=1$ ($i=1,2,3$). Moreover, as one can see in Fig. 3, there are six bond angles (with respect to the strain direction) $\theta_{i}$ per unit cell. Three of them correspond to atom A and three to atom B. Since the bond vectors of atom A and B are the same, the angles $\theta_{i}$ corresponding to the bonds of atom A are the same with those corresponding to atom B. Consequently, only three of those six angles can be considered as independently deformed, and $m_{ij}=2$. Moreover, due to symmetry reasons, $k_{s,i}=k_{s1}$, $k_{s,ij}=k_{s2}$ and $k_{b,ij}=k_{b}$.

Thus, the energy per unit cell in the original stick and spiral model (according to Eq. (6)) is

[TABLE]

where $k_{b}^{\prime}=k_{b}/a_{0}^{2}$.

In the unit cell of graphene shown in Fig. 3, there are six second nearest neighbor interatomic distances, namely $r_{12}$, $r_{23}$, $r_{31}$, $r_{45}$, $r_{56}$ and $r_{64}$, where $r_{45}=r_{12}$, $r_{56}=r_{23}$ and $r_{64}=r_{31}$. Consequently, there are only three second nearest neighbor interatomic distances, which can be deformed independently and $U_{2}$ in Eq. (22) is

[TABLE]

where $\delta r_{ij}$ are given by (25), and therefore, the energy per atom $U$ in the modified model is $U=U_{1}+U_{2}$.

V.2 The strain constrain

As a function of the independently deformed bond vectors $\mathbf{r}_{i}$, the unit cell vectors $\mathbf{a}$ and $\mathbf{b}$ can be written as

[TABLE]

Thus, if $\mathbf{L}_{0}=n\mathbf{a}+m\mathbf{b}$ defines the strain direction, then $\mathbf{L}_{0}=(n+m)\mathbf{r}_{3}-n\mathbf{r}_{2}-m\mathbf{r}_{1}$, and consequently the $q_{i}$s in (14) are $q_{3}=n+m$, $q_{2}=-n$ and $q_{1}=-m$. As we show in the Sec. D,

[TABLE]

where

[TABLE]

and consequently, (as shown in the same Appendix), the strain constraint of Eq. (14) takes the form

[TABLE]

while (17) becomes

[TABLE]

respectively, where $i=1$, or 2, or 3.

V.3 Solving for the deformations $\delta l_{i}$ and $\delta\theta_{i}$

As we show in the Sec. (E), Eqs. (23) and (IV) give

[TABLE]

and

[TABLE]

The solution of these equations, (as shown in the same appendix), is of the form

[TABLE]

and

[TABLE]

where $\xi_{1}^{\prime}=8k_{b}^{\prime}\lambda/(9a_{0}^{2}K^{\prime})$, $\xi_{2}^{\prime}=k_{s2}\lambda(k_{s1}-18k_{b}^{\prime})/[9a_{0}^{2}K^{\prime}(k_{s1}+6k_{s2})]$, $\xi_{3}^{\prime}=2\lambda k_{s1}/(9a_{0}^{2}K^{\prime})$ and $K^{\prime}=k_{s1}k_{s2}+(4k_{s1}+6k_{s2})k_{b}^{\prime}$. For these expressions of $\delta l_{i}$ and $\delta\theta_{j}-\delta\theta_{i}$, Eqs. (33) and (34) yield

[TABLE]

(see Sec. E for details). Consequently,

[TABLE]

where $K_{0}=k_{s1}^{2}+9k_{s1}k_{s2}+18(k_{s1}+3k_{s2})k_{b}^{\prime}$ and therefore,

[TABLE]

Thus,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

Using Eq. (116), (44) gives

[TABLE]

Obviously, Eq. (39) leads to

[TABLE]

The former shows that $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ are not independent.

Moreover, according to the relations between $\phi_{ij}$ and $\theta_{i}$ shown in Sec. A, the relations between the $\phi_{ij}$ and $\theta_{i}$ angles of graphene, shown in Fig. 3 are

[TABLE]

Thus, the bond angle deformations $\delta\phi_{ij}$ are

[TABLE]

Due to the symmetry of the unit cell, the results we find for strain angle $\theta_{0}$, will be the same for strain angles $n\pi/3\pm\theta_{0}$, $n=0,1,2,3,4,5$. Thus, without loss of generality, we may assume that $0\leq\theta_{0}\leq\pi/6$.

V.4 Energy, Young’s modulus and Poisson’s ratio

According to Eq. (21), the deformation energy per unit cell is $U=\lambda\varepsilon/2$. For graphene, $\lambda$ is given by (41), and consequently,

[TABLE]

where

[TABLE]

As for the Young’s modulus $E$, it is easy to show that $E=2U/(V\varepsilon^{2})$, where $V$ is the volume of the unit cell ($V=3\sqrt{3}a_{0}^{2}d_{0}/2$) and $d_{0}$ is the hypothetical depth of the graphene layer, which is assumed to be equal to the graphite interlayer separation ($d_{0}=3.34$Å), in order to direct compare the Young’s modulus values of two dimensional (2D) carbon structures with the known values for three dimensional (3D) systems, like graphite PCCP_Fthenakis . Thus, for the above expression for $A$,

[TABLE]

Moreover, in Sec. F we show that the Poisson’s ratio $\nu$ is

[TABLE]

which for the $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ expressions of (45), (46) and (47) becomes

[TABLE]

As one can see from the above expressions, $U$, $E$ and $\nu$ are independent of the strain angle $\theta_{0}$, and consequently, graphene is isotropic.

V.5 Relations between $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$

with $\xi_{1}$, $\xi_{2}$, $\xi_{3}$ and $A$

One would have thought that Eqs. (45), (46) and (47), which form a $3\times 3$ system of equations, would provide solutions for $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$ as functions of $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$. However, as shown in Eq. (49), $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ are not independent, and therefore, these equations can not provide relations for $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$ as functions of $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$. On the other hand, $A$, which is independent of $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$, is also a function of $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$. Therefore, $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$ could be written as functions of $\xi_{1}$, $\xi_{2}$, $\xi_{3}$ and $A$.

As we show in the Sec. G,

[TABLE]

and

[TABLE]

[TABLE]

and

[TABLE]

V.6 The original stick and spiral model

The corresponding results for the original stick and spiral model (i.e. not including second nearest neighbor interactions for stretching) can be obtained by setting $k_{s2}=0$. Thus, the solution of Eqs. (19) and (20) have again the form of (42), with $\lambda_{i}$ and $\mu_{ij}$ given again by Eqs. (43) and (44), but now

[TABLE]

The first of the Eqs. (49) becomes $9\xi_{1}+2\xi_{3}=4$, while the second remains the same. The energy and the Young’s modulus are again given by (52) and (54), respectively, but now

[TABLE]

and the Poisson’s ratio is

[TABLE]

Moreover, the relations between $k_{s1}$ and $k_{b}^{\prime}$, with $\xi_{1}$, $\xi_{3}$ and $A$ are

[TABLE]

VI Force constants from DFT results and discussion

VI.1 Details of our DFT calculations

For our DFT calculations we used the Quantum Espresso quantesp code at the level of GGA/PBE functional functional and adopted an ultra-soft pseudopotential for Carbon rrkj ; pseudopotential . The two unit cells are shown in Fig. 5. For the rectangular unit cell of Fig. 5(a) we used a 12$\times$12 k-point mesh, while for the unit cell of Fig. 5(b) a 12$\times$6 (12 along the small real space direction). In addition, we used cut-offs 50 and 500 Ryd for the wave functions and charge density, respectively, and occupation smearing of 5 mRyd. As in Ref. PCCP_Fthenakis, , for non zero uniaxial strain, the unit cells were extended in the strain direction while all the atoms in the cell as well as the vertical cell dimension were fully relaxed.

VI.2 Results

As a first step, we want to calculate the parameters $\lambda_{i}$ and $\mu_{ij}$, which depend on the strain direction, as well as $A$, which is independent. To calculate the $\lambda_{i}$ and $\mu_{ij}$ values, we fit the deformations $\delta l_{i}$ and $\delta\phi_{ij}$ in the strain range $[-0.05,0.05]$ to a quadratic form, considering that the coefficient of the linear term represent the corresponding $3a_{0}\lambda_{i}$ and $\mu_{ij}$ values in Eq. (42), respectively. For the calculation of $A$, we fit the corresponding energy per atom values to a fourth order polynomial, considering that $(3a_{0})^{2}A$ is the coefficient of the quadratic term.

Although in real world, graphene sheet bends for negative strains, computationally it is possible to perform calculations for negative strains without bending of the structure. Fitting a curve to the deformations $\delta l_{i}$, $\delta\phi_{ij}$ and $U$ for both negative and the positive strain values, we expect a better estimation of $\lambda_{i}$, $\mu_{ij}$ and $A$ values, than using an extrapolation of $\delta l_{i}$, $\delta\phi_{ij}$ and $U$ at $\varepsilon=0$, which can be obtained from a fitting of the deformation values of $\delta l_{i}$, $\delta\phi_{ij}$ and $U$ for positive strain values only.

Using the DFT method presented above, we calculated the deformations $\delta l_{i}$ and $\delta\phi_{ij}$, $i,j=1,2,3$, and the deformation energy per atom $U$, for uniaxial strain along the high symmetry arm chair and zig-zag directions, as well as the directions along the vectors $\mathbf{L}=2\mathbf{a}+\mathbf{b}$ and $\mathbf{L}_{\perp}=4\mathbf{a}-5\mathbf{b}$, which are perpendicular to each other, and randomly selected. We increase the strain gradually with a 0.01 strain step in the range between $\varepsilon=-0.1$ and $\varepsilon=0.25$. The results are presented in Figs. 4 and 6, respectively. The fitting functions are presented in the Supplementary Data.

The values of $\lambda_{i}$ and $\mu_{ij}$ obtained from the fits for the four strain directions are presented in Table 1, while the corresponding $A$ values are shown in the legends of Fig. 6. Although $A$ was expected to be independent of the strain direction, the values of $A$ shown in Fig. 6 does not seem to agree with this prediction. However, this discrepancy is due to numerical errors introduced from the different unit cells used. The total energy per atom difference between the equilibrium graphene geometries at $\varepsilon=0$ obtained using the two unit cells of Fig. 5 is $2.3\times 10^{-4}$ eV/atom. As one can show, this difference is enough to produce such a discrepancy in $A$, (i.e. of the order of $10^{-3}$ eV/Å2). It is worth noting, however, that the difference between the two $A$ values, corresponding to the two perpendicular strain directions of the same unit cell, is of the order of $10^{-4}$ eV/Å2. For our calculations we will adopt the value $A=3.046$ eV/Å2, which corresponds to an average of the obtained values.

The second step is to calculate the values of $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ using the $\lambda_{i}$ and $\mu_{ij}$ values of Table 1 and Eqs. (43) and (44). According to these equations, $\xi_{1}$, $\xi_{2}$ and $\xi_{3}$ can be obtained using a linear fitting of the $\lambda_{i}$ values as a function of $\cos^{2}\theta_{0i}$ and the $-\mu_{i^{\prime}j^{\prime}}/\sqrt{3}$ values as a function of $\cos 2\theta_{0k^{\prime}}$. The values of $\lambda_{i}$ as a function of $\cos^{2}\theta_{0i}$ and the values of $-\mu_{i^{\prime}j^{\prime}}/\sqrt{3}$ as a function of $\cos 2\theta_{0k^{\prime}}$, as well as the corresponding fitting lines are shown in Fig. 7(a). The smoothness of the fitting is obvious. These fitting lines are

[TABLE]

and

[TABLE]

Thus, $\xi_{1}=0.278921$, $\xi_{2}=-0.002272$ and $\xi_{3}=0.758145$. Using these values, the value of $A$, and Eqs. (57) - (60), we can calculate the values of $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$, as well as the ratios $k_{s2}/k_{s1}$ and $k_{b}^{\prime}/k_{s1}$. Thus, $k_{b}^{\prime}/k_{s1}=0.091975$, $k_{s2}/k_{s1}=0.0096665$, $k_{s1}=41.972$ eV/Å, $k_{s2}=0.40572$ eV/Å and $k_{b}^{\prime}=3.8604$ eV/Å. Therefore, roughly speaking $k_{b}^{\prime}\approx 0.1k_{s1}$ and $k_{s2}\approx 0.01k_{s1}$, which qualitatively provides the relative strength of each deformation mode. Moreover, according to (54) and (55), $E=1012$ GPa and $\nu=0.1744$, in agreement with the results of our previous work PCCP_Fthenakis obtained fitting the stress $\sigma$ and the the transverse strain $\varepsilon_{\perp}$ values as a function of strain, to a third and second order polynomial, respectively.

Knowing the $k_{s1}$, $k_{s2}$ and $k_{b}^{\prime}$ values, we have the ability to predict any mechanical property related to the in-plane deformations of graphene and not only $E$ and $\nu$. For instance, the corresponding biaxial isotropic modulus $E_{B}=\sigma/\varepsilon$, where $\sigma=\sigma_{xx}=\sigma_{yy}$ and $\varepsilon=\varepsilon_{xx}=\varepsilon_{yy}$, is $E_{B}=4\sqrt{3}A^{\prime}/d_{0}$, where for the biaxial isotropic deformation $U=9a_{0}^{2}A^{\prime}\varepsilon^{2}$. Using (28) and (29), it is easy to show that for biaxial isotropic strain $A^{\prime}=k_{s1}/6+k_{s2}$. Thus, for graphene, $E_{B}=2459$ GPa. A different calculation using the relation $U=k\delta l^{2}/2+k(\delta l+\delta l_{\perp})^{2}/2=k\delta l^{2}(1+\nu+\nu^{2}/2)=2U_{u}(1+\nu+\nu^{2}/2)$, or $A^{\prime}=2(1+\nu+\nu^{2}/2)A$, yields $E_{B}=2408$ GPa. As one can see, the two results are very close to each other.

Obviously, the term $U_{2}$ corresponding to the stretching of the second nearest neighbor interatomic distances is the less important energy contribution, but it is not a term that can be ignored. If this term is ignored, (which is equivalent to set $k_{s2}=0$ or $\xi_{2}=0$), the energy model reduces to the original stick and spiral model, which, according to (19), predicts that any bond which is perpendicular to the strain direction remains undeformed. This, however, is in contrast to what we find from our DFT calculations for the $l_{3}$ bond length under uniaxial strain along the zig-zag direction. Just for comparison, we also calculate the corresponding $\xi_{1}$, $\xi_{3}$, $k_{s1}$ and $k_{b}^{\prime}$ values obtained from the original stick and spiral morel. Obviously, the form of Eq. (44) does not change in the original stick and spiral model and consequently the value of $\xi_{3}$ remains the same as the modified model. However, (43) becomes $\lambda_{i}=\xi_{1}\cos^{2}\theta_{0i}$. The corresponding fit for the $\lambda_{i}$ values of Table 1 as a function of $\cos^{2}\theta_{0i}$ yelds $\xi_{1}=0.275981$. In Fig. 7(b) we show the prediction error $\delta\lambda_{i}$ (i.e. the difference between the $\lambda_{i}$ provided by the fitting equations of $\lambda_{i}$ as a function of $\cos^{2}\theta_{0i}$ and the corresponding $\lambda_{i}$ values of Table 1 for the original and the modified stick and spiral model. As we can see, the error for the modified sick and spiral model is between $\pm 0.001$, while the error for the original model is almost double, ranging between -0.0025 and 0.0017. The values of $k_{s1}$ and $k_{b}^{\prime}$ for the original model, according to (64) are $k_{s1}=44.178$ eV/Å and $k_{b}^{\prime}=4.0177$ eV/Å, i.e. they are overesimated by 5 and 4%, respectively, in comparisson with the corresponding values obtained from the modified model. Thus, the original stick and spiral model can not provide an accurate description for the bond and angle deformations of graphene, or at least, it can not provide such an accurate description as the modified model, which is presented here.

VII Conclusions

In summary, we present a method for the study of the equilibrium deformations of 2D planar materials under uniaxial strain. The method is based on the stick and spiral model including angle bending energy terms and either only 1st nearest neighbors bond stretching terms (case 1) or both 1st and 2nd nearest neighbors terms (case 2). The method can be generalized to describe structural deformations not only under uniaxial strain, but also under any strain conditions. We present analytic expressions/equations for the structure deformations under strain, namely the equilibrium angle bending and bond stretching deformations for both case 1 (equations (19) and (20)) and case 2 (equations (23) and (IV)). We then focus on graphene in order to assess the applicability of our method for which we perform DFT calculations for several values of strain in 4 different directions. We find that the original stick and spiral model (case 1) decouples the equations yielding $\delta l_{i}$ from those yielding $\delta\theta_{i}$ and for graphene, it predicts that the vertical to the strain bonds are not modified. This is in contrast with the DFT results. The inclusion of 2nd nearest neighbors stretching terms (case 2) results in the coupling of $\delta l_{i}$ and $\delta\theta_{i}$, improves the model significantly and brings the results in close agreement with DFT. Our method provides a simple and solid method to study the structural deformations of Graphene in the case of uniaxial strain on any direction in the elastic regime. The elastic properties of graphene under strain are very accurately reproduced by our method. Although this first application concerns graphene, our method can be applied to any 2D planar material and it would be interesting to assess its accuracy on different structures and materials like Graphene planar allotropes, h-BN, Si3B, Si2BN, CdS, etc.

Acknowledgements

NNL acknowledges support from the Hellenic Ministry of Education (through ESPA) and from the GSRT through “Advanced Materials and Devices” program (MIS:5002409).

Appendix A Relation between $\phi_{ij}$ and $\theta_{i}$s

Let us define, for each atom of the unit cell, a local anti-clockwise frame of coordinates with its origin at the position of that atom and its x-axis along the strain direction, as shown in Fig. 1(c). Let us denote as $\mathbf{r}_{1}$, $\mathbf{r}_{2}$ and $\mathbf{r}_{3}$ the three bond vectors, which have their tail on atom $i$ and by $\theta_{1}$, $\theta_{2}$ and $\theta_{3}$ the corresponding angles between these bond vectors with the strain direction, respectively, as shown in Fig. 1(c).

Obviously, $\mathbf{r_{i}}\mathbf{r_{j}}=r_{i}r_{j}\cos\phi_{ij}$, where $\phi_{ij}$ is the angle formed by the bonds $i$ and $j$, and $\mathbf{r}_{i}=r_{i}\cos\theta_{i}\hat{\mathbf{i}}+r_{i}\sin\theta_{i}\hat{\mathbf{j}}$, $i=1,2,3$. Thus, the dot product $\mathbf{r_{i}}\mathbf{r_{j}}$ can be written as

[TABLE]

and consequently,

[TABLE]

If $\phi_{0ij}$, $\theta_{0i}$ and $\theta_{0j}$ are the values of the corresponding $\phi_{ij}$, $\theta_{i}$ and $\theta_{j}$ angles at equilibrium for $\varepsilon=0$, then using a first order Taylor expansion around these values, Eq. (68) yields

[TABLE]

where $\phi_{ij}=\phi_{0ij}+\delta\phi_{ij}$, $\theta_{i}=\theta_{0i}+\delta\theta_{i}$ and $\theta_{j}=\theta_{0j}+\delta\theta_{j}$ are the corresponding angles at $\varepsilon\neq 0$. Thus, the derivative of $\delta\phi_{ij}$ with respect to $\delta\theta_{i}$ is

[TABLE]

Imposing that $0<\phi_{ij}\leq\pi$, (68) gives

[TABLE]

If $\theta_{i}$s, $i=1,2,3$ are defined inside the same unit circle (e.g. $0\leq\theta_{i}<2\pi$ or $-\pi<\theta_{i}\leq\pi$), then $-2\pi<\theta_{i}-\theta_{j}<2\pi$. However, according to (71), $\theta_{i}-\theta_{j}$ is out of the range $(-2\pi,2\pi)$, for $k\neq 0$ or $1$, and therefore only $k=0$ and $k=1$ should be considered. Consequently, (i) for $k=0$ (or $0<|\theta_{i}-\theta_{j}|\leq\pi$, according to (71)), $\phi_{ij}=|\theta_{i}-\theta_{j}|$ and (ii) for $k=1$ (or $\pi\leq|\theta_{i}-\theta_{j}|<2\pi$, according to (71)), $\phi_{ij}=2\pi-|\theta_{i}-\theta_{j}|$. Thus, for any case, $\delta\phi_{ij}=\pm(\delta\theta_{i}-\delta\theta_{j})$, which leads to (15).

If $\mathbf{r}_{i}$s, $i=1,2,3$, have their tail at the position of an atom A, then they have their head at the position of the atoms which form bonds with atom A. Assume B is such an atom, which forms a bond with another atom C (different than A), and $\mathbf{r}_{1}$ and $\mathbf{r}_{4}$ are the bond vectors corresponding to the bonds A-B and B-C, respectively. There are two options for the direction of $\mathbf{r}_{4}$: either its head is on the position of atom B and its tail on the position of atom C, or the opposite. In the former case, the relations between the bond angle $\phi_{ij}$ and the bond angle $\theta_{i}$ with respect to the strain direction are the same with those presented above, since $\mathbf{r}_{1}\mathbf{r}_{4}=r_{1}r_{4}\cos\phi_{14}$. However, in the later case, $\mathbf{r}_{1}\mathbf{r}_{4}=r_{1}r_{4}\cos\omega_{14}$, where the bond angle $\phi_{14}$ is $\phi_{14}=\pi-\omega_{14}$. Thus, for this case, the relations presented above will be valid if $\phi_{ij}$ is replaced by $\pi-\phi_{ij}$. Thus, (68), should be replaced by

[TABLE]

[TABLE]

and

[TABLE]

If $0\leq\pi-\phi_{ij}<\pi$, then $0<\phi_{ij}\leq\pi$. For $\phi_{ij}$ in this range, (72) yields (i) if $0<|\theta_{i}-\theta_{j}|\leq\pi$, then $\phi_{ij}=\pi-|\theta_{i}-\theta_{j}|$ and (ii) if $\pi<|\theta_{i}-\theta_{j}|\leq 2\pi$, then $\phi_{ij}=|\theta_{i}-\theta_{j}|-\pi$. Obviously, therefore, for this case, $\delta\phi_{ij}$ is also $\delta\phi_{ij}=\pm(\delta\theta_{i}-\delta\theta_{j})$ and consequently, (15) is also valid.

Appendix B The physical meaning of $\lambda$

Obviously, $\Lambda$ is parametrically dependent on $\varepsilon$, i.e. $\Lambda=\Lambda(\{\delta l_{i}\},\{\delta\theta_{i}\},\lambda;\varepsilon)$. If $\Lambda$ is minimized for $\delta l_{i}=\delta l_{i}^{*}$, $\delta\theta_{i}=\delta\theta_{i}^{*}$, and $\lambda=\lambda^{*}$, where $\delta l_{i}^{*}$s, $\delta\theta_{i}^{*}$s and $\lambda^{*}$ are specific values of $\delta l_{i}$s, $\delta\theta_{i}$s and $\lambda$, respectively, then $\Lambda_{min}=\Lambda(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\},\lambda^{*};\varepsilon)=\Lambda_{min}(\varepsilon)$, where $\Lambda_{min}$ is the minimum of $\Lambda$.

For $\delta l_{i}=\delta l_{i}^{*}$ and $\delta\theta_{i}=\delta\theta_{i}^{*}$, the strain $\varepsilon$ is $\varepsilon=\delta L(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\})/L_{0}$ and $U$ is minimized subject to the constrain $\varepsilon=\delta L/L_{0}$. Thus, if $U_{min}$ is the minimum of $U$ subject to the constrain $\varepsilon=\delta L/L_{0}$, then $U_{min}=U(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\})$ and (according to (14)), $U(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\})=\Lambda(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\},\lambda^{*};\varepsilon)$, or $U_{min}(\varepsilon)=\Lambda_{min}(\varepsilon)$.

According to (19) and (20), for the minimized $\Lambda$, $\delta l_{i}^{*}$ and $\delta\theta_{i}^{*}$ depend linearly on $\lambda^{*}$, and therefore, according to (14), $\lambda^{*}$ should depend linearly on $\varepsilon$. Thus, $\delta l_{i}^{*}=\delta l_{i}^{*}(\varepsilon)$ and $\delta\theta_{i}^{*}=\delta\theta_{i}^{*}(\varepsilon)$, and consequently, $U_{min}=U_{min}(\varepsilon)$. On the other hand, $U_{min}$ is quadratically dependent on $\delta l_{i}^{*}$ and $\delta\theta_{i}^{*}$, and consequently $U_{min}$ should depend quadratically on $\varepsilon$. Therefore we can write $U_{min}(\varepsilon)=K\varepsilon^{2}$, where $K=K(\{k_{si}\},\{k_{bij}\})$.

Obviously, $\partial\Lambda(\{\delta l_{i}\},\{\delta\theta_{i}\},\lambda;\varepsilon)/\partial\varepsilon=\lambda$, and consequently, $d\Lambda_{min}/d\varepsilon=d\Lambda(\{\delta l_{i}^{*}\},\{\delta\theta_{i}^{*}\},\lambda^{*};\varepsilon)/d\varepsilon=\lambda^{*}$. Since, $\Lambda_{min}=U_{min}$, we have $d\Lambda_{min}/d\varepsilon=dU_{min}/d\varepsilon=2K\varepsilon$. Thus, $2K\varepsilon=\lambda^{*}$, which leads to (21).

Appendix C $\delta r_{ij}$ as a function of bond length and bond angle deformations

Let us assume that atoms A, B and C belong to the same planar 2D structure and atom A forms bonds with atoms B and C. Let us also assume that $\mathbf{r}_{0i}$ and $\mathbf{r}_{0j}$ are the bond vectors corresponding to the bonds A-B and A-C at equilibrium for $\varepsilon=0$, having both their tails (or their heads) at the position of atom A. Then the interatomic distance $r_{0ij}$ between atoms B and C is the length of the vector $\mathbf{r}_{0ij}=\mathbf{r}_{0j}-\mathbf{r}_{0i}$, for which

[TABLE]

where $l_{0i}$ and $l_{0j}$ are the lengths of $\mathbf{r}_{0i}$ and $\mathbf{r}_{0j}$, respectively, and $\phi_{0ij}$ the bond angle between bonds A-B and A-C. If at the equilibrium state under strain, $l_{0i}$, $l_{0j}$, $r_{0ij}$ and $\phi_{0ij}$ are deformed to $l_{i}=l_{0i}+\delta l_{i}$, $l_{j}=l_{0j}+\delta l_{j}$, $r_{ij}=r_{0ij}+\delta r_{ij}$ and $\phi_{ij}=\phi_{0ij}+\delta\phi_{ij}$, respectively, then

[TABLE]

For $\delta r_{ij}<<r_{0ij}$, $r_{ij}^{2}\approx r_{0ij}^{2}+2r_{0ij}\delta r_{ij},$ and consequently, (76) leads to

[TABLE]

Therefore, $\delta r_{ij}$ is a function of the deformations of $\delta l_{i}$, $\delta l_{j}$, $\delta\theta_{i}$ and $\delta\theta_{j}$, (see Sec. A).

The derivatives of $\delta r_{ij}$ with respect to $\delta l_{i}$ and $\delta\theta_{i}$ are

[TABLE]

and

[TABLE]

Using (68), (69) and (70) the above equations give

[TABLE]

and

[TABLE]

[TABLE]

However, if the head of $\mathbf{r}_{i}$ and the tail of $\mathbf{r}_{j}$ (or vice versa) are at the position of atom A, then we have to use (72), (73) and (74) instead of (68), (69) and (70) (see Sec. A), and thus, (77), (78) and (79) give

[TABLE]

and

[TABLE]

[TABLE]

Commuting $i$ with $j$ in (78), (79), (81), (82), (84) and (85), we obtain the corresponding relations for $\partial\delta r_{ij}/\partial\delta\theta_{j}$ and $\partial\delta r_{ij}/\partial\delta l_{j}$.

Appendix D Derivation of Eqs. (31), (33) and (34)

If $\mathbf{L}_{0}=n\mathbf{a}+m\mathbf{b}$ defines the strain direction, then $\mathbf{L}_{0}=(\sqrt{3}/2)(\sqrt{3}(n+m)\hat{{\mathbf{i}}}+(n-m)\hat{\mathbf{j}})a_{0}$, and consequently, $\cos\theta_{0}=3(n+m)a_{0}/(2L_{0})$ and $\sin\theta_{0}=\sqrt{3}(n-m)a_{0}/(2L_{0})$, where $\theta_{0}$ is the angle of the strain direction with respect to the x-axis. Solving these two equations with respect to $n$ and $m$, we obtain, $n=2L_{0}/(3a_{0})((1/2)\cos\theta_{0}+(\sqrt{3}/2)\sin\theta_{0})=-2L_{0}/(3a_{0})\cos\theta_{02}$ and $m=2L_{0}/(3a_{0})((1/2)\cos\theta_{0}-(\sqrt{3}/2)\sin\theta_{0})=-2L_{0}/(3a_{0})\cos\theta_{01}$, and consequently, $n+m=2L_{0}/(3a_{0})\cos\theta_{0}=2L_{0}/(3a_{0})\cos\theta_{03}$, which lead to (31). In Sec. H we present useful relations between the trigonometric functions of these angles, which will be used here.

Bearing in mind that in graphene $l_{01}=l_{02}=l_{03}=a_{0}$, and using (31), (14) becomes

[TABLE]

Using (105) for $k=2$ of Sec. H, the above equation leads to Eq. (33).

Moreover, (17) becomes

[TABLE]

which leads to (34). In the last step of the above equation we used (H) of the Sec. H.

Appendix E Derivation of Eqs. (35), (36), (37),

(38) and (39)

As we can see in Fig. 3, the tails of the bond vectors $\mathbf{r}_{1}$, $\mathbf{r}_{2}$ and $\mathbf{r}_{3}$ are at the position of atom A, while the heads of the bond vectors $\mathbf{r}_{4}$, $\mathbf{r}_{5}$ and $\mathbf{r}_{6}$ are at the position of atom B. Therefore, to apply (25), (26) and (27) to (23) and (IV), we have to use the upper signs among $\pm$ and $\mp$. Moreover, $l_{0i}=a_{0}$, $r_{0ij}=\sqrt{3}a_{0}$, $\cos(\theta_{0j^{\prime}}-\theta_{0i^{\prime}})=\cos(2\pi/3)=-1/2$ and $\sin(\theta_{0j^{\prime}}-\theta_{0i^{\prime}})=\sin(2\pi/3)=\sqrt{3}/2$. Consequently, (25), (26) and (27) yield

[TABLE]

$\partial\delta r_{ij}/\partial\delta l_{i}=\sqrt{3}/2$ and $\partial\delta r_{i^{\prime}j^{\prime}}/\partial\delta\theta_{i^{\prime}}=-\partial\delta r_{i^{\prime}j^{\prime}}/\partial\delta\theta_{j^{\prime}}=-a_{0}/2$, respectively. Thus, (23) gives

[TABLE]

which leads to (35), and (IV) gives

[TABLE]

which leads to (36). Summing up the three equations (35) (i.e. for $(i^{\prime},j^{\prime},k^{\prime})=(1,2,3)$, $(2,3,1)$ and $(3,1,2)$), and using (H) we obtain

[TABLE]

Substituting $\sum_{i=1}^{3}\delta l_{i}$ in (35) we take

[TABLE]

Subtracting by parts equations (36) (two at a time) leads to

[TABLE]

The solution of the system of (90) and (91) are (37) and (38).

Using the expressions of (37) and (38) for $\delta l_{i}$ and $\delta\theta_{j}-\delta\theta_{i}$, and (105), (H), (H), (112) and (114), (33) and (34) become

[TABLE]

and

[TABLE]

respectively, leading to (39).

Appendix F Poisson’s ratio

In order to find the Poisson’s Ratio $\nu$, ($\nu=-\varepsilon_{\perp}/\varepsilon$), we need to find the transverse strain $\varepsilon_{\perp}=\delta L_{\perp}/L_{\perp 0}$, where $L_{\perp 0}$ is a length of the material perpendicular to the strain direction and $\delta L_{\perp}$ its deformation upon tensile strain $\varepsilon$. If $\mathbf{L}_{\perp 0}=t_{a}\mathbf{a}+t_{b}\mathbf{b}$ is a lattice vector, which is perpendicular to the vector $\mathbf{L}_{0}=n\mathbf{a}+m\mathbf{b}$, which defines the strain direction, then

[TABLE]

For convenience we may select $t_{a}$ and $t_{b}$ to be $t_{a}=2m+n$ and $t_{b}=-(2n+m)$. Using (30), $\mathbf{L}_{\perp 0}$ becomes $\mathbf{L}_{\perp 0}=(m-n)\mathbf{r}_{1}+(2n+m)\mathbf{r}_{2}-(2m+n)\mathbf{r}_{3}$. The projection of the deformation of a bond vector normal to the strain direction is given by (10). Thus, the deformation $\delta L_{\perp}$ of $\mathbf{L}_{\perp 0}$ is

[TABLE]

where $q_{\perp 1}=m-n=2L_{0}/(3a_{0})(\cos\theta_{03}-\cos\theta_{02})$, $q_{\perp 2}=2n+m=(n+m)+n=2L_{0}/(3a_{0})(\cos\theta_{01}-\cos\theta_{03})$ and $q_{\perp 3}=-(2m+n)=-m-(n+m)=2L_{0}/(3a_{0})(\cos\theta_{02}-\cos\theta_{01})$. Using (117) we have

[TABLE]

and consequently (using (42), (43), (49) and (H))

[TABLE]

The magnitude $L_{\perp 0}$ of the vector $\mathbf{L}_{\perp 0}$ is

[TABLE]

Using (105) and (H)

[TABLE]

Thus,

[TABLE]

and consequently,

[TABLE]

which leads to (55).

Appendix G Derivation of Eqs. (57), (58),

(59), (60)

The first of (57) can be directly obtained if we divide by parts (45) and (47). Using that equation, (46) becomes

[TABLE]

(47) can also be written as

[TABLE]

Dividing (100) and (101) by parts we obtain

[TABLE]

Using the first of (49), this equation leads to the second of (57).

From the expression $K_{0}=k_{s1}^{2}+9k_{s1}k_{s2}+18(k_{s1}+3k_{s2})k_{b}^{\prime}$, it is obvious that the expression $k_{s1}k_{s2}+2(2k_{s1}+3k_{s2})k_{b}^{\prime}$, which appears in (53), is $k_{s1}k_{s2}+2(2k_{s1}+3k_{s2})k_{b}^{\prime}=(K_{0}-k_{s1}^{2}+18k_{s1}k_{b}^{\prime})/9$. Thus, using (100) and the second of (57),

[TABLE]

Thus,

[TABLE]

Using (49) (i.e. $2-\xi_{3}=9\xi_{1}/2+6\xi_{2}$), we find

[TABLE]

Solving this equation with respect to $k_{s1}$ we get (58).

Using the expression in (58) for $k_{s1}$ and (57), the derivation of (59) and (60) is obvious.

Appendix H Useful relations between trigonometric functions of $\theta_{0i}$

of graphene

Some relations, which are used in the present study, between the trigonometric functions of the angles $\theta_{0i}$ defined by (32), are presented here.

As we have already seen in Sec. V.2, $q_{1}=-m$, $q_{2}=-n$ and $q_{3}=n+m$. Thus, $q_{1}+q_{2}+q_{3}=0$, and consequently, $\cos\theta_{01}+\cos\theta_{02}+\cos\theta_{03}=0$, where $\theta_{0i}=\theta_{0i}(\theta_{0})=2\pi i/3-\theta_{0}$, $i=1,2,3$. Obviously, (i) $2\theta_{01}(\theta_{0})=4\pi/3-2\theta_{0}=\theta_{02}(2\theta_{0})$, $2\theta_{02}(\theta_{0})=8\pi/3-2\theta_{0}=2\pi+\theta_{01}(2\theta_{0})$ and $2\theta_{03}(\theta_{0})=4\pi-2\theta_{0}=2\pi+\theta_{03}(2\theta_{0})$, and (ii) $4\theta_{01}(\theta_{0})=8\pi/3-4\theta_{0}=2\pi+\theta_{01}(4\theta_{0})$, $4\theta_{02}(\theta_{0})=16\pi/3-2\theta_{0}=4\pi+\theta_{02}(4\theta_{0})$ and $4\theta_{03}(\theta_{0})=8\pi-2\theta_{0}=6\pi+\theta_{03}(4\theta_{0})$. Consequently, for $k=1$, or 2, or 4,

[TABLE]

The first derivative of the above equation with respect to $\theta_{0}$ gives

[TABLE]

Using (104) for $k=2$ or $k=4$, and the relation $\cos 2\theta=2\cos^{2}\theta-1$ we obtain

[TABLE]

Then, using the relation $\sin^{2}\theta=1-\cos^{2}\theta$, we obtain

[TABLE]

Moreover, using the relation $\sin 2\theta=2\sin\theta\cos\theta$, (105) for $k=2$ yields

[TABLE]

Using (104) for $k=1$ and (H), we obtain

[TABLE]

In turn, using (105) for $k=1$ and (H), we obtain

[TABLE]

Thus,

[TABLE]

Consequently, (H) gives

[TABLE]

and in turn, (H) gives

[TABLE]

Taking the first derivative of (112), we obtain

[TABLE]

Moreover, let us consider the trigonometric identity

[TABLE]

For the angles $\theta_{0i}$ in (32) we have $(\theta_{0j}-\theta_{0i})/2=(j-i)\pi/3$ and $(\theta_{0j}+\theta_{0i})/2=(i+j)\pi/3-\theta_{0}$. For $(i,j,k)=(1,2,3)$, or $(2,3,1)$, or $(3,1,2)$, the sum $i+j+k=6$, and consequently, $i+j=6-k$. Thus, $(i+j)\pi/3-\theta_{0}=(6-k)\pi/3-\theta_{0}=2\pi-k\pi/3-\theta_{0}=2\pi-k\pi+\theta_{0k}$. Consequently, $\cos[m(\theta_{0i}+\theta_{0j})/2]=\cos[m(\theta_{0k}-k\pi)]$ and $\sin[m(\theta_{0j}-\theta_{0i})/2]=\sin[m(j-i)\pi/3]$. Thus, for $m=1$

[TABLE]

and for $m=2$

[TABLE]

Taking the first derivative of (115) with respect to $\theta_{0}$ we obtain

[TABLE]