## Abstract

We establish the most general form of the discrete elasticity of a two-dimensional triangular lattice embedded in three dimensions, taking into account up to next-nearest-neighbour interactions. Besides crystalline system, this is relevant to biological physics (e.g. red blood cell cytoskeleton) and soft matter (e.g. percolating gels, etc.). In order to correctly impose the rotational invariance of the bulk terms, it turns out to be necessary to take into account explicitly the elasticity associated with the vertices located at the edges of the lattice. We find that some terms that were suspected in the literature to violate rotational symmetry are, in fact, admissible.

## 1. Introduction

Since Born & von Kármán, the construction of microscopic elastic models of crystals has formed an important part of solid-state physics, aimed to determine elastic constants and lattice dynamical properties (Born & von Kármán 1912; Born 1914). Imposing global translational and rotational invariance is a delicate matter, although the required relationships are well established (Born & Huang 1954). In fact, inconsistencies associated with rotational invariance in the earlier models were pointed out by Lax (1965) and Keating (1966*a*,*b*). More recently, the importance of triangular lattices has emerged with the study of soft matter systems, such as percolating gels (Feng & Sen 1984; Schwartz *et al*. 1985; Arbabi & Sahimi 1993), membrane with crystalline order (Kantor & Nelson 1987; Seung & Nelson 1988; Gompper & Kroll 1997), biological membranes (Discher *et al*. 1997; Coughlin & Stamenovic 2003), red blood cell cytoskeleton (Saxton 1990; Lim *et al*. 2002; Marcelli *et al*. 2005) and finite element modelling (Gusev 2004).

In this paper, we consider a triangular lattice embedded in three dimensions and we determine the most general expression of the harmonic elastic energy up to next-nearest-neighbour interactions. We formulate a systematic procedure to generate all the energy terms allowed by symmetry. This procedure could be applied to other types of lattices. We show that in order to correctly ensure rotational invariance, it is necessary to consider the boundary of the lattice, as first suggested by Lax (1965). Our central result, equation (3.2), gives *all* the allowed energy terms, with no particular connection to mechanical elements, such as springs, elastic wedges, etc. This contrasts with the above-cited papers that postulate energy terms in a non-exhaustive way. In particular, the most general energy contains a non-central nearest-neighbour term (with coefficient ), which was believed to violate rotational invariance in percolating gels (Feng & Sen 1984) on the grounds of Keating's argumentation (1966*a*,*b*). The situation is actually delicate since in a later paper, Keating (1968) restricted the generality of some assertions of its earlier papers.

Our paper is organized as follows. In §2, we formulate the elastic energy of triangular lattices. In §2*a*, we write down all the scalar energy terms quadratic in the displacements of the lattice vertices. In §2*b*, we retain only the terms satisfying the sixfold symmetry of an infinite triangular lattice. In §2*c*, we also consider the edges: we study the symmetries relative to them and we retain the allowed edge terms. In §2*d*, we impose global translational and rotational invariance, taking explicitly into account the reduced symmetry of the boundary, and we retain the allowed terms. In §2*e*, we collect all the contributions and we extract the resulting bulk terms. In §3, we analyse and discuss the results obtained in §2. In §3*a*, we give a possible interpretation of our energy in terms of mechanical elements. In §3*b*, we rewrite the energy into the generic expression (3.2), which constitutes our central result. In §3*c*, we establish the connection between our model and its continuum version. Finally, in §4, we give a brief summary and we detail the origin of the rotational invariance problem.

## 2. Discrete elasticity model

### (a) General definitions

Let us consider a bounded triangular lattice. In its undistorted state, we assume that it lies in a plane and that its vertices (i.e. its sites) are evenly separated by a distance (figure 1). We denote by , the position of a generic vertex *α*. We also assume that the elastic properties of the lattice are symmetric with respect to .

In an arbitrary distorted state, becomes , where is a three-dimensional displacement. The most general quadratic expansion of the lattice elastic energy, *F*, is given by(2.1)In this expression, we have separated the terms coupling identical vertices from those coupling different vertices. We have excluded the terms linear in , because they must vanish for the undistorted state to be a minimum of the energy. Since *F* is scalar, the tensors and must satisfy and , where the superscript ‘t’ means transpose. Note that, in general, .

We fix the *interaction range* by assuming that all the vanish except when *α* and *β* are nearest or next-nearest-neighbours.

### (b) Bulk terms

To begin with, we consider the *bulk* vertices, i.e. the vertices that are far away enough from the boundary (this point will be found precisely in §2*c*), so that the triangular symmetry may be assumed. In the bulk, the elastic properties are assumed to be spatially homogeneous. Let be a unitary vector normal to . Let be two three-dimensional vectors forming a *local* orthonormal basis in , allowed to vary from site to site. The tensors and can be decomposed into the complete basis , where the symbol represents the tensorial product. Elements linear in , such as , have been excluded, since the undistorted lattice is symmetric with respect to .

To construct , we take along any one of the six nearest-neighbour directions of the undistorted lattice and . Because of the triangular symmetry, each term in the decomposition of must be even in and in , and the coefficients of and must be equal (this condition is necessary for to be invariant under an in-plane rotation). Hence,(2.2)where is the identity tensor in the plane . By virtue of the lattice homogeneity, the coefficients and should not depend on *α*.

We now construct . Let be the unit vector in the direction going from to . We take and . Without loss of generality, the decomposition of over can be rewritten as(2.3)By virtue of the lattice homogeneity, the coefficients take only two sets of values: either if *α* and *β* are nearest neighbours, or if *α* and *β* are next-nearest-neighbours. The in-plane mirror symmetry with respect to requires that must be even in , which leads to(2.4)

#### (i) Collecting bulk terms

We now anticipate that we shall clearly distinguish between the and belonging to the bulk and those belonging to the edges. Denoting by the part of the energy *F* which involves only the *bulk* terms, we arrive, from equations (2.1)–(2.4), at(2.5)where denotes the normal component of and , the projection of on . The symbol indicates the summation over both nearest and next-nearest-neighbour pairs of vertices. As already mentioned, we assign to the former the suffix *i*=1 and to the latter the suffix *i*=2 (see figure 2 where these numbers are indicated as I and II). The function takes the value 1 if the pair {*α*, *β*} is of type *i* and zero otherwise. We shall hereafter call such functions *indicatrix functions*.

### (c) Edge terms

With arbitrary coefficients and , the energy form (2.5) is not invariant under a global translation or rotation of the lattice. In order to correctly ensure rotational invariance, it is necessary to take into account the *edges* of the lattice. This seems counterintuitive, because one expects the contribution of the edges to be subdominant with respect to the expensive contribution of the bulk. However, as it will more precisely be shown in §4, this is not the case, because in a rotation the displacements of the edges are proportional to the size of the system.

#### (i) Neighbour clusters and vertices types

Since our interest lies in the bulk terms, we choose the simple shape of a regular hexagon of edge length (figure 1) as our system. Consistently with the energy expansion, the specificity of the edges must be taken into account up to the next-nearest-neighbour distance. For this purpose, we group the vertices into *neighbour clusters*, in which all the inter-distances are at most equal to the next-nearest-neighbour distance . There are seven types of such clusters, which we label from I to VII (see figure 2). A vertex is said to participate in a neighbour cluster of type Y if it can be grouped with its neighbour vertices into a neighbour cluster of type Y. Note that vertices near the edges participate in a smaller number of neighbour clusters than those in the bulk, because of their lacking neighbours.

To each vertex *α*, we now assign a *type* depending on which neighbour clusters it participates in: two vertices are assigned the same type if and only if their positions within the set of neighbour cluster in which they participate are indistinguishable. Close inspection shows that the edge vertices are classified into five types, labelled from 1 to 5, as shown in figure 3. The vertices of type 1 are the bulk vertices.

We apply a similar procedure to assign a type to the oriented nearest and next-nearest-neighbour vertex pairs, which we call *vertex doublets*. We distinguish the vertex doublets from , since generally . Two vertex doublets are assigned the same type if and only if their positions and orientations within the set of neighbour clusters in which they participate are indistinguishable. We thus determine eight types of vertex doublets, labelled from 1 to 8 (figure 3). The vertex doublets of type 1 are the bulk nearest-neighbour vertex doublets and those of type 2 are the bulk next-nearest-neighbour vertex doublets. We also define types for the vertex *pairs* in the following way. When the two doublets associated with a given pair are of the same type, the pair is assigned that type (this holds for the types 1, 2, 5, 6), while the pairs corresponding to doublet types 3 and 7 are assigned the type 3 and those corresponding to doublet types 4 and 8 are assigned the type 4.

#### (ii) Symmetry allowed edge terms

In order to take into account the local symmetry up to the next-nearest-neighbour distance, we consider the *symmetries of the set of neighbour clusters*, in which the vertices or vertex doublets participate. Let us first examine the coefficients of the tensors for vertices of types 2–5 (edge vertices). Since each one possesses an axis of symmetry, we choose along this axis and as indicated in figure 3. Then, as must be even in , we arrive at(2.6)where may either be or . For definiteness, we choose for vertices of type 2 or 4 and for vertices of type 3 or 5.

We now construct the tensor for the vertex doublets of types 3–8 (edge types). The most general expression for is given by an expression similar to (2.3):(2.7)We choose and as indicated in figure 3. The vertex doublets of types 4 and 8 have a mirror symmetry with respect to , hence must be even in , which implies . Other relationships can be obtained by comparing the coefficients of and those of . First, note that implies , , , and . Then, since the vertex doublets of types 5 and 6 do not change type when exchanging *α* and *β*, their coefficients must satisfy , , , and . Therefore, for types 5 and 6. Finally, we find(2.8)where , and . In the following, we define as the tensor that performs a counterclockwise rotation in the plane ; depending on the relative orientations of and , is either equal to or to . Note that, in general, we do find , contrary to what is the case in the bulk.

#### (iii) Collecting edge terms

The total elastic energy can now be written as , with given by (2.5) and obtained by collecting the edge terms defined above. We regroup the terms and in order to sum on the pairs and not on the doublets. We arrive then at(2.9)In this expression, and are the indicatrix functions of the vertices of type *i* and of the vertex pairs of type *i*, respectively (these types are defined in §2*c*(i)), and the symbol indicates the summation over all pairs of vertices. The third line in (2.9) requires more explanations. We define a closed contour by grouping the vertex pairs of type 3 (see figure 3), another closed contour by grouping the vertex pairs of type 5, and three separate closed contours by grouping the vertex pairs of type 6. The symbol denotes the sums over all these contours followed *counterclockwise*. The symbol is equal to 1 or −1 according to whether the direction from *α* to *β* goes outwards or inwards, respectively. Finally, in order to simplify the notation, we have renamed , , , and .

### (d) Rotational and translational invariance

Our task is now to ensure that the elastic energy remains unchanged if the lattice is translated or rotated as a whole. Although the general relationships are well established (Born & Huang 1954), it is not an easy task to implement them in such a way as to put the elastic energy in a satisfactory form. Our strategy is to rearrange *F* in order to put forward a number of terms that are obviously invariant, then to apply deliberately chosen translations or rotations in order to show that the remaining terms should vanish.

We first separate in the contribution depending on the components from the contribution depending on the components . We thus have with(2.10)(2.11)

#### (i) Translation invariance

Putting forward obviously translationally invariant terms in the form of differences, can be rewritten as(2.12)where the are linear combinations of the coefficients and (the coefficients in (2.12) still form an independent set). In appendix A*a*, we show that the invariance of *F* under any translation along requires(2.13)

Similarly, can be written as obviously translationally invariant terms plus remainders:(2.14)Whereas the first line was clearly obtained by putting forward squares and discarding the remainders in the last line, the other lines deserve detailed explanations. Here and afterwards, the symbol denotes the sum over all the triangular neighbour clusters of type Y, with *α*, *β* and *γ* taken in counterclockwise order. For type IV, the vertex *β* is supposed to be the central one, as in figure 2. The symbol is the indicatrix function of the neighbour clusters of type Y lying within the closed contour made by the vertex pairs of type *i* ( or 6). More precisely, for *i*=6, it means to lie within at least one of the three closed contours made by the vertex pairs of type 6. In the following, we shall often encounter the translationally and rotationally invariant expression (as in the second line):(2.15)where and or are the angles made by the vertices triplet , with *β* being the apex, after and before an elastic distortion of the lattice, respectively. In appendix A*b*(i), we show how the terms involving in (2.11) were rewritten with the help of (2.15), and we give the expressions of the coefficients and that were modified by this procedure. Finally, the lines involving were obtained as detailed in appendix A*b*(ii). The coefficients and are linear combinations of the coefficients , , , , and and the new coefficients in (2.14) still form an independent set. In appendix A*b*(iii), we show that the invariance of *F* under in-plane translations requires(2.16)

Consequently, the total energy can be written in the following form, obviously invariant under any translation of the system:(2.17)where is the first term and is the sum of all the other terms.

#### (ii) Rotation invariance

Putting forward rotationally invariant terms, , as obtained from (2.17), can be rewritten as(2.18)The symbol denotes the sum over all the neighbour clusters of type Y=VI or VII with the vertices placed as shown in figure 2. The angle denotes the wedge angle between the planes defined by the vertices and . In appendix A*c*(i), we show how the wedge angles are related to the normal displacements . The coefficients are linear combinations of , and , and the new coefficients still form an independent set. We show in appendix A*c*(ii) that the invariance of *F* with respect to out-of-plane rotations requires(2.19)Finally, we put forward in seven rotationally invariant terms, obtained by considering various energy terms associated with deviations of lengths or angles with respect to their equilibrium values:(2.20)The abbreviation ‘c.p.’ stands for ‘circular permutations’ of and is the distance between the vertices *α* and *β*. The detail of how (2.20) and (2.17) are related is given in appendix A*d*(i). The coefficient *G* is a linear combination of and the other 's and the coefficients are linear combinations of the corresponding and the 's . The new coefficients still form an independent set. In appendix A*d*(ii), we show that the invariance of *F* with respect to in-plane rotations implies(2.21)

### (e) Final expression of the bulk terms

From §2*d*, it follows that the most general elastic energy *F* of an hexagonal triangular lattice, invariant with respect to translation and rotations, is given by the sum of (2.18) and (2.20), with all the 's, 's and *G* equal to zero.

Putting forward the *bulk terms*, and expressing them in terms of the displacements instead of angles and lengths (using (2.15) and the formulae of appendix A*c*(i),*d*(i)), we finally obtain(2.22)where and and the ellipsis represents terms involving only edge vertex pairs.

## 3. Analysis

### (a) Comments on the energy form (2.22)

As shown in figure 4, each term of (2.22) can be given a simple interpretation in terms of linear or angular springs. The terms with coefficients , corresponding to elastic wedges, tend to keep the lattice flat and yield a form of curvature energy. The terms with coefficients *k* and represent elastic bonds acting between nearest and next-nearest-neighbours, respectively. The terms with moduli , and are angular springs acting at the bond joints.

The decomposition (2.22) is in many respects non-unique. For instance, the term with coefficient could be omitted if we are only interested in the bulk terms, for it has the same decomposition on the bulk vertex pairs as it has on the term with coefficient ; this can be easily checked by comparing the coefficients of and in the equations of appendix A*c*(i). Note, however, that the expressions of these two energies in terms of clusters of four vertices differ, as we can see in the first two terms of (2.22). As another example, we could have included in (2.22) energy terms associated with the angles in the bulk neighbour clusters of type IV; this would bring nothing new, since they can be written as a linear combination of other terms that are already present:(3.1)Nonetheless, although the *form* of (2.22) depends on the choices we made during the various transformations on the energy, our procedure guarantees that it does correspond to the most general elastic energy of a triangular lattice up to next-nearest-neighbour interactions.

### (b) Compact expression of the elastic energy

With the help of the formulae given in appendix A*c*(i),*d*(i), the energy (2.22) can be rewritten in a more compact form, in terms of the normal, transversal and longitudinal components of the relative displacements between nearest and next-nearest-neighbours:(3.2)where and the symbols and indicate the sum over all the neighbour clusters of types I and II, i.e. over the pairs of nearest and next-nearest vertices, respectively. The normal coefficient, which corresponds to curvature energy, is given by . The transverse coefficients and are given by and , and the longitudinal coefficients and are given by and . These new coefficients are independent linear combinations of the former. Moreover, the five functions of (proportional to the coefficients , , , and ), appearing in the form (3.2), can be shown to be linearly independent in their functional space. Therefore, the seven (or more) rigidities which one could associate with the mechanical elements drawn in figure 4 reduce to only five independent elastic constants.

A number of comments follows. The angular coefficients of (2.22), i.e. , and , are present both in the tranverse and the longitudinal coefficients; this is because in a triangular lattice one cannot change the angles between bonds without changing the distances between the vertices. Conversely, the coefficients *k* and appear only in the longitudinal coefficients. Also, as explained above, the two coefficients and collapse into a unique normal coefficient . Most importantly, note that the non-central term with coefficients , which was understood from Keating (1966*a*,*b*) to violate rotational invariance, is actually allowed (even if we restrain the model to nearest-neighbour interaction only). We shall comment on this in §4.

Finally, although (3.2) is expressed as a sum of squared terms, the stability issue is not trivial, since these terms are not independent. For instance, the stability of the undistorted lattice does not require all the five coefficients of (3.2) to be positive. For example, if all the coefficients in (2.22) are positive, which is obviously a sufficient condition for stability, and if , we may have . Further comments on the stability issue will be given in §4.

### (c) Link with the continuous theory

In a continuous description, the displacement of the point on the plane is described by a three-dimensional continuous function . In order to link the discrete description to the continuous one, we assume that the displacements vary very slowly at the scale of the lattice spacing and we identify and . Hence, we transform (3.2) with the help of the Taylor expansion:(3.3)where , with for nearest-neighbour vertices and for next-nearest-neighbour vertices, and . We take an orthonormal coordinate system , with the *z*-axis parallel to and the *x*-axis along one of the three nearest-neighbour directions of the undistorted lattice. Then, for the in-plane part of (3.2), we obtain , with(3.4)in which and and is the derivative of the *i*th component of with respect to the *j*th component of (note that , whereas ). Using integrations by parts, we can rewrite the energy density in the standard, isotropic, form , with . The Lamé coefficients are given by(3.5)(3.6)In terms of the spring constants of (2.22), this gives and . For a system composed only of springs between nearest neighbours, we recover the classical relation (Kantor & Nelson 1987). As expected, the area compressibility modulus depends only on the spring constants and not on the angular constants: .

For the normal part of the energy, it is necessary to add the second-order terms in the Taylor expansion because the first-order vanishes. With , we obtain the isotropic energy density,(3.7)with the bending constant . One recognizes the bending energy associated with the curvature of the surface .

## 4. Summary and discussion

In this paper, we have determined the most general form for the elastic energy of a triangular lattice when both first and second neighbour couplings are taken into account. We have found that, among the six coefficients originating from the decomposition of the energy in terms of the normal, transversal and longitudinal components of the relative displacements between nearest and next-nearest-neighbours, only five are independent (3.2). Consequently, only five combinations of the elastic parameters given, for instance, in figure 4 are measurable. We have also found that only three coefficients survive in the large-scale limit, the Lamé coefficients and the bending modulus; therefore, our discrete description introduces two new elastic coefficients. From the following well-known stability criteria: *κ*>0, *μ*>0 and , it is possible to deduce some necessary stability conditions: , and . A complete stability analysis can be performed by Fourier transform. However, the additional stability conditions on the coefficients being rather intricate, we shall not describe them here.

One of the main issue in the paper was the necessity to take into account the boundary elastic energy in order to correctly impose rotational invariance. This is well exemplified by the following energy term, which represents the sum of the energies of the angular springs between nearest bonds:(4.1)where *i* runs over all the pairs of adjacent bonds, making an angle , and runs over the interior bonds when (*B*) is indicated, or over the boundary bonds when is indicated. In the last sum, the bonds must be taken counterclockwise. Indeed, the left-hand side is manifestly rotationally invariant and, under an infinitesimal in-plane rotation , the last term of the right-hand side (r.h.s.), which is a boundary term, can be shown to yield an extensive contribution , where *A* is the area of the lattice, just as the first term of the r.h.s, which is a bulk term. It follows that: (i) this first term is not rotationally invariant by itself, (ii) boundary terms are therefore necessary to correctly ensure rotational invariance and, as announced, (iii) the term with coefficient is indeed allowed, even if we restrain the model to nearest-neighbour interaction only.

In conclusion, the method developed in this paper—which consists in defining neighbour clusters based on the interaction range, in considering explicitly the edges of the lattice, in considering the local symmetries defined in terms of cluster environment to distinguish the various edge terms, etc.—provides a safe basis to further implement rotational invariance. This method could be applied to a large variety of lattices and to higher dimensions.

## Footnotes

- Received November 9, 2005.
- Accepted February 13, 2006.

- © 2006 The Royal Society