Delusive chirality and periodic strain pattern in moiré systems

Geometric phase analysis (GPA) is a widely used technique for extracting displacement and strain fields from scanning probe images. Here, we demonstrate that GPA should be implemented with caution when several fundamental lattices contribute to the image, in particular in twisted heterostructures featuring moiré patterns. We find that in this case, GPA is likely to suggest the presence of chiral displacement and periodic strain fields, even if the structure is completely relaxed and without distortions. These delusive fields are subject to change with varying twist angles, which could mislead the interpretation of twist angle-dependent properties.


INTRODUCTION
Geometric phase analysis (GPA), the analysis of the spatial variation of the phase of a periodic signal, like in a high-resolution image obtained with a scanning tunneling microscope (STM), has been widely applied to gain insight into the physics of correlated electron systems. Specific examples include phase-sensitive identification of a d -form factor density wave in cuprates [1], the observation of a multiband character of the charge density wave (CDW) in a transition metal dichalcogenide [2], and discommensuration and topological defects in ordered electronic phases [3][4][5][6]. GPA is also extensively deployed to correct image distortions [7,8], to extract chiral CDW displacement fields [9], to quantify local moiré lattice heterogeneity [10], spatial variation of strain [11][12][13], and periodic strain modulations in van der Waals heterostructures [14].
GPA has been implemented in different ways: using spatial lock-in (SL), also known as coarse-grained phase field extraction or Lawler-Fujita method [7,10,19,20] (see also Supplementary Material Section III), using realspace fitting [6], and combining Fourier-filtering with intensity thresholding [9]. These methods all extract the  [15,16], (b) monolayer MoS2 on 2H-NbSe2 (30°, at 77 K(> TCDW ) (adapted with permission from [17], copyright 2018 American Chemical Society), and (c) monolayer MoS2 on Au(111) (7.7°) [18]. Red, blue, and orange circles highlight lattice, satellite, and moiré peaks, respectively. spatial variation of the phase with respect to a reference signal corresponding to the center of an atomic or a CDW Bragg peak in Fourier space. They all involve a filtering step performed either in Fourier space, by masking a region of interest and inverse transforming, or in real space, by convoluting the raw data with a smoothing function. These filterings are equivalent and have a decisive impact on the phase field obtained.
Here, we discuss the limitations and potential pitfalls of GPA, focusing on deceptive deformation fields and non-existent chiral distortions that can appear. We show that ambiguity arises when the analyzed image is formed by several non-overlapping fundamental lattices. This is particularly relevant to STM images of moiré patterns.

EXPERIMENTS AND MODEL
Moiré patterns are commonly observed by STM when imaging two periodic structures with different lattice orientations and/or different lattice constants [10,14,15,17,[21][22][23][24][25]. In Figure 1, we show the Fourier transforms (FTs) of three high-resolution STM images of moiré patterns obtained on different heterostructures: graphene on WSe 2 (Figure 1(a)), monolayer MoS 2 on 2H-NbSe 2 (Figure 1 (b)) and monolayer MoS 2 on Au(111) (Figure 1(c)). A common feature of all three FTs is a high degree of wavevector mixing, with many intense peaks besides the fundamental lattice peaks (q Li ). These additional qvectors correspond to a linear combination of the constituent lattices wavevectors: q = i n i q Li , where n i are integers. Ultimately, we can generate all the FT components by choosing specific sets of {n i } coefficients. The peaks corresponding to the lowest order wavevector mixing, i.e. with the most zero components in {n i }, will typically appear brightest. In Figure 1, example lattice peaks corresponding to q Li are marked by red circles, while the shortest moiré wavevectors q M j , which correspond to the differences between the nearest q Li of the two lattices generating the moiré pattern, are out-ii lined by orange circles. Blue circles highlight example satellite or moiré replica peaks, whose wavevectors correspond to first-order linear combinations of q M and q L components. Since q M j are the shortest wavevectors in the system (in the absence of other periodicities than the constituent lattices), these satellite peaks are the closest to the lattice peaks. Note that one of the satellite peaks always coincides with one of the lattice peaks and that the lattice peaks corresponding to the top layer in a heterostructure are typically more intense. For clarity, we only highlight some of the peaks in Figure 1.
To investigate possible artifacts of the GPA technique, we design numerically generated real-space images whose FT intensity maps contain all the features discussed above. The simplest model consists in a sum of plane waves with properly defined wavevectors. To simulate the three-fold symmetry of the heterostructures shown in Figure 1, we define the q-vectors by considering three non-equivalvent (2π/3 rotated) lattice peaks q Li (i = 1, 2, 3) and their satellites q ij = q Li + q M j (j = 1, 2...N ). For a uniform notation, we define the lattice peaks as q i0 = q Li , indeed q M 0 = 0. In the following simulations, we allow up to six satellite peaks (N = 6) placed symmetrically around each lattice peak to best mimic our data. Note that in general, any number of (satellite) peaks could be included. However, we are only interested in those that are the closest to the lattice peaks as all the others are filtered out when applying GPA to real data. It is for the same reason that we do not consider the peaks (q M j ) close to the origin (q=0). Once the wavevectors of the first lattice peak and of its six satellites have been defined, e.g. q 1j , we obtain the others by 2π/3 rotations around the origin to account for the spatial symmetry of the systems we are simulating. The corresponding topographic signal can then be calculated as: (1) where T i (r) is the signal from the i-th lattice peak and its satellites. To simplify the discussion we set ϑ ij = 0 without impacting our final conclusions. It allows a convenient description in terms of peak intensities in the FT (which are determined by the amplitudes A ij of the plane waves in real space). Considering ϑ ij ̸ = 0 can give a more general description with additional complexity that may be worth investigating in specific systems. Using the identities given in Supplementary Material Section VIII, we obtain A ij (r) sin(q Li r + q M j r) = A i (r) sin(q Li r + φ i (r)). (2) The terms A i (r) and φ i (r) satisfy: and where ∆q M kl ≡ q M k − q M l . This is an amplitude-modulated signal, which is a spatial analog of the beating phenomenon well-known from acoustics. In real space images, it appears as a periodic domain structure defined by the q M j wavevectors, similar to a moiré pattern (Supplementary Material Section I). ). We can explicitly calculate this phase field and the corresponding displacement field, which then allows us to extract a strain field as described in Supplementary Material Section II. In the next section, we do this in two ways. First, we use the phase field obtained from the three main lattice peaks and their satellites to compute the displacement field, solving an over-determined system of equations in the least square sense (Supplementary Material Section II). Second, we show what happens if we only use pairs of two lattice peaks and their satellites, which is in principle sufficient to describe any real distortion in the two-dimensional STM image plane.

DISCUSSION
The key elements influencing the GPA are the satellite peaks, which ultimately determine the nature of possible displacement and strain fields. When all satellite peaks have the same intensity, the corresponding topography shows domains due to the amplitude beating (Equation 3), but the displacement field is strictly zero (Equation 4, Suppl. Fig. 3). However, for any asymmetry in the satellite peak intensities, we not only find an amplitude beating that shapes a domain structure in real space, but also periodically varying displacement and strain fields as shown in Fig. 2 and in Suppl. Figs. 4 and 5. These fields are modulated with a periodicity of 2π/q M and their fine structure depends on the details in the satellite peak intensities as discussed below.
To give an insight into the variety of possible displacement and strain fields depending on the satellite peak structure, we consider a simple case where all the satellite peaks have the same intensity, except for one (the k-th) which is stronger: A ij = C < A ik ∀i and ∀j ̸ = k (j, k ̸ = 0). When the wavevector corresponding to the stronger satellite peak (q ik ) is parallel with q Li , Equation 4 yields a radial displacement field. This field is pointing away from (towards) the center of the beating domain when the satellite peak wavevector is shorter (longer) than q Li , corresponding to tensile (compressive) strain within each domain (see Suppl. Fig. 4).
Whenever q ik is not parallel with q Li , Equation 4 yields a chiral displacement field whirling around the center of each domain (Fig. 2). The left-handed or righthanded vorticity depends on the angle between q ik and q Li (see Fig. 2 and Suppl. Fig. 5). Note that this chiral displacement field is zero in the middle of the domain and increases toward the domain edges, exactly as observed for the chiral CDW displacements in TaS 2 [9].
The analysis so far was purely analytical, evaluating Equation 4 for a given set of A ij amplitudes and q ij wavevectors. It shows that whenever the satellite peak structure has an inhomogeneous structure, delusive displacement, and strain fields appear in a system that is perfectly relaxed and undistorted by construction. Such inhomogeneities can be present for several reasons in real STM images of moiré structures, but in any case, because all satellite peaks are not equivalent since one of them always coincides with a lattice peak. Consequently, there is a significant risk of misinterpreting experimental data in terms of displacement fields and strain which are not present in the actual material.
To illustrate possible strain misinterpretations further, we perform GPA by applying the popular SL technique used to restore and analyse experimental STM data to our numerically generated images. This approach enables us a perfect control over the input parameters and assess their impact on the GPA. We indeed observe the same effects, finding finite periodic displacement and strain fields (Suppl. Fig. 6) although the input structure does not have any of them. Unsurprisingly, since the displacement and strain fields depend on the satellite peaks, the response depends on the filter size applied to the data for the SL procedure. Fig. 3 shows how the displacement field progressively disappears when reducing the Fourier space included in the filter around the lattice peaks at q Li . Likewise, the displacement field will progressively vanish if the satellite peaks shift further away from the lattice peak while using a constant filter size. The moiré satellite peaks move away from the lattice peaks in heterostructures with increasing twist angle, and the resulting vanishing displacement field can be mistaken for a twist angle dependent reduction of strain [14].
An additional analytical element we consider in the following is that two linearly independent in-plane vectors should uniquely define any two-dimensional displacement. Hence, we should find the same displacement field for any choice of lattice peaks pair and their satellites (see Supplementary Material Section II). However, comparing (either the SL or the analytical) results obtained with the three possible pairs of q-vectors for our undistorted model heterostructure, we observe three very different patterns with symmetry-lowering strain and displacement fields in Fig. 4 (and in Suppl. Fig. 8), providing additional evidence that they are not real.
Let us finally consider the case of an STM image with a real spatial distortion. To this end, we use a bespoke periodic chiral displacement field to distort a perfect lattice as described in Supplementary Material Section VI. This chirality is hardly seen in the corresponding topography ( Fig. 5(a)), which is not surprising, as a pure distortion acts on the phase of the signal while the peakto-peak amplitude remains intact. This remains true even when we apply a distortion field so large that the atomic displacements become visually perceptible. However, the periodic chiral distortion is unambiguously seen in the corresponding FT ( Fig. 5(b)) as satellite peaks surrounding the lattice peaks. Their positions correspond to the q Li ± q M j wavevectors, where q Li and q M j are the wavevectors of the undistorted lattice and of the periodic iv displacement field, respectively. The magnified lattice peak region in the inset of Fig. 5(b) reveals asymmetric satellite peaks, as expected in the presence of a periodic displacement field.
Applying GPA to this distorted image, we recover precisely the displacement field we used to generate the distorted image (Fig. 5(c) and we can extract the corresponding strain (Fig. 5(d)). This result demonstrates the suitability of GPA to access genuine strain fields. Since we are dealing here with real strain and not an artifact, the result should not depend on whether we solve the over-constrained problem using all three lattice peaks as in Fig. 5(c) and (d), or whether we solve it using any pair of lattice peaks. Indeed, the results in Fig. 5(e) and (f), obtained using only two lattice peaks, are exactly the same, and this is true for any pair of lattice peaks as shown in Suppl. Fig. 9.
The last example suggests a possible method to disentangle real strain from artifacts by comparing the GPA results obtained with two different pairs of lattice peaks. Unfortunately, this may not work in experimental data of moiré systems, as the satellite peaks can have a combined origin. They may result from a combination of moiré wavevector mixing and periodic distortion. Without a suitable strategy -which we are not yet aware of-to dis-entangle these contributions, GPA-based strain analysis is ambiguous in moiré systems: one may find a periodic strain where there is absolutely no strain present in the system. The only reliable statement possible is in the case of pure deformation, where GPA analysis gives exactly the same result irrespective of the choice of lattice peaks for the analysis. Whenever different results are found, it is impossible to make any strong claim about the presence or absence of strain. This problem is formally very similar to the one discussed in the context of scanning transmission electron microscopy of composite materials, where a variation in the unit cell composition leads to an apparent strain at the interface between two compounds although there is no strain [26].
A further upshot of our analysis is to highlight a risk often overlooked when applying GPA-based distortion corrections (notably the popular Lawler-Fujita method [7]) to scanning probe images. GPA has to be deployed with care and is only applicable when there is no contamination in the vicinity of the lattice peaks, which otherwise leads to deceiving apparent distortions; correcting for them is likely to introduce spurious features that are not genuinely representative of the sample under investigation.

CONCLUSION
We have analyzed possible pitfalls of geometric phase analysis to extract periodic deformation and displacement fields from scanning tunneling microscopy images when several atomic-scale lattices contribute to the imaging contrast. While our focus is mainly on the moiré patterns observed on heterostructures, the conclusions are general. We demonstrate the possible misinterpretation of high resolution images whose contrast is determined by mixed signals from at least two lattices, one of which may also be a charge density wave. They can mimic deformation, chirality and other displacement and strain fields, despite the system under investigation being perfectly undistorted and relaxed. FT of (a) showing the lattice peaks and their satellites arising due to the periodic distortion. The numbers identify three lattice peaks available for the GPA. The inset shows a zoom into one of the lattice peaks with its satellites. (c) Chiral displacement field and (d) biaxial strain field extracted by GPA using all three lattice peaks in the FT. (e) Chiral displacement field and (f) biaxial strain field extracted by GPA using only the lattice peaks #1 and #2. Note that the results do not depend on the choice of lattice peaks combination (see also Suppl. Fig. 9).

METHODS
The samples and data presented in Fig. 1 were prepared and measured at the University of Geneva as described in [15,16], [17], [18] for panel (a), (b) and (c) respectively.

I. AGGREGATED SIGNAL OF THE SUM OF PLANE WAVES
In this section, we examine simple cases to enlighten the meaning of Equation (3) and Equation (4). First, we consider only two q-vectors: q 10 = q L and q 11 = q L + q M ; and with amplitudes: A 10 (r) = αA 0 and A 11 (r) = A 0 . By substituting these into Equations (3) and (4), we find: which clearly shows that Equation (3) and Equation (4) describe an amplitude-modulated signal with spatially evolving phase, given by φ(r) with respect to a T 0 (r) = sin(q L r) reference signal. The high-amplitude domains of this spatial beating signal appear when the constituent waves interfere constructively, i.e. when they are in phase: which shows that the period of the amplitude and the phase modulations is given by q M .
We plot the sum of two one-dimensional sine waves showing the amplitude-modulated beating signal in Suppl. Fig. 1(a). The phase relation is highlighted in Suppl. Fig. 1(b) where we plot the beating signal together with the sin(q L x) reference signal whose amplitude is matched for a better comparison.
Suppl. Fig. 2 shows the phenomenon of spatial beating in two dimensions. In panel

II. DISPLACEMENT AND STRAIN FIELD
In this section, we discuss how to extract a displacement and strain field from topographic images. We consider a topography (T (r)) where the imaging distortions and local strain lead to a change in the original positions: r ′ = r + u(r).
We assume that the topography is composed of the sum of harmonic modulations: T (r ′ ) = i T i (r ′ ) = i sin(q i r ′ ). Taking only two modulations we can write: T 2 (r ′ ) = sin(q 2 r + q 2 u(r)) = sin(q 2 r + φ 2 (r))). with Writing component by component we get, which can be written in a matrix equation form as: where The solution of this system of equations is: The equations above show that we need to consider (at least) two peaks in the FT in order to extract the displacement field, otherwise, we get an under-determined system of equations.
On the other hand, we have not specified which two peaks. As long as they correspond to structural signals (atomic Bragg-peaks), we assume that the obtained displacement field represents structural strain and imaging artifacts.
After obtaining U(r) (either by the spatial lock-in or the real space fitting method [1][2][3][4][5]) we can correct the distortions of the topographic signal by mapping T (r ′ ) to T (r ′ − u(r)) which can be done e.g by the imwarp Material Section VII we show that when the obtained phase variation is not (purely) due to real distortion, using different peaks lead to very different results, e.g. artificial symmetry lowering, see Fig. 4 and Suppl. Fig. 8.
A modified method has been proposed (though, neither implemented nor compared to the original) by Benschop, Jong, Stepanov et al. [2], which involves using the phase of all three independent lattice peaks and solving the obtained over-determined system of equations in a least-square minimization sense. Mathematically, we start with: Writing component by component we get, q 1x u x (r) + q 1y u y (r) = φ 1 (r) q 2x u x (r) + q 2y u y (r) = φ 2 (r) q 2x u x (r) + q 2y u y (r) = φ 2 (r) (12) which can be written in a matrix equation form as: where Supplementary equation (13) is formally the same as the previous supplementary equation (8) for the displacement field. However, this time ⃗ Q is not invertible and we cannot find the solution as U(r) = Q −1 ⃗ φ(r). Instead, at each pixel, we can search for the U(r) which minimizes the Euclidean norm of the QU(r) − ⃗ φ(r) vector, i.e. a solution in the minimum of least-squares sense. This is easily implemented in Matlab using the mldivide (or backslash) [7] operation pixel by pixel.
The displacement vector field U(r) can be used to find the local distortion of the lattice and to define the strain tensor at each position whose elements are defined by the derivatives of the displacement field, as follows [8,9]: where U x (r) and U x (r) are the x and y component of the displacement field at position r.
The normal strains in x and y directions interfere with each other and lead to an isotropic strain which is defined by ϵ bi = (s xx + s yy )/2, commonly referred to as biaxial strain in the literature [10,11], whereas the anisotropic (uniaxial) strain is defined as ϵ uni = (s xx − s yy )/2.

IV. VARIOUS CASES WITHIN THE SUM OF PLANE WAVES MODEL
A. Numerical evaluation of the analytical form using three phases to determine the displacement field Suppl. Fig. 3. We find an amplitude-modulated topography (a) in the case when all the satellite peaks have the same intensity (b) (FT of (a)), but the displacement field (c) is identically zero.
Suppl. Fig. 4. If one satellite peak is stronger than the others and is parallel with the lattice peak there is a radial displacement field: pointing away from (towards) the middle of the domains. In In all cases we find the same as with the evaluation of the analytical formula, apart from the borders where the Gaussian kernel is not yet entirely within the image. In the displacement vector field and biaxial strain maps, the dashed black squares delimit the region 3σ (in real space) away from the edges of the images, where the filtering kernel is entirely within the image and edge effects disappear.
Suppl. Fig. 6. Displacement and biaxial strain fields extracted using the spatial-locking method in the same cases that we studied analytically. (a)-(d) One peak is stronger and parallel with the lattice peak, radial strain (Suppl. Fig. 4). (e) and (f) One peak is stronger (Fig. 2).

V. ANALYTICAL RESULT FOR THE CASE DISCUSSED IN FIGURE 3
Suppl. Fig. 7. Strain field calculated using the analytical form of the phase for the case discussed in Figure 3 of the main text.

VI. GENERATION OF CHIRAL DISPLACEMENT FIELDS
In order to test the effect of a periodic chiral displacement field on a perfect single lattice, we first generate a chiral displacement field and then distort the lattice accordingly. In practice, we first create a base lattice as the sum of three plane waves with q-vectors rotated 2π/3 with respect to each other. For the displacement field, we start again from the sum of three plane waves, however this time with q-vectors which are shorter than for the base lattice; their length and direction can be selected on demand. Then we take the gradient of this field to obtain a vector field. Finally, we locally rotate the vectors of this field around the out-of-plane axis by a γ angle to obtain a chiral displacement field. With this method, we can continuously tune the displacement field from fully radial (pointing toward/away from the domain center) through chiral (whirling around the domains) to entirely tangential (perpendicular to the direction toward the domain center) by choosing the appropriate γ.
Finally, we can apply this displacement field (e.g using the imwarp [6] function in Matlab) to the base lattice to generate the distorted lattice.

PLACEMENT FIELD
A. In case of phase fields obtained analytically Suppl. Fig. 8. Even if the phase fields are obtained by evaluating the analytical formula, the obtained displacement and strain fields depend on the pair of lattice peaks chosen for the calculation.
In each row, left to right: FT of the real-space topography with orange circles marking the peaks whose phase was used to calculate the displacement field, the displacement vector field, and the biaxial strain field that we obtained in the analysis.