{'fraction_non_alphanumeric': 0.037418180545346354, 'fraction_numerical': 0.009240785466764675, 'mean_word_length': 4.615172955974843, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mapping local heterogeneity in open twisted bilayer graphene devices

We introduce a new method to continuously map inhomogeneities of a moir\'e lattice and apply it to open-device twisted bilayer graphene (TBG). We show that the variation in the twist angle, which is frequently conjectured to be the reason for differences between devices with a supposed similar twist angle, is about 0.04{\deg} over areas of several hundred nm, comparable to devices encapsulated between hBN slabs. We distinguish between an effective twist angle and local anisotropy and relate the latter to heterostrain. Our results suggest that the lack of evidence for superconductivity in open devices is not a consequence of higher heterogeneity in the twist angle, but possibly due to the absence of interaction with a top hBN layer. Furthermore, our results imply that for our devices, twist angle heterogeneity has a roughly equal effect to the electronic structure as local strain. The method introduced here is applicable to results from different imaging techniques, and on different moir\'e materials.

The kinetic energy of the electrons changes rapidly as the twist angle is varied, especially around the magic angle, therefore the fabrication of devices with just the right angle is key in making them superconducting. But getting the right angle might not even be the most challenging aspect of fabricating high-quality superconducting TBG devices: contaminations, internal stress, and heterogeneities of the twist angle are difficult to avoid. This is in part because the magic angle is not the lowest energy configuration and in part because of the strong forces associated with the tear-andstack technique. Internal stress and heterogeneities are often conjectured to limit the quality of the devices and are attributed as the main causes for the variability between devices (13). This holds especially for open devices that lack the hBN top layer; notably such devices have never been found to superconduct. Measuring, visualizing, and characterizing heterogeneity in the twist angle and strain in TBG is thus crucial to understand and improve devices.
Probably the most complete visualization of inhomogeneity thus far has been obtained using scanning SQUID-on-tip microscopy (SOT) (14). SOT measures the Landau levels as a function of location and thus has access to the density of states at the Fermi level. On encapsulated devices, SOT has been used to visualize heterogeneity on length scales of a few micron with a resolution of several tens of nanometers, demonstrating that the twist angle varies by less than 4% (14). While being a very precise measure of the density of states, SOT cannot differentiate between inhomogeneities of the chemical potential, the twist angle, the local magnetic screening and other factors that influence the density of states. Other techniques to access homogeneity are Nano-ARPES (15)(16)(17), which can image the full electronic structure in reciprocal space with a spatial resolution of circa 600 nm, low energy electron microscopy (LEEM) (15), which can image structural inhomogeneities at twist angles lower and higher than the magic angle, conductive atomic-force microscopy (AFM) (18), nano-photocurrent mapping (19), which can measure the twist angle with a resolution on the order of ~20 nm, and scanning single electron transistors (20), which can map the twist angle by measuring the inverse local compressibility. Finally, scanning tunneling microscopy (STM), the probe used in this study, has been used to measure both the topography and the local density of states of TBG, including the emergence of correlations at the magic angle (21)(22)(23)(24)(25)(26). In previous STM studies, two different methods have been used to determine the local twist angle. First, one can determine the twist angle using three neighboring moiré lattice sites in real space.
The distances between each lattice site, λ1 , λ2 , λ3 are fit to a set of equations that yield the twist angle at a per-triangle resolution (figure 1a) (22), and, using a model with assumptions about the strain distribution in the two layers, an estimate for the heterostrain ε.
A second method to determine the twist angle, uses the Fourier transform of a real space topography to determine the moiré wavelengths λj in the three directions of the moiré lattice (in principle, two directions are fully determining the lattice, but often all three are used for a better signal-to-noise ratio). The twist angle is determined using = is the lattice constant of graphene. Using the Fourier transform is generally more accurate than fitting three moiré lattice peaks, because it averages over the whole field of view, but this also limits its spatial resolution to the full field of view.
In this work, we introduce an alternative method of quantifying and visualizing the heterogeneity in open devices, with sub-moiré lattice cell resolution over length scales of hundreds of nanometres.
We develop a spatial lock-in method that enables one to map as a function of spatial location the local twist angle θ*(r), the local moiré anisotropy κ(r), and the anisotropy direction ψ(r), as defined below. We then apply our method to determine the heterogeneity in open TBG devices.
We fabricate our devices using the tear and stack method with a special focus on avoiding contamination to ensure the large clean areas needed for this study. A single graphene flake is precut in halves with an AFM tip, ensuring initial crystallographic alignment between them. The first half is subsequently picked up with a hBN flake, mechanically exfoliated on a SiO2/Si chip and adhered to a PDMS/PC stamp at ~100°C. The second half of graphene is manually rotated to a target twist angle of 1.5° -2.0° and consequently picked up by the hBN/graphene stack on PC. In the next step, the PC layer is carefully peeled off of the initial PDMS stamp and placed on another PDMS stamp up-side down. The sacrificial polycarbonate (PC) layer is then removed in 1-Methyl-2-Pyrrolidone.
Subsequently, the TBG/hBN heterostructure is transferred on a target SiO2/Si substrate with a prepatterned navigation structure, two gold electrodes and a graphite gate contacting one of them within the measurement area. We carefully align the TBG/hBN stack with the local graphite gate to avoid short circuiting. The second pre-patterned gold electrode is used to electrically contact TBG using another graphite piece. The devices are inserted into our ultra-high-vacuum setup and annealed at 350° for 12h before inserting them into the low-temperature STM operating at 4.2K. The TBG samples are located using a capacitive navigation scheme (27).  field-of-views are well suited for spectroscopy studies, we require large field of views that encompass many moiré cells for the heterogeneity study using spatial lock-in detection presented here. Figure 2c shows an example. The general method of spatial lock-in is illustrated in figure 3a for the one-dimensional case: the "measured" signal S(x), a not-quite periodic signal, is multiplied with a reference signal, a perfectly periodic complex plane wave Sref(x). The phase of the resulting signal, when low pass-filtered, corresponds to the local phase of the original wave. To obtain the local variations in wavelength λ(x) of the original wave, shown on the bottom, one calculates the derivative of the local phase. Spatial lock-in algorithms like this have been used previously in electron microscopy studies (known as geometric phase analysis) (28)(29)(30) and optical metrology (31). In the context of STM, the most wellknown application is known as the Lawler-Fujita algorithm (32). Lawler, Fujita et al. have, based on earlier work by Slezak et al. (33), introduced a lock-in algorithm to correct topographic images for drift by calculating the displacement field, i.e. the vectors that connect the coordinates of the measured images with the points of an ideal reference lattice. Our motivation here is different: we do not need to correct an imperfect image, but want to extract heterogeneities of the lattice.

S(x)
To do so, we start with defining three reference plane waves * ( ) = , # ⋅ , ∈ {1,2,3}, where the reference wavevectors qj are determined by simultaneously fitting six gaussians to the Bragg peaks in the Fourier transform of the topography (figure 4b). In order to measure deviations from an isotropic triangular lattice, we force the reference wavevectors to be of equal magnitude and 60 degrees with respect to each other (although see supporting information on choice of reference vectors). The reference lattice is then defined as the real part of the sum of the reference plane waves, i.e. 0 ( ) = The transformation between the measured lattice, Tm(r), and this perfectly periodic, hexagonal reference lattice, Tr(r) can approximately be parametrized as the shifts between points in the moiré lattice and corresponding points in the reference lattice. That is, by approximating our topography, including heterogeneity, by introducing the displacement field, u(r), in the following manner: To extract the displacement field from our data, each reference signal is multiplied with the original topographic image and low-pass filtered with a gaussian window. This operation corresponds to convolution of the original topographic image with the plane wave encompassed by a gaussian, calculating the relevant wave vector component of the ultimately small window 2D Fourier transform.
The window of the gaussian filter needs to be chosen large enough (small enough in frequency space), in order to exclude larger frequencies, but simultaneously small enough (big enough in frequency space) to maintain good spatial resolution (see supporting information). In practice a filter width of a few periods is used, as illustrated by the circle in figure 4a. The local phase of the result of this operation corresponds to the local shift between the real image and the reference wave, or more This local phase is 2π periodic and needs to be phase-unwrapped to remove discontinuities. After phase unwrapping, the displacement field u(r) can be extracted from two of the phase maps by pixelwise multiplication with Q -1 , the inverse of a matrix containing the used wave vectors (Although not applied here, using all three wave vectors is more involved but can be beneficial for low signal-tonoise ratio situations, as detailed in the supporting information).
In a second step, we decompose the obtained displacement field, u(r) into the local effective twist angle, θ*(r) and the local moiré anisotropy magnitude and direction, κ(r) and ψ(r) respectively.
To that end, we consider the Jacobian of the transformation, = + , which is the displacement   The local anisotropy parameter κ(r) discussed here can be related to heterostrain, following the model of Kerelsky et al. (22). Here, it is assumed that one of the graphene sheets is strained with a uniaxial strain ε(r), while the other one is unaffected and only undergoes a rotation. To connect to our measurements, we note that for small average twist angles, the displacement field of the moiré lattice is related to relative displacement of the constituting layers by the following formula: where ⟨J⟩ is the Jacobian corresponding to the average angle between the layers and ~( ) is the relative displacement field experienced between the two sheets (see supporting information). The relative displacement field can be decomposed in the same way as before, where the angle corresponding to W now corresponds to the deviation of the twist angle between the two sheets from the average twist angle, and the local anisotropy κ(r) and ψ(r) obtained from the resulting scaling matrix indicate the relative strain between the layers.
Furthermore, from the resulting scaling matrix elements, we can calculate the magnitude of the strain applied to the deformed sheet, ε(r) (see supporting information). We show the resulting ε(r) in figure   4f. On average, we find that ε = 0.14% with a standard deviation of 0.09%.
It is interesting to compare the numbers for strain and twist angle heterogeneity, and their respective influence on the electronic structure of TBG. Calculations using a continuum model have considered both strain and twist angle changes in TBG samples close to the magic angle (34). It was shown that a heterostrain of ε ≈ 0.1% results in a splitting of the van Hove singularities of approximately 5 meV, comparable to the effect of variations in the twist angle of about 0.03° (35).
Furthermore, stress can cause strong qualitative changes to the electronic structure including new van Hove singularities for ε ≈ 0.5%. If we compare these numbers with our measurements, we conclude a roughly equal effect of the observed strain and twist angle inhomogeneity, suggesting that both have to be taken into account when fabricating samples, as both effects significantly alter the electronic structure compared to a perfect lattice.
Before concluding, we want to address one potential challenge of the method introduced here: it is also sensitive to piezo drift. Piezo drift occurs in STM experiments due to thermal fluctuations that influence the piezo, due to piezo relaxation after a change of field of view, or due to the piezo relaxation from the movement necessary to take the topography.

S1.1 Deformations of a lattice
We perform a lock-in measurement on an image that clearly displays a periodic lattice. In STM, this implies we can use any topography of sufficient quality that displays the crystal lattice. The idea is to use a lock-in measurement in order to find a transformation of coordinates between the measured, distorted image and its pristine, undeformed equivalent (in this work, a perfect triangular lattice). Defining the measured and pristine image as T m (r), T r (r ) respectively, both with measurement coordinates r = (x, y) ∈ R 2 and lattice coordinates r = (x , y ) ∈ R 2 , the following relation holds: T m (r) = T r (r + u(r)) = T r (f (r)) = T r (r ) = T m f −1 (r ) where the transformation from measurement coordinates to lattice coordinates is given by f (r) = r + u(r) = r . Here, u(r) is called the displacement field, as is well-established in continuum mechanics. For convenience, we also define the inverse displacement: Note that by substitution, we have the following relation between forward and inverse displacement: With this, we can express the pristine image at lattice coordinates in terms of 13 the measured image: Therefore, if we can determine u(r), and therefore u (r ), we can reconstruct the pristine image. This is the idea of the Lawler-Fujita reconstruction algorithm [1] 1 .

S1.2 Properties of the deformation
The lattice distortion u(r) as defined above, fully describes the deformation of the lattice, but does not directly provide insight into the relevant properties.
To that end, we first define the Jacobian of the transformation f : , where ∇u is the Jacobian of the displacement field, in continuum mechanics terminology the deformation gradient tensor, and in canonical terms defined as follows: In order to fully characterise the deformation of the lattice, we decompose J in its polar form: where W is the rotation matrix corresponding to the rotation of the full lattice and the matrix P describes the local anisotropy and scaling. P can be further decomposed in the rotation matrix V indicating the orientation of the axis of anisotropy (i.e. the axis of largest scaling, with the axis of smallest scaling perpendicular to it) and the diagonal scaling matrix D = d 1 0 0 d 2 , where by convention and implementation d 1 ≥ d 2 holds for any position r. The geometric mean of these directional scaling factors is equal to the square root of the determinant of D and therefore of J: √ d 1 d 2 = det(J). We use this to quantify the local twist angle, as it corresponds to the local scaling of the wavelength of the moiré lattice: Where |q j | is the length of the chosen reference vectors. This local wavelength is then converted to a local twist angle using the well-known expression: , where a = 2.46Å is the lattice constant of graphene and θ(r) the local twist angle.
A quantification of the local anisotropy is given by the ratio κ = d 1 /d 2 and the angle between the anisotropy axis and the x-axis is finally calculated from V : ψ = arctan Vxy Vxx . In our practical implementation, the singular value decomposition (SVD) is used to obtain the decomposition in equation 1 for each point in the image, and Matlab's atan2 is used to find the right quadrant of the angles from the signs of V xx and V xy .

S1.3 Determination of the displacement field u(r)
In order to determine u(r) for a certain image, we perform a lock-in measurement. To clarify, we can represent any (nearly) periodic image as: where φ j = q j · u(r) is the position-dependent phase of the lattice. The summation runs over the reciprocal lattice vectors q j (j ∈ {1, 2, 3} for a hexagonal lattice), T 0 is the constant indicating the amplitude of the modulation and u(r) is again the displacement field.
The phase is measured using standard lock-in procedure: The existing image is mixed with a reference image containing a specific plane wave. If we choose the periodicity of this reference wave equal to that of the lattice itself, we can then low-pass filter the mixed image and end up with a phase map for a specific wave. For clarification: where the cosine in the first term denotes the (real-valued) measured image, whereas the complex exponential denotes the reference wave and → denotes lowpass filtering in order to get rid of the last term between brackets, corresponding to a rotating wave approximation. Alternatively, for a gaussian low-pass filter, this corresponds to a real space gaussian integration window of the lock-in. By taking the (pointwise) angle of the complex, filtered product image, we end up with the phase map. In particular, this phase map contains information about the displacements of each pixel in the measured image T m (r) with respect to the pristine reference lattice T r (r). This procedure is repeated for at least one extra reciprocal lattice vector. The two phase maps are then used to find the displacement field u(r). Earlier, we defined the distorted coordinates as r = r + u(r). Multiplying this equation by the reciprocal lattice vectors, we get a system of equations expressing the projection of the distortion onto the reciprocal lattice vectors: q j · r = q j · r + φ j , j ∈ 1, 2, 3 Selecting only j ∈ {1, 2}, we have in matrix notation: such that we can write for φ = φ1 φ2 : Multiplying by Q −1 , we find r = r + Q −1 φ, and therefore u(r) = Q −1 φ(r).

S1.4 Additional notes on choice of reference vectors
Selecting two reference vectors To obtain u(r) as described above, we only used the phase of the lock-in signal of two reference vectors. For a triangular/hexagonal lattice, a priori three possible choices of which two reference vectors to use are possible from the three linear independent references vectors as fitted to the FFT of the image. To select which two vectors to use for the reconstruction of u(r), we either selected the ones with the largest average lockin amplitude, or by inspecting the phase-unwrapped images and selecting the ones where no remaining phase slips occured.
Using more than two reference vectors In principle, information is lost when only selecting the phase of the lock-in signal of two reference vectors to obtain u(r). Although not used in this work, in low signal-to-noise ratio situations, it could be beneficial to use all the information. Equation 2 also holds for more than two phases and reference vectors. Although Q is not square in this case, a solution can be obtained for each pixel using linear least squares minimization of the following equivalent equation: Where additionally the amplitude of the lock-in signals can be used as weights to the minimization problem.
Isotropy Enforcing the reference lattice to be isotropic can be done either in advance by enforcing isotropic reference wavevectors as applied in this work, or alternatively after the initial lock-in step, by adding an additional linear phase ∆φ j = ∆q j · r to the obtained phase, where ∆q j is the difference between the used reference wavevector and the isotropic wavevector. The advantage of the latter method would be a slightly improved signal-to-noise ratio, as the smoothing window can be centered around the actual average wavevector occuring in the image instead of the ideal, equal-length, 60 degree rotated ones.

S2 Relation of moire lattice to relative displacement
For a non-homogeneous bilayer, we can fully describe the system by two displacement fields u ↑ (r), u ↓ (r) of respectively the top and bottom layer compared to an undistorted system.
where T r (r) denote the atomic lattices, r ↓ , r ↑ the lattice coordinates of both lattices and ⊕ denotes the (as of now, unspecified) operation of the combination of two lattices into one image.
We can express the deformation of one atomic lattice w.r.t the coordinates of the other: Assert u ↓ (r) = J ↓ r + v ↓ (r), i.e. a rotation and/or scaling plus local variations. Note that here, J ↓ is constant 2-by-2 matrix corresponding to a mean ∇u, and therefore corresponding to J − I in terms of the J defined in the previous section. In this case, we have: For two real lattice plane waves T r (r ) = cos(q j · r ) and taking pointwise product for the ⊕ operator, we have: = cos(q j r ↑ ) cos(q j r ↑ + δ(r)) = 1 2 cos(2q j r ↑ + δ(r)) + 1 2 cos(−δ(r)) = 1 2 cos(2q j r ↑ + δ(r)) + 1 2 cos(+δ(r)) For the modulation δ(r) the following holds: Substituting r ↑ = r + u ↑ (r) and u ↑ (r ↑ ) = −u ↑ (r): where u ∼ (r) denotes the relative displacement between the upper layer and the rotated lower layer. Substituting back in T m : Note that for a 2D lattice consisting of the sum of 2 or more cosines, each with it's own q j , this construction can be made for each q j separately, nevertheless resulting in a single, joint u moire (r) (as expected).
For a small twist angle θ between two equal lattices, e.g. magic angle twisted bilayer graphene, we have: Therefore, in this case the topography T m (r) consists of a sum of a cosine with approximately twice the atomic frequency and a cosine with approximately θ times the atomic frequency: the moire frequency. As expected, this lattice is rotated by 90 degrees plus half the angle of the original rotation, i.e. angled halfway inbetween both lattices.

S2.1 Relation to uniaxial strain models
Graphene has a Poisson ratio δ = 0.17, so if a strain is applied in one direction, it shrinks in the perpendicular direction by δ . By applying the decomposition into θ(r), κ(r) and ψ(r) as described in Section S1.1, to the relative displacement between the layers u ∼ (r) and assuming the relative strain is dominated by the strain of one layer, we can calculate that strain (r). For uniaxial strain, we have with these assumptions in terms of the decomposition into relative displacement: and therefore we can express the strain of a single layer as follows: which can then be related to other measurements and models [2,3]. Note that the measured quantity u moire (r) is related to the relative displacement by a multiplication of J −1 . For small twist angles, J −1 ≈ 1 θ (with θ in radians, i.e. for θ = 1.05 • we have J −1 ≈ 55), strongly amplifying effects of small relative displacement.

S3 Phase unwrapping & singularities
In this work, phase unwrapping of a periodic phase is needed in two separate places: unwrapping the lock-in phase φ j (r) before reconstructing u(r), and to obtain a single valued anisotropy angle ψ(r). Our code is written in Matlab, which has a built-in function for one-dimensional phase unwrapping. The phase is unwrapped in both directions of the image. The order in which this is done Figure S1: Phase maps of the data shown in main text figure 4, figure S3a. a,b,c correspond to the phase maps of the Bragg peak labeled q1,q2 and q3 respectively (see main text figure 2b). Because the map corresponding to q2 shows some phase singularities, we use φ1 and φ2 for determining the displacement field.
usually does not matter, provided there are no phase singularities present in the image. We occasionally encountered some phase singularities in one of the three phase maps (figure S1), but we worked around this simply by using the other two phase maps in order to find the displacement field.
In case this is not an option, for example when applying this technique to a square lattice, and/or when phase slips are present in all phase maps, there are more sophisticated algorithms for phase unwrapping available: [4][5][6]. Some of these phaseslips were present in the ψ(r) maps, for example the one displayed in the main text, figure 4e. Here, we used a Matlab implementation of a least-squares based phase unwrapping algorithm [4,7].

S4 Device overview
A schematic of the devices studied in this work is presented in figure S2. More information about the actual fabrication process can be found in the main text.

S5 Validity check with more data
In this work, we claim that the contribution of piezo drift to the output of our algorithm is negligible. To verify this, we apply it to multiple topographies, all sequentially measured on the same area. All of them are measured with a scan speed of 54 nm/s, whereas the last one (figure S3f) is measured at 65 nm/s. Because piezo drift changes with time and scan speed, comparing these datasets provides us with insight to which degree the algorithm output is affected by this effect. The algorithm output for these measurements is displayed in figure S3. Figure S3a corresponds to the data shown in the main text. For completeness, we also show ξ, the angle corresponding to the matrix W (see section S1.1). Comparing these results from different scans, we observe that almost all deformations are reproduced, in particular the vertical line-like feature on the right and the two minima in κ(r). The only features not reproduced are horizontal 'creases', corresponding to line-to-line scan artefacts. Additionally, no significant difference is observed for figure S3f with the deviating scan speed comapred to the rest. From this, we conclude that most observed deformations are intrinsic to the sample.   Figure S3: Spatial lock-in output for sequentially measured topographies in the same field of view. Figure f was measured at 65nm/s, whereas a-e were measured with a scan speed of 54nm/s. The setup condition was kept constant between measurements: V = 250 mV, I = 20 pA.

S6 Accuracy of the algorithm
As an additional consistency check, we used the Lawler-Fujita algorithm to reconstruct the undistorted image [1], and then applied the algorithm on the corrected data in order to extract the residual displacement field and compare it to the previously extracted displacement field. Here, a perfectly performing and consistent algorithm would extract a zero residual displacement field. Therefore, this gives an indication of the error of the quantities extracted by the algorithm. Since we decompose the displacement field, for an almost zero displacement field, we expect the effective twist angle map to become more centered around the average twist angle (in this case, 2.02 • ). Furthermore, we expected that most of the anisotropy is gone i.e., κ → 1 and → 0. We check this using the data topography presented in the main text (figure 4) and in figure S3a, and show the results in figure S4. Aside from edge effects in the corner indeed, the residual anisotropy is less than one percent and the residual variations in the twist angle are also below a percent, both more than an order of magnitude smaller than the originally obtained values..

S7 Heterogeneity comparison with other devices
We measured 2 additional devices, with average twist angles of 2.16 • and 2.01 • . The output of the spatial lock-in algorithm for these topographies is displayed in figure S5 and figure S6. Calculating the standard deviation for the twist angle maps, we find 0.03 and 0.06 respectively which is consistent with the result presented in the main text.

S8 Data processing
Regarding data pre-processing and post-processing, we made the following manipulations: • Topographies are obtained from the measured data by subtracting a polynomial background up to 8 th order. It was verified that this did not significantly influence the extracted displacement fields.