🔭 astrophysics

Fast partial-sky spherical harmonic transforms

This paper presents a novel, efficient algorithm for partial-sky spherical harmonic transforms that achieves significant speed-ups (3–10x) on small, high-resolution sky regions, thereby overcoming the limitations of arbitrary pixelizations for demanding Cosmic Microwave Background analyses like those from SPT-3G and BICEP3.

Original authors: Julien Carron, Martin Reinecke

Published 2026-03-19

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Julien Carron, Martin Reinecke

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Picture: Mapping the Sky Without the Whole Map

Imagine you are an astronomer trying to study a tiny, specific patch of the night sky—maybe a small cluster of galaxies or a specific region where the Cosmic Microwave Background (the afterglow of the Big Bang) is very interesting.

Usually, to analyze data on a sphere (like the sky), scientists use a mathematical tool called a Spherical Harmonic Transform (SHT). Think of this as a way to translate a 3D map of the sky into a list of musical notes (frequencies) that describe the patterns on that map.

The Problem:
Standard methods for doing this translation are like trying to listen to a single violin in a symphony orchestra by recording the entire orchestra, the whole concert hall, and every single instrument playing at once. Even if you only care about that one violin, the computer has to process the data for the whole globe. This is incredibly slow and wasteful, especially when you need to do this calculation thousands of times for high-precision science.

The Solution:
Julien Carron and Martin Reinecke have invented a new, faster way to do this. They figured out how to "zoom in" on just the small patch of sky you care about, stretch it out to fill a virtual screen, and do the math only on that small area.

The Core Idea: The "Sticky Note" Analogy

Here is how their new algorithm works, broken down into three simple steps:

1. The "Stretching" Trick (Magnifying the Patch)

Imagine you have a tiny, crumpled piece of paper with a drawing on it (your small sky patch). If you try to analyze the drawing while it's crumpled, it's hard.

Old Way: You try to analyze the drawing while it's still crumpled, but you have to look at the whole table it's sitting on to understand the context.
New Way: You take that tiny crumpled paper and stretch it out until it fills a whole large poster board. Now, the tiny details of your drawing are huge and easy to see.

In the paper's math, they "stretch" the small patch of the sky so it looks like a bigger, simpler shape. Because the features look bigger, the computer doesn't need to use as many "notes" (frequencies) to describe them. It's like describing a giant mountain vs. a tiny pebble; the mountain needs fewer complex details to be recognized than the pebble does when you are looking at it from far away.

2. The "Pole" Advantage (The Ice Skater)

Once they stretch the patch, they rotate the sky so that the center of their patch sits right at the North Pole.

The Analogy: Imagine an ice skater spinning. If they spin near the center of the ice (the pole), they only have to move a tiny distance to go all the way around. If they spin near the edge of the rink (the equator), they have to skate a huge circle.
The Benefit: By placing the patch at the pole, the "circles" the computer has to calculate are tiny. This means the computer does much less work. The paper calls this "Polar Optimization." It's like realizing you only need to walk around a small coffee table instead of a giant stadium to get your exercise.

3. The "Soft Edge" (The Window)

There's a catch: If you just cut out a square piece of the sky and stretch it, the edges look jagged and messy, which creates mathematical "noise" or errors.

The Fix: They use a special mathematical "window" (a filter) that acts like a soft fade-out. Imagine taking a photo and using a blur tool on the very edges so the image fades gently into blackness instead of cutting off sharply.
The Result: This ensures the math stays smooth and accurate. They use a very specific type of mathematical curve (based on something called "Slepian functions") that is the perfect balance between being smooth and keeping the data accurate.

Why Does This Matter? (The Real-World Impact)

The authors tested this on data from the South Pole Telescope, which is looking for evidence of primordial gravitational waves (ripples from the very beginning of the universe).

The Speed Boost: On the large areas they studied, the new method was 3 to 10 times faster. On the smallest, most detailed patches, it was up to 40 times faster.
The "No-Compromise" Benefit: Usually, to get this speed, scientists have to make approximations (like pretending the sky is flat, which isn't true). This new method is fast and perfectly accurate. It removes the need to force data into rigid grids or make "flat-sky" guesses.
The Outcome: Because the calculations are so much faster, scientists can run complex simulations on regular computers instead of needing massive, expensive supercomputers. This makes high-precision cosmology accessible to more researchers.

Summary

Think of this paper as inventing a specialized microscope for the sky. Instead of scanning the entire universe to find a tiny detail, they figured out how to isolate that detail, magnify it, and analyze it with a fraction of the effort, all while keeping the picture crystal clear.

In one sentence: They found a way to do complex 3D sky math on small patches by stretching the patch, spinning it to the pole, and using a soft mathematical filter, making the process 3 to 40 times faster without losing any accuracy.

1. Problem Statement

The paper addresses the computational bottleneck associated with Spherical Harmonic Transforms (SHTs) when applied to partial-sky data (small, arbitrary regions of the celestial sphere).

Context: Modern cosmological surveys, such as the South Pole Telescope (SPT-3G) and BICEP3, target small sky fractions to study Cosmic Microwave Background (CMB) lensing and primordial gravitational waves (B-modes). These analyses require high-resolution, numerically intensive iterative methods (e.g., Wiener filtering, quadratic estimators) involving repeated forward and inverse SHTs.
Limitations of Existing Methods:
- Full-sky algorithms: Standard efficient SHTs (like those based on the Fourier-sphere method) evaluate maps on a full set of $\ell_{\text{max}} + 1$ iso-latitude rings. When the data is confined to a small patch, calculating the entire sphere is computationally wasteful.
- Flat-sky approximations: While fast for small patches, these approximations break down at high resolutions or for large patches, introducing systematic errors in precision cosmology.
- Arbitrary pixelizations: Existing methods often struggle with arbitrary point distributions or require specific pixelizations (like HEALPix) that may not align with the observation geometry.

The goal is to develop an algorithm that performs SHTs on arbitrary locations within small sky patches with speed-ups of 3x to 10x (or more) compared to full-sky methods, while maintaining high numerical accuracy.

2. Methodology

The authors propose a novel algorithm that combines the Fourier-sphere method (Reinecke et al., 2023) with a local stretching and apodization technique.

Core Concept: Stretching and Periodicity

Instead of evaluating the map on the full sphere, the algorithm treats the small patch as a periodic function on a "stretched" Fourier sphere.

Coordinate Transformation: The patch is re-centered so its center lies at the North Pole ( $\theta = 0$ ).
Stretching: The physical region of interest (radius $\theta^*$ ) is mathematically "stretched" to cover a larger spherical cap (radius $\theta_{\text{max}}$ ). This effectively magnifies the features within the patch, reducing their effective spatial frequency.
Band-limit Reduction: Because the features are magnified, the effective band-limit required to represent the map decreases from $\ell_{\text{max}}$ to a lower value $L_{\text{max}}$ .
Apodization (Windowing): To ensure the stretched map remains periodic and smooth (avoiding discontinuities at the new boundary), a weight function $w(\Theta)$ is applied. This function is unity within the region of interest and smoothly decays to zero at the edge of the stretched cap.

Algorithmic Steps

Define Parameters: Choose an apodization length ( $\epsilon_{\text{apo}}$ ) and a window function based on the desired error tolerance ( $\epsilon$ ).
Construct Weighted Function: Define a new function $\tilde{f}(\Theta, \phi) = w(\Theta) f(\theta(\Theta), \phi)$ , where $\Theta$ is the stretched coordinate.
Reduced Ring Evaluation: Compute the SHT of $\tilde{f}$ using the standard Fourier-sphere method, but only on $L_{\text{max}} + 1$ rings (where $L_{\text{max}} < \ell_{\text{max}}$ ).
Polar Optimization: By centering the patch at the pole, the algorithm exploits the fact that rings near the pole are physically shorter. This allows for polar optimization, where high- $m$ modes (longitudinal frequencies) that exceed the physical resolution of the small ring are discarded, further reducing computational cost.
Reconstruction: The result is mapped back to the original coordinates.

Window Function Selection

The authors utilize Slepian functions (or prolate spheroidal wave functions) and their approximation, the Kaiser-Bessel window (specifically an "exponential of semi-circle" kernel). These functions are chosen because they offer the optimal trade-off between compact support in real space (minimizing the apodization region) and rapid decay in frequency space (minimizing spectral leakage).

3. Key Contributions

Novel Partial-Sky Algorithm: A method that eliminates the need for full-sky evaluation or flat-sky approximations for small patches, handling arbitrary pixelizations naturally.
Theoretical Framework for Error Control: The paper derives explicit formulas for the error introduced by the stretching and apodization process. It relates the window function parameters to the desired numerical tolerance ( $\epsilon$ ), allowing users to tune the algorithm for specific precision requirements.
Polar Optimization Integration: The method explicitly leverages the geometry of the patch (centering it at the pole) to maximize the efficiency of the Legendre transforms, a feature often underutilized in general-purpose SHT codes.
Implementation Efficiency: The algorithm is implemented as a minor modification to the existing high-performance ducc0 library, preserving all low-level CPU vectorizations and optimizations.

4. Results and Performance

The authors benchmarked the algorithm on three specific sky regions relevant to CMB research:

SPT-3G Main Field: ~1500 square degrees.
BICEP3 Deep Footprint: ~800 square degrees.
Euclid Deep Field South (EDFS): ~57 square degrees (very small patch).

Performance Metrics (Speed-up factors compared to standard full-sky ducc0.sht):

SPT-3G: Speed-up of 3.2x – 3.9x.
BICEP3: Speed-up of 6.4x – 9.0x.
EDFS (Small Patch): Speed-up of 23x – 43x.
Adjoint Operations: Similar speed-ups were observed for the adjoint (inverse) transforms.

Key Observations:

The method outperforms even truncated full-sky transforms (which only compute rings intersecting the patch) by a factor of ~1.5 to 3x, due to the combined benefits of reduced band-limit and polar optimization.
For very small patches (like EDFS), the method is competitive with or faster than 2D Fast Fourier Transforms (FFT) on flat-sky grids, but without the geometric distortions of flat-sky approximations.
The error is controllable and scales with the chosen tolerance parameter; for a tolerance of $10^{-7}$ , the errors are negligible for cosmological applications.

5. Significance

This work is significant for the future of high-precision cosmology, particularly for:

CMB B-mode Searches: Experiments like SPT-3G and future CMB-S4 require iterative lensing reconstruction to remove foregrounds and isolate primordial gravitational wave signals. The speed-up allows these computationally expensive iterations to be performed on standard workstations rather than requiring massive High-Performance Computing (HPC) clusters.
Flexible Analysis: It enables rigorous spherical analysis on arbitrary, non-standard survey footprints without resorting to approximations that could bias cosmological parameter estimation.
Software Accessibility: The method is integrated into the lenspyx software package, making these high-performance capabilities accessible to the broader astrophysics community.

In summary, Carron and Reinecke have developed a mathematically rigorous and computationally efficient solution that bridges the gap between the accuracy of full-sky spherical harmonics and the speed of flat-sky approximations, specifically tailored for the small, high-resolution patches characteristic of modern CMB experiments.