# Simulated TEM imaging of a heavily irradiated metal

Daniel R. Mason<EMAIL_ADDRESS>UK Atomic Energy Authority, Culham
Centre for Fusion Energy, Oxfordshire OX14 3DB, United Kingdom  Max
Boleininger<EMAIL_ADDRESS>UK Atomic Energy Authority, Culham
Centre for Fusion Energy, Oxfordshire OX14 3DB, United Kingdom  Jack Haley
<EMAIL_ADDRESS>UK Atomic Energy Authority, Culham Centre for Fusion
Energy, Oxfordshire OX14 3DB, United Kingdom  Eric Prestat
<EMAIL_ADDRESS>UK Atomic Energy Authority, Culham Centre for Fusion
Energy, Oxfordshire OX14 3DB, United Kingdom  Guanze He
<EMAIL_ADDRESS>Department of Materials, University of Oxford,
Parks Road, Oxford, OX1 3PH, United Kingdom  Felix Hofmann
<EMAIL_ADDRESS>Department of Engineering, University of Oxford,
Parks Road, Oxford, OX1 3PJ, United Kingdom  Sergei L. Dudarev
<EMAIL_ADDRESS>UK Atomic Energy Authority, Culham Centre for Fusion
Energy, Oxfordshire OX14 3DB, United Kingdom

###### Abstract

We recast the Howie-Whelan equations for generating simulated transmission
electron microscope (TEM) images, replacing the dependence on local atomic
displacements with atomic positions only. This allows very rapid computation
of simulated TEM images for arbitrarily complex atomistic configurations of
lattice defects and dislocations in the dynamical two beam approximation.
Large-scale massively-overlapping cascade simulations performed with molecular
dynamics, are used to generate representative high-dose nanoscale irradiation
damage in tungsten at room temperature, and we compare the simulated TEM
images to experimental TEM images with similar irradiation and imaging
conditions. The simulated TEM shows ‘white-dot’ damage in weak-beam dark-field
imaging conditions, in line with our experimental observations and as expected
from previous studies, and in bright-field conditions a dislocation network is
observed. In this work we can also compare the images to the nanoscale lattice
defects in the original atomic structures, and find that at high dose the
white spots are not only created by small dislocation loops, but rather arise
from nanoscale fluctuations in strains around curved sections of dislocation
lines.

##  I Introduction

The thermomechanical properties of structural materials are strongly dependent
on the material microstructure, yet materials proposed for advanced nuclear
fission and fusion must tolerate exposure to irradiation, which introduces
nanoscale defects changing the microstructure _in service_ [1, 2]. There is,
therefore, a real interest in characterizing the evolution of microstructure
as a function of irradiation dose and temperature. A popular and long-
established tool for investigating nanoscale defects characteristic of
irradiated microstructures is conventional TEM, as it offers a direct window
onto the microstructure at the micrometre scale with sub-nanometre
resolution[3, 4]. Many features of irradiated microstructure can be identified
unambiguously with TEM– voids[5, 6], bubbles[7, 8], large dislocation loops[9,
10], dislocation lines[11, 12], stacking faults[13], stacking fault
tetrahedra[14, 15] and second-phase particles[16]. But small features of a few
nanometres in size are much more difficult to characterize. Notwithstanding
the issue that a small cluster of defected atoms may image too faintly to be
readily detected, by eye or otherwise, those small features which do image
brightly typically appear as featureless spots, offering no real insight into
their true atomic nature. These small features, characteristic of radiation
damage, are often referred to as ‘black-dot’ damage when observed in bright-
field conditions (or ‘white-dot’ in dark-field)[17, 18, 19]. The small
features are of fundamental importance to understanding microstructural
evolution, for as well as being the characteristic state for low temperature
irradiation, quasi-independent prismatic dislocation loops are the first
defects formed during irradiation[20, 21], and they also form the building
blocks for defect annealing and coarsening.

Image simulations suggest a single, isolated black-dot is consistent with a
single, perfect, isolated dislocation loop. Experimentally, we sometimes find
dots that become invisible at a particular set of imaging $\bf{g}$-vectors,
suggesting that $\bf{g}\cdot\bf{b}=0$ in these conditions, and hence
indicating a loop with a single Burgers vector[4]. But this invisibility
criterion is only really true in the limit of long-range linear elasticity,
and cannot be relied on in cases where the defect is much smaller than this
elastic limit[4, 22]. This identification is further complicated by the fact
that a single, perfect, isolated dislocation loop is often the one thing the
spot cannot possibly be: molecular dynamics simulations show that small
dislocation loops in Body-Centre Cubic (BCC) metals are extremely mobile, with
a diffusion constant order 0.01 µm2/s[23], and so prismatic loops would find
their way to the TEM foil surface before observation is possible[24]. This
contradiction has been explained previously by assuming the loop is pinned by
impurity atoms[25, 26], or by elastic interactions with other lattice
defects[27]. Molecular dynamics simulations of overlapping cascades have
offered another possibility– that the defects are not simple loops with one
Burgers vector, but may be complex objects[28, 29]. High-dose, massively
overlapping cascade simulations have suggested that the transition from simple
loop-like defects to complex network dislocation microstructure may happen at
a dose of 0.01 to 0.1 dpa [30, 31, 32]. The true nature of the black dots is
therefore uncertain.

Rather than adopting the standard approach of simulating TEM images of known,
isolated defects, and checking they image in a way consistent with
experimental observations[33], in this work we perform image simulations of
characteristic irradiated microstructure and view them dispassionately. While
this approach is similar in motivation to that employed by Schäublin et al.
previously[34]; our work has important computational differences allowing us
to explore larger simulation sizes with a wider range of defect environments.
In section II, we rederive the equations for simulating dynamical two-beam
imaging. We recast the partial differential equations, replacing the
dependence on an atomistic displacement field with a dependence only on atomic
position. This is more physically meaningful for arbitrary atomic
configurations typically generated using molecular dynamics simulations where
the reference lattice may be difficult to define, and is computationally very
efficient, generating a simulated image of a million atom configuration in
seconds.

We emphasize that more sophisticated (and therefore more computationally
expensive) models for generating simulated TEM or STEM images exist. ABTEM
[35] and PRISTMATIC [36] consider the dispersion of electrons during their
propagation through the material, requiring a multi-slice simulation
technique. While more accurate, this added level restricts the number of atoms
which can be handled on a desktop computer, and requires more user input in
the form of testing for convergence in the solution.

In section III we describe large-scale MD simulations used to produce
characteristic high-dose irradiated microstructures. Similar simulations have
been shown previously to have a maximum hydrogen retention capacity in
agreement with hydrogen plasma-loading experiments[37], a qualitative strain
response in agreement with X-ray diffraction measurements[31], ion-beam mixing
in agreement with Rutherford Backscattering (channelling) experiments[38, 39,
40], and thermal conductivity in agreement with transient grating spectroscopy
experiments[41]. For this work we use simulation boxes of 21 M atoms; from a
computational perspective, this size is necessary to minimise elastic periodic
image effects from the defects generated, but more importantly the cell side
length (70 nm) is directly comparable to the typical thickness of a TEM
transparent foil, and the simulation boxes contain sufficient defects to
enable a non-trivial statistical image analysis. An image calculation of this
kind takes a couple of minutes on a desktop computer, and has good parallel
scaling (currently to 128 cores) on a computer cluster.

Finally, in section IV we use a JEOL 2100 TEM equipped with a LaB6 source to
image tungsten foils irradiated to 1 dpa with 20 MeV self ions at room
temperature. These irradiation experiments are, as close as possible, a match
to the simulations. We show that there is a very good qualitative agreement
between simulated and experimental TEM images. But when we look closer at the
_known_ nanoscale defects in the microstructure responsible for the simulated
image, we see a rather weak correlation between the position of brightly
imaging spots and the dislocation loops. Rather we find that, in dark field
conditions, the brightly imaging regions are due to any strain fields which
fluctuate at the nanometre scale. These strain fields are generated by all
defects, be they simple or complex loops or parts of the dislocation line
network. We conclude that black-dot damage is consistent with a range of
complex microstructural features with a characteristic nanometre scale. And so
for characterizing a high dose microstructure, we show by direct simulation
that it is necessary to consider both dynamical and weak beam conditions.

##  II Howie-Whelan approximation for generating simulated TEM images

In this section we rederive the Howie-Whelan equations for dynamical two-beam
imaging in a form which allows us to use atomic positions only, rather than a
strain or displacement field, computed in the elasticity theory approximation.
We start with the Schrödinger equation to describe the propagation of high-
energy electrons:

$-{\hbar^{2}\over 2m}{\partial^{2}\over\partial{\bf r}^{2}}\Psi({\bf
r})+U({\bf r})\Psi({\bf r})={\hbar^{2}{\bf k}^{2}\over 2m}\Psi({\bf r}),$ (1)

where $\Psi({\bf r})$ is the one-electron wave function of a high-energy
electron, as a function of position $\bf{r}$, ${\bf k}$ is the wave vector of
incident electrons, and $U({\bf r})$ is the potential energy of interaction
between the high-energy electron and the atoms. Although the above equation is
taken in the non-relativistic form, in the treatment of electron diffraction
the relativistic effects can be accounted for by replacing $m$ with the
relativistic electron mass [42].

We look for the solution in the form of a sum of propagating and diffracted
beams, where the amplitude of each varies slowly as a function of spatial
coordinates

$\Psi({\bf r})=\Phi_{0}({\bf r})\exp(i{\bf k}\cdot{\bf r})+\Phi_{\bf g}({\bf
r})\exp[i({\bf k}+{\bf g})\cdot{\bf r}],$ (2)

where ${\bf g}$ is a reciprocal lattice vector used for imaging.

A dark field TEM image of microstructure is given by the intensity
distribution $I_{\bf g}(x,y)$ computed as $|\Phi_{\bf g}(x,y,L)|^{2}$ at the
exit surface of the foil at $z=L$.

Substituting (2) into (1), we find

$\displaystyle{\hbar^{2}{\bf k}^{2}\over 2m}\exp(i{\bf k}\cdot{\bf
r})\Phi_{0}({\bf r})+{\hbar^{2}({\bf k}+{\bf g})^{2}\over 2m}\exp[i({\bf
k}+{\bf g})\cdot{\bf r}]\Phi_{\bf g}({\bf r})$ $\displaystyle-i{\hbar^{2}\over
m}\exp(i{\bf k}\cdot{\bf r})\left({\bf k}\cdot{\partial\over\partial{\bf
r}}\right)\Phi_{0}({\bf r})-i{\hbar^{2}\over m}\exp[i({\bf k}+{\bf
g})\cdot{\bf r}]\left(({\bf k}+{\bf g})\cdot{\partial\over\partial{\bf
r}}\right)\Phi_{\bf g}({\bf r})$ $\displaystyle+U({\bf
r})\left\\{\Phi_{0}\exp(i{\bf k}\cdot{\bf r})+\Phi_{\bf g}\exp[i({\bf k}+{\bf
g})\cdot{\bf r}]\right\\}$ $\displaystyle=$ $\displaystyle{\hbar^{2}{\bf
k}^{2}\over 2m}\left\\{\Phi_{0}\exp(i{\bf k}\cdot{\bf r})+\Phi_{\bf
g}\exp[i({\bf k}+{\bf g})\cdot{\bf r}]\right\\}$

Separating terms associated with either $\exp(i{\bf k}\cdot{\bf r})$ and
$\exp(i[{\bf k}+{\bf g})\cdot{\bf r}]$, and choosing the direction of the $z$
axis in the direction of ${\bf k}$, where $|{\bf k}|\gg|{\bf g}|$, we arrive
at a system of coupled equations for the amplitudes of the transmitted and
diffracted beams

$\displaystyle{\partial\over\partial z}\Phi_{0}({\bf r})$ $\displaystyle=$
$\displaystyle-i{U_{-{\bf g}}\over\hbar v}\exp[i{\bf g}\cdot{\bf u}({\bf
r})]\Phi_{\bf g}({\bf r}),$ $\displaystyle{\partial\over\partial z}\Phi_{\bf
g}({\bf r})$ $\displaystyle=$ $\displaystyle-i{\epsilon_{\bf g}\over\hbar
v}\Phi_{\bf g}({\bf r})-i{U_{{\bf g}}\over\hbar v}\exp[-i{\bf g}\cdot{\bf
u}({\bf r})]\Phi_{0}({\bf r}),$

where $v=\hbar k/m$ is the velocity of electrons, ${\bf u}({\bf r})$ is the
field of atomic displacements at ${\bf r}$, and $U_{\bf g}$ is the Fourier
component of the periodic potential in an ideal crystal with no atomic
distortions, ie

$U({\bf r})=\sum_{\bf h}U_{\bf h}\exp(i{\bf h}\cdot{\bf r}),$

where summation over ${\bf h}$ is performed over reciprocal lattice vectors.
Parameter

$\epsilon_{\bf g}=\frac{\hbar^{2}({\bf k}+{\bf
g})^{2}}{2m}-\frac{\hbar^{2}{\bf k}^{2}}{2m}$

characterises the deviation of the orientation of the incident electron beam
from the exact Bragg condition. In applications, it is often advantageous to
use imaging conditions where the magnitude of $\epsilon_{\bf g}$ is
substantial, even though this implies that the amplitude of $\Phi_{\bf g}({\bf
r})$ is relatively small and the overall intensity of the diffraction image is
lower.

The derivation of equations (LABEL:eqn:two_beam1) assumes that the notion of
the Fourier component of the potential $U_{\bf g}$ is still well defined, and
the field of displacements ${\bf u}({\bf r})$ varies on the scale much larger
than the size of the unit cell. This condition is often satisfied even in the
core regions of defects and dislocations [43, 44, 45].

In many cases of practical relevance $U_{{\bf g}}=U_{-{\bf g}}<0$, and it is
convenient to define the so-called extinction distance $\xi_{\bf g}$

$\xi_{\bf g}=\frac{\pi\hbar v}{|U_{\bf g}|},$

characterising the spatial scale of variation of the solution describing
diffraction of high-energy electrons in the crystal. At the exact Bragg
condition $\epsilon_{\bf g}=0$, and distance $\xi_{\bf g}$ corresponds to a
half the period of oscillation of solutions of (LABEL:eqn:two_beam1).

Finally, as a purely technical note, we observe that in literature [13, 46],
the energy parameter $\epsilon_{\bf g}$ is often replaced by another parameter
$s_{\bf g}$ that also characterizes the deviation of the direction of the
incident beam from the exact Bragg diffraction condition. This parameter is
defined as

$s_{\bf g}=-\frac{1}{2\pi}\frac{\epsilon_{\bf g}}{\hbar v}.$

Using the above notation, and also adding the term describing the effect of
the average constant crystal potential, equations (LABEL:eqn:two_beam1) can be
written in the form [47]

$\displaystyle{\partial\over\partial z}\Phi_{0}({\bf r})$ $\displaystyle=$
$\displaystyle{i\pi\over\xi_{0}}\Phi_{0}({\bf r})+{i\pi\over\xi_{\bf
g}}\exp[i{\bf g}\cdot{\bf u}({\bf r})]\Phi_{\bf g}({\bf r}),$
$\displaystyle{\partial\over\partial z}\Phi_{\bf g}({\bf r})$ $\displaystyle=$
$\displaystyle\left({i\pi\over\xi_{0}}+i2\pi s_{\bf g}\right)\Phi_{\bf g}({\bf
r})+{i\pi\over\xi_{{\bf g}}}\exp[-i{\bf g}\cdot{\bf u}({\bf r})]\Phi_{0}({\bf
r}).$

Note that in our notations, the reciprocal lattice vectors already include the
factor of $2\pi$ that in earlier literature on electron microscopy is still
written separately in the relevant formulae [47, 46].

To integrate equations (II), we need smoothly varying fields of atomic
displacements. Consider atoms at positions $\mathbf{R}_{j}$, where
$j=\\{1,2,\ldots,N\\}$. If we imagine perfect reference positions for these
atoms, $\mathbf{R}^{(0)}_{j}$, then small displacements
$\mathbf{u}_{j}=\mathbf{R}_{j}-\mathbf{R}^{(0)}_{j}$ are well defined. But
near a dislocation core, displacements are not small, and even the reference
lattice site for a given atom may be difficult to uniquely identify. To solve
this problem, note that to integrate equation II we actually only need the
_phase factor_ $\exp\left[i\mathbf{g}\cdot\mathbf{u}(\bf r)\right]$.

If our imaging vector $\bf{g}$ is a reciprocal lattice vector satisfying

$\exp\left(i\mathbf{g}\cdot\mathbf{R}^{(0)}_{j}\right)=1,$ (6)

we find that at atom positions
$\exp\left(i\mathbf{g}\cdot\mathbf{u}_{j}\right)=\exp\left(i\mathbf{g}\cdot\mathbf{R}_{j}\right)$.
Therefore, to provide an unambiguous answer for the smoothly varying field of
displacements needed to integrate equation II, we can interpolate the phase
factor from atomic positions themselves, and not the atomic displacements, by
interpolating the value of $\exp\left(i\mathbf{g}\cdot\mathbf{R}_{j}\right)$,

$x(\mathbf{r})=\sum_{j}\kappa(\mathbf{r};\mathbf{R}_{j})\exp\left(i\mathbf{g}\cdot\mathbf{R}_{j}\right),$
(7)

where $\kappa$ is a suitable local interpolation kernel normalised such that
$|x(\mathbf{r})|=1$. We use the Gaussian form

$\kappa(\mathbf{r};\mathbf{R}_{j})\sim\exp\left(-\frac{|\mathbf{r}-\mathbf{R}_{j}|^{2}}{2\sigma^{2}}\right)H\left(3\sigma-|\mathbf{r}-\mathbf{R}_{j}|\right),$

where $H(r)$ is the Heaviside function and the lengthscale $\sigma$ is taken
to be approximately half the unit cell parameter.

One can then reformulate equations (II) using a gauge transformation defined
by the formulae

$\displaystyle\Phi_{0}({\bf r})$ $\displaystyle=$ $\displaystyle\phi_{0}({\bf
r}),$ $\displaystyle\Phi_{\bf g}({\bf r})$ $\displaystyle=$ $\displaystyle
x^{-1}(\mathbf{r})\phi_{\bf g}({\bf r}).$ (8)

This does not affect the images since $|\Phi_{\bf g}({\bf r})|^{2}=|\phi_{\bf
g}({\bf r})|^{2}$. As a final detail, we note that in a void region electron
waves propagate freely. If we write an atomic density function $\rho(\bf{r})$
defined as equal to one in the crystal and zero in void regions, then we find
that equations (II) acquire the final form to be numerically integrated for
each pixel in the output image:

${\partial\over\partial z}\left(\begin{array}[]{c}\phi_{0}(z)\\\
\phi_{\bf{g}}(z)\end{array}\right)=i\pi\left[\left(\begin{array}[]{c c}0&0\\\
0&2s_{\bf{g}}\end{array}\right)+\rho(\bf{r})\left(\begin{array}[]{c
c}1/\xi_{0}&1/\xi_{\bf{g}}\\\
1/\xi_{\bf{g}}&1/\xi_{0}\end{array}\right)+\rho(\bf{r})\left(\begin{array}[]{c
c}0&0\\\ 0&x^{-1}(\mathbf{r}){\partial\over\partial
z}x(\mathbf{r})\end{array}\right)\right]\left(\begin{array}[]{c}\phi_{0}(z)\\\
\phi_{\bf{g}}(z)\end{array}\right)$ (9)

As an aside, we note that we can re-express the phase-factor term using

$x^{-1}(\mathbf{r}){\partial\over\partial
z}x(\mathbf{r})=i\mathbf{g}\cdot\left[\mathbf{F}(\bf{r})\,\hat{\mathbf{z}}\right],$

where $\mathbf{F}(\bf{r})$ is the deformation gradient tensor [48] at position
$\bf{r}$ with elements $F_{ij}=\partial{u}_{i}/\partial{x}_{j}$, and
$\hat{\mathbf{z}}$ is a unit vector pointing along $z$. This shows explicitly
that it is the deformation gradient tensor projected onto the imaging vector
${\bf g}$ and the direction of propagation of the electron beam $\hat{\bf z}$,
rather than atomic displacements as such, that is responsible for the observed
variation of intensity in TEM images of defects.

We can recognise the three matrix terms in equation 9 as the partial
differential equation coefficients for propagation through vacuum, propagation
through perfect crystal, and the effective change in the local lattice
orientation due to the strain field. Lattice strain and rotation enters the
equations, through the lattice deformation gradient tensor, via the third term
in (9) that is analogous to the deviation $s_{\bf g}$ from the exact Bragg
diffraction conditions entering the first term in the same equation. The
effect of local lattice distortions on the propagation of high-energy
electrons is therefore exactly equivalent to the effect of variation of the
local Bragg condition.

Furthermore, our phase field $x(\mathbf{r})$ has an associated multi-valued
coordinate-dependent displacement field, $\mathbf{\cal U}({\bf r})$, defined
by $x(\mathbf{r})=\exp[i\mathbf{g}\cdot\mathbf{\cal U}(\mathbf{r})]$. Near
atom sites, this displacement field satisfies

$\mathbf{g}\cdot\mathbf{\cal U}(\mathbf{r})=\mathbf{g}\cdot\mathbf{\cal
U}(\mathbf{R}_{j}+\bm{\delta})=\mathbf{g}\cdot\mathbf{u}_{j}+\mathbf{g}\cdot\left(\mathbf{F}\,{\bm{\delta}}\right)+2n\pi,$

where $n$ is an integer. This shows that by introducing a single-valued phase
field $x(\mathbf{r})$, we are able to truly circumvent the problem of having a
multi-valued displacement field near a dislocation core.

We use a fourth-order Runge-Kutta integration scheme with four samples of
$x(z)$ per unit cell length. The extinction distances are computed from the
crystal structure factors $|U_{\bf{g}}|$ in the Doyle-Turner approximation
[49, 50] including a finite temperature Debye-Waller factor [51].

More technical details about our implementation of these equations are given
in the Appendix. We validate the correct functioning of the code by comparing
to TEMACI [52] in the supplementary material, section X.1. Choosing a foil
orientation for a good two-beam condition is discussed in section X.2. The
geometry of how the imaging space relates to the frame represented by the
input atomic positions is in section X.3.

##  III Generation of high-dose simulated microstructure

To be able to compare simulated TEM with ground truth high dose network
dislocation microstructures, we need to generate representative atomic
configurations. Recently it has been shown that massively overlapping MD
cascade simulations [38, 53, 39, 40, 54] produce high dose microstructures
with characteristic properties consistent with several independent
experiments: Zhang et al. [55] and Markelj et al. [56] have shown simulations
agree with Rutherford Back Scattering Channeling measurements; Mason et al.
have shown that changes in lattice strain in high dose simulations of tungsten
are consistent with micro-Laue X-ray diffraction [31], the vacancy content in
high dose simulations of tungsten is consistent with the observed saturated
deuterium retention [37], and the thermal diffusivity consistent with
transient grating spectroscopy experiments [41]. But these are all
calculations of ensemble properties of the simulated microstructures, and here
we need to consider the spatial variation in the microstructure.

We therefore use here the simulated irradiation microstructures generated for
an earlier paper by some of the authors, Ref [57], where full details of the
simulation methodology can be found. Briefly, a large single crystal box of 21
M tungsten atoms was subjected to multiple 10 keV PKA collision cascades using
the MD code LAMMPS [58], with traction-free periodic boundary conditions,
until an NRT dose [59] of 1 dpa was reached at order $6\times 10^{5}$ cascades
total.

These simulation boxes have a good size for our purpose, as the simulation
cell side of 70 nm is comparable to TEM foil thickness. Sub-optimally they
have periodic boundary conditions in all directions, rather than a free
surface, so there are no image forces on the microstructure [60], and no loop
loss to the surface [61]. To generate the images, we take multiple periodic
replicas in the plane normal to the beam direction, and a single replica
parallel to the beam, essentially creating a new, unrelaxed surface. The
geometry of the image generation is discussed in section X.3.

In these and similar massively overlapping cascade simulations in bcc metals,
the microstructure at low dose (order 0.01 dpa) shows small, isolated,
interstitial character dislocation loops. The size and the Burgers vector of
these loops is somewhat potential dependent [54]– here we used the potential
of Ref [62]. At a critical dose (order 0.1 dpa), the density of the loops gets
so high that they are no longer separated, but rather form complex networks.
At high dose, the loops have become so large that they span the simulation
cell, and can be better described as planes of perfect crystal cut by
dislocation lines. At all doses we see some dislocation loops, but we note
that MD simulations have a short timescale, order tens of nanoseconds, and so
these isolated loops might be able to diffuse and further coalesce if evolved
over a longer timescale or at higher temperature.

##  IV Experimental imaging of high-dose microstructure

TEM imaging of dislocation loops was carried out on a set of tungsten samples
irradiated to 0.01,0.1, and 1.0 dpa. These samples were electro-polished and
then irradiated by 20 MeV self ions at room temperature. Images were taken
using weak beam dark field conditions on a JEOL-2100 TEM (accelerating voltage
200 kV, LaB6 source). All images were taken with $\mathbf{g}=[200]$ and zone
axis $[001]$ diffraction conditions. A detailed description can be found in
[63].

##  V Results and Discussion

In this section we show images of simulated high-dose microstructure and
compare to overlays of the ground-truth dislocations and to experimental TEM
images.

The principal result of this work is shown in figure 1. This shows the
evolution of the microstructure from dislocation loops through a network to
dislocation lines. Kinematical bright-field and weak-beam dark-field images
were generated with $s_{\bf g}\approx 0.03\rm{nm}^{-1}$ and $s_{\bf g}\approx
0.20\rm{nm}^{-1}$ respectively, corresponding to $n_{\bf g}=1.50$ and $n_{\bf
g}=6.25$. A description of how the foil orientation is chosen to find these
deviation parameters is given in appendix X.2. In these, and subsequent
simulated TEM images, we use an accelerating voltage of 200 keV, and assume a
sample temperature of 300K for the extinction distances. The bright-field
images can be compared directly to the images of high-dose, room-temperature
tungsten irradiation by El-Atwani et al. [6, 64] who observed spatially
ordered structures at the tens of nanometre lengthscales. An image from Ref
[64] is reproduced in figure 2. We see the same qualitative pattern of a high
density of curved dislocation lines in both experiment, fig 2 and the
corresponding dose in simulation, fig 1. El Atwani et al. attributed the
spatial ordering to loop rafting [65], where mobile loops with the same
Burgers vector orient and align themselves to reduce their mutual elastic
interaction energy. Here, we have found the complex dislocation network is
generated at low temperature in such a way as to lower elastic energy density,
which leads to orientation and alignment of parts of the the network without
the need for assuming the existence of independent mobile dislocation loops.
This is in better agreement with the network growth mechanism suggested by
Wang et al. [66]. Note that our simulations are low-temperature irradiation,
and so would not show the loop rafting mechanism where the defect density is
lower at high temperatures.

Figure 1: Snapshots of massively overlapping cascade simulations with zone
axis close to the [001] direction. At each dose, the three panels are Left:
simulated kinematical bright-field two-beam image with $\mathbf{g}=[110]$
horizontal; Centre: weak-beam dark field image; Right: dislocation network
shown computed with Ovito[67] ( green lines $1/2\langle 111\rangle$, pink
lines $\langle 100\rangle$ ).

The image size of the spots seen at low dose in weak beam conditions in figure
1 is very similar to the dislocation loop size. But at high dose it is clear
that the spot size in weak beam conditions is correlated with the curved
sections of dislocation lines. Individual small dislocation loops are present
in the high-dose simulation images but image much more faintly in comparison.
A discussion of brightly and faintly imaging defects in otherwise perfect
crystal is given in section X.1.

Figure 2: Bright-field experimental image of tungsten, ion-irradiated at room
temperature to 0.25 dpa, reproduced from Ref [64]. Note the simulated images
in this work are about 100 nm across, the size of the scale bar in this image.

To understand the difference in the bright-field and dark-field images, we
consider a high-dose microstructure snapshot at a range of deviation
parameters. Figure 3 shows a series of images with the same $\mathbf{g}=[200]$
imaging vector and zone axis close to [001], with varying fine tilt in the
g-vector direction of the order of a few degrees to change the value of
$s_{g}$. The images go from dynamical two beam condition to weak beam
condition. The dose for this sequence of images is 0.1 dpa, where the
microstructure is one of densely packed complex network dislocations and
dislocation loops. From equation 9, we see that in kinematical conditions,
$|x^{-1}{\partial x}/\partial z|\gg|s_{\bf g}|$, image intensity associated
with the long range elastic fields is seen. As the foil is tilted to weak-beam
condition, $|s_{\bf g}\gg|x^{-1}{\partial x}/\partial z|$, we only see short
range intensity peaks at the positions of the highest strain. The spatial
correlation between the simulated TEM image and the ground truth dislocation
configuration is equally good in each image, we simply highlight different
magnitudes of the strain field.

Figure 3: Snapshots of overlapping cascade simulations at 0.1 dpa, with zone
axis close to the [001] direction, computed for different deviation parameters
$s_{g}$. Top dynamical two-beam conditions, rotating to bottom weak beam
conditions. Dislocation network computed with Ovito[67] overlaid coloured by
Burgers vector: green ${}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\langle
111\rangle$ type and pink $\langle 100\rangle$. All images shown with inverted
contrast. Note we have exploited the simulation box periodic boundary
conditions to extend the images.

The sequence of images in figure 3 demonstrates how the weak beam image offers
a correct but incomplete view of the microstructure, and interpreting a high
dose microstructure as consisting of quasi-independent dislocation loops from
the weak beam image alone can be misleading. Figure 4 makes a direct
comparison between experimental weak beam TEM images of ion-irradiated
tungsten with the corresponding simulated TEM images. As the experimental
images cover a much larger area with a long spatial correlation length
compared to the simulations, we artificially increase the apparent simulated
area by stitching together six symmetry-related ${\bf g}$-vector - zone axis
orientations of the same simulation.

Figure 4: A direct comparison between experimental and simulated weak beam
dark field images of irradiated tungsten, shown in inverted contrast with same
scale bar and ${\bf g}$=[200] imaging vector.

Comparing these image shows the same qualitative behaviour in experiment and
simulation. At low dose (0.01 dpa) the image is one of isolated black dots,
apparently homogeneously distributed. At intermediate dose (0.1 dpa), the dots
show emergent spatial correlation, appearing clustered rather than
homogeneous. At high dose (1.0 dpa), the spatial arrangement is of strings of
dots. With the ground-truth dislocation microstructure, we therefore identify
the low dose as dislocation loops, the intermediate as complex dislocation
networks, and the high dose as curved sections of dislocation lines. This is
the same conclusion as that drawn by Wang et al. [66] looking at dynamical
two-beam images.

We note that the spot size in our simulations is smaller than seen in
experiment. As noted above, this is partly potential dependent, though it may
also be related the fact that loops at room temperature will have non-zero
mobility, even when self-trapped by elastic interactions [27], and so be able
to coalesce to some extent at timescales beyond those accessible with MD.

A slightly better conventional TEM image can be generated by using the
precession method of Haley et al. [10]. This is a microscopy technique similar
in intent to the convergent weak-beam method of Prokhodtseva et al. [18], in
that both methods sample a small range of deviation parameters $s_{g}$ by
rocking the incident beam angle over a small angle. Where Prokhodtseva used a
linear tilt, Haley precesses the beam in a circle. Both suggest averaging over
order 10 images with the angle chosen so that the diffracted spot remains
within the objective aperture. In simulation we use a tilt angle of 5 mrad and
10 images. Taking more images does not improve the result to the naked eye.
The result of this procedure is shown in figure 5. Two ${\bf g}$-vectors are
shown here, both in weak beam conditions, $s_{\bf g}\sim 0.2\rm{nm}^{-1}$. For
both imaging vectors the contrast oscillations due to variations in strain in
the background and the beading seen in dislocation lines, exaggerated by high
deviation parameter, are smoothed out by averaging over several images with
varying $s_{\bf g}$ [68]. This has the effect of joining the black spots into
more clear lines. However the $\bf{g}\cdot\bf{b}$ invisibility criterion [4]
is still clearly in effect, with approximately half the dislocation network
visible at each ${\bf g}$-vector.

Figure 5: A snapshot of a massively overlapping cascade simulations at 0.10
dpa, with zone axis close to the [001] direction. Simulated TEM with g-vectors
${\mathbf{g}}=[110]$, $n_{\bf g}$ = 6.25 in images a) and b) and
${\mathbf{g}}=[200]$, $n_{\bf g}$ = 4.0 in c) and d). Images b) and d) use the
precession technique. Images are shown in reverse contrast. Right: dislocation
network coloured by $\bf{g}\cdot\bf{b}$ ( red $\bf{g}\cdot\bf{b}=+1$, white
$\bf{g}\cdot\bf{b}=0$, blue $\bf{g}\cdot\bf{b}=-1$ ).

This observation leads us to conclude that it will be difficult to establish
the true extent of a complex dislocation network using TEM, unless care is
taken to superimpose the images from a range of $\bf{g}$-vectors. A dramatic
example of this can be seen in Ref [6], where elastic energy minimisation
leads to the majority of large loops having the same Burgers vector at high
dose.

##  VI Conclusions

In this paper we have described a simple transformation of the Howie-Whelan
equations which enables their quick and efficient evaluation over arbitrarily
complex atomic configurations without the need to find displacement or strain
fields. This is done by constructing the complex phase
$x({\bf{r}})=\exp(i\bf{g}\cdot\bf{r})$ at atom sites, and interpolating to
form a continuous field.

Partly for didactic purposes, we showed how established results in expert TEM
lab use transfer to simulated TEM image generation. We showed how changing the
deviation parameter highlights different parts of the strain field, and how
each images produced show a strong spatial correlation to the ground truth
dislocation structure, and how our simulated microstructure is a good
qualitative fit to both dynamical and weak beam images in the literature. We
conclude that the generation and evolution of our low-temperature high dose
microstructure in simulation is a match to the description of network
formation in Wang et al. [66], and does not require the diffusion of mobile
loops. Importantly here we were able to demonstrate unambiguously, by
comparing to ground-truth dislocation microstructure, how the network of
dislocations produced by high dose irradiation is almost invisible in weak-
beam dark-field imaging due to a combination of the suppression of the image
intensity of small strain fields and $\bf{g}\cdot\bf{b}$ invisibility.

We conclude that attempts to model the irradiated microstructure with simple
dislocation loop objects, such as used in object kinetic Monte Carlo [69, 70,
71] or cluster dynamics [72], must fail at the high-dose limit where the
formation of a complex dislocation network is driven by fluctuating elastic
stresses rather than diffusing point-like objects. Conversely, attempts to
characterize high-dose microstructure in terms of a size-frequency
distribution of dislocation loops will also be misleading, without first
demonstrating that a complex dislocation network has not formed.

##  VII Acknowledgements

The authors would like to thank Hongbing Yu and Chris Grovenor for helpful
insights and Andrew Warwick for code testing.

This work has been carried out within the framework of the EUROfusion
Consortium, funded by the European Union via the Euratom Research and Training
Programme (Grant Agreement No. 101052200 - EUROfusion), and by the RCUK Energy
Programme, Grant No. EP/W006839/1. To obtain further information on the data
and models underlying the paper please contact<EMAIL_ADDRESS>The views and opinions expressed herein do not necessarily reflect those of
the European Commission. The authors acknowledge the use of the Cambridge
Service for Data Driven Discovery (CSD3) and associated support services
provided by the University of Cambridge Research Computing Services [73] that
assisted the completion of this study.

##  VIII Data availability

_The code used to generate the dynamical two-beam images will be made public
on acceptance of this paper._

##  IX Author Contributions

S.L.D. and D.R.M. developed the concept and the equations used to solve for
simulated TEM images. D.R.M. implemented the code and ran the image
simulations. M.B. ran the high dose microstructure simulations. J.H. and E.P.
ensured the simulation methodology matched that of expert TEM use. F.H and
G.H. provided high dose experimental TEM images.

## References

  * You  _et al._ [2016] J. H. You, E. Visca, C. Bachmann, T. Barrett, F. Crescenzi, M. Fursdon, H. Greuner, D. Guilhem, P. Languille, M. Li, S. McIntosh, A. V. Müller, J. Reiser, M. Richou, and M. Rieth, European DEMO divertor target: Operational requirements and material-design interface, Nuclear Materials and Energy 9, 171 (2016).
  * Hasegawa  _et al._ [2011] A. Hasegawa, T. Tanno, S. Nogami, and M. Satou, Property change mechanism in tungsten under neutron irradiation in various reactors, Journal of Nuclear Materials 417, 491 (2011), proceedings of ICFRM-14.
  * Kirk  _et al._ [2015] M. A. Kirk, X. Yi, and M. L. Jenkins, Characterization of irradiation defect structures and densities by transmission electron microscopy, Journal of Materials Research 30, 1195–1201 (2015).
  * Jenkins and Kirk [2001] M. L. Jenkins and M. A. Kirk,  _Characterization of Radiation Damage by Transmission Electron Microscopy_ , Series in Microscopy in Materials Science (IOP, Bristol, 2001).
  * Allen  _et al._ [2006] T. Allen, J. Cole, C. Trybus, D. Porter, H. Tsai, F. Garner, E. Kenik, T. Yoshitake, and J. Ohta, The effect of dose rate on the response of austenitic stainless steels to neutron radiation, Journal of Nuclear Materials 348, 148 (2006).
  * El-Atwani  _et al._ [2018a] O. El-Atwani, E. Aydogan, E. Esquivel, M. Efe, Y. Wang, and S. Maloy, Detailed transmission electron microscopy study on the mechanism of dislocation loop rafting in tungsten, Acta Materialia 147, 277 (2018a).
  * El-Atwani  _et al._ [2014] O. El-Atwani, J. Hinks, G. Greaves, S. Gonderman, T. Qiu, M. Efe, and J. Allain, In-situ TEM observation of the response of ultrafine- and nanocrystalline-grained tungsten to extreme irradiation environments, Scientific Reports 4, 4716 (2014).
  * El-Atwani  _et al._ [2015] O. El-Atwani, K. Hattar, J. Hinks, G. Greaves, S. Harilal, and A. Hassanein, Helium bubble formation in ultrafine and nanocrystalline tungsten under different extreme conditions, Journal of Nuclear Materials 458, 216 (2015).
  * Yao  _et al._ [2013] B. Yao, D. Edwards, and R. Kurtz, TEM characterization of dislocation loops in irradiated bcc Fe-based steels, Journal of Nuclear Materials 434, 402 (2013).
  * Haley  _et al._ [2017] J. C. Haley, S. A. Briggs, P. D. Edmondson, K. Sridharan, S. G. Roberts, S. Lozano-Perez, and K. G. Field, Dislocation loop evolution during in-situ ion irradiation of model FeCrAl alloys, Acta Materialia 136, 390 (2017).
  * Hirsch  _et al._ [1956] P. B. Hirsch, R. W. Horne, and M. J. Whelan, Direct observations of the arrangement and motion of dislocations in aluminium, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 1, 677 (1956).
  * Wu and Schäublin [2019] W. Wu and R. Schäublin, TEM diffraction contrast images simulation of dislocations, Journal of Microscopy 275, 11 (2019).
  * Hirsch  _et al._ [1967] P. B. Hirsch, A. Howie, R. Nicholson, D. W. Pashley, and M. J. Whelan,  _Electron Microscopy of Thin Crystals_, 2nd ed. (Krieger Publishing Company, Malabar, Florida, 1967) Chap. 10.
  * Silcox and Hirsch [1959] J. Silcox and P. B. Hirsch, Direct observations of defects in quenched gold, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 4, 72 (1959).
  * Loretto  _et al._ [2015] M. Loretto, P. Phillips, and M. Mills, Stacking fault tetrahedra in metals, Scripta Materialia 94, 1 (2015).
  * Baluc  _et al._ [2011] N. Baluc, J. L. Boutard, S. L. Dudarev, M. Rieth, J. Brito Correia, B. Fournier, J. Henry, F. Legendre, T. Leguey, M. Lewandowska, R. Lindau, E. Marquis, A. Muñoz, B. Radiguet, and Z. Oksiuta, Review on the EFDA work programme on nano-structured ODS RAF steels, Journal of Nuclear Materials 417, 149 (2011).
  * Yao  _et al._ [2008] Z. Yao, M. Hernández-Mayoral, M. L. Jenkins, and M. A. Kirk, Heavy-ion irradiations of Fe and Fe–Cr model alloys Part 1: Damage evolution in thin-foils at lower doses, Philosophical Magazine 88, 2851 (2008).
  * Prokhodtseva  _et al._ [2013] A. Prokhodtseva, B. Décamps, A. Ramar, and R. Schäublin, Impact of He and Cr on defect accumulation in ion-irradiated ultrahigh-purity Fe(Cr) alloys, Acta Materialia 61, 6958 (2013).
  * Yi  _et al._ [2016] X. Yi, M. L. Jenkins, M. A. Kirk, Z. Zhou, and S. G. Roberts, In-situ TEM studies of 150 keV W+ ion irradiated W and W-alloys: Damage production and microstructural evolution, Acta Materialia 112, 105 (2016).
  * Sand  _et al._ [2013] A. E. Sand, S. L. Dudarev, and K. Nordlund, High-energy collision cascades in tungsten: Dislocation loops structure and clustering scaling laws, EPL (Europhysics Letters) 103, 46003 (2013).
  * Yi  _et al._ [2015] X. Yi, A. E. Sand, D. R. Mason, M. A. Kirk, S. G. Roberts, K. Nordlund, and S. L. Dudarev, Direct observation of size scaling and elastic interaction between nano-scale defects in collision cascades, EPL (Europhysics Letters) 110, 36001 (2015).
  * Mason  _et al._ [2019] D. Mason, D. Nguyen-Manh, M.-C. Marinica, R. Alexander, A. Sand, and S. L. Dudarev, Relaxation volumes of microscopic and mesoscopic irradiation-induced defects in tungsten, J. Appl. Phys 126, 075112 (2019).
  * Derlet  _et al._ [2011] P. M. Derlet, M. R. Gilbert, and S. L. Dudarev, Simulating dislocation loop internal dynamics and collective diffusion using stochastic differential equations, Phys. Rev. B 84, 134109 (2011).
  * Wu  _et al._ [2012] W. Wu, R. Schäublin, and J. Chen, General dislocation image stress of anisotropic cubic thin film, Journal of Applied Physics 112, 093522 (2012).
  * Arakawa  _et al._ [2007] K. Arakawa, K. Ono, M. Isshiki, K. Mimura, M. Uchikoshi, and H. Mori, Observation of the one-dimensional diffusion of nanometer-sized dislocation loops, Science 318, 956 (2007).
  * Castin  _et al._ [2019] N. Castin, A. Dubinko, G. Bonny, A. Bakaev, J. Likonen, A. De Backer, A. Sand, K. Heinola, and D. Terentyev, The influence of carbon impurities on the formation of loops in tungsten irradiated with self-ions, Journal of Nuclear Materials 527, 151808 (2019).
  * Mason  _et al._ [2014] D. R. Mason, X. Yi, M. A. Kirk, and S. L. Dudarev, Elastic trapping of dislocation loops in cascades in ion-irradiated tungsten foils, Journal of Physics: Condensed Matter 26, 375701 (2014).
  * Sand  _et al._ [2018] A. Sand, J. Byggmästar, A. Zitting, and K. Nordlund, Defect structures and statistics in overlapping cascade damage in fusion-relevant bcc metals, Journal of Nuclear Materials 511, 64 (2018), special Section on ”18th International Conference on Fusion Reactor Materials”.
  * Byggmästar  _et al._ [2019] J. Byggmästar, F. Granberg, A. E. Sand, A. Pirttikoski, R. Alexander, M.-C. Marinica, and K. Nordlund, Collision cascades overlapping with self-interstitial defect clusters in Fe and W, Journal of Physics: Condensed Matter 31, 245402 (2019).
  * Derlet and Dudarev [2020] P. M. Derlet and S. L. Dudarev, Microscopic structure of a heavily irradiated material, Phys. Rev. Materials 4, 023605 (2020).
  * Mason  _et al._ [2020] D. R. Mason, S. Das, P. M. Derlet, S. L. Dudarev, A. J. London, H. Yu, N. W. Phillips, D. Yang, K. Mizohata, R. Xu, and F. Hofmann, Observation of transient and asymptotic driven structural states of tungsten exposed to radiation, Phys. Rev. Lett. 125, 225503 (2020).
  * Warwick  _et al._ [2021] A. R. Warwick, M. Boleininger, and S. L. Dudarev, Microstructural complexity and dimensional changes in heavily irradiated zirconium, Phys. Rev. Materials 5, 113604 (2021).
  * Head  _et al._ [1973] A. Head, P. Humble, L. Clarebrough, A. Morton, and C. Forwood,  _Computed Electron Micrographs and Defect Identification_ (North-Holland Publishing Company, Amsterdam, 1973).
  * Schäublin  _et al._ [2002] R. Schäublin, M.-J. Caturla, M. Wall, T. Felter, M. Fluss, B. Wirth, T. Diaz de la Rubia, and M. Victoria, Correlating TEM images of damage in irradiated materials to molecular dynamics simulations, Journal of Nuclear Materials 307-311, 988 (2002).
  * Madsen and Susi [2020] J. Madsen and T. Susi, abTEM: Ab initio transmission electron microscopy image simulation, Microscopy and Microanalysis 26, 448 (2020).
  * Rangel DaCosta  _et al._ [2021] L. Rangel DaCosta, H. G. Brown, P. M. Pelz, A. Rakowski, N. Barber, P. O’Donovan, P. McBean, L. Jones, J. Ciston, M. Scott, and C. Ophus, Prismatic 2.0 – simulation software for scanning and high resolution transmission electron microscopy (STEM and HRTEM), Micron 151, 103141 (2021).
  * Mason  _et al._ [2021a] D. R. Mason, F. Granberg, M. Boleininger, T. Schwarz-Selinger, K. Nordlund, and S. L. Dudarev, Parameter-free quantitative simulation of high-dose microstructure and hydrogen retention in ion-irradiated tungsten, Phys. Rev. Materials 5, 095403 (2021a).
  * Granberg  _et al._ [2016] F. Granberg, K. Nordlund, M. W. Ullah, K. Jin, C. Lu, H. Bei, L. M. Wang, F. Djurabekova, W. J. Weber, and Y. Zhang, Mechanism of radiation damage reduction in equiatomic multicomponent single phase alloys, Phys. Rev. Lett. 116, 135504 (2016).
  * Velişa  _et al._ [2017] G. Velişa, M. W. Ullah, H. Xue, K. Jin, M. L. Crespillo, H. Bei, W. J. Weber, and Y. Zhang, Irradiation-induced damage evolution in concentrated Ni-based alloys, Acta Mater. 135, 54 (2017).
  * Granberg  _et al._ [2020] F. Granberg, J. Byggmästar, and K. Nordlund, Defect accumulation and evolution during prolonged irradiation of Fe and FeCr alloys, J. Nucl. Mater. 528, 151843 (2020).
  * Mason  _et al._ [2021b] D. R. Mason, A. Reza, F. Granberg, and F. Hofmann, Estimate for thermal diffusivity in highly irradiated tungsten using molecular dynamics simulation, Phys. Rev. Materials 5, 125407 (2021b).
  * Peng  _et al._ [2004] L.-M. Peng, S. L. Dudarev, and M. J. Whelan,  _High Energy Electron Diffraction and Microscopy_ (Oxford University Press, 2004).
  * Dudarev [2003] S. L. Dudarev, Coherent motion of interstitial defects in a crystalline material, Philosophical Magazine 83, 3577–3597 (2003).
  * Boleininger  _et al._ [2018] M. Boleininger, T. D. Swinburne, and S. L. Dudarev, Atomistic-to-continuum description of edge dislocation core: Unification of the Peierls-Nabarro model with linear elasticity, Physical Review Materials 2, 083803 (2018).
  * Boleininger and Dudarev [2019] M. Boleininger and S. L. Dudarev, Continuum model for the core of a straight mixed dislocation, Physical Review Materials 3, 093801 (2019).
  * Arakawa and Short [2020] K. Arakawa and M. P. Short, Experimental validation of models: In situ TEM for radiation damage, in  _Handbook of Materials Modelling_, edited by W. Andreoni and S. Yip (Springer, 2020) pp. 2503–2516.
  * Howie and Whelan [1961] A. Howie and M. J. Whelan, Diffraction contrast of electron microscope images of crystal lattice defects. II. The development of a dynamical theory, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 263, 217 (1961).
  * Roters  _et al._ [2010] F. Roters, P. Eisenlohr, L. Hantcherli, D. D. Tjahjanto, T. R. Bieler, and D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications, Acta Materialia 58, 1152 (2010).
  * Dudarev  _et al._ [1995] S. Dudarev, L.-M. Peng, and M. Whelan, On the Doyle-Turner representation of the optical potential for RHEED calculations, Surface Science 330, 86 (1995).
  * Peng  _et al._ [1996a] L. Peng, G. Ren, S. L. Dudarev, and M. J. Whelan, Robust parameterization of elastic and absorptive electron atomic scattering factors, Acta Crystallographica Section A 52, 257 (1996a).
  * Peng  _et al._ [1996b] L. Peng, G. Ren, S. L. Dudarev, and M. J. Whelan, Debye–waller factors and absorptive scattering factors of elemental crystals, Acta Crystallographica Section A 52, 456 (1996b).
  * Zhou  _et al._ [2006] Z. Zhou, M. L. Jenkins, S. L. Dudarev, A. P. Sutton, and M. A. Kirk, Simulations of weak-beam diffraction contrast images of dislocation loops by the many-beam Howie–Basinski equations, Philosophical Magazine 86, 4851 (2006).
  * Byggmästar  _et al._ [2018] J. Byggmästar, F. Granberg, and K. Nordlund, Effects of the short-range repulsive potential on cascade damage in iron, J. Nucl. Mater. 508, 530 (2018).
  * Granberg  _et al._ [2021] F. Granberg, J. Byggmästar, and K. Nordlund, Molecular dynamics simulations of high-dose damage production and defect evolution in tungsten, Journal of Nuclear Materials 556, 153158 (2021).
  * Zhang  _et al._ [2017] S. Zhang, K. Nordlund, F. Djurabekova, F. Granberg, Y. Zhang, and T. S. Wang, Radiation damage buildup by athermal defect reactions in nickel and concentrated nickel alloys, Mater. Res. Lett. 5, 433 (2017).
  * Markelj  _et al._ [2023] S. Markelj, X. Jin, F. Djurabekova, J. Zavašnik, E. Punzón-Quijorna, T. Schwarz-Selinger, M. L. Crespillo, G. G. López, F. Granberg, E. Lu, K. Nordlund, A. Šestan, and M. Kelemen, Unveiling the radiation-induced defect production and damage evolution in tungsten using multi-energy rutherford backscattering spectroscopy in channeling configuration, Acta Materialia , 119499 (2023).
  * Boleininger  _et al._ [2022] M. Boleininger, S. L. Dudarev, D. R. Mason, and E. Martínez, Volume of a dislocation network, Phys. Rev. Mater. 6, 063601 (2022).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117, 1 (1995).
  * Norgett  _et al._ [1975] M. J. Norgett, M. T. Robinson, and I. M. Torrens, A proposed method of calculating displacement dose rates, Nuclear Engineering and Design 33, 50 (1975).
  * Fikar  _et al._ [2017] J. Fikar, R. Gröger, and R. Schäublin, Interaction of irradiation-induced prismatic dislocation loops with free surfaces in tungsten, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 393, 186 (2017), Computer Simulation of Radiation effects in Solids Proceedings of the 13 COSIRES Loughborough, UK, June 19-24 2016\.
  * Zheng and Han [2020] R.-Y. Zheng and W.-Z. Han, Comparative study of radiation defects in ion irradiated bulk and thin-foil tungsten, Acta Materialia 186, 162 (2020).
  * Mason  _et al._ [2017] D. R. Mason, D. Nguyen-Manh, and C. S. Becquart, An empirical potential for simulating vacancy clusters in tungsten, Journal of Physics: Condensed Matter 29, 505501 (2017).
  * He  _et al._ [2023] G. He, H. Yu, P. Karamched, J. Liu, and F. Hofmann, Elastic strain associated with irradiation-induced defects in self-ion irradiated tungsten, Scripta Materialia 237, 115687 (2023).
  * El-Atwani  _et al._ [2018b] O. El-Atwani, E. Esquivel, M. Efe, E. Aydogan, Y. Wang, E. Martinez, and S. Maloy, Loop and void damage during heavy ion irradiation on nanocrystalline and coarse grained tungsten: Microstructure, effect of dpa rate, temperature, and grain size, Acta Materialia 149, 206 (2018b).
  * Dudarev  _et al._ [2014] S. L. Dudarev, K. Arakawa, X. Yi, Z. Yao, M. L. Jenkins, M. R. Gilbert, and P. M. Derlet, Spatial ordering of nano-dislocation loops in ion-irradiated materials, Journal of Nuclear Materials 455, 16 (2014), Proceedings of the 16th International Conference on Fusion Reactor Materials (ICFRM-16).
  * Wang  _et al._ [2023] S. Wang, W. Guo, T. Schwarz-Selinger, Y. Yuan, L. Ge, L. Cheng, X. Zhang, X. Cao, E. Fu, and G.-H. Lu, Dynamic equilibrium of displacement damage defects in heavy-ion irradiated tungsten, Acta Materialia 244, 118578 (2023).
  * Stukowski [2009] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Cockayne [1973] D. J. H. Cockayne, The principles and practice of the weak-beam method of electron microscopy, J. Microscopy 98, 116 (1973).
  * Domain  _et al._ [2004] C. Domain, C. Becquart, and L. Malerba, Simulation of radiation damage in Fe alloys: an object kinetic Monte Carlo approach, Journal of Nuclear Materials 335, 121 (2004).
  * Becquart and Domain [2009] C. Becquart and C. Domain, An object kinetic Monte Carlo simulation of the dynamics of helium and point defects in tungsten, Journal of Nuclear Materials 385, 223 (2009), Nuclear Materials {IIIProceedings} of the E-MRS 2008 Spring Meeting: Third Symposium N on Nuclear Materials.
  * Martin-Bragado  _et al._ [2013] I. Martin-Bragado, A. Rivera, G. Valles, J. L. Gomez-Selles, and M. J. Caturla, MMonCa: An object kinetic Monte Carlo simulator for damage irradiation evolution and defect diffusion, Computer Physics Communications 184, 2703 (2013).
  * Marian and Bulatov [2011] J. Marian and V. V. Bulatov, Stochastic cluster dynamics method for simulations of multispecies irradiation damage accumulation, Journal of Nuclear Materials 415, 84 (2011).
  * [73] www.csd3.cam.ac.uk.
  * Howie and Basinski [1968] A. Howie and Z. S. Basinski, Approximations of the dynamical theory of diffraction contrast, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 17, 1039 (1968).
  * Shoemake [1985] K. Shoemake, Animating rotation with quaternion curves, SIGGRAPH85 19, 245 (1985).

##  X Supplementary Material

###  X.1 Validation by comparison to previous work

We validate the formulation of the Howie-Whelan equations based on atomic
positions only (equations 9) by comparing to the results from an existing TEM
simulation code, TEMACI[52] based on analytic formulae for the displacement
field in the elastic limit. We consider four typical irradiation-induced
defects which image brightly in weak-beam dark-field conditions, and four
which are faint. The bright defects are prismatic interstitial loops with
$1/2[111]$ and $[100]$ Burgers vectors, with diameters of 2 nm and 6 nm,
imaged with $\mathbf{g}=2\pi/a_{0}[200]$, and using $s_{\bf
g}=0.18$nm${}^{-1},\xi_{0}=13.4,\xi_{\bf g}=36.6$ nm. The defects were placed
at a depth 25 nm in a 50 nm thick foil, and relaxed with fixed zero strain
boundary conditions using LAMMPS[58] using the MNB empirical tungsten
potential[62]. The simulation box was $80\times 80\times 160$ conventional
cubic unit cells (2 M atoms).

The faintly imaging defects were prismatic interstitial loops with $[010]$
Burgers vectors, so that $\mathbf{g}\cdot\mathbf{b}=0$, and voids. Again these
had diameters 2 nm and 6 nm and relaxed as above. The bright defect simulated
TEM images are shown in figure 6, and the faint in figure 7. For the bright
defects the saturated white level is set to $|\phi_{\bf g}|^{2}=0.15$, for the
faint defects the white level is at $|\phi_{\bf g}|^{2}=0.0625$.

The TEMACI images technically solve the related Howie-Basinski equations[74].
In the columnar approximation, and with two beams only, the Howie-Basinski
equations reduce to the Howie-Whelan equations that we use here. The elastic
strain fields used were computed analytically by approximating the loops as a
hexagon of dislocation line segments.

a)  b)
c)  d)
e)  f)
g)  h)

Figure 6: A comparison between simulated TEM images produced from atomistic
data (a,c,e,g) and elastic strain (b,d,f,h). Top to bottom: 2 nm 1/2 [111]
loop, 6 nm 1/2 [111] loop, 2 nm [100] loop, 6 nm [100] loop at a depth 25 nm
in a 50 nm foil. Weak beam dark field imaging conditions $s_{g}=0.1777$ nm-1,
equivalent to $(\mathbf{g},4g)$ with $\mathbf{g}=[200]$. All images have equal
intensity range.

a)  b)
c)  d)
e)  f)
g)  h)

Figure 7: A comparison between simulated TEM images produced from atomistic
data (a,c,e,g) and elastic strain (b,d,f,h). Top to bottom: 2 nm [010] loop, 6
nm [010] loop, 2 nm void, 6 nm void at a depth 25 nm in a 50 nm foil. Weak
beam dark field imaging conditions $s_{g}=0.1777$ nm-1, equivalent to
$(\mathbf{g},4g)$ with $\mathbf{g}=[200]$. All images have equal intensity
range, with maximum intensity level set to one quarter that of figure 6. Note
the relaxed atomistic simulation cell has double the linear dimension of the
image, to minimise periodic boundary effects.

It is clear by eye that the two sets of images are very similar, clearly
showing the same extent and symmetries. To make a more systematic comparison,
we find the correlation between the images. We can find the similarity between
two images, whose intensities at 2d position $\mathbf{x}$ are given by
$f(\mathbf{x})$ and $g(\mathbf{x})$, using a simple correlation function, $c$,
defined by

$c=\frac{\int\left(f(\mathbf{x})-\bar{f}\right)\left(g(\mathbf{x})-\bar{g}\right)\mathrm{d}^{2}\mathbf{x}}{\sqrt{\int\left(f(\mathbf{x})-\bar{f}\right)^{2}\mathrm{d}^{2}\mathbf{x}\,\int\left(g(\mathbf{x})-\bar{g}\right)^{2}\mathrm{d}^{2}\mathbf{x}}},$
(10)

where $\bar{f},\bar{g}$ are the average intensities in each image
respectively. This correlation function is invariant to affine transformations
of the intensity ( brightness and contrast shifts ) and weakly dependent on
noise providing the signal-to-noise ratio is large. A correlation value of one
means that two images have identical (relative) spatially varying intensities,
a value of minus one means they form a pair of negatives.

Table 1 shows the correlation factors between the eight images. If we read
across each row looking for the closest match between an atomistic TEM image
and a strain field TEM image, we find a very highly correlated like-for-like
match in most cases, with the exception being that the atomistic image for the
‘invisible’ $\mathbf{b}=[010]$ loops could still be mistaken for a small
$\mathbf{b}=[100]$ loop.

| TEMACI strain formulae
---|---
atomic |  $[111]$(2nm) |  $[111]$(6nm) |  $[100]$(2nm) |  $[100]$(6nm) |  $[010]$(2nm) |  $[010]$(6nm) | void (2nm) | void (6nm)
$[111]$(2nm) | 0.91 | 0.26 | 0.83 | 0.55 | 0.51 | 0.30 | 0.71 | 0.23
$[111]$(6nm) | 0.13 | 0.94 | 0.40 | 0.58 | 0.53 | 0.15 | 0.00 | 0.69
$[100]$(2nm) | 0.72 | 0.46 | 0.90 | 0.74 | 0.71 | 0.34 | 0.72 | 0.00
$[100]$(6nm) | 0.51 | 0.57 | 0.82 | 0.94 | 0.75 | 0.38 | 0.17 | 0.62
$[010]$(2nm) | 0.32 | 0.37 | 0.74 | 0.65 | 0.74 | 0.24 | 0.04 | 0.49
$[010]$(6nm) | 0.25 | 0.21 | 0.54 | 0.52 | 0.42 | 0.50 | 0.04 | 0.31
void (2nm) | 0.87 | 0.04 | 0.35 | 0.06 | 0.18 | 0.04 | 0.92 | -0.05
void (6nm) | 0.07 | 0.66 | 0.22 | 0.77 | 0.30 | 0.45 | -0.03 | 0.97
Table 1: Image correlation factors between simulated TEM images generated
using atomic position data (figures 6,7 a,c,e,g) and strain formulae using
TEMACI (same figures b,d,f,h). For each row ( each atomic figure ) the best
correlated TEMACI figure is highlighted.

###  X.2 Finding good two-beam diffraction conditions

In this section we derive equations required for finding good dynamical two
beam conditions. Start with an incident high energy electron beam with wave
vector given by

$\mathbf{k}=\frac{mv}{\hbar}\hat{\mathbf{k}},$ (11)

where $v$ is the velocity of the electrons, and $\hat{\mathbf{k}}$ a unit
vector in the direction of the beam.

We have a set of atoms in general position defining the foil. For good
diffraction contrast, we hope most of the atoms are in a lattice, so we encode
the conventional unit cell of the atomic lattice as the matrix $\mathbf{C}$,
where the component $C_{ij}$ is the $i^{th}$ Cartesian component of the
$j^{th}$ lattice vector in lab frame. For a cubic unit cell, oriented with a
zone axis along $[001]$, this matrix might be simply
$\mathbf{C}=a_{0}\mathbf{I}$, where $\mathbf{I}$ is the $3\times 3$ identity
matrix, but more generally the columns of $\mathbf{C}$ are the three
conventional lattice vectors $\mathbf{c}_{1},\mathbf{c}_{2}$ and
$\mathbf{c}_{3}$, ie

$\mathbf{C}=\left(\begin{array}[]{ccc}\mathbf{c}_{1}&\mathbf{c}_{2}&\mathbf{c}_{3}\\\
|&|&|\\\ |&|&|\end{array}\right).$ (12)

$\mathbf{C}$ will be determined by the positions of the atoms in the input
file, which may have any orientation or strains present, and not by the
bounding simulation box.

Note that this equates the coordinate frame of the input file with the lab
frame, while the crystal frame is allowed to vary.

The reciprocal lattice vectors by convention are given by
$\mathbf{b}_{1}=2\pi\,\mathbf{c}_{2}\times\mathbf{c}_{3}/(\mathbf{c}_{1}\cdot(\mathbf{c}_{2}\times\mathbf{c}_{3}))$
etc, which can be encoded in the matrix $\mathbf{B}=2\pi\mathbf{C}^{-1}$, so
that the $[hkl]$ reciprocal lattice vector is a vector in the crystal frame
given by

$\mathbf{g}_{hkl}=\mathbf{B}^{T}\left(\begin{array}[]{c}h\\\ k\\\
l\end{array}\right).$ (13)

If we solve the dynamical two beam (Howie-Whelan) equations using the imaging
reciprocal vector $\mathbf{g}$, we can define the (energy) deviation parameter

$\epsilon_{\mathbf{g}}=\frac{\hbar^{2}({\mathbf{k}}+{\mathbf{g}})^{2}}{2m}-\frac{\hbar^{2}{\mathbf{k}}^{2}}{2m}=\frac{\hbar^{2}(2\mathbf{k}+\mathbf{g})\cdot\mathbf{g}}{2m}.$
(14)

Combining equations 11,13,14, we find the Bragg condition is met at the
$[hkl]$ reciprocal lattice vector when

$\frac{2mv}{\hbar}\,\hat{\mathbf{k}}\cdot\mathbf{B}^{T}\left(\begin{array}[]{c}h\\\
k\\\ l\end{array}\right)+\left|\mathbf{B}^{T}\left(\begin{array}[]{c}h\\\ k\\\
l\end{array}\right)\right|^{2}=0.$ (15)

This condition is trivially met for $[hkl]=[000]$, and can be solved for other
reciprocal vectors by rotating the incident beam vector $\hat{\mathbf{k}}$. We
will look for a rotation where this condition is met for a particular $[hkl]$,
but not met for other reciprocal lattice vectors
$[h^{\prime}k^{\prime}l^{\prime}]$.

We start by defining an orthonormal basis set
$\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}}$. The user selects a zone
axis direction, $[uvw]$, with a direction in the lab frame given by
$\hat{\mathbf{z}}\sim\mathbf{C}[uvw]^{T}$. We choose a second unit vector
$\hat{\mathbf{y}}$ to have no projection along the desired $g$-vector
direction. To align the crystal with the incident beam, we therefore first
rotate the crystal by $\mathbf{U}$ given by

$\mathbf{U}=\left(\begin{array}[]{ccc}\hat{\mathbf{x}}&\rule[2.15277pt]{10.76385pt}{0.5pt}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\hat{\mathbf{y}}&\rule[2.15277pt]{10.76385pt}{0.5pt}&\rule[2.15277pt]{10.76385pt}{0.5pt}\\\
\hat{\mathbf{z}}&\rule[2.15277pt]{10.76385pt}{0.5pt}&\rule[2.15277pt]{10.76385pt}{0.5pt}\end{array}\right),$
(16)

and then use a virtual ‘tilt stage’ to fine-tune the orientation of the
crystal through a rotation matrix $\mathbf{R}$ in order to get a good two-beam
diffraction pattern. Under these two rotations a vector $\mathbf{x}$ rotates
to $\mathbf{R}\mathbf{U}\mathbf{x}$. The $\mathbf{g}-$vector rotates to
$\mathbf{g}^{\prime}=\mathbf{R}\mathbf{U}\mathbf{g}$ so that for an atom at
position $\mathbf{x}$, we see
$\mathbf{g}^{\prime}\cdot\mathbf{x}^{\prime}=\mathbf{g}\cdot\mathbf{x}$, as
required.

If our tilt stage first rotates by $\theta$ about the x-axis, then by $\psi$
about the y-axis, the rotation matrix describing the transformation is

$\mathbf{R}=\left(\begin{array}[]{ccc}\cos\psi&\sin\theta\sin\psi&\cos\theta\sin\psi\\\
0&\cos\theta&-\sin\theta\\\
-\sin\psi&\sin\theta\cos\psi&\cos\theta\cos\psi\end{array}\right).$ (17)

The Bragg condition, equation 15, after rotation reads
$\hat{\mathbf{k}}\cdot\mathbf{g}^{\prime}=\alpha$, where
$\alpha=-\hbar|\mathbf{g}|^{2}/(2mv)$ is a constant, so if the components of
$\mathbf{g}$ are

$\left(\begin{array}[]{c}g_{1}\\\ g_{2}\\\
g_{3}\end{array}\right)=\mathbf{U}\mathbf{B}^{T}\left(\begin{array}[]{c}h\\\
k\\\ l\end{array}\right),$ (18)

this gives the vector equation

$\left(\begin{array}[]{c}-\sin\psi\\\ \cos\psi\sin\theta\\\
\cos\psi\cos\theta\end{array}\right)\cdot\left(\begin{array}[]{c}g_{1}\\\
g_{2}\\\ g_{3}\end{array}\right)=\alpha$ (19)

Equation 19 has solutions $\\{\psi\\}$ for each choice of $\theta$, given by

$\displaystyle\psi$ $\displaystyle=$ $\displaystyle\pm\rm{acos}\left[u\pm
v\right]+2n\pi\quad\quad(n\in\mathcal{Z}),$ $\displaystyle u$ $\displaystyle=$
$\displaystyle\frac{\alpha(g_{2}\sin\theta+g_{3}\cos\theta)}{g_{1}^{2}+(g_{2}\sin\theta+g_{3}\cos\theta)^{2}}$
$\displaystyle v$ $\displaystyle=$
$\displaystyle\frac{g_{1}\sqrt{g_{1}^{2}+(g_{2}\sin\theta+g_{3}\cos\theta)^{2}-\alpha^{2}}}{g_{1}^{2}+(g_{2}\sin\theta+g_{3}\cos\theta)^{2}}.$
(20)

We therefore need an algorithm to decide which pair $\\{\psi,\theta\\}$
represents best diffraction conditions.

For the perfect crystalline foil thickness $L$, the equations to solve for the
propagation of incident and diffracted electron beams (equation 9) reduce to:

${\partial\over\partial z}\left(\begin{array}[]{c}\phi_{0}(z)\\\
\phi_{\bf{g}}(z)\end{array}\right)=i\pi\left(\begin{array}[]{c
c}1/\xi_{0}&1/\xi_{\bf{g}}\\\
1/\xi_{\bf{g}}&1/\xi_{0}+2s_{\bf{g}}\end{array}\right)\left(\begin{array}[]{c}\phi_{0}(z)\\\
\phi_{\bf{g}}(z)\end{array}\right),$ (21)

This can be solved analytically with the boundary conditions
$\phi_{0}(0)=1,\phi_{\bf{g}}=0$, to give the solution at the back face of the
foil, $z=L$,

$\displaystyle\left|\phi_{0}(L)\right|^{2}$ $\displaystyle=$
$\displaystyle\frac{s_{\bf{g}}^{2}+\xi_{\bf{g}}^{-2}\cos^{2}\left[\pi\sqrt{s_{\bf{g}}^{2}+\xi_{\bf{g}}^{-2}}L\right]}{s_{\bf{g}}^{2}+\xi_{\bf{g}}^{-2}}$
$\displaystyle\left|\phi_{\bf{g}}(L)\right|^{2}$ $\displaystyle=$
$\displaystyle\frac{\xi_{\bf{g}}^{-2}\sin^{2}\left[\pi\sqrt{s_{\bf{g}}^{2}+\xi_{\bf{g}}^{-2}}L\right]}{s_{\bf{g}}^{2}+\xi_{\bf{g}}^{-2}}.$
(22)

The foil thickness $L$ is found in the imaging space, defined in section X.3.
Note that both $s_{\bf{g}}$ and $L$ depend on the virtual tilt rotation
$\mathbf{R}$,

We can now define a score for a rotation $\mathbf{R}$,

$S=\sum_{\mathbf{g}}\left|\phi_{\bf{g}}(L;\mathbf{R})\right|^{2},$ (23)

where the sum is over a range of reflections excluding $[000]$ and the
requested reflection $[hkl]$. The best rotation for dynamic two beam imaging
has the smallest value $S$, this corresponds to the rotation where unwanted
reflections have minimum intensity. For our score, we consider the set of
permitted reflections $[hkl]$ where $h^{2}+k^{2}+l^{2}\leq 10^{2}$.

If we set the $n_{g}^{\rm{th}}$ diffraction spot, ie $n_{g}\times[hkl]$
bright, and others faint, we then place the ‘aperture’ at $[hkl]$. (Noting
that in the Howie-Whelan approximation the aperture is a delta function. ).
Large $n_{g}$ corresponds to a large deviation parameter $s_{g}$. In this work
we have also used fractional values for $n_{g}$, which are more appropriate in
experimental work, as strongly-excited systematic reflections can cause
extinction effects in the image if integer values are chosen [4]. To find a
fractional value, we find the quaternions for the integer values of $n_{g}$,
and spherical interpolate between as appropriate [75].

###  X.3 Mapping atoms to the image space

In this section we consider how the atoms in the input file are placed into an
imaging space, whose length and breadth define the extend of the pixels in the
output image. This enables us to define the foil thickness.

Figure 8: Cartoon illustrating the geometry of the periodic input cell and the
imaging space. The input file contains atoms in the $[000]$ periodic replica,
illustrated by the dark shaded box. Periodic replicas inside the imaging space
( heavy black line box ) are kept.

We start by stating that the atoms are in a periodic supercell whose repeat
vectors, of length 100 nm or more, are given by the columns of the matrix
$\mathbf{A}$, where

$\mathbf{A}=\left(\begin{array}[]{ccc}\mathbf{a}_{1}&\mathbf{a}_{2}&\mathbf{a}_{3}\\\
|&|&|\\\ |&|&|\end{array}\right).$ (24)

If periodicity in one direction is not desired (to represent a foil), and it
is not present in the original file, then the user can add padding in that
direction. In figure 8, we show an input file with a free surface in one
direction, represented by the heavy shaded box in the centre.

After an affine transformation consisting of a rotation $\mathbf{U}$, and
translation $\bm{\delta}$, the position of an atom in the $[ijk]$ periodic
replica of $\mathbf{A}$ transforms to

$\mathbf{x}^{\prime}=\mathbf{U}\left(\mathbf{x}+\mathbf{A}\left(\begin{array}[]{c}i\\\
j\\\ k\end{array}\right)\right)+\bm{\delta}.$ (25)

The periodic replicas of the atoms in the input file after rotation and
translation are shown with light shaded boxes in figure 8. The atom replicas
are placed in in _imaging space_ , defined as a new orthorhombic box, whose
lattice vectors are given by the columns of the matrix $\mathbf{D}$, given by

$\mathbf{D}=\left(\begin{array}[]{ccc}d_{1}&0&0\\\ 0&d_{2}&0\\\
0&0&d_{3}\end{array}\right),$ (26)

where $d_{1}$ and $d_{2}$ are the size of the output image (in nm) and $d_{3}$
is the integration depth total, chosen by the user. The original periodic
boundaries across replicas of $\mathbf{A}$ will be respected within the
transformed cell, but there is no periodicity across replicas of $\mathbf{D}$–
rather this transformed cell should be viewed as a cuboidal crystallite whose
sides have lengths $d_{1},d_{2},d_{3}$, suspended in vacuum. The imaging space
is represented by the heavy lined box in figure 8, with its 3-axis aligned
with the electron beam.

The vector offset for the affine transformation in equation 25, $\bm{\delta}$,
is given by

$\bm{\delta}=\mathbf{D}\left(\begin{array}[]{c}{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\\\
{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\\\
{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\end{array}\right)-\mathbf{U}\mathbf{A}\left(\begin{array}[]{c}{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\\\
{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\\\
{}^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\end{array}\right),$ (27)

so that an atom in the centre of the first periodic replica of the input file
is transformed to the centre of the image.

Succinctly, the $[ijk]$ replica of an atom at $\mathbf{x}$ in the input file
is within the imaging space if for each Cartesian component
$\alpha\in\\{1,2,3\\}$,

$\mathrm{int}\left(\mathbf{U}\left(\mathbf{x}+\mathbf{A}\left[i-^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}},j-^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}},k-^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\right]^{T}\right)_{\alpha}/d_{\alpha}+^{\scriptstyle{1}}\\!\\!/\\!_{\scriptstyle{2}}\right)=0.$
(28)

In figure 8, some atoms from the central replica $[ijk]=[000]$, and its
neighbours $[\bar{1}00]$,$[100]$ are included in the imaging space.

The thickness of the foil, $L$, can then be determined from the positions of
the atoms in the imaging space, as shown in figure 8. Note that without the
periodic replicas, the edges of the image do not have the same thickness:
periodic replicas are required to see the full information contained in the
atomic input file. We find the thickness by conceptually dividing the atoms in
$\mathbf{D}$ into columns aligned with the 3-axis, and finding the maximum and
minimum height of atoms in that column. The foil thickness for the column is
given by the difference plus half a lattice parameter, ie $a_{0}/2$. The mean
foil thickness $L$ is the average thickness of the columns, and the mean
normal $\hat{\mathbf{n}}$ to the surface from the average spatial variation of
the maxima and of the minima. Under a small tilt stage rotation $\mathbf{R}$,
the foil thickness changes to

$L^{\prime}=\frac{[001]\hat{\mathbf{n}}}{[001]\mathbf{R}\hat{\mathbf{n}}}L.$
(29)