Moiré, Euler and self-similarity – the lattice parameters of twisted hexagonal crystals

A real-space approach for the calculation of the moiré lattice parameters for superstructures formed by a set of rotated hexagonal 2D crystals such as graphene or transition-metal dichalcogenides is presented. Apparent moiré lattices continuously form for all rotation angles, and their lattice parameter to a good approximation follows a hyperbolical angle dependence. Moiré crystals, i.e. moiré lattices decorated with a basis, require more crucial assessment of the commensurabilities and lead to discrete solutions and a non-continuous angle dependence of the moiré-crystal lattice parameter. In particular, this lattice parameter critically depends on the rotation angle, and continuous variation of the angle can lead to apparently erratic changes of the lattice parameter. The solutions form a highly complex pattern, which reflects number-theoretical relations between formation parameters of the moiré crystal. The analysis also provides insight into the special case of a 30° rotation of the constituting lattices, for which a dodecagonal quasicrystalline structure forms.

1. Introduction

Recently Cao et al. (2018) demonstrated that stacked graphene layers with relative rotations can have drastically different properties to their regularly aligned counterparts. They observed that a relative rotation of about 1.1° leads to superconductivity in double-layer graphene. This corresponds to one of the `magic angles' previously predicted by Bistritzer & MacDonald (2011), who calculated the band structure of twisted double-layer graphene and found that the narrowing of bands at small angles is non-continuous, and at 1.05° and other distinct angles the Dirac-point velocity vanishes. These publications, and several others dealing with rotation-controlled band-structure modifications in double-layer systems (see e.g. Suárez Morell et al., 2010; Moon & Koshino, 2012; Fang & Kaxiras, 2016; Trambly de Laissardière et al., 2012; and Kerelsky et al., 2019) opened up the field of `twistronics' which currently attracts considerable scientific interest. Also at large rotation angles, twisted bilayers bear striking effects, such as the formation of a dodecagonal quasicrystalline phase at a rotation angle of 30°.

In order to fully understand the physical effects in twistronics and to allow precise device design, it is imperative to understand the relation of the relative rotation angle of the constituting layers and the resulting structure of the moiré crystal. This is particularly important owing to the fact that both the physical properties – viz. the occurrence of distinct `magic angles' – as well as the lattice parameter of the moiré crystal, and therewith the local order and the length scale of its modulation, are critically dependent on the rotation angle. Recently the field evolved and, apart from graphene bilayers, dichalcogenides such as MoS₂ and WS₂ (Chen et al., 2019; Lu et al., 2019) as well as heterogeneous systems, e.g. stacked BN and graphene layers (Finney et al., 2019) or MoS₂ on graphite (Chen et al., 2019), are now also considered.

In this paper a straightforward and presumption-free framework for the direct calculation of the angle dependence of the moiré-crystal lattice parameter is provided. The paper has two aims: First, to establish a set of equations that describe the moiré lattice parameter as a function of the rotation angle of the constituting lattices that can straightforwardly be used for spreadsheet calculation or be implemented in compact computer code for everyday use, e.g. in the design of twisted-bilayer devices or the interpretation of high-resolution transmission electron micrographs of twisted bilayers. For example, the equations can be used to identify angles of particularly critical angle dependence, or for back-calculation of the rotation angle by measurement of lattice parameters of actual devices in the transmission electron microscope for quality control or reproducibility checks. Second, to point out the salient relations between the pattern of solutions (hereafter the `solution pattern') of moiré-crystal lattice parameters and mathematical number theory, which may imply a corresponding relation between number theory and the actual properties of twisted bilayer devices.

2. Phenomenology

2D hexagonal lattices can be represented by the periodic arrangement of a rhombic unit cell in a plane, with two lattice vectors a and b oriented at an angle of 120° (Fig. 1). This can be regarded as the projection of the hexagonal Bravais lattice along the sixfold axis. When two such hexagonal lattices are superposed and rotated relative to each other (Fig. 2), moiré patterns occur. These patterns, at angles below and above 30°, have the same rotational symmetry and consist of area-filling rhombuses as the constituting hexagonal lattices, but they are rotated with respect to the latter and the periodicity length, in the following referred to as the moiré lattice parameter, is larger.

Figure 1 2D hexagonal lattice with its unit cell (gray) and the lattice vectors a and b.

Figure 2 Moiré patterns at different angles and dodecagonal quasicrystal at 30° formed by relative rotation of two hexagonal lattices.

Starting at low angles of relative rotation, the pattern develops in the form of broad zones of high and low net-plane density. These zones are broad at low angles, because the net planes are almost parallel in the constituting lattices. The rotation center appears bright because at the common point of rotation the net planes converge and do not fill up the space in between. Further bright zones appear on a line through the rotation center approximately perpendicular to the a axis of one of the lattices, in the following referred to as the reference lattice, and approximately perpendicular to the b axis of the other lattice, referred to as the rotated lattice. The angle between these lines, in the following referred to as moiré lines, is 60°, corresponding to the narrow angle in the rhombic cell used for the description of a hexagonal structure. The symmetry of the underlying lattices demands that further moiré lines parallel to the lines crossing the rotation center occur, which consequently form a rhombic moiré pattern.

As the rotation proceeds, the bright low-density zones move towards the rotation center and new bright zones occur far out and move towards the axes of the constituting lattices, thus forming a moiré lattice of increasingly smaller rhombuses. Careful inspection reveals that the angle of the moiré line sets, with respect to the normal of the a and b axis of the reference and rotated lattice, respectively, increases by half the rotation angle of the underlying lattices and as a result the rhombuses of the moiré pattern are always oriented with the narrow angle towards the rotation center.

Upon further rotation, the lattice parameter of the moiré pattern decreases continuously, until, at angles approaching 30°, which corresponds to half the periodic angle, it approaches the length scale of the underlying lattices. Around 30° no discernible moiré pattern is visible, but highly complex structures form. This is the case in a range of about 8° below and above 30°, i.e. from about 22° to 38°. At precisely 30° a dodecagonal quasicrystalline pattern occurs. At angles larger than 38°, a moiré pattern becomes visible again, developing from smaller to larger moiré lattices with increasing rotation angle. The angle of the moiré line sets with respect to the normal of the a and b axis of the reference and rotated lattice, respectively, is now given by half the rotation angle of the underlying lattices plus 30°, and therefore the moiré pattern still appears as lattice of rhombuses with the narrow edge oriented towards the rotation center.

The lattice parameter of the moiré lattice can be determined by calculating the distance from the rotation center to the center of the first bright zone on the line crossing the rotation center. We will refer to this first bright zone as the first-order moiré point. The second bright zone will be referred to as the second-order moiré point etc. Owing to the discrete nature of the lattices, the lines of the underlying lattices do not cross perfectly at one point in all the moiré points.

We will refer to moiré points, in which the lines cross perfectly (or very close to perfectly) as `clean' moiré points. For an overview of the nomenclature and coordinate systems used, see Fig. 3.

Figure 3 Coordinate system and nomenclature for reference and rotated hexagonal lattices (here rotated by α = 4.3°). The vector pointing at the first moiré point is shown in green.

3. Calculation of the moiré lattice parameter

As a basis of the reference lattice we choose the vectors

and for the rotated lattice

The moiré lattice resulting from the rotation has basis vectors a_M, b_M. Let us choose the direction for a_M as perpendicular to a₁ at α = 0, so that in the lattice section shown in Fig. 3, vectors a_M can be represented by a positive linear combination of a₁, b₁, and a₂, b₂.

In order to find the moiré lattice parameter, we thus have to find the first point on the moiré line, which is common to both the reference lattice and the rotated lattice. Let us express a vector to a point of the moiré lattice in terms of the reference lattice and the rotated lattice, respectively, as

where n, m, n′ and m′ are integers. We are only interested in vectors to moiré points on the moiré line crossing the origin, i.e. having very specific length and direction. This narrows down the set of coefficients n, m, n′ and m′ to a subset that fulfils specific interrelations. Careful analysis of rotated lattices at various angles establishes that moiré points only form for coefficients fulfilling the following conditions simultaneously:

and

where p is the order of the moiré point. These empirical results will be verified below, when all coefficient interrelations have been established.

Note that relation (6) holds for angles α < 30°, and using it restricts the results to this angular range. At angles α > 30°, it becomes 2m − n = p, and all following calculations can be carried out analogously for that angular range.

Using (3) and (6), the length of the moiré lattice parameter can be calculated directly. For the first moiré point p = 1 holds, and we obtain

In order to establish further relations between the coefficients let us now consider the angle β between the vector a_M, expressed in terms of the reference lattice, and the a₁ lattice vector. Calculating the scalar product , we obtain the relation

This relation can immediately be used to identify candidates for moiré lattice points: for a given angle β we have to find combinations of integers n and m fulfilling (8) within a defined error margin, which then, using (3), can be used to calculate the corresponding vectors a_M in the basis of the reference lattice.

Analogously, we consider the angle of the vector , expressed in terms of the reference lattice, and the a₁ lattice vector. Using the basis (2), we calculate the scalar product and obtain the relation

Since both scalar products calculated relate to the a₁ vector of the reference lattice, they both describe the same angle. Therefore (8) and (9) can be equated, which gives us a condition for the identification of moiré points as lattice points common to both the reference and the rotated lattice. Because the denominators are independent of the rotation angle, they have to be equal by themselves for a solution of the whole equation to exist. Equating the denominators gives us the condition n² + m² − nm = n′² + m² − n′m. This can be fulfilled only if

Further solving the equation then reduces to equating the numerators, which finally gives us a relation between the rotation angle and the coefficient n and m:

For first-order moiré points, using (6) with p = 1, we then obtain

Combining this result with (7), the angular dependence of the length of the moiré lattice parameter is obtained as

For small angles α, the sine can be very well approximated by its argument, which leaves us with the expression (α in radians) or, in degrees, .

Note that with (10), all interrelations between the coefficients of (3) and (4) are established. It can straightforwardly be shown that (5), (6) and (10) solve the equation of the magnitudes of a_M and , which confirms retrospectively that the empirically deduced conditions (5) and (6) are valid.

For a given rotation angle, we can now use (12) to calculate m, then calculate n using (6) under the condition that both m and n are integers. Then (3) is used to calculate the basis vector a_M of the moiré lattice and its length, the moiré lattice parameter.

Table 1 displays integer solutions for moiré lattice parameters, calculated according to this procedure accepting a deviation of ±0.005 for n and m from being integers and sampling the rotation angle at 10⁻⁵°. The full set of solutions between 0.1 and 30° is displayed in Fig. 4.

Table 1
Integer solutions of equation (12) and lengths of moiré lattice vectors

α (°)	n	m	aM/a1
0.21	158	315	272.80
1.09	31	61	52.83
2.13	16	31	26.85
2.87	12	23	19.92
3.15	11	21	18.19
5.09	7	13	11.27
6.01	6	11	9.54
9.43	4	7	6.08
13.17	3	5	4.36
21.79	2	3	2.65

Figure 4 Moiré lattice parameter (in units of the reference lattice parameter) as a function of α. Red: exact solution; blue: approximation aM = 180/πα; black: length of actual moiré-lattice basis vectors for which m and n are integers (Table 1).

Note that the last entry involves the indices n, m = 2, 3, implying n′ = 1, which is the lowest possible combination of positive integers ≠ 0 and corresponds to an angle of 21.8°. This is in agreement with our previous observation that in the range between 22° and 30° no moiré lines can be seen.

Fig. 4 compares this continuous description with the length of actual basis vectors of the moiré lattice, for which m and n are integers. This graph can be understood in the following way: the continuous line according to (13) describes the lattice parameter of all moiré lattices, as we apparently see them developing in a set of rotated lattices. If we take a closer look, however, we see that for most angles the first-order moiré points are not clean – the net planes of the reference and the rotated lattice do not meet in a single point, but just get more or less close.

The black circles in Fig. 4, corresponding to the solutions of (12) for which m is integer, are the moiré-lattice parameters of clean moiré points, in which the underlying net planes meet in a single point. The angle dependence of moiré points is thus continuous, while the angle dependence of clean moiré points is discrete.

For the sake of completeness, let us introduce the b-lattice vector, b_M, for the moiré lattice. We choose the basis such that b_M is perpendicular to the b axis of the rotated lattice, b₂, for α = 0 (Fig. 5). This gives us a basis that is not right-handed, but naturally suits the resulting moiré lattice. In particular, the moiré points corresponding to positive linear combinations of a₁, b₁, and a₂, b₂ also have positive indices in terms of the moiré-lattice basis.

Figure 5 The moiré lattice spanned by the basis aM, bM (green arrows) and its indexing using the coefficients p and q.

The b-axis vector of the moiré lattice in terms of the reference lattice is then given by

which together with (3) forms the basis of the moiré lattice.

All points of the moiré lattice A_M can then be expressed using two indices p and q, where p is the previously introduced order of the moiré points along the a axis, and q is the corresponding parameter along the direction of the b axis, as

The moiré lattice for a given rotation angle of the constituting lattices can thus fully be described by the following procedure: for a given rotation angle α, use (12) to calculate m, then calculate n using (6) for p = 1. With (3) and (14) the moiré-lattice basis for the given angle is obtained, which, using the indices p and q, spans the full moiré lattice (Fig. 5).

4. Higher-order moiré points and moiré crystals

Let us now consider moiré crystals, i.e. we add one or more atoms to the unit cell of the constituting lattices (i.e., in crystallographic terms, we add a basis). If we want to consider moiré crystals we necessarily have to take into account higher-order moiré points, as can be seen from the following example.

Fig. 6(a) displays a moiré lattice at α = 6.6°. The moiré points of order p = 0, 1, 2, 3 can clearly be seen. A closer look reveals, however, that while for the second-order moiré point (blue circle) the lines of the lattices cross in a single point, allowing an unambiguous identification of its exact position, this is not the case for the first-order point (red circle). For the latter, the crossings of the net planes along the a and b directions of the rotated lattices do not coincide anywhere. This means that the second-order moiré point is clean but the first-order moiré point is not.

Figure 6 Moiré lattice (a) and graphene moiré crystal at a rotation angle of 6.6°. The circles display the lattice planes and local atomic arrangement at the first-order (red) and second-order (blue) moiré point.

This has a striking effect if we consider crystals rather than mere lattices: in Fig. 6(b) a moiré crystal resulting from two graphene layers rotated by the same angle of 6.6° is displayed. In the first-order moiré point (red circle) a completely different local atomic arrangement to that in the second-order point (blue circle) now appears. The atomic arrangement in the second-order moiré point corresponds to that of the origin, and therefore the distance between these constitutes the moiré lattice parameter for the given rotation angle. The moiré crystal unit cell, correspondingly, is given by the green rhombus in Fig. 6(b), and has twice as long a lattice parameter than a first look at the lattice alone would suggest.

The example illustrates a general property of moiré crystals: at a given rotation angle between the constituting lattices, the lattice parameter of the resulting moiré crystal is given by the clean moiré point of lowest order at that angle.

In order to describe the moiré crystal we thus have to calculate the higher-order moiré points, and then figure which of those are clean. For this we combine (11) and (6) which yields

The length of the vectors to higher-order moiré points is calculated directly from (3) and using (6) we obtain

and combining (16) and (17) yields the dependence on the rotation angle α,

for a given order p. We will refer to a_M,p as the length of the higher-order moiré lattice vectors, the angle dependence of which is now represented by a discrete family of curves. The curve for p = 1 corresponds to the line in Fig. 4.

The lattice parameter of the moiré crystal for a given angle is determined by finding the solutions of (16) for integers m, which, for most angles, will yield more than one solution for different orders p. The solutions for 1 to 30° taking into account the first 30 orders is shown in Appendix A, Fig. 9.

From these solutions, the one with the lowest order corresponds to the moiré lattice parameter; all higher orders are mere multiples of the latter and can be neglected. With the so-determined values for m and p, the length of the moiré lattice parameter is calculated using (17), and the moiré lattice vectors by (3) and (14) using the value for n obtained through (6).

Fig. 7 displays the so-determined angle dependence of the moiré-crystal lattice parameter for angles between 1 and 30° up to the 30th order. Also shown is the family of curves (red lines in Fig. 7) representing the continuous solutions of (18). The positions of the higher-order multiples of the moiré-crystal lattice parameters are indicated as breaks in the red curves.

Figure 7 Moiré-crystal lattice parameters (in units of the reference lattice parameter) as a function of α for p = 1 to 30. The continuous angle dependence of the moiré lattice points are shown as guide to the eye. The lines are broken at the positions of higher-order multiples.

Each dot in the solution pattern in Fig. 7 represents a moiré crystal, which forms at the respective angle with the shown lattice parameter. The green circle, for example, marks a point at the angle 11.98°, which means that here a solution is found for a moiré crystal with a lattice parameter of 19.16a₁ at a moiré point of fourth order. At this angle, there exists no clean solution at a moiré point of lower order that would lead to a shorter lattice parameter. A slightly lower rotation angle of 11.64° leads to a moiré crystal with a lattice parameter of 14.80a₁ (at a third-order point), and a slightly larger angle of 12.20° to a moiré crystal with a lattice parameter of 23.52a₁ (at a fifth-order point). The blue circle marks the example shown in Fig. 6, a moiré crystal at 6.6° at a second-order point.

The diagram shows that moiré crystals exist for all angles, but for some of them moiré points of very high order have to be taken into account. These high-order points correspond to very large moiré lattice parameters that may exceed the size of the flake of sample material, which implies an ultimate upper limit for the relevance of the high orders. For graphene, which has a lattice parameter a of about 2.5 Å, a typical consistent flake size of 1 µm corresponds to solutions of about 80th to 100th order.

The diagram also shows that the moiré lattice parameter critically depends on the rotation angle. In some ranges, variation of α leads to small variations of the moiré lattice parameter, but more often small angle variations lead to considerable changes in lattice parameter. This is for example the case for the angles for which first-order solutions exist, e.g. at 22°. A small change in angle leads to a change to a large lattice parameter of the highest order considered. Accordingly, the sequence of moiré-crystal lattice parameters upon variation of the rotation angle can be seemingly erratic.

Fig. 8(a) depicts the evolution of the moiré lattice parameter for the case of graphene around the magic angle of 1.1° (blue line), taking into account orders of up to 8. Within the small angular range of 0.035° the moiré-crystal lattice parameter takes on various, greatly different values between 130 and 1050 Å. Even if only second- or third-order solutions are taken into account, the lattice parameter can take values varying by a factor of almost 2 or 3, respectively. The figure also shows that at the magic angle of 1.1° no first-order solution exists. The closest solution is of 7th order and corresponds to a lattice parameter of about 910 Å. The closest first-order solution, on the other hand, occurs at an angle of 1.085° and corresponds to a lattice parameter of about 130 Å. This comes close to the value mentioned by Cao et al. (2018). Fig. 8(b) is a similar depiction of the situation at higher angles, 13° to 20°, taking into account orders of up to 16. The critical angle dependence, which includes changes of the lattice parameter of up to 160 Å upon minor changes of the rotation angle, is clearly seen.

Figure 8 Evolution of graphene moiré-crystal lattice parameters as a function of rotation angle, (a) around the magic angle 1.1° (blue line) taking into account orders up to 8, and (b) in the range from 13 to 20° taking into account orders of up to 16.

The solution pattern of the moiré-crystal lattice parameters depicted in Fig. 7 is separated into intervals which are defined by the angle positions of the first-order solutions. The most obvious interval in the figure extends from 13.2 to 21.8°, narrower ones extend from 9.4 to 13.2°, from 7.3 to 9.4° etc., and an incomplete interval extends from 21.8° to higher angles. At the limits of the intervals the solutions apparently diverge, which can be described by a set of pole functions (Section B2 of Appendix B). Note that each interval contains an identical but scaled version of the same solution pattern. The whole angular range can hence be described by repeated, scaled versions of a single-interval pattern. Conversely, the features identified for the pattern in one interval are valid for all other intervals as well. In the following we will use the pattern between 13.2 and 21.8° for further analysis.

A distinct feature of the solution pattern is that it displays several features that reflect a close relationship between mathematical number theory and the structural properties of moiré crystals, which in turn may be critical for their physical properties. This is due to the fact that the discrete series of solutions along the continuous curves of equation (18) (see Appendix A, Fig. 9) are interrupted at the limits of the intervals, i.e. at those angles at which a lower-order solution exists. We can thus immediately understand that those red lines, the order of which is a number with many factors, contain a lower density of solutions. On the other hand, along those lines, the order of which is a prime number, the density of solutions is highest. This is a direct consequence of the fact that for orders corresponding to prime numbers no lower-order solutions exist, except p = 1. Since in each interval defined by the presence of first-order solutions, p solutions of pth order exist (see Appendix A), we find sequences of n_p − 1 solutions along these lines, where n_p is the respective prime number. The lines corresponding to an order which is a square of a prime have a number of gaps corresponding to that prime, and generally those corresponding to powers r of a prime number n_p have n_p^r−1 gaps.

For lines the order of which does not correspond to a prime number, there are gaps in the sequence of solutions which correspond to solutions of their factors and solutions of the factors of their factors. For example, each interval for the line corresponding to the sixth order contains two solutions, since the two solutions being multiples of the two third-order solutions and the multiple of the second-order solutions are omitted. Lines corresponding to an order that has many factors, e.g. 10, 12, 24 etc., hence have particularly low densities of solutions. For a given order p, for all solutions that have a common divisor with p, a solution of lower order exists. Therefore, the number of solutions in each interval corresponds to the number of integers for which the greatest common divisor with p is equal to 1. Or – in other words – the number of solutions in each interval is equal to the number of integers that are co-prime with p. The number of solutions in each interval as a function of the order p is thus given by Euler's totient function φ(p) (see Appendix B, Fig. 11). Indeed, if one counts the solutions on for each order, the characteristic sequence 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, … is found. The properties of the solution pattern are discussed in more detail in Section B2 of Appendix B. Notably, in the solution pattern of the moiré lattice parameters, which itself follows Euler's totient function, each solution defines a subset of solutions, which again follows Euler's totient function and thus establishes a salient type of self-similarity. This is explained in more detail in Section B3 of Appendix B.

5. 30° case

For the special case α = 30° the resulting pattern is not a moiré lattice but a dodecagonal quasicrystal. The quasicrystal is nonperiodic and thus has no finite lattice parameter. This corresponds to the fact that no exact solution for a moiré-crystal lattice parameter exists for the 30° angle. If solutions for orders of up to 100 are calculated, the closest one deviates about 3 × 10⁻³% from 30° and has the order 97, and for the case of graphene would correspond to a lattice parameter of 8.7 µm. The lowest-order solution that falls well below a 1% margin off 30° is of 26th order and has a lattice parameter of 125 Å in graphene.

While the quasicrystal has no translation symmetry, it does possess a scaling symmetry, referred to as inflation symmetry (Socolar, 1989). This means that the structure can be reproduced by scaling with a certain factor, i.e. it possesses a certain type of self-similarity. For the decagonal case the scaling factor is . A second smaller scaling factor for this lattice exists, which leads to a reproduction of the structure when an additional rotation of 15° is taken into account. This second scaling factor is given by .

The scaling factor λ alternatively can be expressed as 1/tan(30°/2). Therefore the coefficient m for the first-order moiré point in the 30° case (12) equals . With this value, we can calculate a hypothetical moiré lattice parameter for the 30° case using (7) and obtain . The hypothetical lattice parameter of the moiré lattice at 30° thus corresponds to the smallest scaling factor of the dodecagonal quasicrystal lattice.

6. Discussion and conclusions

In this paper a real-space approach for the calculation of lattice parameters of moiré crystals formed by the relative rotation of two constituting hexagonal crystal layers is worked out. This can be applied to the case of graphene or metal-dichalcogenide crystals. A closed and consistent framework for the description of moiré crystals and their structural parameters is provided, with solutions that are straightforward to implement. In the literature, preceding papers including aspects of the present work are available. These, however, mainly focus on band-structure calculation and the interlayer electronic coupling (Shallcross et al., 2009, 2010) and restrict themselves to first-order moiré crystals. The focus of the present paper, on the other hand, is rather on the discussion of the higher-order moiré crystals and the structure of their solution pattern. One early paper focusing on geometric aspects of moiré crystals was presented by Rong & Kuiper (1993), who carried out scanning tunneling microscopy (STM) of [0001] graphite surfaces. They observed a region with a top layer rotated by 2.1° with respect to the bulk, and identified it as a moiré superlattice with a period of 66. Indeed this result corresponds to the clean moiré point at 2.13° with a lattice parameter of 65.78 Å, listed in the third row of Table 1 (scaled using a lattice parameter of 2.45 Å for graphite as in the reference). This pioneering paper also includes an expression for the continuous angle dependence of the first order moiré points, which corresponds to equation (13) in the present paper. In more recent STM work on twisted WS₂ bilayers (Chen et al., 2019) a relation between twist angle and moiré period is established, which compares well with our results for first-order moiré crystals. Higher-order moiré lattices were considered by Lopes dos Santos et al. (2007), who calculate the Fourier components of the hopping amplitudes and show that in the low-angle limit only first-order solutions are relevant for the corresponding physical properties.

The results of the present paper imply that the angle dependence of the moiré lattice parameter is critical. This directly relates to the required fabrication precision for the rotation angle of twin-layer devices. The current state-of-the-art for graphene bilayer systems is about 0.02° (Hill, 2019), and within this margin the resulting device may have considerably varying lattice parameters. If the exact value of the lattice parameter is critical for the functionality of the device, this may explain the low reproducibility of working devices, which is quoted as `… in three months of trying, just 2 of the 30 devices worked' (Hill, 2019). On the other hand, the critical angle dependence allows for a highly accurate determination of the rotation angle in existing bilayer systems by means of transmission electron microscopy, via measurement of the lattice parameter, which can thus be used for quality control and reproducibility checks.

In summary, the present work allows for the following conclusions:

(i) The discussion of superstructures formed by rotated hexagonal structures requires discrimination between moiré lattices, which apparently possess a continuous angle dependence, and moiré crystals. For the latter, a basis is taken into account in addition to the constituting hexagonal lattice, and requires a clear identification of `clean' moiré points.

(ii) For moiré crystals, the local atomic order in a non-clean moiré point can be significantly different from the reference point, i.e. the rotation center at the origin.

(iii) The angle dependence of the moiré-crystal lattice parameters has discrete solutions and higher-order moiré points have to be taken into account. Solutions exist for all angles, but small changes of the rotation angle can lead to drastic changes of the resulting lattice parameter.

(iv) For the upcoming field of twistronics, i.e. the deliberate fabrication of twisted bilayer devices with tailored physical properties, it is important to know the exact relation between the rotation angle and the structural parameters, in particular the lattice parameter, of the resulting devices.

(v) The angle dependence of the moiré-crystal lattice parameter has to be considered critical, which has consequences for the required fabrication precision for bilayer devices.

(vi) The low reproducibility in the production of magic-angle devices may be a direct consequence of this critical angle dependence. Investigating functioning magic-angle devices by transmission electron microscopy to determine the actual angle is highly recommended.

(vii) The solution pattern of moiré-crystal lattice parameters is highly complex and reflects a close relation between mathematical number theory and the formation parameters of the moiré crystal, and thus potentially also the physical properties of twisted-bilayer devices.

(viii) The formation of a twelvefold quasicrystalline structure at 30° corresponds to the case of a hypothetical first-order moiré lattice parameter taking the value of the smallest inflation factor of the dodecagonal lattice.

(ix) In the intervals defined by the positions of the first-order solutions the number of solutions as a function of the order p follows Euler's totient function φ(p).

(x) A novel type of self-similarity is found: in the solution pattern of the moiré lattice parameters, which itself follows Euler's totient function, each solution defines a subset of solutions, which again follows Euler's totient function.

Finally it should be noted that the present treatment only considers the purely geometric case of rigid and non-interacting constituting lattices, in which any rotation of the lattice translates to exactly defined atom positions uniquely determined by the angle. In a real bilayer system relaxations may take place and locally lead to more favorable relative atom positions and thus an overall lower energy of the system. Potential relaxation pathways were recently discussed by Kalinin (2020), and scanning transmission electron microscopy investigations on homo- and heterogeneous bilayers of MoS₂ and WS₂ (Weston et al., 2020) and on MoSe₂ and WSe₂ heterostructures (Rosenberger et al., 2020) indicate that atomic reconstruction indeed takes place. Local rearrangements may also lead to locking of favorable states under rotation and thus to a less continuous angle dependence than in the purely geometric case. In order to investigate such effects further analysis, e.g. by means of density-functional theory, is required.

Figure 9 Length of higher-order moiré lattice vectors (in units of the reference lattice parameter) as a function of the rotation angle α for p = 1 to 30; black: length of actual higher-order moiré lattice vectors for which m and n are integer.

APPENDIX A

Structure of the diagram of higher-order moiré points

The angle dependence of the higher-order moiré lattice vectors calculated according to (18) and the length of the vectors pointing at higher-order moiré points are shown in Fig. 9. This figure corresponds to the (differently scaled) Fig. 3(a) in Lopes dos Santos et al. (2007), where the first six orders are displayed.

Fig. 9 is the basis for the determination of the lattice parameters of moiré crystals, for which at each angle only the solution of lowest order is regarded.

The general structure of the solution pattern in Fig. 9 can be understood as follows: the first-order solutions are solutions for all higher orders as well, because the lattice parameters of the higher-order patterns are multiples of those of the first order. These sets of solutions make up the vertical sets of points, dividing the α axis into intervals. The interval from 13 to 20° is shown in Fig. 10, and similar intervals, becoming narrower at lower and wider at higher angles, exist between all first-order solutions. Within each interval, there exists one solution for the first order, two solutions for the second order, three for the third and so forth. Since all solutions correspond to periodicities of moiré crystals, all second-order solutions are solutions for the fourth order, sixth order etc. Generally all pth order solutions are solutions of the n × pth order as well (blue arrows in Fig. 10). These rules define the general structure of the pattern.

Figure 10 Section of the pattern in Fig. 9, displaying the coefficient m for each point and arrows indicating the relation between lower-order points and their higher-order multiples.

APPENDIX B

Structure of the solution pattern for moiré crystals

B1. Euler's totient function

The solution pattern for moiré crystals is subdivided into intervals defined by the angles at which first-order solutions exist. It is argued in the text that within each of these intervals the number of solutions is given by Euler's totient function φ(p). This is illustrated in Fig. 11, displaying the number of solutions in the interval from 13.2 to 21.8° as a function of the order p for the first 80 orders in the typical representation of φ(p). The red line represents the prime numbers n_p for which φ(n_p) = n_p − 1.

Figure 11 The number of solutions of a given order as a function of the order p, within an interval between two first-order solutions. The line represents the prime numbers. The number of solutions is described by Euler's totient function φ(p).

B2. Pole functions

In Fig. 7 it can be seen that at the angle of each first-order solution an apparent pole exists, at which the angle dependence of the higher-order solutions diverges. This is very obvious at e.g. 9.4°, 13.2° and 21.8°, but at lower angles narrower poles also exist for each first-order solution. We will refer to the angle dependence of the sequences of solutions forming these apparent poles as pole functions. They can be analytically described in a very similar way as the continuous description of the moiré lattice vectors (18) and its approximation by a function of the form of an hyperbola. In particular, the pole functions can be described as

where α_p is the angle at which a solution of order p exists and s is an integer. The equation thus does not only describe the pole functions corresponding to first-order solutions, but the more general family of curves corresponding to any solution.

The red line in Fig. 12(a) represents the pole function calculated according to (19) for the case s = 1, formed at the first-order solution at α_p = 21.8°, which is located on the right-hand edge of the diagram and is marked by a red circle. The blue and black lines represent the pole function for the second-order solution at 16.4° (blue circle) and the fifth-order solution at 14.3° (black circle), respectively, both for s = 1. The cyan lines represent the pole function for s = 2 of the third-order solution at 15.2° (cyan circle).

Figure 12 Pole functions for the moiré crystal solution pattern in the interval from 13.2 to 21.8°. (a) Examples red: p = 1, s = 1; blue: p = 2, s = 1; cyan: p = 3, s = 2; black: p = 5, s = 1. (b) Family of pole functions for p = 1 (red and green), p = 2 (blue) and s = 1 to 30.

In Fig. 12(b) the families of pole functions for s = 1 to 30 corresponding to the first-order solutions at 13.2° (green) and 21.8° (red), and the second-order solution at 16.4° (blue), are shown. It is obvious that the solutions are all located on crossing points of the curves, but not all crossing points entail a solution. This would be the case if instead of solutions for moiré crystals (cf. Fig. 7), solutions for moiré lattice vectors (cf. Fig. 9) were considered.

B3. Self-similarity

For each family of pole functions the sequence of solutions on a given pole curve reflects the divisibility of the coefficient s of the family member. Again, Euler's totient function is hidden in this graph, which can be seen in Fig. 13(a): if you count the number of solutions in the subdivision for each s on the pole functions corresponding to the first-order solution at 13.2° (red lines and point), you obtain the sequence 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, …, i.e. the number of solutions follows Euler's totient function φ(s). In the same graph, the same is found for for the pole functions corresponding to the first-order solution at 21.8° (green lines and point). This recurrence of Euler's totient function is also found for higher-order solutions: Fig. 13(b) shows the family of pole functions corresponding to the second-order solution at 16.43° (blue lines and circle), where the number of solutions in the subdivisions and is given by φ(s). The remaining solutions in the subdivisions and fall on the family of pole functions corresponding to the first-order solutions (red and green), and for both individually the number of solutions is given by φ(s). Fig. 13(c) shows the same for the third-order solutions at 15.2° and 17.9° (cyan lines and circles). In the four corresponding subdivisions, the solutions individually follow φ(s), and the remaining subdivisions are covered by lower-order solutions (here second- and first-order), each individually following φ(s). This can be continued ad infinitum for each solution in the pattern – each solution defines a family of pole functions and a set of angular subdivisions, in each of which the solutions follow φ(s), and the rest of the pattern is subdivided into ranges covered by solutions of the same and lower orders, for which the same holds.

Figure 13 Self-similarity and Euler's totient function (see text).

We thus find a novel type of self-similarity in the solution pattern of the moiré lattice parameters: in each interval defined by the poles at the first-order solutions the number of solutions for each order p is given by φ(p), and each of those solutions subdivides the angular interval into further sub-intervals, in which the number of solutions for a given s of the corresponding pole-function family again is given by an Euler totient function φ(s).