#  A differentiable forward model for the concurrent, multi-peak Bragg
coherent x-ray diffraction imaging problem

S. Maddali  Materials Science Division, Argonne National Laboratory Lemont, IL
60439 (USA)  T. D. Frazer  Materials Science Division, Argonne National
Laboratory Lemont, IL 60439 (USA)  N. Delegan  Materials Science Division,
Argonne National Laboratory Lemont, IL 60439 (USA)  Center for Molecular
Engineering, Argonne National Laboratory, Lemont, IL 60439 (USA)  K. J. Harmon
Materials Science Division, Argonne National Laboratory Lemont, IL 60439 (USA)
S. E. Sullivan  Materials Science Division, Argonne National Laboratory
Lemont, IL 60439 (USA)  Center for Molecular Engineering, Argonne National
Laboratory, Lemont, IL 60439 (USA)  M. Allain  Aix-Marseille Univ, CNRS,
Centrale Marseille, Institut Fresnel, Marseille, France  W. Cha  X-ray Science
Division, Argonne National Laboratory, Lemont, IL 60439 (USA)  A. Dibos
Center for Nanoscale Materials, Argonne National Laboratory, Lemont, IL 60439
(USA)  I. Poudyal  Materials Science Division, Argonne National Laboratory
Lemont, IL 60439 (USA)  S. Kandel  X-ray Science Division, Argonne National
Laboratory, Lemont, IL 60439 (USA)  Y. S. G. Nashed  Mathematics & Computer
Science Division, Argonne National Laboratory, Lemont, IL 60439 (USA)  F. J.
Heremans  Materials Science Division, Argonne National Laboratory Lemont, IL
60439 (USA)  Center for Molecular Engineering, Argonne National Laboratory,
Lemont, IL 60439 (USA)  Pritzker School of Molecular Engineering, University
of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637 (USA)  H. You
Materials Science Division, Argonne National Laboratory Lemont, IL 60439 (USA)
Y. Cao  Materials Science Division, Argonne National Laboratory Lemont, IL
60439 (USA)  Pritzker School of Molecular Engineering, University of Chicago,
5640 South Ellis Avenue, Chicago, IL 60637 (USA)  S. O. Hruszkewycz  Materials
Science Division, Argonne National Laboratory Lemont, IL 60439 (USA)
Corresponding author<EMAIL_ADDRESS>
###### Abstract

We present a general analytic approach to spatially resolve the nano-scale
lattice distortion field of strained and defected compact crystals with Bragg
coherent x-ray diffraction imaging (BCDI). Our approach relies on fitting a
differentiable forward model simultaneously to multiple BCDI datasets
corresponding to independent Bragg reflections from the same single crystal.
It is designed to be faithful to heterogeneities that potentially manifest as
phase discontinuities in the coherently diffracted wave, such as lattice
dislocations in an imperfect crystal. We retain fidelity to such small
features in the reconstruction process through a Fourier transform -based
resampling algorithm designed to largely avoid the point spread tendencies of
commonly employed interpolation methods. The reconstruction model defined in
this manner brings BCDI reconstruction into the scope of explicit optimization
driven by automatic differentiation. With results from simulations and
experimental diffraction data, we demonstrate significant improvement in the
final image quality compared to conventional phase retrieval, enabled by
explicitly coupling multiple BCDI datasets into the reconstruction loss
function.

##  1 Introduction

Bragg coherent diffraction imaging (BCDI) is a synchrotron-based
characterization technique that utilizes coherent x-ray illumination and 3D
phase retrieval algorithms to interrogate structural heterogeneities in
crystalline materials on the scale of tens of nanometers [1, 2, 3]. The Bragg
diffraction geometry is sensitive to the distortion field of the atomic
lattice, components of which are encoded in the complex-valued phase of the
scattered wave. Conventional BCDI seeks to spatially resolve individual
projections of the lattice distortion field by applying phase retrieval
methods to a single three-dimensional Bragg diffraction pattern. This
diffraction pattern corresponds to interfering x-rays reflected from a single
family of atomic planes within a crystal. The corresponding component of the
elastic strain tensor may be subsequently obtained by spatial differentiation
of the retrieved phase. This technique has greatly benefited research that
non-destructively interrogates crystalline materials in which properties and
behavior are influenced by lattice strain, such as charging cycle -induced
degradation of solid-state batteries [4, 5], catalytic reactions at metal
nanoparticle interfaces [6] and quantum sensors [7, 8]. Numerically, phase
retrieval typically involves fixed point iterative projections (FPIPs) acting
on an initial guess object. The solution is alternately constrained to an
assumed support mask in real space and the measured signal in Fourier space
[9, 10]. The initial support estimate is intermittently refined using a
‘shrinkwrap’ algorithm [11]. The iterations proceed until the solution is
deemed converged. A given combination of FPIP algorithms for a particular
reconstruction task, or ‘recipe’, is determined by trial-and-error and is
typically validated only in hindsight upon successful reconstruction. This
method of recipe design becomes increasingly difficult when reconstructing
defected crystals, for which convergence is far less stable, and sometimes
even elusive.

Recent efforts aim to directly retrieve the 3D lattice distortion field by
concurrently solving the phase retrieval problems for multiple independent
BCDI datasets and then fitting the lattice distortion to the reconstructed
phase profiles. Proof-of-concept demonstrations so far have used combinations
or variants of established FPIPs in their reconstruction algorithms [12, 13,
14, 15]. The majority of these are purely simulation-based with highly
idealized scattering crystals and distortion fields. This variety of multi-
reflection BCDI (MR-BCDI) has only recently been demonstrated with
experimental data [15].

In this paper we present an gradient descent -based solution for the general
MR-BCDI problem which offers the following major advantages over current
practices:

  1. 1.

Feature and data fidelity: In our reconstructions we maintain fidelity to
small features within the scattering crystals ($\sim$ 1 pixel), as well as to
the noise statistics in the diffraction data itself, by implementing a Fourier
transform -based re-gridding method in real space, and avoiding data re-
gridding altogether. Grid interpolation either in real or Fourier space is
typically necessary in BCDI to reconcile signal digitization with a chosen
reconstruction grid. Current methods that employ straightforward multi-linear
or polynomial interpolation are susceptible to point spread artifacts, which
our method greatly minimizes. The significance of our Fourier-based approach
is the ability to reconstruct crystals with higher defect and dislocation
densities than the usual FPIP approach.

  2. 2.

Flexible convergence: Stochastic gradient descent implemented over the whole
body of BCDI data affords us a flexible path to convergence to the desired
solution, in the manner of training of present-day machine learning models
[16, 17, 18]. In contrast, FPIP approaches using trial-and-error recipes
commit to an arbitrary rigid progression in solution space [19], punctuated by
uncontrolled ‘resets’ due to shrinkwrap. The stochastic optimization approach
is shown to be far more reliable for convergence than a typical trial-and-
error recipe.

  3. 3.

Reconstruction quality: We demonstrate an improvement in the reconstruction
quality enabled by explicit coupling of the BCDI data sets into the loss
function. This forward model approach is fundamentally different from
intermittent or _post hoc_ bootstrapping of the concurrently calculated phases
through the underlying lattice distortion [12, 15].

This paper is organized as follows. Section 2 describes the MR-BCDI loss
function over which the optimization is performed, similar to the well-known
Gaussian loss function used in phase retrieval [20]. It is formulated in terms
of the lattice distortion field suitably digitized to be consistent with each
of the BCDI signals. Section 3 demonstrates reconstruction results from
simulated as well as experimental data, the latter obtained from a specially
prepared sub-micron sized silicon carbide crystal for quantum information and
sensing applications. Section 3 also compares the results to conventional
phase retrieval reconstructions of the same data. Section 4 discusses the
place of our algorithm in the current landscape of coherent diffraction
imaging techniques and its impact on future BCDI experiments, along with an
assortment of issues like stability, support update, interpretability and
optimization constraints. Finally, the Methods section 5 details the various
Fourier transform-based operations employed in the multi-reflection forward
model, as well as specifics of the silicon carbide sample preparation.

##  2 The multi-reflection forward model

To describe the the distortion field we wish to reconstruct, we nominally
adopt the right-handed ortho-normal laboratory frame
$[\hat{\boldsymbol{s}}_{1}~{}\hat{\boldsymbol{s}}_{2}~{}\hat{\boldsymbol{s}}_{3}]$
commonly used at the Advanced Photon Source (Figure 1). The following
formalism can be extended generally to lab frames at other synchrotron
facilities. Here, $\hat{\boldsymbol{s}}_{3}$ is a unit vector pointing
downstream along the incident beam and $\hat{\boldsymbol{s}}_{2}$ is
vertically upward.

Figure 1:  BCDI setup at Beamline 34-ID-C of the Advanced Photon Source
showing the laboratory and detector frames
$[\hat{\boldsymbol{s}}_{1}~{}\hat{\boldsymbol{s}}_{2}~{}\hat{\boldsymbol{s}}_{3}]$
and
$[\hat{\boldsymbol{k}}_{1}~{}\hat{\boldsymbol{k}}_{2}~{}\hat{\boldsymbol{k}}_{3}]$
for a single Bragg condition.

In this frame, we define the following quantities:

  * •

The scalar field $\mathcal{A}:\mathbb{R}^{3}\rightarrow\mathbb{R}$ and the 3D
vector field $\boldsymbol{u}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$,
respectively denoting the local electron density and the vector displacements
of the atoms from their equilibrium positions, at location
$\boldsymbol{x}\in\mathbb{R}^{3}$. For optimization purposes we assume
$\mathcal{A}(\boldsymbol{x})\in[0,1]$ in arbitrary units.
$\mathcal{A}(\boldsymbol{x})=0$ and
$\boldsymbol{u}(\boldsymbol{x})=\boldsymbol{0}$ is enforced outside a user-
defined bounding box $\mathcal{V}$. $\mathcal{A}$ and $\boldsymbol{u}$ are
discretized on a user-defined grid within $\mathcal{V}$, given by a voxel size
$s_{0}$ along each of the Cartesian axes
$[\hat{\boldsymbol{s}}_{1}~{}\hat{\boldsymbol{s}}_{2}~{}\hat{\boldsymbol{s}}_{3}]$.

  * •

The reciprocal vectors of the crystal lattice corresponding to the BCDI data
sets ($M$ in number) are denoted by the set
$\\{\boldsymbol{G}_{i}\\}_{i=1}^{M}$. The physical units of the
$\boldsymbol{G}_{i}$ are chosen to be the inverse of those chosen for
$\boldsymbol{u}(\boldsymbol{x})$, _i.e._ ,
$\left\lVert\boldsymbol{G}_{i}\right\rVert=1/d_{i}$ (crystallographers’
convention, as opposed to the physicists’ convention of $2\pi/d_{i}$), where
$d_{i}$ is the corresponding atomic plane spacing.

  * •

The $3\times 3$ rotation matrices $\\{\mathcal{R}_{i}\\}_{i=1}^{M}$, to bring
the crystal lattice into the $\boldsymbol{G}_{i}$ Bragg condition (hereafter
referred to as the $i$’th _mounting position_).

  * •

The $3\times 3$ matrices
$\left\\{\boldsymbol{B}_{\text{real}}^{(i)}\right\\}_{i=1}^{M}$ whose columns
form a basis of 3D real-space sampling vectors, corresponding to the $i$’th
mounting position. These are determined by the diffraction geometry specific
to each Bragg reflection [21].

  * •

Dataset-dependent scaling factors
$\left\\{\chi_{i}\left|\chi_{i}\in\mathbb{R}\right.\right\\}_{i=1}^{M}$ for
the electron density $\mathcal{A}$ that controls total energy from each of the
$i$ diffracted waves that is incident upon the detector over the course of the
scan. We note that for a measured BCDI dataset, this total energy depends not
only on the structure factor of the scattering crystal but also the
acquisition time for each image in the scan.

These definitions lead to the following deductions:

  1. 1.

The complex-valued object wave consistent with the projection of
$\boldsymbol{u}(\boldsymbol{x})$ along $\boldsymbol{G}_{i}$, _prior to
rotating into the Bragg condition_ , is simply:
$\psi^{(i)}(\boldsymbol{x})=\mathcal{A}(\boldsymbol{x})\exp[\iota
2\pi\boldsymbol{G}_{i}^{T}\boldsymbol{u}(\boldsymbol{x})]$. After the crystal
is mounted into the $i$’th position, the _rotated_ version of this complex-
valued object, _i.e._ ,
$\psi^{(i)}(\mathcal{R}_{i}^{-1}\boldsymbol{x})\equiv\mathcal{A}(\mathcal{R}_{i}^{-1}\boldsymbol{x})\exp[\iota
2\pi\boldsymbol{G}_{i}^{T}\boldsymbol{u}(\mathcal{R}_{i}^{-1}\boldsymbol{x})]$,
represents the 3D object wave to be eventually propagated to the far field. We
note that the phase of $\psi^{(i)}(\mathcal{R}_{i}^{-1}\boldsymbol{x})$
already accounts for the rotation $\mathcal{R}_{i}$ being applied to the
crystal lattice and therefore to $\boldsymbol{G}_{i}$, in the following
manner:

$\displaystyle\boldsymbol{G}_{i}$ $\displaystyle\xrightarrow[]{\text{rotated
by }\mathcal{R}_{i}}\mathcal{R}_{i}\boldsymbol{G}_{i}$ (Transformation of
single vector) $\displaystyle\boldsymbol{u}(\boldsymbol{x})$
$\displaystyle\xrightarrow[]{\text{rotated by
}\mathcal{R}_{i}}\mathcal{R}_{i}\boldsymbol{u}\left(\mathcal{R}_{i}^{-1}\boldsymbol{x}\right)$
(Transformation of vector field)
$\displaystyle\therefore\angle\left[\psi^{(i)}\left(\mathcal{R}_{i}^{-1}\boldsymbol{x}\right)\right]$
$\displaystyle=2\pi\left(\boldsymbol{G}_{i}^{T}\mathcal{R}_{i}^{T}\right)\left[\mathcal{R}_{i}\boldsymbol{u}\left(\mathcal{R}_{i}^{-1}\boldsymbol{x}\right)\right]$
$\displaystyle=2\pi\boldsymbol{G}_{i}^{T}\boldsymbol{u}\left(\mathcal{R}_{i}^{-1}\boldsymbol{x}\right)$
(since $\mathcal{R}_{i}^{T}\mathcal{R}_{i}=\mathds{1}$, the identity matrix)

This last expression denotes the rotated phase field within the bounding box
$\mathcal{V}$.

  2. 2.

With the definition of $\boldsymbol{B}_{\text{real}}^{(i)}$ above, the $i$’th
_digital_ object wave may be indexed by an integer column vector
$\boldsymbol{m}\in\mathbb{Z}^{3}$:

$\displaystyle\psi_{\boldsymbol{m}}^{(i)}=\psi^{(i)}\left(\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{real}}^{(i)}\boldsymbol{m}\right)=A\left(\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{real}}^{(i)}\boldsymbol{m}\right)\exp\left[\iota
2\pi\boldsymbol{G}_{i}^{T}\boldsymbol{u}\left(\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{real}}^{(i)}\boldsymbol{m}\right)\right]$
(1)

Assuming negligible cyclic aliasing, the discrete Fourier transform of
$\psi_{\boldsymbol{m}}^{(i)}$ from Eq. (1) is guaranteed to return the
continuous Fourier transform $\Psi^{(i)}$ of
$\psi^{(i)}(\mathcal{R}^{-1}_{i}\boldsymbol{x})$, sampled on the conjugate
grid $\boldsymbol{B}^{(i)}_{\text{recip}}\boldsymbol{n}$, where
$\boldsymbol{n}\in\mathbb{Z}^{3}$.

$\psi^{(i)}_{\boldsymbol{m}}=\psi^{(i)}\left(\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{real}}^{(i)}\boldsymbol{m}\right)\underset{\mathcal{F}^{-1}}{\stackrel{{\scriptstyle\mathcal{F}}}{{\rightleftharpoons}}}\Psi^{(i)}\left(\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{recip}}^{(i)}\boldsymbol{n}\right)=\Psi^{(i)}_{\boldsymbol{n}}$
(2)

This assumption of cyclic aliasing is discussed briefly in Section 4 For
correctly computed $\boldsymbol{B}_{\text{real}}^{(i)}$ (typically by
calculating $\boldsymbol{B}_{\text{recip}}^{(i)}$ first, see Ref. [21]), the
measured coherent diffraction data in the far field will be correctly
represented by $\Psi_{\boldsymbol{n}}^{(i)}$. The relation between each pair
of $\boldsymbol{B}^{(i)}_{\text{real}}$ and
$\boldsymbol{B}^{(i)}_{\text{recip}}$ depends on the specific diffraction
geometry and the size of the signal-space array chosen to encompass the signal
[21]. This conjugate basis formulation serves to conveniently bring some of
the capabilities of existing tools such as _xrayutilities_ [22, 23] into the
scope of automatic differentiation.

In our algorithm we choose to solve for the primary quantities of interest,
namely $\mathcal{A}$ and $\boldsymbol{u}$, in the laboratory frame, but apply
our Fourier-based re-gridding algorithm to $\psi^{(i)}_{\boldsymbol{m}}$ in
each of the $i$ detector frames in turn (_i.e._ , each of the $M$ coordinate
frames
$[\hat{\boldsymbol{k}}_{1}~{}\hat{\boldsymbol{k}}_{2}~{}\hat{\boldsymbol{k}}_{3}]_{i}$
from Figure 1), after appropriate transformation to that frame. We note that
Eq. (1) makes no explicit reference to the coordinate frame in which these
quantities are represented and therefore affords us this flexibility. In
Section 5, we switch to a particular detector frame to describe Fourier-based
interpolation.

  3. 3.

The global scaling factor $\chi_{i}$ in single-reflection BCDI is trivially a
constant determined by the signal strength, but cannot be ignored in the MR-
BCDI forward model. In our algorithm we treat the $\chi_{i}$ as additional
parameters to be optimized, to mathematically reconcile the different signal
energies with the constraint that $\mathcal{A}\in[0,1]$ (since they are
related through the Fourier transform).

The operation defined by
$T_{i}\equiv\mathcal{R}_{i}^{-1}\boldsymbol{B}_{\text{real}}^{(i)}$ in Eq. (1)
denotes a re-gridding transformation of the object wave achieved by the
combination of (i) rotation into the $i$’th Bragg condition followed by (ii)
resampling to achieve the effect of $\boldsymbol{B}_{\text{real}}^{(i)}$,
starting from the _discretized_ $\mathcal{A}(\boldsymbol{x})$ and
$\boldsymbol{u}(\boldsymbol{x})$ fields. This is achieved by our Fourier
interpolation method described in Section 5. With this in mind, the formal MR-
BCDI inverse problem is set up in terms of a loss function
$\mathcal{L}_{\text{multi}}$ defined by:

$\displaystyle\mathcal{L}_{\text{multi}}[A,\boldsymbol{u}]$
$\displaystyle=\sum\limits_{i=1}^{M}\frac{1}{N}_{\text{vox}}\sum\limits_{\boldsymbol{n}}\left[\left\lVert\mathcal{F}\left\\{\chi_{i}\mathcal{A}(T_{i}\boldsymbol{m})e^{\iota
2\pi\boldsymbol{G}_{i}^{T}\boldsymbol{u}(T_{i}\boldsymbol{m})}\right\\}\right\rVert_{\boldsymbol{n}}-\left(I_{\boldsymbol{n}}^{(i)}\right)^{1/2}\right]^{2}$
(3)
$\displaystyle\mathcal{A}^{\star},\boldsymbol{u}^{\star},\boldsymbol{\chi}^{\star}$
$\displaystyle=\arg\min\limits_{\mathcal{A},\boldsymbol{u},\boldsymbol{\chi}}\mathcal{L}_{\text{multi}}\left[\mathcal{A},\boldsymbol{u},\boldsymbol{\chi}\right]$
(4)

where $\boldsymbol{\chi}$ represents the set of $M$ scaling factors.
$I^{(i)}_{\boldsymbol{n}}$ is the measured diffraction data corresponding to
the $i$’th Bragg reflection and $\boldsymbol{n}\in\mathbb{Z}^{3}$ ranges over
the pixels in each $I^{(i)}_{\boldsymbol{n}}$. Eq. 3 represents the mean error
in signal amplitude per BCDI data set, aggregated over the number of data
sets. Here, $N_{\text{vox}}$ is the number of voxels in a single BCDI data set
($128^{3}=2097152$ throughout this paper). In our stochastic gradient descent
scheme, the outermost summation in Eq. (3) was evaluated either over all the
$M$ BCDI data sets, or a subset thereof. It is an extension of the familiar
phase retrieval loss function typically minimized by FPIP [20], based on the
square root transformation of the measured data for robustness. Further, for
Eq. (4) to be over-determined in $\boldsymbol{u}(\boldsymbol{x})$, we require
$M\geq 3$ [12, 14, 13, 24]. To analytically constrain $\mathcal{A}\in[0,1]$ in
practice, we define an activation function
$\alpha:\mathbb{R}^{3}\rightarrow\mathbb{R}$ such that:

$\mathcal{A}(\boldsymbol{x})\equiv\mathcal{A}[\alpha(\boldsymbol{x})]=\frac{1}{2}\left[1+\tanh\left(\frac{\alpha(\boldsymbol{x})}{\alpha_{0}}\right)\right]$
(5)

With this definition, we see that $\mathcal{A}\simeq 0$ for
$\alpha(\boldsymbol{x})\ll 0$ and $\mathcal{A}\simeq 1$ for
$\alpha(\boldsymbol{x})\gg 1$. The transition between $0$ and $1$ in the
neighborhood of $\alpha(\boldsymbol{x})=0$ is controlled by the user-defined
positive constant $\alpha_{0}$.

Figure 2:  Activation function $\alpha(\boldsymbol{x})$ used to constrain
$\mathcal{A}(\boldsymbol{x})\in[0,1]$, for different activation parameters
$\alpha_{0}$.

Qualitatively, a larger value of the hyperparameter $\alpha_{0}$ results in a
smaller gradient descent step for $\mathcal{A}$ at $\boldsymbol{x}$ (as is
seen by differentiating Eq. (5) with respect to $\alpha$), and a more gradual
and controlled approach of $\mathcal{A}$ to its intended value. Smaller
$\alpha_{0}$ results in $\mathcal{A}$ converging through sporadic switching
between values close to $0$ and $1$, and as such is less controlled.
$\alpha_{0}$ is therefore interpreted as parameterizing how strongly the
gradient descent solution oscillates around the local minimum in the subspace
of $\alpha(\boldsymbol{x})$ as it approaches convergence. We found that
$\alpha_{0}=1$ best served our global convergence rate, and have used this for
all reconstructions in this paper. Given this, the actual optimization problem
to be solved in practice is:

$\displaystyle\alpha^{\star},\boldsymbol{u}^{\star},\boldsymbol{\chi}^{\star}$
$\displaystyle=\arg\min\limits_{\alpha,\boldsymbol{u},\boldsymbol{\chi}}\mathcal{L}_{\text{multi}}\left[\mathcal{A}[\alpha],\boldsymbol{u},\boldsymbol{\chi}\right]$
(6) $\displaystyle\mathcal{A}^{\star}$
$\displaystyle=\mathcal{A}[\alpha^{\star}]$ (7)

We also note that this optimization procedure as described does away with the
use of FPIP algorithms altogether. Significantly, it also avoids shrinkwrap
for support update, provided that the initial bounding box $\mathcal{V}$ is
chosen judiciously.

The size of $\mathcal{V}$ relative to the buffered array is an important
factor in overdetermining the optimization problem. This is seen as follows.
We assume $\mathcal{V}$ to have an edge size of $n$ voxels, within a buffered
array of size $N$ (with $n<N$). Therefore we have $M+4n^{3}$ unknowns
($\mathcal{A}$, $\boldsymbol{u}$ and $\boldsymbol{\chi}$ combined), where $M$
is the number of BCDI data sets. The number of measurements (photon counts) is
given by: $MN^{3}$. The MR-BCDI problem is therefore over-determined by a
factor of $MN^{3}/(M+4n^{3})$. In the case of our dislocation-free crystal,
$M=5$, $N=128$ and $n=46$ and therefore the overdetermination factor is $\sim
27$ which is significantly higher than the minimum required $1$ and therefore
highly desirable. Further, since our MR-BCDI reconstruction uses the fast
Fourier transform, the issue of signal oversampling must be addressed for a
truly well-conditioned inverse problem. This is achieved by ensuring that the
bounding box $\mathcal{V}$, _after transformation into each of the $M$
detector frames_, spans $\leq N/2$ voxels along each array axis. The
transformed $\mathcal{V}$ is no longer cubic, but in general sheared into a
parallelopiped.

Lastly, we note that the optimal solution in Eq. (4) does not preclude the
_universal twin_ solution, akin to the twin image degeneracy in single-
reflection BCDI [25]. In other words, the crystal defined by the
transformations
$\mathcal{A}(\boldsymbol{x})\rightarrow\mathcal{A}(-\boldsymbol{x})$ and
$\boldsymbol{u}(\boldsymbol{x})\rightarrow-\boldsymbol{u}(-\boldsymbol{x})$ is
also a solution to the MR-BCDI problem, in that it results in the twinned
version of the object wave for _all_ Bragg conditions, and therefore the same
inferred diffraction patterns. In this sense, MR-BCDI is still susceptible to
a reconstruction degeneracy. A related benefit of this coupled reconstruction
approach over traditional FPIP methods is that the global optimization will
converge to either the physically correct structure or the universal twin
structure automatically, without requiring manual intervention to align
individual phase retrieval instances.

##  3 Results

In this section we present three reconstruction results with the MR-BCDI
technique described above. Two of these are with numerically synthesized
digital nanocrystals of gold, with coherent diffraction signal simulated in
the far field. In these numerical studies, one crystal has a slowly varying
continuous distortion field (_i.e._ , strained with no line defects), and the
other has a discontinuous, winding displacement field consistent with two
orthogonal screw dislocations within the crystal. The third reconstruction
result is from experimental diffraction data of a silicon carbide (SiC)
nanocrystal consisting of six Bragg reflections measured using BCDI at
Beamline 34-ID of the Advanced Photon Source. In each case, we present the MR-
BCDI reconstruction of the local electron density
$\mathcal{A}(\boldsymbol{x})$ and the lattice deformation
$\boldsymbol{u}(\boldsymbol{x})$. The MR-BCDI reconstruction software was
written as a Python module for GPU hardware using the Adam optimizer routine
bundled in PyTorch [26] and is available for general use [27]. All
reconstructions were performed using a single Nvidia Tesla P100 GPU with 16 GB
of RAM.

We have adopted the matrix coordinate convention for all the cross-section
plots in this paper. As an example, a cross-section image labeled as $X-Y$
implies that the X-axis is directed from top to bottom along the image, and
the Y-axis from left to right. The third axis (in this example, the Z-axis)
emerges out of the plane of the figure in order to maintain right-handedness.
The $X$, $Y$ and $Z$ axes correspond to the laboratory frame of reference
(_i.e._ , the $\hat{\boldsymbol{s}}_{1}$, $\hat{\boldsymbol{s}}_{2}$ and
$\hat{\boldsymbol{s}}_{3}$ directions in Figure 1 respectively). With this
convention, we ensure that the reference frame of the reconstructed crystal is
consistent with the programmed order of the array axes in the multi-
dimensional FFT routines in Python and PyTorch.

###  3.1 Simulated crystal without dislocations

A synthetic crystal with arbitrary facets was generated on a Cartesian grid of
size $128\times 128\times 128$ voxels, with a voxel size of $s_{0}=12$ nm
along each axis in the laboratory frame
($[\hat{\boldsymbol{s}}_{1}~{}\hat{\boldsymbol{s}}_{2}~{}\hat{\boldsymbol{s}}_{3}]$
in Figure 1). A face-centered cubic (fcc) gold lattice was assumed with a
lattice constant of $a_{0}=4.078$ Å. A slowly-varying internal field
$\boldsymbol{u}(\boldsymbol{x})$ was generated by first generating uniform
random samples for the three components of $\boldsymbol{u}(\boldsymbol{x})$
($-0.1a_{0}\leq u_{i}(\boldsymbol{x})\leq 0.1a_{0}$) and then retaining the
slow variations by convolving each $u_{i}(\boldsymbol{x})$ with a low-pass
filter. The crystal spanned $(39,39,40)$ voxels along the laboratory frame
axes, and therefore a cubic bounding box $\mathcal{V}$ of size $46\times
46\times 46$ voxels was chosen for reconstruction purposes. This highly
deliberate choice of the bounding box size is for demonstration purposes only.
In practice it may be chosen based on the size of the Patterson function of
the diffraction signal [28]. An arbitrary initial lattice orientation was
assigned to this digital crystal in the laboratory frame. Assuming an x-ray
photon energy of $9$ keV, the crystal was actively rotated about
$\hat{\boldsymbol{s}}_{2}$ into several Bragg conditions. The resultant
complex-valued object wave was passively transformed into the corresponding
detector frame
$[\hat{\boldsymbol{k}}_{1}~{}\hat{\boldsymbol{k}}_{2}~{}\hat{\boldsymbol{k}}_{3}]$
using the method described in 5.3. This two-step composite rotation is denoted
by $\mathcal{R}_{i}$ in Eq. (1). The real-space sampling grid of the rotated
crystal was then transformed to the span of the basis
$\boldsymbol{B}_{\text{real}}^{(i)}$ using the methods described in Sections
5.1 (bulk) and 5.2 (shear). The rocking direction encoded in the third column
of $\boldsymbol{B}_{\text{recip}}^{(i)}$ (_i.e._ , $\boldsymbol{q}_{k}$) was
chosen as either a rotation about the $\theta$ or $\phi$ axis (see Figure 1)
based on whether the basis vectors
$[\boldsymbol{q}_{i}~{}\boldsymbol{q}_{j}~{}\boldsymbol{q}_{k}]$ were
maximally orthogonal to each other. Appendix A describes the customized
‘mutual orthogonality’ (MO) metric used in this choice. For the purpose of
this demonstration, five scans were selected based on the highest reciprocal-
space MO, corresponding to the $[1\bar{1}\bar{1}]$, $[\bar{1}\bar{1}\bar{1}]$,
$[\bar{2}\bar{2}0]$, $[20\bar{2}]$ and $[0\bar{2}\bar{2}]$ lattice planes in
their respective Bragg conditions. A straightforward 3D FFT of the resampled
complex-valued object wave in each case automatically gave the diffracted far-
field wave consistent with the correct sampling of a real-world BCDI scan. The
peak intensity for each far-field wave was set to $10^{5}$ photon counts.
Subsequent addition of Poisson noise to this rescaled diffraction pattern was
sufficient to emulate a high-quality signal collected with a synchrotron-based
BCDI instrument. Thus, the total number of unknowns to solve for was: $4\times
46^{3}+5=389349$. The prefactor $4$ comes from the electron density
$\mathcal{A}$ and the three components of the vector field
$\boldsymbol{u}(\boldsymbol{x})$ at each point $\boldsymbol{x}$. Further,
there are $5$ global scaling factors $\chi_{i}$ corresponding to each
simulated Bragg peak. During the MR-BCDI reconstruction, the crystal was
initialized to a constant-amplitude cube ($\mathcal{A}=1$) occupying the
entire $46\times 46\times 46$ bounding box with no interior lattice distortion
($\boldsymbol{u}=0$). The former was achieved by setting
$\alpha(\boldsymbol{x})$ from Eq. (5) to 2 and $\alpha_{0}=1$. Each of the
global scaling factors $\chi_{i}$ was initialized for the object wave to match
the total energy of the corresponding simulated signal. An Adam optimizer [18]
with an initial learning rate of 0.005 was employed. A total variation (TV)
regularizer was applied to $\alpha(\boldsymbol{x})$ with a coupling constant
of $10^{-5}$. During optimization, $\boldsymbol{u}(\boldsymbol{x})$ was
constrained to lie within the lattice plane separations $\pm d_{hkl}/2$
corresponding to $[hkl]=[100],[010]$ and $[001]$. This constraint may be
employed without loss of generality and is discussed further in Section 4.

The optimization was carried out sequentially over increasingly larger
‘minibatches’ of the BCDI data consisting of randomly selected scans, with the
final stretch of optimization being performed on all 5 scans. Table 1 in the
Appendix summarizes the sizes of the randomized minibatches and optimization
iterations therewith. After the execution of this plan, the voxels in the
crystal interior alone (for which $\mathcal{A}>0.2$) were optimized for
another 1000 iterations over all 5 scans, with the other voxels held constant
at their reconstructed values. This was done to refine the solution even
further within the smaller space consisting of the most relevant optimization
dimensions. Additionally, after each optimization epoch, a median filter with
a kernel size of $3\times 3\times 3$ voxels was applied to the reconstructed
lattice distortion field $\boldsymbol{u}(\boldsymbol{x})$ in order to remove
spurious isolated discontinuities induced by phase wraps. The median filter
was found to be better suited to the task of smoothing
$\boldsymbol{u}(\boldsymbol{x})$ than the TV regularizer used for
$\alpha(\boldsymbol{x})$. The entire optimization took $\sim 2$ hours on an
Nvidia Tesla P100 GPU.

Figure 3:  Ground truth (top) and reconstructed (bottom) electron density
cross-sections ($\mathcal{A}(\boldsymbol{x})$) of the synthetic strained
crystal with no internal dislocations. The array size represents the original
bounding box chosen for reconstruction, _i.e._ , 46 pixels along each lab
frame axis. The spatial resolution was estimated by fitting an error function
to the reconstructed electron density along the dashed line (see Figure 5(c)).
Figure 4:  Orthogonal cross sections of (a) ground truth and (b) reconstructed
lattice distortion components $u_{1}(\boldsymbol{x})$, $u_{2}(\boldsymbol{x})$
and $u_{3}(\boldsymbol{x})$, for the dislocation-free synthetic crystal. The
length scale is the same as in Figure 3.  Figure 5:  (a) Trend in multi-
reflection loss function given by Eq. (3) for the dislocation-free synthetic
crystal. The red lines denote the beginning of a new optimization epoch. (b)
Histogram of the simulated and reconstructed components of the vector
$\boldsymbol{u}(\boldsymbol{x})$, along with point-to-point residuals of each
vector component. (c) Error function fit to the reconstructed profile along
the dashed line in Figure 3 (X-Y slice). The spatial resolution was estimated
from this fit to be $\simeq 1.86$ pixels, or $\simeq 23$ nm.  Figure 6:
Cross-sections of the noisy simulated (top) and reconstructed (bottom)
coherent diffraction patterns from the five Bragg peaks chosen for the
dislocation-free synthetic crystal. The reconstructed diffraction patterns are
shown here without the scaling effect of the $\chi_{i}$, but rather scaled to
match the corresponding simulated intensity peaks. To fully match the color
scales, each reconstructed diffraction was clipped below to the smallest
nonzero photon count in the corresponding signal (_i.e._ , 1 photon).

Figure 3 shows the ground truth and the reconstructed electron density
$\mathcal{A}$, on the original grid of step size $s_{0}=12$ nm. The spatial
resolution of the reconstruction was estimated by fitting an error function to
the density profile along the dashed line shown in Figure 3. This error
function is shown in Figure 5(c) and has a width of $\sigma=0.788$ pixels,
which is the standard deviation of the underlying Gaussian. From this, the
spatial resolution was estimated to be $\sqrt{8\ln 2}\sigma\simeq 1.86$
pixels. Figure 4 shows three orthogonal cross-sections of the three components
of $\boldsymbol{u}(\boldsymbol{x})$ (ground truth and reconstructed). Figure
5(a) shows the progression of the loss function
$\mathcal{L}_{\text{multi}}[\mathcal{A},\boldsymbol{u}]$ over the entire
optimization process while Figure 5(b) shows the histograms of the voxel-by-
voxel ground truth and reconstructed vector components of
$\boldsymbol{u}(\boldsymbol{x})$ (all three vector components in the same
histogram), as well as the corresponding point-to-point error (residuals). The
dashed lines in Figure 5(a) indicate the beginning of each new optimization
epoch. Within each epoch, the rapid oscillations are attributed to the abrupt
change in the loss function landscape due to each new randomized minibatch.
For example, in the first epoch, the loss function is successively optimized
over 400 sets of of 2 randomly selected BCDI scans for 6 iterations each (from
Table 1). Whenever a new randomized set of scans is used, the energy function
landscape abruptly changes from its previous state, and the gradient descent
is now, in general, forced along a completely different slope.

We note that that the majority of the residuals in Figure 5(b) are
concentrated around $0$, and the second peak in residuals originates from the
voxels within the margin of mismatch between the estimated object edge with
respect to the ground truth. Prior to calculating these residuals, the
reconstruction was spatially aligned with the ground truth by aligning their
respective centroids in $\mathcal{A}$, to minimize the effect of any
translational offsets. The residuals for voxels at which $\mathcal{A}<0.5$
were ignored in the histogram.

Figure 6 shows the comparison of the simulated and reconstructed diffraction
patterns. The top row depicts the simulated detector measurement at the Bragg
peak of each scan. The reconstructed diffraction patterns shown here do not
contain the scaling effect of their respective $\chi_{i}$. We note the good
agreement between the simulated and reconstructed diffraction, which is
quantified in Figure 5(a). This shows $\mathcal{L}_{\text{multi}}\sim 1$,
which is the estimate of the aggregate error in signal amplitude per voxel
over all the BCDI datasets.

Figure 7:  (a) Phase retrieval reconstructions of the object wave amplitudes
($\left\lVert\psi^{(i)}\right\rVert$), corresponding to the 5 individual BCDI
scans from the dislocation-free crystal. The voxel sizes in nanometers are the
magnitudes of the columns of the corresponding
$\boldsymbol{B}_{\text{real}}^{(i)}$ matrix. We note that these basis vectors
are not orthogonal in general, and therefore the slices shown in each column
do not correspond to orthogonal slices. (b) Magnified view of one of the
reconstructions, with a line profile extracted along the dashed line. (c)
Error function fit to the line profile in (b), with an estimated spatial
resolution of $1.3\sqrt{8\ln 2}=3.06$ pixels, or $1.3\sqrt{8\ln 2}\times
13=39.8$ nm. The multiplier of $13$ comes for the $\Delta k$ pixel pitch in
the corresponding image in (a).

We now compare these reconstructions to those of conventional phase retrieval.
Figure 7 shows the results of individual phase retrieval reconstructions of
the 5 BCDI scans simulated for this crystal (object wave amplitude). If the
pixels of the final reconstruction $\psi_{ijk}$ in each case are indexed by
integers $i$, $j$ and $k$, then each row in Figure 7(a) corresponds to the
central slice along the $i-j$, $k-i$ and $j-k$ planes. These slices are in
general not orthogonal, since the real-space sampling directions encoded in
each $\boldsymbol{B}_{\text{real}}$ are not orthogonal. The phase retrieval
was performed with the open-source Python module Phaser [29]. The phase
retrieval recipe chosen was a combination of ER, HIO and SF (solvent flipping)
with intermittent shrinkwrap-based update of the object support. The exact
recipe is given in Section C of the Appendix. Each reconstruction took 0.5
seconds with an Nvidia Tesla P100 GPU. An error function fit to the line
profile in 7(b) is shown in 7(c), with an estimated spatial resolution of
$39.8$ nm. We see that the reconstructed electron density from Figure 3, and
therefore any subsequent object wave amplitude, is more uniform and exhibits
better spatial resolution than those obtained obtained from conventional phase
retrieval.

###  3.2 Simulated crystal with screw dislocations

A synthetic gold crystal ($a_{0}=4.078$ Å) with arbitrary facets was generated
on the same Cartesian grid with the same voxel size as the dislocation-free
case ($s_{0}=12$ nm), with an arbitrary initial crystal orientation. The
$\boldsymbol{u}(\boldsymbol{x})$ field at each voxel was consistent with the
presence of two spatially separated screw dislocations within the crystal
volume, whose cores were along the orthogonal $[111]$ and $[2\bar{2}0]$
crystallographic directions. In this paper we favor the well-known continuum
model of lattice distortion for screw dislocations [30, 31] for our simulation
over the more sophisticated, discontinuity-free atomistic simulations that use
known embedded atomic potentials [32]. The former introduces an unphysical
discontinuity in the lattice distortion field, which nevertheless generates
the correct outgoing object wave at the desired discretization of $\sim 10$
nm. The Burgers vector magnitude $\left\lVert\boldsymbol{b}\right\rVert$ for
each screw dislocation was set to the corresponding lattice-plane spacing.
Namely, $\left\lVert\boldsymbol{b}\right\rVert_{111}=a_{0}/\sqrt{3}$ and
$\left\lVert\boldsymbol{b}\right\rVert_{2\bar{2}0}=a_{0}/\sqrt{8}$, where
$a_{0}$ is the lattice constant of gold. The net displacement field was
modeled as the vector sum of the individual displacement fields due to each
screw dislocation. The simulation procedure of the BCDI data sets, including
the selection criterion for Bragg reflections, rocking directions and
introduction of Poisson noise, was identical to that for the dislocation-free
crystal. Four peaks with the best reciprocal-space mutual orthogonality (MO)
were chosen for reconstruction, corresponding to the $[200]$, $[002]$, $[202]$
and $[2\bar{2}0]$ reciprocal lattice vectors. The simulation box $\mathcal{V}$
was chosen to be $40\times 40\times 40$ voxels in size. The number of
optimization variables in this case was therefore $4\times 40^{3}+4=256004$.
The Adam optimizer was initialized with a learning rate of 0.02, with the
chosen sequential optimization plan shown in the Appendix Table 2.

The final optimization stretch for the interior voxels over all 4 scans was
carried out for 5000 iterations. As with the dislocation-free crystal, a
median filter of size $7\times 7\times 7$ voxels was applied to the components
of $\boldsymbol{u}(\boldsymbol{x})$ after each optimization epoch. This
optimization took $\sim 3.5$ hours on an Nvidia Tesla P100 GPU.

Figure 8:  Ground truth (top) and reconstructed (bottom) electron density
cross-sections of the synthetic strained crystal with two orthogonal screw
dislocations. The array size represents the original bounding box chosen for
reconstruction, _i.e._ , 40 pixels along each lab frame axis. The arrows
indicate the points of intersection of the dislocation cores with the figure
plane.  Figure 9:  Orthogonal cross sections of (a) ground truth and (b)
reconstructed lattice distortion components $u_{1}(\boldsymbol{x})$,
$u_{2}(\boldsymbol{x})$ and $u_{3}(\boldsymbol{x})$. The length scale is the
same as Figure 3. The regions of discontinuity in the lattice distortion match
with the drop in electron density in Figure 8. The circled regions show
imperfect reconstruction at locations where the lattice discontinuity meets
the crystal surface.  Figure 10:  (a) Trend in multi-reflection loss function
given by Eq. (3) as a function of Adam optimizer iteration (initial learning
rate = 0.02) The red lines denote the beginning of a new optimization epoch.
(b) Histogram of the simulated and reconstructed components of the vector
$\boldsymbol{u}(\boldsymbol{x})$, along with point-to-point residuals of each
vector component.  Figure 11:  Cross-sections of the noisy simulated (top) and
reconstructed (bottom) coherent diffraction patterns from the five Bragg peaks
chosen for the gold crystal with two orthogonal screw dislocations.

We again note the good agreement between the simulated and reconstructed
crystals (Figures 8 and 9). In Figure 8, we note the drop in reconstructed
electron density along the lines where the half-plane of discontinuous lattice
distortion intersects the plane of the figure. The terminal points of these
“trenches” of low density are the locations of the dislocation cores (marked
by arrows). These trenches are discussed further in Section 4. In Figure 9 we
note the regions circled in red where the discontinuity in the reconstructed
$\boldsymbol{u}(\boldsymbol{x})$ intersects the crystal surface (_i.e._ , an
abrupt drop in $\mathcal{A}$) and results in an imperfect reconstruction of
$\boldsymbol{u}(\boldsymbol{x})$. Both these effects may be attributed to the
inherently band-limited nature of the reconstruction, due to Poisson noise
-induced information loss.

Figure 12:  (a) Individual amplitude reconstructions
$\left\lVert\psi^{(i)}\right\rVert$ from phase retrieval applied to the BCDI
scan the crystal with screw dislocations (complex amplitude), using the recipe
in Appendix Section C. The cross sections shown correspond to the same slices
in Figure 7.

Figure 10 shows the reconstruction metrics for this simulated structure
similar to those in Figure 5 for the dislocation-free simulated crystal.
Figure 11 compares the simulated and reconstructed diffraction patterns from
the four chosen Bragg peaks, on the same color scale. Figure 12 shows the
individual phase retrieval reconstructions of the BCDI scans from the
dislocated crystal lattice. The recipe used was the same for the dislocation-
free crystal in Figure 7, detailed in the Appendix section C. We note the
abject failure of the conventional phase retrieval approach in the presence of
dislocations, compared to the reconstruction in Figures 8 and 9.
Reconstructing dislocated crystals with FPIP algorithms has thus far most
successfully been achieved with parallelized reconstruction threads in a
guided (_i.e._ , genetic algorithm) approach [33, 14].

###  3.3 Silicon carbide nanocrystal

Having demonstrated MR-BCDI on synthetic datasets, we now move on to
experimental data acquired from real crystals. We chose silicon carbide (SiC)
as a case study, which is a promising platform for quantum sensing and quantum
information technologies [34, 35, 36, 37, 38, 39]. The crystal structure is a
4H polymorph with lattice parameters of $a=3.073$ Å and $c=10.053$ Åwith
quasi-hexagonal and quasi-cubic lattice planes. Arrays of sub-$\mu$m sized,
‘D’-shaped SiC crystals were fabricated out of an _i-p-n_ SiC wafer procured
from Norstel ($[11\bar{2}0]$ major flat, $4^{\circ}$ miscut from the
crystallogrphic surface normal$[0001]$, see Figure 16(d)) and transferred onto
a silicon wafer. The SiC crystals had known positions, nominally aligned
orientations, and separations substantially larger than the x-ray beam
footprint. This enabled us to collect BCDI data from 6 independent Bragg
reflections from an isolated SiC crystal at the 34-ID-C BCDI instrument at the
Advanced Photon Source. These reflections corresponded to the $[10\bar{1}1]$,
$[10\bar{1}0]$, $[01\bar{1}2]$, $[01\bar{1}0]$, $[01\bar{1}1]$ and
$[10\bar{1}3]$ lattice planes. The fabrication process is described in detail
in Section B of the Appendix. The following pre-processing steps were
performed on this experimental data set prior to the MR-BCDI optimization:

  1. 1.

The raw BCDI data from each Bragg reflection was collected as a sequence of
$256\times 256$ images with $80$ steps along the rocking curve. The data sets
were clipped to a size of $128\times 128$ pixels in the imaging plane and
zero-padded to $128$ steps along the rocking curve to give a cubic simulation
volume $128\times 128\times 128$ voxels in size. In addition to reducing
memory requirements, this uniform array size ensured that the FFT-based
implementation of the shear operator described in Section 5.2 does not suffer
from signal normalization issues inherent to applying the FFT along individual
array dimensions, while leaving the other dimensions untouched. Future
iterations of this reconstruction tool will be generalized to non-uninform
array dimensions.

  2. 2.

Unlike with the simulations, the precise lab frame crystal orientation was
known only approximately prior to reconstruction and had to be refined. This
was done by a simple least squares optimization of the crystal orientation
(expressed as a 3-parameter rotation vector) over the $\delta$, $\gamma$ and
$\theta$ motor positions (see Figure 1) of all 6 BCDI scans with the Nelder-
Mead optimizer [40] available in the SciPy software package.

The maximum pixel intensities in each of the 6 measured Bragg peaks were
respectively observed to be: 40980, 16072, 32661, 14577, 34123 and 7805
photons.

Figure 13:  Loss function trend for the SiC crystal reconstruction, for the
optimization plan shown in Table 3.  Figure 14:  Final estimated SiC
nanocrystal electron density $\mathcal{A}$, from a highly permissive mask
threshold of $0.1$. The red line is the contour at a threshold of
$\mathcal{A}=0.65$, representing the approximate crystal surface.  Figure 15:
Reconstructed components of the lattice distortion within the SiC crystal
(centered at zero volume-averaged distortion in each dimension), within the
estimated crystal surface at $\mathcal{A}=0.65$.  Figure 16:  (a), (b), (c)
Isosurface plots of the reconstructed $u_{1}(\boldsymbol{x})$,
$u_{2}(\boldsymbol{x})$ and $u_{3}(\boldsymbol{x})$ respectively in the SiC
crystal (color scale in nanometers), with an isovalue of $0.65$ in electron
density. (d) Schematic of the nanoparticle, with the ‘D’-shaped frustum cross
section, also visible in the cross section in Figure 14(c). Also indicated are
the $[0001]$ and $[11\bar{2}0]$ crystallographic planes relative to the
facets. (e) SEM image of the nanoparticle structure on the substrate, with the
diffracting SiC portion highlighted within the red circle.  Figure 17:
Measured (top) and reconstructed (bottom) diffraction patterns of SiC.

Figure 13 depicts the progression of the loss function for the optimization
scheme described in Table 3. Figure 14 shows the electron density after
optimization of the interior voxels, ‘interior’ being defined in this case as
above the highly permissive threshold of $\mathcal{A}>0.1$. Also shown is the
contour at $\mathcal{A}=0.65$, depicting the approximate surface of the
crystal. Figure 15 shows the reconstructed components of
$\boldsymbol{u}(\boldsymbol{x})$ along with the same contour. Figure 16 shows
the corresponding 3D isosurface plots of the components of $\boldsymbol{u}$,
with the color scale denoting the spatially varying lattice distortion in
nanometers. We see from the absence of discontinuities that the particle
appears to lack dislocations or sharply varying strain fields. Figure 17 shows
the measured and inferred diffraction patterns, with the reconstructed
patterns transformed to the ranges of the measured signals.

We see that the known ‘D’-shaped cross-section of the nanoparticle in Figure
16(d) is discernable in the electron density and lattice distortion
reconstructions shown in Figures 14 and 15, despite the loss of spatial
resolution as a result of clipping the BCDI data sets to a $128\times
128\times 128$ -pixel array. The relatively uniform frustum of the nanocrystal
in the 16(a), (b) and (c) isosurface plots is corroborated by the SEM image in
16(e). We note that the SEM image shows the highly porous _p_ -type tail from
the fabrication process. This volume did not diffract appreciably compared to
the _i_ -type section, so only the latter was reconstructed.

##  4 Discussion

In this paper we present a completely differentiable forward model for the MR-
BCDI reconstruction problem for vector-valued lattice distortion fields,
tailored to non-uniformly strained and defective crystalline nanoparticles. We
have shown reconstruction results from simulated and experimental coherent
diffraction data. By using a gradient descent-based global optimization over
randomized minibatches, we have obviated the need for trial-and-error FPIP
recipes and, more significantly, support update via the shrinkwrap algorithm.
We have demonstrated consistently better reconstructions for both defect-free
and defected crystals with dislocations that are isolated on the scale of BCDI
resolution. The reconciliation of different signal gridding geometries with
the Fourier transform is crucial for reliable imaging of such topologically
intricate structural heterogeneities within the crystalline nanoparticles.
This is particularly true for heterogeneities that manifest nominally as phase
wraps and discontinuities of the coherent diffracted wave, such as edge and
screw dislocations. We have demonstrated that BCDI data sets from independent
Bragg reflections, when of sufficiently high quality and explicitly coupled
together, can ‘conspire’ to give a more uniform electron density profile and
potentially better spatial resolution (as evidenced by the estimated $23$ nm
and $39.8$ nm spatial resolutions for the respective MR-BCDI and phase
retrieval reconstructions of the dislocation-free crystal). This is consistent
with previous experimental findings that data redundancy from coupled
diffraction data sets improves spatial resolution in BCDI [15].

On the topic of shrinkwrap [11], we note that this currently prevalent method
is essentially a ‘blind’ process contingent on reasonable prior convergence to
a compact shape during phase retrieval. Failing this, it becomes irrelevant
and potentially misleading, as our examples with the screw dislocations
indicate (Figure 12). Our global optimization approach avoids this
inconvenience entirely by simply making the support one of the products of the
fitting process. The only choices involved are the size of the initial
bounding box $\mathcal{V}$, and the hyperparameter $\alpha_{0}$ from Eq. (5).
In contrast, shrinkwrap involves two hyperparameters, namely the width of the
blurring Gaussian kernel (usually denoted $\sigma$) and an arbitrary threshold
below which the amplitude is suppressed to zero. In addition, there is the
more subtle choice of how often to employ shrinkwrap in a recipe, a largely
trial-and-error process when reconstructing heavily dislocated crystals.

This formulation of MR-BCDI offers another potential advantage in that it
permits the analytical interrogation of the high-dimensional solution space in
the neighborhood of the true solution ($4\times 46^{3}+5=389349$ dimensions in
the case of our dislocation-free synthetic crystal). Without the explicit
evaluation of a loss function, FPIP methods offer little to no insight in this
regard, preventing rigorous BCDI error analysis. In contrast, our analytical
prescription is a potential first step in developing rigorous error
propagation and uncertainty quantification capabilities, which is currently
lacking in current FPIP-based BCDI and will be in high demand as expected
improvement in data quality from fourth-generation synchrotrons pushes the
limits of spatial resolution.

In Section 3.1 we noted that we restricted $\boldsymbol{u}(\boldsymbol{x})$ to
lie within the lattice plane spacings defined by the $[100]$, $[010]$ and
$[001]$ directions of the crystal. These directions are not orthogonal in
general, as is seen from the equivalent $[10\bar{1}0]$, $[01\bar{1}0]$ and
$[0001]$ for our SiC crystal, in which the last two directions are separated
by $120^{\circ}$. This constraint is not a loss of generality when we consider
the conditions under which phase wrapping takes place for a particular Bragg
reflection. The phase at a point $\boldsymbol{x}$ corresponding to a
reflection $\boldsymbol{G}_{hkl}$ is given by: $\phi(\boldsymbol{x})\equiv
2\pi\boldsymbol{G}_{hkl}\cdot\boldsymbol{u}(\boldsymbol{x})$. Therefore a
phase wrap is associated with an excess displacement $\boldsymbol{d}_{hkl}$
given by:
$\phi(\boldsymbol{x})+2m\pi=2\pi\boldsymbol{G}_{hkl}\cdot[\boldsymbol{u}(\boldsymbol{x})+m\boldsymbol{d}_{hkl}]$,
where $m\in\mathbb{Z}$. Here, $\boldsymbol{d}_{hkl}$ is the vector in the
direction of $\boldsymbol{G}_{hkl}$ and whose magnitude is equal to the $hkl$
lattice plane spacing, with
$\left\lVert\boldsymbol{G}_{hkl}\cdot\boldsymbol{d}_{hkl}\right\rVert=1$ in
crystallographers’ units ($2\pi$ in physicists’ units). We see that for a
given $\boldsymbol{G}_{hkl}$, $\boldsymbol{u}(\boldsymbol{x})$ is determined
up to an integer multiple of $\boldsymbol{d}_{hkl}$ in the same sense that the
phase is only determined up to an integer multiple of $2\pi$. We therefore
deduce that any $\boldsymbol{u}(\boldsymbol{x})$ may be reduced to a
‘fundamental zone’ (terminology borrowed from crystallography) between the
parallel lattice planes located at
$\boldsymbol{u}(\boldsymbol{x})=\boldsymbol{0}$ and
$\boldsymbol{u}(\boldsymbol{x})=\boldsymbol{d}_{hkl}$ (or equivalently,
$\pm\boldsymbol{d}_{hkl}/2$), by addition or subtraction of an appropriate
number of $\boldsymbol{d}_{hkl}$ [24]. The reduced
$\boldsymbol{u}(\boldsymbol{x})$ is guaranteed to generate the same phase in
the outgoing object wave. Thus, as a purely computational convenience, we
constrain $\boldsymbol{u}(\boldsymbol{x})$ at each point to lie within the
intersection of the three largest ‘fundamental zones’ of the crystal lattice,
which correspond to the three largest (close packed) lattice plane separations
(_i.e._ , $\left\lVert\boldsymbol{d}_{100}\right\rVert$,
$\left\lVert\boldsymbol{d}_{010}\right\rVert$ and
$\left\lVert\boldsymbol{d}_{001}\right\rVert$). This volume of intersection is
rectangular for a cubic or orthorhombic crystal, but in general is a sheared
parallelopiped. This constraint serves the following purposes:

  1. 1.

It plays a role in preventing spurious discontinuities in the reconstructed
$\boldsymbol{u}(\boldsymbol{x})$ resulting from excessively large gradient
descent steps, which could cause the $\boldsymbol{u}(\boldsymbol{x})$ at
adjacent voxels to converge to different but phase-equivalent values.

  2. 2.

It retains faithful representation of the _genuine_ discontinuities in
$\boldsymbol{u}(\boldsymbol{x})$ from lattice defects that are equal in size
to a closed-packed lattice plane spacing, such as edge and screw dislocations.

We also note that the conjugate basis formalism of
$\boldsymbol{B}_{\text{real}}$ and $\boldsymbol{B}_{\text{recip}}$ is _not_
the most general parameterization of the signal space. In fact, it immediately
breaks down when imaging sufficiently small crystals, for which sufficient
fringe information is obtained over a much larger volume of Fourier space,
whose size is comparable to $\left\lVert\boldsymbol{G}_{hkl}\right\rVert$.
Under these circumstances, the curvature of the Ewald sphere can no longer be
ignored when rocking the sample. Accounting for this aspect would require the
design of better forward models, which is currently an active field of study.
Further, the use of the fast Fourier transform (FFT) leaves us open to the
problem of cyclic aliasing of the inferred scattering. This imposes a limit on
the spatial resolution of MR-BCDI that is not consistent with the amount of
Fourier-space information collected in a scan. Perversely, this effect is more
pronounced with a higher signal-to-noise ratio measurement, which is expected
at a fourth-generation synchrotron. Research to address this issue is ongoing
in the BCDI community [41].

All this taken into account, this work represents a significant step in the
imaging of defected materials on the $\sim 10$ nm scale, while eliminating the
need for _ad hoc_ phase retrieval recipe development and shrinkwrap-based
support updates. In addition, it is another step in bringing coherent
diffraction methods into the scope of global optimization techniques [42, 43],
thereby enabling the use of highly optimized software packages capable of
handling large data volumes and running on seemingly ever-improving high-
performance computing hardware.

##  5 Methods: Fourier transform-based resampling

We start with the discrete object $\psi(s_{0}\boldsymbol{m})$, which denotes
the object $\psi(\boldsymbol{x})$ sampled along the orthogonal directions
$[\hat{\boldsymbol{k}}_{1}~{}\hat{\boldsymbol{k}}_{2}~{}\hat{\boldsymbol{k}}_{3}]$
of the detector frame, in steps of size $s_{0}$. We assume the digitized
scatterer to be inside a cubic array of size $N_{1}\times N_{2}\times N_{3}$
pixels. In the MR-BCDI problem, this is one of $i$ mounting configurations
corresponding to a single Bragg condition. We now describe the basic
operations of (i) bulk deformation, (ii) shear deformation and (iii) rigid-
body rotation of the sampling grid.

The basis matrix $\boldsymbol{B}_{\text{real}}$ for this demonstrative
configuration is based on an actual BCDI measurement described in Ref. [21]:

$\boldsymbol{B}_{\text{real}}=\left[\begin{matrix}12.72&0&0\\\ 0&12.72&0\\\
3.21&2.70&11.87\end{matrix}\right]~{}\text{nanometers}$ (8)

The resampling method
$\psi(s_{0}\boldsymbol{m})\rightarrow\psi(\boldsymbol{B}_{\text{real}}\boldsymbol{m})$
consists of two interpolation operations acting on
$\psi(s_{0}\boldsymbol{m})$, performed in order: (i) ‘bulk’ resampling along
the orthogonal directions in _new_ steps of $B_{11}$, $B_{22}$ and $B_{33}$,
where $B_{ii}$ denote the diagonal elements of $\boldsymbol{B}_{\text{real}}$,
followed by (ii) a ‘shear’ resampling into each of the three independent
orthogonal directions, corresponding to the off-diagonal elements of
$\boldsymbol{B}_{\text{real}}$. We introduce for simplicity $\mathcal{F}_{i}$
as the 1-D Fourier transform along the axis $i$ (where $i=1,2,3$).

Further, we assign $s_{0}=(\det\boldsymbol{B}_{\text{real}})^{1/3}=12.43$ nm.
This choice allows us to demonstrate both bulk stretching and compression,
since is it greater than one of the $B_{ii}$ but less than the other two. For
an MR-BCDI experiment, $s_{0}$ is the size of the reconstruction grid within
the bounding box $\mathcal{V}$.

###  5.1 Bulk resampling

Real-space resampling along the orthogonal array axes of
$\psi(s_{0}\boldsymbol{m})$ is achieved by truncating or padding the numerical
array in Fourier space. To this end, we define operations $\mathcal{P}$ and
$\mathcal{T}$ that symmetrically pad or truncate the input array along the
specified axis:

  * •

$\mathcal{P}(\hat{\boldsymbol{k}}_{i},n)\psi_{\boldsymbol{m}}$ denotes the
operator $\mathcal{P}$ acting on an array $\psi_{\boldsymbol{m}}$ that returns
an array zero-padded on either side by $n$ pixels. This results in $2n$ more
pixels along the $\hat{\boldsymbol{k}}_{i}$ axis.

  * •

$\mathcal{T}(\hat{\boldsymbol{k}}_{i},n)\psi_{\boldsymbol{m}}$ denotes the
operator $\mathcal{T}$ acting on an array $\psi_{\boldsymbol{m}}$ that returns
an array truncated/clipped on either side by $n$ pixels. This results in $2n$
fewer pixels along the $\hat{\boldsymbol{k}}_{i}$ axis.

We seek to resample the real-space object array $\psi_{\boldsymbol{m}}$ along
axis $\hat{\boldsymbol{k}}_{i}$ in steps of $B_{ii}$ (different from the
original sampling steps size $s_{0}$). The number $n$ of pixels by which to
zero-pad/truncate is easily computed from the relation between the real-space
step size and the Fourier-space numerical aperture: $s_{0}\propto 1/N$. If
$N_{i}^{\text{(new)}}$ is the number of pixels along
$\hat{\boldsymbol{k}}_{i}$ after zero-padding/truncation, then we have:

$\displaystyle N_{i}^{\text{(new)}}$
$\displaystyle=\left(\frac{s_{0}}{B_{ii}}\right)N_{i}$ (9)
$\displaystyle\Longrightarrow n$ $\displaystyle=\frac{1}{2}\left\lvert
N_{i}^{\text{(new)}}-N\right\rvert=\frac{1}{2}\left\lvert
N\left(1-\frac{s_{0}}{B_{ii}}\right)\right\rvert$ (10)

Here $\left\lvert\cdot\right\rvert$ denotes the absolute value and the factor
$1/2$ comes from the requirement to pad/truncate symmetrically. Armed with
this notation, the resampling operator
$\mathcal{B}(\hat{\boldsymbol{k}}_{i}\left|s_{0}\rightarrow B_{ii}\right.)$ is
given by the following sequence of operations:

$\mathcal{B}(\hat{\boldsymbol{k}}_{i}\left|s_{0}\rightarrow
B_{ii}\right.)=\left\\{\begin{matrix}\mathcal{T}\left(\hat{\boldsymbol{k}}_{i},n\right)\mathcal{F}_{i}^{-1}\mathcal{P}\left(\hat{\boldsymbol{k}}_{i},n\right)\mathcal{F}_{i}&\text{if}&B_{ii}<s_{0}\\\
\mathcal{P}\left(\hat{\boldsymbol{k}}_{i},n\right)\mathcal{F}_{i}^{-1}\mathcal{T}\left(\hat{\boldsymbol{k}}_{i},n\right)\mathcal{F}_{i}&\text{if}&B_{ii}>s_{0}\end{matrix}\right.$
(11)

The bulk-resampled array is denoted as:
$\psi^{\prime}_{\boldsymbol{m}}=\mathcal{B}(\hat{\boldsymbol{k}}_{i}\left|s_{0}\rightarrow
B_{ii}\right.)\psi_{\boldsymbol{m}}$. The value of $n$ in Eq. (11) is computed
from Eq. (10), in practice rounded to the nearest integer. The final, bulk-
resampled object along all three orthogonal axes is given by:

$\displaystyle\psi^{\prime}_{\boldsymbol{m}}$
$\displaystyle\equiv\psi([B_{11}~{}B_{22}~{}B_{33}]\boldsymbol{m})$
$\displaystyle=\underbrace{\mathcal{B}(\hat{\boldsymbol{k}}_{1}\left|s_{0}\rightarrow
B_{11}\right.)\mathcal{B}(\hat{\boldsymbol{k}}_{2}\left|s_{0}\rightarrow
B_{22}\right.)\mathcal{B}(\hat{\boldsymbol{k}}_{3}\left|s_{0}\rightarrow
B_{33}\right.)}_{\text{total bulk resampling operator}}\psi_{\boldsymbol{m}}$
(12)

Bulk resampling of a synthetic object originally sampled in steps of $s_{0}$
is shown in Fig. 18. The first two rows show the amplitudes of the original
and bulk-resampled object. The third and fourth rows show the corresponding
complex phases. The winding phase profile is indicative of a screw dislocation
passing through the center of the crystalline bulk (third row). We note the
fidelity of the resampled phase to the original object, specifically the
discontinuity at the dislocation core. This could not have been achieved with
a simple interpolation within a voxel, as seen in present-day methods [12, 13,
15].

Figure 18:  Fourier transform -based bulk resampling demonstrated on a
synthetic crystalline volume with a screw dislocation, marked by a winding
phase at the dislocation core. Row 1: orthogonal slices through center of the
amplitude profile of the original crystal. The arrows denote the bulk
compression/expansion along the principal axes with respect to the original
object, as a result of resampling. Row 2: central cross-sections of the
resampled amplitude. Row 3: corresponding slices of the original phase, with
the screw dislocation indicated by the phase discontinuity. Row 4: resampled
phase profile. The plots show phase in radians.

Further, there is little difference in the sizes of the original and resampled
objects because of the proximity of each $B_{ii}$ to the chosen value of
$s_{0}$. For demonstrative purposes, the fidelity of the bulk-resampling
procedure to high-frequency phase features even with greatly exaggerated
under/over-sampling steps is shown in the Appendix Section D.

We note the following:

  1. 1.

The concept of upsampling with the padding operator $\mathcal{P}$ has been
explored in the context of image registration [44].

  2. 2.

The operator $\mathcal{B}$ in Eq. (11) returns a new array of size
$N_{1}\times N_{2}\times N_{3}$ pixels which contains the original object
resampled along axis $i$ by a factor $B_{ii}/s_{0}$.

  3. 3.

The upsampling enabled by Eq. (11) when $B_{ii}<s_{0}$ does not actually
result in better sampling of any sub-pixel features in the physical object.
The smallest feature dictated by the available Fourier content is merely
sampled more finely. This is owing to the fact that no new Fourier information
was added in the process. However, on truncation ($B_{ii}>s_{0}$), high-
frequency information in Fourier space is truly discarded. As expected, the
downsampling operation results in a loss of spatial resolution.

  4. 4.

While truncation and zero-padding are not differentiable operations by
themselves, _every pixel value_ in the output array
$\psi^{\prime}_{\boldsymbol{m}}$ is still related analytically to the pixels
in $\psi_{\boldsymbol{m}}$ through the Fourier transform. In the MR-BCDI
problem, this analytical dependency extends back to the original unknowns,
namely $A_{\boldsymbol{m}}$ and $\boldsymbol{u}_{\boldsymbol{m}}$ through
multiplication and exponentiation, as well as each of the inferred diffraction
patterns $\left\lVert\Psi^{(i)}_{\boldsymbol{m}}\right\rVert^{2}$ and the loss
function in Eq. (3) itself. Multiplication, exponentiation and the Fourier
transform are all differentiable functions and have highly optimized
implementations in modern auto-differentiation software. In the language of
auto-differentiation, the computational graph is unbroken from
$\mathcal{A}_{\boldsymbol{m}}$ and $\boldsymbol{u}_{\boldsymbol{m}}$ to
$\mathcal{L}_{\text{multi}}$, even in the presence of a “destructive”
operation like truncation. Therefore this formulation lends itself easily to
auto-differentiation, which has already started to underlie computational
imaging applications with x-rays [42, 45, 46, 47].

  5. 5.

The $\mathcal{B}$ operators in Eq. (12) commute with each other since the
resampling along each of the three axes are independent operations.

###  5.2 Shear resampling

The bulk resampling method described in Section. 5.1 was achieved with non-
analytical operations such as zero-padding and truncation, which nevertheless
connect the unknowns $A$ and $\boldsymbol{u}$ analytically to the inferred
signal. In contrast, the shear transformations corresponding to the off-
diagonal elements of $\boldsymbol{B}_{\text{real}}$ are truly differentiable
in nature. As we show in this subsection, they can be achieved through a
combination of 1-D FTs and phase ramps. The 2D shear of a general 3D function
$\xi(\boldsymbol{x})\equiv\xi(x,y,z)$ in the $yz$-plane is given by the
following sequence of operations:

$\displaystyle\xi(x+\alpha y+\beta z,y,z)$ $\displaystyle=$
$\displaystyle\underbrace{\int_{\mathbb{R}}dk~{}e^{\iota 2\pi
kx}}_{\mathcal{F}_{1}^{-1}}$ $\displaystyle\underbrace{e^{\iota 2\pi k(\alpha
y+\beta z)}}_{\text{phase shift}}\underbrace{\int_{\mathbb{R}}dx~{}e^{-\iota
2\pi kx}}_{\mathcal{F}_{1}}\xi(x,y,z)$ (13)
$\displaystyle=\underbrace{\mathcal{F}_{1}^{-1}\Phi_{23}(k\left|\alpha,\beta\right.)\mathcal{F}_{1}}_{\text{off-
diagonal shear}}\xi(x,y,z)$ (14)

Here, $\Phi_{23}(k\left|\alpha,\beta\right.)\equiv\exp[\iota 2\pi k(\alpha
y+\beta z)]$ is the phase ramp applied to axes 2 and 3 ( in this case, $y$ and
$z$ respectively). Similar phase ramps $\Phi_{31}(k|\alpha,\beta)$ and
$\Phi_{12}(k|\alpha,\beta)$ are constructed with the 1D Fourier transforms
$\mathcal{F}_{2}$ and $\mathcal{F}_{3}$ to achieve shears along the $zx$\- and
$xy$-planes. These are given by:

$\displaystyle\Phi_{31}(k|\alpha,\beta)$ $\displaystyle=\exp[\iota 2\pi(\alpha
z+\beta x)]$ (15) $\displaystyle\Phi_{12}(k|\alpha,\beta)$
$\displaystyle=\exp[\iota 2\pi(\alpha x+\beta y)]$ (16)

The use of the sampling matrix $\boldsymbol{B}_{\text{real}}$ from Ref. [21]
lends itself easily to the compact formulation of the shear operator in Eq.
(14). To see this, we consider the desired sampling step along
$\hat{\boldsymbol{k}}_{1}$, namely $B_{11}l+B_{12}m+B_{13}n$ (where
$\boldsymbol{m}\equiv[l~{}m~{}n]^{T}\in\mathbb{Z}^{3}$). This may be obtained
from $\psi^{\prime}_{\boldsymbol{m}}$ in Eq. (12) by:

$\displaystyle\psi^{\prime\prime}_{\boldsymbol{m}}\equiv\psi^{\prime}([B_{11}~{}$
$\displaystyle B_{22}~{}B_{33}]^{T}\boldsymbol{m})\longrightarrow$
$\displaystyle\psi^{\prime}([B_{11}l+B_{12}m+B_{13}n~{}B_{22}~{}B_{33}]^{T}\boldsymbol{m})$
$\displaystyle=\underbrace{\left[\mathcal{F}_{1}^{-1}\Phi_{23}\left(k\left|\frac{B_{12}}{B_{11}},\frac{B_{13}}{B_{11}}\right.\right)\mathcal{F}_{1}\right]}_{\text{$yz$
shear operator}}\psi^{\prime}_{\boldsymbol{m}}$ (17)

where the $yz$-shear operator acts not on the original objet
$\psi(s_{0}\boldsymbol{m})$, but the _bulk-resampled_ object
$\psi^{\prime}_{\boldsymbol{m}}$ from Eq. (12). Consequently, the complete
shear operation is given by the composition of the shear operators along the
independent Cartesian directions:

$\displaystyle\psi^{\prime\prime}_{\boldsymbol{m}}$
$\displaystyle\equiv\psi(\boldsymbol{B}_{\text{real}}\boldsymbol{m})$
$\displaystyle=\left[\mathcal{F}_{3}^{-1}\Phi_{12}\left(k\left|\frac{B_{31}}{B_{33}},\frac{B_{32}}{B_{33}}\right.\right)\mathcal{F}_{3}\right]\circ\left[\mathcal{F}_{2}^{-1}\Phi_{31}\left(k\left|\frac{B_{23}}{B_{22}},\frac{B_{21}}{B_{22}}\right.\right)\mathcal{F}_{2}\right]\circ\left[\mathcal{F}_{1}^{-1}\Phi_{23}\left(k\left|\frac{B_{12}}{B_{11}},\frac{B_{13}}{B_{11}}\right.\right)\mathcal{F}_{1}\right]\psi^{\prime}_{\boldsymbol{m}}$
(18)

where $\circ$ denotes operator composition. Fig. 19 shows the effect of the
shear resampling operation in Eq. (18) applied to the bulk-resampled object
(Fig. 18, row 2). Once again we note the fidelity to the original
discontinuity, which would have not been possible with a simple real-space
interpolation within the space of a single voxel.

Figure 19:  Row 1: amplitude cross-section of the _bulk-resampled_ digital
object in Fig. 18, before application of the shearing operator in Eq. (18).
The arrows denote the direction of shearing. The $\otimes$ and $\odot$ symbols
denote arrows entering and exiting the plane of the figure respectively. Row
2: shear-resampled object amplitude cross-sections. Row 3,4: corresponding
phase profiles, indicating the screw dislocation. All phase plots are in
radians.

Eq. (12) and Eq. (18) together denote the sequence of analytic operations that
take the unknown object $\psi(s_{0}\boldsymbol{m})$ to the appropriately
sampled object $\psi(\boldsymbol{B}_{\text{real}}\boldsymbol{m})$, for direct
use in Eq. (3). The complete resampling operation
$\psi(s_{0}\boldsymbol{m})\longrightarrow\psi(\boldsymbol{B}_{\text{real}}\boldsymbol{m})$
_differentiably_ transforms the original object (Fig. 18, row 1) to the
appropriately resampled object (Fig. 19, row 2), while preserving the fidelity
to high-frequency features within the object (in this case, phase
discontinuities).

###  5.3 Resampling by grid rotation

Building upon the shear-resampling methodology in Section 5.2, we describe our
last fundamental resampling operation in which the discrete grid of the object
is arbitrarily rotated. This is of crucial importance in the MR-BCDI problem,
since the scatterer has to be rotated into different Bragg conditions in turn.
The method decomposes the desired rotation into a sequence of shear
operations, and therefore has all the advantages of the shear-resampling
method, including differentiability and maximum fidelity to small features.

As an example, we wish to actively rotate the discretized object
$\psi(s_{0}\boldsymbol{m})$ by an angle $\theta$ about the axis
$+\hat{\boldsymbol{k}}_{3}$ (in a right-handed manner). Practically, we rotate
the Cartesian sampling grid by $-\theta$ about $+\hat{\boldsymbol{k}}_{3}$.
According to Refs. [48, 49, 50], Eq. (14) may be employed to achieve this by
three successive shear operations:

$\displaystyle\mathcal{R}(\theta,\hat{\boldsymbol{k}}_{3})\psi(\boldsymbol{x})$
$\displaystyle\equiv\psi\left(x\cos\theta+y\sin\theta,-x\sin\theta+y\cos\theta,z\right)$
$\displaystyle=\left[\mathcal{F}_{1}^{-1}\Phi_{12}\left(k\left|0,\tan\frac{\theta}{2}\right.\right)\mathcal{F}_{1}\right]\circ\left[\mathcal{F}_{2}^{-1}\Phi_{12}\left(k\left|-\sin\theta,0\right.\right)\mathcal{F}_{2}\right]\circ\left[\mathcal{F}_{1}^{-1}\Phi_{12}\left(k\left|0,\tan\frac{\theta}{2}\right.\right)\mathcal{F}_{1}\right]\psi(\boldsymbol{x})$
(19)

Similar formulations to Eq. (19) hold for rotations
$\mathcal{R}(\theta,\hat{\boldsymbol{k}}_{1})$ and
$\mathcal{R}(\theta,\hat{\boldsymbol{k}}_{2})$ about the
$\hat{\boldsymbol{k}}_{1}$ and $\hat{\boldsymbol{k}}_{2}$ directions. Using
this method, we can express the grid resampling of an _arbitrary_ rotation
$\mathcal{R}(\theta,\hat{\boldsymbol{n}})$ as a composition of three
rotations:

$\mathcal{R}(\theta,\hat{\boldsymbol{n}})=\mathcal{R}(\alpha_{3},\hat{\boldsymbol{k}}_{3})\circ\mathcal{R}(\alpha_{2},\hat{\boldsymbol{k}}_{2})\circ\mathcal{R}(\alpha_{1},\hat{\boldsymbol{k}}_{1})$
(20)

Here the $\alpha_{i}$ denote the Euler angles of rotation in the $XYZ$
convention (_i.e._ , rotation about $X$ followed by $Y$ and then $Z$). We note
that regardless of the Euler convention, an arbitrary rotation may be achieved
by a triplet of successive rotations about the principal Cartesian axes.

The rotation of the example crystal from the previous examples by an angle of
$45^{\circ}$ about the $Y$-axis is shown in Fig. 20. Again we note the
fidelity the resampled object to the discontinuity in the bulk. The Appendix
Section F contains an example that contrasts Fourier-based rotation resampling
with regular grid interpolation in a crystal containing multiple bulk
discontinuities.

Figure 20:  Row 1: Mutually orthogonal cross sections of the original
scatterer amplitude. The arrows indicate the direction of the $45^{\circ}$
rotation about the $Y$-axis viewed from different perspectives. Row 2:
scatterer amplitude after it has been rotated within the grid. Row 3, 4:
corresponding complex phase cross sections. All phase plots are in radians.

We note here in our global optimization scheme, $\mathcal{A}$ and
$\boldsymbol{u}$ are defined and reconstructed only within an estimated
bounding box $\mathcal{V}$ of the object. However, the transformations
described in this section are applied to a zero-padded version of
$\mathcal{V}$. In this paper, all buffered arrays are of size $128\times
128\times 128$. While our interpolation approach ensures that phase
discontinuities are well-preserved, we note that object amplitude suffers from
Gibbs phenomenon -induced oscillations at sharp edges and discontinuities.
This is evident on comparing the first and second rows in Figures 18 and 19.
However, this problem is automatically addressed in our global optimization
scheme through the use of the TV regularizer. As we have shown, this results
in smoothed edges in the reconstructions.

Further, we note that a large rotation angle about a principal axis (for
example, $\theta$ about the $Z$-axis from Eq. (19)), when implemented with the
fast Fourier transform, potentially causes the transformed object to be split
across the periodic array boundary, due to an excessively large shear. Whether
or not this splitting takes place depends on the size of $\theta$, and the
ratio of the crystal-to-buffer sizes in each array dimension. A smaller ratio
implies a larger $\theta$ can be accommodated before splitting commences. This
undesirable effect can be addressed by performing a large principal rotation
in a sequence of smaller rotations. In this work, we have conservatively set
this smaller rotation size to $45^{\circ}$. With this convention, a rotation
of $210^{\circ}$ about the $Z$-axis would be decomposed as four rotations of
$45^{\circ}$ followed by a $30^{\circ}$ rotation. While this certainly
lengthens the chain of analytic computations for each optimization step, we
have found this to introduce negligible overhead when implemented on the GPU.

##  Appendix A Mutual orthogonality

Consider a $3\times 3$ matrix
$\boldsymbol{M}\equiv[\boldsymbol{v}_{1}~{}\boldsymbol{v}_{2}~{}\boldsymbol{v}_{3}]$
whose columns represent 3 linearly independent vectors spanning 3D Euclidean
space $\mathbb{R}^{3}$. Here each of the $\boldsymbol{v}_{i}\in\mathbb{R}^{3}$
is a $3\times 1$ column vector. The _mutual orthogonality_ (MO) of
$\boldsymbol{M}$ is defined as:

$\displaystyle\mathcal{O}(\boldsymbol{M})$
$\displaystyle\equiv\frac{\boldsymbol{v}_{1}\cdot\boldsymbol{v}_{2}\times\boldsymbol{v}_{3}}{\left\lVert\boldsymbol{v}_{1}\right\rVert\left\lVert\boldsymbol{v}_{2}\right\rVert\left\lVert\boldsymbol{v}_{3}\right\rVert}$
(21)
$\displaystyle=\frac{\det(\boldsymbol{M})}{\left\lVert\boldsymbol{v}_{1}\right\rVert\left\lVert\boldsymbol{v}_{2}\right\rVert\left\lVert\boldsymbol{v}_{3}\right\rVert}$
(22)

Here $\left\lVert\cdot\right\rVert$ denotes the Euclidean distance or the
$L_{2}$-norm. We see that $\mathcal{O}(\boldsymbol{M})$ is invariant under
even permutations of the $\boldsymbol{v}_{i}$’s and changes sign but not
magnitude under odd permutations. This follows simply from the same properties
of the determinant $\det(\boldsymbol{M})$. A positive (negative) sign of
$\mathcal{O}(\boldsymbol{M})$ indicates that the columns of $\boldsymbol{M}$
are right (left) -handed in order.

If the basis vectors are mutually orthogonal (_i.e._ ,
$\boldsymbol{v}_{1}\perp\boldsymbol{v}_{2}\perp\boldsymbol{v}_{3}$), then we
see that $\mathcal{O}(\boldsymbol{M})=\pm 1$, indicating a rectangular basis
for $\mathbb{R}^{3}$. If they all lie in the same plane,
$\mathcal{O}(\boldsymbol{M})=0$ and the ‘basis’ does not in fact span
$\mathbb{R}^{3}$. Thus, $\mathcal{O}(\boldsymbol{M})\in[-1,1]$ contains
information about the relative independence of the basis vectors of
$\boldsymbol{M}$, as well as their handedness.

We refer to Fig. 21.

Figure 21: Schematic elucidating angle parameterization of MO for a basis
matrix
$\boldsymbol{M}=[\boldsymbol{v}_{1}~{}\boldsymbol{v}_{2}~{}\boldsymbol{v}_{3}]$.

Here, $\mathcal{O}(\boldsymbol{M})$ can be parameterized in terms of the
angles between different vectors from the magnitudes of the cross and dot
products:

$\displaystyle\mathcal{O}(\boldsymbol{M})$
$\displaystyle=\frac{\left\lVert\boldsymbol{v}_{3}\right\rVert\left\lVert\boldsymbol{v}_{1}\times\boldsymbol{v}_{2}\right\rVert\cos\theta_{3}}{\left\lVert\boldsymbol{v}_{1}\right\rVert\left\lVert\boldsymbol{v}_{2}\right\rVert\left\lVert\boldsymbol{v}_{3}\right\rVert}$
(23)
$\displaystyle=\frac{\left\lVert\boldsymbol{v}_{3}\right\rVert\left(\left\lVert\boldsymbol{v}_{1}\right\rVert\left\lVert\boldsymbol{v}_{2}\right\rVert\sin\phi_{12}\right)\cos\theta_{3}}{\left\lVert\boldsymbol{v}_{1}\right\rVert\left\lVert\boldsymbol{v}_{2}\right\rVert\left\lVert\boldsymbol{v}_{3}\right\rVert}$
(24) $\displaystyle=\sin\phi_{12}\cos\theta_{3}$ (25)

We see that the MO of a basis of vectors does not depend on their magnitudes,
but only their relative orientations. This parameterization is useful in
determining the relation between the orthogonalities of mutually Fourier-
conjugate bases, as we see below.

The formulation of MO as given above is particularly useful in BCDI. The
measurement consists of rectilinear (but not orthogonal) samples of x-ray
intensity in reciprocal (Fourier) space, given by a $3\times 3$ sampling basis
matrix $\boldsymbol{B}_{\text{recip}}$ [21, 51]. In other words, Fourier space
is sampled at points $\boldsymbol{q}$ obtained from integer combinations of
the basis vectors.

$\displaystyle\boldsymbol{q}$
$\displaystyle=l\boldsymbol{v}_{1}+m\boldsymbol{v}_{2}+n\boldsymbol{v}_{3}$
(26) $\displaystyle=\boldsymbol{B}_{\text{recip}}\left[\begin{matrix}l\\\ m\\\
n\end{matrix}\right]$ (27) $\displaystyle\text{where }l,m,n$
$\displaystyle\in\mathbb{Z}.$

In practice, the integers $l,m,n$ range over $N_{1},N_{2},N_{3}$ values
respectively: $l=0,1,\ldots,N_{1}-1$, and so on for $m$ and $n$ as well. Thus,
the Fourier window is of size $N_{1}\times N_{2}\times N_{3}$ pixels.

The columns of $\boldsymbol{B}_{\text{recip}}$ are determined by the
diffraction geometry and the chosen manner of crystal rotation (‘rocking’).
This in turn unambiguously determines the conjugate basis
$\boldsymbol{B}_{\text{real}}$ whose columns denote the real-space samplings
steps of the object wave (_i.e._ , the object itself). The real-space object
sampled using $\boldsymbol{B}_{\text{real}}$ is related to the propagated wave
sampled in the far field using $\boldsymbol{B}_{\text{recip}}$ through the
discrete Fourier transform. The relation between the conjugate bases is given
by [21]:

$\displaystyle\boldsymbol{B}_{\text{real}}=\boldsymbol{B}_{\text{recip}}^{-T}{\underbrace{\left[\begin{matrix}N_{1}&&\\\
&N_{2}&\\\
&&N_{3}\end{matrix}\right]}_{\mathcal{D}=\text{diag}(N_{1},N_{2},N_{3})}}^{-1}$
(28)

where we have denoted $\mathcal{D}\equiv\text{diag}(N_{1},N_{2},N_{3})$. Here,
$(\cdot)^{-T}$ equivalently denotes the inverse of the transpose or the
transpose of the inverse. As we can see, a larger Fourier window (greater
$N_{1}$, $N_{2}$, $N_{3}$) results in a better real-space resolution _limit_
(smaller sampling steps).

The rotating crystal geometry in BCDI ensures that the diffraction signal is
always sampled not in a well-behaved, rectangular manner, but with a linear
shear. In other words,
$\left\lVert\mathcal{O}(\boldsymbol{B}_{\text{recip}})\right\rVert<1$. This
naturally introduces a shear into the real-space sampling basis
$\boldsymbol{B}_{\text{real}}$ through Eq. (28). In the light of this, we
would like to determine the relation between
$\mathcal{O}(\boldsymbol{B}_{\text{recip}})$ and
$\mathcal{O}(\boldsymbol{B}_{\text{real}})$. Ideally, we would prefer
$\mathcal{O}(\boldsymbol{B}_{\text{recip}})$ for each scan in a multi-
reflection BCDI measurement to be as close to unity as possible, and therefore
one must arrange their measurements accordingly. In our simulations, we choose
to exercise this flexibility rather coarsely through one of two rocking
directions, namely $\theta$ and $\phi$ from Figure 1. These correspond to
motor stages at the 34-ID-C coherent diffraction instrument at the Advanced
Photon Source.

We seek an expression for $\mathcal{O}(\boldsymbol{B}_{\text{real}})$ in terms
of $\mathcal{O}(\boldsymbol{B}_{\text{recip}})$. In order to do this, we first
define
$\boldsymbol{\beta}_{\text{recip}}\equiv\boldsymbol{B}_{\text{recip}}\mathcal{D}^{1/2}$
and
$\boldsymbol{\beta}_{\text{real}}\equiv\boldsymbol{B}_{\text{real}}\mathcal{D}^{1/2}$.
This implies that Eq. (28) can be rewritten more simply as:

$\boldsymbol{\beta}_{\text{real}}^{T}\boldsymbol{\beta}_{\text{recip}}=\mathds{1}$
(29)

where $\mathds{1}$ is the $3\times 3$ identity matrix. From the above
definitions, we see that $\boldsymbol{B}_{\text{real}}$
($\boldsymbol{B}_{\text{recip}}$) differs from
$\boldsymbol{\beta}_{\text{real}}$ ($\boldsymbol{\beta}_{\text{recip}}$) only
in the magnitudes of the columns, and therefore we have:
$\mathcal{O}(\boldsymbol{\beta}_{\text{recip}})=\mathcal{O}(\boldsymbol{B}_{\text{recip}})$
and
$\mathcal{O}(\boldsymbol{\beta}_{\text{real}})=\mathcal{O}(\boldsymbol{B}_{\text{real}})$.
It is therefore sufficient to determine the relation between
$\mathcal{O}(\boldsymbol{\beta}_{\text{recip}})$ and
$\mathcal{O}(\boldsymbol{\beta}_{\text{real}})$.

If
$\boldsymbol{\beta}_{\text{recip}}\equiv[\boldsymbol{k}_{1}~{}\boldsymbol{k}_{2}~{}\boldsymbol{k}_{3}]$
and
$\boldsymbol{\beta}_{\text{real}}\equiv[\boldsymbol{r}_{1}~{}\boldsymbol{r}_{2}~{}\boldsymbol{r}_{3}]$,
then Eq. (29) tells us that:

$\displaystyle\boldsymbol{r}_{l}$
$\displaystyle=\frac{\epsilon_{lmn}\boldsymbol{k}_{m}\times\boldsymbol{k}_{n}}{\det(\boldsymbol{\beta}_{\text{recip}})}=\frac{\epsilon_{lmn}\boldsymbol{k}_{m}\times\boldsymbol{k}_{n}}{\boldsymbol{k}_{1}\cdot\boldsymbol{k}_{2}\times\boldsymbol{k}_{3}}~{}\text{for
}l,m,n=1,2,3$ (30)
$\displaystyle\Longrightarrow\left\lVert\boldsymbol{r}_{l}\right\rVert$
$\displaystyle=\frac{1}{\left\lVert\boldsymbol{k}_{l}\right\rVert\cos\theta_{l}}$
(from the geometry in Fig. 21)

Here $\epsilon_{lmn}$ is the Levi-Civita symbol denoting the cyclic relation
between the indices $l$, $m$ and $n$. We therefore have:

$\displaystyle\mathcal{O}(\boldsymbol{\beta}_{\text{real}})$
$\displaystyle=\frac{\det\left(\boldsymbol{\beta}_{\text{real}}\right)}{\left\lVert\boldsymbol{r}_{1}\right\rVert\left\lVert\boldsymbol{r}_{2}\right\rVert\left\lVert\boldsymbol{r}_{3}\right\rVert}$
(31)
$\displaystyle=\frac{\left\lVert\boldsymbol{k}_{1}\right\rVert\left\lVert\boldsymbol{k}_{2}\right\rVert\left\lVert\boldsymbol{k}_{3}\right\rVert}{\det\left(\boldsymbol{\beta}_{\text{recip}}\right)}\cos\theta_{1}\cos\theta_{2}\cos\theta_{3}$
($\because\det\boldsymbol{\beta}_{\text{real}}=1/(\det\boldsymbol{\beta}_{\text{recip}})$
from Eq. (29)) $\displaystyle\Longrightarrow$
$\displaystyle\boxed{\mathcal{O}(\boldsymbol{\beta}_{\text{real}})=\frac{\cos\theta_{1}\cos\theta_{2}\cos\theta_{3}}{\mathcal{O}(\boldsymbol{\beta}_{\text{recip}})}}$
(32)

We note that a high MO in reciprocal space indicates a similarly high MO in
real space. These MOs respectively pertain to how the signal and the scatterer
are discretized. Hence choosing $\mathcal{O}(\boldsymbol{B}_{\text{recip}})$
as close to $1$ as possible by judicious choice of crystal rotation (in
particular, rocking about the $\theta$ or $\phi$ axis) facilitates the best
orthogonal discretization of the reconstructed object wave, owing to
interpolation from a minimally sheared sampling basis
$\boldsymbol{B}_{\text{real}}$.

##  Appendix B Silicon carbide nanoparticle fabrication

The SiC nanoparticles were fabricated from a single 4H-SiC wafer with an _i-p-
n_ structure procured from Norstel. The wafer comprised a $400$ nm-thick
intrinsic SiC layer with $<10^{15}$ cm-3 intrinsic defect content, a
$2~{}\mu$m-thick _p_ -type SiC layer with $<10^{19}$ cm-3 of aluminum atoms,
and a $0.5$ mm-thick _n_ -type handle with $<10^{18}$ cm-3 of nitrogen atoms.
The wafer was diced into $5\times 5$ mm2 chips, patterned by electron beam
lithography (bi-layer PMMA 495K-A3/950K-A6; $1000\mu$C cm-2 dose, 3-minute
development in 1:3 MIBK:IPA at room temperature), and metallized (10 nm of Ti
and 60 nm of Cr via e-beam evaporation) to transfer the asymmetric, D-shaped
arrays onto the SiC using plasma etching. During synthesis, the D shape flat
was aligned to the crystallographic $[11\bar{2}0]$ direction (see Figure
16(d)). The SiC was etched (SF6-Ar ICP/RIE dry-etching) with the Ti/Cr acting
as a hard mask. The $1.2~{}\mu$m etch depth exposed the _p_ -type layer for
subsequent photo-electrochemical (PEC) etching of the nanoparticles. The
plasma etch-defined nanopillars were then PEC etched (0.2M KOH, -0.3V biased
etch under 365nm UV illumination at $<500$ mW) to partially and selectively
remove the _p_ -type SiC to allow for nanoparticle formation and subsequent
detachment. For this particular experimental run, the _p_ -type SiC was
underetched, resulting in a more substantial, weakly scattering, porous _p_
-type tether remaining. This resulted in the nanoparticle having a larger
perceived size as seen in Figure 16(e). These partially undercut
nanopillar/nanoparticle arrays were stamp-transferred to a PMMA-coated silicon
substrate. The wafers were baked at the PMMA glass transition temperature,
which locked the nanoparticles on the Si substrate. The trasnferred
nanoparticles were arranged on the Si wafer in $100~{}\mu$m pitched arrays
bordered by macroscopic fiducials consisting of tightly packed nanoparticles.
Finally, the PMMA was removed by oxygen plasma, leaving pristine SiC on Si
nanoparticle arrays that were then conformally covered with alumina (22 nm
thick) via ALD.

##  Appendix C Phase retrieval recipe

Given below is the phase retrieval recipe used in all the reconstructions in
this paper. All these algorithms are reviewed in Ref. [19].

  1. 1.

$150$ iterations of ER, with shrinkwrap every $30$ iterations.

  2. 2.

$300$ iterations of HIO.

  3. 3.

$100$ iterations of solvent-flipping (SF), with shrinkwrap every $25$
iterations.

  4. 4.

$300$ iterations of HIO.

  5. 5.

$450$ iterations of ER, with shrinkwrap every $90$ iterations.

##  Appendix D Feature fidelity with exaggerated resmpling

Visual demonstration of feature fidelity with exaggerated under- and over-
sampling (Figure 22).

Figure 22:  Orthogonal cross sections of a discontinuous 3D phase object
subjected to uniform bulk transformation
$\mathcal{B}(\hat{\boldsymbol{k}}_{i}\left|s_{0}\rightarrow\xi
B_{ii}\right.)$. For simplicity, $s_{0}=B_{ii}=1$ for all $i=1,2,3$ (along all
three axes), with (top row) $\xi=1$ (original), (middle row) $\xi=1.5$ and
(bottom row) $\xi=0.75$.

##  Appendix E Sequential optimization plans

Table 1:  Sequential optimization plan for dislocation-free crystal. The BCDI scans within each minibatch were randomly selected.  Epoch | Number of minibatches | Scans in minibatch | Iterations per minibatch
---|---|---|---
1 | 400 | 2 | 6
2 | 320 | 3 | 12
3 | 240 | 3 | 25
4 | 160 | 3 | 50
5 | 80 | 3 | 100
6 | 40 | 4 | 200
7 | 20 | 4 | 400
8 | 1 | 5 (full dataset) | 2000
Table 2:  Optimization plan for crystal with orthogonal screw dislocations along $[111]$ and $[2\bar{2}0]$. The BCDI scans within each minibatch were randomly selected.  Epoch | Number of minibatches | Scans in minibatch | Iterations per minibatch
---|---|---|---
1 | 800 | 2 | 3
2 | 400 | 2 | 6
3 | 320 | 2 | 12
4 | 240 | 3 | 25
5 | 160 | 3 | 50
6 | 80 | 3 | 100
7 | 40 | 3 | 200
8 | 20 | 3 | 400
9 | 1 | 4 (full dataset) | 2000
Table 3:  Optimization plan for the silicon carbide nanocrystal. The BCDI scans within each minibatch were randomly selected.  Epoch | Number of minibatches | Scans in minibatch | Iterations per minibatch
---|---|---|---
1 | 800 | 4 | 6
2 | 640 | 4 | 12
3 | 480 | 4 | 25
4 | 320 | 5 | 50
5 | 160 | 5 | 100
6 | 80 | 5 | 200
7 | 40 | 5 | 400
8 | 1 | 6 (full dataset) | 2000

##  Appendix F Comparison of interpolation methods

Figure 23:  (a) 2D slice of 3D discontinuous phase object prior to resampling
by grid rotation. (b) Grid rotation by simple multi-linear interpolation
between neighboring voxels. (c) Grid rotation by Fourier transform -based
interpolation (Eq. (19)). The inset plots highlight the difference in fidelity
to small, pixel-scale features.

##  Appendix G Funding

The development of the analytical MR-BCDI model and the Fourier transform-
based interpolation method, proof-of-concept demonstration and design and
fabrication of the SiC nanoparticles was supported by the US Department of
Energy (DOE), Office of Science, Basic Energy Sciences, Materials Science and
Engineering Division. Additional support for materials preparation came from
the Q-NEXT Quantum Center, a U.S. Department of Energy, Office of Science,
National Quantum Information Science Research Center, under Award Number DE-
FOA-0002253. Silicon carbide deterministic nanoparticle fabrication and SEM
characterization work was performed under proposals 72483 and 775514 in the
Center for Nanoscale Materials clean room. Work performed at the Center for
Nanoscale Materials, a U.S. Department of Energy Office of Science User
Facility, was supported by the U.S. DOE, Office of Basic Energy Sciences,
under Contract No. DE-AC02-06CH11357. Refinement of the geometric,
computational and optimization concepts was supported by the European Research
Council (European Union’s Horizon H2020 research and innovation program grant
agreement No. 724881). Generation of the simulated structures and the BCDI
data acquisition was supported by the Laboratory Directed Research and
Development (LDRD) funding from Argonne National Laboratory, provided by the
Director, Office of Science, of the U.S. Department of Energy under Contract
No. DE-AC02-06CH11357. This research uses the resources of the Advanced Photon
Source, a US DOE Office of Science User Facility operated for the DOE Office
of Science by Argonne National Laboratory under contract No. DE-
AC02-06CH11357.

##  Appendix H Acknowledgments

The authors gratefully acknowledge numerous valuable discussions with Drs.
Anthony Rollett, Robert Suter and Matthew Wilkin (Carnegie Mellon University),
Nicholas Porter and Dr. Richard Sandberg (Brigham Young University), Dr. Ross
Harder (Argonne National Laboratory) and Dr. Anastasios Pateras (DESY).

##  Appendix I Disclosures

The authors declare no conflicts of interest. Data and analysis scripts beyond
those available in the code repository will be made available upon reasonable
request.

## References

  * [1] IA Vartanyants, JA Pitney, JL Libbert, and IK Robinson.  Reconstruction of surface morphology from coherent x-ray reflectivity.  Physical Review B, 55(19):13193, 1997.
  * [2] I. K. Robinson, I. A. Vartanyants, G. J. Williams, M. A. Pfeifer, and J. A. Pitney.  Reconstruction of the shapes of gold nanocrystals using coherent x-ray diffraction.  Phys. Rev. Lett., 87:195505, Oct 2001.
  * [3] Jianwei Miao, Tetsuya Ishikawa, Ian K. Robinson, and Margaret M. Murnane.  Beyond crystallography: Diffractive imaging using coherent x-ray light sources.  Science, 348(6234):530–535, 2015.
  * [4] Andrew Ulvestad, Andrej Singer, Jesse Clark, Josh Cho, Jong Woo Kim, Ross Harder, Jorg Maser, Shirley Meng, and Oleg Shpyrko.  Topological defect dynamics in operando battery nanoparticles.  Science, 348(6241):1344–1347, 2015.
  * [5] A. Singer, M. Zhang, S. Hy, D. Cela, C. Fang, T. A. Wynn, B. Qiu, Y. Xia, Z. Liu, A. Ulvestad, N. Hua, J. Wingert, H. Liu, M. Sprung, A. V. Zozulya, E. Maxey, R. Harder, Y. S. Meng, and O. G. Shpyrko.  Nucleation of dislocations and their dynamics in layered oxide cathode materials during battery charging.  Nature Energy, 3(8):641–647, Aug 2018.
  * [6] Tomoya Kawaguchi, Thomas F. Keller, Henning Runge, Luca Gelisio, Christoph Seitz, Young Yong Kim, Evan R. Maxey, Wonsuk Cha, Andrew Ulvestad, Stephan O. Hruszkewycz, Ross Harder, Ivan A. Vartanyants, Andreas Stierle, and Hoydoo You.  Gas-induced segregation in pt-rh alloy nanoparticles observed by in situ bragg coherent diffraction imaging.  Phys. Rev. Lett., 123:246001, Dec 2019.
  * [7] H. Kraus, V. A. Soltamov, F. Fuchs, D. Simin, A. Sperlich, P. G. Baranov, G. V. Astakhov, and V. Dyakonov.  Magnetic field and temperature sensing with atomic-scale spin defects in silicon carbide.  Scientific Reports, 4(1):5303, Jul 2014.
  * [8] S. O. Hruszkewycz, S. Maddali, C. P. Anderson, W. Cha, K. C. Miao, M. J. Highland, A. Ulvestad, D. D. Awschalom, and F. J. Heremans.  Strain annealing of sic nanoparticles revealed through bragg coherent diffraction imaging for quantum technologies.  Phys. Rev. Materials, 2:086001, Aug 2018.
  * [9] R. W. Gerchberg.  A practical algorithm for the determination of the phase from image and diffraction plane pictures.  Optik, 35:237–246, 1972.
  * [10] J. R. Fienup.  Phase retrieval algorithms: a comparison.  Appl. Opt., 21(15):2758–2769, Aug 1982.
  * [11] S. Marchesini, H. He, H. N. Chapman, S. P. Hau-Riege, A. Noy, M. R. Howells, U. Weierstall, and J. C. H. Spence.  X-ray image reconstruction from a diffraction pattern alone.  Phys. Rev. B, 68:140101, Oct 2003.
  * [12] Marcus C. Newton.  Concurrent phase retrieval for imaging strain in nanocrystals.  Phys Rev B, 102:014104, Jul 2020.
  * [13] Yuan Gao, Xiaojing Huang, Hanfei Yan, and Garth J. Williams.  Bragg coherent diffraction imaging by simultaneous reconstruction of multiple diffraction peaks.  Phys Rev B, 103:014102, Jan 2021.
  * [14] Felix Hofmann, Nicholas W. Phillips, Suchandrima Das, Phani Karamched, Gareth M. Hughes, James O. Douglas, Wonsuk Cha, and Wenjun Liu.  Nanoscale imaging of the full strain tensor of specific dislocations extracted from a bulk sample.  Phys. Rev. Materials, 4:013801, Jan 2020.
  * [15] Matthew J. Wilkin, Siddharth Maddali, Stephan O. Hruszkewycz, Anastasios Pateras, Richard L. Sandberg, Ross Harder, Wonsuk Cha, Robert M. Suter, and Anthony D. Rollett.  Experimental demonstration of coupled multi-peak bragg coherent diffraction imaging with genetic algorithms.  Phys. Rev. B, 103:214103, Jun 2021.
  * [16] David Saad.  Online algorithms and stochastic approximations.  Online Learning, 5:6–3, 1998.
  * [17] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of machine learning research, 12(7), 2011.
  * [18] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [19] S. Marchesini.  Invited article: A unified evaluation of iterative projection algorithms for phase retrieval.  Review of Scientific Instruments, 78(1):011301, 2007.
  * [20] Pierre Godard, Marc Allain, Virginie Chamard, and John Rodenburg.  Noise models for low counting rate coherent diffraction imaging.  Optics express, 20(23):25914–25934, 2012.
  * [21] S. Maddali, P. Li, A. Pateras, D. Timbie, N. Delegan, A. L. Crook, H. Lee, I. Calvo-Almazan, D. Sheyfer, W. Cha, F. J. Heremans, D. D. Awschalom, V. Chamard, M. Allain, and S. O. Hruszkewycz.  General approaches for shear-correcting coordinate transformations in Bragg coherent diffraction imaging. Part I.  Journal of Applied Crystallography, 53(2), Apr 2020.
  * [22] S. Maddali, M. Allain, P. Li, V. Chamard, and S. O. Hruszkewycz.  Detector tilt considerations in bragg coherent diffraction imaging: A simulation study.  Crystals, 10(12), 2020.
  * [23] Dominik Kriegner, Eugen Wintersberger, and Julian Stangl.  xrayutilities: a versatile tool for reciprocal space conversion of scattering data recorded with linear and area detectors.  Journal of Applied Crystallography, 46(4):1162–1170, Aug 2013.
  * [24] Kandel, Saugat.  Using automatic differentiation for coherent diffractive imaging applications.  PhD thesis, 2021.
  * [25] J. R. Fienup and C. C. Wackerman.  Phase-retrieval stagnation problems and solutions.  J. Opt. Soc. Am. A, 3(11):1897–1907, Nov 1986.
  * [26] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [27] S. Maddali.  _mrbcdi_ : Differentiable multi-reflection bragg coherent diffraction imaging for lattice distortion fields in crystals (github), 2022.
  * [28] Henry N. Chapman, Anton Barty, Stefano Marchesini, Aleksandr Noy, Stefan P. Hau-Riege, Congwu Cui, Malcolm R. Howells, Rachel Rosen, Haifeng He, John C. H. Spence, Uwe Weierstall, Tobias Beetz, Chris Jacobsen, and David Shapiro.  High-resolution ab initio three-dimensional x-ray diffraction microscopy.  J. Opt. Soc. Am. A, 23(5):1179–1200, May 2006.
  * [29] S. Maddali.  _Phaser_ : Python-based bcdi phase retrieval for cpu and gpu (github), 2020.
  * [30] Rob Phillips.  Crystals, Defects and Microstructures: Modeling Across Scales.  Cambridge University Press, 2001.
  * [31] A. Ulvestad, M. Menickelly, and S. M. Wild.  Accurate, rapid identification of dislocation lines in coherent diffractive imaging via a min-max optimization formulation.  AIP Advances, 8(1):015114, 2018.
  * [32] Ke Xu, Liang-Liang Niu, Shuo Jin, Xiaolin Shu, Hongxian Xie, Lifang Wang, and Guang-Hong Lu.  Atomistic simulations of screw dislocations in bcc tungsten: From core structures and static properties to interaction with vacancies.  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 393:174–179, 2017.  Computer Simulation of Radiation effects in Solids Proceedings of the 13 COSIRES Loughborough, UK, June 19-24 2016.
  * [33] A. Ulvestad, Y. Nashed, G. Beutier, M. Verdier, S. O. Hruszkewycz, and M. Dupraz.  Identifying defects with guided algorithms in bragg coherent diffractive imaging.  Scientific Reports, 7(1):9920, Aug 2017.
  * [34] F. Fuchs, V. A. Soltamov, S. Väth, P. G. Baranov, E. N. Mokhov, G. V. Astakhov, and V. Dyakonov.  Silicon carbide light-emitting diode as a prospective room temperature source for single photons.  Scientific Reports, 3(1):1637, Apr 2013.
  * [35] A. Lohrmann, N. Iwamoto, Z. Bodrog, S. Castelletto, T. Ohshima, T. J. Karle, A. Gali, S. Prawer, J. C. McCallum, and B. C. Johnson.  Single-photon emitting diode in silicon carbide.  Nature Communications, 6(1):7783, Jul 2015.
  * [36] Corey J. Cochrane, Jordana Blacksberg, Mark A. Anders, and Patrick M. Lenahan.  Vectorized magnetometer for space applications using electrical readout of atomic scale defects in silicon carbide.  Scientific Reports, 6(1):37077, Nov 2016.
  * [37] Mete Atatüre, Dirk Englund, Nick Vamivakas, Sang-Yun Lee, and Joerg Wrachtrup.  Material platforms for spin-based photonic quantum technologies.  Nature Reviews Materials, 3(5):38–51, May 2018.
  * [38] Kevin C. Miao, Alexandre Bourassa, Christopher P. Anderson, Samuel J. Whiteley, Alexander L. Crook, Sam L. Bayliss, Gary Wolfowicz, Gergo Thiering, Peter Udvarhelyi, Viktor Ivady, Hiroshi Abe, Takeshi Ohshima, Adam Gali, and David D. Awschalom.  Electrically driven optical interferometry with spins in silicon carbide.  Science Advances, 5(11):eaay0527, 2019.
  * [39] Alexandre Bourassa, Christopher P. Anderson, Kevin C. Miao, Mykyta Onizhuk, He Ma, Alexander L. Crook, Hiroshi Abe, Jawad Ul-Hassan, Takeshi Ohshima, Nguyen T. Son, Giulia Galli, and David D. Awschalom.  Entanglement and control of single nuclear spins in isotopically engineered silicon carbide.  Nature Materials, 19(12):1319–1325, Dec 2020.
  * [40] John A Nelder and Roger Mead.  A simplex method for function minimization.  The computer journal, 7(4):308–313, 1965.
  * [41] Pierre Godard.  On the use of the scattering amplitude in coherent x-ray bragg diffraction imaging.  Journal of Applied Crystallography, 54(3):797–802, 2021.
  * [42] S. Kandel, S. Maddali, M. Allain, S. O. Hruszkewycz, C. Jacobsen, and Y. Nashed.  Using automatic differentiation as a general framework for ptychographic reconstruction.  Opt. Express, 27(13):18653–18672, Jun 2019.
  * [43] Saugat Kandel, S. Maddali, Youssef S. G. Nashed, Stephan O. Hruszkewycz, Chris Jacobsen, and Marc Allain.  Efficient ptychographic phase retrieval via a matrix-free levenberg-marquardt algorithm.  Opt. Express, 29(15):23019–23055, Jul 2021.
  * [44] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup.  Efficient subpixel image registration algorithms.  Opt. Lett., 33(2):156–158, Jan 2008.
  * [45] Alexander Scheinker and Reeju Pokharel.  Adaptive 3d convolutional neural network-based reconstruction method for 3d coherent diffraction imaging.  Journal of Applied Physics, 128(18):184901, November 2020.
  * [46] Mathew J. Cherukara, Tao Zhou, Youssef Nashed, Pablo Enfedaque, Alex Hexemer, Ross J. Harder, and Martin V. Holt.  AI-enabled high-resolution scanning coherent diffraction imaging.  Applied Physics Letters, 117(4):044103, July 2020.
  * [47] Ross Harder.  Deep neural networks in real-time coherent diffraction imaging.  IUCrJ, 8(1):1–3, January 2021.
  * [48] Michael Unser, Philippe Thevenaz, and Leonid Yaroslavsky.  Convolution-based interpolation for fast, high-quality rotation of images.  IEEE Transactions on image processing, 4(10):1371–1381, 1995.
  * [49] Kieran G Larkin, Michael A Oldfield, and Hanno Klemm.  Fast fourier method for the accurate rotation of sampled images.  Optics communications, 139(1-3):99–106, 1997.
  * [50] Philippe Thévenaz, Thierry Blu, and Michael Unser.  Interpolation revisited [medical images application].  IEEE Transactions on medical imaging, 19(7):739–758, 2000.
  * [51] P. Li, S. Maddali, A. Pateras, I. Calvo-Almazan, S.O. Hruszkewycz, W. Cha, V. Chamard, and M. Allain.  General approaches for shear-correcting coordinate transformations in Bragg coherent diffraction imaging. Part II.  Journal of Applied Crystallography, 53(2), Apr 2020.