Number of equidistant neighbors on honeycomb lattice

A convenient scheme is presented for calculating potential energy of van der Waals interacting bilayer graphene and other similar 2D compounds. It is based on the notion of the existence of two types of local symmetry of carbon atoms ordering, a 3- and 6-fold one. Potential energy of an atom is expressed as a sum of contributions from rings of equidistant atoms on neighboring layer. Methods are described to compute the radius of rings of equidistant atoms and number of atoms they contain. Exact positions of atoms are found as well, allowing to apply the introduced method in modelling of anisotropic potentials and to be used when twisting between layers is present.

1 Van der Waals potential energy of an atom on graphene structure.
After the experimental realisation of single graphene layers [1] a growing interest arose in artificial 2-dimensional structures. Other newly discovered materials, named van der Waals heterostructures [2] become the subject of broad research. The reason of that is richness of physical phenomena observed and high potential of their possible applications [3]. Novel materials include graphene/silicene, graphene/M oS 2 , silicene/M oS 2 systems [3], as well silicene [4] or M oS 2 alone [5].
When two layers of graphene are rotated [6], [7] Moiré patterns result. They are related to formation of spatial patterns due to singularities in density of electronic states. In bilayer graphene twisted at about 1.06 • an unconventional superconductivity was found with critical temperature of 1.7K [8], [9], [10], [11]. The discovery led to formation of a new research direction named twistronics, where physical material properties are controlled by rotation of layers along with traditional methods such as chemical doping, electrical and magnetic field or pressure. Coexistence of insulating states mixed with superconductivity has been controlled by electric field and twist angle, offering a chance for theoretical bridging novel materials with high-T c superconductors [12], [13], [14].
For a proper understanding of the role of interlayer van der Waals forces in these 2D materials it is desirable to expand available theoretical concepts, for the use along the existing methods such as for instance these in large scale modelling of huge collections of atoms in molecular dynamics simulations [22]. With right computational algorithms there is no difficulty finding the number and positions of neighbors to carbon atom on graphene lattice. An insight into the role of symmetries, possible broader implications may be overlooked however when purely computer solutions are at use.
There are two types of relative ordering of two layers of graphene that are of special importance in any theoretical considerations. These are named AA ordering, when atoms of one layer are on top of atoms of another layer, and AB stacking when some atoms of one layer are located over the centers of hexagon cells of another layer. AB ordering is energetically stable and common in nature while AA order is observed in some rear, special cases only. When there is some disorder in arrangements of layers, or there is a rotation between their lattices, locally both orderings are found.
In case of AA stacking we will call the local ordering of atoms on neighboring layers a type I, and for AB stacking -a type II, as shown in Fig. 1 (type I is also present for AB planes while type II can be found for AB structures, only). Figure 1: Rings of equidistant atoms on neighboring planes of graphene. On the left is type I ordering, present for AA and AB stacking and type II (right) present in AB stacking, only.
Potential energy Φ I (z) for an atom in position AA at a distance z from the plane is given by an infinite sum of contributions from equidistant atoms of type I: where the number of neighbors N (i) of equal distance depends on distance projection ρ(i) measured within the plane and a is the unit cell length for in-plane atoms (a = 1.42Å in case of graphene). The potential energy between two single atoms, U , is any of potentials discussed often in literature (Lennard-Jones, Mie [17], Kolmogorov-Crespi [18], [19], Lebedeva's [20], [21], as reviewed in [15]).
In case of AB planes, there is equal number of atoms that are in ordering of type I and type II. Equation on energy for type II atoms, Φ II (z), is similar to Eq. 1 but the number of neighbors at equal distances is different. In case of AB ordering we must take an average of Φ I (z) and Φ II (z) as an average energy of an atom.
A quick convergence was found [15] of potential energy Φ as a function of the total number of atoms N in a symmetric molecule for type I and type II molecules: an approximately 1/N 2 scaling is found (or 1/R, where R is molecule radius), when the rings of equidistant neighbors are added in Eq. 1.
A convenient scheme may be created of calculating potential energy for different stackings by taking advantage of local symmetries of neighboring atoms. For that we need projection (in plane) of distances between an atom position over the plane together with number of atoms at these distances (i.e. atoms in rings of equidistant atoms, as these in Fig.  1). For symmetrical atom-atom interaction potentials there is no need for other information. When potentials of KC and Lebedeva are considered, relative orientation of atoms in x − y plane does not play a role either. When an interaction potential depends on orientation between atomic orbitals within x − y plane, like a one introduced recently by Carr et al. [23], exact coordinates of atoms are needed in computation, as provided by the method we propose.
2 Radius of rings of equidistant atoms and number of their atoms.

Radius of rings.
Positions (x, y) of carbon atoms in a graphene lattice of C-C spacing a = 1.42Å may be written as 1 : for any pair of integer numbers, (m, n) ∈ Z, except these when a Boolean function g(m, n) is true: g(m, n) ≡ (n mod 3 = 0) ∧ (m mod 3 = 1) ∨ (n mod 3 = 1) ∧ (m mod 3 = 2) ∨ (n mod 3 = 2) ∧ (m mod 3 = 0) = 0. (3) The condition (3) is equivalent to: Equations 3 and 4 are used for finding indexes (m, n) of positions that are in centers of hexagonal cells not belonging to honeycomb lattice. The condition that (m, n) belongs to graphene lattice is G(m, n): In Equation 2, the distance from the center (0, 0) of a point defined by indexes (m, n) is given by: Graphene structure should be considered, in general, as an object possessing a 3-fold rotational symmetry. However, locally, it acquires symmetry 3-or 6-fold, depending on which position is assumed as the center of rotation axis. Also, 6-fold symmetry is preserved when performing any lattice transformations that do not take into account lattice missing points described by Eq. 5.

Number of atoms in rings of equidistant atoms
We observe that the distance expressed by Eq. 6 has the following properties: If a pair (m, n) fulfills one of conditions given by Eq. 3 than one needs to check if other pairs (m ′ , n ′ ) fulfill Eq. 3, such that: since m 2 + mn + n 2 = n 2 + nm + m 2 . Also we find that the following pair needs to be checked for conditions (3): We observe also that replacing both signs of (m, n) at the same time also leads to the same distance in Eq. 6.
Let us define an identity transformation matrix 2 , M 0 : The above rules (7,8,9) can be expressed, respectively (10,11,12), in the following way, as matrix transformations of (m, n) vector,  Indexes (m, n) are displayed along atoms. Circles centered at (0, 0) pass through atoms belonging to rings of equidistant atoms. This lattice, with an atom located in rotational center has the symmetry properties of type I ordering defined in Fig. 1: a 3-fold rotational symmetry with a broken parity symmetry (a lack of mirror reflection symmetry, e.g. along 0Y axis). Lines and arrows are used for explanation of parity symmetry of some atoms.
The above rules may be combined. One finds that M 3 is commutative with M 1 and M 2 while the last two, M 1 and M 2 , do not commutate together. Hence, we may generate additional unique transformations. By defining M ji as M j × M i , we find: There are also tertiary transformations. Of these, M 312 = −M 12 , and M 321 = −M 21 , and we have two more left: The set of transformations (Eqs 10 -19) is complete, in the sense that any transformation of these within that set of any other transformation from within there leads to result that is within the same set of transformations.
Hence, since there are 12 unique transformation matrices, we conclude that there are maximum 12 only possible pairs of (m, n) with the same value of ρ, as given by Eq. 6 and summarized in Table 1. We may consider pairs (m ′ , n ′ ) that result by applying transformations M i to (m, n), as these listed in the second column of Table 1, as a unique and complete set of elements (group) of a set (group) of transformations we name M m,n .
Another approach of studying the number of atoms in a ring could go through classical rotations. For counterclockwise 60 • rotation the transformation matrix R xy 1 has the form (we are transforming coordinates of a point (x, y) of positions indexed by (m, n); We distinguish between matrix R xy 1 acting in x-y space and R 1 acting on indexes (m, n)): If to write x ′ = m ′ + n ′ /2, y ′ = √ 3n ′ /2, one obtains R xy 1 in a form acting on indexes (m, n): We find that the matrix R 1 of counterclockwise rotation for π/3 is the same as M 21 . In a similar way we may find that matrix of clockwise rotation for π/3, is the same as M 12 . Obviously, a matrix for 120 • rotation, R 2 , may be written as: For completeness, let us have explicitly form of rotation matrices for other multiplicities of 60 • : There is an equivalence between matrices M and R. Two subsets of transformations M i may be created: Every next element of subset A or B in Eq. 24, in order as defined there, is related to its previous element by transformation R 1 . The last elements transform to first ones. Transition between elements of subsets A and B is possible by transformation of indexes (m, n) ⇒ (n, m), that is by M 1 .
Therefore, we find one more way of creating a set of transformations of initial point (m, n) into points of the same distance, by using 60 • rotations on points (m, n) and (n, m). Figure 3 shows how angles (with respect to 0X axis) and difference between them change for points (m, n) and (m ′ , n ′ ), where (m ′ , n ′ ) = M 1 × (m, n), as a function of n/m ratio.
The parallel lines in Fig. 2 are drawn for a geometrical explanation of parity symmetry and for finding relation between indexes (m, n) for which parity holds.
The red parallel lines (30 • slope with respect to 0X axis) pass through points belonging to graphene lattice, only. The lines are shifted with respect to each other by vector of length 3 along the 0X axis (as shown by horizontal arrows). Any parallel line shifted for multiplicity of 3 along 0X passes through points belonging to graphene lattice, only. For all points lying on the red line passing through (0, 0) (all of them have the same m and n indexes) the parity transformation is valid: if (m, n) belongs to graphene then (m ′ , n ′ ) = M 3 × (m, n) belongs to graphene lattice. Hence, for any value of n, points with m = n belong to graphene lattice and fulfill the condition of parity symmetry. The right-shifted red line (arrows in Fig. 2 are shift vectors) passes through points shifted right for 3 by m and it passes points (3, 1), (3,2), etc.
Hence, any points (m, n) fulfilling the condition: where k ∈ Z, belong to graphene lattice and have a parity counterpart (m ′ , n ′ ) = (−m, −n) belonging to graphene lattice. Figure 2 shows lines parallel to these in red, in blue and light blue colors, passing through points m = n + 3k + 1 and m = n + 3k + 2. They are also periodic in m, with period 3, however periodicity is present for positive only or negative only values of m, i.e. there is a reflection symmetry for these lines with respect to value m = 0. Any point lying on a line m = n + 3k + i under transformation m ↔ −m and n ↔ −n finds its image on a line m = n − 3k − i, where i = 1, 2. However, points on blue lines m = n + 3k + 1 do not belong to graphene lattice while all their transform images belong to. In a similar and reverse way, all points on light blue lines, m = n + 3k + 2, belong to graphene lattice while their transformed images do not.
It is convenient now to change the notation. Let us name elements of sets defined by Eq. 24 by using rotation transformations: where R i M is i times repeated rotation operation R 1 on M 1 . When Eq. 25 is valid, all 12 available transformations (26) produce a 12-element ring of 6-fold rotational symmetry. In cases when parity symmetry is broken (m = n + 3k + 1 or m = n + 3k + 2), a 3-fold rotation symmetry holds. That means that from between 6 available rotation matrices R i in Eq. 26 three only transformations are allowed, and from 6 matrices R i M as well three only. Again, analysis of positions of lattice points in Fig. 2 may help in resolving the question on which transformations between graphene atoms are valid: points on m = n + 3k + 1 line are not valid graphene lattice points. However, points that are obtained by applying R 1 transformation to them are valid lattice points. Then two other points must be valid as well that are rotated by 120 • and 2 × 120 • . Therefore valid transformations in this case are R 1 , R 3 , R 5 . We see also that point R 0 M (i.e. M 1 ) is a valid lattice point and deduce in a similar way that transformations R 0 M, R 2 M, R 4 M are all valid. By using similar reasoning we find points fulfilling condition m = n + 3k + 2. Let us name these sets as C 1 and C 2 : Table 2: Pairs of (m, n) with the same distance ρ from the ring center, for any m, n < 25.
(m, n) (m ′ , n ′ ) ρ 2 (5,3) (7,0) 49 (6,5) ( There are special cases of (m, n). When m = n, there is no difference between points (m, n) and (n, m. This is seen also in Fig. 3, where α m,n = α m ′ ,n ′ = 30 • . Therefore, transformations R i M in Eq. 26 are the same as R i alone and the set of available 6 transformations is A given by Eq. 26. Another special case is when n = 0. One may check then (by using for instance simple operations on matrices) that in this case the results of some of transformations of subset A are identical to some transformations of subset B defined in Eq. 26. For instance, R 1 = R 0 M, R 3 = R 2 M, etc. In result, in that case we should distinguish between cases when m = 3k + 1 and m = 3k + 2. The allowed operations are subsets of these defined by Eqs. 27 and 28:

Multiply pairs of (m,n) for equidistant rings of atoms.
There are points (m ′ , n ′ ) which have the same distance ρ as (m, n) (Eq. 6) but do not belong to the set of operations M m,n . For instance: (7, 0) / ∈ M 5,3 and (7, 7) / ∈ M 11,2 . For ρ < 25 (Table 2), the only other pairs with the same ρ are these (16, 1) and (11,8), (19,0) and (16,5), (13,10) and (17,5). Another coincidence of that kind are pairs that are related by a trivial relation: if (m, n) generates an equidistant ring than (q · m, q · n) forms a ring that has radius q times larger, where q is any natural number. In Table 2 such a case is found for pairs (21,0) and (15,9), since they are related to (7, 0) and (5, 3) by common divisor of 3. At larger values of (m, n) such coincidences are found often. A few of them are triplets of pairs and many are quartets of pairs. When to consider pairs of co-prime numbers only, the first triplet is found for (80, 19), the first quartet for (25,23).
A'posteriori analysis of the data indicates on existence of a few rules governing the occurrence of multiply pairs of rings with the same values of ρ(m, n), as summarized in Table 3 ( Figure 4 demonstrates some of restrictions on allowed values of (m ′ , n ′ )).
Equation 6 may be written in another form: For pairs of (m, n) and (m ′ , n ′ ) of equal-distance ρ we will have: which may be written as:  1. The sum of indices m and n for rings of the same radius ρ fulfills the condition: The sum of indices m and n, and m ′ and n ′ , for rings of the same radius, (m + n) + (m ′ + n ′ ), differs by integer values p: , and maximum value of p, p max , fulfills the condition: p max < (2/ √ 3 − 1)ρ (Fig. 4). 4. The sum of multiplicities of indices m · n, and m ′ · n ′ , for different rings of the same radius, is less than 2/3 times ρ 2 : m · n + m ′ · n ′ < 2/3 · ρ 2 . We have also: m · n < 1/3 · ρ 2 and m ′ · n ′ < 1/3 · ρ 2 .
In Fig. 4 a) the quantity |(m + n) − (m ′ + n ′ )| is displayed as a function of ρ: discrete values of |(m + n) − (m ′ + n ′ )| = p are found, with p ≥ 1 and p ≤ p max , where p max is integer part of (2/ √ 3 − 1)ρ and lower than that value. The above condition may be used for finding allowed pairs (m ′ , n ′ ) by the method of trial and error, for any reasonably low value of p.

Transformations and graphs.
Many aspects of discussed transformation properties may be well illustrated in graphs. Figure 5 shows examples of transformation paths between points on equidistant rings.
In a) the path is from (3, −1) to (−2, 3), in the same order as in listing in Table 1. In b) it is in an ordered manner, as in set T in Eq. 26. In general, the transformation paths and number of visited points depend on the order of transformations, points (and paths) may repeat itself.

Type II symmetry.
In case of type I symmetry, the parity symmetry is broken (there is no symmetry of positions of atoms after reflection along OY axis, or change x ⇒ −x). In case of type II symmetry, with restored parity symmetry, one more geometrical degree of freedom is available. In result, the number of atoms in rings of the same radius as in case of type I symmetry may become two times larger.  Table 1. Regardless of the starting point of transformation, all points on the ring are visited. In general, however, the transformation paths and number of visited graph vertices depend on the order of transformations. In b) transformations are ordered, i.e. each subsequent is a rotation of the previous one by 60 • , except the path between points (3, −2) ⇒ (1, 2) (red arrow) which uses (m, n) ⇒ (n, m) transform between two subsets of 60 • rotations. Table 4: Summary of available transformation operations on co-prime pairs of numbers (m, n) (including n = 0 and n = m: m ≥ n, n ≥ 0) preserving the invariant m 2 + mn + n 2 . N is the number of available transformations (for type I ordering; In parentheses -for type II). Here, k is an integer number, k ≥ 0. Indexes (m, n) refer to type I ordering.

Conditions
Type I Type II N Eqs m = n + 3k, n > 0 (A, B) -12(−) 26 m = n + 3k + 1, n > 0 C 1 (A, B) 6(12) 27 m = n + 3k + 2, n > 0 C 2 (A, B) 6(12) 28 m = n (R 0 , · · · , R 5 ) -6(−) 26 m = 3k, n = 0 (R 0 , · · · , R 5 ) -6(−) 26 m = 3k + 1, n = 0 C 1,n=0 (R 0 , · · · , R 5 ) 3(6) 29 m = 3k + 2, n = 0 C 2,n=0 (R 0 , · · · , R 5 ) 3(6) 30 Figure 6 illustrates that a transition between type I and type II surroundings may be achieved by introducing a displacement operator D acting on indexes (m, n). It increases m by 1, while a reverse operator, D −1 , decreases m by 1: Any transformation operations on indexes of atoms discussed so far for type I symmetry may be carried out on by transformations DM i D −1 on atoms of type II symmetry. The only difference in results between I and II symmetries is due to conditions on belonging to graphene lattice (Equations 3-5). These conditions must be checked again in case of operations on atoms of type II. Table 4 summarizes results on available transformation operations derived for both types of ordering. A column lists there the number of possible graphene atoms in rings, N . In parentheses is the number of atoms in type II ordering. For any ring of type II of the same radius as a ring of type I the number of graphene atoms in it is either the same, it is two times larger, or it is none.  In-plane projection of distances between an atom on one graphene plane and nearest neighbors on next plane, for type I of neighbors, as in Fig. 1 (left), found in AA and AB stacking of layers, and for type II ordering, as in Fig. 1  3 Summary and conclusions.
A scheme is presented for calculating potential energy of van der Waals interacting bilayer graphene and other similar 2D compounds. It is based on the notion of the existence of two types of local symmetry of carbon atoms ordering, a 3-and 6-fold one (type I and type II). Potential energy of an atom is expressed as a sum of contributions from rings of equidistant atoms on neighboring layer. Methods are described to compute radiuses of rings of equidistant atoms and number of atoms they contain. The methods are based on analysis of transformation properties of equation on atoms arrangement on equidistant rings and on symmetries of underlying crystallographic lattice. Table 4 summarizes available transformation operations and number of atoms on rings for both types of orderings for pairs of co-prime numbers (m, n) indexing atoms positions in lattice. Table 5 lists the data obtained for equidistant rings of radius up to about 12Å. Not only the number of neighbors is found but their exact coordinates as well (parameterized by indexes (m, n)). That allows to apply the introduced method to modelling anisotropic potentials, when potential energy depends on angles between bonds, and in particular it allows to compute properties of structures formed by twisted bilayers of graphene (Fig. 7). For any atom its neighbors may be described with an infinite series of concentric rings. For atoms located at superstructure lattice nodes the centers of rings of atoms from both layers coincide, with arrangement of atoms rotated between layers. When to neglect distinction between atoms on different layers than a 6-fold symmetry of Moiré pattern is observed. There are however differences between arrangement of atoms on equidistant rings (blue and green concentric circles in Fig. 7) formed around vertices of unit cells, restoring the 3-fold symmetry of the structure.  , 1)). Full green circles represent atoms of unrotated layer and blue empty ones of rotated one. Arrows indicate basis vectors forming Moiré superstructure. When to neglect distinction between atoms on different layers than a 6-fold symmetry of Moiré pattern is observed. There are however differences between arrangement of atoms on equidistant rings (blue and green concentric circles) formed around vertices of unit cells, restoring the 3-fold symmetry of the structure.