Considerations for extracting moiré-level strain from dark field intensities in transmission electron microscopy Open Access

Moiré materials, formed by stacking atomically thin materials with a slight interlayer twist or lattice mismatch, exhibit spatially modulated potentials that can localize excitons and charge carriers, serving as a basis for many novel devices.^1–11 With both highly tunable and proposed to be a physical realization of the Hubbard and other correlated models with experimentally accessible phase transitions, moiré materials are also a preferred platform to investigate and understand correlated electron physics outside of computationally tractable regimes.^12–34 However, as these materials locally reconstruct^35–40 in a manner that alters their electronic behavior,^41–46 a detailed understanding of the reconstruction mechanism in these systems is crucial to realize practical devices and further our understanding of their rich physics.

We discuss an interferometric methodology based on four-dimensional scanning transmission electron microscopy (4D-STEM)^47,48 that permits the direct measurement of inter-layer strain fields.^38,49 This approach involves fitting the intensity interference⁵⁰ of overlapping Bragg disks in the diffraction patterns to a known functional form to extract structural details.^38,49,51 Alternative 4D-STEM based methods typically fit the locations of Bragg disks^52–56 or have analyzed a combination of both intensity and disk locations.⁵⁷ However, determining the locations of these large converged beam electron diffraction (CBED) disks to the precision needed for revolving moiré-level strain is difficult without specialized apertures⁵⁸ and has not yet been feasible in marginally twisted structures. Previous works have also analyzed the Bragg disk overlap intensities in single-shot CBED patterns;^59–62 both for the focused probes, we consider and with large defocus values where holographic fringes encoding out-of-plane information are seen. However, these works had not yet fit these intensities to extract a local structure within a scanning setup, as is more characteristic of ptychography.

Other ptychographic techniques used in electron microscopy similarly analyze changes in the scattering intensities when scanning a converged probe and can provide much higher resolution,^63–66 but are currently arduous to acquire and analyze over the moiré length scale. Bragg interferometry (BI) avoids the same level of computational complexity by only fitting select areas in the acquired diffraction patterns of anticipated interest for determining moiré stacking to simple physically derived analytical expressions, omitting the central beam and averaging over the sample within the illuminated region to provide a coarse-grained picture of reconstruction on the moiré-scale. We refer to this heavily constrained analog of dark field ptychography, realized in x-ray diffraction as Bragg ptychography,^67–70 as Bragg interferometry³⁸ to reflect the lack of an attempt to reconstruct the probe wavefunction.

We outline how this interferometric methodology can measure strain from bilayer interfaces in generic structures and discuss and present solutions to various numerical challenges encountered in its practical application. Numerical approaches for the fitting and extraction of strain are compared, and various considerations are presented for improved strain resolution and analyzed using simulated multi-slice data and through re-processing published^38,49 data.

Likewise, we define $u x y bottom=( r x y bottom− r x y ref)$⁠. These displacements can be analyzed using small displacement theory (also known as infinitesimal strain theory),⁷¹ in which the strain tensors are linearized. This framework neglects higher-order terms from changes to the coordinate system, omitting the distinction between deformed and reference axes. Such an approximation is valid when the displacements change relatively slowly over the field of view such that $|∇ u x y top | ∞≪1$⁠. This necessitates that the relaxation is not too dramatic and that the lattice scale $a 0$ is much smaller than the moiré-wavelength, coinciding with a small interlayer twist. More precisely, the local deformation gradient within a rigid structure is on the order of $sin(θ/2)$ as $u x y top$ varies between 0 and $a 0/4$ over a length scale $a 0/sin(θ/2)$⁠. The following analysis, therefore, assumes $sin(θ/2)≪1$⁠.

The full derivation is shown in the supplemental information of Refs. 38 and 49 for a centro-symmetric homo-bilayer and a general bilayer, respectively. Briefly, a series of strict assumptions (that the beam is completely focused and that scattering takes place within a single plane and is well represented by the weak phase object approximation) are used to yield the simple analytic expression above. We also note that an equivalent expression for centrosymmetric materials may be obtained by taking the $Δf=0$ limit of functional forms derived for in-line CBED holography.^59–62 One possible avenue for increased generalization of these fitting functions may use the anticipated z-dependence of these expressions for $Δf≠0$ to additionally measure corrugations; however, we note that large defocus values strongly impact the subsequent discussion of resolution; therefore, such an approach is difficult in practice.⁵¹ 

The accuracy of this expression can be investigated through comparison to simulated electron diffraction patterns, obtained using the multi-slice algorithm.⁷³ In Fig. 2, we present multi-slice results obtained at a series of van der Waals (VdW) gap choices and multi-slice thicknesses, parameterizing the number of scattering events. Each multi-slice simulation is performed on a large uniform bilayer graphene sheet of a given stacking order to investigate the accuracy of Eq. (2) without resolution considerations, which will be discussed in Sec. IV. An $α$ of 4 mrad and a perfectly focused aberration-free lens are assumed. We note that the primary effect of changing $α$ in the absence of resolution considerations will be to increase the impact of small aberrations not included here. Aberrations will also change the internal structure of the overlap intensity, effects from which can be mitigated through averaging over the whole overlap region. The implications of a subtle internal structure in aberrated probes for non-spherical virtual apertures is outside the scope of this paper. The defocus also will have little impact when resolution considerations are ignored so long as the entire overlap region intensity is averaged.⁵¹ 

Figure 2(d) shows a slight disagreement between Eq. 2 and the multi-slice results. This discrepancy presents as a slight decrease in the multi-slice intensity of the second-order reflections in the AA-type stacking region and a small offset of the cosine origin. Increasing the VdW gap of the material leads to a greater deviation from the expected expression, as it is derived under the assumption that the free-space propagation present in a multi-slice between the two layers plays a minor role. However, the expression remains accurate up to at least gaps of 15 Å. We note that these discrepancies are far too small to be seen in the experimental data due to the effects of noise and finite resolution for materials with relatively small VdW gaps, but may need to be considered in bilayer materials with many intermediate layers between.

The slight discrepancy between Eq. (2) and the multislice results obtained using a slice thickness greater than the VdW gap (and, therefore, a single scattering event) can be attributed to the weak phase object assumption. The goodness of fit for the data simulated with a 4 Å slice thickness reflects that the weak phase object assumption is well motivated for bilayer graphene and slight deviations from Eq. (2) are primarily driven by multiple scattering. We conclude that the presented expression is sufficient for describing the intensity variations in graphene bilayers, and inaccuracies arising from finite spatial resolution, post-processing, and experimental considerations, such as sample purity, will dominate.

One major resolution limiting consideration in previous studies was the effect of beam-width biasing, as previous works reported stacking order and strain averaged over illumination regions of approximately 1–1.5 nm. These probe sizes are also larger than those seen in other electron microscopy techniques,^63–66 which often use probe widths on the order of 0.2-0.5 nm and also analyze the sample variability within illumination regions to obtain sub-angstrom resolution. Unlike many microscopy techniques, the minimum viable spatial resolution of BI is dependent on the sample morphology. To understand why we first note that the full-width at half maximum of the point spread function, or the modulo square of Fourier transform of $exp(iχ)$ where $χ$ is the Scherzer aberration function,⁷⁴ and the $α$-dependent phase shift introduced by the lens, yields a typical measure for probe width.⁷⁵ For many applications, the maximum tolerable phase error is then used to determine the maximal semi-convergence angle (⁠ $α$⁠) and its associated probe width of approximately 0.61 $λ/α$ where $λ$ is the wavelength of the electron probe [Fig. 1(e)].⁷⁶ We will not consider the effect of this phase error on the CBED pattern signal to noise as it is beyond the scope of this work, instead focusing on $α<10$ mrad where it is deemed tolerable with sufficient aberration correction⁷⁶ and consider the number of pixels we can obtain at each $α$ the practical limitations associated with their extraction. We note that in this discussion, we consider the information limit associated with fitting the average intensity within overlapping Bragg disks, instead of the information limit set by the detector pixel size, which may be realized in principle by considering the internal structure of these regions. However, we note that dark field intensities do not provide all of the phase information contained in the bright field. A larger resolution may also be obtained in theory through oversampling.

In BI, we often cannot realize the optimal $α$ as the large convergence angle will lead to an overlap of additional Bragg disks and prohibit the selective investigation of desired interfaces. To see how these considerations translate to the obtained spatial resolution, we first note that the probe wavelength we use is associated with a relatively low 80 kV accelerating voltage, which increases the probe width. The use of lower accelerating voltages is, however, more necessary in thin samples to reduce knock-on damage.^77–79 The resulting probe width provides us with a measure of our spatial resolution, where we can assess changes in average displacement down to a minimum spatial resolution of roughly d if not slightly smaller.⁸⁰ This can be used to determine the minimum $α$ needed at each twist angle to resolve stacking order domain size changes to a given precision. For this, we define the AA stacking domains to include stacking configurations of $| u x y|=0±0.25 a 0$ so that a rigidly twisted structure has AA domains of width $0.25 a 0/sin(θ/2)$⁠. The $α$ needed to resolve domain contractions to a given percent of the rigid value is shown in Fig. 3(a). We note that the BI methodology is most suitable for relatively small twist angles where dramatic changes to the AA domains are seen.

For many sample morphologies, $α$ must be smaller than 10 mrad. To see this, we show a series of anticipated overlap regions [Fig. 3(b)], obtained using the relations illustrated in Figs. 1(a) and 1(b). We can then straightforwardly compute the number of overlap pixels (in units of the camera length squared) at a given choice of sample morphology and $α$⁠. We also note that small overlap regions are difficult to isolate in practice. This currently limits the spatial resolution and the range of samples that can be investigated with this technique. However larger $α$ and, therefore, higher spatial resolution are possible than those used in previous works (⁠ $α≈2−3$  mrad)^38,49 if care is taken to define precise virtual apertures of variable shape [Fig. 3(c)] instead of avoiding overlap with encapsulating layers entirely. We note that for many samples, this will require virtual apertures to be defined separately for each pixel (or at least local regions of real space) instead of from diffraction patterns averaged over larger domains as was done previously^38,49,81 and that precise automated determination of these small overlap regions is difficult for the same reason that the precise disk locations in these individual CBED patterns are hard to extract. This limitation may also be addressed by fitting regions associated with a greater number of overlapping disks to more complex expected expressions so that larger virtual apertures can be used, although such an approach requires improvements in the computational post-processing to be discussed.

We compare multiple fitting procedures to obtain the displacements $u x y$ prescribed by Eq. (2). One approach we use involves fitting $u x y$ and the coefficients $A i, B i$ sequentially such that $u x y$ at each pixel location $(x,y)$ and the coefficients for each of the twelve disks can each be fit independently and parallelized. In addition to being much faster and allowing for the coefficients to be optimized using linear least squares, we found that this use of alternating optimization⁸² was also more robust to experimental noise and permitted the use of many initial guesses at each iterative step, helping to avoid local minima.

The fitting of $A i , B i$ and $u x y$ is then iterated until convergence. We compare this staged optimization to a joint optimization in which the coefficients and $u x y$ are optimized simultaneously following one iteration of the same multi-start u $x y$ optimization with fixed $A i=1, B i=0$⁠. The mean — $u x y$— obtained over the full set of pixels is shown for both approaches in Fig. 4(a) for artificial data corresponding to a rigidly twisted structure with disk intensities that exactly obey Eq. (2) with independent Poisson noise associated with acquisition included for each disk intensity. Representative experimental datasets are shown in Fig. 4(c), showing that the extent of noise included is comparable or greater than that expected of weakly scattering samples, such as graphene, and much greater than that expected of transition metal dichalcogenides (both acquired using $α=1.7$ mrad). Figure 5(a) shows that the iterative approach consistently achieves a lower error in the displacements $u x y$ [Fig. 4(a)]; however, this may not be robust when considering subtle differences in the initial conditions. Despite this, the increased speed and ability to use many initial guesses at each interactive step are still beneficial. Avoidance of local extrema and noise mitigation are both active areas of research, and alternative algorithms, such as genetic algorithms,⁸³ could prove better than this iterative multi-start approach, but were not required for experimental data of this nature.

The displacements obtained in the fitting procedure cannot yet be differentiated to assess local strain. This is because the obtained $u x y$ are not uniquely determined by the data since adding integer multiples of the lattice vectors $a 1$ and $a 2$ to $u x y$ results in the same set of intensities due to the sample’s symmetry. For centro-symmetric materials, flipping the sign of $u x y$ also results in the same pattern. Owing to these ambiguities, the raw displacement vectors obtained are not smoothly oriented and cannot be naïvely differentiated, even without experimental noise or imperfections in the fitting procedure. To resolve this, we can determine the sign and lattice vector offsets needed to obtain smoothly varying displacements. While not specific to BI-obtained displacements, such a procedure is not needed when analyzing simulated data with uniform moiré domains as the effect of lattice symmetry can be accounted for by easily removing low-frequency modes. Additionally, this degeneracy does not feature in more conventional strain-mapping procedures that monitor gradual shifts in the Bragg disk locations. The problem will, however, generally apply when attempting to differentiate experimentally obtained displacements given that the extent of noise is large enough that we cannot naïvely choose locally optimal offsets.

We note that this physical degeneracy is from a phase ambiguity in the projected electrostatic potential and that phase ambiguities, such as these, feature prominently across a wide array of physical techniques, from remote sensing to magnetic resonance imaging. A great deal of previous research has attempted to resolve these degeneracies through phase unwrapping procedures, often but not always cast as optimizing a non-local cost function.^84–90 These existing strategies for phase unwrapping we are aware of primarily find a smoothly varying scalar quantity $ϕ x y u n w r a p= ϕ x y$ $+2π n x y$ given $ϕ x y$ with $n x y∈Z$⁠, inverting $tan − 1⁡(sin(x)/cos(x))$ to resolve the argument degeneracy of $cos(x)$ or similar.

We seek instead to resolve a degeneracy in $cos 2⁡( x⋅ g ^ i)$⁠, therefore unwrapping a vector quantity along directions set by normalized reciprocal lattice vectors $g ^ i$⁠. Importantly, there is also an additional sign degeneracy in a squared cosine; therefore, the strategies developed for interpreting phases proportional to $u⋅ g ^ i$ in geometric phase analysis⁹¹ or general 2D phase unwrapping techniques^92–94 cannot be used directly. This sign degeneracy is a manifestation of a more general issue involving a loss of phase information in overlapping dark field reflections compared to the bright-field Fourier spectrum.

So since the degeneracy here is a more complex manifestation of the lattice symmetry, we develop an application-specific approach for which we can use physical intuition regarding sample morphology to reduce computational cost.^38,49 Future work adapting more general phase unwrapping approaches to this application may be promising. We also note that since we only care for the derivative of the displacement, it may be possible in principle to avoid unwrapping altogether by modifying the finite difference stencil as done in essentially non-oscillatory schemes.^95,96 We found this strategy sensitive in practice to noise.

The liberal use of integer quadratic programming is prohibitively expensive; therefore, we considerably pre-process the data so that it needs only to be applied in small regions. We first pre-processed the data using physical intuition regarding its geometry, partitioning the displacements into zones of constant lattice vector offsets $( n r, m r)$⁠, such that $n x y, m x y= n r, m r$ for all pixels $(x,y)$ in region r.⁴⁹ 

One partitioning approach we used involves the watershed algorithm.⁹⁷ In this method, several regions are specified [in this case the centers of the $| u x y|≈0$ (AA-type) stacking regions] and pixels of increasing $| u x y|$ are successively assigned. We denote this as watershed-AA. The watershed-AA approach proved effective on the data sets with large, dispersed regions of $| u x y|≈0$ stacking separated by thin boundaries, as was seen for anti-parallel aligned twisted transition metal dichalcogenides.⁴⁹ We also investigate an alternative approach based on a Voronoi partition,⁹⁸ where region boundaries are drawn halfway between the AA locations. This method was successful for highly noisy data but failed when the moiré wavelength varied even marginally over the field of view.⁴⁹ Last, we investigate an approach where the watershed algorithm is applied to the angle of $u x y$ instead of its magnitude to select Bernal-type stacking regions for structures that minimize $| u x y|≈0$ like large-angle twisted graphene and parallel aligned twisted transition metal dichalcogenides. We denote this new approach as watershed-AB, which was not used in our prior applications,⁴⁹ but is conceptually similar.³⁸ After use of any of these three segmentation approaches, the values $n r, m r$ are determined based on the connectivity of neighboring regions and $s x y$ is chosen to maximize a local curl. A quadratic integer program can then be applied to the pre-processed data.⁴⁹ 

We compare these three segmentation approaches on artificial $u x y unwrap$ data corresponding to a rigidly twisted structure with Gaussian-distributed noise to model the effects of sample defects and contamination, which is then converted into $u x y$ following Eq. 4. Figure 4(b) shows the fraction of incorrectly assigned pixels after each segmentation approach. These results illustrate that the Watershed-AA approach performs best for the rigidly twisted domains, but its performance declines more rapidly with increasing noise than the Voronoi approach. The Watershed-AB approach is more variable to subtle details in the applied noise, which can impact the angle of $u x y$ more easily than its magnitude. For applications to experimental data, the most appropriate choice will depend heavily on the morphology of the moiré pattern.

Following the unwrapping procedure, $u x y unwrap$ can then be differentiated numerically, for which we use a centered finite difference stencil. We note that the measured displacements are cell-averaged quantities, which some other stencil choices would need to account for and a factor of one half is introduced as we consider the intralayer strain $∇ u x y top$⁠. In practice, there is often a rotation between the real-space coordinate system and that used to define $u x y$⁠, which need be accounted for.^38,49

Using the techniques described here, we can directly measure the moiré-scale strain in bilayer systems. Experimental strain distributions are shown in Fig. 5 for MoS $2$ with a twist angle of $1.77°$ using a series of coordinate choices. The sheer strain, curl, and divergence maps are shown in Figs. 5(c)–5(e), expected to be uniformly 0, $θ t$⁠, and 0 in a rigid structure, respectively. The clear local structure in 5(c) and 5(d) are consequences of atomic reconstruction, showing that MoS $2$ distorts to increase the local twist angle within the $| u x y|≈0$ [shown black in Fig. 5(a)], reducing the size of these high-energy stacking regions. Simultaneously, the distortion decreases the local twist angle within the low-energy stacking regions [Fig. 5(d)]. These two rotational effects act in a manner that accumulates shear strain along their boundaries [Fig. 5(c)], while the dilational strain associated with uniform stretching of the lattice remains negligible [Fig. 5(e)].⁴⁹ 

We have outlined and analyzed various numerical and instrumental limitations for the BI technique and presented solutions to generalize its applicability to a wider range of bilayer systems. Future work will aim to expand this method for use in the limit of smaller converged probe sizes and more general samples. This may involve any of the discussed strategies, including the use of more precisely defined non-spherical virtual apertures, to enable smaller CBED probes, fitting disk locations in tandem with their intensities, more involved post-processing techniques along the lines of those discussed for interpreting multi-layer overlap regions, and improvements to the algorithms used.

We acknowledge helpful conversations with C. Groshner and formative contributions from N. Kazimerzak, K. Bustillo, H. Brown, J. Ciston, T. Taniguchi, and K. Watanabe in the referenced initial realization of this approach. This material is based upon work supported by the U.S. National Science Foundation Early Career Development Program (CAREER), under Award No. 2238196 (D.K.B). I.M.C. acknowledges a pre-doctoral fellowship award under Contract No. FA9550-21-F-0003 through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program, sponsored by the Air Force Research Laboratory (AFRL), the Office of Naval Research (ONR) and the Army Research Office (ARO). Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

Isaac M. Craig: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). Madeline Van Winkle: Investigation (supporting); Methodology (equal); Writing – review & editing (supporting). Colin Ophus: Conceptualization (supporting); Methodology (supporting); Writing – review & editing (supporting). D. Kwabena Bediako: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Resources (lead); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available within the article. A Python implementation of the discussed approach is available at github.com/bediakolab/pyInterferometry, which imports some functionality from py4DSTEM⁹⁹ among other open-source modules, such as GEKKO.¹⁰⁰