Accurate, rapid identification of dislocation lines in coherent diffractive imaging via a min-max optimization formulation

Defects such as dislocations impact materials properties and their response during external stimuli. Defect engineering has emerged as a possible route to improving the performance of materials over a wide range of applications, including batteries, solar cells, and semiconductors. Imaging these defects in their native operating conditions to establish the structure-function relationship and, ultimately, to improve performance has remained a considerable challenge for both electron-based and x-ray-based imaging techniques. However, the advent of Bragg coherent x-ray diffractive imaging (BCDI) has made possible the 3D imaging of multiple dislocations in nanoparticles ranging in size from 100 nm to1000 nm. While the imaging process succeeds in many cases, nuances in identifying the dislocations has left manual identification as the preferred method. Derivative-based methods are also used, but they can be inaccurate and are computationally inefficient. Here we demonstrate a derivative-free method that is both more accurate and more computationally efficient than either derivative- or human-based methods for identifying 3D dislocation lines in nanocrystal images produced by BCDI. We formulate the problem as a min-max optimization problem and show exceptional accuracy for experimental images. We demonstrate a 260x speedup for a typical experimental dataset with higher accuracy over current methods. We discuss the possibility of using this algorithm as part of a sparsity-based phase retrieval process. We also provide the MATLAB code for use by other researchers.

In materials, crystallographic imperfections such as dislocations often dictate performance and properties. For example, dislocation cores can act as fast diffusion sites [1][2][3][4], mitigate strain and plasticity during structural phase transformations [5][6][7], and govern crystal growth [8][9][10]. Increasingly, Bragg coherent x-ray diffractive imaging (BCDI) is being utilized at synchrotron and x-ray free electron laser [11,12] sources to address this challenge of understanding and optimizing materials properties via tuning of lattice distortions by nondestructively imaging the 3D lattice distortion field under in situ and operando conditions [13][14][15][16][17][18][19]. The typical resolution of the technique is roughly 10 −4 in strain sensitivity and 10-20 nm spatial resolution at 5-10 minute temporal resolution. Recent studies have revealed the 3D dislocation line dynamics in an individual nanocrystal during reactive processes such as hydrogen uptake [5,6], battery charging [20], crystal growth and dissolution [9,21], and grain growth in polycrystalline materials [22]. In all these studies, accurate identification of the dislocation line in the reconstructed image is essential to understanding the underlying physics and the dislocation impact. While these studies have shown great promise, the breadth of 3D BCDI dislocation dynamics measurements and techniques could expand substantially if accurate, robust, and rapid methods existed to determine the 3D dislocation line structure. For example, the datasets generated at diffraction-limited storage rings will likely be too large for existing derivative-based and human-in-the-loop-based methods [23], and the "dislocation basis" could potentially be used as a sparse basis to circumvent constraints in phase retrieval [24]. Here we present such a method by reformulating the dislocation core identification problem as a min-max optimization problem that can be solved rapidly.
It is counterintuitive that an imaging experiment with 10-20 nm spatial resolution is sensitive to atomic-scale defects such as dislocations. To understand how this is possible, consider the relationship between the continuum representation of the crystal, ρ(r), and the diffraction intensity, I(q) in the far field under a perfectly coherent illumination and in the kinematical scattering approximation [13,25,26]: Here, r and q are the real and reciprocal space coordinates, respectively, F is the Fourier transform that provides the map between these two spaces, Q is the measured reciprocal lattice vector (e.g., the 111 Bragg peak for a face-centered cubic crystal lattice), and u(r) is the vector displacement field that is a continuum description of how the atoms are displaced from their equilibrium positions. If we consider a cubic crystal with a screw dislocation along the z direction, then the displacement field is given by u x = u y = 0 and where b is the Burgers vector and θ measures the angle around the dislocation core [27]. If the Bragg peak (reciprocal lattice vector) measured is the 001 peak, then Q · u(r) = 2π θ. The Burgers vector in this case is equal to the lattice spacing in the z direction, and so Q · u(r) = θ. Thus, the signature of a dislocation in the complex image is a point around which the phase varies from 0 to 2π, or −π to π depending on the chosen convention. This is the maximum "signal" possible in the image phase and is a direct consequence of dislocations introducing large displacement fields. Note that this argument can be extended to other Bragg peaks that are not parallel to the displacement field induced by the dislocation.
To demonstrate a BCDI study of a dislocation, we show in The 001 Bragg peak is used after Gaussian blurring to simulate finite resolution and, assuming perfect phase retrieval, to reconstruct a BCDI image of the crystal. The particle shape is represented as an isosurface of the amplitude, which is proportional to the Bragg electron density, while the color projected onto the shape shows the atomic displacement field. (c) The 2D amplitude cross-section corresponding to the location shown in (b). The dislocation core signature can be seen as a low-amplitude region. (d) The 2D phase crosssection; the dislocation manifests itself as a phase vortex, or a region where the phase varies from −π to π. Note that the dislocation core is unique but the phase discontinuity is not. To see why, consider (e) and (f), which have two different global (constant) phase offsets applied; (e) shows a global phase offset of e iπ/2 while (f) shows a global phase offset of e iπ . Both images produce the same diffraction pattern according to Eq. 1. Thus, the reconstructed image is not sensitive to the absolute phase. Different global phase offsets move the phase discontinuity around the dislocation core. This fact is used in the derivative-based dislocation identification method. Fig. 1 the atomic and BCDI description of a screw dislocation in a cubic crystal. A periodic array of atoms is created to fill a particular particle shape, and the screw dislocation displacement field is applied (Fig. 1a). The 001 Bragg peak, blurred by a Gaussian function to simulate a finite resolution, is used to reconstruct the particle image at 15 nm resolution assuming perfect phase retrieval [28][29][30][31]. The reconstructed BCDI image is complex. The amplitude corresponds to the Bragg, or diffracting, electron density [32], while the phase corresponds to the displacement of the atoms from their equilibrium positions. The amplitude is used to draw the isosurface while the phase is projected onto the isosurface as the colormap (Fig. 1b). While the atoms cannot be identified in Fig. 1b because of the finite resolution, there is a key signature in the phase in the form of a phase vortex, or a region where the phase varies from −π to 0 to π. The 2D cross-sections through the center of the 3D image are shown in Fig. 1c-f. There is a signature of the dislocation core in the amplitude cross-section (Fig. 1c), but this is often difficult to distinguish from the spatial amplitude variation in real experimental data (see, e.g., Fig. 3b). The dislocation signature in the phase is much stronger (Fig. 1d). Again, the reason is that dislocations are large displacements of the atoms from their equilibrium values. However, care must be taken in interpreting the displacement field image. In this case the dislocation core is unique. However, if periodicity is not appropriately accounted for, then the spatial location of the displacement field discontinuity, defined as the jump from −π to π, is not unique, and so this discontinuity cannot be used directly for dislocation identification. To understand why, consider applying a global (constant) phase offset of e iπ/2 to the reconstructed image (Fig. 1e). The diffraction pattern does not change because the factor is a constant and |e iφ0 | = 1 for all φ 0 (see Eq. 1). However, the phase discontinuity line shifts around the dislocation core. An additional example is shown in Fig. 1f for a different phase offset. Thus, the phase discontinuity's spatial location is not unique, but the dislocation core's spatial location is unique. All phase offsets shift the discontinuity around the dislocation core, which is exploited in the standard derivative-based method for dislocation core identification described next.
In previous work, dislocation cores were identified by eye or with a derivative-based method that considered a range of phase offsets as shown in Figs. 1d-f. The derivative-based method was previously detailed in [5,9]; only the key steps are repeated here. A derivative of the displacement field is taken in all three orthogonal image directions, thresholding is applied to determine what constitutes a large derivative value, and the locations of these large derivatives are stored. The process is repeated for 360 different global phase offsets in 1 • steps. The intersection (across these 360 offsets) of the locations with large derivatives is used to identify the dislocation core. This process is computationally and memory intensive (three derivatives for each voxel for all 360 phase offsets need to be computed and stored) and requires tuning multiple thresholds. We now define the new algorithm and show how it is both more efficient and more accurate.
We denote byx the integer-valued spatial coordinates of a pixel for which we have obtained the complex-valued signal from phase retrieval, where we denote by ρ(x) and θ(x) the amplitude and phase values at that pixel, respectively. We further assume that, in preprocessing, a binary mask has been applied to θ so that a set of pixelsx with low amplitude (i.e., |ρ(x)| < τ for some specified threshold τ > 0) are not considered. We denote the set of all masked pixelsx by M and its set complement byM . Such a mask has, for example, been applied to Figs. 3b-c to mask the pixels that have nearzero amplitude. For arbitraryx and α ≥ 0, denote the αneighborhood ofx by B α (x) = {ŷ : x −ŷ ∞ ≤ α}. Consider the α-parameterized function of phase offset φ 0 and spatial coordinatesx given by where we have used ≡ to denote modulo (e.g., [2π + 1 ≡ 2π] = 1). We adopt the convention that if B α (x) ∩M = ∅ (i.e., B α (x) ⊂ M ), then for any value of φ 0 , f α (x, φ 0 ) = 0. Based on our previous discussion, large-magnitude values of the function signal thatx is potentially at or near a dislocation core, since Ψ α (x) represents the least possible length (in radians) of a radius 1 arc containing all the phases, modulo 2π, of pixels in the neighborhood B α (x). We are thus interested in solving the min-max problem We note that Eq. 3 is posed as a max-min problem, but we conform with the standard optimization terminology of "minmax" since any max-min problem has an equivalent min-max formulation obtained by negating the objectives. The interpretation of Eq. 3 is that it seeks the largest values of phase differences under the best possible value of the phase offset φ 0 . Because φ 0 is a continuous parameter andx is a discrete pixel location, Eq. 3 is a potentially challenging mixed-integer nonlinear robust optimization problem [33]. In our case, however, Ψ α (x) is straightforward to evaluate. We employ the finite list of the phases associated with all unmasked neighbors of a pixelx. If Θ(x) is empty or composed of a single element, or ifx ∈ M , then we follow the convention that Ψ α (x) = 0. Provided thatx is not masked and that Θ(x) has at least two elements, we run the following procedure.

Return the least entry of A.
In practice, we typically set α = 1, as in Fig. 2. We remark that in two dimensions, there are at most 9 phases in B 1 (x) ∩M ; in three dimensions, there are at most 27 phases in B 1 (x) ∩M . Combining this information with the fact that anyx ∈ M satisfies Ψ α (x) = 0, we can evaluate Ψ α (x) at every unmasked pixel in a given 2D or 3D dataset in time that is linear in the number of unmasked pixels. The min-max problem in Eq. 3 is then trivially solved by taking the maximum of Ψ α (x) over all such pixelsx, a task that can be performed in an online fashion (i.e., without explicitly storing the Ψ α (x) values).
We can gain insight into the algorithm by considering the three cases shown in Fig. 2. Figure 2a shows the sorted experimental phase values in the neighborhood of a pixel in a region of minimal phase variation. The MATLAB range has been shifted from [−π, π] to [0, 2π] by applying +π to all values. These values were all close to zero in the reconstructed image, and thus they are all close to π in the plot. Then, we compute 2π minus the pairwise differential (A i in step 2). This computation leads to values close to 2π for most of the pairs, since the differences are near zero. When computing the difference between the last point and the first point, the value is simply the difference, not 2π minus the difference (A n in step 2). The minimum of this set of values (step 3) is then close to zero, since point 27 and point 1 are close in value. This point is thus not likely to be near a dislocation core. Figure 2b shows the sorted experimental phase values in the neighborhood of a pixel in the dislocation core region. Again, we compute 2π minus the pairwise differential (A i in step 2) and the pairwise differential for point 27 and point 1 (A n in step 2). In this case, the minimum (step 3) in the 2π minus differential list comes from the difference between point 18 and point 17 (a difference of roughly 2), and is thus roughly 4.2. Pixel values above π/2 indicate potential proximity to a dislocation core, with this likelihood growing as the values approach the 2π upper bound on Ψ α . Fig. 2c shows the sorted experimental phase values for a pixel outside of, but close to, the phase discontinuity region. In this case, the minimum in the 2π minus differential list comes from the difference between point 9 and point 10. The differential is approximately 3π/2; since 2π − 3π/2 = π/2, the minimum value is approximately 1.62.
We compare our new method with the derivative-based method for identifying dislocations in images reconstructed from experimental data in Fig. 3. Figure 3a shows the experimental reconstruction of a silver nanoparticle after a dissolu- Identifying the dislocation line in images reconstructed from experimental data in which a nanocrystal underwent electrochemicalinduced dissolution. (a) The particle shape represented as an isosurface of the amplitude, which is proportional to the Bragg electron density, while the color projected onto the shape shows the atomic displacement field. The same color scale from −π to π is used as in Fig. 1. (b) The 2D amplitude cross-section corresponding to the central vertical slice. The dislocation core signature can somewhat be seen as a low amplitude region. However, there are two dislocation cores and potentially three low amplitude regions so the identification is ambiguous. (c) The 2D phase cross-section. The two dislocation cores manifest themselves as phase vortices. (d) The 3D dislocation line (yellow) along with the particle shape (semi-transparent red) using the derivative-based method. (e) A close-up view of the dislocation line identified by using the derivative-based method. (f) A close-up view of the line identified using the min-max algorithm. The min-max algorithm is more accurate in identifying the line over all space and more computationally efficient. tion step [21]. The particle shape is represented as an isosurface of the amplitude, which is proportional to the Bragg electron density, while the color projected onto the shape shows the atomic displacement field. The large displacement field on the particle surface is due to the termination of two dislocation lines. A cross section of the amplitude (Fig. 3b) and phase (Fig. 3c) at the particle center show the signature of two dislocation lines. We used the derivative-based method to identify the dislocation line and show the line as a pair of yellow points superimposed onto the particle shape, which is shown as the semi-transparent red isosurface (Fig. 3d). To better visualize the line, we show just the line in Fig. 3e. The line clearly has a horseshoe-like shape, which is indicative of a dislocation transition between edge and screw character [27]. The dislocation line shows some artifacts of the derivativebased method, namely, that some parts of the line appear to be larger or smaller than the voxels. Figure 3f shows the line identified with our new algorithm. The line is more clearly identified, and subtle changes in direction are visible. The new algorithm is also 150x faster than the derivative-based method for this 128x128x64 array size (430s versus 3s). The scaling with larger array size is also much better. When tested on a 256x256x96 array with 83,890 nonzero elements, the difference is 1,065s versus 4.7s, a speedup of 260x. This speedup makes incorporating dislocation identification into existing phase retrieval algorithms feasible. It could be used to perform the transformation to the dislocation basis, which tends to be sparse [24]. Sparsity using the dislocation basis could be used to circumvent traditional constraints in BCDI, ptychography, and other methods that rely on phase retrieval. One may even be able to use subpixel resolution of angles to further refine dislocation core regions by using more sophisticated min-max algorithms such as that in [? ].
We developed a new algorithm that is based on minimum differentials in a local neighborhood for identifying dislocation cores in BCDI images. Using experimentally determined images, we demonstrated that this algorithm is both more accurate and computationally efficient. Using the new algorithm, we identified additional geometric features of the dislocation line. The computational speedup opens the possibility of incorporating the dislocation line into the phase retrieval algorithm. For example, the image is sparse in the dislocation basis, and this feature could be exploited to develop new and improved algorithms for phase retrieval and new experimental techniques in which traditional constraints are relaxed. We expect this algorithm will find immediate use in identifying dislocations in reconstructed experimental images, and we provide the algorithm in the supplemental material.