Aberrated Image Recovery of Ultracold Atoms

Learning Objectives

By the end of this project you should be able to:

• explain how ultracold atomic gases are produced and imaged, and why their density is often Gaussian;

• describe the resolution limits of an imaging system: the point spread function, the Airy disk, the depth of field and pixelization;

• treat an image as a high-dimensional data element, and explain why its pixels are highly correlated;

• estimate a physical quantity (the atom number) from noisy images, and use K-means clustering to try to resolve discrete values; and

• use the singular value decomposition (SVD) to compress images and identify which features survive compression.

Overview

Over the past 20–30 years the groundwork has been laid for precise experimental control of atomic gases at ultracold temperatures. Ultracold atom experiments exploit the quantum mechanics of their underlying atomic systems for a diverse set of applications, ranging from the simulation of computationally difficult physics problems [1] to searches for new physics beyond the Standard Model [2]. Most experiments in this field rely on directly imaging the ultracold gas, to extract the size and shape of the atomic number density distribution.

The goal of this project is to introduce you to some relevant image-processing techniques, and to get familiar with an image as a data element. As you will see, images are a kind of data with a very large number of features, almost all of which are highly correlated within a region of interest. Of particular interest is how a real imaging system distorts the information contained in an image, and how those effects can be unfolded from the data to recover information about the true density distribution.

Capturing all possible kinds of aberrations and noise in real images of ultracold atoms is outside the scope of the simulated data in this project. Instead, we limit ourselves to a few key effects:

• optical aberrations, in the form of defocus and primary spherical aberration;

• pixelization from the finite resolution of the detector; and

• in some datasets, a simple model of number noise.

Data Sources

• Course mirror: https://courses.physics.illinois.edu/phys498mlp/sp2025/data/projects/AberratedImages/data.tgz (1 MB). The code below downloads and extracts it.

The archive contains two folders of simulated images, stored as TIFF files of floating-point values, plus the notebook ImageSimulation.ipynb that generated them (also available in the course repository). You are encouraged to read that notebook to see how the simulation works, and even to generate your own datasets.

All images simulate a thermal gas of \(^{87}\)Rb atoms imaged with \(\lambda = 780\) nm light onto a camera with pixel size \(\text{PS} = 13\) µm. The true parameters of each image are encoded in its file name, e.g. HiNA.atom_num=100000.mu=(0.0...0.0...0.0).sigma=(1.0...1.0...20.0).Z04=0.2.seed=0.tiff: the imaging system, the atom number \(N\), the centre \(\mu\) and widths \(\sigma = (\sigma_x, \sigma_y, \sigma_z)\) of the Gaussian density in µm (with \(z\) the imaging axis), the spherical-aberration coefficient Z04, and the random seed of the noise.

Dataset	Imaging system	Images	Image size	Gas widths \(\sigma\) (µm)	Atom numbers
Number_estimate	Low NA: NA = 0.2, magnification M = 4, no spherical aberration (Z04 = 0)	60	31 × 31 pixels	(10, 10, 10)	9,000 / 10,000 / 11,000 and 99,000 / 100,000 / 101,000; 10 noise seeds each
Heavy_tails	High NA: NA = 0.4, M = 20, Z04 = 0.2	18	185 × 185 pixels	\(\sigma_{xy}\) ∈ {1, 2.5, 5, 10, 20}, with \(\sigma_z = \sigma_{xy}\) or \(\sigma_z = 20\)	100,000; 2 noise seeds each

Noise model. Every image is normalized so that its pixels sum to the atom number \(N\). Then noise is added to each pixel with a standard deviation equal to the square root of the pixel value (see ImSim in ImageSimulation.ipynb). This is a toy model for fluctuations in the measured density. It can produce small negative pixel values, which were deliberately left in.

Field of view. The images cover ±50 µm (Number_estimate) and ±60 µm (Heavy_tails) in the object plane, and the simulation includes the gas over about ±3σ along the imaging axis. Even the largest gases are therefore fully contained in the images. (The datasets were regenerated in 2026 with a larger field of view than the original ±10 µm.)

%matplotlib inline
import glob
import os
import re
import tarfile
import urllib.request

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import skimage.io

DATA_URL = "https://courses.physics.illinois.edu/phys498mlp/sp2025/data/projects/AberratedImages/data.tgz"

# Download and extract the data once
if not os.path.isdir("Number_estimate"):
    urllib.request.urlretrieve(DATA_URL, "data.tgz")
    with tarfile.open("data.tgz") as tar:
        try:
            tar.extractall(".", filter="data")   # safe extraction (Python >= 3.12, and recent 3.8-3.11)
        except TypeError:
            tar.extractall(".")

PATTERN = re.compile(r"(?P<system>\w+)\.atom_num=(?P<N>\d+)\."
                     r"mu=\((?P<mu_x>[-\d.]+?)\.\.\.(?P<mu_y>[-\d.]+?)\.\.\.(?P<mu_z>[-\d.]+?)\)\."
                     r"sigma=\((?P<sigma_x>[\d.]+?)\.\.\.(?P<sigma_y>[\d.]+?)\.\.\.(?P<sigma_z>[\d.]+?)\)\."
                     r"Z04=(?P<Z04>[\d.]+)\.seed=(?P<seed>\d+)\.tiff")

def load_dataset(dataset_name):
    """Load all images in a dataset folder into a DataFrame: one row per image,
    with the true parameters parsed from the file name and the image itself in column 'image'."""
    rows = []
    for path in sorted(glob.glob(os.path.join(dataset_name, "*.tiff"))):
        m = PATTERN.fullmatch(os.path.basename(path))
        row = {k: float(v) for k, v in m.groupdict().items() if k != "system"}
        row["system"] = m["system"]
        row["image"] = skimage.io.imread(path)
        rows.append(row)
    df = pd.DataFrame(rows)
    return df.astype({"N": int, "seed": int})

number_df = load_dataset("Number_estimate")
tails_df = load_dataset("Heavy_tails")
print(len(number_df), "Number_estimate images of shape", number_df.image[0].shape)
print(len(tails_df), "Heavy_tails images of shape", tails_df.image[0].shape)
tails_df.drop(columns="image").head()

Questions

Answer each question in a markdown cell directly below it (add code cells where needed). For the background questions, aim for one to two well-organized paragraphs, and cite your sources. Your notebook should run from top to bottom without errors.

Part A: Background

Question 01

What is an ultracold atomic gas? How cold are we talking? How do experiments bring atoms down to these temperatures?

When you research these experiments, you will find that alkali atoms (Rb, K, Na, etc.) are the most commonly used species. Why these atoms? (Hint: recall high-school chemistry. What do these atoms have in common with everyone’s favourite quantum system, the hydrogen atom?)

A helpful resource is Foot’s Atomic Physics [3].

Question 02

To measure the number density distribution of an ultracold gas, we typically release the gas from its trap and let it expand in free fall before taking an image after some time of flight (TOF). The density distribution in such an image often has a Gaussian shape. Why is that? (Hint: think statistical mechanics. What is the velocity distribution of a thermal gas, and how does the cloud size grow with TOF?) How can a TOF measurement be used to determine the temperature of the gas?

Question 03

Real imaging systems have a finite resolution, set by the diffraction limit. Consider a point source emitting a spherical electromagnetic wave. A real lens propagates the field from the object plane to the image plane, producing a pattern called the imaging response function, or point spread function (PSF), of the system. Any object of finite size can be thought of as a collection of point sources, so its image is the sum of the PSFs of all its points, i.e. a convolution.

For a circular lens, the PSF at best focus, for a point source on the optical axis, is the Airy disk, \(I(\rho) = I_0 \left[2J_1(\rho)/\rho\right]^2\). Here \(\rho = 2\pi\,\text{NA}\, r/\lambda\) is a dimensionless radial coordinate, and its first zero is at \(\rho \approx 3.83\).

1. Calculate the fraction of the total irradiance (power) contained in the central peak, inside the first zero.

2. Approximate the Airy disk by a Gaussian, \(I_0 \exp(-2\rho^2/w^2)\), with the same peak irradiance and the same total power. How is the first zero of the Airy disk related to the \(1/e^2\) radius \(w\) of the Gaussian? Compare with a least-squares Gaussian fit to the central peak.

3. Convert the radius of the first zero to µm for the two imaging systems used in this project (NA = 0.4 and NA = 0.2, \(\lambda = 780\) nm).

In addition to the links above, a useful resource is Goodman’s Introduction to Fourier Optics [4].

Question 04

The resolution in Question 03 describes the size of the image of a point source in the plane transverse to the optical axis. What about along the optical axis? This axial resolution is often called the depth of field (DOF), although different sources define its value differently.

• How are the DOF and the transverse resolution related to each other, and how does each scale with NA? No derivation is necessary, although you may give one if you restrict yourself to the axial PSF along \(\rho = 0\).

• Estimate the DOF for the two imaging systems in this project.

• Each camera pixel corresponds to \(\text{PS}/M\) in the object plane. For each imaging system, which limits the resolution of the images: diffraction or pixel size?

• Which of the gases in the Heavy_tails dataset are thicker along \(z\) than the DOF? What do you expect their images to look like?

---

Part B: Images as data

The remaining questions use the Number_estimate and Heavy_tails datasets described in the Data Sources section.

• With the Number_estimate dataset, you will explore how noise affects your ability to estimate the total number of atoms in an image.

• With the Heavy_tails dataset, you will test the ability of a simple linear decomposition (SVD) to reduce the dimensionality of the data, and check whether the “compressed” images still retain the relevant features of the gas. Gases with \(\sigma_z \gtrsim\) DOF acquire a noticeable tail in the image, which departs from the true Gaussian density distribution. Of interest is how this tail is carried through to the compressed and reconstructed images.

Question 05

Display a selection of images from both datasets, with axes in µm in the object plane (use the pixel size \(\text{PS}/M\)). Plot horizontal line profiles through the centre of each Heavy_tails image on a logarithmic scale.

• Check that the sum over the pixels of each image is close to its atom number \(N\).

• Check that every gas is fully contained in its image, by comparing the pixel values at the edge of each image with its peak. Why is it important for this project that the field of view is much larger than the gas?

• Compare the images with \(\sigma_{xy} = 1\) µm and \(\sigma_z = 1\) or 20 µm. What differences do you see, and why?

Question 06

The data in this project are images: 2D arrays of greyscale values representing the irradiance of the image of the gas. If each image is treated as a single sample, what is the dimension of the data for each dataset?

Suppose the question we are interested in is: is the number of atoms in an image greater than some value \(X\)? If none of the pixels were correlated, roughly how many samples would we need to sample the space of possible images? Assume, for example, that each pixel can take one of only 10 distinguishable values. Why is the true number of degrees of freedom of these images much smaller?

Question 07

Take the images in the Number_estimate dataset with a true atom number of order \(N \sim 10^4\). Estimate the atom number in each image by summing over all its pixels, and draw an appropriately binned histogram of these estimates. Repeat this for the \(N \sim 10^5\) subset of the dataset.

Note: in a real experiment, images of ultracold gases are not recorded directly in terms of the number density, but in terms of a quantity proportional to the column density called the optical depth (OD). See ImageSimulation.ipynb for details.

How are these two histograms different from each other? Compare the spread of the estimates for each true \(N\) with what you expect from the noise model described in the Data Sources section. How does the relative uncertainty of the atom number scale with \(N\)? With only 10 images per atom number, your estimate of the spread is itself uncertain. Estimate that uncertainty, e.g. with a bootstrap, before you conclude whether the data agree with the noise model.

Question 08

Use K-means clustering on the 1D atom-number estimates from Question 07 to try to recover the true atom number of each image. Treat the \(N \sim 10^4\) and \(N \sim 10^5\) images separately, with 3 clusters each.

• Compare the cluster assignments with the true atom numbers (e.g. with a confusion matrix, or sklearn.metrics.adjusted_rand_score).

• Are you able to clearly resolve the atom number in each subset? Explain the difference between the subsets, in terms of the spacing between the true atom numbers and the noise you found in Question 07.

• Roughly how many images per atom number, or how much less noise, would you need to resolve the \(N \sim 10^5\) subset as well as the \(N \sim 10^4\) one?

Question 09

Use the singular value decomposition (SVD) to compress the Heavy_tails images. Write each 185 × 185 image as \(M = U S V^T\) and keep only the \(k\) largest singular values to form a rank-\(k\) approximation \(M_k\).

1. Plot the singular values of each image on a logarithmic scale. A 2D Gaussian with independent \(x\) and \(y\) widths is an exactly rank-1 matrix. Show why. What, then, do the smaller singular values represent?

2. Reconstruct the images for a handful of \(k\) (e.g. 1, 2, 3, 5, 10), and compare line profiles (on a log scale) of the reconstructions with the original images.

3. Focus on pairs of images with the same \(\sigma_{xy}\) but different \(\sigma_z\) (\(\sigma_z = \sigma_{xy}\) and \(\sigma_z = 20\) µm). Up to what \(k\) do the two images in each pair appear to have the same-sized tail in the compressed representation? At what \(k\) does the difference become visible? Compare the singular-value spectra of the two images in each pair. For which values of \(\sigma_{xy}\) is there a visible difference between the two images at all? Explain why, by comparing \(\sigma_z\) with the DOF and \(\sigma_{xy}\) with the resolution of the imaging system (Questions 03–04).

4. Optional: instead of compressing each image separately, flatten all 18 images into vectors of 185 × 185 = 34,225 pixels and perform an SVD (PCA) of the whole dataset. What do the leading components look like?

References

[1] S. S. Kondov, W. R. McGehee, J. J. Zirbel, B. DeMarco, “Three-dimensional Anderson localization of ultracold matter,” Science 334, 66 (2011). doi:10.1126/science.1209019

[2] W. B. Cairncross et al., “Precision measurement of the electron’s electric dipole moment using trapped molecular ions,” Phys. Rev. Lett. 119, 153001 (2017). doi:10.1103/PhysRevLett.119.153001

[3] C. J. Foot, Atomic Physics (Oxford University Press, 2005).

[4] J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company, 2005).

Acknowledgements

• Initial version: Max Gold with some guidance from Mark Neubauer

• Revised: 2026

© Copyright 2026