Non-Abelian topological defects and strain mapping in 2D moiré materials

Transmission electron microscopy (TEM) provides an experimental route to characterizing local atomic configurations in real space, including detecting the change in ${\bf u}$ across a domain wall, known as the Burgers vector. Changes in the stacking order are distinguishable by a dark field imaging technique that consists of inserting an aperture into the diffraction plane around a single Bragg position and recording the resulting filtered real space image. Depending on the choice of diffraction peak, contrast between domains (see Fig. 2(a)) or the partial dislocations that form the domain walls (see Fig. 2(b)) in the bilayer stacking are visible [20]. The Burgers vectors of the dislocations are exactly determinable, as the vector perpendicular to the diffraction peak for which their contrast vanishes in the dark field image [20, 26]. The Burgers vector of a dislocation in the bilayer is equal to $\Delta{\bf u}$, the change in shift vector across the boundary.

We note that the dislocations form a network with distinct topology. Twisted bilayer graphene (Fig. 2(b)) has a structure where six dislocation lines meet at a node. Despite being a different material with different symmetries, MoS₂ slightly twisted from a 3R-like stacking (close to 0^∘ twisting) appears to share the same dislocation network topology as graphene (Fig. 2(c)). In contrast, MoS₂ twisted from a 2H configuration (close to 180^∘ twisting) has a structure where three lines meet at a node rather than six (Fig. 2(d)). As we will explore later, the topology of the graphene and 3R-like case is defined by a punctured torus, originating from the single energy maximum in the stacking fault energy as a function of configuration (and two degenerate minima). However, for twisted 2H MoS₂ there are two energy maxima and one minimum, leading the configuration space to be described as a twice-punctured torus and generating a different topology. While in this paper we focus on the topology generated by the energy profile of bilayer graphene, we are motivated by the realization that a topological description of each system’s configuration space should distinguish the different possible connectivities of the network, including those arising from materials of different lattice symmetry, such as a square lattice.

While constructing the tools to analyze the topological defects of the moiré superlattice in configuration space, we also attain the ability to analyze the configuration distribution and extract the overall strain profile in the twisted bilayer. In order to do these analyses, we need to consider several important quantities of the moiré lattice: the displacement gradient matrix, the moiré vectors, and the Burgers vectors. All of these quantities give important information that allows us to reconstruct ${\bf{u}}$.

A displacement vector field $\bf{u}$ which describes the shift in positions of a lattice $({\bf{x}_{t}})$ compared to a reference lattice $({\bf{x}_{b}})$ such that ${\bf{x}_{t}}={\bf{x}_{b}}+{\bf{u}(\bf{x}_{b})}$, is closely related to strain. Typically, ${\bf{x}_{t}}$ is a 2D vector corresponding to the positions in the strained lattice, and ${\bf{x}_{b}}$ is a 2D vector corresponding to the positions in the intrinsic lattice, prior to application of forces. In the case of moiré materials, presented in Fig. 1, the reference lattice (with coordinates ${\bf{x}_{b}}$) is the bottom layer of the material, and the displaced lattice (${\bf{x}_{t}}$) is the top layer. Thus the vector field that produces a moiré pattern is in a sense analogous to the relative strain (heterostrain) between the two layers.

The linear strain matrix, used in the modeling of strain energies, is defined by taking derivatives of the ${\bf{u}}$ field and symmetrizing:

Note that this linear strain tensor removes rotation contributions to first order in angle, since rotating the lattice does not cost energy and thus should not be counted towards the total strain energy.

Without symmetrizing, the gradient of the ${\bf{u}}$ field is known as the displacement gradient matrix, $\overline{d}$, obtained from

The vector ${\bf{u}}$ can be recovered by integrating $\overline{d}$ over distance. If we assume that $\bf{u}(0)=0$ and that the strain is spatially uniform from $\bf{0}$ to $\bf{x}_{b}$, we obtain

Thus, ${\bf{x}_{t}}={\bf{x}_{b}}+\overline{d}{\bf{x}_{b}}$ or

Any arbitrary displacement of layer 2 relative to layer 1, formed by a combination of twist and strain, can be written in terms of a displacement gradient matrix, $\overline{d}$, with four independent components. It is useful to define the four components as follows:

where $\theta$ is a linear approximate of the twist (measured in radians), $\alpha$ is isotropic strain, and $\beta$ and $\gamma$ are uniaxial and shear strain, respectively [27].

The moiré length is the distance over which the lattice ${\bf{x}_{t}}$ has been shifted by a unit vector with respect to lattice ${\bf{x}_{b}}$, creating a local return to the starting configuration.

Consider the shift vector ${\bf{u}}$, as defined in Fig. 1. If we twist our lattices about a point where two atoms are stacked on top of each other, the origin has ${\bf{u}}=0.$ Traveling away from the origin in a direction ${\bf{x}}$, the displacement ${\bf{u}}$ increases until we reach the point where the lattices have diverged by a whole unit cell (${\bf{u}}=\pm{\bf{a}}_{1}$ or $\pm{\bf{a}}_{2}$, where ${\bf{a}}_{1}$ and ${\bf{a}}_{2}$ are two lattice vectors of the unstrained lattice) and thus are aligned (${\bf{u}}=0)$ again, resulting in moiré periodicity.

A general interface of 2D lattices with heterostrain, ${\bf{x}_{t}}=(1+\overline{d}){\bf{x}_{b}}$, can be characterized by a pair of moiré vectors, ${\bf{m}}_{i}$ ($i=$1, 2), which describes the moiré periodicity in the 2D space. Considering the corresponding pair of two coincident lattice points in the upper and lower layers, ${\bf{x}_{t}}_{,i}$ and ${\bf{x}_{b}}_{,i}$, respectively, the two moiré vectors can be expressed

where the constant vectors ${\bf{s}}_{i}$ can each be $\pm{\bf{a}}_{1}$ or $\pm{\bf{a}}_{2}$, depending on the coincident lattice condition for the moiré superlattice.

If ${\bf{s}}_{1}$ and ${\bf{s}}_{2}$ are collinear (without ${\bf{m}}_{1}$ and ${\bf{m}}_{2}$ being collinear), then the matrix $\overline{d}$ has determinant zero and the moiré pattern is 1D rather than 2D. We will ignore this case for now, and assume that ${\bf{a}}_{1}$ and ${\bf{a}}_{2}$ are each used once.

We can relate the lattice constants $\bf{a_{i}}$, the matrix $\overline{d}$ and the moiré vectors, by ${\bf{m}}_{i}={\bf{x}_{t}}_{,i}=(1+\overline{d}){\bf{x}_{b}}_{,i}.$ Thus if $\overline{d}$ is invertible, ${\bf{x}_{b}}_{,i}=\overline{d}^{-1}{\bf{s}}_{i}$ and

Putting Eq. 7 together with the corresponding equation for the other moiré vector, ${\bf{m}}_{i}$, we obtain a matrix equation that can be solved for $\overline{d}$ if we have measured $\overline{m}$ and $\overline{s}$:

where $\overline{m}$ and $\overline{s}$ are $2\times 2$ matrices formed by horizontally concatenating the ${\bf{m}}_{i}$ and ${\bf{s}}_{i}$ column vectors, respectively. Note that if the two ${\bf{m}}_{i}$’s are not linearly independent, that is again the $|\overline{d}|=0$ case.

The shift vector ${\bf{u}}$ is also related to the Burgers vector of the dislocations that form in the relaxed moiré system. To define the Burgers vector, first one considers a closed path along the lattice points of a perfect crystal, sometimes called a Burgers circuit. When a dislocation is introduced within this path the circuit becomes broken, and the vector connecting the now separated start and end points of the circuit is called the Burgers vector. In a moiré superlattice, the twisted interface acts as a dislocation, and we can obtain Burgers vectors by considering circuits that traverse the interface. The simplest closed circuits in the aligned, untwisted structure take the following form: travel from a lattice point $\bf{x}_{0}$ to $\bf{x}_{1}$ along the top layer, then down vertically into the bottom layer, then along $\bf{x}_{1}$ to $\bf{x}_{0}$ in the bottom layer, and finally return to $\bf{x}_{0}$ in the top layer by moving vertically upwards. To connect the Burgers vector $\bf{b}$ to the shift vector $\bf{u}$, it is easiest to think of obtaining the twisted geometry by keeping the bottom layer fixed and introducing a dislocation via rotation of the top layer. In this case, the failure of the circuit to close after applying the twist is equal to the relative change in positions of the lattice points that corresponded to $\bf{x}_{0}$ and $\bf{x}_{1}$ in the top layer. However this is simply the change in the local shift between the two points, e.g.

where the non-integral form is only true if one does not map $\bf{u}$ into the compact unit-cell torus.

While this expression was obtained under the assumption of a uniformly twisted system, the result is quite general. For the relaxed systems that occur in experimental devices, one finds that the Burgers vector is only non-zero if a “dislocation line” (that is, a domain wall) is enclosed in the circuit. Each dislocation line between a pair of AA nodes can be associated with a specific Burgers vector, as was done experimentally using DF TEM in Fig. 2. The pair of AA nodes also provides a moiré vector ${\bf{m}}_{i}$, which can be linked to ${\bf{s}}_{i}$ obtained from the Burgers vectors. Therefore, knowledge of the moiré vectors ${\bf{m}}_{i}$ and the Burgers vectors across the dislocations associated with ${\bf{m}}_{i}$ is sufficient to solve for the displacement gradient matrix using Eq. 8. As the moiré vectors and Burgers vectors can be measured from experimental dark field and diffraction images, it is possible to obtain a map of $\overline{d}$ from the information provided by DF TEM images, as we will discuss in Sec. IV.2.

Although a torus is topologically nontrivial, a small loop around the AA site (${\bf u}=0$) on the torus can be contracted to a point and thus does not inherit any nontrivial topological properties from the torus. The topological defects associated with winding around the holes of the torus are the dislocations [24], but this description alone provides no constraint on the manner in which the dislocations meet at the nodes. In the case of a graphene moiré superlattice, the three dislocation lines with different Burgers vector, colored red($R$), green($G$), and blue($B$), converge at the AA defect and diverge out again. What we need is a proper mathematical description of the topological nature of the AA nodes consistent with the experimental findings.

We start by investigating the distribution in configuration space of the order parameter for unrelaxed and relaxed moiré systems. Fig. 3(a) shows an unrelaxed twisted bilayer of a honeycomb lattice, and Fig. 3(d) shows the configuration space, colored via a scheme that is used consistently throughout this work. A Gaussian intensity distribution of a chosen color is centered around each key point in configuration space. The region of configuration space centered around the AA point is colored white, and those centered around AB (BA) are dark gray (light gray). Red, green, and blue regions are placed on the midpoint of the three equidistant shortest paths between AB and BA, two of which require crossing the unit cell boundaries. The coloring in real space in Fig. 3(a-c) is determined for each point, by determining the local shift vector, finding it as a point in configuration space, and adopting the color corresponding to that point.

Fig. 3(b-c) depict real space structures that have been modulated by a periodic lattice distortion to mimic relaxation, with two different amplitudes. As the amplitude of relaxation is increased from Fig. 3(a) to (b) to (c), light or dark gray regions, corresponding to nearly AB or BA stacking, take up increasing areas in real space, while AA regions shrink, and red, green, and blue regions evolve into lines. Because the red, green, and blue color tell us which path on the torus was used to get between BA and AB, i.e. $\Delta{\bf u}$, the color tells us the Burgers vector of that line. The corresponding distributions in configuration space are shown in Fig. 3(d-f). As relaxation strength increases, decreasing point density around the AA configuration reveals that fewer points in real space correspond to AA stacking, while configurations on the red, green, and blue lines, and especially AB and BA sites, become more numerous. Thus, the atomic lattice relaxation process in real space can be viewed as emptying out most of configuration space and populating only the AB and BA points and dislocation lines connecting them.

A similar emptying of AA and concentration at AB, BA, and the colored lines between AB and BA, occurs for the three strain types: isotropic, uniaxial, and shear. Fig. 4(a-c) show the configurations formed by ${\bf{x}_{t}}=(1+\overline{d}){\bf{x}_{b}}$ with a spatially constant displacement gradient matrix $\overline{d}$, corresponding to only one nonzero component $\alpha,\beta,$ or $\gamma$ in Eq. 5. Fig. 4(d-f) show corresponding structures after applying the modulation function to mimic atomic scale relaxation. We note that the order and orientation of red, green, and blue lines differs in real space for each strain component, as well as twist. In this sense, the colors (encoding local order parameter) provide information about which strain components are present.

To understand the algebraic structure of paths in the theta space (Fig. 5(c)), it is useful to consider the same path through configuration space in both the unit-cells shown in Fig. 1(e). Let us consider the full encirclement of the AA node shown in Fig. 5(d). We begin by observing that the BA $\rightarrow$ AB $\rightarrow$ BA transition through the $R$ and $G^{-1}$ dislocations in Fig. 5(d) becomes, in the parallelogram scheme shown in Fig. 6(a), a non-contractible loop about the torus’ periodic boundary conditions. The other BA $\rightarrow$ AB $\rightarrow$ BA transition through the $B$ and $R^{-1}$ arrows in the hexagonal scheme becomes another non-contractible loop in the parallelogram unit cell shown in Fig. 6(b). We label the first and the second transitions as $a$ and $b$, which will soon be identified with two generators of the free group $F_{2}$ [28, 29]. The third BA $\rightarrow$ AB $\rightarrow$ BA move through $G$ and $B^{-1}$ arrows of the hexagon can be decomposed as the inverse of $a$ followed by the inverse of $b$, or $a^{-1}b^{-1}$ [Fig. 6(c)]. A complete loop around the edges of the hexagonal unit cell is equivalent to the algebraic operation $aba^{-1}b^{-1}\equiv[a,b]$. Our convention is to perform the operation appearing on the left side of the product first.

The commutator $[a,b]$ in the language of the free group with two generators $a,b$ represents a “vortex” centered about the AA defect. The vorticity of this topological defect can naturally defined as the commutator $[a,b]$ (more detailed discussion in Appendix A). Similar to the conventional vortex defined in $S^{1}$, this vorticity defined for the AA node is non-trivial in the sense that it is not contractible to an identity, as illustrated in Fig. 6(d). After cancelling out the paths that are traversed both ways, the overall path for $[a,b]$ becomes equivalent to four partial loops around the four corners of the parallelogram, equal to a full loop round ${\bf u}=0$. On a torus such a loop can be contracted to zero and become trivial, but not for a punctured torus. The anti-vortex has the algebraic representation $[b,a]=bab^{-1}a^{-1}=[a,b]^{-1}$. Geometrically, this amounts to starting from the same BA point on the upper left corner of the hexagonal unit cell in Fig. 5(d) and making a complete counter-clockwise loop. Appendix A gives a more complete account of the vortex structures in the language of free groups with two generators $a,b$. One can find group-theoretic representations for vortex dipoles (vortex + antivortex) and vortex quadrupoles (two vortices and two antivortices) as well.

As we described in Section IIA, dark field TEM imaging can identify the dislocation lines with given Burgers vectors and the AB and BA domains separated by them. Each dislocation line converging on a given AA node can then be color-coded as $R$, $G$, $B$ or one of their inverses $R^{-1}$, $G^{-1},B^{-1}$ considering Burgers vector and the neighboring AB/BA domains in the TEM measurement. The free group language of the previous subsection gives a mathematically complete account of the vortex and antivortex structures, but it is helpful to translate the same statement to the more tangible and experimentally measurable $RGB$ scheme according to

By direct substitution we obtain the commutator

which is a product of transition vectors over the six domain walls in succession. In the same scheme we have the antivortex commutator

Any cyclic permutation of the six letters gives rise to the equivalent vortex or antivortex.

While the sign of the exponent in $R$, $G$, and $B$ operators can be obtained considering the order of the neighboring AB and BA domains by combining the DF TEM images with the first and second order Bragg peaks as shown in Fig. 2, there is a simpler scheme to assign the sign of the operators, considering AB/BA domains are always complementary. For example, if the six dislocation lines converging on an AA node appear, for instance, in the order of $RGBRGB$ while going clockwise around it, it ought to be interpreted as $RG^{-1}BR^{-1}GB^{-1}$ given in Eq. (11) and classified as a vortex. If the colors appear as $RBGRBG$, it is an anti-vortex according to Eq. (12). One only needs to keep in mind that the sequence of colors is to be understood as one color letter followed by the inverse of another color letter, and vice versa. Generalizations of the RGB scheme to vortex-antivortex dipole and/or vortex quadrupole structure are discussed in the Appendix B.

While moiré patterns and commensurated domain systems with vortex-type nodes have been studied extensively, those with antivortex-type nodes have not been demonstrated. We postulate this is due to the energy required to maintain sufficient strain, whereas twist and lattice constant mismatch can create vortex-type moiré without global strain. Nonetheless, we observe a line of antivortex nodes along a boundary of non-uniformly strained moiré superlattice.

Fig. 7(a) shows combined DF TEM images of a twisted bilayer graphene sample that contains a $\sim$1 $\mu$m sized bubble formed underneath the sample. Near the boundary of the bubble, non-uniform relative strain builds up in the moiré superlattice, which in turn creates the various strain components discussed in Eq. 5. The antivortices (vortices) can be identified by the $RBG$ ($RGB$) order in which the dislocations occur in a clockwise loop. The antivortices, which form along the top edge of a closed-loop dislocation, are each capable of annihilating with a nearby vortex, keeping the net winding number $w$ constant within a fixed-boundary region. In Fig. 7(b), loop A surrounds a vortex-antivortex pair, with $w=0$. Loop D, which surrounds the entire closed-loop dislocation, also has net $w=0$ as the entire feature could annihilate if the local strain were removed, in which case it would become similar to loop E. When circling a single antivortex (B) or vortex (C), the winding number is nonzero.

Existence of antivortices is a measure of the fact that anisotropic strains (uniaxial and shear) are dominating over the isotropic and twist components. Anisotropic strains alter the band structure and can produce pseudomagnetic fields [30]. We can quantify the various strain components from the displacement gradient matrix. Computation of the displacement gradient matrix from a DF TEM image is discussed in Appendix C. In brief, one component of the order parameter is known at every colored line, and an elastic model is used to interpolate in between, obtaining an estimate of the order parameter in the continuum. The displacement gradient matrix can be obtained by differentiating.

The density of vortices or antivortices is computed in Fig. 7(c), by taking the determinant of the displacement gradient matrix, $|\overline{d}|$, which by the definition in Eq. 5 is equal to $(\alpha^{2}+\theta^{2})-(\beta^{2}+\gamma^{2})$ in terms of the strain and twist components. If $|\overline{d}|>0$, vortices are present and if $|\overline{d}|<0,$ antivortices are present (see Appendix C). Thus, if the isotropic components (twist and isotropic scaling) outweigh the anisotropic components (shear and uniaxial strain), vortices form, and if the opposite, antivortices form.

We further use our estimated displacement gradient matrix to create strain maps of the three strain components, plus twist. Note that the moiré pattern only provides information on the heterostrain, or difference in strain between the two lattices. Unlike other methods to estimate heterostrain from the spatial structure of a moiré pattern [31], this DF TEM method includes knowledge of the lattice orientation and Burgers vector information, avoiding the need to make additional assumptions. Still, the need to interpolate within the moiré cell means that information smaller than the moiré scale is not deterministic from the data. In Fig. 8, strain maps are shown for a heterostructure of MoSe₂ and WSe₂. The intrinsic lattice constant mismatch of 0.3% should show up in the isotropic component. However, isotropic mismatch smaller than 0.3% is measured, indicating that the lattice globally strains to achieve closer to epitaxial matching, whereas global lattice mismatch is often assumed to be fixed when calculating moiré lengths. This type of strain mapping can reveal useful information about the phenomenology of moiré materials.