Single-pass STEM-EMCD on a zone axis using a patterned aperture: progress in experimental and data treatment methods

Measuring magnetic moments in ferromagnetic materials at atomic resolution is theoretically possible using the electron magnetic circular dichroism (EMCD) technique in a (scanning) transmission electron microscope ((S)TEM). However, experimental and data processing hurdles currently hamper the realization of this goal. Experimentally, the sample must be tilted to a zone-axis orientation, yielding a complex distribution of magnetic scattering intensity, and the same sample region must be scanned multiple times with sub-atomic spatial registration necessary at each pass. Furthermore, the weak nature of the EMCD signal requires advanced data processing techniques to reliably detect and quantify the result. In this manuscript, we detail our experimental and data processing progress towards achieving single-pass zone-axis EMCD using a patterned aperture. First, we provide a comprehensive data acquisition and analysis strategy for this and other EMCD experiments that should scale down to atomic resolution experiments. Second, we demonstrate that, at low spatial resolution, promising EMCD candidate signals can be extracted, and that these are sensitive to both crystallographic orientation and momentum transfer.


Introduction to this script and the supplementary information
This document guides the reader through all of the processing steps necessary to reach the conclusions drawn in the paper it supplements. It was written using the Matlab Live Script utility in version R2018b. The purpose is to introduce the code used in the paper and reproduce the entire processing workflow necessary to generate the figures. The functions used in this script are provided alongside the data and are available through the Zenodo website at the following URL or using the digital object identifier listed below: https://zenodo.org/record/3361582 10.5281/zenodo.3361582

Description of the files
There are four data files that accompany this script. Two of the files are the raw data as captured in the TEM and two of them are timestamp metadata. The two raw data files are saved in Gatan's proprietary .dm4 format, while the two timestamp files (which allow for synchronizing the data frame to the probe position) are saved in Gatan's proprietary .gtg persistent tag structure format. The timestamp files are used to arrange the data into a 4D datacube, which subsequently allows for the STEM-DP and 4D-EELS datasets to be correlated through their common spatial dimensions.
• EELS4D_Movie.dm4: Image stack of all of the 2D EELS frames that were captured from the Ultrascan camera running in "view" mode using the hook-up script written by Linus Schönström and a modification of the frame recorder written by Vincent Hou (see the manuscript for details). • DP_Movie.dm4: Image stack of all of the CBED frames that were captured from the Orius camera running in "view" mode using the hook-up script written by Linus Schönström and a modification of the frame recorder written by Vincent Hou (see the manuscript for details). • STEMDP_TimeStamps.gtg: Collection of time stamps for the STEM-DP experiment. It goes along with the file DP_Movie.dm4. • EELS4D_TimeStamps.gtg: Collection of time stamps for the 4D EELS experiment. It goes along with the file EELS4D_Movie.dm4.

Dependencies
The first public release of this script and the accompanying code is v1.0. This script assumes that this release is being used. Later releases may break the code.
The use of this script relies on a number of third party software packages. First and foremost, both a recent version of Matlab as well as the optimization toolbox is required. These are commercial packages that must be purchased from The Mathworks, Inc. This script was written in version R2018b.
The raw data were recorded and are thus provided in a proprietary format developed by Gatan, Inc. Andreas Korinek's dmread.m script, which, at the time of writing, can be downloaded from the Matlab file exchange, is used to read and import this format into Matlab. It is included in the package for this work in unmodified form aside from a small bug fix that is commented.
The divergent colormaps come from a Matlab implentation BrewerMap, which is based off of the colormaps developed by Cynthia Brewer (http://colorbrewer.org). A copy of this library is included in the repository with this code.

Explore the STEM-DP datacube
The code provided with this script does not include a Matlab function to explore 4D Datasets effectively. We recommend that one exports the datacube to Hyperspy and uses the PyXem package. For now, the datacube can be programatically explored by generating virtual HAADF images, to ensure that the patterns were properly imported. Orientation masks using virtual dark field samples will then also be created. Proper placement of the virtual apertures was empirically determined by interactively exploring the data in PyXem. Finally, a "grain" mask that represents the entire iron grain under investigation will also be produced. These masks reproduce the Region Of Interest (ROI) figure presented in the paper. In this section, we use the STEM-DP data to generate orientation masks. If one browses the data (for example using PyXem), one sees that some regions of the iron grain are closer to the [1 0 0] zone axis orientation than others.
To estimate this, we look at a Kikuchi line that varies quite a bit as the orientation changes. This will give us an estimate of a zone-axis orientation vs. an orientation that is tilted further away.

CBED patterns for different orientations
The average CBED patterns for the different orientation masks are generated here, reproducing part of the ROI figure from the paper.

4D EELS pretreatment
The EELS pretreatment is performed before the stack of 2D EELS frames is assembled into a 4D datacube.

X-Ray spike identification and removal
In this section, x-ray spikes are identified and removed. The x-ray spikes are determined by looking for datapoints that deviate from a moving average by more than 5 sigma. The area surrounding these datapoints is then dilated and they are replaced by a linear interpolation from the nearest datapoints.

Rough energy alignment (correcting for the drift-tube shift)
This first step corrects for the sawtooth drift tube offset that was applied during the data acquisition. Since the offset was acquired in a separate thread, the exact values are unknown. Cross-correlation with a reference spectrum is required to retrieve them.
The EELS data are first vertically summed (integration along q y ). The assumption here is that this is what happens under standard EELS acquisition. Also, each 2D EELS frame was acquired with the same drift tube offset, and this integration helps with the signal to noise ratio. Once this is done, a good looking spectrum that can be used as a reference is recorded to variable ref.

Raw energy drift tube shift
The sawtooth function is clearly visible. In regions where the Fe L 2,3 edges were too weak (such as over the background), the cross-correlation failed, giving a very large erroneous value.

Correct for outliers
Erroneous drift tube offset datapoints are easily corrected for by selecting extreme values. Extrapolation of the sawtooth function towards the end of the ESI acquisition is achieved by modelling the sawtooth function (although this is admittedly not necessary, as these datapoints don't contain iron and, thus, will not be used in the EMCD analysis).

Apply the rough ISOmap and visualize it
The offset correction is now applied to the full 4D EELS datacube by iterating through each q y vector. colormap inferno (1024) ; 8 axis normal 9 catch 10 fprintf (1 , Ensure that the GUI Layout Toolbox is installed ) ; 11 end 19/36

Fine energy drift correction
As discussed in the manuscript, the above energy drift correction only corrects (approximately) for the drift tube offset, providing gain averaging. However, close inspection of the EELS data reveals that there is some energy drift dependence on q y as well. This is corrected for here by using a peak shifting methodology based on regressing the spectral mean and its first derivation against the raw data. The input for this step is the roughly-aligned dataset.
1 mask . Ref = 1500:1615; % channels . Only select Fe L3 2 mask . meanEELSSpectra = mean (B , 1) ; % Restrict to iron grain 3 mask . meanEELSSpectra = squeeze ( mask . meanEELSSpectra ) ; 4 mask . BG = mask . meanEELSSpectra < 28; % Based on 4 D summation 5 6 % taylorShift is designed for typical 3 D hyperspectral EELS datacubes , which 7 % we currently have , since we have not interpolated to an (x , y ) datacube . The energy offset of Fe L 3 as a function of both scan index (pixel position) and q y can now be visualized. A noticable energy shift between the middle region (low values of q y ) and the outer regions is observed. Thus, this correction step is crucial when comparing chiral spectra that are derived from q y regions of opposite sign. We now return to the original (raw) dataset and apply the full energy drift correction. This includes both the drift tube offset as well as the q y -dependent correction. The result is saved in variable C, which is ready for interpolation into a 4D datacube (allowing it to be correlated to the STEM-DP datacube through the x and y spatial dimensions).

Interpolation to 4D datacube
For the final preprocessing stage, the stack of 2D EELS frames (dimensions ∆E × q y × N px ) is assembled via interpolation into a 4D datacube (dimensions ∆E × q y × x × y). This datacube is moreover cropped to the Fe L 2,3 edges and extreme values of q y falling outside of the passthrough region of the aperture are removed to cut down on the data size. To write the 2D EELS frames into a 4D datacube, the following for loop iterates through each row. The timestamps (in variable ts) assign the frame index to the appropriate row. While most probe positions can be assigned a single 2D EELS frame, some will be missing and some will have excess frames. This is mitigated through the use of a nearest neighbor interpolation for each row, which is what the code below achieves. The 4D EELS datacube is now reshaped into the four dimensions and permuted for convenience.

DPEELS = reshape ( DPEELS , length ( mskE ) , length ( goodQ ) , nCols , nRows ) ; 2 DPEELS = permute ( DPEELS , [2 4 3]) ;
A comparison with the virtual HAADF image from the CBED datacube can be used to ensure that the spatial registration is reasonable. This is achieved in the 4D EELS datacube here by integrating along both the ∆E and q y dimensions, yielding an x, y map that represents the total differential inelastic scattering cross sections.

Candidate EMCD signal extraction:
Necessary function handles q is a simple function handle that selects the opposing q-values for the summed 2D EELS spectra. This assumes that the center of the summation lies at q = 0. This was determined by careful placement of the patterned aperture during experimental acquisition, and the centering we belive to be most accurate is applied here. If you want to see the effect of changing this center value, then re-run the section above with different values for goodQ. sp is a function handle that sharpens (profile matches) one of the two chiral spectra. The approach is a standard removal of the second derivative. Validation of this approach can be performed by integrating the second derivative. If there is a zero-crossing between the L 3 and L 2 peaks, then the area under the peaks will not change and, hence, this will not directly bias the calculation of m L /m S (although this is not explicitly studied in this manuscript due to the low SNR). This profile matching step is necessary since there is a very slight energy blur that arises from spectrometer errors. The matching routine is parameterized by multiplying the second derivative by a scalar value (k) and this is what is allowed to vary in the optimization code for EMCD signal extraction. This function handle serves to set an initial starting parameter, which is occasionally needed to avoid reaching a local minimum. This is something that you can play around with to see the result to get a feel for what it does to the data. The starting paremeter used below was empirically determined, but typically was in the range of -1 to -2.
1 sp = @( Sp , k ) Sp -k .* gradient ( gradient ( smooth ( Sp ) ) ) ; To ensure that the sp function does not change the area under the curve: The following sections run the optimization function optimMLMS.m for EMCD signal extraction. This function is outlined in the manuscript but more details are available within the function code itself. In short, it takes two chiral 27/36 spectra as input (output from q above), determines starting parameters (prepEMCDoptim.m), writes these to variable fi (only accessible in the function workspace), passes this to Matlab's fmincon function along with the objective function objMLMS.m. The parameters in fi are listed in the table in the manuscript. Some of these parameters are used to generate a fit function to the EMCD signal (which is the difference between background-subtracted, profile-matched, properly normalized chiral spectra). This fit function is subtracted from the raw EMCD signal leaving residual errors. The sum of the square of these errors is minimized with fmincon. The resulting optimal parameters are returned as variable fo and written to the active Matlab workspace. A struct containing all of these parameters as well as the processed data is returned as fStruct in the Matlab workspace and can then be inspected.
As input parameters, a number of different models can be used. The standard inverse power law and the linear combination of power laws (Cornell Spectrum Imager) are options for the pre-edge background model. The post-edge linearization can either use the standard integration approach or the linear regression approach discussed in the manuscript. The residual errors that are minimized can be weighted (not performed here) and decoupled as needed. In this case, the pre-edge regions for both chiral spectra are "decoupled" and optimized in serial with the residuals from the EMCD signal, which helps keep the backgrounds more reasonable. Finally, the EMCD signal model can be chosen. The standard double-gaussian approach can be chosen, but this fails to capture the extended tails of the EMCD signal. The pseudo-Voigt function approach is called with the 'GLmix' input argument and is what is performed here. An additional skewed Gaussian model is available for playing with. Details are provided within the function itself.
The data are visualized by calling the function plotEMCDfitCum.m, which accepts the struct output from optimMLMS.m as an argument and returns a formatted plot that is used in the manuscript. Unfortunately, at the time of writing, Matlab's Live Script environment does not allow for interactive or updatable figures. Consequently, it is not possible to watch the progression of this optimization routine during the iterations. This can be accomplished by copying the code below that calls optimMLMS (along with all its additional arguments) but adding a fourth output argument (for example, flm) and executing this in the Matlab command line directly. This will both allow the graph to update at every iteration, but will also save each frame (iteration) as an entry in a large struct variable in the workspace (flm). Then, the included makeMovie.m function can be called as follows: makeMovie ( flm , ' FileName . mp4 ' , 1 0 ) ; This will create and save an .mp4 file into the current path with the name 'FileName.mp4' and, in this case, a speed of 10 frames per second.

Using mask.orient01 (close to zone axis) and restricted qy
The first EMCD signal is extracted from the sample region where the iron grain was oriented close to a zone axis. This is the ROI called "Orient01" in the paper. The values of q y that are used for this signal extraction are limited to the "outer" q y region. The reason for this is clear from below when the q y maps are generated. This reproduces figure 5 in the manuscript.

Using mask.orient02 (tilted away from zone axis) and restricted qy
The second EMCD signal is extracted from the sample region where the iron grain was mistilted. This is the ROI called "Orient02" in the paper. The values of q y are the same as before, in order to compare the effect of crystallographic orientation on the EMCD signal quality  This reproduces figure 6 in the manuscript.

Using mask.grain and all qy values
The third EMCD signal is a summation over all crystallographic orientations and q y values.  This reproduces figure 7 in the manuscript.

Using mask.grain (integral over all orientations) and restricted qy
The effect of q y can be assessed by integrating over a restricted range of q y values. This EMCD signal is extracted using the "outer" 4 rows, as with the orientation masks above.

EMCD q-maps
The q − E maps for the EMCD signal can be created by running optimMLMS.m on each individual q y value. This is achieved by iterating over all 12 available conjugate q y vectors in the dataset. All three of the ROI masks are tested with this approach. These final three maps are composited together to make figure 9 in the manscript.