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

Isaac M. Craig  Department of Chemistry, University of California, Berkeley
Madeline Van Winkle  Department of Chemistry, University of California,
Berkeley  Colin Ophus  The Molecular Foundry, Lawrence Berkeley National
Laboratory, Berkeley  D. Kwabena Bediako<EMAIL_ADDRESS>Department of
Chemistry, University of California, Berkeley  Chemical Sciences Division,
Lawrence Berkeley National Laboratory

###### Abstract

Bragg interferometry (BI) is an imaging technique based on four-dimensional
scanning electron microscopy (4D-STEM) wherein the intensities of select
overlapping Bragg disks are fit or more qualitatively analyzed in the context
of simple trigonometric equations to determine local stacking order. In
4D-STEM based approaches, the collection of full diffraction patterns at each
real-space position of the scanning probe allows the use of precise virtual
apertures much smaller and more variable in shape than those used in
conventional dark field imaging, such that even buried interfaces marginally
twisted from other layers can be targeted. A coarse-grained form of dark field
ptychography, BI uses simple physically derived fitting functions to extract
the average structure within the illumination region and is therefore viable
over large fields of view. BI has shown a particular advantage for selectively
investigating the interlayer stacking and associated moiré reconstruction of
bilayer interfaces within complex multi-layered structures. This has enabled
investigation of reconstruction and substrate effects in bilayers through
encapsulating hexagonal boron nitride and of select bilayer interfaces within
trilayer stacks. However, the technique can be improved to provide a greater
spatial resolution and probe a wider range of twisted structures, for which
current limitations on acquisition parameters can lead to large illumination
regions and the computationally involved post-processing can fail. Here we
analyze these limitations and the computational processing in greater depth,
presenting a few methods for improvement over previous works, discussing
potential areas for further expansion, and illustrating the current
capabilities of this approach for extracting moiré-scale strain.

††preprint: AIP/123-QED

##  I Introduction

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. Huang  _et al._ (2022); Mak and Shan (2022); Zhang  _et al._
(2018); Seyler  _et al._ (2019); Alexeev  _et al._ (2019); Jin  _et al._
(2019); Tran  _et al._ (2019); Liu  _et al._ (2021); Dandu  _et al._ (2022);
Naik  _et al._ (2022); Susarla  _et al._ (2022) Both highly tunable and
proposed to be a physical realization of the Hubbard model with experimentally
accessible transition temperatures, moiré materials are also a preferred
platform to investigate correlated effects.Balents  _et al._ (2020); Huang
_et al._ (2022); Mak and Shan (2022); Lau  _et al._ (2022); Tang  _et al._
(2020); Wang  _et al._ (2020); Regan  _et al._ (2020); Xu  _et al._ (2020); Li
_et al._ (2021a); Huang  _et al._ (2021); Xu  _et al._ (2022) However as these
materials locally reconstruct Yoo  _et al._ (2019); Weston  _et al._ (2020);
Rosenberger  _et al._ (2020); Kazmierczak  _et al._ (2021); Li  _et al._
(2021b); Sung  _et al._ (2022) in a manner that alters their electronic
behavior Shabani  _et al._ (2021); Naik  _et al._ (2020); Li  _et al._
(2021c); Naik and Jain (2018); Enaldiev  _et al._ (2020); Ferreira  _et al._
(2021), 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 Ophus  _et al._ (2014); Latychevskaia  _et
al._ (2018); Ophus (2019) that permits the direct measurement of inter-layer
strain fields. Kazmierczak  _et al._ (2021); Van Winkle  _et al._ (2023) This
approach involves fitting interference fringes in the diffraction patterns to
a known functional form to extract structural details. Kazmierczak  _et al._
(2021); Zachman  _et al._ (2021); Van Winkle  _et al._ (2023) Alternative
4D-STEM based methods typically fit the locations of Bragg disks Mukherjee
_et al._ (2020); Yuan, Zhang, and Zuo (2019); Munshi  _et al._ (2022); Sari
_et al._ (2023); Mahr  _et al._ (2021), or have analyzed a combination of both
intensity and disk locations Shi  _et al._ (2023). 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 Zeltmann  _et al._ (2020) and has not yet been feasible
in marginally twisted structures. Other ptychographic techniques similarly
analyze scattering intensities and can provide much higher resolution Jiang
_et al._ (2018); Chen  _et al._ (2020); Yang  _et al._ (2024); Pelz  _et al._
(2017), but are currently arduous to acquire and analyze over the moiré length
scale. 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 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
Kazmierczak  _et al._ (2021); Van Winkle  _et al._ (2023) data.

##  II Strain Fields in Moiré Systems

We begin with a simple picture of strain in moiré materials. We can first
assume that there is negligible hetero-strain within the moiré such that the
lattice vectors of both layers are identical apart from a small interlayer
twist angle $\theta_{t}$ (Figure 1a), resulting in the following relationship
between the atomic positions of the top layer $\textbf{r}_{ij}^{\rm{top}}$,
the bottom layer $\textbf{r}_{ij}^{\rm{bottom}}$, and the reference layer
$\textbf{r}_{ij}^{\rm{ref}}$ we take as halfway between them at each pixel
location $(i,j)$. The details of the $\textbf{r}_{ij}$ coordinates do not
factor into this analysis as we consider the strain of the locally averaged
offset between the lattice vectors of each layer.

$\displaystyle\textbf{r}_{ij}^{\rm{top}}=\begin{bmatrix}\cos(\frac{\theta_{t}}{2})&-\sin(\frac{\theta_{t}}{2})\\\
\sin(\frac{\theta_{t}}{2})&\cos(\frac{\theta_{t}}{2})\\\
\end{bmatrix}\textbf{r}_{ij}^{\rm{ref}}$

$\displaystyle\textbf{r}_{ij}^{\rm{bottom}}=\begin{bmatrix}\cos(\frac{\theta_{t}}{2})&\sin(\frac{\theta_{t}}{2})\\\
-\sin(\frac{\theta_{t}}{2})&\cos(\frac{\theta_{t}}{2})\\\
\end{bmatrix}\textbf{r}_{ij}^{\rm{ref}}$

The expression is straightforwardly adapted to account for an interlayer
lattice mismatch originating either in hetero-strain or from the use of two
chemically distinct layers. The intralayer displacement
$\textbf{u}^{\rm{top}}$ (Fig 1b) associated with the displacement of each atom
in the top layer away from the reference is therefore given by
$\textbf{u}^{\rm{top}}_{ij}=(\textbf{r}_{ij}^{\rm{top}}-\textbf{r}_{ij}^{\rm{ref}})$,
yielding the following relation.

$\displaystyle\textbf{u}_{ij}^{\rm{top}}=\begin{bmatrix}\cos(\frac{\theta_{t}}{2})-1&-\sin(\frac{\theta_{t}}{2})\\\
\sin(\frac{\theta_{t}}{2})&\cos(\frac{\theta_{t}}{2})-1\\\
\end{bmatrix}\textbf{r}_{ij}^{\rm{ref}}$

Likewise, we define
$\textbf{u}^{\rm{bottom}}_{ij}=(\textbf{r}_{ij}^{\rm{bottom}}-\textbf{r}_{ij}^{\rm{ref}})$.
These displacements can be analyzed using small displacement theory (also
known as infinitesimal strain theory) Slaughter (2002), 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
$|\nabla\textbf{u}^{\rm{top}}|_{\infty}\ll 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(\theta/2)$ as $\textbf{u}^{\rm{top}}_{ij}$ varies between 0
and $a_{0}/4$ over a length scale $a_{0}/\sin(\theta/2)$. The following
analysis therefore assumes $\sin(\theta/2)\ll 1$.

The strain tensor $\nabla\textbf{u}_{ij}^{\rm{top}}$ then consists of partial
derivatives of $\textbf{u}_{ij}^{\rm{top}}$ with respect to the reference
coordinates. For a rigidly twisted structure that perfectly obeys this
relation at each pixel, $\nabla\textbf{u}_{ij}^{\rm{top}}$ is independent of
pixel location and obeys the following uniformly.

$\displaystyle\nabla\textbf{u}^{\rm{top}}_{ij}=\begin{bmatrix}\cos(\frac{\theta_{t}}{2})-1&-\sin(\frac{\theta_{t}}{2})\\\
\sin(\frac{\theta_{t}}{2})&\cos(\frac{\theta_{t}}{2})-1\\\ \end{bmatrix}$ (1)

This rigid displacement field corresponds to a uniform local curl and
divergence of $\nabla\times\textbf{u}^{\rm{top}}_{ij}=2\sin(\theta_{t}/2)$ and
$\nabla\cdot\textbf{u}^{\rm{top}}_{ij}=2\cos(\theta_{t}/2)-2$, equal to
$\theta_{t}$ and 0 respectively to first order. Atomic reconstruction will
manifest deviations from this uniform curl and divergence, which we will
denote as rotations and dilations respectively. The use of small deformation
theory further permits us to consider reconstruction-driven changes to the
local rotation and dilation as additive changes to what is expected of a rigid
moiré. In this framework, we can decompose $\nabla\textbf{u}^{\rm{top}}$ as
the sum of the curl-free and divergence-free matrices $\epsilon$ and $\omega$,
encoding dilational and rotational strain respectively.

$\displaystyle\nabla\textbf{u}^{\rm{top}}=\underbrace{\frac{1}{2}\nabla\textbf{u}^{\rm{top}}+\frac{1}{2}(\nabla\textbf{u}^{\rm{top}})^{T}}_{\epsilon}+\underbrace{\frac{1}{2}\nabla\textbf{u}^{\rm{top}}-\frac{1}{2}(\nabla\textbf{u}^{\rm{top}})^{T}}_{\omega}$

The eigenvalues and eigenvectors of $\epsilon$ are termed the principal
stretches and their directions, and the sum and difference of these
eigenvalues are the total intralayer dilation (the first invariant, equivalent
to $\nabla\cdot\textbf{u}^{\rm{top}}$) and maximum intralayer engineering
shear strain $\gamma$ respectively. The curl of $\omega$ corresponds to the
local rotation $\nabla\times\textbf{u}^{\rm{top}}$, consisting of the sum of
the rigidly imposed twist angle and a local reconstruction rotation. For a
homo-bilayer structure, the deformation of the top and bottom layers are equal
and opposite and all reconstruction information can be determined through
consideration of a single layer. In the discussion that follows, we use u to
denote $\textbf{u}^{\rm{total}}$.

Figure 1:  (1a) Schematic illustrating the anticipated interlayer
displacements utotal within a rigid moiré structure at a twist angle of
$\theta_{t}$, for which the curl and divergence of the displacement utop
(equal to one half utotal) are uniformly $\theta_{t}$ and 0 respectively. (1b)
Definition of the total interlayer displacement utotal, where utotal = utop
and utop \- ubottom in terms of the intralayer displacements of each layer,
utop and ubottom. |utop| $\approx$ |ubottom| such that utotal $\approx$ 2utop
is assumed. (1c) Legend illustrating the coloring scheme used to denote the
atomic stacking, where AA denotes regions where the two layers are aligned
out-of-plane and AB/BA denote regions where the two layers are aligned Bernal
out-of-plane. SP1, SP2, and SP3 denote the various symmetrically distinct
saddle point stacking regions.

##  III Interference Region Intensities for Displacement Measurements

Figure 2:  (2a) Schematics illustrating the 4D-STEM Bragg interferometry
approach where a series of converged beam electron diffraction patterns are
collected at various real-space positions as the beam is scanned across the
sample. Various relevant parameters including the camera length L, convergence
semi-angle $\alpha$, and electron probe width d are noted. (2b-c) Multi-slice
CBED patterns at select u choices. Corresponding approximate real space
stackings are below. $10\bar{1}0$ and $\bar{2}110$ reflections are highlighted
in white, corresponding to those plotted in 2d. (2d) Variation in intensity of
the $10\bar{1}0$ and $1\bar{2}10$ reflections (red) and their Friedel pairs
(blue) over a densely-sampled scan of u choices. These reflections were
selected for ease as they align with the chosen cartesian axis orientation.
Both intensities show good agreement with a fit to equation 2 (solid lines)
with the need for a small offset from the origin, equal and opposite for each
Friedel pair. (2e) Corresponding root mean squared errors and offsets added to
the cosine argument are shown for the $10\bar{1}0$ line-cut (left column of
2d) in units of $\pi$ for a series of acquisition parameters. Unless otherwise
noted, multi-slice simulations were carried out using $\alpha$ = 4 mrad, a
slice thickness of 3 Å, perfect aberration correction, and the VdW gap
expected for graphene of $3.55$Å using AbTEM Madsen and Susi (2021) Figure 3:
(3a) Convergence semi-angle $\alpha$ needed to resolve AA-type domain
contractions to a given precision. The probe size is assumed to be 0.61
$\lambda/\alpha$, AA region defined as $|u|=0\pm 0.25a_{0}$, and domain size
changes are taken to be observable when larger than the beam width. The
$\alpha$ needed will be over-estimated as observed beam widths are slightly
smaller than those predicted and many independent measurements are used. (3b)
To-scale illustrations of the overlap regions expected at a series of $\alpha$
choices for a hBN, bilayer graphene hetero-structure with inter-layer twists
of 2 degrees, 15 degrees, and 17 degrees between the two graphene layers and
each graphene layer and the hBN respectively. (3c) Number of viable overlap
pixels in units of the camera length squared for a series of twist angles and
$\alpha$ choices. The two smallest of the three twist tables are labeled, with
the third implicit angle being the sum of these other two.

To extract these displacements u and associated strain $\nabla\textbf{u}$, we
look at the intensity of the overlap regions between Bragg disks arising from
the scattering off of desired layers in the material. The overlap intensities
$I_{j}$ at reciprocal lattice positions $\textbf{g}_{j}$ can be shown to obey
the following equation, where u is the (in-plane projection of the) interlayer
displacement vector dictating stacking order.

$\displaystyle I_{j}$
$\displaystyle=A_{j}+B_{j}cos^{2}(\pi\textbf{g}_{j}\cdot\textbf{u})$ (2)

The full derivation is shown in the supplemental information of Kazmierczak
_et al._ (2021) and Van Winkle  _et al._ (2023) 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.

The accuracy of this expression can be investigated through comparison to
simulated electron diffraction patterns, obtained using the multi-slice
algorithm. Cowley and Moodie (1957) In figure 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 equation (2) without
resolution considerations, which will be discussed in a later section. An
$\alpha$ of 4 mrad and a perfectly focused aberration-free lens are assumed.
We note that the primary effect of changing $\alpha$ 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 subtle internal structure in
abberated 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 Zachman  _et al._ (2021).

Fig. 2d shows a slight disagreement between Eqn 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 multi-slice
between the two layers plays a minor role. However, the expression remains
accurate up to at least gaps of 15 angstroms. 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 equation (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 equation (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.

##  IV Resolution-Driven Restrictions on Sample Morphology

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 to 1.5 nm. These probe
sizes are also larger than those seen in other electron microscopy techniques
Jiang  _et al._ (2018); Chen  _et al._ (2020); Yang  _et al._ (2024); Pelz
_et al._ (2017), 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\chi)$ where
$\chi$ is the Scherzer aberration function Scherzer (1949), and
$\alpha$-dependent phase shift introduced by the lens, yields a typical
measure for probe width. Kirkland (1998) For many applications, the maximum
tolerable phase error is then used to determine the maximal semi-convergence
angle ($\alpha$) and its associated probe width of approximately 0.61
$\lambda/\alpha$ where $\lambda$ is the wavelength of the electron probe
(Figure 1e). Weyland and Muller (2020) 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 $\alpha<10$ mrad where it is deemed tolerable
with sufficient aberration correction Weyland and Muller (2020) and consider
the number of pixels we can obtain at each $\alpha$ the practical limitations
associated with their extraction.

In BI, we often cannot realize the optimal $\alpha$ 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
Kretschmer  _et al._ (2020); Bell  _et al._ (2014); Egerton (2012). 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 James and Browning (1999).
This can be used to determine the minimum $\alpha$ 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
$|\textbf{u}|=0\pm 0.25a_{0}$ so that a rigidly twisted structure has AA
domains of width $0.25a_{0}/\sin(\theta/2)$. The $\alpha$ needed to resolve
domain contractions to a given percent of the rigid value are shown in Figure
3a. 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, $\alpha$ must be smaller than 10 mrad. To see
this, we show a series of anticipated overlap regions (Figure 3b), obtained
using the relations illustrated in Figure 1a-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 $\alpha$. 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 $\alpha$ and therefore higher spatial
resolution are possible than those used in previous works ($\alpha\approx 2-3$
mrad) Kazmierczak  _et al._ (2021); Van Winkle  _et al._ (2023) if care is
taken to define precise virtual apertures of variable shape (Figure 3c)
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
Kazmierczak  _et al._ (2021); Van Winkle  _et al._ (2023); Craig  _et al._
(2024) 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.

Figure 4:  (4a) Mean |u| errors (in units of $a_{0}$) obtained through fitting
intensities in the form of equation 2 (with A=1, B=0) using the iterative and
joint optimizations described. Poisson noise is included by drawing each disk
intensity independently from the Poisson distribution $P(k;\lambda I)=(\lambda
I)^{k}\exp(-\lambda I)/k!$ where I is the noise-free intensity. A range of
specified $\lambda$ choices are shown, which in practice will depend on the
electron dose. Noise mitigation is expected when averaging over many pixels to
obtain overlap regions, increasing the effective $\lambda$. The resulting
signal-to-noise ratio of the interferometry pattern is shown using a colored
virtual dark field in the inset. Error bars represent the standard deviation
using 25 datasets randomly generated in this manner. The result of the joint
optimization is shown as a solid line and its corresponding standard deviation
as the greyed region. (4b) The fraction of incorrectly assigned pixels using
each of the three described unwrapping procedures applied to the displacement
fields expected of a rigidly twisted structure (shown in 4d) with Gaussian
noise of zero mean and the specified standard deviation $\sigma_{u}$ (in units
of $a_{0}$). Error bars represent the standard deviation using 25 datasets
randomly generated in this manner. We note that pixels on the edge of the
field of view can be incorrectly assigned even without noise in these
approaches. (4c) Colored virtual dark fields (see Fig 4e) and displacement
field maps (visualized using the color scheme defined in Fig 1c) illustrating
the signal-to-noise ratios and fit efficacy of representative experimental
data. (4d) Colored virtual dark fields and displacement field maps (visualized
using the color scheme defined in Fig 1c) for artificial data assuming a
rigidly twisted structure and given Gaussian noise, as used in the comparison
of the unwrapping approaches in (4b). (4e) Color legend for virtual dark
fields illustrating qualitative atomic stacking classifications. Craig  _et
al._ (2024)

##  V Intensity Fitting To Determine Inter-layer Stacking

We compare multiple fitting procedures to obtain the displacements u
prescribed by equation 2. One approach we use involves fitting ujk and the
coefficients $A_{i},B_{i}$ sequentially, such that ujk at each pixel location
$(j,k)$ 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 Bezdek and Hathaway (2003) was also more
robust to experimental noise and permitted the use of many initial guesses at
each iterative step, helping to avoid local minima.

In this iterative approach, we first determine the optimal displacement vector
ujk independently at each pixel (j,k) provided the fixed coefficients
${A_{i},B_{i}}$. The experimental intensities are normalized, and
$A_{i}=0,B_{i}=1$ are assumed. We used a dense uniform grid of 36 initial
starting conditions to decrease the chance of obtaining local minima and
constrain the values of u to reside within a single unit cell such that
$\textbf{u}=c_{1}\textbf{a}_{1}+c_{2}\textbf{a}_{2}$ with $|c_{1}|\leq
1/2,|c_{2}|\leq 1/2$ in terms of the separately measurable in-plane lattice
vectors $\textbf{a}_{1}$ and $\textbf{a}_{2}$. This results in the following
over-determined equation with an easily obtainable Jacobian, which we solved
using Newton non-linear least squares independently for each pixel.

$\displaystyle\begin{bmatrix}I_{10}\\\ ...\\\ I_{hk}\\\ \end{bmatrix}$
$\displaystyle=\begin{bmatrix}A_{10}+B_{10}cos^{2}(\pi c_{1})\\\ ...\\\
A_{hk}+B_{hk}cos^{2}(\pi(hc_{1}+kc_{2}))\\\ \end{bmatrix}$ (3)

Using these ujk values, we determine the optimal coefficients ${A_{i},B_{i}}$
independently for each of the twelve Bragg disk intensities $I_{ijk}$ at each
pixel $(j,k)$ using linear least squares. Fitting the coefficients at a fixed
choice of u results in the following linear equation for each diffraction disk
i. We note that the coefficients $A_{i}$, $B_{i}$ are assumed to be
independent of pixel location over the N by N pixel field of view to avoid
over-fitting, and can be assumed to be approximately 1 and 0 respectively
following normalization of $I_{i}$.

$\displaystyle\underbrace{\begin{bmatrix}cos^{2}(\pi\textbf{g}_{i}\cdot\textbf{u}_{00})&1\\\
cos^{2}(\pi\textbf{g}_{i}\cdot\textbf{u}_{10})&1\\\ ...\\\
cos^{2}(\pi\textbf{g}_{i}\cdot\textbf{u}_{NN})&1\\\
\end{bmatrix}}_{A_{coef}}\cdot\begin{bmatrix}A_{i}\\\ B_{i}\\\
\end{bmatrix}=\underbrace{\begin{bmatrix}I_{i00}\\\ I_{i10}\\\ ...\\\
I_{iNN}\\\ \end{bmatrix}}_{b_{coef}}$ (4)

The fitting of ${A_{i},B_{i}}$ and ujk are then iterated until convergence. We
compare this staged optimization to a joint optimization in which the
coefficients and u are optimized simultaneously following one iteration of the
same multi-start ujk optimization with fixed $A_{i}=1,B_{i}=0$. The mean |u|
obtained over the full set of pixels is shown for both approaches in Fig 4a
for artificial data corresponding to a rigidly twisted structure with disk
intensities that exactly obey equation 2 with independent Poisson noise
associated with acquisition included for each disk intensity. Representative
experimental datasets are shown in Fig 4c, showing that the extent of noise
included is comparable or greater than that expected of weakly scattering
samples like graphene and much greater than that expected of transition metal
dichalcogenides (both acquired using $\alpha$ = 1.7 mrad). Fig 5a shows that
the iterative approach consistently achieves a lower error in the
displacements u (Fig 4a). Avoidance of local extrema and noise mitigation are
both active areas of research and alternative algorithms such as genetic
algorithms Mirjalili and Mirjalili (2019) could prove better than this
iterative multi-start approach, but were not required for experimental data of
this nature.

Figure 5:  (5a) Displacement field map (visualized using the color scheme
defined in Fig 1c) for a representative experimental dataset of hBN
encapsulated bilayer MoS2 with a 1.77 degree twist, fit using the described
iterative procedure. The field of view is 100 by 100 pixels, corresponding to
50 by 50 nm. (5b) Components of the strain tensor
$\nabla\textbf{u}_{ij}^{\rm{top}}$ obtained by numerically differentiating the
displacement field following unwrapping via the Watershed-AA approach. Values
provided are unit-less after accounting for the difference in scale for the
displacements ($a_{0}$ = 0.315 nm for MoS2) and pixels (step size of 0.5 nm).
Results are shown using a coordinate system where the y-axis is taken to lie
along $SP_{1}$, necessary to obtain $\nabla\cdot\textbf{u}^{\rm{top}}\approx
0$ and $\nabla\times\textbf{u}^{\rm{top}}\approx\theta_{t}$ in rigidly twisted
structures. (5c-e) The corresponding engineering sheer strain, local rotations
$\nabla\times\textbf{u}^{\rm{top}}$ (shown in degrees), and local dilation
$\nabla\cdot\textbf{u}^{\rm{top}}$

##  VI Displacement Unwrapping for Strain Extraction

The displacements obtained in the fitting procedure cannot yet be
differentiated to assess local strain. This is because the obtained u are not
uniquely determined by the data, since adding integer multiples of the lattice
vectors $\textbf{a}_{1}$ and $\textbf{a}_{2}$ to u results in the same set of
intensities due to the sample’s symmetry. For centro-symmetric materials,
flipping the sign of u also results in the same pattern. Owing to these
ambiguities, the raw displacement vectors obtained are not smoothly oriented
and cannot be naively 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 naively choose locally
optimal offsets.

We therefore seek to find the optimal $s_{jk},n_{jk},m_{jk}$ at each pixel
location $(j,k)$ given $\textbf{u}_{jk}$ such that the (squared) euclidean
distance of $\textbf{u}^{\rm{unwrap}}_{jk}$ from all four of its neighbors is
minimized, thus minimizing a quadratic cost function to variations in integer
parameters.

$\displaystyle\textbf{u}^{\rm{unwrap}}_{jk}=s_{jk}\textbf{u}_{jk}+n_{jk}\textbf{a}_{1}+m_{jk}\textbf{a}_{2}$
(5) $\displaystyle n_{jk},m_{jk}\in Z,s_{jk}=\pm 1$ (6)

Liberal use of integer quadratic programming is prohibitively expensive, so we
considerably pre-process the data so that it need only 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,m)$, such that $n_{jk},m_{jk}=n_{i},m_{i}$ for all pixels
$(j,k)$ in region i.

One partitioning approach we used involves the watershed algorithm. Kornilov
and Safonov (2018) In this method, several regions are specified (in this case
the centers of the $|\textbf{u}|\approx 0$ (AA-type) stacking regions) and
pixels of increasing $|\textbf{u}|$ are successively assigned. We denote this
as watershed-AA. The watershed-AA approach proved effective on the data sets
with large, dispersed regions of $|\textbf{u}|\approx 0$ stacking separated by
thin boundaries, as was seen for anti-parallel aligned twisted transition
metal dichalcogenides Van Winkle  _et al._ (2023). We also investigate an
alternative approach based on a Voronoi partition Okabe  _et al._ (2009),
where region boundaries are drawn halfway between the $x_{i}$ locations. This
method was successful for highly noisy data but failed when the moiré
wavelength varied even marginally over the field of view. Van Winkle  _et al._
(2023) Lastly, we investigate an approach where the watershed algorithm is
applied to the angle of u instead of its magnitude to select Bernal-type
stacking regions for structures that minimize $|\textbf{u}|\approx 0$ like
large-angle twisted graphene and parallel aligned twisted transition metal
dichalcogenides. We denote this approach as watershed-AB. After use of any of
these three segmentation approaches, the values $n_{i},m_{i}$ are determined
based on the connectivity of neighboring regions and $s_{ij}$ is chosen to
maximize a local curl. A quadratic integer program can then be applied to the
pre-processed data. Van Winkle  _et al._ (2023)

We compare these three segmentation approaches on artificial
$\textbf{u}_{jk}^{\rm{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_{jk}$ following Eqn
4. Fig 4b 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 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, uunwrap 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 $\nabla\textbf{u}^{\rm{top}}$.
In practice there is often a rotation between the real-space coordinate system
and that used to define u, which need be accounted for. Kazmierczak  _et al._
(2021); Van Winkle  _et al._ (2023)

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 MoS2 with a twist angle of 1.77 degrees using a series of coordinate
choices. The sheer strain, curl, and divergence maps are shown in Fig 5c-e,
expected to be uniformly 0, $\theta_{t}$, and 0 in a rigid structure
respectively. The clear local structure in 5c and 5d are consequences of
atomic reconstruction, showing that MoS2 distorts to increase the local twist
angle within the |u| $\approx 0$ (shown black in Fig 5a), reducing the size of
these high-energy stacking regions. Simultaneously, the distortion decreases
the local twist angle within the low-energy stacking regions (Fig 5d). These
two rotational effects act in a manner that accumulates shear strain along
their boundaries (Fig 5c) while the dilational strain associated with uniform
stretching of the lattice remains negligible (Fig 5e). Van Winkle  _et al._
(2023)

##  VII Conclusion

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.

###### Acknowledgements.

This material is based upon work supported by the US 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
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.

## Conflict of Interest

The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are available within the
article.

## References

  * Huang  _et al._ (2022) D. Huang, J. Choi, C.-K. Shih, and X. Li, “Excitons in semiconductor moiré superlattices,” Nat. Nanotechnol. 17, 227–238 (2022).
  * Mak and Shan (2022) K. F. Mak and J. Shan, “Semiconductor moiré materials,” Nat. Nanotechnol. 17, 686–696 (2022).
  * Zhang  _et al._ (2018) N. Zhang, A. Surrente, M. Baranowski, D. K. Maude, P. Gant, A. Castellanos-Gomez, and P. Plochocka, “Moiré intralayer excitons in a $\mathrm{MoSe_{2}}$/$\mathrm{MoS_{2}}$ heterostructure,” Nano Lett. 18, 7651–7657 (2018).
  * Seyler  _et al._ (2019) K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao, and X. Xu, “Signatures of moiré-trapped valley excitons in $\mathrm{MoSe_{2}}$/$\mathrm{WSe_{2}}$ heterobilayers,” Nature 567, 66–70 (2019).
  * Alexeev  _et al._ (2019) E. M. Alexeev, D. A. Ruiz-Tijerina, M. Danovich, M. J. Hamer, D. J. Terry, P. K. Nayak, S. Ahn, S. Pak, J. Lee, J. I. Sohn, _et al._ , “Resonantly hybridized excitons in moiré superlattices in van der waals heterostructures,” Nature 567, 81–86 (2019).
  * Jin  _et al._ (2019) C. Jin, E. C. Regan, A. Yan, M. Iqbal Bakti Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, _et al._ , “Observation of moiré excitons in $\mathrm{WSe_{2}}$/$\mathrm{WS_{2}}$ heterostructure superlattices,” Nature 567, 76–80 (2019).
  * Tran  _et al._ (2019) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, _et al._ , “Evidence for moiré excitons in van der waals heterostructures,” Nature 567, 71–75 (2019).
  * Liu  _et al._ (2021) E. Liu, E. Barré, J. van Baren, M. Wilson, T. Taniguchi, K. Watanabe, Y.-T. Cui, N. M. Gabor, T. F. Heinz, Y.-C. Chang, _et al._ , “Signatures of moiré trions in $\mathrm{WSe_{2}}$/$\mathrm{MoSe_{2}}$ heterobilayers,” Nature 594, 46–50 (2021).
  * Dandu  _et al._ (2022) M. Dandu, G. Gupta, P. Dasika, K. Watanabe, T. Taniguchi, and K. Majumdar, “Electrically tunable localized versus delocalized intralayer moiré excitons and trions in a twisted $\mathrm{MoS_{2}}$ bilayer,” ACS Nano 16, 8983–8992 (2022).
  * Naik  _et al._ (2022) M. H. Naik, E. C. Regan, Z. Zhang, Y.-H. Chan, Z. Li, D. Wang, Y. Yoon, C. S. Ong, W. Zhao, S. Zhao, _et al._ , “Intralayer charge-transfer moiré excitons in van der waals superlattices,” Nature 609, 52–57 (2022).
  * Susarla  _et al._ (2022) S. Susarla, M. H. Naik, D. D. Blach, J. Zipfel, T. Taniguchi, K. Watanabe, L. Huang, R. Ramesh, F. H. da Jornada, S. G. Louie, _et al._ , “Hyperspectral imaging of excitons within a moiré unit-cell with a sub-nanometer electron probe,” Preprint at arXiv:2207.13823 (2022).
  * Balents  _et al._ (2020) L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, “Superconductivity and strong correlations in moiré flat bands,” Nat. Phys. 16, 725–733 (2020).
  * Lau  _et al._ (2022) C. N. Lau, M. W. Bockrath, K. F. Mak, and F. Zhang, “Reproducibility in the fabrication and physics of moiré materials,” Nature 602, 41–50 (2022).
  * Tang  _et al._ (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, _et al._ , “Simulation of hubbard model physics in $\mathrm{WSe_{2}}$/$\mathrm{WS_{2}}$ moiré superlattices,” Nature 579, 353–358 (2020).
  * Wang  _et al._ (2020) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, _et al._ , “Correlated electronic phases in twisted bilayer transition metal dichalcogenides,” Nat. Mater. 19, 861–866 (2020).
  * Regan  _et al._ (2020) E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, _et al._ , “Mott and generalized wigner crystal states in $\mathrm{WSe_{2}}$/$\mathrm{WS_{2}}$ moiré superlattices,” Nature 579, 359–363 (2020).
  * Xu  _et al._ (2020) Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, “Correlated insulating states at fractional fillings of moiré superlattices,” Nature 587, 214–218 (2020).
  * Li  _et al._ (2021a) H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, _et al._ , “Imaging two-dimensional generalized wigner crystals,” Nature 597, 650–654 (2021a).
  * Huang  _et al._ (2021) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, _et al._ , “Correlated insulating states at fractional fillings of the $\mathrm{WS_{2}}$/$\mathrm{WSe_{2}}$ moiré lattice,” Nat. Phys. 17, 715–719 (2021).
  * Xu  _et al._ (2022) Y. Xu, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, “A tunable bilayer hubbard model in twisted $\mathrm{WSe_{2}}$,” Nat. Nanotechnol. 17, 934–939 (2022).
  * Yoo  _et al._ (2019) H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, _et al._ , “Atomic and electronic reconstruction at the van der waals interface in twisted bilayer graphene,” Nat. Mater. 18, 448–453 (2019).
  * Weston  _et al._ (2020) A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, _et al._ , “Atomic reconstruction in twisted bilayers of transition metal dichalcogenides,” Nat. Nanotechnol. 15, 592–597 (2020).
  * Rosenberger  _et al._ (2020) M. R. Rosenberger, H.-J. Chuang, M. Phillips, V. P. Oleshko, K. M. McCreary, S. V. Sivaram, C. S. Hellberg, and B. T. Jonker, “Twist angle-dependent atomic reconstruction and moiré patterns in transition metal dichalcogenide heterostructures,” ACS Nano 14, 4550–4558 (2020).
  * Kazmierczak  _et al._ (2021) N. P. Kazmierczak, M. Van Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe, and D. K. Bediako, “Strain fields in twisted bilayer graphene,” Nat. Mater. 20, 956–963 (2021).
  * Li  _et al._ (2021b) H. Li, S. Li, M. H. Naik, J. Xie, X. Li, J. Wang, E. Regan, D. Wang, W. Zhao, S. Zhao, _et al._ , “Imaging moiré flat bands in three-dimensional reconstructed $\mathrm{WSe_{2}}$/$\mathrm{WS_{2}}$ superlattices,” Nat. Mater. 20, 945–950 (2021b).
  * Sung  _et al._ (2022) S. H. Sung, Y. M. Goh, H. Yoo, R. Engelke, H. Xie, K. Zhang, Z. Li, A. Ye, P. B. Deotare, E. B. Tadmor, _et al._ , “Torsional periodic lattice distortions and diffraction of twisted 2d materials,” arXiv preprint arXiv:2203.06510 (2022).
  * Shabani  _et al._ (2021) S. Shabani, D. Halbertal, W. Wu, M. Chen, S. Liu, J. Hone, W. Yao, D. N. Basov, X. Zhu, and A. N. Pasupathy, “Deep moiré potentials in twisted transition metal dichalcogenide bilayers,” Nat. Phys. 17, 720–725 (2021).
  * Naik  _et al._ (2020) M. H. Naik, S. Kundu, I. Maity, and M. Jain, “Origin and evolution of ultraflat bands in twisted bilayer transition metal dichalcogenides: Realization of triangular quantum dots,” Phys. Rev. B 102, 075413 (2020).
  * Li  _et al._ (2021c) E. Li, J.-X. Hu, X. Feng, Z. Zhou, L. An, K. T. Law, N. Wang, and N. Lin, “Lattice reconstruction induced multiple ultra-flat bands in twisted bilayer $\mathrm{WSe_{2}}$,” Nat. Commun. 12, 1–7 (2021c).
  * Naik and Jain (2018) M. H. Naik and M. Jain, “Ultraflatbands and shear solitons in moiré patterns of twisted bilayer transition metal dichalcogenides,” Phys. Rev. Lett. 121, 266401 (2018).
  * Enaldiev  _et al._ (2020) V. Enaldiev, V. Zólyomi, C. Yelgel, S. Magorrian, and V. Fal’ko, “Stacking domains and dislocation networks in marginally twisted bilayers of transition metal dichalcogenides,” Phys. Rev. Lett. 124, 206101 (2020).
  * Ferreira  _et al._ (2021) F. Ferreira, S. Magorrian, V. Enaldiev, D. Ruiz-Tijerina, and V. Fal’ko, “Band energy landscapes in twisted homobilayers of transition metal dichalcogenides,” Appl. Phys. Lett. 118, 241602 (2021).
  * Ophus  _et al._ (2014) C. Ophus, P. Ercius, M. Sarahan, C. Czarnik, and J. Ciston, “Recording and using 4d-stem datasets in materials science,” Microscopy and Microanalysis 20, 62–63 (2014).
  * Latychevskaia  _et al._ (2018) T. Latychevskaia, C. R. Woods, Y. B. Wang, M. Holwill, E. Prestat, S. J. Haigh, and K. S. Novoselov, “Convergent beam electron holography for analysis of van der waals heterostructures,” Proceedings of the National Academy of Sciences 115, 7473–7478 (2018).
  * Ophus (2019) C. Ophus, “Four-dimensional scanning transmission electron microscopy (4d-stem): From scanning nanodiffraction to ptychography and beyond,” Microsc. Microanal. 25, 563–582 (2019).
  * Van Winkle  _et al._ (2023) M. Van Winkle, I. M. Craig, S. Carr, M. Dandu, K. C. Bustillo, J. Ciston, C. Ophus, T. Taniguchi, K. Watanabe, A. Raja, _et al._ , “Rotational and dilational reconstruction in transition metal dichalcogenide moiré bilayers,” Nature communications 14, 2989 (2023).
  * Zachman  _et al._ (2021) M. J. Zachman, J. Madsen, X. Zhang, P. M. Ajayan, T. Susi, and M. Chi, “Interferometric 4d-stem for lattice distortion and interlayer spacing measurements of bilayer and trilayer 2d materials,” Small 17, 2100388 (2021).
  * Mukherjee  _et al._ (2020) D. Mukherjee, J. T. Gamler, S. E. Skrabalak, and R. R. Unocic, “Lattice strain measurement of core@ shell electrocatalysts with 4d scanning transmission electron microscopy nanobeam electron diffraction,” ACS Catalysis 10, 5529–5541 (2020).
  * Yuan, Zhang, and Zuo (2019) R. Yuan, J. Zhang, and J.-M. Zuo, “Lattice strain mapping using circular hough transform for electron diffraction disk detection,” Ultramicroscopy 207, 112837 (2019).
  * Munshi  _et al._ (2022) J. Munshi, A. Rakowski, B. H. Savitzky, S. E. Zeltmann, J. Ciston, M. Henderson, S. Cholia, A. M. Minor, M. K. Chan, and C. Ophus, “Disentangling multiple scattering with deep learning: application to strain mapping from electron diffraction patterns,” npj Computational Materials 8, 254 (2022).
  * Sari  _et al._ (2023) B. Sari, S. E. Zeltmann, C. Zhao, P. M. Pelz, A. Javey, A. M. Minor, C. Ophus, and M. C. Scott, “Analysis of strain and defects in tellurium-wse2 moiré heterostructures using scanning nanodiffraction,” ACS nano 17, 22326–22333 (2023).
  * Mahr  _et al._ (2021) C. Mahr, K. Müller-Caspary, T. Grieb, F. F. Krause, M. Schowalter, and A. Rosenauer, “Accurate measurement of strain at interfaces in 4d-stem: A comparison of various methods,” Ultramicroscopy 221, 113196 (2021).
  * Shi  _et al._ (2023) C. Shi, N. Mao, K. Zhang, T. Zhang, M.-H. Chiu, K. Ashen, B. Wang, X. Tang, G. Guo, S. Lei, _et al._ , “Domain-dependent strain and stacking in two-dimensional van der waals ferroelectrics,” Nature communications 14, 7168 (2023).
  * Zeltmann  _et al._ (2020) S. E. Zeltmann, A. Müller, K. C. Bustillo, B. Savitzky, L. Hughes, A. M. Minor, and C. Ophus, “Patterned probes for high precision 4d-stem bragg measurements,” Ultramicroscopy 209, 112890 (2020).
  * Jiang  _et al._ (2018) Y. Jiang, Z. Chen, Y. Han, P. Deb, H. Gao, S. Xie, P. Purohit, M. W. Tate, J. Park, S. M. Gruner, _et al._ , “Electron ptychography of 2d materials to deep sub-ångström resolution,” Nature 559, 343–349 (2018).
  * Chen  _et al._ (2020) Z. Chen, M. Odstrcil, Y. Jiang, Y. Han, M.-H. Chiu, L.-J. Li, and D. A. Muller, “Mixed-state electron ptychography enables sub-angstrom resolution imaging with picometer precision at low dose,” Nature communications 11, 2994 (2020).
  * Yang  _et al._ (2024) W. Yang, H. Sha, J. Cui, L. Mao, and R. Yu, “Local-orbital ptychography for ultrahigh-resolution imaging,” Nature Nanotechnology , 1–6 (2024).
  * Pelz  _et al._ (2017) P. M. Pelz, W. X. Qiu, R. Bücker, G. Kassier, and R. D. Miller, “Low-dose cryo electron ptychography via non-convex bayesian optimization,” Scientific reports 7, 9883 (2017).
  * Slaughter (2002) W. S. Slaughter, “Kinematics,” in  _The Linearized Theory of Elasticity_ (Birkhäuser, Boston, MA, 2002) pp. 97–155.
  * Madsen and Susi (2021) J. Madsen and T. Susi, “The abtem code: transmission electron microscopy from first principles,” Open Research Europe 1 (2021).
  * Cowley and Moodie (1957) J. M. Cowley and A. F. Moodie, “The scattering of electrons by atoms and crystals. i. a new theoretical approach,” Acta Crystallographica 10, 609–619 (1957).
  * Scherzer (1949) O. Scherzer, “The theoretical resolution limit of the electron microscope,” Journal of Applied Physics 20, 20–29 (1949).
  * Kirkland (1998) E. J. Kirkland,  _Advanced computing in electron microscopy_ , Vol. 12 (Springer, 1998).
  * Weyland and Muller (2020) M. Weyland and D. A. Muller, “Tuning the convergence angle for optimum stem performance,” arXiv preprint arXiv:2008.12870 (2020).
  * Kretschmer  _et al._ (2020) S. Kretschmer, T. Lehnert, U. Kaiser, and A. V. Krasheninnikov, “Formation of defects in two-dimensional mos2 in the transmission electron microscope at electron energies below the knock-on threshold: the role of electronic excitations,” Nano letters 20, 2865–2870 (2020).
  * Bell  _et al._ (2014) D. C. Bell, M. Mankin, R. W. Day, and N. Erdman, “Successful application of low voltage electron microscopy to practical materials problems,” Ultramicroscopy 145, 56–65 (2014).
  * Egerton (2012) R. Egerton, “Mechanisms of radiation damage in beam-sensitive specimens, for tem accelerating voltages between 10 and 300 kv,” Microscopy research and technique 75, 1550–1556 (2012).
  * James and Browning (1999) E. James and N. Browning, “Practical aspects of atomic resolution imaging and analysis in stem,” Ultramicroscopy 78, 125–139 (1999).
  * Craig  _et al._ (2024) I. M. Craig, M. Van Winkle, C. Groschner, K. Zhang, N. Dowlatshahi, Z. Zhu, T. Taniguchi, K. Watanabe, S. M. Griffin, and D. K. Bediako, “Local atomic stacking and symmetry in twisted graphene trilayers,” Nature Materials , 1–8 (2024).
  * Bezdek and Hathaway (2003) J. C. Bezdek and R. J. Hathaway, “Convergence of alternating optimization,” Neural, Parallel & Scientific Computations 11, 351–368 (2003).
  * Mirjalili and Mirjalili (2019) S. Mirjalili and S. Mirjalili, “Genetic algorithm,” Evolutionary algorithms and neural networks: Theory and applications , 43–55 (2019).
  * Kornilov and Safonov (2018) A. S. Kornilov and I. V. Safonov, “An overview of watershed algorithm implementations in open source libraries,” Journal of Imaging 4, 123 (2018).
  * Okabe  _et al._ (2009) A. Okabe, B. Boots, K. Sugihara, and S. N. Chiu, “Spatial tessellations: concepts and applications of voronoi diagrams,” (2009).


# Atomistic origins of the phase transition mechanism in Ge2Sb2Te5

Juarez L. F. Da Silva,1 Aron Walsh,1 Su-Huai Wei,1 and Hosun Lee2 1National
Renewable Energy Laboratory, 1617 Cole Blvd., Golden, CO 80401, USA
2Dept. of Applied Physics, Kyung Hee University, Suwon 446-701, South Korea

###### Abstract

Combined static and molecular dynamics first-principles calculations are used
to identify a direct structural link between the metastable crystalline and
amorphous phases of Ge2Sb2Te5. We find that the phase transition is driven by
the displacement of Ge atoms along the rocksalt [111] direction from the
stable-octahedron to high-energy-unstable tetrahedron sites close to the
intrinsic vacancy regions, which give rise to the formation of local 4-fold
coordinated motifs. Our analyses suggest that the high figures of merit of
Ge2Sb2Te5 are achieved from the optimal combination of intrinsic vacancies
provided by Sb2Te3 and the instability of the tetrahedron sites provided by
GeTe.

Crystalline-amorphous phase transition, mechanism, density functional theory,
Ge2Sb2Te5

###### pacs:

61.43.-j,61.50.Ks,71.15.Nc

††preprint: GST compounds

Ternary (GeTe)m(Sb2Te3)n materials, in particular the Ge2Sb2Te5 (GST)
composition, have been considered as the most natural candidates for non-
volatile memory applications through exploiting the fast and reversible
resistance change between a metastable (m-GST) crystalline phase (low
resistivity) and an amorphous (a-GST) phase (high resistivity). Ovshinsky
(1968); Yamada et al. (1991); Wuttig and Yamada (2007); Chong et al. (2008)
However, the mechanism of the phase transition is still under intense debate.
The existing models,Kolobov et al. (2004, 2006); Welnic et al. (2006); Hegedüs
and Elliott (2008) have provided a preliminary understanding of the transition
mechanism, but fail to provide a clear and direct structural link between the
m-GST and a-GST phases, which play a key role in the understanding of the
reservible transition at an atomistic level. The m-GSTYamada and Matsunaga
(2000); Njoroge et al. (2002); Matsunaga and Yamada (2004); Park et al.
(2005); Sun et al. (2006); Wuttig et al. (2007); Shportko et al. (2008) phase
crystallizes in a rocksalt-type (RS) structure, in which the Te atoms occupy
the anion sites and Ge, Sb, and the naturally occurring intrinsic vacancies
from Sb2Te3 (20% in GST) occupy the cation sites. It has been suggested that
a-GST is characterized by the presence of 4-fold coordinated Ge atoms,Kolobov
et al. (2004); Baker et al. (2006); Akola and Jones (2007); Caravati et al.
(2007); Hegedüs and Elliott (2008); Jóv’ari et al. (2008) in which the sum of
the occurrences GeTe4, Ge(SbTe3), and Ge(GeTe3) is about 66%.Jóv’ari et al.
(2008) Ge$-$Ge and Ge$-$Sb bonds are found in those motifs, which is assumed
to be due to disorder effects, since they are not present in the crystalline
phases. Furthermore, a-GST shows a volume expansion of $6-7$% compared with
the m-GST phase, Weidenhof et al. (1999); Njoroge et al. (2002) and it has a
higher energy ($28-40$ meV/atom) with respect to m-GST.Kalb et al. (2003)
Theoretically, first-principles molecular dynamics (MD) starting from a liquid
phase with slow cooling rates have been used to generate a metastable phase
(RS-type structure), however, no direct transition path was identified to link
the proposed m-GST and a-GST phases. Thus, a new approach that connects the
two phases at the atomistic level becomes highly desirable.

Figure 1:  (Color online) Structure models of the GST phases. (a) Metastable
crystalline GST (m-GST). (b) m-GST (shift) structure, in which the Ge atoms
occupy the 4-fold tetrahedron sites with lowest energy. (c) Amorphous GST
obtained at zero temperature (am-GST) using modified m-GST structures (m-GST
with Ge shift), in which the tetrahedron Ge sites were initially occupied. (d)
Amorphous GST obtained by high temperature molecular dynamics DFT calculations
(a-GST). The Ge, Sb, and Te atoms are indicated in green, blue, and red,
respectively.

In this work, using first-principles methods, we will address the following
open questions: Is there a dominant structure link between both phases? What
are the roles of GeTe and Sb2Te3 in GST? We will show in this Letter that
combined static (zero temperature) and MD (high-temperature) first-principles
calculations can explain the phase transition mechanism between the m-GST and
a-GST phases. Moreover, our study shows that generating the amorphous phase
from a known crystalline phase provides a better understanding of the
structural relationship between both phases. Thus, it provides a new avenue
for further study of amorphous materials phase change transitions.

Our static total energy and MD calculations are based on the all-electron
projected augmented wave (PAW) methodBlöchl (1994); Kresse and Joubert (1999)
and density functional theory (DFT) within the generalized gradient
approximation (GGA-PBE)Perdew et al. (1996) as implemented in VASP. Kresse and
Hafner (1993); Kresse and Furthmüller (1996) To represent the metastable phase
(RS-type structure), we employ a hexagonal $(2{\times}2{\times}1)$ unit cell,
in which the Te atoms are stacked along of the [0001] direction. Da Silva et
al. (2008) The MD calculations were performed employing cubic and hexagonal
cells with 108 to 126 atoms. The total energies and equilibrium volumes for
all structures in both crystalline and amorphous phases were obtained by full
relaxation of the volume, shape, and atomic positions of the unit cell to
minimize the quantum mechanical stresses and forces.

To understand the phase transition, we first established the crystal structure
of the m-GST phase,Da Silva et al. (2008) as shown in Fig. 1. The obtained
structure is consistent with experimental results and provides new insights
into m-GST.Yamada and Matsunaga (2000); Matsunaga and Yamada (2004); Park et
al. (2005); Kolobov et al. (2004); Matsunaga et al. (2007) In this layered-
structure the ordered intrinsic vacancies separate the building block units
(GST), in which the Ge and Sb atoms are intermixed in planes. All Ge atoms are
6-fold coordinated in m-GST. However, it has been reported that up to one
fifth of the Ge atoms are 4-fold coordinated with Te atoms in a-GST (GeTe4),
while the remaining Ge atoms form 4-fold motifs with combined Ge, Sb, and Te
atoms.Jóv’ari et al. (2008) We notice that tetrahedral Ge atoms can be
obtained by shifting the octahedral Ge atoms in m-GST along the hexagonal $c$
direction, i.e., there are two tetrahedron sites for each Ge atom, Fig. 1. In
order to identify the lowest energy tetrahedron sites, we calculated the
energetics for the occupation of each site by Ge atoms. The lowest energy
sites are located in the intrinsic vacancy regions, while the highest energy
sites are located in the center of the GST building blocks, i.e., there is a
strong preference for the four-fold Ge atoms to be located in or near
intrinsic vacancy regions. Assuming that all the Ge atoms are shifted from
their octahedron sites and occupy the lowest energy tetrahedron sites, we find
that 50% of Ge will shift from the octahedra to tetrahedra along the [0001]
direction, while the remaining 50% shift along of the opposite direction. The
m-GST (shift) structure in which all Ge atoms occupy the tetrahedral sites
according to the distribution of intrinsic vacancies and energy barriers is
shown in Fig. 1b, which leads to the formation of Ge$-$Ge bonds. This
configuration is highly unstable, and the system will relax without energy
barrier to a lower energy phase (see am-GST structure in Fig. 1c).

Figure 2:  Pair-correlation functions of various GST phases. Amorphous GST obtained by molecular dynamic calculations (a-GST, black lines). Amorphous GST obtained from occupation of tetrahedron sites in m-GST and complete relaxation (am-GST, red lines). Meta-stable GST phase (m-GST, blue lines, scaled by 0.50). Table 1:  Bond lengths (in Å) of Ge2Sb2Te5 (GST) in the amorphous and crystalline phases.  | Amorphous GST | Meta-stable GST
---|---|---
| a-GST | am-GST | Exp. | m-GST | Exp.
Ge$-$Te | 2.74 | 2.79 |  $2.60-2.63$111Exp. Reference Jóv’ari et al., 2008; Kolobov et al., 2004; Baker et al., 2006. | $2.87-3.24$ |  $2.83-3.15$222Exp. Reference Kolobov et al., 2004.
Sb$-$Te | 2.91 | 2.96 |  $2.82-2.85$111Exp. Reference Jóv’ari et al., 2008; Kolobov et al., 2004; Baker et al., 2006. | $2.96-3.30$ | 2.91222Exp. Reference Kolobov et al., 2004.
Te$-$Te | 4.16 | 4.28 |  | $4.14-4.38$ | 4.26222Exp. Reference Kolobov et al., 2004.
Ge$-$Ge | 2.63 | 2.64 |  $2.47-2.48$111Exp. Reference Jóv’ari et al., 2008; Kolobov et al., 2004; Baker et al., 2006. | $4.27-4.62$ |
Ge$-$Sb | 2.79 | 2.78 | 2.69111Exp. Reference Jóv’ari et al., 2008; Kolobov et al., 2004; Baker et al., 2006. | $4.23-4.53$ |
Sb$-$Sb | 2.93 | 2.92 |  | $4.27-4.62$ |

To provide a more direct structural link between the m-GST and a-GST phases,
we first generated a-GST structures using first-principles MD simulations at
high temperatures, $T$, using the same approach adopted in previous a-GST
studies.Akola and Jones (2007); Sun et al. (2007); Caravati et al. (2007);
Hegedüs and Elliott (2008) Secondly, we generated several amorphous structures
from modified m-GST structures, in which a percentage of the Ge atoms (100%,
75%, 50%, 25%) are shifted to the tetrahedron sites from the lower energy
octahedron sites (m-GST with shift Ge). The goal is to show that amorphous
structures obtained in this way (am-GST in Fig. 1) can preserve most of the
structural features present in the a-GST generated by conventional MD
calculations (a-GST in Fig. 1), and therefore provide a direct structural link
and solid evidence to support the mechanisms that determines the phase
transition from m-GST to a-GST. The amorphous structures obtained by both
approaches are shown in Fig. 1. To quantify our analysis, we calculated the
pair correlation (PC) functions, which are shown in Fig. 2. For the a-GST
structures, the PC functions were averaged over five structures, while the PC
functions of the am-GST structures were calculated for ten structures with
different initial occupation of the Ge tetrahedron sites; the structure that
provided the best agreement with the a-GST PC functions is shown in Fig. 2.

Figure 3:  Potential energy path for atomic displacements of Ge atoms along of
the rocksalt (RS) [111] direction of GeTe. (a) Distorted RS structure. (b)
Perfect RS structure. (c) Long Ge$-$Te bonds zincblende (ZB) structure. (d)
Graphite like structure. (e) Perfect ZB structure. (f) RS-layer structure.

Our PC function analysis shows that for all the am-GST structures, the one in
which 50% of the Ge atoms are shifted from the octahedron to the tetrahedron
sites along the hexagonal [0001] direction and the rest 50% moves along the
$[000\bar{1}]$ direction reproduces almost all features present in the PC
functions of a-GST, although some minor differences still exist. Furthermore,
even minor features are well-described by both structures, with the formation
of Ge$-$Ge bonds and cavity regions, both of which have been identified as key
characteristics of a-GST.Baker et al. (2006); Akola and Jones (2007); Caravati
et al. (2007); Hegedüs and Elliott (2008) We observe that am-GST structures in
which less than half of the Ge atoms are moved to the tetrahedron sites do not
yield PC functions similar to the a-GST structures, instead they show strong
similarity to the PC function calculated for m-GST (see Fig. 2). Furthermore,
we observed that only the am-GST structures in which the Ge atoms initially
occupy four-fold sites in or near the intrinsic vacancies lead to structure
properties in good agreement with the calculated MD a-GST structures. Thus, it
suggests that the location of the intrinsic vacancies plays an important role
in the phase transition, which can be explained by the lower energy barriers
for Ge displacements close to intrinsic vacancy regions. For the lowest energy
m-GST structures, in which the intrinsic vacancies are ordered in a plane
perpendicular to $c$. However, at high temperature or under non-equilibrium
growth conditions the intrinsic vacancies may distribute more randomly among
the cation sites, which is expected to play an important role in the pattern
of shifted Ge atoms from their stable octahedra.

Our predicted results are in good agreement with available experimental data.
For example, using the calculated equilibrium volumes for both phases, we
obtained a density of 5.89 g/cm3 (m-GST) and 5.35 g/cm3 (a-GST and am-GST),
i.e. the amorphization gives rise to a volume expansion, which decreases the
density by about 9.20%. The experimentally observed expansion is on the order
of 6.4%.Njoroge et al. (2002) The volume expansion upon amorphization is a
consequence of the Ge atoms moving to the lower coordination sites in the
a-GST structures. Therefore, the smaller volume deformation observed in the
experimental sample may indicate that the amorphization process in not
complete, Yamada et al. (1991); Jóv’ari et al. (2008) i.e. not all the Ge
atoms are moved away from their stable octahedron sites.

Comparison of the total energies reveals that the a-GST structure is about
$140-182$ meV/atom higher in energy than the lowest energy m-GST structure,
which corresponds to the energy limit between the fully amorphized (100% shift
of the Ge atoms) and the ordered m-GST structure. Differential scanning
calorimetry measurements obtained $28-42$ meV/atom.Kalb et al. (2003) We found
that the calculated energy differences decrease by about 30 meV/atom if the
intrinsic vacancies become disordered in m-GST. Furthermore, the energy
difference could be much smaller (e.g. about 50 meV/atom) if only a fraction
of the Ge atoms undergo site transitions. This again suggests that full scale
amorphorization of GST or a complete ordering of Ge, Sb, and intrinsic
vacancies in m-GST may not be typical in the GST phases.

The averaged bond lengths calculated for both phases are summarized in Table 1
along with available experimental results. Jóv’ari et al. (2008); Kolobov et
al. (2004); Baker et al. (2006) The calculated bond lengths deviate by about
$3-6$% from the experimental results; however, most of the error is due to the
use of GGA in our calculations, which systematically overestimates the lattice
constants by about 3%. Furthermore, it is important to notice that the
nearest-neighbor distances are spread over a large range of values, e.g.
Ge$-$Te is from 2.67 to 2.94 Å and Sb$-$Te is from 2.86 to 3.23 Å. We found a
contraction in the averaged Ge$-$Te bond lengths in a-GST of up to 10%
compared with m-GST, e.g. Ge$-$Te decreases from $2.87-3.24$ Å (m-GST) to
$2.67-2.94$ Å (a-GST), while experimental measurements obtained a decrease of
about 12%. Similar trends exist for Sb$-$Te.

To understand the relaxation effects introduced by the shift of Ge atoms from
octahedron to tetrahedron sites, we calculated the potential energy path along
the RS [111] direction for GeTe as a function of Ge shift from octahedron
(perfect RS) to the tetrahedron (zincblende) sites. The results are shown in
Fig. 3. As expected, the distorted RS structure has the lowest energy (54 meV
lower than the perfect RS structure), in which the distortion is driven by
Peierls-type level repulsion near the band edge. Unexpectedly, the zincblende
(ZB) structure in which the Ge atoms occupy the tetrahedron sites with bond
angles (Ge$-$Te$-$Ge) of $109.47^{\circ}$, is not a local minimum as would be
expected based on the general trends for binary semiconductors. In fact, we
found that the ZB structure relaxes without energy barrier to the ‘graphite-
like’ or to the ‘long-Ge$-$Te’ structures, which have lower energies than the
high-symmetry ZB phase. Thus, Ge at ideal tetrahedral sites are intrinsically
unstable in GeTe, which drives the Ge atoms at tetrahedral sites in GST to
move away and adopt a variety of lower symmetry coordination environments.

The variety of coordination environments found in the GeTe energy surface is
remarkable. From Ge site occupation of (0.25,0.25,0.25) to (0.40,0.40,0.40),
three structures have similar energies, i.e. ‘long-Ge$-$Te’, layered-ZB, and
layered-RS. In ‘long-Ge$-$Te’, the Ge atoms form three short bonds (2.77 Å)
and one long bond (4.68 Å) with the Te atoms. However, Ge is only three-fold
coordinated in the layered structures with bond lengths of 2.76 Å, which is
2.90% (14.51%) smaller than the short (longer) Ge$-$Te bond lengths in the
distorted RS structure. As the layered GeTe structures are lower in energy
than the graphite-like phase and only about 100 meV/f.u. higher than the
distorted RS structure, it indicates a strong tendency of Ge atoms to form
four-fold motifs with three short Ge$-$Te bonds (about 2.76 Å) and bond angles
of about $90^{\circ}$. Similar results, e.g. short bond lengths and average
bond angles of about $90^{\circ}$, are observed by our calculations for a-GST,
which is also consistent with previous MD results for a-GST,Akola and Jones
(2007); Caravati et al. (2007); Hegedüs and Elliott (2008) as well as by
experimental observations.Kolobov et al. (2004, 2006); Jóv’ari et al. (2008)
Therefore, the inherent instability of Ge at the tetrahedral sites, low
displacement energy, and unique coordination preferences of GeTe plays an
important role in the formation of a-GST.

In summary, using first-principles calculations, we obtained a direct
structural link between the meta-stable and amorphous GST phases, as well as
the role of the parent compounds. The Sb2Te3 provides intrinsic lattice
vacancies, while GeTe contributes its RS-type structure in which Ge
displacements along the RS [111] direction can be realized at low energy cost.
The instability at the tetrahedral sites leads to the generation of disordered
GST structures in which the Ge atoms are mostly four-fold coordinated with
three short Ge$-$Te bond lengths. As the displacement has the lowest energy
near intrinsic vacancy sites, our analysis suggests that a high degree of
amorphization can be achieved most easily when the system has a composition of
(GeTe)2(Sb2Te3), i.e., is consistent with the observation that GST has the
highest figure of merit of all Ge$-$Sb$-$Te compounds. Moreover, we show that
generating amorphous materials directly from its crystalline counterpart
provides a better approach to understand these type of phase transitions
present in phase change materials.

## References

  * Ovshinsky (1968) S. R. Ovshinsky, Phys. Rev. Lett. 21, 1450 (1968).
  * Yamada et al. (1991) N. Yamada, E. Ohno, K. Nishiuchi, N. Akahira, and T. Takao, J. Appl. Phys. 69, 2849 (1991).
  * Wuttig and Yamada (2007) M. Wuttig and N. Yamada, Nature Mater. 6, 824 (2007).
  * Chong et al. (2008) T. C. Chong, L. P. Shi, X. Q. Wei, R. Zhao, H. K. Lee, P. Yang, and A. Y. Du, Phys. Rev. Lett. 100, 136101 (2008).
  * Kolobov et al. (2004) A. V. Kolobov, P. Fons, A. I. Frenkel, A. L. Ankudinov, J. Tominaga, and T. Uruga, Nature Mater. 3, 703 (2004).
  * Kolobov et al. (2006) A. V. Kolobov, J. Haines, A. Pradel, M. Ribes, P. Fons, J. Tominaga, Y. Katayama, T. Hammouda, and T. Uruga, Phys. Rev. Lett. 97, 035701 (2006).
  * Welnic et al. (2006) W. Welnic, A. Pamungkas, R. Detemple, C. Steimer, S. Blügel, and M. Wuttig, Nature Mater. 5, 56 (2006).
  * Hegedüs and Elliott (2008) J. Hegedüs and S. R. Elliott, Nature Mater. 7, 399 (2008).
  * Yamada and Matsunaga (2000) N. Yamada and T. Matsunaga, J. Appl. Phys. 88, 7020 (2000).
  * Njoroge et al. (2002) W. K. Njoroge, H.-W. Wöltgens, and M. Wuttig, J. Vac. Sci. Technol. A 20, 230 (2002).
  * Matsunaga and Yamada (2004) T. Matsunaga and N. Yamada, Phys. Rev. B 69, 104111 (2004).
  * Park et al. (2005) Y. J. Park, J. Y. Lee, M. S. Youm, Y. T. Kim, and H. S. Lee, J. Appl. Phys. 97, 093506 (2005).
  * Sun et al. (2006) Z. Sun, J. Zhou, and R. Ahuja, Phys. Rev. Lett. 96, 055507 (2006).
  * Wuttig et al. (2007) M. Wuttig, D. Lüsebrink, D. Wamwangi, W. Welnic, M. Gillessen, and R. Dronskowski, Nature Mater. 6, 122 (2007).
  * Shportko et al. (2008) K. Shportko, S. Kremers, M. Woda, D. Lencer, J. Robertson, and M. Wuttig, Nature Mat. 7, 653 (2008).
  * Baker et al. (2006) D. A. Baker, M. A. Paesler, G. Lucovsky, S. C. Agarwal, and P. C. Taylor, Phys. Rev. Lett. 96, 255501 (2006).
  * Akola and Jones (2007) J. Akola and R. O. Jones, Phys. Rev. B 76, 235201 (2007).
  * Caravati et al. (2007) S. Caravati, M. Bernasconi, T. D. Kühne, M. Krack, and M. Parrinello, Appl. Phys. Lett. 91, 171906 (2007).
  * Jóv’ari et al. (2008) P. Jóv’ari, I. Kaban, J. Steiner, B. Beuneu, A. Schöps, and M. A. Webb, Phys. Rev. B 77, 035202 (2008).
  * Weidenhof et al. (1999) V. Weidenhof, I. Friedrich, S. Ziegler, and M. Wuttig, J. Appl. Phys. 86, 5879 (1999).
  * Kalb et al. (2003) J. Kalb, F. Spaepen, and M. Wuttig, J. Appl. Phys. 93, 2389 (2003).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 48, 13115 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Da Silva et al. (2008) J. L. F. Da Silva, A. Walsh, and H. Lee, Phys. Rev. B 78, 224111 (2008).
  * Matsunaga et al. (2007) T. Matsunaga, R. Kojima, N. Yamada, K. Kifune, Y. Kubota, and M. Takata, Appl. Phys. Lett. 90, 161919 (2007).
  * Sun et al. (2007) Z. Sun, J. Zhou, and R. Ahuja, Phys. Rev. Lett. 98, 055505 (2007).