Constitutive modeling of some 2D crystals: graphene, hexagonal BN, MoS 2 , WSe 2 and NbSe 2
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n onstitutive modeling of some 2D crystals: graphene,hexagonal BN, MoS , WSe and NbSe . D. Sfyris, G.I. Sfyris, C. GaliotisOctober 15, 2018
Abstract
We lay down a nonlinear elastic constitutive framework for the modeling of some2D crystals of current interest. The 2D crystals we treat are graphene, hexagonalboron nitride and some metal dichalcogenides: molybdenium disulfide (MoS ), tung-sten selenium (WSe ), and niobium diselenide (NbSe ). We first find their arithmeticsymmetries by using the theory of monoatomic and diatomic 2-nets. Then, by con-finement to weak transformation neighbourhoods and by applying the Cauchy-Bornrule we are able to use the symmetries continuum mechanics utilizes: geometric sym-metries. We give the complete and irreducible representation for energies dependingon an in-plane measure, the curvature tensor and the shift vector. This is done for thesymmetry hierarchies that describe how symmetry changes at the continuum level: C ν → C ν → C for monoatomic 2-nets and C ν → C ν → C for diatomic two nets.Having these energies at hand we are able to evaluate stresses and couple stresses foreach symmetry regime. Keywords: graphene, hexagonal BN, MoS , WSe , NbSe , nonlinear elasticity. Recently, strictly 2D atomic crystals have been isolated from three dimensional layeredmaterials. Novoselov et al ([17]) report free standing atomic crystals that can be viewedas individual atomic planes pulled out of the bulk crystal. Using micromechanical cleavagethese authors study single layers of graphene, hexagonal boron nitride and some metaldichalcogenides (such as MoS , NbSe , and WSe ).Motivated by this work, we here lay down a constitutive framework for studying such2D crystalline materials suitable for the nonlinear elasticity theory. We view grapheneas a monoatomic 2-net ([26, 27, 28]), while hexagonal BN, MoS , NbSe , and WSe areviewed as diatomic 2-nets. The arithmetic symmetry of such 2-nets is well reported in theworks of Fadda-Zanzotto ([9, 11]) based on the earlier works of Ericksen ([4, 5, 6, 7]), Parry([18, 19, 20]), Pitteri ([21, 22]) on the definition of arithmetic symmetry. The classification2f symmetry in arithmetic classes offers a more stringent classification than crystallogrpahicpoint groups since conjugacy is taken within the general linear group and not within theorthogonal group ([23, 12]).Monoatomic 2-nets consist of two simple Bravais lattices which have indistinguishableatoms. Graphene belongs to this category ([26, 27, 28]), since it is made of two hexago-nal Bravais lattices with carbon atoms occupying atomic positions. On the other hand,diatomic 2-nets consist of two simple Bravais lattices the atomic positions of which areoccupied by different atoms. Boron nitride is an example that belongs to this class: onehexagonal Bravais lattice consist of boron atoms only, while the second lattice consists ofnitride atoms only. The metal dichalcogenides MoS , NbSe , and WSe belong to the samecategory and are also treated in the analysis.We here confine the analysis to weak transformation neighbourhoods ([25]) and use theCauchy-Born rule ([8]) to lay down the complete and irreducible representation ([34, 35]) foran energy depending on three arguments. The first argument is the surface right Cauchy-Green deformation tensor which is a measure of the in-plane deformations of the 2D crystal.The second argument is the curvature tensor, which introduces out-of-plane deformationsand is motivated by the work of Steigmann-Ogden ([31]) and earlier approaches on thetopic ([2, 14, 16]). The third argument that the energy depends on is the shift vector; thisis the vector that relates the two lattice. Within these limits arithmetic and geometricsymmetries become equivalent which is of particular interest since geometric symmetriesare used in continuum mechanics.Having the complete and irreducible representation for the energy, we are able to evalu-ate the surface stress tensor and the surface couple stress tensor. The first being responsiblefor the in-plane motions, the second for the out-of-plane motions. These measures par-ticipate to the momentum and the moment of momentum equations which are the fieldequations for this problem. From the physical point of view, momentum equation is theforce balance for the surface, while moment of momentum equation is the couple balancefor the surface. To these field equations one should add the equation ruling the shift vec-tor. Form the physical point of view this equation says that the shift vector adjust so thatequilibrium is reached.The analysis carries over to cases where symmetry changes for monoatomic 2-netsaccording to the hierarchies C ν → C ν → C . For the diatomic 2-nets symmetry hierarchieschange as C ν → C ν → C . These groups are the geometric symmetry groups whichone derives if the analysis is confined to weak transformation neighbourhoods startingfrom the arithmetic symmetries. The suitable geometric symmetry group is found byevaluating the eigenvalues for matrices of the arithmetic symmetry and corresponding themto appropriate generators of a geometric symmetry group. We lay down the complete andirreducible representation for the energies for these cases without studying what happensto the transition regime. This is a work in progress in line with similar fundamental workson zirconia ([32]).The paper is structured as follows. Section 2 deals with the definition of monoatomicand diatomic 2-nets as well as their symmetries and symmetry hierarchies. Section 3presents the limitations for the proper transition to the classical continuum viewpoint.3ection 4 describes the basic kinematics for a surface energy depending on a surface mea-sure, the curvature tensor and the shift vector. Section 5 desrcibes the way materialsymmetry should be viewed for the present framework and also gives the field equations.Then, Section 6 gives the complete and irreducible representation for the energy undera specific symmetry group. These symmetry groups are the above mentioned geometricsymmetry groups that describe how symmetry breaks for such materials. Surface stressand couple stress tensor can then be evaluated. The article ends up at Section 7 with someconcluding remarks and future directions. The importance in the difference between arithmetic and geometric symmetries for crystalsstem form the fundamental work of Ericksen ([4]). We refer to the book of Pitteri-Zanzotto([25], and references therein) for a nice exposition of the topic.A three dimensional simple lattice in R is defined as (see e.g. [25]) L ( e a ) = { x ∈ R : x = M a e a , a = 1 , , , M a ∈ Z} , (1)with e a being the lattice vectors and Z the space of positive integers. The geometricsymmetry group of L is ([25, 5]) P ( e a ) = { Q ∈ O : L ( Qe a ) = L ( e a )= { Q ∈ O : Qe a = m ba e b , m ∈ GL (3 , Z ) } , (2)where GL (3 , Z ) and O are the general linear, and the orthogonal group, respectively.Essentially, this group gives all orthogonal transformations Q that map L to itself. Hereone finds, for a three dimensional lattice, the 7 crystal systems and the 28 crystallographicpoint groups continuum mechanics utilizes.Due to the fact that conjugacy in GL (3 , Z ) is more stringent than conjugacy in O ,a finer description of symmetry is given by the arithmetic symmetry which is defined as([25, 5]) L ( e a ) = { m ∈ GL (3 , Z ) : m ba e b = Qe a , Q ∈ P ( e a ) } . (3)This group gives all distinct types of lattices that are compatible with a given geometricgroup. Essentially, this is a finer description of symmetry that can differentiate betweenBravais lattice types within the same crystal system.A multilattice is a generalization of a simple lattice in the sense that it is the finiteunion of translates of some suitable simple lattice M ( p i , e a ) = ∪ n − i =0 L ( p i , e a ) . (4)The particular case of a 2-lattice is the union of two simple lattices M ( e a , p ) = L ( e a ) ∪ { p + L ( e a ) } , (5)4 ( e a ) being the simple lattice generated by the basis e a , while the vector p = p e + p e (6)is called the shift vector and gives the separation of the two simple lattices constituting M .The geometric and arithmetic symmetry groups of multilattices are defined in a similarfashion with the corresponding definitions of simple lattices and we refer to [24, 9, 10, 12,11, 25] for more informations.A particular class of a multilattice is the case of monoatomic 2-lattice (net). In this case,two simple lattices constitute the 2-net and atoms at lattice points are indistinguishable,in the sense that they belong to the same species. For the 2 dimensional case the unitcell of monoatomic 2-nets are depicted in Figure 1 ([9]). There are five Bravais lattice forFigure 1: Bravais lattices for monoatomic 2-net and their symmetry hierarchies (Figuretaken from [9]). Atoms at lattice points of both lattices belong to the same species.monoatomic 2-nets: (1) oblique, (2) side-rectangular, (3) axis-rectangular, (4) rhombic and(5) hexagonal. 5raphene is a monoatomic 2-lattice consisting of carbon atoms in all lattice positionsthat also has a hexagonal unit cell with a shift vector that is depicted in Figure 2. WhenFigure 2: A schematic representation of the unit cell of a hexagonal 2-net, depicting thelattice vectors and the shift vector.suitable loading is applied symmetry changes following the scheme (5) → (4) → (1) that isshown in Figure 1. For the hexagonal case arithmetic symmetries are given by the matrices − − −
11 0 00 0 1 , , − − −
10 1 00 0 1 , (7) − − −
10 0 1 , , − − −
10 0 1 . (8)To these matrices one should add their counterparts which have -1 at the (3, 3) componentof the matrix, since monoatomic 2-nets admit the central inversion as a symmetry opera-tion. The eigenvalues of these matrices are 1 , − , e i π , e − i π . Thus, these matrices describe:the identity transformation, central inversion and rotations by 60 , −
60 degrees. The planegroup corresponding to this case is p6mm ([9]).6or the rhombic case the matrices of the arithmetic symmetry read , . (9)The eigenvalues of these matrices are -1 and 1, so they describe central inversion andidentity transformation. As above, to these matrices one should add their counterpartswhich have -1 at the (3, 3) component of the matrix. The plane group for this case isc2mm ([9]).For the oblique case, the arithmetic symmetry is described by the matrix , (10)together with its counterpart with -1 in the (3, 3) component. The corresponding planegroup is p2.Diatomic nets are another class of a multilattice: they consist of 2 simple latticesbut now the lattice points are occupied by atoms belonging to different species (examplesinclude MoS , NbSe and WSe and hexagonal BN). In Figure 2 a diatomic 2-net is depictedwith solid and hollow circles referring to different atom species. For the planar case, theunit cell of diatomic 2-nets and their symmetry hierarchies are seen in Figure 3 ([11]) Thematerials we study here (hexagonal BN, MoS , WSe and NbSe ) are of the hexagonaltype (class (10) of Figure 3). Symmetry breaks following the scheme (10) → (8) → (2) ofFigure 3. Diatomic 2-nets differ from their 1-net counterparts since they lack the centralinversion. Thus, matrices with -1 at the (3, 3) component of the arithmetic matrices and-1 eigenvalues should be excluded, since they describe central inversion. For expressing the discrete nature of a lattice to the continuum scale, we use the Cauchy-Born rule (see e.g. [8] and references therein). Limitations of this rule can be found in [13](see also [3]), but in the present work validity of this rule is enforced.According to this rule, for multilattices, atomic motion of the lattice agrees with thegross deformation, while the shift vector is free to adjust so as to reach equilibrium. Weassume the existence of a stored energy function φ for the multilattice that depends on thecurrent lattice vectors and the shift vector φ = φ ( e a , p ) . (11)The Cauchy-Born rule then states that the reference, e a , and the current lattice vectors, e a , are related according to the formula ([25]) e a = Fe a , (12)7igure 3: Bravais lattices for diatomic 2-nets and their symmetry hierarchies (Figure takenfrom [11]). Solid circles denote atoms of one species, hollow circles refer to atoms belongingto another species.where F = ∂ x ∂ X is the well known deformation gradient of continuum mechanics. The shiftvector is adjusted so that equilibrium is reached. Using minimization arguments one mayshow that ([25]) ∂φ∂ p = 0 . (13)This equation plays the role of the field equation for the evaluation of the shift vector p .Applying the Cauchy-Born rule to the stored energy we write φ = φ ( F , p ) . (14)8lassical invariance of the energy then would require φ ( F , p ) = φ ( QFH T , pH T ) , H ∈ L ( e a ) , Q ∈ Q , F ∈ N e a , (15)but when confined to a weak transformation neighborhood this requirement becomes φ ( F , p ) = φ ( QFH T , pH T ) , H ∈ P ( e a ) , Q ∈ Q , F ∈ N . (16)The action of H is due to material symmetry, while the action of Q is due to frame indiffer-ence. Reduction of the arithmetic symmetry group to the geometric can be accomplishedwhen one is confined to weak transformation neighborhoods. In short terms, for simple lat-tices say that B is the 9 dimensional space of all basis e a . Then one can prove ([6, 21, 22])that there exists a neighborhood N ⊃ B such that the action of L ( e a ) coincides withthat of P ( e a ). Namely, when one is confined to such neighborhoods there is no need todistinguish between geometric and arithmetic symmetry groups. Similar arguments holdfor multilattices as well ([25]).So, confined to this neighborhood, material symmetry uses the geometric symmetrygroup, i.e. the crystal systems continuum mechanics uses. Frame indifference then leadsto φ = φ ( C , p ) , (17)where C = F T F is the right Cauchy-Green deformation tensor. So, the Cauchy-Born ruleallows the transition from a lattice to its continuum counterpart (see Section 6.2, [25]).What is new for multilattices is the dependence on the shift vector as well ([25, 3]). Themotivation for the transformation rule for the shift vector under the action of the symmetrygroup in eqs. (15, 16) is given in [26].To sum up, there are two crucial assumptions that are necessary for the validity of ourmodel: first is the enforcement of the Cauchy-Born rule and second the confinement of theanalysis to weak transformation neighbourhoods. This enables us to work with an energythat has the form φ = φ ( C , p ) , (18)augmented properly to take into account bending effects, while symmetry is the one con-tinuum mechanics uses. Arguments of Section 3 pertain to classical three dimensional (bulk) elasticity. In thisSection we lay down the kinematics of surface elasticity ([31, 33]). We assume in thereference configuration a smooth surface A , described by Y = Y ( θ , θ ), where Y is theposition vector for the point of the surface from the origin and the parameters θ α arecurvilinear coordinates on the surface. Covariant and contravariant vectors are defined by A α = Y ,α , A α · A β = δ αβ , (19)9here δ αβ is the kronecker delta in two dimensions. Deformation of the surface bringspoint Y to point y ( θ , θ ) on surface A, in the current configuration. The covariant andcontravariant vectors related with the surface A are then defined a α = y ,α , a α · a β = δ αβ , (20)The linear mapping that maps vectors on the tangent plane of A to those of A, is thesurface deformation gradient defined by F s = a α ⊗ A α . (21)The surface right Cauchy-Green deformation tensor is then defined as C s = F Ts · F s . (22)The surface curvature tensors in the reference and the current configuration can be ex-pressed as b = b αβ A α ⊗ A β , (23) b = b αβ a α ⊗ a β . (24)Surface divergence is defined by ∇ s () = ∇ () − n ( n · ∇ ()) , (25)where ∇ is the common divergence operator in three dimensions, while n and N are theoutward unit normals on surface A and A , respectively.A curvature dependent surface energy is an energy of the form ([31, 2, 16]) W = W ( C s , b ) . (26)So, aside from the in-surface strain measure C s , dependence on b is assumed. Thisdependence allows the modeling of bending effects since it takes into account out-of-planedeformations. For a monolayer 2D crystal this energy should be augmented with thedependence on the shift vector which gives W = W ( C s , b , p ) , (27)where p is the shift vector. The need for the dependence on the shift vector stems fromthe fact that monoatomic and diatomic 2-nets are multilattices. Material symmetry for curvature dependent surface energies is a subject tackled elegantlyby Steigmann and Ogden ([31] Section 6) based on the earlier work of Murdoch and Cohen10[16]). For such energies, elements of the symmetry group are pairs which leave the responseof the surface invariant to superposed deformations. Steigmann and Ogden ([31]) concludedthat for surfaces with constant non-negative curvature (namely, for planes and spheres)amenability to available representation theory is possible when the symmetry group reads { R , } while for the energy it then holds W ( C s , b ) = W ( RC s R T , Rb R T ) , RR T = I , det R = 1 . (28)The form of the energy should be augmented by taking into account dependence on theshift vector. So, collectively, the action of the symmetry group for a curvature dependentsurface energy reads W ( C s , b , p ) = W ( RC s R T , Rb R T , pR T ) . (29)The motivation for the transformation rule of the shift vector is given in [26].Field equations for materials described by curvature dependent surface energies arestudied in [1]. According to these authors, the momentum equation for the static case andin the absence of body forces reads σ bulk · n + ∇ S T S = 0 , (30)where T S is the surface first Piola-Kirchhoff stress tensor defined by ([31]) T S = ∂ ¯ W∂ F S (31)when W = ¯ W ( F S , b , p ), while σ bulk is the Cauchy stress tensor for the bulk materialsurrounding the surface. Here since we speak about free standing surfaces, there is no bulkmaterial, so the bulk stress tensor should be set equal to zero, σ bulk = . The surface firstPiola-Kirchhoff stress tensor is related to the second Piola-Kirchhoff surface stress tensor, S S , according to the formula S S = F − S · T S . (32)The second Piola-Kirchhoff stress tensor can also be written as S S = ∂W∂ C S . (33)The momentum equation in the absence of body forces and inertia can be expressed usingthe second Piola-Kirchhoff stress tensor as¯ ∇ S S S = , (34)where ¯ ∇ S () = ∇ C S − N ( N · ∇ C ()) is the surface gradient when the bulk gradient ∇ C istaken with respect to the metric C of the bulk. In this respect, the bulk divergence in eq.(25) is taken with respect to the metric G of the bulk reference configuration. Essentially,11q. (34) is the surface analog of the momentum equation written using the second Piola-Kirchhoff stress tensor in classical elasticity ([15, 29, 30]). From the physical point of view,the momentum equation is the force balance for the surface.The moment of momentum balance, in the absence of body couples and inertia reads([1]) when setting σ bulk = ∇ S m S − ∇ S ( F − S · S S × y ) = , (35)where the surface couple stress tensor is defined as ([31]) m S = ∂W∂ C S . (36)The symbol × in eq. (35) denotes the cross product of the three dimensional space. Themoment of momentum equation renders the couple balance for the surface.For the shift vector the field equation reads ([25, 3]) ∂W∂ p = . (37)Form the physical point of view, this equation says that the shift vector adjusts in suchaway that equilibrium is reached. To obtain the appropriate geometric symmetry group we first evaluate eigenvalues for thematrices of the arithmetic symmetry group. We then correspond these eigenvalues toappropriate generators of the geometric symmetry group. This is done for the hierarchies(5) → (4) → (1) for monoatomic 2-nets and (10) → (8) → (1) for diatomic 2-nets. For hexagonal lattices the arithmetic symmetry group is given by the matrices of eq. (7, 8)for monoatomic 2-nets. Evaluation of the eigenvalues gives 1 , − , e i π , e − i π which describeidentity transformation, central inversion and rotation by -60, 60 degrees. The space groupis p6mm. For their diatomic counterparts central inversion is excluded. The correspondingcrystallographic point group has generators R ( π ) and R j and belongs to class 10 accordingto the classification of Zheng ([34, 35]) and is denoted by C ν . R ( θ ) denotes a rotation ofangle θ and R j is a two-dimensional reflection transformation. Diatomic 2-nets have thesame geometric symmetry group, C ν , since this group does not admit central inversion.The structure tensor for C ν is denoted by P and defined by ([34]) P = Re ( a + i a ) , (38)12r equivalently as P = M ⊗ M ⊗ M − ( N ⊗ M ⊗ N + N ⊗ N ⊗ M ) , (39)where M = a ⊗ a − a ⊗ a , N = a ⊗ a − a ⊗ a , a , a , { a , a } an orthonormal basisvector. It can also be written as P = Re [ e i θ ( c + i c ) ] , (40)where c = cos ( θ ) a + sin ( θ ) a , c = − sin ( θ ) a + cos ( θ ) a . This tensor is an irreducibletensor since C ν is compact. In two dimensions it has only two independent components([35]) P = cos(6 θ ) , P = sin(6 θ ) . (41)These two components introduce the anisotropy and they model the zig-zag and armchairdirection. Since θ = π we have P = cos (2 π ) = 1 , P = sin (2 π ) = 0.Thus, for this case we have using curvature dependent surface energy W = W anisotropic ( C S , b , p ) , (42)where the anisotropy stems from the fact that the symmetry group is not the full orthogonalgroup. Reduction to an isotropic function is done through the use of the principle ofisotropy of space ([35]), that gives W = W anisotropic ( C S , b , p ) = W isotropic ( C S , b , P , p ) . (43)Thus, we take an isotropic function at the expense of using the structure tensor (thatdescribes the anisotropy) as an addittional argument.The complete and irreducible representation of such a scalar function under the group C ν consists of the following quantities ([35, 34]) I = tr C S , I = det C S , I = tr(Π C S C S ) , I = tr( C S b ) ,I = tr(Π b b ) , I = tr b , I = det b , I = p · C S p ,I = p · b p , I = p · π p , I = p · p , I = tr(Π p C S ) , I = tr(Π p b ) . (44)Thus, in general for such a model we have the following expression for the energy W = ˜ W ( I i ) , i = 1 , , ..., . (45)The term Π A for a symmetric tensor of second order A is defined as ([35], indices rangingfrom 1 to 2) Π A = P ijklmn A kl A mn c i ⊗ c j , (46)and renders a second order tensor. The term π x with respect to the vector z is defined as π z = P ijklmn z j z k z l z m z n c i , (47)13hile for Π z we have Π z = P ijklmn z k z l z m z n c i ⊗ c j . (48)The material parameters related to I , I describe pure bending effects since det b , tr b are the mean and the Gaussian curvature of the surface, respectively. The term relatedwith I describes the effect of the armchair and the zig-zag direction at bending. Theparameters related with I , I are related to pure stretching, while those related with theterm I describe the anisotropy (zig-zag, armchair directions) to stretching. The parameterrelated to I describes coupling between bending and stretching response. The terms I , I describe the effect of the in plane and the out of plane deformations, respectively, on theshift vector. The term I describes the way anisotropy affects the shift vector, while I describe changes related with the shift vector solely. Terms I , I are related to couplingof anisotropy with the shift vector for the in plane and the out of plane deformations,respectively.The surface stress tensor and the surface couple stress tensor are evaluated by deter-mining the derivatives S S = ∂ ˜ W∂ C S , m S = ∂ ˜ W∂ b . (49)Also, the field equation for the shift vector is ∂ ˜ W∂ p as shown in eq. (37). Using the expressionsof eq. (44) in eq. (49), after some calculations we obtain S S = ∂ ˜ W∂I I + ∂ ˜ W∂I [tr( C S ) − C S ] + 3 ∂ ˜ W∂I P : ( C S ⊗ C S ) + ∂ ˜ W∂I b + ∂ ˜ W∂I p ⊗ p + ∂ ˜ W∂I Π p , (50) m S = ∂ ˜ W∂I I + ∂ ˜ W∂I [tr( b ) − b ] + 3 ∂ ˜ W∂I P : ( b ⊗ b ) + ∂ ˜ W∂I C S + ∂ ˜ W∂I p ⊗ p + ∂ ˜ W∂I Π p , (51) ∂W∂ p = 2 ∂ ˜ W∂I C S p + 2 ∂ ˜ W∂I b p + 6 ∂ ˜ W∂I P • ( p ⊗ p ⊗ p ⊗ p ⊗ p ) + ∂ ˜ W∂I p +4 ∂ ˜ W∂I [ P : ( p ⊗ p ⊗ p )] : C S + 4 ∂ ˜ W∂I [ P : ( p ⊗ p ⊗ p )] : b . (52)By making the simplest possible assumption that ˜ W is linear with respect to the invariants I i , i = 1 , , , ...,
13 we take S S = α I + β [tr( C S ) − C S ] + 3 γ P : ( C S ⊗ C S ) + δ b + θ p ⊗ p + ρ Π p , (53) m S = ǫ I + ζ [tr( b ) − b ] + 3 η P : ( b ⊗ b ) + δ C S + ι p ⊗ p + τ Π p , (54) ∂W∂ p = θ C S p + ι b p + 6 λ P • ( p ⊗ p ⊗ p ⊗ p ⊗ p ) + ξ p +4 ρ [ P : ( p ⊗ p ⊗ p )] : C S + 4 τ [ P : ( p ⊗ p ⊗ p )] : b . (55)14he Greek letters α, β, γ, δ, θ, ρ, ǫ, ζ , η, ι, τ, λ, ξ are material parameters to be determinedby experiments.Written with respect to indices ranging from 1 to 2 these formulas are S S AB = αδ AB + β [ tr ( C S ) δ AB − C S AB ] + 3 γP ABCDEF C S EF C S CD + δd AB + θp A p B + ρP ABCDEF p C p D p E p F , (56) m S AB = ǫδ AB + ζ [ tr ( b ) δ AB − b AB ] + 3 ηP ABCDEF b EF b CD + δC S AB + ιp A p B + τ P ABCDEF p C p D p E p F , (57) ∂W∂p A = θC S AB p A + ιb AB p B + 6 λP ABCDEF p B p C p D p E p F + ξp A +4 ρP ABCDEF p D p E p F C S BC + 4 τ P ABCDEF p D p E p F b BC . (58)The elasticities of this model are given by the following fourth order tensors A = ∂ W∂ C S , B = ∂ W∂ b , C = ∂ W∂ C S ∂ b . (59)Quantities of the first term are related to the in-plane motion, the second to the out-of-plane while the third to the coupling between in-plane and out-of-plane motions. Allin all, modeling of the hexagonal lattice at the continuum level for monoatomic and di-atomic 2-nets requires specification of 13 material parameters; these are the Greek letters α, β, γ, δ, θ, ρ, ǫ, ζ , η, ι, τ, λ, ξ . Breaking of symmetry is dictated by Figures 1 and 3 which describe the symmetry hierar-chies for monoatomic and diatomic 2-nets. We again follow the approach of evaluating theeigenvalues of the matrices of the arithemtic symmetry group and finding the correspond-ing geometric symmetry groups with generators describing the same symmetry operation.For the first braking of symmetry one has to distinguish between monoatomic and diatomic2-nets.
For monoatomic 2-nets the rhombic unit cell has arithmetic symmetry group describedby eq. (9), and the corresponding eigenvalues are 1 and -1. The corresponding geometricsymmetry group is the group C ν (no. 4 in the classification adopted in [34, 35]) with gen-erators R ( π ) , R j . This case corresponds to orthotropy in two dimensions and the structuretensor in this case is the tensor M = a ⊗ a − a ⊗ a .Thus, in this case we have W anisotropic ( C S , b , p ) = W isotropic ( C S , b , p , M ) . (60)15he complete and irreducible representation of such a scalar function under the group C ν consists of the following quantities ([34, 35]) I = tr C S , I = det C S , I = tr( MC S ) , I = tr( C S b ) ,I = tr( Mb ) , I = tr b , I = det b , I = p · p ,I = p · Mp , I = p · C S p , I = p · b p . (61)Thus, for the energy we obtain W = ˜ W ( I i ) , i = 1 , , ..., . (62)The surface stress and surface couple stress that correspond to this energy are calculatedas S S = ∂ ˜ W∂I I + ∂ ˜ W∂I [tr( C S ) − C S ] + ∂ ˜ W∂I M : ( C S ⊗ C S ) + ∂ ˜ W∂I b + ∂ ˜ W∂I p ⊗ p . (63) m S = ∂ ˜ W∂I I + ∂ ˜ W∂I [tr( b ) − b ] + ∂ ˜ W∂I M : ( b ⊗ b ) + ∂ ˜ W∂I C S + ∂ ˜ W∂I p ⊗ p . (64)For the term ∂W∂ p we have ∂W∂ p = ∂ ˜ W∂I p + ∂ ˜ W∂I Mp + ∂ ˜ W∂I C S p + ∂ ˜ W∂I b p . (65)By making the simplest possible assumption that ˜ W is linear with respect to the invariants I i , i = 1 , , , ...,
11 we take S S = α G + β [tr( C S ) − C S ] + γ M : ( C S ⊗ C S ) + δ b + θ p ⊗ p , (66) m S = ǫ G + ζ [tr( b ) − H S ] + η M : ( b ⊗ b ) + δ C S + ι p ⊗ p (67)and ∂W∂ p = λ p + µ Mp + θ C S p + ι b p . (68)Now the effect of anisotropy at the level of the constitutive law is introduced throughtthe terms where the structure tensor, M , is present. Namely the material parameters γ, η, µ measure the effect of anisotropy at the in-plane motion, the out-of-plane motionand the shift vector, respectively. For monoatomic 2-nets the number of independentmaterial constants to be observed and measured in experiments are 11: the Greek letters α, β, γ, δ, θ, ǫ, ζ , η, ι, λ, µ . Diatomic 2-nets lack central inversion. Thus, we assume that the corresponding geomet-ric symmetry group is the group C ν (no.3 in the classification adopted in [34, 35]) with16enerators R j . In this case the structure tensor is the vector a . So, for this case energyreads W anisotropic ( C S , b , p ) = W isotropic ( C S , b , p , a ) . (69)The complete and irreducible representation of such a scalar function under the group C ν consists of the following quantities ([34]) I = tr C S , I = tr( MC S ) , I = tr C S , I = p · a , I = p · p I = tr( C S b ) , I = p · C S a , I = tr b , I = tr( Mb ) ,I = tr b , I = p · b a . (70)Thus, for the energy we obtain W = ˜ W ( I i ) , i = 1 , , ..., . (71)The surface stress tensor for this case then reads S = ∂W∂ C S = ∂ ˜ W∂I I + ∂ ˜ W∂I M + ∂ ˜ W∂I [(tr C S ) I − C S ] + ∂ ˜ W∂I b + ∂ ˜ W∂I p ⊗ a = α I + β M + γ [(tr C S ) I − C S ] + δ b + ǫ p ⊗ a , (72)and the energy is assumed to be linear with respect to the invariants. For the surfacecouple stress under the same energy assumption we then have m = ∂W∂ b = ∂ ˜ W∂I C S + ∂ ˜ W∂I I + ∂ ˜ W∂I M + ∂ ˜ W∂I [(tr b ) I − b ] + ∂ ˜ W∂I p ⊗ a = δ C S + ζ I + η M + θ [(tr b ) I − b ] + ι p ⊗ a . (73)For the term ∂W∂ p we find ∂W∂ p = ∂ ˜ W∂I a + ∂ ˜ W∂I p + ∂ ˜ W∂I ( C S · a ) + ∂ ˜ W∂I ( b · a )= κ a + λ p + ǫ ( C S · a ) + ι ( b · a ) (74)The material parameters related with I , I describe pure stretching, while those relatedwith I describe the effect of anisotropy to pure stretching. Pure bending characteristicsare introduced throught terms I , I , while the effect of anisotropy to pure bending isdescribed by I . Combined stretching and bending effects are given by I , while theeffect of anisotropy to the shift vector is throught the term I . Terms I , I describe howanisotropy affect the combined response of the shift vector with stretching and bending,respectively. Finally, I is a term for the shift vector solely. For diatomic 2-nets the numberof independent material constants to be observed and measured in experiments are 11: theGreek letters α, β, γ, δ, θ, ǫ, ζ , η, ι, κ, λ . 17 .3 Second breaking of symmetry In this case both monoatomic and diatomic 2-nets have as geometric symmetry group thegroup C (no. 1 in the classifiaction adopted by [34, 35]). The structure tensors for thisgroup are the vectors a , a . So, for the energy we have W anisotropic ( C S , b , p ) = W isotropic ( C S , b , p , a , a ) . (75)Invariants for this energy are then I = tr C S , I = tr( MC S ) , I = tr( NC S ) , I = p · a , I = p · a I = tr b , I = tr( Mb ) , I = tr( Nb ) . (76)Thus, for the energy we obtain W = ˜ W ( I i ) , i = 1 , , ..., . (77)By considering that the energy is linear with respect to the invariants we find thesurface stresses S S = ∂W∂ C S = ∂ ˜ W∂I I + ∂ ˜ W∂I M + ∂ ˜ W∂I N = α I + β M + γ N . (78)For the surface couple stress tensor we have m S = ∂W∂ b = ∂ ˜ W∂I I + ∂ ˜ W∂I M + ∂ ˜ W∂I N = δ I + ǫ M + ζ N . (79)For the term related with the shift vector we evaluate ∂W∂ p = ∂ ˜ W∂I a + ∂ ˜ W∂I a = η a + θ a (80)Terms I , I describe pure stretching and bending, respectively, effects. The effect ofanisotropy to stretching and bending is described by terms I , I and I , I , respectively.The effect of anisotropy to the shift vector is introduced throught terms I , I . All inall, the second breaking of symmetry requires specification of the 8 material parameters: α, β, γ, δ, ǫ, ζ , η, θ . We present a nonlinear elastic constitutive framework for the modeling of 2D crystals ofcurrent interest such as graphene, hexagonal BN, MoS , WSe , and NbSe . We use the18heory of monoatomic and diatomic 2-nets to find their arithmetic symmetries. Confinedto weak transformation neighbourhoods and using the Cauchy-Born rule we are able towork with geometric symmetries. For finding the geometric symmetry group we evaluatethe eigenvalues of the arithmetic symmetries and find the corresponding generators amongthe crystallographic point groups. We then apply the theory of invariants for an energyfunction depending on the surface Cauchy-Green deformation tensor, the curvature tensorand the shift vector. We lay down the expression for the stress tensor, the couple stresstensor as well as a term related with the shift vector. This is done for all case wheresymmetry changes due to applied loads.Future directions of our line of research are at two levels. Firstly, one important aspect isto find the symmetry breaking and symmetry preserving deformations for such materialsin line with the approach of [25]. Secondly, we currently work on utilizing non-convexenergies that can capture the phenomena at the transition regime, i.e., when symmetrybreaks, in line with fundamental works in zirconia ([32]); this might be expanded to thenonlinear case. D. Sfyris sincerely thanks G. Zanzotto (Padova, Italy) for some very important clarifica-tions concerning differences between monoatomic and diatomic 2-nets. This research hasbeen co-financed by the European Union (European Social Fund - ESF) and Greek na-tional funds through the Operational Program ”Education and Lifelong Learning” of theNational Strategic Reference Framework (NSRF) - Research Funding Program: ERC-10”Deformation, Yield and Failure of Graphene and Graphene-based Nanocomposites”. Thefinancial support of the European Research Council through the projects ERC AdG 2013(Tailor Graphene) is greatfully acknowledged.
References [1] P. Chhapadia, P. Mohammadi, P. Sharma, Curvature dependent surface energy andimplications for nanostructures. J. Mech. Phys. Sol. 59 (2011) 2103-2115.[2] H. Cohen, C. N. DeSilva, Nonlinear theory of elastic surfaces. J. Math. Phys. 7 (1966)246-253.[3] W. E, P. Ming, Cauchy-Born rule and the stability of crystalline solids: static problem.Arch. Rat. Mech. Anal. 183 (2007) 241-297.[4] J. L. Ericksen, Nonlinear elasticity of diatomic crystals. Int. J. Sol. Struct. 6 (1970)951-957.[5] J. L. Ericksen, On the symmetry of deformable crystals. Arch. Rat. Mech. Anal. (1979)72, 1-13. 196] J.L. Ericksen, Some phase transitions in crystals, Arch. Rat. Mech. Anal. 73 (1980)99-124.[7] J. L. Ericksen, Multi-valued strain energy functions for crystals. Int. J. Sol. Struct. 18(1982) 913-916.[8] J.L. Ericksen, On the Cauchy-Born rule, Math. Mech. Sol. 13 (2008) 199-220.[9] G. Fadda, G. Zanzotto, The arithmetic symmetry of monoatomic 2-nets. Acta Crystal.A, 56 (2000) 36-48.[10] G. Fadda, G. Zanzotto, Symmetry breaking in monoatomic 2-lattices. Int. J. Non-Linear Mechanics 36 (2001) 527-547.[11] G. Fadda, G. Zanzotto, The arithmetic symmetry of colored crystals: classification of2-color 2-lattices. J. Appl. Cryst. 37 (2004) 1-7.[12] G. Fadda, G. Zanzotto, On the arithmetic classification of crystal structures. ActaCryst. A, 57 (2001) 492-506.[13] G. Friesecke, F. Theil, Validity and failure of the Cauchy-Born hypothesis in a twodimensional mass spring lattice. J. Nonl. Sci. 12 (2002) 445-478.[14] M.E. Gurtin, A. I. Murdoch, A continuum theory of elastic material surfaces, Arch.Ration. Mech. Anal. 57 (1975) 291-323.[15] J.E. Marsden, T.J.R. Hughes, Mathematical foundations of elasticity. Dover Publica-tions (1994).[16] A.I. Murdoch, H. Cohen, Symmetry consideration for material surfaces. Arch. Rat.Mech. Anal. 72 (1979) 61-97.[17] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotkevick, S.V. Morozov,A.K. Geim, Two dimensional atomic crystals. PNAS 102 (2005) 10451-10453.[18] G. Parry, On diatomic crystals. Int. J. Sol. Struct. 14 (1978) 281-187.[19] G. P. Parry, On phase transitions involving internal strain. Int. J. Sol. Struct. 17(1981) 361-378.[20] G. P. Parry, On internal variable models of phase transitions. J. Elast. 17 (1987) 63-70.[21] M. Pitteri, Reconciliation of local and global symmetries of crystals, J. Elast. 14 (1984)175-190.[22] M. Pitteri, On νν