key: cord-1019191-25qex9sj
authors: Jong, Tobias A. de; Benschop, Tjerk; Chen, Xingchen; Krasovskii, Eugene E.; Dood, Michiel J.A. de; Tromp, Rudolf M.; Allan, Milan P.; Molen, Sense Jan van der
title: Imaging moir'e deformation and dynamics in twisted bilayer graphene
date: 2021-07-30
journal: Nat Commun
DOI: 10.1038/s41467-021-27646-1
sha: 3825bac1149a5449a7d4def57ad46f7da0477e7b
doc_id: 1019191
cord_uid: 25qex9sj

In twisted bilayer graphene (TBG) a moir'e pattern forms that introduces a new length scale to the material. At the 'magic' twist angle of 1.1{deg}, this causes a flat band to form, yielding emergent properties such as correlated insulator behavior and superconductivity [1-4]. In general, the moir'e structure in TBG varies spatially, influencing the local electronic properties [5-9] and hence the outcome of macroscopic charge transport experiments. In particular, to understand the wide variety observed in the phase diagrams and critical temperatures, a more detailed understanding of the local moir'e variation is needed [10]. Here, we study spatial and temporal variations of the moir'e pattern in TBG using aberration-corrected Low Energy Electron Microscopy (AC-LEEM) [11,12]. The spatial variation we find is lower than reported previously. At 500{deg}C, we observe thermal fluctuations of the moir'e lattice, corresponding to collective atomic displacements of less than 70pm on a time scale of seconds [13], homogenizing the sample. Despite previous concerns, no untwisting of the layers is found, even at temperatures as high as 600{deg}C [14,15]. From these observations, we conclude that thermal annealing can be used to decrease the local disorder in TBG samples. Finally, we report the existence of individual edge dislocations in the atomic and moir'e lattice. These topological defects break translation symmetry and are anticipated to exhibit unique local electronic properties.

In twisted bilayer graphene (TBG) a moiré pattern forms that introduces a new length scale to the material. At the 'magic' twist angle θ m ≈ 1.1 • , this causes a flat band to form, yielding emergent properties such as correlated insulator behavior and superconductivity [1] [2] [3] [4] . In general, the moiré structure in TBG varies spatially, influencing the local electronic properties [5] [6] [7] [8] [9] and hence the outcome of macroscopic charge transport experiments. In particular, to understand the wide variety observed in the phase diagrams and critical temperatures, a more detailed understanding of the local moiré variation is needed [10] . Here, we study spatial and temporal variations of the moiré pattern in TBG using aberration-corrected Low Energy Electron Microscopy (AC-LEEM) [11, 12] . The spatial variation we find is lower than reported previously. At 500 • C, we observe thermal fluctuations of the moiré lattice, corresponding to collective atomic displacements of less than 70 pm on a time scale of seconds [13] , homogenizing the sample. Despite previous concerns, no untwisting of the layers is found, even at temperatures as high as 600 • C [14, 15] . From these observations, we conclude that thermal annealing can be used to decrease the local disorder in TBG samples. Finally, we report the existence of individual edge dislocations in the atomic and moiré lattice. These topological defects break translation symmetry and are anticipated to exhibit unique local electronic properties.

In charge transport experiments, a percolative average of the microscopic properties is measured. Therefore, any local variation in twist angle and strain in TBG will influence the result of such experiments. However, imaging such microscopic variations is non-trivial. A myriad of experimental techniques has been applied to the problem [16] [17] [18] [19] [20] , each only resolving part of the puzzle due to practical limitations (capping layer or device substrate, surface quality or measurement speed).

Here, we use (AC-)LEEM, which measures an image of the reflection of a micron-sized beam of electrons at a landing energy E 0 (0-100 eV, referenced to the vacuum energy) in real space, in reciprocal space (diffraction), or combinations thereof. This allows us to perform large-scale, fast and non-destructive imaging of TBG. Additionally, spectroscopic measurements, yielding information on the material's unoccupied bands can be done by varying E 0 [21, 22] : Device-scale imaging of TBG. a, Local spectra used to determine the graphene layer count. b, Stitched composite bright field overview of a sample using 4 eV, 8 eV and 17 eV as imaging energies in red, green and blue respectively (see main text for color interpretation). Visible defects include folds in black, tears, where the monolayer or even bare hBN shines through, bubbles (bright) and some polymer residue in the lower and upper right (dark speckles). c, Ab-initio calculations of LEEM spectra of different relative stackings of bilayer graphene, 37 eV indicated. d, Stitched bright field overview of the same sample imaged at E0 = 37 eV, for optimal stacking contrast, revealing the moiré patterns. e-h, Crops of different twist angles areas from d, inset shows Fourier transforms and the detected moiré peaks. All measurements have been done on the same sample, featured in b.

In Figure 1a , LEEM spectra are shown, taken on several locations of a TBG sample. These LEEM spectra are directly related to layer count, as described in refs. [22] [23] [24] [25] ; on the one hand via interlayer resonances in the 0-5 eV range, on the other hand via the gradual disappearance of a minimum at 8 eV. Here, more graphene layers (having a band gap at 8 eV) are progressively masking an hBN band underneath. This allows us to determine the local graphene layer count for each point on our sample. To visualize that, we choose three characteristic energies, i.e. E 0 = 4 eV (red), E 0 = 8 eV (green) and E 0 = 17 eV (blue) (see Figure 1a) , and combine stitched overviews at these energies into a single false-color image (Figure 1b) . This overview confirms that the sample consists of large TBG areas (darker green in Figure 1b ) surrounded by monolayer graphene (pink), on an hBN flake (blue/purple) on silicon (black). Stripes of brighter green indicate areas of 2-on-2 graphene layers (lower stripe), 2-on-1 (upper stripe), and 1-on-2 (wedge on the lower right). The relatively homogeneous areas are themselves separated by folds, appearing as black lines. The folds locally combine in larger dark nodes (confirmed by AFM, see Supplementary Information A). A few folds, however, have folded over and appear as lines of higher layer count. Hence, Figure 1b provides a remarkable overview of a larger-scale sample, with detailed local information. Increasing E 0 beyond 25 eV, stacking boundaries and AA-sites become visible [23, 26] . This is consistent with ab-initio calculations of LEEM spectra for different relative stackings, as presented in Figure 1c . Therefore, imaging at E 0 = 37 eV (indicated in Figure 1c ) yields a precise map of the moiré lattice over the full TBG area (see Figure 1d ). We find that separate areas, between folds, exhibit different moiré periodicities and distortion [27] . This allows us to study different moiré structures on a single sample. Figure 1e -h shows full resolution crops of areas indicated in Figure 1d . The observed twist angles on this sample range from < 0.1 • to 0.7 • . For smaller angles, we observe local reconstruction towards Bernal stacking within the moiré lattice, consistent with literature [17, 20] . The best resolution was reached on another sample with a twist angle of 1.3 • (See Supplementary Figure 11 ).

The moiré patterns shows distortions, corresponding to local variations in twist angle and (interlayer) strain. Near folds, for instance, the strain increases resulting in strongly elongated triangles, for example in the lower right corner of Figure 1d [28] . Despite their relative homogeneity, the moiré areas in Figure 1e -h also show subtle distortions. As structural variations correlate directly with local electronic properties, we quantify them in detail [5, 6] . For this, we use adaptive geometric phase analysis (GPA), extending our earlier work on STM data of moiré patterns in TBG (see Supplementary Information) [29] [30] [31] [32] [33] . This method, illustrated in Figure 2a -c, yields the displacement field with respect to a perfect lattice, by multiplying the original image with complex reference waves followed by low-pass filtering to obtain the GPA phase differences, which are then converted to the displacement field. This field fully describes the distortion of the moiré lattice and allows us to extract key parameters such as the local twist angle θ * ( r) (see Figure 2d ), and heterostrain magnitude ( r) and direction (see Figure 2e ). [31, 34] The distortions of the moiré pattern correspond directly to distortions of the atomic lattices, magnified by a factor 1/θ and rotated by 90 • + θ/2 [31, 35] .

The extracted variation in twist angle and heterostrain for various regions of the sample, including those in Figure 1e -h, is summarized in Figure 2f ,g, respectively. The twist angle variation within each domain is much smaller than the variation in twist angle between the separate, foldbounded areas. Within domains, standard deviations range from 0.005 • to 0.015 • , i.e. significantly smaller (by a factor 3-10) than previously reported [18, 20, 31] . The strain observed is around a few tenths of a percent, which is considerable. In some domains, we find an average strain of the atomic lattice of up to = 0.4 %. According to earlier theoretical work, such values are high enough to locally induce a quantum phase transition [8] .

The variation of , as for θ * , within the domains is significantly lower than in earlier studies. We do note that the use of GPA introduces a point spread function (PSF) that is broader than the PSF of the instrument, resulting in a lower displacement field frequency response at small scales and therefore a somewhat reduced variation. Nevertheless, the combined PSF of instrument and analysis is still comparable to other techniques that do not image the unit cell directly, allowing for a direct comparison.

We hypothesize that the difference in variations with literature stems from the relatively high temperature at which we annealed and measured the sample, combined with the relatively long averaging time of this measurement (≥ 16 s for all data in Figure 2 ). The high temperature induces thermal fluctuations of the lattice (as demonstrated below), allowing the system to approach a lower energy, more homogeneous, state.

So far, we have discussed structural properties varying on the moiré length scale. However, the moiré magnification of deformations is general and extends to atomic edge dislocations (visualized in Figure 3a ). This type of topological defect stems from a missing row of atomic unit cells and is characterized by an in-plane Burgers vector (in red) [36, 37] . The addition of a second (twisted) atomic layer magnifies (and rotates) the defect to an edge dislocation in the moiré lattice (illustrated in Figure 3d ,e) [35] . In all cases, the defect can be characterized using GPA, pinpointing its location and Burgers vector (Figure 3b ,e) [38] .

In our sample, we find a few such defects in the moiré lattice (see SI). In Figure 3f ,g we show an edge dislocation in a topographically flat region with θ = 0.63 • (see SI). Contrary to TEM observations on single-layer graphene [38] , we do not observe creation or annihilation of edge dislocation pairs in the microscope, even at elevated measurement temperatures (500 • C) and under low-energy electron irradiation. Furthermore, the mobility of the defects is low. One edge dislocation did move, over several moiré cells between measurements, after which it remained at the same position even after a month at room temperature and reheating (see supplementary information). This stability suggests that the moiré lattice itself plays a role in stabilizing these defects, via a minimum of the local stacking fault energy within the moiré unit cell.

These topological dislocations break translational symmetry, which may lead to singular electronic properties on the local scale [39] [40] [41] . Specifically, a phase difference will appear between electron paths encircling the defect clockwise and counterclockwise.

All measurements presented so far were performed at 500 • C, to minimize hydrocarbon contamination under the electron beam. In literature, there is concern about the graphene layers untwisting at such temperatures, due to energy differences between different rotations [14, 15] . We, however, see no sign of that. The twist angles within the domains are stable from 100 • C up to 600 • C for all samples studied. However, we did observe a more subtle thermal influence on the moiré pattern. (Figure 4d,e) . Moreover, we can quantify these fluctuations via the difference in displacement field with respect to the image at t = 0 s using GPA (Figure 4f,g) . Interestingly, these involve the collective movement of millions of atoms, but only over very small distances.

We stress that a translation of the domain boundary by 4 nm, as observed, corresponds to a shift of less than half the width of the domain boundary itself [17, 42] . As the relative shift of the layers over the full domain boundary is a single carbon bond length, the corresponding atomic translations are less than half of that, i.e. less than 70 pm. Hence, the 'moiré magnification' makes it possible to detect these sub-Angstrom changes in TBG in real time using LEEM. Our data suggest that domain boundary displacement follows a random pattern of forward and backward steps. This indicates a possible source for the twist angle disorder observed in low(er) temperature experiments [10, 18, 20, 31] : frozen-in thermal fluctuations of the moiré lattice. The thermal fluctuations found, corresponding to ±0.005 • for twist angle and ±0.02% for strain, are smaller than the extracted static deformations, though not negligible. Note that these values are damped by the intrinsic broadening of GPA and the time integration. Future experiments will focus on deducing the detailed statistics of the domain boundary dynamics versus temperature. Following these local collective excitations in time, will yield quantitative information on the energy landscape of these atomic lattice deformations within the moiré lattice. This will be important to answer the question if moiré lattices can be relaxed and homogenized using controlled annealing. If so, this would yield higher-quality magic-angle TBG devices in which charge transport is not limited by percolative effects and higher critical temperatures are reached.

Our quantitative LEEM study on TBG reveals a wide variety in twist angles and strain levels in a single sample. We show that spontaneous changes in global twist angle do not occur, even during high-temperature annealing, but that local collective fluctuations do take place. This suggests that high-temperature annealing causes relaxation of the local moiré lattice, reducing lattice disorder. Vice versa, this points to frozen-in thermal fluctuations as a possible source of the (significant) short-range twist angle disorder observed previously. Furthermore, this potentially offers insight into energetic aspects of the atomic lattice deformation within the moiré lattice.

We also report the observation of stable topological defects, i.e. edge dislocations, in the moiré lattice of two Van der Waals layers. Combining our methods with other techniques that can access the electronic structure, such as STS, nanoARPES, and even in-situ potentiometry [43] , will allow for a systematic study of the electronic properties around these defects. Finally, the methods we describe here extend beyond TBG, to any type of twisted system. Therefore, our work introduces a new way of studying deformations of moiré patterns and of connecting these to the (local) electronic properties of this exciting class of materials.

The twisted bilayer graphene sample was fabricated using the standard tear-and-stack method [1, 44] . The monolayer graphene was first exfoliated with scotch tape on to a SiO 2 /Si substrate. A polycarbonate (PC)/polydimethylsiloxane (PDMS) stamp was used for the transfer process, where the PC covered only half of the PDMS surface. After the first half of the graphene flake was successfully torn and picked up, it was rotated by 1.0 • . The flake was then overlapped with the bottom half and used to pick it up. The stack was then stamped on a moderately thick (∼140nm) hBN flake, priorly exfoliated with PDMS on to a silicon substrate, along with the PC layer. Part of the graphene flake is deliberately put in contact with the Silicon surface for electrical contact purposes, i.e. to absorb the beam current. The whole substrate is then left in chloroform for 3 hours to dissolve the PC. All flakes were exfoliated from crystals, commercially bought from HQ Graphene and the fabrication process was performed using the Manual 2D Heterostructure Transfer System sold by the same company.

All LEEM measurements where performed in the ESCHER LEEM, based on the SPECS P90 [11, 12, 45] . Samples were loaded into the ultrahigh-vacuum (base pressure better than 1.0×10 −9 mbar) LEEM main chamber and slowly heated to 500 • C (as measured by pyrometer and confirmed by IR-camera) and left to anneal overnight to get rid of any (polymer) residue. All measurements were conducted at elevated temperatures, 450 • C -500 • C, unless specified otherwise. The sample was located on the substrate using photo-emission electron microscopy (PEEM) with an unfiltered mercury short-arc lamp, by comparing to optical microscopy images taken beforehand. Spectra were taken in high-dynamic-range mode and drift corrected and all images were corrected for detector artefacts, as described in ref. [46] . When needed to obtain a sufficient signal-to-noise ratio, multiple 250 ms exposures were accumulated for each image, e.g. 8 exposures (2 seconds) per landing energy in the spectra in Figure 1b and 16 exposures (4 seconds) at each location for the overview at 37 eV in Figure 1c ,e.

To measure the dynamics as presented in Figure 4 , a time series of accumulated 4 × 250 ms = 1 s exposure images was taken back-to-back. After regular detector correction and drift correction, each image was divided by a Gaussian smoothed version of itself (σ = 50 pixels) to get rid of spatial and temporal fluctuations in electron illumination intensity. To further reduce noise, a Gaussian filter with a width of σ = 1 image ∼ 1 s, was applied in the time direction before applying GPA.

To enable high resolution, large field-of-view LEEM imaging, the LEEM sample stage [47] was scanned in a rectangular pattern over the sample, taking an image at each position, leaving sufficient overlap (2 µm steps at a 4.7 µm field-of-view). To obtain meaningful deformation information from this, care needs to be taken to use a stitching algorithm that does not introduce additional deformation, i.e. as faithfully reproducing reality as the constituting images.

To achieve this, a custom stitching algorithm tailored towards such LEEM data, was developed, as described in the Supplementary information and in the implementation [48] .

In addition, for the composite bright field in Figure 1a , minor rotation and magnifications differences due to objective lens focus differences were compensated for. This was done by registering the stitches for different energies using a log-polar transformation based method to obtain relative rotations and magnification. Subsequently, areas where a color channel was missing, were imputed using a k-nearest neighbor lookup in a regularly sliced subset of the area with all color channels present.

To quantify the large deviations in lattice shape due to the moiré magnification of small lattice distortions, we extended the GPA algorithm to use an adaptive grid of reference wave vectors, based on related to earlier work in laser fringe analysis [32] .

The spatial lock-in signal is calculated for a grid of wave vectors around a base reference vector, converting the GPA phase to reference the base reference vector every time. For each pixel, the spatial lock-in signal with the highest amplitude is selected as the final signal. To avoid the problem of globally consistent phase unwrapping, the gradient of each GPA phase was directly converted to the displacement gradient tensor. More details of the used algorithm are given in the Supplementary information.

All image analysis code was written in Python, using Numpy [49] , Scipy [50] , scikit-image [51] and Dask [52] . The core algorithms will be made available as an open source Python package [33] . Throughout the development of the algorithms and writing of the paper, matplotlib [53, 54] was extensively used for plotting and figure creation.

The theoretical reflectivity spectra are obtained with the ab-initio Bloch-wave-based scattering method described in ref. [55] . Details of the application of this method to stand-alone twodimensional films of finite thickness can be found in ref. [56] . The underlying all-electron Kohn-Sham potential was obtained with a full-potential linear augmented plane-wave method within the local density approximation, as explained in ref. [57] . Inelastic scattering is taken into account by an absorbing imaginary potential −iV i , which is taken to be spatially constant (V i = 0.5 eV) over a finite slab (where the electron density is non-negligible) and to be zero in the two semiinfinite vacuum half spaces. In addition, a Gaussian broadening of 1 eV is applied to account for experimental losses.

We thank Marcel Hesselberth and Douwe Scholma for their indispensable technical support. We thank D. K. Efetov and J. Aarts for getting us started on twisted bilayer graphene. We would like to thank J. Jobst and D.N.L. Kok for scientific discussions and L. Visser and G. Stam in particular for their input on the stitching algorithm. We thank Federica Galli for both scientific discussion as well as technical support. This work was supported by the Dutch Research Council (NWO) as part of the Frontiers of Nanoscience programme. It was also supported by the Spanish Ministry of Science, Innovation and Universities (Project No. PID2019-105488GB-I00).

TAdJ performed the LEEM measurements and wrote the LEEM data analysis code. XC fabricated the sample and performed AFM measurements. TB contributed to the analysis. EK performed the reflectivity calculations. TAdJ and SJvdM took the initiative for the writing of the manuscript, with contributions from several others. All authors contributed to the scientific discussion. To further characterize the surface properties of the sample, an AFM (JPK, NanoWizard 3) measurement was performed in AC tapping mode following the LEEM measurements. Predominantly, the results show a very flat and clean graphene surface between folds, indicating annealing at 500 • C in UHV had successfully removed the polymer residue left on the surface.

In profile 1, the terrace height sees a difference of 0.3 nm, demonstrating the graphene layer count goes down by one at this location. This corresponds to the layer counts extracted from the LEEM spectra.

Profile 2 to 4 shows 3 different kinds of defects in the bilayer graphene region. The ridge at location 2 seems to be a neat folding of both the bilayer graphene flake (1.5 nm in height, four layers of graphene), whereas profile 4 shows wrinkles that are up to 120 nm tall. This is also reflected by the distinct patterns in the LEEM bright field overview image, respectively. While the wrinkles merely appear black, the ridge resembles more like a unique layer count domain. Profile 3 shows two tears within the one layer of graphene, corresponding nicely to the defect region observed in LEEM where monolayer graphene shines through.

The zoomed-in small scale measurements marked by the red and yellow box shows the topography on top of two dislocations observed in LEEM. As shown in Figure 6 , no distinctive feature was observed at either dislocation. The topography, however, shows an exceptionally flat surface with a height variation (peak-to-peak) of less than 1 nm.

Regular GPA is limited in the wave vector deviations (with respect to the reference wave vector) it can measure, due to the limitations in spectral leakage. This is no problem when applied to atomic lattices, as the expected deviations are very small there. However, due to the moiré magnification of small lattice distortions, it does become a limiting factor when applying GPA to small twist angle moiré lattices.

To overcome this limitation, we extended the GPA algorithm to use adaptive reference wave vectors, based on the combination of two ideas and related to earlier work in laser fringe analysis [32] : First, a GPA phase calculated with respect to one reference vector can always be converted to the GPA phase with respect to another reference vector by adding a phase corresponding to the phase difference between the reference vectors. Second, a larger lock-in amplitude corresponds to a better fit between the reference vector and the data.

The adaptive GPA algorithm therefore works as follows: The spatial lock-in signal is calculated for a grid of wave vectors around a base reference vector, converting the GPA phase to reference the base reference vector every time. For each pixel, the spatial lock-in signal with the highest amplitude is selected as the final signal.

It was realized that to deduce the deformation properties, reconstruction to a globally consistent phase (requiring 2D phase unwrapping), as reported previously [31] , is not strictly needed, making it possible to circumvent the problems associated with 2D phase unwrapping. Instead, the gradient of each GPA phase was calculated, requiring only local 1D phase unwrapping (i.e. assuming the derivative of the phase in both the x and y direction will never be more than π per pixel, an assumption in practice always met). Subsequently, these three GPA gradients are converted to the displacement gradient tensor (in real space coordinates), estimating the transformation via weighted least squares, using the local spatial lock-in amplitudes as weights.

As an added benefit, this entire procedure is local, i.e. not depending on pixels beyond nearest neighbors in any way except for the initial Gaussian convolution in determining the GPA. This reduces the effect of artefacts in the image to a minimum local area around each artefact (where for in the 2D phase unwrapping they have a global influence on the phases).

However, when the gradient is computed based on phase values stemming from two different GPA reference vectors, i.e. at the edge of their valid/optimal regions, artefacts appear due to their relatively large absolute error. To prevent this, the local gradient of the phase with the highest lock-in magnitude is stored alongside the lock-in signal itself in the GPA algorithm. This way, the gradient is calculated based on a single reference phase, propagating only the much smaller relative/derivative error between the two signals instead of the absolute error.

As mentioned in the main text, even adaptive GPA has its limits. In particular, too large deviations from the base reference vector can not be resolved correctly, causing an erroneous, lower, extracted deviation, as is visible in the lower right of Figure 2e in the main text). As the deformation becomes too large, e.g. towards the folds in the TBG, the highest lock-in amplitude will occur at a different moiré peak or at the near-zero components of the fourier transform, causing an incorrect value to be extracted.

To achieve stitching of images without inducing any additional deformation, a custom stitching algorithm tailored towards such LEEM data, was developed, working as follows:

To compensate sample stage inaccuracy, nearest neighbor (by sample stage coordinates) images are compared, finding their relative positions by cross-correlation. Using an iterative procedure, calculating cross correlations of overlapping areas at each step, the absolute positions of all images are found. Images are then combined in a weighted fashion, with the weight sloping to zero at the edges of each image, to smooth out any mismatch due to residual image warping. The full stitching algorithm is implemented in Python, available as a Jupyter Notebook [48] .

the valid area of the image defined by mask/radius, take an area as close to the intended area but still within the valid area as possible.

ii. Find the location of the maximum in the cross-correlation. This corresponds to the correction to the estimate of the difference vector between the corresponding image position pair. iii. Calculate the weight of the match by dividing the maximum in the cross-correlation by the square root of the maximum of the auto-correlations.

3. Compute a new estimate of the difference vectors by adding the found corrections.

Reconvert to a new estimate of pixel coordinates by minimizing the squared error in the system of equations for the positions, weighing by modified weights, either:

i. w mod = w − w min for w > w min , w = 0 else, with w min the maximum lower bound such that the graph of nearest neighbours with non-zero weights is still connected ii. Only use the 'maximum spanning tree' of weights, i.e. minus the minimum spanning tree of minus the weights, such that only the n best matches are used.

6. (Optional) Refine the estimate of the transformation matrix, using all estimated difference vectors with a weight better than w minest and restart from step 3.

7. Repeat step 4. and 5. until sigma is satisfactory small. Optional repeat a final time with the original data if the signal to noise of the original data permits.

8. Select only the images for stitching where the average of the used weights (i.e. where w > w min ) is larger than q thresh for an appropriate value of q thresh .

9. (Optional) For those images, match the intensities by calculating the intensity ratios between the overlap areas of size fftsize*fftsize and perform a global optimization. 10 . Define a weighting mask, 1 in the center and sloping linearly to zero at the edges of the valid region, over a width of bandwidth pixels, as illustrated in Figure 8 .

11. Per block of output blocksize*blocksize, select all images that have overlap with the particular output block, multiply each by the weighting mask and shift each image appropriately. Divide by an equivalently shifted stack of weighting masks. As such information at the center of images gets prioritized, and transitions get smoothed.

For square grids with a decent amount of overlap, it makes sense to put n neighbors to 5 (including the image itself), however, for larger overlaps or datasets where an extra dimension is available (such as landing energy), it can be appropriate to increase the number of nearest neighbors to which each image is matched. Parameters and intermediate results of the iteration are saved in an xarray and saved to disk for reproducibility.

Using dask, the following steps are parallelized:

• step 5B , where each pair of images can be treated independently. In practice parallelization is performed over blocks of subsequent images with their nearest neighbours. This could be improved upon in two ways: firstly by treating each pair only once, and secondly by making a nicer selection of blocks of images located close to each other in the nearest neighbor graph. This would most likely require another (smarter) data structure than the nearest neighbour indexing matrix used now.

• Step 6 is quite analogous to 5B and is parallelized similarly. • Step 11 is parallelized on a per block basis. To optimize memory usage, results are directly streamed to a zarr array on disk.

• The minimizations are parallelized by scipy natively. 

Full movie showing larger Field of View compared to Figure 4 , in real space data, difference data and GPA-extracted displacement field.

Unconventional superconductivity in magic-angle graphene superlattices

Moire bands in twisted double-layer graphene

Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene

Observation of flat bands in twisted bilayer graphene

Twisted graphene bilayer around the first magic angle engineered by heterostrain

Strain impacts on commensurate bilayer graphene superlattices: Distorted trigonal warping, emergence of bandgap and direct-indirect bandgap transition". en

Designing flat bands by strain

Strain-Induced Quantum Phase Transitions in Magic-Angle Graphene

Flat band carrier confinement in magic-angle twisted bilayer graphene

Graphene bilayers with a twist

A new aberration-corrected, energy-filtered LEEM/PEEM instrument. I. Principles and design

A new aberration-corrected, energy-filtered LEEM/PEEM instrument II. Operation and results

Moiré phonons in twisted bilayer graphene

Electronic correlations in twisted bilayer graphene near the magic angle

Correlated insulator behaviour at half-filling in magic-angle graphene superlattices

Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene

Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene

Mapping the twist-angle disorder and Landau levels in magic-angle graphene

Moiré metrology of energy landscapes in van der Waals heterostructures

Strain fields in twisted bilayer graphene". en

Intensity-voltage low-energy electron microscopy for functional materials characterization

Quantifying electronic band interactions in van der Waals materials using angle-resolved reflected-electron spectroscopy

Growth and electronic transport properties of epitaxial graphene on SiC

Low-energy electron reflectivity from graphene

Measuring the Local Twist Angle and Layer Arrangement in Van der Waals Heterostructures

Intrinsic stacking domains in graphene on silicon carbide: A pathway for intercalation

Rotational Disorder in Twisted Bilayer Graphene

Evidence of flat bands and correlated states in buckled graphene superlattices". en

Quantitative measurement of displacement and strain fields from HREM micrographs

Intra-unit-cell electronic nematicity of the high-Tc copper-oxide pseudogap states

Measuring local moiré lattice heterogeneity of twisted bilayer graphene

Two-dimensional windowed Fourier transform for fringe pattern analysis: Principles, applications and implementations

Maximized electron interactions at the magic angle in twisted bilayer graphene

Moiré pattern as a magnifying glass for strain and dislocations in van der Waals heterostructures". en

Geometrical considerations concerning the structural irregularities to be assumed in a crystal

Mathematical Theory of Dislocations and Fracture

Dislocation-Driven Deformations in Graphene

Berry phase of dislocations in graphene and valley conserving decoherence

Direct observation of magneto-electric Aharonov-Bohm effect in moir\'escale quantum paths of minimally twisted bilayer graphene

Aharonov-Bohm Oscillations in Minimally Twisted Bilayer Graphene

Structural and electronic transformation in lowangle twisted bilayer graphene". en

Low-Energy Electron Potentiometry: Contactless Imaging of Charge Transport on the

Tunable moiré bands and strong correlations in small-twist-angle bilayer graphene

Low-energy electron microscopy and spectroscopy with ESCHER: Status and prospects

Quantitative analysis of spectroscopic Low Energy Electron Microscopy data: High-dynamic range imaging, drift correction and cluster analysis

A versatile ultra high vacuum sample stage with six degrees of freedom

TAdeJong/LEEM-analysis: LEEM-analysis

Array programming with NumPy

SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python

scikit-image: image processing in Python

Dask Development Team. Dask: Library for dynamic task scheduling

Matplotlib: A 2D graphics environment

matplotlib/matplotlib: REL: v3.3.2. 2020

Augmented-plane-wave approach to scattering of Bloch electrons by an interface

Ab Initio Theory of Photoemission from Graphene

Augmented Fourier components method for constructing the crystal potential in self-consistent band-structure calculations

Numba

It is designed for use with ESCHER LEEM images. For those images, their positions are known approximately in terms of stage coordinates, i.e. the positions as reported by the sensors in the sample stage. It should however generalize to any set of overlapping images where relative positions of the images are known in some coordinate system which can approximately be transformed to coordinates in terms of pixels by an affine transformation (rotation, translation, mirroring).The algorithm consists of the following steps:1. Using the stage coordinates for each image, obtain a nearest neighbour graph with the nearest n neighbors neighbouring images for each image.2. Obtain an initial guess for the transformation matrix between stage coordinates and pixel coordinates, by one of the following options:1. Copying a known transformation matrix from an earlier run of a comparable dataset.2. Manually overlaying some nearest neighbor images from the center of the dataset, either refining the estimate, or making a new estimate for an unknown dataset 3. Calculate an initial estimate of the pixel coordinates of the images by applying the corresponding transformation to the stage coordinates 4. Apply a gaussian filter with width sigma to the original dataset and apply a magnitude sobel filter. Optionally scale down the images by an integer factor z in both directions to be able to reduce fftsize by the same factor, without reducing the sample area compared.5. Iterate the following steps until the calculated image positions have converged to within sigma:1. Obtain a nearest neighbour graph with per image the nearest n neighbors neighbouring images from the current estimate of the pixel coordinates and calculate the difference vectors between each pair of nearest neighbours.2. For each pair of neighboring images:i. Calculate the cross-correlation between areas estimated to be in the center of the overlap of size fftsize*fftsize of the filtered data. If the estimated area is outside