Multislice simulations

Multislice simulations#

In the previous walkthrough on potentials, we described how to integrate the potential into a series of thin slices \(V_n\) along the \(z\)-axis (see (1)). Now we will describe how the potential slices are used in the multislice algorithm.

Given a wave function of fast electrons at the entrance of the \(n\)’th slice, \(\psi_n(\vec{r})\), and a weak potential slice \(V_n\), the wave function at exit plane of that slice may be written as

(2)#\[ \psi_{n + 1}(\vec{r}) = p(\vec{r}) * \left[t_n(\vec{r}) \psi_n(\vec{r}) \right] \]

where

\[ p(\vec{r}) = \frac{1}{i \lambda \Delta z}\exp\left[\frac{i\pi}{\lambda \Delta z} r^2 \right] \]

is the Fresnel free-space operator for propagation by a distance \(\Delta z\) along the \(z\)-axis, \(*\) is the convolution operator and

\[ t_n(\vec{r}) = \exp\left[i\sigma V_n(\vec{r})\right] \]

is the transmission function that applies a phase shift proportional to the magnitude of the potential slice \(V_n\), where the proportionality constant, \(\sigma\), is called the interaction constant. The derivation of these equations may be found in textbooks or in our appendix.

Using the Fourier convolution theorem, we can write Eq. (2) as

\[ \psi_{n+1}(\vec{r}) = \mathcal{F}^{-1} \left\{P(\vec{k}) \ \mathcal{F}\left[t(\vec{r}) \psi_n(\vec{r})\right] \right\} := \mathcal{M}_n \psi_n(\vec{r}) \quad , \]

where

\[ P(\vec{k}) = \exp(-i \pi \lambda k^2 \Delta z) \]

is the Fresnel free space propagator in reciprocal space and \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) is the Fourier transform and its inverse. For brevity, we have defined a multislice operator, \(\mathcal{M}_n\), acting on a wave function to step it forward through the \(n\)’th potential slice.

Thus, given an initial wave function \(\psi_0\), we can obtain the exit wave function for a potential with slice indices \(n=1,\ldots N\) by sequentially applying these operators

(3)#\[ \psi_{exit}(\vec{r}) := \psi_{N}(\vec{r}) = \mathcal{M}_N \mathcal{M}_{N-1} \ldots \mathcal{M}_{1} \psi_0(\vec{r}) \quad . \]

Warning

Phonon scattering is not included in the following simplified examples, although their inclusions are important. Hence the following results may not correspond to experiment, and more realistic calculations with frozen phonons will be discussed later.

Multislice simulations with plane waves#

Below we create a Potential representing gold with a lattice constant of \(4.08 \ \mathrm{Å}\) in the \(\left<100\right>\) zone axis and use a PlaneWave as the initial wave function.

unit_cell = ase.build.bulk("Au", cubic=True)

atoms = unit_cell * (1, 1, 30)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
abtem.show_atoms(atoms, ax=ax1, title="Beam view")
abtem.show_atoms(atoms, ax=ax2, plane="xz", title="Side view", linewidth=0);

../../_images/b87be55c47e6e7a7dc06c1f248ddfed805cb622e19b4a86afd661c4bdef5bbd1.png

After repeating the structure along \(z\), the simulation cell has a size of \((4.08\times4.08\times122.4) \ \mathrm{Å}\).

atoms.cell

Cell([4.08, 4.08, 122.4])

We set the slice thickness to half the unit cell height \(4.08 / 2 = 2.04 \ \mathrm{Å}\) and the \(xy\)-sampling of the potential to \(0.04 \ \mathrm{Å}\).

Tip

Slice thickness It is beneficial, but not always possible, to make the slice thickness the same as the distance between the atomic planes of the crystal. Doing so typically allows you to use thicker slices.

In all cases, the multislice algorithm is only accurate in the limit of small slice thicknesses, and you should always make sure that your calculation is sufficiently converged with respect to the slice thickness.

potential = abtem.Potential(atoms, slice_thickness=4.08 / 2, sampling=0.04)

plane_wave = abtem.PlaneWave(energy=200e3)

We set up our calculation by calling the multislice method, producing an exit wave function.

exit_wave = plane_wave.multislice(potential)
exit_wave.shape

(102, 102)

As discussed previously, nothing is yet calculated. To execute the simulation we need to call compute.

exit_wave.compute();

Note

Notice that we did not provide the sampling or extent of the wave function above, and so the wave function automatically adopted the grid of the potential. A GridUndefinedError will be thrown if the grid is not defined for both the wave function and potential.

We can then visualize the intensity of the resulting exit wave function.

exit_wave.intensity().show(common_color_scale=True, cbar=True);

../../_images/e829536cf0bc0e4fe088b7492c481cc0dd2d7b16ab09d9f44462fd31356c019e.png

In realistic HRTEM experiments, the wave functions have to be magnified by an objective lens which introduces aberrations and effectively cuts off large scattering angles. How this is handled will be discussed in our introduction to the contrast transfer function.

Here, we apply a defocus of \(-50 \ \mathrm{Å}\) and an objective aperture of \(20 \ \mathrm{mrad}\).

exit_wave.apply_ctf(defocus=-50, semiangle_cutoff=20).intensity().show(cbar=True);

../../_images/d6e150d5846e8fdecf0928cf38fb5a1c2cd4de918b6a8356b0dda7dcb13124e3.png

Electron diffraction patterns#

Instead of an image, we can instead simulate a selected area diffraction (SAD) experiment by using the DiffractionPatterns method. We use block_direct=True to block the direct beam: it typically has a much higher intensity than the scattered beams, and thus it is typically not possible to show it on the same scale.

diffraction_patterns = exit_wave.diffraction_patterns()

diffraction_patterns.block_direct().show(units="mrad", cbar=True)

<abtem.visualize.visualizations.Visualization at 0x29c44165a90>

../../_images/2f8942595bea171cb7202ea240a22a07dd259b0c1175c7e5296fc151d68d5a8a.png

You may wonder; why do the diffraction spots look like squares? This is because the incoming wave function is a periodic and infinite plane wave, hence the intersection with the Ewald sphere is pointlike. However, since we are discretizing the wave function on a square grid (i.e. pixels), the spots can only be as small as single pixels. In real SAD experiments, the spot size is broadened due to the finite extent of the crystal as well instrumental effects.

We can use the index_diffraction_spots method to create a represention of SAD patterns as a mapping of Miller indices to the intensity of the corresponding reflections. The conventional unit cell have to be provided in order to index the pattern, we can provide this as the unit cell of the gold crystal we created earlier, we cannot use the the repeated cell.

indexed_spots = diffraction_patterns.crop(120).index_diffraction_spots(cell=unit_cell)

indexed_spots

<abtem.measurements.IndexedDiffractionPatterns at 0x29c4ce17fd0>

The IndexedDiffractionPatterns facilitates the creation of a visualization that corresponds closer with typical textbook illustrations, where the area of the disks are proportional to the intensity. The disks may be scaled using scale and to limit cluttering miller index annotations maybe excluded by providing a threshold in the annotation_kwargs argument.

visualization = (
    indexed_spots.block_direct()
    .crop(120)
    .show(
        scale=0.02,
        cbar=True,
        annotation_kwargs={"threshold": 0.002, "fontsize": 7},
        figsize=(8, 8),
    )
)

../../_images/21c94cc7dd21aef7e27ec0f86a29082132122208c3eae492a9dc333c453cf82a.png

We see that the \(\{100\}\) reflections are extinguished, as is expected from the selection rules of an F-centered crystal. We can also observe that the \(\left<220\right>\) spots have significantly higher intensity than the \(\left<200\right>\) spots; this is due to dynamical scattering — which is accounted for by the multislice algorithm.

It is possible to obtain a pandas dataframe of the intensity values of the indexed diffraction spots. We only include spots with an intensity of at least \(0.005\) (as a fraction of the incoming beam) using the remove_low_intensity method.

df = indexed_spots.remove_low_intensity(0.005).to_dataframe()
df

	0 0 0	0 4 0	0 -4 0	2 2 0	2 -2 0	4 0 0	-4 0 0	-2 2 0	-2 -2 0
0	0.612	0.005	0.005	0.007	0.007	0.005	0.005	0.007	0.007

Writing exit wave functions#

You can write exit wave function directly to disk, which will also trigger the computation to run.

exit_wave = plane_wave.multislice(potential)

exit_wave.to_zarr("./data/exit_waves.zarr", overwrite=True);

We can read the wave function back in as shown below, and see that they are identical to the calculated exit wave.

imported_wave = abtem.from_zarr("./data/exit_waves.zarr").compute()

assert imported_wave.compute() == exit_wave.compute()

Thickness series#

abTEM easily allows us to obtain the wave function at intermediate steps of the multislice algorithm, thus allowing us to see how the wave function evolves as it passes through the potential. To create such a thickness series we set the exit_planes keyword of the potential. exit_planes may be given as a tuple of slice indices at which to return the wave function, or simply as a single integer to indicate the step between those slice indices.

Below we create a Potential as above, but we set exit_planes=4. When running the multislice simulation we obtain an ensemble of wave functions \(\psi_n(\vec{r})\) at slice indices \(n=0,2,4,\ldots\), ie. after every unit cell.

potential_series = abtem.Potential(
    atoms, slice_thickness=4.08 / 2, sampling=0.05, exit_planes=2
)

exit_wave_series = plane_wave.multislice(potential_series).compute()

The thickness series consists of \(31\) exit waves.

exit_wave_series.shape

(31, 82, 82)

We do not want to show all \(32\) exit waves, hence, we can select every fifth exit wave using standard NumPy indexing. We show the wave function intensity as an exploded plot, enforcing a joint intensity color scale with a given maximum value.

exit_wave_series[::5].show(
    explode=True,
    figsize=(13, 5),
    common_color_scale=True,
    cbar=True,
    vmin=0,
    vmax=6,
);

../../_images/5d84b598699449783e678e4ee5d250d24092a9cd644a118f7aea8d1c03357a19.png

We can naturally get the indexed diffraction patterns for the series.

diffraction_patterns_series = exit_wave_series.diffraction_patterns(max_angle=120)

indexed_spots = diffraction_patterns_series.index_diffraction_spots(
    cell=4.08,
)

For a relatively heavy atom such as gold, the diffraction spots in the first order Laue zone are excited already at \(20 \ \mathrm{Å}\).

indexed_spots[::5].block_direct().show(
    explode=True,
    common_color_scale=True,
    figsize=(16, 5),
    cbar=True,
    scale=0.05,
    annotation_kwargs={"threshold": 1.0, "fontsize": 7},
);

../../_images/4d2c7ced2381918747b8137255470f08dd45bd2be46ea7551903760af38857fb.png

The diffraction intensities can be obtained as a pandas dataframe indexed by the thickness. We set a threshold to only include spots which at some thickness have an intensity above \(0.005\).

df = indexed_spots.remove_low_intensity(0.005).to_dataframe()

We select every sixth diffraction spot and show the dataframe.

df.iloc[::6]

hkl	0 0 0	0 2 0	0 4 0	0 -4 0	0 -2 0	2 0 0	2 2 0	...	-4 2 0	-4 -2 0	-2 0 0	-2 2 0	-2 4 0	-2 -4 0	-2 -2 0
z [Å]
0.000	1.000	0.000	0.000	0.000	0.000	0.000	0.000	...	0.000	0.000	0.000	0.000	0.000	0.000	0.000
24.480	0.698	0.015	0.007	0.007	0.015	0.015	0.011	...	0.006	0.006	0.015	0.011	0.006	0.006	0.011
48.960	0.842	0.007	0.001	0.001	0.007	0.007	0.001	...	0.001	0.001	0.007	0.001	0.001	0.001	0.001
73.440	0.679	0.005	0.006	0.006	0.005	0.005	0.007	...	0.005	0.005	0.005	0.007	0.005	0.005	0.007
97.920	0.683	0.009	0.004	0.004	0.009	0.009	0.005	...	0.003	0.003	0.009	0.005	0.003	0.003	0.005
122.400	0.629	0.004	0.004	0.004	0.004	0.004	0.006	...	0.004	0.004	0.004	0.006	0.004	0.004	0.006

6 rows × 21 columns

We can now leverage all the pandas features, for example, we can select three diffraction spots and plot their intensities together.

ax = df[["0 4 0", "0 2 0", "2 2 0"]].plot()
ax.set_ylabel("intensity [arb. unit]");

../../_images/85cf98ce035ee6ff833a1b41bf810cccd7680eeeddfe4d57b3c3d5f4ce915726.png

Multislice simulation with a probe#

Simulations of a focused electron probe differ only in minor details. We create our initial wave function as a convergent beam with an energy of \(200 \ \mathrm{keV}\) (note the conversion from \(\mathrm{keV}\) to the default units of \(\mathrm{eV}\)) and a convergence semiangle of \(9 \ \mathrm{mrad}\).

probe = abtem.Probe(energy=100e3, semiangle_cutoff=9.0)

There are two additional considerations to take into account when performing multislice simulations with probes: self-interaction errors and probe placement.

Self-interaction#

Self-interaction errors happen when the probe overlaps with its periodic images. We see this happening below: the clearest visual clue is that the probe does not have a circular symmetry, even though we defined no non-symmetric aberrations. (Note that we need to define a grid for the probe to display it, which can be done manually above, or by matching it to the grid of the potential as below).

probe.grid.match(potential)
waves = probe.build().compute()

waves.show(title="probe intensity (small window)");

../../_images/318a3069f38d2bade22c93ee1e76f83415359422de3c126e23fa12270bbc8083.png

We get rid of self-interaction by increasing the extent of the potential until the probe intensity is sufficiently small at the boundary. Note that this should be true at any slice of the potential when running the multislice algorithm.

The potential may be extended by tiling the atoms (see walkthrough on atomic models), but for crystals we can also use the CrystalPotential we introduced earlier. Here we repeat the potential used above by \(4\) times in and \(x\) and \(y\).

tiled_potential = abtem.CrystalPotential(potential, (4, 4, 1))

Now, when the probe is matched to the potential, there is no significant self-interaction. Note that the probe window should be large enough to accomodate the probe through all the slices in a multislice simulation.

probe.grid.match(tiled_potential)
waves = probe.build().compute()

waves.show(title="probe intensity (adequate window)");

../../_images/3a6f13779e348a8f96edae368d082ea88f32312a19d589a132743d108b63f338.png

Probe position#

Above the probe was placed (by default) at the center of the potential. Later we will learn how to scan the probe for STEM simulations, but for our purposes here, we manually select a few positions.

We create an array of \(5\) probe positions across the center of the potential (found by taking half of its full extent in both directions)

center = (tiled_potential.extent[0] / 2, tiled_potential.extent[1] / 2)

displacements = np.array([[0, 0], [0.501, 0], [1.02, 0]])

positions = center + displacements

positions

array([[8.16 , 8.16 ],
       [8.661, 8.16 ],
       [9.18 , 8.16 ]])

Running the multislice algorithm with a scan given by these probe positions, we obtain an ensemble of \(3\) wave functions, one for each of the positions.

exit_wave_probe_positioned = probe.multislice(scan=positions, potential=tiled_potential)

exit_wave_probe_positioned.compute();

exit_wave_probe_positioned.show(
    explode=True,
    figsize=(10, 5),
);

../../_images/9b6dafdbcdc2de127bb986f1f0857c9908dbe7d72e90164445007ad8aacf346d.png

This is the first time we encounter a slightly longer calculation time: in fact, we are running five independent full multislice passes, one for each probe position, which takes more time.

Convergent beam diffraction#

A convergent beam (electron) diffraction (CBED) pattern may be calculated in the same manner as the SAD pattern. Here we select the third probe position by indexing the ensemble created above and calculate its diffraction pattern.

exit_wave_probe_positioned[2].diffraction_patterns(max_angle=30).show(cbar=True);

../../_images/e5f55e863a5605f9160f09940a5cafc38a0e7632b6cda50028e1a20f8062c0b5.png

The diffraction pattern is rather pixelated due to our limited sampling. Improving reciprocal-space sampling is equivalent to tiling in real space. Below the potential is repeated by \(8\) times in \(x\) and \(y\), thus improving the reciprocal space sampling by a factor of \(2\). We also set exit_planes=36 to obtain a series of \(5\) exit waves.

potential_extra_tiled = abtem.CrystalPotential(potential, (8, 8, 1), exit_planes=36)

exit_wave_extra_tiled = probe.multislice(potential_extra_tiled).compute()

The potential is now significantly larger, and thus the simulations again take a little time. We show the resulting series of diffraction patterns below.

exit_wave_extra_tiled.diffraction_patterns(max_angle=30).show(
    explode=True,
    common_color_scale=True,
    cbar=True,
    figsize=(12, 5),
);

../../_images/96698dcab63aa03058b8e8050a5708cd5273ec7b4f6f3cf509eeca7c80b25e75.png

Small-angle beam tilt#

Small amounts of beam tilt (or equivalently sample tilt) may be included by a small modification to the propagator function

\[ P(\vec{k}) = \exp\left[-i \pi \lambda k^2 \Delta z + 2 \pi i \Delta z (k_x \tan\theta_x + k_y \tan\theta_y)\right] \quad , \]

where \(\theta_x\), \(\theta_y\) is the beam tilt in the \(x\), \(y\) directions. This is equivalent to shifting the wave function between slices and is only valid for very small tilts of no more than about \(100 \ \mathrm{mrad}\).