Single sideband imaging ----------------------- High-level overview ^^^^^^^^^^^^^^^^^^^ Single sideband (SSB) imaging extracts phase information from 4D-STEM data through a direct, non-iterative process. Here's how it works in four steps: **1. Collect 4D-STEM data** Scan a focused electron probe across your sample in a grid pattern. At each scan position, record a full 2D diffraction pattern on the detector. This creates a 4D dataset :math:`I(\mathbf{k}; \mathbf{R})` where :math:`\mathbf{R}` is the scan position and :math:`\mathbf{k}` is the detector coordinate. **2. Separate sidebands in Fourier space** Fourier transform over all scan positions: :math:`G(\mathbf{k}, \mathbf{Q}) = \mathcal{F}_{\mathbf{R} \to \mathbf{Q}}[I(\mathbf{k}; \mathbf{R})]`. This separates the signal into three components: - Center disk at :math:`\mathbf{Q} = 0` (no phase information) - Right sideband at :math:`\mathbf{Q} = +\mathbf{k}` (contains object phase) - Left sideband at :math:`\mathbf{Q} = -\mathbf{k}` (conjugate of right sideband) **3. Isolate and divide by probe** Select one sideband using a circular aperture, shift it to the center, then divide by the probe intensity: :math:`\tilde{O}(\mathbf{k}) = S(\mathbf{k}) / |P(\mathbf{k})|^2`. This removes the probe's influence and isolates the object's Fourier transform. The probe :math:`P(\mathbf{k})` must be either measured from a vacuum region or constructed from known aberration coefficients. **4. Transform back to get phase** Inverse Fourier transform the isolated sideband back to real space: :math:`o(\mathbf{r}) = \mathcal{F}_{\mathbf{k} \to \mathbf{r}}^{-1}[\tilde{O}(\mathbf{k})]`. Take the imaginary part to extract the phase map :math:`\phi(x,y) = \text{Im}[o(\mathbf{r})]`, which represents the projected potential of your sample. **Why it's fast**: All operations are linear (Fourier transforms, multiplication, division). No iterative loops, no optimization, no convergence issues. **Key requirement**: Works best for thin samples under the weak-phase approximation where the object transmission is :math:`O(x,y) \approx 1 + i\phi(x,y)`. The following sections explain the physics, mathematics, and implementation details behind each of these steps. .. _ssb_motivation: Why SSB? ^^^^^^^^ Ptychography is a powerful phase retrieval technique, but it has limitations: - **Slow reconstruction**: Iterative algorithms (ePIE, rPIE) require hundreds to thousands of iterations, making reconstruction computationally expensive - **Initialization dependence**: Poor initial guesses for the probe or object can lead to slow convergence or failure - **Probe uncertainty**: The probe function must be recovered simultaneously with the object, adding complexity Single sideband (SSB) imaging addresses these challenges by providing: 1. **Fast phase extraction**: Direct algebraic reconstruction without iteration 2. **Simple probe measurement from reference samples**: With a known reference (vacuum, amorphous film), SSB can quickly measure the probe including aberrations. This measured probe can then be used for unknown samples. 3. **Better initialization for iterative ptychography**: SSB gives approximate solutions quickly. Both the object phase φ(x,y) and the measured probe P(k) can initialize ePIE/rPIE, dramatically accelerating convergence compared to random initialization. SSB is particularly valuable in two scenarios: - **Standalone phase imaging**: For thin, weakly scattering samples where the weak-phase approximation holds, SSB provides quick, high-quality phase maps without iteration. This is often sufficient for analysis. - **Hybrid workflow**: SSB provides fast initial estimates → iterative ptychography refines to higher accuracy. SSB handles the "easy part" (getting close to the solution) in seconds, while ptychography handles the "hard part" (fine-tuning beyond weak-phase approximation, correcting probe variations, improving resolution). This combination achieves both speed and accuracy for challenging samples. Physical setup and approximations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. _weak_phase_approximation: **Weak-phase approximation** The object's transmission function under the weak-phase approximation is: .. math:: O(x,y) = e^{i\phi(x,y)} \approx 1 + i\phi(x,y) This approximation assumes that the object changes the probe primarily by altering its phase rather than its amplitude, and that this modification is mild. It applies to thin specimens (typically < 10 nm for biological materials, < 50 nm for light elements) where multiple scattering is negligible. **Exit wave formation** The exit wave at probe position :math:`\mathbf{R} = (R_x, R_y)` in the scan is: .. math:: \psi_{\text{exit}}(x,y;\mathbf{R}) = P(x - R_x, y - R_y) \cdot O(x,y) Substitution of the weak-phase expansion yields: .. math:: \psi_{\text{exit}}(x,y;\mathbf{R}) \approx P(x - R_x, y - R_y) + i\phi(x,y) P(x - R_x, y - R_y) The first term represents the unmodified probe. The second term contains the object information. .. _sideband_formation: Sideband formation in reciprocal space ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **The 4D-STEM dataset structure** SSB uses the **entire 4D-STEM dataset** to extract phase. The dataset is a 4D array: .. math:: I(x_{\text{scan}}, y_{\text{scan}}, k_x, k_y) where :math:`(x_{\text{scan}}, y_{\text{scan}})` are the probe scan positions and :math:`(k_x, k_y)` are the detector coordinates in reciprocal space. At each scan position, a full 2D diffraction pattern is recorded. The sidebands do not appear in individual diffraction patterns. They emerge only after **Fourier transforming over the scan positions** :math:`(x_{\text{scan}}, y_{\text{scan}})`, combining information from all probe positions across the field of view. This is the key difference between SSB and simply looking at a single diffraction pattern. **Formation of interference terms at one probe position** To understand how the sidebands form, first consider what happens at a single probe position :math:`\mathbf{R} = (R_x, R_y)`. Fourier transformation of the exit wave to reciprocal space shows that the probe term, when shifted, acquires a phase ramp: .. math:: \mathcal{F}\{P(x - R_x, y - R_y)\} = P(k_x, k_y) e^{-i(k_x R_x + k_y R_y)} .. note:: **Phase ramp**: The factor :math:`e^{-i(k_x R_x + k_y R_y)}` is the Fourier shift theorem in action—when a function shifts in real space by :math:`\mathbf{R}`, its Fourier transform acquires a linear phase proportional to :math:`\mathbf{R}`. This encodes the position information. The object term in real space (product of :math:`\phi` and the shifted probe) becomes a convolution in Fourier space: .. math:: \mathcal{F}\{\phi(x,y) P(x - R_x, y - R_y)\} = [\tilde{O}(k_x, k_y) * P(k_x, k_y)] e^{-i(k_x R_x + k_y R_y)} Taking the squared magnitude to get the measured intensity in reciprocal space at probe position :math:`\mathbf{R}`: .. math:: I(\mathbf{k}; \mathbf{R}) = \left| P(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{R}} + [\tilde{O}(\mathbf{k}) * P(\mathbf{k})] e^{-i\mathbf{k} \cdot \mathbf{R}} \right|^2 Expanding this produces four terms: .. math:: I(\mathbf{k}; \mathbf{R}) = |P(\mathbf{k})|^2 + |\tilde{O}(\mathbf{k}) * P(\mathbf{k})|^2 + P(\mathbf{k}) [\tilde{O}^*(\mathbf{k}) * P^*(\mathbf{k})] e^{+i\mathbf{k} \cdot \mathbf{R}} + P^*(\mathbf{k}) [\tilde{O}(\mathbf{k}) * P(\mathbf{k})] e^{-i\mathbf{k} \cdot \mathbf{R}} The four terms are: 1. **Probe autocorrelation**: :math:`|P(\mathbf{k})|^2` (no :math:`\mathbf{R}` dependence) 2. **Object autocorrelation**: :math:`|\tilde{O}(\mathbf{k}) * P(\mathbf{k})|^2` (weak, second-order) 3. **Cross term with** :math:`e^{+i\mathbf{k} \cdot \mathbf{R}}` (positive phase ramp) 4. **Cross term with** :math:`e^{-i\mathbf{k} \cdot \mathbf{R}}` (negative phase ramp) The crucial observation is that the cross terms carry opposite phase ramps that encode the probe position :math:`\mathbf{R}`. As the probe scans across the sample, :math:`\mathbf{R}` changes, making these phase ramps vary from one diffraction pattern to the next. **How the full 4D dataset creates sidebands** Now consider the entire 4D dataset. The measurement consists of :math:`I(\mathbf{k}; \mathbf{R})` for many scan positions :math:`\mathbf{R}`. At each position, the intensity has the same four-term structure shown above, but the phase ramps :math:`e^{\pm i\mathbf{k} \cdot \mathbf{R}}` have different values because :math:`\mathbf{R}` changes. **Why Fourier transform over scan positions?** Look at the cross terms: they contain :math:`e^{+i\mathbf{k} \cdot \mathbf{R}}` and :math:`e^{-i\mathbf{k} \cdot \mathbf{R}}`. These are oscillating functions of :math:`\mathbf{R}` with frequency proportional to :math:`\mathbf{k}`. **Concrete example**: Say one detector pixel has :math:`\mathbf{k} = (0.5, 0)` Å⁻¹. As the probe scans in a line from :math:`R_x = 0` to :math:`R_x = 100` Å, the phase ramp :math:`e^{i k_x R_x} = e^{i \cdot 0.5 \cdot R_x}` completes :math:`0.5 \times 100 / (2\pi) \approx 8` full oscillations. Different detector pixels :math:`\mathbf{k}` have different oscillation rates. The raw data :math:`I(\mathbf{k}; \mathbf{R})` is a messy sum of: - Constant term (no oscillation): :math:`|P|^2` - Fast oscillating term at frequency :math:`+\mathbf{k}`: cross term with :math:`e^{+i\mathbf{k} \cdot \mathbf{R}}` - Fast oscillating term at frequency :math:`-\mathbf{k}`: cross term with :math:`e^{-i\mathbf{k} \cdot \mathbf{R}}` Fourier transforming over :math:`\mathbf{R}` **unmixes these components** by their oscillation frequency, just like Fourier analysis separates a musical chord (440 Hz + 550 Hz + 660 Hz) into individual notes. **Result**: After FT, the term oscillating at frequency :math:`+\mathbf{k}` appears at :math:`\mathbf{Q} = +\mathbf{k}`, the one at :math:`-\mathbf{k}` appears at :math:`\mathbf{Q} = -\mathbf{k}`, and the constant term stays at :math:`\mathbf{Q} = 0`. The four overlapping terms **spatially separate** into three distinct regions that can be isolated with apertures. **The key SSB operation**: Fourier transform over all scan positions :math:`\mathbf{R}` to create a new representation of the data. We define: .. math:: G(\mathbf{k}, \mathbf{Q}) = \int I(\mathbf{k}; \mathbf{R}) e^{-i\mathbf{Q} \cdot \mathbf{R}} d\mathbf{R} .. important:: **This is THE key equation!** :math:`I(\mathbf{k}; \mathbf{R})` contains hidden sidebands as oscillating terms :math:`e^{\pm i\mathbf{k} \cdot \mathbf{R}}`. Applying this Fourier transform **unmixes** those oscillations by frequency, causing the sidebands to **appear** at specific :math:`\mathbf{Q}` locations in :math:`G(\mathbf{k}, \mathbf{Q})`. The sidebands were always present in :math:`I`, just hidden—the FT makes them visible and separable. **What this function does**: - **Input**: :math:`I(\mathbf{k}; \mathbf{R})` - intensity at detector pixel :math:`\mathbf{k}` when probe is at scan position :math:`\mathbf{R}` - **Operation**: For a fixed detector pixel :math:`\mathbf{k}`, take the intensity pattern across all scan positions :math:`\mathbf{R}` and compute its Fourier transform - **Output**: :math:`G(\mathbf{k}, \mathbf{Q})` - frequency spectrum showing which object periodicities (encoded by :math:`\mathbf{Q}`) contribute to that detector pixel :math:`\mathbf{k}` where: - :math:`\mathbf{k}`: Detector pixel (stays as an index) - :math:`\mathbf{Q}`: Object spatial frequency (replaces :math:`\mathbf{R}` after transformation) - :math:`e^{-i\mathbf{Q} \cdot \mathbf{R}}`: Complex exponential that extracts frequency component :math:`\mathbf{Q}` from the scan pattern **Physical intuition for Q**: Consider an object phase :math:`\phi(x,y)` with a periodic stripe pattern with spacing :math:`d`. As the probe scans across these stripes, the detector sees a periodic modulation in intensity with the same period :math:`d`. Fourier transforming over the scan positions :math:`\mathbf{R}` produces a peak at :math:`|\mathbf{Q}| = 2\pi/d`. Thus: - **Q = 0**: Uniform object (no spatial variation) - **|Q| large**: Fine features (small :math:`d`, rapid phase changes) - **|Q| small**: Coarse features (large :math:`d`, slow phase changes) - **Direction of Q**: Orientation of the periodic structure Think of **Q** like musical notes: high **Q** corresponds to high pitch (fast oscillations in space), low **Q** corresponds to low pitch (slow oscillations). .. note:: **What "FT over scan positions" means**: For each detector pixel :math:`\mathbf{k}`, extract how its intensity varies across all scan positions :math:`\mathbf{R}` (a 2D map), then FFT that map. You're transforming **scan patterns**, not diffraction patterns. Example: 100×100 scan with 256×256 detector → FFT 256×256 different 100×100 arrays, one per detector pixel. Applying this Fourier transform to each of the four terms: 1. :math:`|P(\mathbf{k})|^2` has no :math:`\mathbf{R}` dependence → becomes :math:`|P(\mathbf{k})|^2 \delta(\mathbf{Q})` at :math:`\mathbf{Q} = 0` 2. :math:`|\tilde{O} * P|^2` is weak (second-order) → small contribution 3. Cross term with :math:`e^{+i\mathbf{k} \cdot \mathbf{R}}` → shifts to :math:`\mathbf{Q} = +\mathbf{k}` by Fourier shift theorem 4. Cross term with :math:`e^{-i\mathbf{k} \cdot \mathbf{R}}` → shifts to :math:`\mathbf{Q} = -\mathbf{k}` by Fourier shift theorem The result is **three overlapping disks** in the 4D :math:`(\mathbf{k}, \mathbf{Q})` space. The delta functions below are idealizations (in reality the sidebands have finite width): .. math:: G(\mathbf{k}, \mathbf{Q}) = |P(\mathbf{k})|^2 \delta(\mathbf{Q}) + P(\mathbf{k}) [\tilde{O}^*(\mathbf{k}) * P^*(\mathbf{k})] \delta(\mathbf{Q} - \mathbf{k}) + P^*(\mathbf{k}) [\tilde{O}(\mathbf{k}) * P(\mathbf{k})] \delta(\mathbf{Q} + \mathbf{k}) The three terms represent distinct regions in :math:`(\mathbf{k}, \mathbf{Q})` space: - **Center disk** at :math:`\mathbf{Q} = 0`: Contains :math:`|P|^2` (probe autocorrelation), purely real, no phase information - **Right sideband** at :math:`\mathbf{Q} = +\mathbf{k}`: Contains :math:`P^*(\mathbf{k})[\tilde{O}(\mathbf{k}) * P(\mathbf{k})]` - typically selected for SSB - **Left sideband** at :math:`\mathbf{Q} = -\mathbf{k}`: Contains :math:`P(\mathbf{k})[\tilde{O}^*(\mathbf{k}) * P^*(\mathbf{k})]` - complex conjugate .. important:: **Why sidebands contain phase information (KEY CONCEPT):** The center disk :math:`|P|^2` is **purely real** - magnitude squared destroys phase. But the sidebands contain the **complex object-probe product** :math:`\tilde{O}(\mathbf{k}) * P(\mathbf{k})`, which preserves both amplitude AND phase. **The key steps to extract phase**: 1. **Sideband contains**: :math:`P^*(\mathbf{k})[\tilde{O}(\mathbf{k}) * P(\mathbf{k})] \approx |P(\mathbf{k})|^2 \tilde{O}(\mathbf{k})` 2. **Divide by probe**: :math:`\frac{|P(\mathbf{k})|^2 \tilde{O}(\mathbf{k})}{|P(\mathbf{k})|^2} = \tilde{O}(\mathbf{k})` - the **complex** Fourier transform of your object 3. **Inverse FT**: :math:`o(x,y) = \mathcal{F}^{-1}[\tilde{O}(\mathbf{k})] \approx i\phi(x,y)` - complex value encoding phase 4. **Extract phase**: :math:`\phi(x,y) = \text{Im}[o(x,y)]` - the imaginary part is your phase map! **Why this works**: Under weak-phase approximation :math:`O(x,y) \approx 1 + i\phi(x,y)`, the sideband isolates the :math:`i\phi` term (the DC "1" stayed in the center disk). The division by :math:`|P|^2` removes probe effects, leaving pure object information with phase intact. .. note:: **Which sideband to select?** Either sideband works, but **only one** should be selected. The right sideband (at :math:`\mathbf{Q} = +\mathbf{k}`) and left sideband (at :math:`\mathbf{Q} = -\mathbf{k}`) contain complex-conjugate information. Selecting the right sideband gives :math:`\tilde{O}(\mathbf{k})`; selecting the left gives :math:`\tilde{O}^*(\mathbf{k})`. After taking the imaginary part in real space, both yield the same phase :math:`\phi(x,y)` up to a sign. By convention, the right sideband is typically selected. Note: :math:`G(\mathbf{k}, \mathbf{Q})` is a 4D function (2D detector + 2D spatial frequency). For visualization, we typically show 2D slices or projections. **Understanding the overlap condition** The condition :math:`|\mathbf{Q}| < 2k_{\text{max}}` (equivalently :math:`\lambda |\mathbf{Q}| < 2\alpha`, where :math:`\alpha` is convergence angle) determines when the disks overlap. Here: - :math:`\mathbf{Q}` represents object spatial frequency (how fast :math:`\phi(x,y)` varies) - :math:`k_{\text{max}}` is the probe aperture radius When object features are coarse (:math:`|\mathbf{Q}|` small), the sidebands at :math:`\pm \mathbf{k}` remain close to center, creating significant overlap. The **double overlap regions** between center and sidebands contain interference terms preserving phase. The **triple overlap** has highest information but is contaminated. By selecting one sideband with an aperture (preferably where center disk overlap is minimal), object information is isolated for direct reconstruction. The diagram below shows how the three regions separate in reciprocal space: .. plot:: :include-source: False import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle fig, ax = plt.subplots(figsize=(10, 6)) # Create three circles with significant overlap to show triple overlap region disk_radius = 1.0 # All disks same size separation = 0.8 # Very close spacing for significant triple overlap center_circle = Circle((0, 0), disk_radius, color='lightgray', alpha=0.5, label='0-order disk') left_circle = Circle((-separation, 0), disk_radius, color='#FF6B6B', alpha=0.6, label=r'$-\Delta q$ sideband') right_circle = Circle((separation, 0), disk_radius, color='#4A90E2', alpha=0.6, label=r'$+\Delta q$ sideband') ax.add_patch(center_circle) ax.add_patch(left_circle) ax.add_patch(right_circle) # Add labels ax.text(-separation, 0, r'$\tilde{O}(k_x, k_y) * P(k_x, k_y)$' + '\nTransfer', ha='center', va='center', fontsize=9, weight='bold', color='darkred') ax.text(separation, 0, r'$\tilde{O}(k_x, k_y) * P(k_x, k_y)$' + '\nTransfer', ha='center', va='center', fontsize=9, weight='bold', color='darkblue') # Mark overlap regions ax.text(-0.5, 0.8, 'double overlap', ha='center', fontsize=9, style='italic', color='gray') ax.text(0.5, 0.8, 'double overlap', ha='center', fontsize=9, style='italic', color='gray') ax.text(0, 1.05, 'triple overlap', ha='center', fontsize=9, style='italic', color='gray') # Add arrows showing separation arrow_y = -1.5 ax.annotate('', xy=(-separation, arrow_y), xytext=(0, arrow_y), arrowprops=dict(arrowstyle='<->', color='black', lw=1.5)) ax.text(-separation/2, arrow_y - 0.25, r'$\Delta q$', ha='center', fontsize=11) ax.annotate('', xy=(separation, arrow_y), xytext=(0, arrow_y), arrowprops=dict(arrowstyle='<->', color='black', lw=1.5)) ax.text(separation/2, arrow_y - 0.25, r'$\Delta q$', ha='center', fontsize=11) # Axis settings ax.set_xlim(-3.5, 3.5) ax.set_ylim(-2.2, 2.2) ax.set_aspect('equal') ax.axis('off') ax.set_title('Spatial separation of sidebands in reciprocal space', fontsize=13, weight='bold', pad=20) # Add legend ax.legend(loc='upper right', fontsize=10, framealpha=0.9) plt.tight_layout() .. important:: **This is NOT a diffraction pattern!** Individual diffraction patterns recorded at each scan position do NOT show these three disks. At one scan position, only the interference pattern :math:`I(\mathbf{k}; \mathbf{R})` is visible - a single 2D intensity distribution on the detector. These diagrams show the structure of :math:`G(\mathbf{k}, \mathbf{Q})` - a mathematical object created by Fourier transforming over **all scan positions**. Since :math:`G(\mathbf{k}, \mathbf{Q})` is 4D, we visualize it by showing 2D slices or projections. The "three disks" represent different components that emerge in this transformed space: the center disk (autocorrelation at :math:`\mathbf{Q} = 0`), and the two sidebands at :math:`\mathbf{Q} = \pm \mathbf{k}` (containing object information). The sidebands exist only after processing the entire 4D dataset - they are not visible in any single measurement. .. _phase_extraction: Phase extraction from single sideband ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Why only one sideband is needed** The two sidebands (at :math:`\mathbf{Q} = +\mathbf{k}` and :math:`\mathbf{Q} = -\mathbf{k}`) contain the same information, just complex-conjugated. Keeping both causes interference and introduces redundancy. Selecting only one sideband produces a clean, linear signal. This is the reason for the term "single sideband." **The complete SSB workflow** What is actually detected is the total intensity pattern—the sidebands at :math:`\pm \Delta q` are not directly visible as separate entities. They exist mathematically as components of the Fourier decomposition of the recorded intensity, mixed together with the autocorrelation terms. The spatial separation discussed earlier refers to how these terms separate in reciprocal space after Fourier analysis. **SSB reconstructs the entire field of view simultaneously** Because SSB requires Fourier transforming over all scan positions, it reconstructs the **entire field of view** simultaneously, not individual spots. SSB cannot be performed on a single diffraction pattern—the full 4D dataset is required. The phase map :math:`\phi(x,y)` that emerges covers the complete scanned area. **Step 1: Identify minimal sideband overlap** When the probe scan step size is large enough, sidebands separate cleanly with minimal overlap: .. plot:: :include-source: False import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle fig, ax = plt.subplots(figsize=(10, 6)) # Create three circles with minimal overlap disk_radius = 1.0 separation = 1.5 # Closer separation to show minimal overlap center_circle = Circle((0, 0), disk_radius, color='lightgray', alpha=0.5, label='0-order disk') left_circle = Circle((-separation, 0), disk_radius, color='#FF6B6B', alpha=0.6, label=r'$-\Delta q$ sideband') right_circle = Circle((separation, 0), disk_radius, color='#4A90E2', alpha=0.6, label=r'$+\Delta q$ sideband') ax.add_patch(center_circle) ax.add_patch(left_circle) ax.add_patch(right_circle) # Add labels ax.text(-separation, 0, r'$\tilde{O}(k_x, k_y) * P(k_x, k_y)$' + '\nTransfer', ha='center', va='center', fontsize=9, weight='bold', color='darkred') ax.text(separation, 0, r'$\tilde{O}(k_x, k_y) * P(k_x, k_y)$' + '\nTransfer', ha='center', va='center', fontsize=9, weight='bold', color='darkblue') # Add arrows showing separation arrow_y = -1.5 ax.annotate('', xy=(-separation, arrow_y), xytext=(0, arrow_y), arrowprops=dict(arrowstyle='<->', color='black', lw=1.5)) ax.text(-separation/2, arrow_y - 0.25, r'$\Delta q$', ha='center', fontsize=11) ax.annotate('', xy=(separation, arrow_y), xytext=(0, arrow_y), arrowprops=dict(arrowstyle='<->', color='black', lw=1.5)) ax.text(separation/2, arrow_y - 0.25, r'$\Delta q$', ha='center', fontsize=11) # Draw aperture around one sideband aperture = Circle((separation, 0), disk_radius * 1.3, fill=False, edgecolor='green', linewidth=3, linestyle='--', label='Aperture selection') ax.add_patch(aperture) # Axis settings ax.set_xlim(-4.5, 4.5) ax.set_ylim(-2.5, 2.5) ax.set_aspect('equal') ax.axis('off') ax.set_title('Minimal sideband overlap - ideal for clean phase extraction', fontsize=13, weight='bold', pad=20) # Add legend ax.legend(loc='upper right', fontsize=10, framealpha=0.9) plt.tight_layout() This shows the :math:`G(\mathbf{k}, \mathbf{Q})` space with minimal sideband overlap (achieved when scan step size is large enough). A circular aperture (shown in green) is drawn around one sideband. The aperture size trades off resolution (smaller aperture, lower resolution) versus signal-to-noise ratio (larger aperture, more signal). This aperture selection happens in the transformed :math:`(\mathbf{k}, \mathbf{Q})` space after Fourier transforming over all scan positions, not on individual diffraction patterns. **Step 2: Crop and shift sideband to center** .. plot:: :include-source: False import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import Circle, FancyArrowPatch import matplotlib.patches as mpatches fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6)) # Left plot: Show aperture selection disk_radius = 1.0 separation = 1.5 # Left panel - before cropping center_circle1 = Circle((0, 0), disk_radius, color='lightgray', alpha=0.5) left_circle1 = Circle((-separation, 0), disk_radius, color='#FF6B6B', alpha=0.6) right_circle1 = Circle((separation, 0), disk_radius, color='#4A90E2', alpha=0.6) ax1.add_patch(center_circle1) ax1.add_patch(left_circle1) ax1.add_patch(right_circle1) # Draw aperture around right sideband aperture = Circle((separation, 0), disk_radius * 1.2, fill=False, edgecolor='green', linewidth=3, linestyle='--') ax1.add_patch(aperture) ax1.text(separation, 0, r'$\tilde{O} * P$', ha='center', va='center', fontsize=11, weight='bold', color='darkblue') ax1.set_xlim(-3.5, 3.5) ax1.set_ylim(-2.5, 2.5) ax1.set_aspect('equal') ax1.axis('off') ax1.set_title('Step 1: Select one sideband', fontsize=13, weight='bold', pad=20) # Right panel - after cropping and shifting cropped_circle = Circle((0, 0), disk_radius, color='#4A90E2', alpha=0.7) ax2.add_patch(cropped_circle) ax2.text(0, 0, r'$\tilde{O} * P$' + '\n(centered)', ha='center', va='center', fontsize=11, weight='bold', color='darkblue') # Add coordinate axes ax2.arrow(-1.8, 0, 3.3, 0, head_width=0.15, head_length=0.15, fc='black', ec='black', alpha=0.3) ax2.arrow(0, -1.8, 0, 3.3, head_width=0.15, head_length=0.15, fc='black', ec='black', alpha=0.3) ax2.text(1.7, -0.3, r'$k_x$', fontsize=11, alpha=0.5) ax2.text(0.2, 1.7, r'$k_y$', fontsize=11, alpha=0.5) ax2.set_xlim(-2.5, 2.5) ax2.set_ylim(-2.5, 2.5) ax2.set_aspect('equal') ax2.axis('off') ax2.set_title('Step 2: Shift to center for reconstruction', fontsize=13, weight='bold', pad=20) # Add arrow between panels fig.text(0.48, 0.5, '→', fontsize=40, ha='center', va='center', color='green', weight='bold') plt.tight_layout() These diagrams show operations in :math:`G(\mathbf{k}, \mathbf{Q})` space (after Fourier transforming over all scan positions). Left: aperture selection isolates one sideband. Right: the selected sideband shifted to center for reconstruction. After Fourier transforming over scan positions, the sideband at :math:`\mathbf{Q} = +\mathbf{k}` contains the object-probe convolution. An aperture :math:`A(\mathbf{Q})` in :math:`(\mathbf{k}, \mathbf{Q})` space isolates this sideband: .. math:: S(\mathbf{k}, \mathbf{Q}) = A(\mathbf{Q} - \mathbf{k}) \cdot G(\mathbf{k}, \mathbf{Q}) With sufficient separation (minimal overlap case), the aperture primarily captures: .. math:: S(\mathbf{k}, \mathbf{Q}) \approx P^*(\mathbf{k}) [\tilde{O}(\mathbf{k}) * P(\mathbf{k})] \delta(\mathbf{Q} - \mathbf{k}) Shifting the sideband to center (setting :math:`\mathbf{Q}' = \mathbf{Q} - \mathbf{k}` for the right sideband at :math:`\mathbf{Q} = +\mathbf{k}`): .. math:: S_{\text{centered}}(\mathbf{k}) \approx P^*(\mathbf{k}) [\tilde{O}(\mathbf{k}) * P(\mathbf{k})] \approx |P(\mathbf{k})|^2 \tilde{O}(\mathbf{k}) .. warning:: **Critical approximation**: The second :math:`\approx` above is a **strong assumption** that requires: - **Localized probe in k-space**: The object's spectrum :math:`\tilde{O}(\mathbf{k})` varies slowly over the support of the probe kernel, so inside the aperture :math:`\tilde{O}(\mathbf{k})` is almost constant on the scale of the probe kernel - **Convolution behaves like multiplication**: :math:`\tilde{O}(\mathbf{k}) * P(\mathbf{k}) \approx \tilde{O}(\mathbf{k}) |P(\mathbf{k})|^2` This is where SSB can break. When NOT true (sharp object features, strong phase variations, thick samples), the approximation fails and SSB produces artifacts, contrast errors, or frequency-dependent distortion. Your derivation is correct **given this extra assumption**. Treat this in your head as "this is where SSB can break." **Step 3: Divide by the probe** Now centered at the origin, the sideband contains approximately :math:`|P(\mathbf{k})|^2 \tilde{O}(\mathbf{k})`. Dividing by the probe intensity isolates the object: .. math:: \frac{|P(\mathbf{k})|^2 \tilde{O}(\mathbf{k})}{|P(\mathbf{k})|^2} = \tilde{O}(\mathbf{k}) The probe :math:`P(\mathbf{k})` must be known (measured from vacuum or constructed analytically). This division recovers the object's Fourier content without iteration. **Step 4: Inverse Fourier transform back to real space** After dividing by the probe, we have :math:`\tilde{O}(\mathbf{k})` in Fourier space. But remember, we're still in the mixed :math:`(\mathbf{k}, \mathbf{Q})` representation. We need to **inverse Fourier transform back over** :math:`\mathbf{Q}` to return to real space coordinates :math:`\mathbf{r} = (x, y)`: .. math:: o(\mathbf{r}) = \int \tilde{O}(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{r}} d\mathbf{k} This recovers the **phase component** of the object transmission function in real space. **What do we actually get?** Under the weak-phase approximation :math:`O(x,y) \approx 1 + i\phi(x,y)`, the sideband isolation removed the constant term (which stayed in the center disk at :math:`\mathbf{Q} = 0`). The reconstructed quantity :math:`o(x,y)` contains only the phase-dependent part: .. math:: o(x,y) \approx i\phi(x,y) Taking the imaginary part gives the phase directly: .. math:: \phi(x,y) = \text{Im}[o(x,y)] This is the object's phase map. SSB directly recovers :math:`\phi(x,y)` without iteration or integration. The weak-phase approximation makes this simple division-based reconstruction possible. .. note:: **SSB vs DPC**: Do not confuse SSB with differential phase contrast (DPC). DPC measures deflection angles and requires integration to get phase. SSB directly extracts phase via sideband isolation and division. They are fundamentally different techniques, though both use 4D-STEM data. .. note:: **Caveats on recovered phase**: - **DC content lost**: The constant term (overall phase offset and slow background) is lost when the center disk is discarded. Absolute offset and very low-frequency behavior are not well constrained, similar to DPC. - **Regularization in practice**: Practitioners often apply regularization or filtering in frequency space to improve signal-to-noise ratio and stability. **Summary: Why SSB is fast** All four SSB steps are linear operations: 1. **Fourier transform over scan positions** :math:`\mathbf{R} \to \mathbf{Q}` (separates sidebands in :math:`G(\mathbf{k}, \mathbf{Q})`) 2. **Aperture selection** in :math:`(\mathbf{k}, \mathbf{Q})` space (isolates one sideband, suppresses conjugate) 3. **Divide by probe** :math:`|P(\mathbf{k})|^2` (removes probe modulation to get :math:`\tilde{O}(\mathbf{k})`) 4. **Inverse Fourier transform** :math:`\mathbf{k} \to \mathbf{r}` (get :math:`o(\mathbf{r})`, extract :math:`\phi = \text{Im}[o]`) No iterative solver. No nonlinear optimization. No thousands of ePIE iterations. Just linear algebra. This simplicity makes SSB fast and robust to initialization errors. .. _probe_specification: Probe specification and aberration measurement ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **What probe information does SSB need?** SSB reconstruction requires knowing the probe function :math:`P(k_x, k_y)` to perform the division step in reciprocal space. The probe function consists of two parts: .. math:: P(k_x, k_y) = A(k_x, k_y) e^{i\chi(k_x, k_y)} where :math:`A(k_x, k_y)` is the aperture modulus (amplitude, typically a sharp circular cutoff at :math:`k_{\text{max}}`) and :math:`\chi(k_x, k_y)` is the aberration phase function encoding optical aberrations. **Simplified probe assumption for highly corrected microscopes** For modern aberration-corrected microscopes (Cs-corrected systems), most higher-order aberrations are negligible. In this case, a highly simplified probe model can be used: .. math:: \chi(k_x, k_y) \approx \pi C_1 \lambda k^2 where: - :math:`C_1` is the defocus (typically the only significant aberration) - :math:`\lambda` is the electron wavelength - :math:`k = \sqrt{k_x^2 + k_y^2}` is the radial coordinate in reciprocal space SSB requires accurate knowledge of the convergence angle to properly define :math:`k_{\text{max}}`. Fortunately, the convergence angle is usually stable and easily obtained from the microscope settings. **General shape of the probe aperture** A circular, uniform aperture is usually a good approximation: .. math:: A(k) = \begin{cases} 1 & |k| < k_{\text{max}} \\ 0 & \text{else} \end{cases} Unless your condenser system is misaligned or your aperture is damaged, this is safe. **Near-zero defocus** For near-zero defocus conditions: .. math:: C_1 \approx 0 This is one of the most important parts. SSB is extremely sensitive to defocus error, but if yours really is near zero, guessing "0" is fine. **Small higher-order aberrations** If your microscope is well-corrected (Cs-corrector), then C3, C5... are small enough that: .. math:: \chi(k) \approx 0 \quad \text{as a first guess} This is not perfect, but usually good enough to do SSB without horrible artifacts. **When more detailed probe measurement is needed** For uncorrected microscopes or when higher accuracy is required, the full aberration phase :math:`\chi(k_x, k_y)` must be measured. This includes defocus, spherical aberration, astigmatism, coma, and higher-order terms. **Why aberrations appear in the sideband** Because the selected sideband :math:`S(\mathbf{k})` depends directly on the probe: .. math:: S(\mathbf{k}) \approx |P(\mathbf{k})|^2 \tilde{O}(\mathbf{k}) it includes both the probe amplitude :math:`|P(\mathbf{k})|` and the aberration phase :math:`\chi(\mathbf{k})` (contained in :math:`P = A e^{i\chi}`). The phase and shape of :math:`S(\mathbf{k})` directly encode :math:`\chi(\mathbf{k})`. **Breaking the chicken-and-egg problem** SSB faces a fundamental issue: to extract phase, the probe must be known, but to measure the probe, the object must be known. SSB breaks this cycle using **two separate measurements**: 1. **First, measure the probe** on a simple reference (vacuum, amorphous film) where the object is known 2. **Then, use that measured probe** to extract phase from unknown samples For a reference with known :math:`\tilde{O}(k_x, k_y)`, the probe is recovered from the sideband: .. math:: P(k_x, k_y) = \frac{S(k_x, k_y)}{\tilde{O}(k_x, k_y)} Once measured, this probe can be used for all subsequent reconstructions in the same imaging session (assuming the probe doesn't change). **Required conditions for accurate measurement** SSB provides clean probe aberration measurements under specific conditions: - **Simple object**: Vacuum edge, uniform amorphous carbon, or ultrathin film allow accurate estimation of :math:`\tilde{O}(k_x, k_y)` for division - **High signal-to-noise ratio**: Phase extraction requires low noise to resolve fine details of :math:`\chi(k_x, k_y)` - **Well-separated sidebands**: The sideband must be far from the center disk and noise to avoid contamination from other terms - **Coherent probe**: Partially coherent or incoherent probes broaden the sideband and degrade phase resolution - **Weak-phase approximation valid**: Thick samples or strong-scattering materials violate the approximation and introduce errors - **Negligible multiple scattering**: Thick or crystalline samples scatter multiple times, corrupting the sideband When all these conditions hold, division by the object term yields a clean :math:`P(k_x, k_y)`, and the aberrations can be reliably extracted. Parameters that cannot be reliably estimated ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Even with a highly corrected microscope, certain parameters introduce artifacts if guessed incorrectly. These are the danger areas: **Residual astigmatism** Even tiny astigmatism distorts the sideband shape and can introduce anisotropic artifacts. This cannot be blindly guessed. **Residual coma** Coma shifts the probe and skews the sideband. If small, SSB survives; if moderate, it breaks. **Finite source size and partial coherence** If the illumination is not coherent enough, the sideband contrast decreases and blurs. Guessing wrong means SSB incorrectly amplifies noise on division. **Defocus drift** If the defocus changes during the scan due to stage drift or thermal instability, the probe varies across the field of view. SSB assumes a constant probe, so this introduces systematic errors that cannot be corrected by a single probe measurement. .. _practical_guide: Practical guide for using SSB ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Required microscope parameters** SSB reconstruction requires these experimental parameters: - **Convergence angle** (α): Defines :math:`k_{\text{max}}` for the probe aperture (units: mrad) - **Acceleration voltage**: Determines the electron wavelength :math:`\lambda` (units: keV → Å) - **Camera length**: Converts detector pixels to reciprocal space units (mrad) These are standard microscope parameters that should be stable and easy to obtain from your instrument. .. note:: **Unit conversions**: - Convergence angle in mrad → :math:`k_{\text{max}}` in Ų: :math:`k_{\text{max}} = \alpha / \lambda` where :math:`\alpha` is in radians - Acceleration voltage in keV → wavelength in Å: :math:`\lambda = \frac{12.398}{\sqrt{E(E + 1022)}}` where :math:`E` is in keV - Example: 60 keV → λ = 0.0487 Å, 34 mrad → k_max = 0.698 Ų **Method 1: Use a vacuum region inside your scan** This is the easiest and most robust method. If your scan area includes vacuum (no sample), the SSB sideband over vacuum directly measures the probe: .. math:: S_{\text{vacuum}}(k_x, k_y) \approx P(k_x, k_y) because over vacuum, :math:`O(x,y) = 1` (no object), so :math:`\tilde{O} = \delta(k)` and the sideband equals the probe. **Procedure:** 1. Collect 4D-STEM data with part of the scan over vacuum (sample edge or hole) 2. Process the vacuum region separately to extract :math:`P(k_x, k_y)` 3. Use this measured probe for all sample positions in the same dataset This provides a real probe measurement without estimation. The vacuum measurement captures all aberrations, including residual astigmatism, coma, and defocus, as they actually exist at the time of data acquisition. **Method 2: Thin amorphous reference film** If no vacuum is available, use a uniform thin amorphous carbon or silicon nitride film where :math:`\phi(x,y) \approx` constant. The sideband over this region approximates the probe with only a weak, slowly varying phase offset. **Method 3: Simple analytical probe (for quick tests)** For rapid prototyping or when accuracy is not critical, construct an analytical probe from microscope parameters: .. math:: P(k_x, k_y) = A(k) e^{i\chi(k)} where: - :math:`A(k) = 1` for :math:`k < k_{\text{max}}`, 0 otherwise (from convergence angle) - :math:`\chi(k) = \pi C_1 \lambda k^2` (defocus only, if near zero set :math:`C_1 = 0`) This works for highly corrected microscopes but ignores residual aberrations. **Processing workflow** 1. **Load 4D-STEM dataset**: :math:`I(\mathbf{k}; \mathbf{R})` where :math:`\mathbf{R}` are scan positions, :math:`\mathbf{k}` are detector pixels 2. **Measure or construct probe** :math:`P(\mathbf{k})` using vacuum region, reference film, or analytical model 3. **Fourier transform over scan positions**: Transform the 4D data :math:`I(\mathbf{k}; \mathbf{R})` over the :math:`\mathbf{R}` dimensions to create :math:`G(\mathbf{k}, \mathbf{Q})`, where :math:`\mathbf{Q}` is spatial frequency 4. **Select one sideband**: Apply circular aperture :math:`A(\mathbf{Q} - \mathbf{k})` in :math:`(\mathbf{k}, \mathbf{Q})` space to isolate sideband at :math:`\mathbf{Q} = \pm \mathbf{k}` 5. **Shift sideband to center**: Re-center by setting :math:`\mathbf{Q}' = \mathbf{Q} \mp \mathbf{k}` to get :math:`S(\mathbf{k})` 6. **Divide by probe intensity**: :math:`\tilde{O}(\mathbf{k}) = S(\mathbf{k}) / |P(\mathbf{k})|^2` 7. **Inverse Fourier transform**: Transform :math:`\tilde{O}(\mathbf{k})` back to real space to get :math:`o(x,y) \approx i\phi(x,y)` 8. **Extract phase**: :math:`\phi(x,y) = \text{Im}[o(x,y)]` The vacuum region method (Method 1) for measuring the probe is strongly recommended because it captures actual experimental conditions without guesswork. Advantages and limitations ^^^^^^^^^^^^^^^^^^^^^^^^^^^ **Advantages:** - Direct, non-iterative phase extraction via linear operations - Simple workflow once probe is known - Works well for thin samples under weak-phase approximation - Easy to measure probe aberrations from reference samples - Entire field of view reconstructed simultaneously **Limitations:** - Requires weak-phase approximation (thin samples, :math:`\phi \ll 1`) - Resolution limited by sideband separation and signal-to-noise ratio - Probe knowledge critical—errors in aberration coefficients cause artifacts - Lower redundancy exploitation compared to iterative ptychography - Phase contrast transfer function falls off at high spatial frequencies Implementation in quantem ^^^^^^^^^^^^^^^^^^^^^^^^^^ The following sections describe how the conceptual framework above is implemented in the quantem library. The implementation follows the same physics but uses computational notation and provides practical tools for real data analysis. Probe representation ~~~~~~~~~~~~~~~~~~~~ The probe function :math:`P(\mathbf{k})` from :ref:`probe_specification` is implemented as: .. math:: P(\mathbf{k}) = A(\mathbf{k}) \cdot \exp\left[-i\chi(\mathbf{k})\right] where :math:`A(\mathbf{k})` is the aperture function (soft or hard-edged circular) and :math:`\chi(\mathbf{k})` is the aberration phase function. Aberration function expansion ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The aberration function :math:`\chi(\mathbf{k})` is parameterized using Zernike-like polynomial expansion: .. math:: \chi(\mathbf{k}) = \frac{2\pi}{\lambda} \sum_{n,m} \frac{1}{n+1} C_{nm} k^{n+1} \cos[m(\phi - \phi_{nm})] where: - :math:`\lambda` is the electron wavelength - :math:`C_{nm}` are the aberration coefficients (e.g., :math:`C_{10}` = defocus, :math:`C_{12}` = astigmatism, :math:`C_{30}` = spherical aberration) - :math:`k = |\mathbf{k}|` and :math:`\phi = \arctan(k_y/k_x)` are polar coordinates - :math:`\phi_{nm}` are azimuthal angles for non-rotationally symmetric aberrations The most important aberration terms are: .. math:: \chi(\mathbf{k}) \approx \pi C_{10} \lambda k^2 + \frac{\pi}{2} C_{12} \lambda k^2 \cos[2(\phi - \phi_{12})] + \frac{\pi}{2} C_{30} \lambda^3 k^4 + \ldots where :math:`C_{10}` is defocus, :math:`C_{12}` is two-fold astigmatism with angle :math:`\phi_{12}`, and :math:`C_{30}` is third-order spherical aberration. These are the same aberration coefficients discussed in the "Probe specification" section earlier. SSB transfer function ^^^^^^^^^^^^^^^^^^^^^^ The implementation uses a transfer function :math:`\Gamma(\mathbf{Q}, \mathbf{k})` (sometimes called the deconvolution kernel) to perform sideband isolation and probe division. This implements the conceptual steps from :ref:`phase_extraction` (aperture selection, shift, divide) in a single mathematical operation. Three deconvolution kernel options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The ``deconvolution_kernel`` parameter selects different reconstruction methods: **1. 'none' (Bright-field imaging)** No deconvolution applied. Simply returns the measured bright-field intensity without phase extraction. Fast but provides no phase information. **2. 'quadratic' (Parallax method)** .. math:: \Gamma(\mathbf{Q}, \mathbf{k}) \approx \mathbf{Q} \cdot \mathbf{k} A simplified approximation that's faster to compute but less accurate than full SSB. Useful for hyperparameter optimization where speed matters. **3. 'full' (Complete SSB - default)** .. math:: \Gamma(\mathbf{Q}, \mathbf{k}) = P(\mathbf{Q} - \mathbf{k}) P^*(\mathbf{k}) - P^*(\mathbf{Q} + \mathbf{k}) P(\mathbf{k}) This implements the complete SSB method from :ref:`phase_extraction`: - First term: selects the :math:`+\mathbf{k}` sideband - Subtraction: suppresses the :math:`-\mathbf{k}` conjugate sideband (removes twin-image artifacts) - The form includes implicit probe division The asymmetric form (subtraction) gives better noise properties than the symmetric form (addition). Normalization ~~~~~~~~~~~~~ To prevent numerical instabilities from division by near-zero values, the operator is normalized: .. math:: \Gamma_{\text{norm}}(\mathbf{Q}, \mathbf{k}) = \frac{\Gamma^*(\mathbf{Q}, \mathbf{k})}{|\Gamma(\mathbf{Q}, \mathbf{k})| + \epsilon} where :math:`\epsilon = 10^{-8}` is a regularization parameter that prevents division by zero when the probe is weak. Reconstruction workflow ^^^^^^^^^^^^^^^^^^^^^^^^ The implementation follows the conceptual workflow from :ref:`phase_extraction`: Data preparation ~~~~~~~~~~~~~~~~ 1. Load 4D-STEM dataset :math:`I(\mathbf{k}; \mathbf{R})` with shape ``(scan_x, scan_y, detector_x, detector_y)`` 2. Identify bright-field disk using center-of-mass or fit_probe_circle 3. Extract BF intensities within aperture mask Probe specification ~~~~~~~~~~~~~~~~~~~ Define the probe :math:`P(\mathbf{k})` using one of the three methods from :ref:`practical_guide`: .. dropdown:: Code example: Probe specification :icon: info :color: info .. code-block:: python # Method 1: From vacuum region (recommended) # Method 2: From thin amorphous film # Method 3: Analytical from microscope parameters cmplx_probe = evaluate_probe( alpha, phi, # Polar coordinates in k-space semiangle_cutoff, # Convergence angle [mrad] wavelength, # Electron wavelength [Å] aberration_coefs={ # Aberration coefficients [Å] 'C10': 0, # Defocus (near zero) 'C12': 10, # Astigmatism magnitude 'phi12': 0.5, # Astigmatism angle [rad] } ) Applying the transfer function ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For each scan position, apply the transfer function :math:`\Gamma(\mathbf{Q}, \mathbf{k})` which implements the steps from :ref:`phase_extraction`: 1. **Fourier transform over scan positions**: Transform the 4D dataset :math:`I(\mathbf{k}; \mathbf{R})` over the scan coordinate :math:`\mathbf{R}` to create :math:`G(\mathbf{k}, \mathbf{Q})` from :ref:`sideband_formation`: .. math:: G(\mathbf{k}, \mathbf{Q}) = \mathcal{F}_{\mathbf{R} \to \mathbf{Q}}[I(\mathbf{k}; \mathbf{R})] In the implementation, this is ``torch.fft.fft2(vbf_stack)`` where ``vbf_stack`` is the bright-field intensity at each scan position. 2. **Apply** :math:`\Gamma(\mathbf{Q}, \mathbf{k})` to isolate and deconvolve one sideband: .. math:: \tilde{O}(\mathbf{Q}) = G(\mathbf{k}, \mathbf{Q}) \cdot \Gamma(\mathbf{Q}, \mathbf{k}) This multiplication in Fourier space implements the aperture selection, shift, and probe division. 3. **Inverse Fourier transform** to get object :math:`o(\mathbf{r})`: .. math:: o(\mathbf{r}) = \mathcal{F}_{\mathbf{Q} \to \mathbf{r}}^{-1}[\tilde{O}(\mathbf{Q})] In code: ``torch.fft.ifft2(fourier_factor)`` 4. **Extract phase**: :math:`\phi(\mathbf{r}) = \text{Im}[o(\mathbf{r})]` or equivalently ``corrected_stack.angle()`` The final result is averaged over all positions for improved signal-to-noise. Scan-detector rotation correction ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ As mentioned in the practical guide, there can be a rotation :math:`\theta` between the scan coordinate system (probe motion) and the detector coordinate system (diffraction pattern recording). The implementation corrects for this: .. math:: \begin{pmatrix} k_x' \\ k_y' \end{pmatrix} = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} k_x \\ k_y \end{pmatrix} where :math:`\theta` is in **radians** for the mathematical operations. .. note:: **Angle units**: Internal calculations use **radians**, but the user-facing API accepts **degrees** for convenience. The function parameters are: - ``rotation_angle_deg``: User input in degrees (e.g., 0.38°) - ``rotation_angle_rad``: Internal computation in radians (e.g., 0.0066 rad) Conversion: :math:`\theta_{\text{rad}} = \theta_{\text{deg}} \times \pi / 180` This rotation is applied to the reciprocal space grid before probe evaluation, ensuring the reconstructed object aligns with the scan coordinate system. The rotation angle can be optimized (see next section). Optimizing probe parameters ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Motivation As discussed in :ref:`probe_specification`, accurate probe knowledge is critical. Without a vacuum region or reference sample (see :ref:`practical_guide`), the probe parameters (aberration coefficients and rotation angle) can be optimized directly from the data. Self-consistency loss ~~~~~~~~~~~~~~~~~~~~~ The optimization uses a variance-based loss function: .. math:: \mathcal{L} = \frac{1}{N} \sum_{n=1}^{N} \|\phi_n(\mathbf{r}) - \bar{\phi}(\mathbf{r})\|^2 where: - :math:`\phi_n(\mathbf{r})`: Phase reconstructed from scan position :math:`n` alone - :math:`\bar{\phi}(\mathbf{r})`: Mean phase across all :math:`N` positions **Physical intuition**: If the probe parameters are correct, all scan positions should reconstruct the same object (since the sample doesn't move between measurements). The variance measures how much the individual reconstructions disagree—lower variance means better probe parameters. **How it works**: For each trial set of aberration coefficients and rotation angle: 1. Reconstruct phase :math:`\phi_n(\mathbf{r})` from each scan position independently 2. Compute the mean reconstruction :math:`\bar{\phi}(\mathbf{r})` across all positions 3. Calculate variance (disagreement) between individual and mean reconstructions 4. Try new parameters to minimize this variance **Important limitation**: This optimization requires **sufficient detector sampling**. With only 12×12 detector pixels, there's insufficient information in each diffraction pattern to constrain the parameters well. The optimization may converge but won't produce meaningful results. Reliable optimization typically requires ≥64×64 detector pixels, ideally 128×128 or more. Optimization algorithm ~~~~~~~~~~~~~~~~~~~~~~ The implementation uses **Optuna** with Tree-structured Parzen Estimator (TPE) sampling: 1. **Parameter Space**: For aberration-corrected microscopes (recommended starting ranges): .. math:: \mathbf{\Theta} = \{C_{10} \in [-100, 100]\text{ Å}, \, C_{12} \in [0, 50]\text{ Å}, \, \phi_{12} \in [-\pi/2, \pi/2], \, \theta \in [-10°, 10°]\} For uncorrected microscopes or initial exploration, use wider ranges: :math:`C_{10} \in [-400, 400]` Å, :math:`\theta \in [-180°, 180°]`, and increase :math:`C_{12}` range to [0, 100] Å. 2. **Objective Function**: .. code-block:: python def objective(trial): # Sample parameters C10 = trial.suggest_float('C10', -400, 400) C12 = trial.suggest_float('C12', 0, 100) phi12 = trial.suggest_float('phi12', -np.pi/2, np.pi/2) theta = trial.suggest_float('rotation_angle_rad', -np.pi, np.pi) # Reconstruct with these parameters reconstruct(aberration_coefs={'C10': C10, 'C12': C12, 'phi12': phi12}, rotation_angle_rad=theta) # Return variance loss return variance_loss() 3. **Optimization Loop**: Run for :math:`N_{\text{trials}}` iterations (typically 50-100) to find: .. math:: \mathbf{\Theta}^* = \arg\min_{\mathbf{\Theta}} \mathcal{L}(\mathbf{\Theta}) Code organization ^^^^^^^^^^^^^^^^^ Module structure ~~~~~~~~~~~~~~~~ The SSB algorithm is implemented across two main modules: **1. complex_probe.py** Contains core mathematical operators: - ``gamma_factor()``: Computes the SSB transfer function :math:`\Gamma(\mathbf{Q}, \mathbf{k})` - ``evaluate_probe()``: Evaluates complex probe with aberrations - ``aberration_surface()``: Computes aberration phase surface - ``aperture()``: Soft/hard aperture functions - ``spatial_frequencies()``: Generates Fourier grids with rotation **2. direct_ptychography.py** Orchestrates the reconstruction pipeline: - ``DirectPtychography`` class: Main interface - ``from_dataset4d()``: Constructor from 4D-STEM data - ``reconstruct()``: Performs SSB reconstruction - ``optimize_hyperparameters()``: Optuna-based optimization - ``_compute_gamma_operator()``: Wrapper for transfer function computation - ``variance_loss()``: Self-consistency metric Key implementation functions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Transfer function computation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. dropdown:: Code example: Transfer function computation :icon: info :color: info .. code-block:: python def gamma_factor( qmks: tuple[Tensor, Tensor], # (q - k) coordinates qpks: tuple[Tensor, Tensor], # (q + k) coordinates cmplx_probe_at_k: Tensor, # P(k) wavelength: float, semiangle_cutoff: float, soft_edges: bool, aberration_coefs: dict, angular_sampling: tuple, asymmetric_version: bool = True, normalize: bool = True, ) -> Tensor: """ Compute the SSB transfer function Γ(Q,k). Returns: Γ(q,k) = P(q-k)P*(k) - P*(q+k)P(k) [asymmetric] or Γ(q,k) = P(q-k)P*(k) + P*(q+k)P(k) [symmetric] """ # Evaluate probe at q-k probe_m = evaluate_probe(q_m * wavelength, phi_m, ...) # Evaluate probe at q+k probe_p = evaluate_probe(q_p * wavelength, phi_p, ...) # Compute transfer function Γ if asymmetric_version: gamma = probe_m * cmplx_probe_at_k.conj() - probe_p.conj() * cmplx_probe_at_k else: gamma = probe_m * cmplx_probe_at_k.conj() + probe_p.conj() * cmplx_probe_at_k # Normalize if normalize: gamma /= gamma.abs().clamp(min=1e-8) return gamma.conj() Probe evaluation ^^^^^^^^^^^^^^^^ .. dropdown:: Code example: Probe evaluation :icon: info :color: info .. code-block:: python def evaluate_probe( alpha: Tensor, # Radial angles [rad] phi: Tensor, # Azimuthal angles [rad] semiangle_cutoff: float, # Aperture size [mrad] angular_sampling: tuple, # Pixel size in angle space wavelength: float, # Wavelength [Å] aberration_coefs: dict, # Aberration coefficients ) -> Tensor: """ Evaluate complex probe function: P(k) = A(k) exp[-iχ(k)] """ # Aperture function (soft or hard edge) aperture = compute_aperture(alpha, semiangle_cutoff, ...) # Aberration phase chi = aberration_surface(alpha, phi, wavelength, aberration_coefs) # Complex probe return aperture * torch.exp(-1j * chi) Reconstruction pipeline ^^^^^^^^^^^^^^^^^^^^^^^ .. dropdown:: Code example: Reconstruction pipeline :icon: info :color: info .. code-block:: python def reconstruct( self, aberration_coefs: dict = None, rotation_angle_rad: float = None, upsampling_factor: int = 1, deconvolution_kernel: str = 'full', # 'full'=SSB, 'quadratic'=parallax max_batch_size: int = None, ): """ Perform SSB/parallax/BF reconstruction. """ # Setup coordinate grids with rotation kxa, kya = spatial_frequencies(self.gpts, self.sampling, rotation_angle) qxa, qya = upsampled_frequencies(upsampling_factor) # Compute complex probe at all k-vectors cmplx_probe = evaluate_probe(alpha, phi, ...) # Process each scan position in batches for batch_idx in batcher: # Get bright-field data in Fourier space vbf_fourier = torch.fft.fft2(self.vbf_stack[batch_idx]) # Compute transfer function Γ(Q,k) for this batch operator = self._compute_gamma_operator( kxa, kya, qxa, qya, aberration_coefs, cmplx_probe, batch_idx ) # Apply transfer function (deconvolution) fourier_factor = vbf_fourier * operator # Inverse FFT to get real-space reconstruction corrected_stack[batch_idx] = torch.fft.ifft2(fourier_factor) # Extract phase for SSB self.corrected_stack = corrected_stack.angle() # Mean across all positions self.mean_corrected_bf = self.corrected_stack.mean(dim=0) return self Usage examples ^^^^^^^^^^^^^^^ These examples demonstrate how to apply the concepts from :ref:`ssb_motivation`, :ref:`sideband_formation`, and :ref:`practical_guide` using the quantem library. Method 1: Extract probe from vacuum region (recommended) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If your scan includes vacuum (sample edge or hole), measure the probe directly: .. dropdown:: Code example: Vacuum probe extraction :icon: info :color: info .. code-block:: python import quantem as em import numpy as np # Load 4D-STEM dataset dataset = em.io.read_emdfile_to_4dstem('data.h5') # Method 1: Identify vacuum region in your scan # (assume vacuum is in the top-left corner, positions 0-10) vacuum_positions = dataset[0:10, 0:10] # Adjust to your data # Extract probe from vacuum region # Over vacuum, O(x,y) = 1, so the sideband directly measures P(k) vacuum_probe = extract_probe_from_vacuum( vacuum_positions, energy=60e3, semiangle_cutoff=34.0, ) # Use measured probe for reconstruction on sample region direct_ptycho = em.diffractive_imaging.DirectPtychography.from_dataset4d( dataset, energy=60e3, semiangle_cutoff=34.0, measured_probe=vacuum_probe, # Use measured probe device='mps', ) result = direct_ptycho.reconstruct(deconvolution_kernel='full') em.visualization.show_2d(result.obj, title='Phase [rad]') Method 3: Analytical probe from microscope parameters ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ For quick tests when vacuum is not available: .. dropdown:: Code example: Analytical probe reconstruction :icon: info :color: info .. code-block:: python import quantem as em # Load 4D-STEM dataset dataset = em.io.read_emdfile_to_4dstem('data.h5') # Initialize with analytical probe from microscope parameters direct_ptycho = em.diffractive_imaging.DirectPtychography.from_dataset4d( dataset, energy=60e3, # 60 keV semiangle_cutoff=34.0, # 34 mrad convergence angle aberration_coefs={'C10': 0}, # Near-zero defocus rotation_angle_deg=0, # Assume aligned device='mps', # 'mps' (Apple), 'cuda' (NVIDIA), or 'cpu' ) # Reconstruct the object phase result = direct_ptycho.reconstruct( deconvolution_kernel='full', # Full SSB deconvolution ) # The result is the phase map φ(x,y) from the conceptual section em.visualization.show_2d(result.obj, title='Phase [rad]') With hyperparameter optimization (recommended) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When you don't have a reference measurement, optimize the probe parameters from your data using the self-consistency loss described above: .. dropdown:: Code example: Hyperparameter optimization :icon: info :color: info .. code-block:: python import optuna import numpy as np # Optimize aberrations and rotation angle # Use narrow ranges around expected values for faster convergence sampler = optuna.samplers.TPESampler(seed=42) direct_ptycho.optimize_hyperparameters( aberration_coefs={ 'C10': (-100, 100), # Search defocus [Å] - narrow for corrected scope 'C12': (0, 50), # Search astigmatism [Å] 'phi12': (-np.pi/2, np.pi/2), # Astigmatism angle [rad] }, rotation_angle_deg=(-10, 10), # Search rotation [degrees] - typically small deconvolution_kernel='quadratic', # Use faster parallax method for optimization sampler=sampler, n_trials=30, # 30 trials usually sufficient with narrow ranges ) # Note: For wider exploration (uncorrected microscope, uncertain defocus): # Use C10: (-400, 400), rotation_angle_deg: (-180, 180), n_trials: 50-100 # Results show optimized parameters print(f"Optimized C10: {direct_ptycho._optimized_parameters['C10']:.1f} Å") print(f"Optimized C12: {direct_ptycho._optimized_parameters['C12']:.1f} Å") print(f"Optimized rotation: {direct_ptycho._optimized_parameters['rotation_angle_deg']:.2f}°") # Reconstruct with optimized parameters result = direct_ptycho.reconstruct_with_optimized_parameters( deconvolution_kernel='full', # Full SSB for final result ) # This is the best estimate of φ(x,y) given your data em.visualization.show_2d(result.obj, title='Optimized phase [rad]') Comparing reconstruction methods ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Compare different approximations discussed in :ref:`phase_extraction`: .. dropdown:: Code example: Method comparison (BF vs Parallax vs SSB) :icon: info :color: info .. code-block:: python # 1. Bright-field: No deconvolution, intensity only bf_result = direct_ptycho.reconstruct_with_optimized_parameters( deconvolution_kernel='none', use_center_of_mass_weighting=False, ) # 2. Parallax: Quadratic approximation to full SSB parallax_result = direct_ptycho.reconstruct_with_optimized_parameters( deconvolution_kernel='quadratic', ) # 3. SSB: Full deconvolution as described conceptually ssb_result = direct_ptycho.reconstruct_with_optimized_parameters( deconvolution_kernel='full', ) # Compare amplitudes and phases em.visualization.show_2d([ [np.abs(bf_result.obj), np.abs(parallax_result.obj), np.abs(ssb_result.obj)], [np.angle(bf_result.obj), np.angle(parallax_result.obj), np.angle(ssb_result.obj)] ], title=[ ['BF amplitude', 'Parallax amplitude', 'SSB amplitude'], ['BF phase', 'Parallax phase', 'SSB phase'] ]) This comparison shows: - **BF**: Fast but limited resolution, no deconvolution - **Parallax**: Faster approximation, good for optimization - **SSB**: Full reconstruction following the conceptual framework References ^^^^^^^^^^^ Key papers on SSB and direct ptychography: 1. Ophus, C. (2019). Four-dimensional scanning transmission electron microscopy (4D-STEM): from scanning nanodiffraction to ptychography and beyond. *Microscopy and Microanalysis*, 25(3), 563-582. 2. Chen, Z., et al. (2021). Electron ptychography achieves atomic-resolution limits set by lattice vibrations. *Science*, 372(6544), 826-831. 3. Rodenburg, J. M., & Bates, R. H. T. (1992). The theory of super-resolution electron microscopy via Wigner-distribution deconvolution. *Philosophical Transactions of the Royal Society A*, 339(1655), 521-553. 4. Maiden, A. M., & Rodenburg, J. M. (2009). An improved ptychographical phase retrieval algorithm for diffractive imaging. *Ultramicroscopy*, 109(10), 1256-1262.