Core-loss

Core-loss#

This tutorial demonstrates how to calculate an “energy-filtered” STEM image using abTEM. We show how to calculate projected overlap integrals following to C. Dwyer (see https://doi.org/10.1016/j.ultramic.2005.03.005 and use them in a dynamical scattering simulation. We use the Density Functional Theory code GPAW for calculating the wave functions for the overlap integrals. Please see our citation guide if you use this code in a publication.

We will first roughly reproduce the results of C. Dwyer, which the implementation is largely based on. In the second half of the tutorial, we present a the typical core-loss energy filtered simulation using abTEM.

Currently Under Construction

Check back for more…

Core-loss transitions#

The first step is to select the transitions. Here we target the K-edge of Silicon i.e. the 1s subshell; quantum numbers $(n, \ell) = (1, 0)$. In a hypothetical experiment, this would be roughly equivalent to only collecting electrons with an energy loss of $10 \ \mathrm{eV}$ above the ionization threshold of $1840 \ \mathrm{eV}$.

Z = 14  # atomic number
n = 1  # principal quantum number
l = 0  # azimuthal quantum number
xc = "PBE"  # exchange-correlation functional

Si_transitions = SubshellTransitions(Z=Z, n=n, l=l, xc="PBE", order=2, epsilon=10)

print("bound electron configuration:", Si_transitions.bound_configuration)
print("ionic electron configuration:", Si_transitions.excited_configuration)

bound electron configuration: 1s2 2s2 2p6 3s2 3p2
ionic electron configuration: 1s1 2s2 2p6 3s2 3p2

There is an infinite number of possible transitions from the bound state to the continuum, however, only transitions with small changes in orbital angular momentum quantum number, $\Delta \ell$, will to contribute meaningfully. Choosing a maximum transition order, $\Delta \ell_{max}$, we include only transitions with continuum orbital angular momentum quantum number: $$ \ell' = \ell, \ell \pm 1, \ell \pm 2, \ldots , \ell \pm \Delta \ell_{max} \quad \ell' \geq 0 $$ This is controlled by setting i.e. order = 0 gives just the monopole term, order = 1 gives also includes all dipole terms and order = 2 includes up to all quadropole terms etc.

Given an orbital angular momentum quantum number $\ell$ (or $\ell'$) the possible magnetic quantum numbers are:

\[ m_{\ell} = -\ell, \ell + 1, \ldots, \ell \quad . \]

Below we print the quantum numbers for the bound and continuum states up quadropole transitions:

for bound_state, continuum_state in Si_transitions.get_transition_quantum_numbers():
    print(f"(n, l, ml) → (l', ml') = {bound_state} → {continuum_state[1:]}")

(n, l, ml) → (l', ml') = (1, 0, 0) → (0, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, -1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (2, -2)
(n, l, ml) → (l', ml') = (1, 0, 0) → (2, -1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (2, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (2, 1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (2, 2)

The wave functions for the bound and continuum states, the atomic wave functions can be written as products of radial and angular parts, i.e $$ \phi_{nm\ell} = \frac{P_{n\ell}(r)}{r} Y_{\ell}^{m}(\hat{\mathbf{r}}) \quad \mathrm{and} \quad \phi_{nm\ell} = \frac{P_{\epsilon_f\ell'}(r)}{r} Y_{\ell'}^{m'}(\hat{\mathbf{r}}) \quad , $$

where $Y_{\ell}^{m}$ represent spherical harmonics. In abTEM the radial wave functions, $P_{n\ell}$ and $P_{\epsilon_f\ell}$, are calculated using a all-electron DFT with using the scalar-relatistically corrected Schrödinger equation. We can get the bound radial wave function, $P_{10}$, using get_bound_wave_function and the continuum wave function, $P_{\epsilon_f 0}$, $P_{\epsilon_f 1}$ and $P_{\epsilon_f 2}$, using get_excited_wave_functions.

bound_state = Si_transitions.get_bound_wave_function()
excited_states = Si_transitions.get_excited_wave_functions()

bound_state, excited_states

(<abtem.inelastic.core_loss.RadialWavefunction at 0x7f4574f99650>,
 [<abtem.inelastic.core_loss.RadialWavefunction at 0x7f4575163390>,
  <abtem.inelastic.core_loss.RadialWavefunction at 0x7f4476c93990>,
  <abtem.inelastic.core_loss.RadialWavefunction at 0x7f4476ab5810>])

Below we plot the the bound and radial wave functions, this is a reproduction of Fig. 1 in Dwyers paper. Note, that there are minor differences, as Dwyer uses a Herman–Skillman wave functions instead of DFT.

The left and right y-axes correspond to $P_{10}$ and $P_{\epsilon_f \ell '}$. The unit of length in the figure is the Bohr radius $a_0 = 0.529 \ \mathrm{Å}$ and the unit of energy is the Rydberg energy $\mathrm{Ry} = 13.6 \ \mathrm{eV}$.

../../_images/0dc982602046dd80baed8cfea86d1873d9c641395a307342d6b32ed84ca1fb01.png

Transition potentials#

We want to calculate the lateral part of the wave function of the inelastically scattered electrons, $\psi_{f}(\mathbf{r}_{\perp}, z)$, just after undergoing scattering at a depth $z$. We first calculate the elastic wave function $\psi_{i}(\mathbf{r}_{\perp}, z)$ at by propagating the wave function using multislice, then the inelastically scattered wave function may be given as a product with a projected transition potential $$ \psi_{f}(\mathbf{r}_{\perp}, z) = V_{fi}(\mathbf{r}_{\perp}) \psi_{i}(\mathbf{r}_{\perp}, z) \quad . $$

The inelastic transition potential is given by $$ V_{fi}(\mathbf{r}_{\perp}) = \int_{-\Delta z / 2}^{\Delta z / 2} \langle \phi_f | V_e | \phi_i \rangle e^{2 \pi i (k_i − k_f) z} \mathrm{d}z $$ with $$ \langle \phi_f | V_e | \phi_i \rangle = \int \phi^*_f(\mathbf{r}) \frac{e^2}{4\pi\epsilon_0 |\mathbf{r} - \mathbf{r}'|} \phi_i(\mathbf{r}') \mathrm{d}\mathbf{r}' \quad , $$ where the operator, $V_e$, is the Coulomb potential mediating the interaction, $\phi_i$ is the bound state of the transition and $\phi_f$ is one of the included continuum states.

We will from now on only include 1st order transitions, hence, we again create the SubshellTransitions for the Silicon K-edge setting order=1.

Si_transitions = SubshellTransitions(Z=Z, n=n, l=l, xc="PBE", order=1, epsilon=10)

for bound_state, continuum_state in Si_transitions.get_transition_quantum_numbers():
    print(f"(n, l, ml) → (l', ml') = {bound_state} → {continuum_state[1:]}")

(n, l, ml) → (l', ml') = (1, 0, 0) → (0, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, -1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 1)

We create the transition potentials for an incoming wave function with an energy of $100 \ \mathrm{keV}$, defined on a grid with $200 \times 200$ grid points and an extent of $2 \ \mathrm{Å}$.

transition_potential = Si_transitions.get_transition_potentials(
    gpts=200, extent=2, energy=100e3
)

We can show the complex transition potentials using domain coloring. This figure is equivalent in content to Fig. 2 in the Dwyer paper.

visualization = transition_potential.show(
    explode=True,
    figsize=(13, 4),
    cbar=True,
    common_color_scale=True,
    vmax=0.42,
)

../../_images/66cb26dc203fd9f6d5ff8f3794402aa80e1183c4fe7bee79fbe8448a80204011.png

Energy differential cross-sections#

Note

The following is not equivalent to a typical abTEM simulation using a transition potential, see next section.

Next, we shall calculate the probabilty per unit energy (or energy differential cross-section) for core-loss scattering at a depth $z$ from a specific atom up to a given scattering angle $k_{cut}$

(10)#\[ \frac{\mathrm{d}\sigma}{\mathrm{d}E}\left(z\right) = \sigma^2 \int_{0}^{k_{cut}} A(k) \sum_{f} |\psi_f(\mathbf{k}_{\perp}, z)|^2 \ \mathrm{d}\mathbf{k} \quad , \]

where $\sigma$ on the rhs is the interaction constant, $V_{fi}$ is the transition potential for a specified transiton of an atom in a specified position, $\psi(\mathbf{k}, z)$ is the reciprocal space wave function at a depth $z$ and $A(k)$ is the aperture function describing the maximum collection angle of the energy filter

\[\begin{split} A(k) = \begin{cases} 1 & \text{if } k \leq k_{cut} \\ 0 & \text{if } k > k_{cut} \end{cases} . \end{split}\]

This equation neglects all elastic scattering after the inelastic scattering event, sometimes called single-channeling. We describe the more accurate double-channeling, in the next section.

The above equation is essentially equivalent Eq. 16 in the Dwyer paper, the main difference being that we integrate in reciprocal space. We shall roughly reproduce the bold full line shown in Fig. 4 of that paper, this is the result obtained by placing the beam on an atomic column and calculating the equation above for inealstic scattering from only atoms in that column.

We create an atomic structure of silicon in the $\left<100\right>$ zone axis with a thickness $500 \ \mathrm{Å}$ and a converged lateral extent. Then we create a potential with exit planes every ten slices.

atoms = ase.build.bulk("Si", cubic=True) * (3, 3, 92)
atoms.center(axis=2)

potential = abtem.Potential(atoms, gpts=256, exit_planes=1, slice_thickness=5.43)

Note

In order to keep the calculation time lower

We create the initial probe placed on the atomic column at $(0,0)$ with the same parameters as Dwyer and calculate the exit_waves for different thicknesses (i.e. $\psi_i(\mathbf{r}_{\perp}, z)$).

probe = abtem.Probe(
    semiangle_cutoff=30, energy=100e3, Cs=-0.053e-3 * 1e10, defocus=-90, C5=50e-3 * 1e10
)

wave = probe.multislice(potential, scan=(0, 0))

exit_waves = wave.compute()

[########################################] | 100% Completed | 960.79 ms

We can use the transition potentials to calculate the scattered waves, $\psi_f(\mathbf{k}_{\perp}, z)$, for scattering from an atoms placed at $(0,0)$ by setting sites=(0,0).

transition_potentials = Si_transitions.get_transition_potentials()

scattered = transition_potentials.scatter(exit_waves, sites=(0, 0))

The scattered wave functions are returned as a 4D array, where the first axis represents the transitions, the second axis represents the exit planes and the last two dimensions represents the $xy$-plane.

scattered.shape

(4, 93, 256, 256)

Next, we calculate the incoherent sum over the transitions in reciprocal space, $\sum_{f} |\psi_f(\mathbf{k}_{\perp}, z)|^2$, by first getting the diffraction_patterns and then summing the first axis.

diffraction_patterns = scattered.diffraction_patterns().sum((0,))

We show every twentieth “diffraction pattern” below.

diffraction_patterns.shape

(93, 171, 171)

visualization = (
    diffraction_patterns[1::20]
    .crop(90)
    .show(
        explode=True, units="mrad", common_color_scale=True, cbar=True, figsize=(12, 5)
    )
)
visualization.axis("lower left", ticks=True)

../../_images/52e22fdcb0b9b671ba57ea15281fda8f84d23fd096feaf27f4061b44692cd8b9.png

To reproduce Fig. 4, we perform integrals with three different integration ranges from $0$ to $5$, $20$ and $50 \ \mathrm{mrad}$. The minor visible differences between Dwyers results and ours could be due to the slightly different atomic wave functions or because we neglected phonons to keep the demonstration simple.

outer = (5, 20, 50)
measurements = diffraction_patterns.integrate_radial(inner=0, outer=outer)

We can use xarray show the data.

da = measurements.to_data_array()

fig, axes = plt.subplots(1, 3, figsize=(13, 4))

for i, ax in enumerate(axes):
    da[i].plot.line(ax=ax)

../../_images/c2afee4e7b7e459ca7143fbe075b5a72212ef6eed10b89302878651153c7b3b2.png

Core-loss filtered#

We will now run a simulation using the standard interface of abTEM. The main difference compared to above is that the contribution from all inelastic scattering events are collected into one single. When multiple exit planes are specified for the potential all events above each exit plane “collects” all previous scattering events.

We also neglected scattering after the inelastric event (single-channeling), by default this is not the case in abTEM, but you can use it by setting double_channel=False. Double-channeling is always more accurate, but it is drastically slower, and as we shall see may not be necessary in some case for modern EELS detectors with large collection angles.

Single channeling requires the same number of multislice step as standard multislice; however, it is still necessary to calculate the scattering probabilities in Eq. (10) which may be expensive.

Core-loss in atomic resolution imaging#

For the purpose of keeping the computational load of this tutorial realtively low, we will use the single-channelling approach, thus these simulations are only accurate for large maximum collection angles.

We will simulate a zone axis simulation of SrTiO₃, we use slightly

from ase.io import read

atoms = read("../walkthrough/data/srtio3_100.cif")

slice_thickness = atoms.cell[2, 2] / 2

atoms *= (4, 4, 25)

potential = abtem.Potential(
    atoms, gpts=256, slice_thickness=slice_thickness, device="cpu"
)

abtem.show_atoms(atoms, legend=True);

../../_images/ed2a1da6706709d94eda8a0c9ff26fa45a08dda75de9279e329b44defe7df740.png

We target the K-edge of oxygen and the L$_{23}$-edge of strontium. As those would most typically be used in an experiment.

O_transitions = SubshellTransitions(Z=8, n=1, l=0, xc="PBE", order=1, epsilon=10)

Sr_transitions = SubshellTransitions(Z=38, n=2, l=1, xc="PBE", order=1, epsilon=10)

print(f"Oxygen K: N_t = {len(O_transitions)}")
for bound_state, continuum_state in O_transitions.get_transition_quantum_numbers():
    print(f"(n, l, ml) → (l', ml') = {bound_state} → {continuum_state[1:]}")

print(f"Strontium L23: N_t = {len(Sr_transitions)}")
for bound_state, continuum_state in Sr_transitions.get_transition_quantum_numbers():
    print(f"(n, l, ml) → (l', ml') = {bound_state} → {continuum_state[1:]}")

Oxygen K: N_t = 4
(n, l, ml) → (l', ml') = (1, 0, 0) → (0, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, -1)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 0)
(n, l, ml) → (l', ml') = (1, 0, 0) → (1, 1)
Strontium L23: N_t = 27
(n, l, ml) → (l', ml') = (2, 1, -1) → (0, 0)
(n, l, ml) → (l', ml') = (2, 1, -1) → (1, -1)
(n, l, ml) → (l', ml') = (2, 1, -1) → (1, 0)
(n, l, ml) → (l', ml') = (2, 1, -1) → (1, 1)
(n, l, ml) → (l', ml') = (2, 1, -1) → (2, -2)
(n, l, ml) → (l', ml') = (2, 1, -1) → (2, -1)
(n, l, ml) → (l', ml') = (2, 1, -1) → (2, 0)
(n, l, ml) → (l', ml') = (2, 1, -1) → (2, 1)
(n, l, ml) → (l', ml') = (2, 1, -1) → (2, 2)
(n, l, ml) → (l', ml') = (2, 1, 0) → (0, 0)
(n, l, ml) → (l', ml') = (2, 1, 0) → (1, -1)
(n, l, ml) → (l', ml') = (2, 1, 0) → (1, 0)
(n, l, ml) → (l', ml') = (2, 1, 0) → (1, 1)
(n, l, ml) → (l', ml') = (2, 1, 0) → (2, -2)
(n, l, ml) → (l', ml') = (2, 1, 0) → (2, -1)
(n, l, ml) → (l', ml') = (2, 1, 0) → (2, 0)
(n, l, ml) → (l', ml') = (2, 1, 0) → (2, 1)
(n, l, ml) → (l', ml') = (2, 1, 0) → (2, 2)
(n, l, ml) → (l', ml') = (2, 1, 1) → (0, 0)
(n, l, ml) → (l', ml') = (2, 1, 1) → (1, -1)
(n, l, ml) → (l', ml') = (2, 1, 1) → (1, 0)
(n, l, ml) → (l', ml') = (2, 1, 1) → (1, 1)
(n, l, ml) → (l', ml') = (2, 1, 1) → (2, -2)
(n, l, ml) → (l', ml') = (2, 1, 1) → (2, -1)
(n, l, ml) → (l', ml') = (2, 1, 1) → (2, 0)
(n, l, ml) → (l', ml') = (2, 1, 1) → (2, 1)
(n, l, ml) → (l', ml') = (2, 1, 1) → (2, 2)

O_transition_potentials = O_transitions.get_transition_potentials(energy=100e3)
O_transition_potentials.grid.match(potential)
O_transition_potentials_array = O_transition_potentials.build()

Sr_transition_potentials = Sr_transitions.get_transition_potentials(energy=100e3)
Sr_transition_potentials.grid.match(potential)
Sr_transition_potentials_array = Sr_transition_potentials.build()

We see that including all first order transitions will require a significant amount of computations, however, we can eliminate some of the transitions on the basis of the corresponding local transition potential of a specific transition, $V_{fi, \ local}(\mathbf{k}_{\perp})$, as a fraction of the local transition potential for all transitions

\[ \hat{V}_{local}(\mathbf{k}_{\perp}) := \sum_{fi} \hat{V}_{fi, \ local}(\mathbf{k}_{\perp}) = \sum_{fi} |\hat{V}_{fi}(\mathbf{k}_{\perp})|^2 A(\mathbf{k}_{\perp}) \]

We can pick a threshold, $t$, and include the smallest set of initial and final states, $fi'$, to satisfy

\[ t < \sum_{fi'} \frac{\int \hat{V}_{fi', \ local}(\mathbf{k}_{\perp}) \mathrm{d}\mathbf{k}_{\perp}}{\int \hat{V}_{local}(\mathbf{k}_{\perp}) \mathrm{d}\mathbf{k}_{\perp}} \quad . \]

The bar plot below shows the integrated local transition potential (left axis) for each transition ordered according to their integrated intensity. The red line indicates the cumulative sum of integrated local transition potentials. From this figure we can deduce that including the seven most important transitions gives us roughly $95 \ \%$ of the local potential, which we shall be satisfied with here.

intensities = Sr_transition_potentials_array.integrated_intensities(60)

order = np.argsort(-intensities)
intensities = intensities[order]
labels = np.array(Sr_transition_potentials_array.ensemble_axes_metadata[0].values)[
    order
]

fig, ax1 = plt.subplots(figsize=(6, 4))

ax1.bar(labels, intensities)
ax1.set_xticks(labels, labels, rotation=90)
ax1.set_ylabel("$\int V_{fi, \ local} \mathrm{d}k \ [eV]$")

ax2 = ax1.twinx()
ax2.axhline(95, color="k", label="95 %")
ax2.text(labels[-1], 94, "95 %", color="k", ha="center", va="top")
ax2.plot(labels, np.cumsum(intensities / intensities.sum()) * 100, "r.-")
ax2.set_ylim([0, 105])
ax2.set_ylabel("Cumulative local potential [%]")

Text(0, 0.5, 'Cumulative local potential [%]')

../../_images/fcb5c55ea6670ec652a1527698fb958638345f3eeb542cd193da196e1055ed5d.png

We can remove the least important transitions using filter_by_intensity, reducing the number of transitions from $27$ to $9$.

filtered_Sr_transition_potentials = Sr_transition_potentials_array.filter_by_intensity(
    threshold=0.95, max_angle=60
)

len(filtered_Sr_transition_potentials)

We create a probe, a scan and a detector as we would for a normal STEM simulation. We use the FlexibleAnnularDetector as it is slightly faster than the PixelatedDetector (unlike standard multislice the detector has to be applied a significant number of times).

probe = abtem.Probe(
    semiangle_cutoff=25,
    energy=100e3,
    device="cpu",
)
probe.grid.match(potential)

scan = abtem.GridScan(
    start=(0, 0),
    end=(1 / 4, 1 / 4),
    fractional=True,
    potential=potential,
    endpoint=False,
)

detector = abtem.FlexibleAnnularDetector(to_cpu=True)

Since both the probe and the transition potentials are localized, we can further cut the computational cost by discounting excitations from atoms, where the transition potentials does not overlap significantly with the wave function. We use a rough estimate to determine the contribution to the total intensity of a site $\mathbf{r}_p$ given a wave function $\psi(\mathbf{r}_{\perp}, z)$

\[ I(\mathbf{r}_p, z) \approx \int V_{local}(\mathbf{r}_{\perp}, \mathbf{r}_{p}) |\psi(\mathbf{r}_{\perp}, z)|^2 \mathrm{d}\mathbf{r}_{\perp} \]

Given the set of all sites $\mathbf{r}_p$ within a given thickness, we will include the smallest subset those sites, $\mathbf{r}_p'$, to satisfay

(11)#\[ t < \frac{\sum_{\mathbf{r_p}'} I(\mathbf{r}_p', z)}{\sum_{\mathbf{r}_p} I(\mathbf{r}_p, z)} \quad . \]

Note that the detailed implementation differs slightly from the above description.

Given the probe wave function at position, $\mathbf{r}_p$, at the center of the unit cell (the white cross) at $z=0$, we show included Oxygen sites below for the three different thresholds. We see that at a threshold of $99 \ \%$ about half the sites in the cell are included, while at $95 \ \%$ only the column and it nearest neighbours are included. The condition in Eq. (11) may be expensive to update at every slice, hence by default it is only updated every 10 slices.

sites_all = atoms[(atoms.numbers == 8) * (atoms.positions[:, 2] < 2)]

position = sites_all[38].position[:2]
waves = probe.build(position).compute(progress_bar=False)

ts = [1.0, 0.99, 0.95]

fig, axes = plt.subplots(1, len(ts), figsize=(13, 4))
for t, ax in zip(ts, axes):
    absolute_threshold = O_transition_potentials_array.absolute_threshold(
        waves, threshold=t
    )
    sites = O_transition_potentials_array.filter_sites(
        waves, sites_all, absolute_threshold
    )
    abtem.show_atoms(sites, ax=ax, title=f"Sites $t={t}$")
    ax.scatter(*position.T, c="w", marker="x", s=60, linewidth=2)

plt.tight_layout()

../../_images/57997299d29a55b5e10b5da9b08856c4c6d08260cd8ad96baa8bb3de365f7b92.png

We run the transition_potential_scan with the Oxygen transition potential. We set the threshold for including a site to $t=0.95$, we use single-channeling here. We set a lower max_batch than normally selected to allow for better parallelization over probe positions. We will also enable task level progress to be displayed, this is only available with the local threads scheduler (the default). The simulation should take slightly more than a minute.

O_measurements = probe.transition_potential_scan(
    scan=scan,
    potential=potential,
    detectors=detector,
    transition_potentials=O_transition_potentials_array,
    double_channel=False,
    max_batch=40,
    threshold=0.95,
)

O_measurements.compute(scheduler="threads", num_workers=6)

[########################################] | 100% Completed | 71.59 ss

<abtem.measurements.PolarMeasurements at 0x7f442be5eed0>

We run an identical simulation for Strontium, which should be slightly faster than the Oxygen simulation.

Sr_measurements = probe.transition_potential_scan(
    scan=scan,
    potential=potential,
    detectors=detector,
    transition_potentials=filtered_Sr_transition_potentials,
    double_channel=False,
    threshold=0.95,
    max_batch=40,
)

Sr_measurements.compute(scheduler="threads", num_workers=6)

[########################################] | 100% Completed | 43.18 ss

<abtem.measurements.PolarMeasurements at 0x7f442bea5310>

We stack the measurements and perform the standard post-processing steps for STEM simulations.

stacked_measurements = abtem.stack((O_measurements, Sr_measurements), ("O", "Sr"))

post_processed_measurements = (
    stacked_measurements.integrate_radial(inner=0, outer=60)
    .tile(repetitions=(4, 4))
    .interpolate(sampling=0.05)
)

visualization = post_processed_measurements.show(
    explode=True, common_color_scale=True, cbar=True
)

../../_images/35222f5e2d243f01a9f7d6af1a01b9963e7dd365dc42dafb14c17d0e5439dcd5.png

The simulations for Oxygen here are made to be similar to those of Fig. 3(c) in DOI:10.1103/PhysRevB.77.184107, we can interpolate a line profile to get a comparable figure. Our simulation has slightly higher intensity presumably because they use and absorptive potential.