An Example XPS Analysis#

In this example analysis, we are going to look at how some elements of PyARPES fit together. We’ll use the example_data.nano_xps dataset.

[1]:

import numpy as np
import arpes
from arpes.io import example_data
from arpes.plotting.dos import plot_core_levels
import matplotlib.pyplot as plt

xps = example_data.nano_xps

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
xps.sum("eV").S.plot(ax=axes[0])
plot_core_levels(xps.sum(["x", "y"]), ax=axes[1])
plt.tight_layout()

Activating auto-logging. Current session state plus future input saved.
Filename       : logs/unnamed_2026-03-24_23-30-08.log
Mode           : backup
Output logging : False
Raw input log  : False
Timestamping   : False
State          : active
../_images/notebooks_full-analysis-xps_1_14.png

Decomposition Analysis#

Let’s look at this XPS data by performing a PCA-based decomposition across the spatial axes. This, in combination with manual exploration, will give us some sense of what’s going on across the sample surface.

[2]:

from arpes.analysis.decomposition import pca_along
import xarray as xr

n_components = 5
data, pca = pca_along(xps.spectrum, ["x", "y"], n_components=n_components)

[3]:

fig, ax = plt.subplots(n_components, 2, figsize=(8, 4 * n_components))

for component in range(n_components):
    data.isel(components=component).S.plot(ax=ax[component, 0])
    ax[component, 0].set_title(f"Component {component}")

    xr.DataArray(pca.components_[component], {"eV": xps.eV.values}, ["eV"]).plot(
        ax=ax[component, 1]
    )

plt.tight_layout()
../_images/notebooks_full-analysis-xps_4_0.png

Let’s parse this. Each row above has a spatial map of the coefficient in the decomposition (left) and the XPS spectrum corresponding to that component (right). Very intense locations in the spatial maps correspond to places where the spectrum on the right is a good representation to the data (red) or to its negative (blue).

Based on the images, we can see that the first four components (0 through 3) explain almost all the variation in the data.

The first component corresponds to a \(\text{WS}_2\) 4f core level spectrum, so regions of high intensity here indicate the presence of \(\text{WS}_2\).

The next component has a wide background, but it still contains the imprint of the core level peaks. What’s going on here? Note that in the Component 0 image, the area surrounding the \(\text{WS}_2\) has a negative coefficient, but if we look at the average below, corresponding to an area to the left of the \(\text{WS}_2\) patch, there are no peaks. The decomposition has learned that it needs to add some core-level like peaks back into the data in order to recover these regions. Clearly, this is not an ideal decomposition, and it points to the need to interpret your analysis. There are other decompositions we could perform which have better biases for spectroscopic data.

Moving to the third component we see some internal structure on the \(\text{WS}_2\) sample. By looking at the corresponding PCA component, this appears to represent a core level shift toward lower binding energy.

The fourth and final component we will interpret (Component 3) again explains internal structure on the \(\text{WS}_2\). Looking at the corresponding spectrum, it looks like these regions correspond to wider core levels as measured by photoemission.

Selecting data using the PCA decomposition#

Suppose we want to continue analysis using only the core levels for the \(\text{WS}_2\) region, as identified by Component 0. We can do this by masking and selecting the data before averaging.

[4]:

ws2_mask = data.isel(components=0) > 500
substrate_mask = data.isel(components=1) > 1000

fig, ax = plt.subplots(2, 2, figsize=(10, 10))
ws2_mask.S.plot(ax=ax[0, 0])
substrate_mask.S.plot(ax=ax[0, 1])

xps.spectrum.where(ws2_mask).mean(["x", "y"]).S.plot(ax=ax[1, 0])
xps.spectrum.where(substrate_mask).mean(["x", "y"]).S.plot(ax=ax[1, 1])

ax[0, 0].set_title("Component 0 Mask")
ax[0, 1].set_title("Component 1 Mask")
ax[1, 0].set_title("Component 0 Masked XPS")
ax[1, 1].set_title("Component 1 Masked XPS")

plt.tight_layout()
../_images/notebooks_full-analysis-xps_7_0.png

Improving our masks#

We can see there’s still some \(\text{WS}_2\) contamination in the substrate_mask. We can improve this by masking according to both of the first two components.

[5]:

ws2_mask = data.isel(components=0) > 500
substrate_mask = (data.isel(components=1) > 1000) & ~ws2_mask

fig, ax = plt.subplots(2, 2, figsize=(10, 10))
ws2_mask.S.plot(ax=ax[0, 0])
substrate_mask.S.plot(ax=ax[0, 1])

xps.spectrum.where(ws2_mask).mean(["x", "y"]).S.plot(ax=ax[1, 0])
xps.spectrum.where(substrate_mask).mean(["x", "y"]).S.plot(ax=ax[1, 1])

ax[0, 0].set_title("Component 0 Mask")
ax[0, 1].set_title("Component 1 Mask")
ax[1, 0].set_title("Component 0 Masked XPS")
ax[1, 1].set_title("Component 1 Masked XPS")

plt.tight_layout()
../_images/notebooks_full-analysis-xps_9_0.png

Looking at the wider vs narrower peak regions#

[6]:

from arpes.config import use_tex

use_tex(enable=False)

ws2_mask = data.isel(components=0) > 500
wide_peak_mask = (data.isel(components=3) > 500) & ws2_mask

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
wide_peak_mask.S.plot(ax=ax[0])

xps.spectrum.where(ws2_mask & ~wide_peak_mask).mean(["x", "y"]).S.plot(ax=ax[1], label="ws2_mask")
xps.spectrum.where(wide_peak_mask).mean(["x", "y"]).S.plot(ax=ax[1], label="wide_peak_mask")

ax[1].legend()

ax[0].set_title("Wide Peak Mask")
ax[1].set_title("XPS Comparison")

plt.tight_layout()
../_images/notebooks_full-analysis-xps_11_0.png

These peaks look very similar, but sure enough the ones captured by wide_peak_mask are wider, especially to high binding energy.

Let’s look at refining this analysis now with some curve fitting.

Curve Fitting#

First off, let’s get a general model working using a single XPS curve.

[7]:

from lmfit.models import GaussianModel, LinearModel

test_curve = xps.spectrum.where(ws2_mask & ~wide_peak_mask).mean(["x", "y"]).sel(eV=slice(-36, -31))

test_model = LinearModel() + GaussianModel(prefix="a_") + GaussianModel(prefix="b_")

target = test_curve - test_curve.min()

x = target.coords["eV"].values
roi1 = x < 45
roi2 = x > -33.5

params = GaussianModel(prefix="a_").guess(target.values[roi1], target.coords["eV"].values[roi1])
params.update(
    GaussianModel(prefix="b_").guess(target.values[roi2], target.coords["eV"].values[roi2])
)


result = target.S.modelfit("eV", test_model, params=params)

result.modelfit_results.item().plot()
result.modelfit_results.item()

[7]:

Fit Result

Model: ((Model(linear) + Model(gaussian, prefix='a_')) + Model(gaussian, prefix='b_'))

../_images/notebooks_full-analysis-xps_13_1.png

This looks reasonably good, but we can improve the simplicity and quality of the model if we just remove an estimated background first. Let’s see how that looks

[8]:

from lmfit.models import GaussianModel
from arpes.analysis.shirley import remove_shirley_background
from arpes.fits.utilities import result_to_hints

mask = ws2_mask | (data.isel(components=3) > 800)
test_curve = xps.spectrum.where(mask).mean(["x", "y"]).sel(eV=slice(-36, -31))
test_curve = remove_shirley_background(test_curve)

test_model = GaussianModel(prefix="a_") + GaussianModel(prefix="b_")
x = test_curve.coords["eV"].values
roi_a = x < -33.5
roi_b = x > -33.5

params = GaussianModel(prefix="a_").guess(test_curve.values[roi_a], x[roi_a])
params.update(GaussianModel(prefix="b_").guess(test_curve.values[roi_b], x[roi_b]))

result = test_curve.S.modelfit("eV", test_model, params=params)

result.modelfit_results.item().plot()
result.modelfit_results.item()

[8]:

Fit Result

Model: (Model(gaussian, prefix='a_') + Model(gaussian, prefix='b_'))

../_images/notebooks_full-analysis-xps_15_1.png

Looking at the residual this is quite a bit better. Now let’s perform fitting across the entire \(\text{WS}_2\) region.

[9]:

# subtract the Shirley background, calculated independently at each (x,y)
bkg_removed_xps = remove_shirley_background(xps.spectrum.sel(eV=slice(-36, -31)))

# first mask
masked_xps = bkg_removed_xps.where(mask)

[10]:

# Performs ~500 curve fits... make a selection if you don't want to wait a few seconds
fit_results = masked_xps.fillna(0).S.modelfit("eV", test_model, params=params, progress=True)

Interepting the fit quality#

Let’s have a look at the five worst fits now.

[11]:

worst = fit_results.modelfit_results.F.worst_fits().values

for i in range(5):
    fig = plt.figure(figsize=(4, 4))
    worst[i].plot(fig=fig)
    for ax in fig.axes:
        ax.set_title("")
../_images/notebooks_full-analysis-xps_20_0.png
../_images/notebooks_full-analysis-xps_20_1.png
../_images/notebooks_full-analysis-xps_20_2.png
../_images/notebooks_full-analysis-xps_20_3.png
../_images/notebooks_full-analysis-xps_20_4.png

The first three fits have quality issues, because the background looks not to be appropriately captured and the peak shape is different from the rest. Three problematic fits out of several hundred is not so bad, so we can continue.

For proper analysis you might want to exclude these points.

[12]:

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
fit_results.modelfit_results.F.p("a_amplitude").T.S.plot(ax=ax[0])
fit_results.modelfit_results.F.p("b_amplitude").T.S.plot(ax=ax[1], vmin=0)
plt.tight_layout()
../_images/notebooks_full-analysis-xps_22_0.png 
[13]:

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
fit_results.modelfit_results.F.p("a_center").T.S.plot(ax=ax[0], vmin=-34.8, vmax=-34.1)
fit_results.modelfit_results.F.p("b_center").T.S.plot(ax=ax[1], vmin=-32.7, vmax=-32.05)
plt.tight_layout()
../_images/notebooks_full-analysis-xps_23_0.png 
[14]:

fig, ax = plt.subplots(1, 2, figsize=(12, 5))
fit_results.modelfit_results.F.p("a_sigma").T.S.plot(ax=ax[0], vmax=0.22, vmin=0.14)
fit_results.modelfit_results.F.p("b_sigma").T.S.plot(ax=ax[1], vmax=0.22, vmin=0.14)
plt.tight_layout()
../_images/notebooks_full-analysis-xps_24_0.png 
[15]:

plt.scatter(
    fit_results.modelfit_results.F.p("a_amplitude").values.ravel(),
    fit_results.modelfit_results.F.p("b_amplitude").values.ravel(),
    color=(0, 0, 0, 0.05),
)
plt.gca().set_aspect(1)
plt.gca().set_xlabel("Lower Peak Intensity")
plt.gca().set_ylabel("Upper Peak Intensity")
plt.gca().set_xlim([0, 1600])
plt.gca().set_ylim([0, 2400])

[15]:

(0.0, 2400.0)
../_images/notebooks_full-analysis-xps_25_1.png 
[16]:

plt.scatter(
    fit_results.modelfit_results.F.p("a_center").values.ravel(),
    fit_results.modelfit_results.F.p("b_center").values.ravel(),
    color=(0, 0, 0, 0.05),
)
plt.gca().set_aspect(1)
plt.gca().set_xlabel("Lower Peak BE")
plt.gca().set_ylabel("Upper Peak BE")

[16]:

Text(0, 0.5, 'Upper Peak BE')
../_images/notebooks_full-analysis-xps_26_1.png 
[17]:

plt.scatter(
    fit_results.modelfit_results.F.p("a_center").values.ravel(),
    fit_results.modelfit_results.F.p("a_sigma").values.ravel(),
    color=(0, 0, 0, 0.1),
)
plt.gca().set_aspect(1)
plt.gca().set_xlabel("Lower Peak BE")
plt.gca().set_ylabel("Upper Peak Width")

[17]:

Text(0, 0.5, 'Upper Peak Width')
../_images/notebooks_full-analysis-xps_27_1.png

The two cohorts in each of the above plots are the piece of \(\text{WS}_2\) on substrate and on plain Si respectively, confirming earlier expectations from the PCA analysis.

Exercises#

  1. Perform a selection of the data according to the fitting values and plot the corresponding mask. What does the collection of points for which Lower Peak BE > -34.4eV look like?

  2. Perform a selection of the data for the decompositions above and plot their corresponding average XPS curves. How do these results compare to the PCA results we found before?

  3. What might you conclude about the sample and experimental conditions given the extracted peak width map?