Multi-Wavelength Spherical Harmonics Maps

• Multi-wavelength spherical harmonics maps are spherical data products that represent angular structure via harmonic coefficients while retaining spectral diversity across frequencies.
• They enable the fusion of disparate data from instruments by employing varied techniques such as spherical Fourier-Bessel transforms, m-mode decomposition, and adaptive reconstruction algorithms.
• Advanced implementations leverage multi-resolution HEALPix schemes and fast spherical harmonic transforms to optimize multi-frequency and multi-messenger sky analyses.

Multi-wavelength spherical harmonics maps are spherical data products in which angular structure is represented on the sphere—typically through spherical harmonic coefficients—and spectral variation is retained across multiple frequencies, wavelengths, or spectral lines. In contemporary practice, the term encompasses several closely related objects: per-band coefficient sets such as $a_{\ell m}$ or $b_{\ell m}(k_0)$, harmonic-space combinations of maps from instruments with different beams and noise properties, spherical Fourier-Bessel representations that keep both angular and radial information, and multi-resolution HEALPix products used to exchange localized sky information across instruments with sharply different angular resolution (Chown et al., 2018, Carozzi, 2015, Liu et al., 2016, Martinez-Castellanos et al., 2021).

1. Mathematical forms and representational regimes

At the scalar level, a spherical map is expanded in the usual harmonic basis. In the interferometric spherical formulation, the sky intensity is written as

$$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$

while the measured visibilities admit a corresponding spherical-wave expansion, with the coefficients related by

$\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$

This yields a direct identification of the spherical harmonic content of the sky and, when repeated at each observing frequency $k_0$, gives a family of wavelength-dependent harmonic maps $b_{\ell m}(k_0)$ rather than a single monochromatic sky representation (Carozzi, 2015).

For three-dimensional line-intensity surveys, the natural generalization is the spherical Fourier-Bessel basis. The temperature field $T(\mathbf r)$ is decomposed as

$$\overline{T}_{\ell m}(k)\equiv \sqrt{\frac{2}{\pi}}\int d\Omega\int dr\,r^2\,j_\ell(kr)\,Y_{\ell m}^*(\hat{\mathbf r})\,T(\mathbf r),$$

and the corresponding power spectrum proxy is $S_\ell(k)$, which separates angular structure indexed by $\ell$ from radial structure indexed by $b_{\ell m}(k_0)$0. This construction is the curved-sky analogue of $b_{\ell m}(k_0)$1, but it remains well defined beyond the narrow-field, flat-sky approximation (Liu et al., 2016).

The scalar formalism is not exhaustive. For arbitrary tensorial systems on $b_{\ell m}(k_0)$2, spin-weighted spherical harmonics provide a unified spectral representation: $b_{\ell m}(k_0)$3 In that setting, the spin-raising and spin-lowering operators act algebraically in spectral space, for example

$b_{\ell m}(k_0)$4

This is the relevant extension for vector, tensor, and polarization-like fields rather than scalar intensity alone (Beyer et al., 2013).

Representation	Core object	Representative use
Spherical harmonic map	$b_{\ell m}(k_0)$5, $b_{\ell m}(k_0)$6	CMB and wide-field imaging
Spherical Fourier-Bessel map	$b_{\ell m}(k_0)$7	Intensity mapping
Spin-weighted spectral map	$b_{\ell m}(k_0)$8	Tensorial fields on $b_{\ell m}(k_0)$9
Multi-resolution spherical map	order-indexed HEALPix pixels	Multi-messenger localization

2. Reconstruction from visibilities and broadband observations

One major route to multi-wavelength spherical harmonics maps begins directly from interferometric measurements. In the Spherical Wave Harmonic Transform, the visibility field satisfies the Helmholtz equation, its eigenfunctions are $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$0, and the coefficient extraction formula for a dataset sampled at positions $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$1 is

$$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$2

The method gives an exact and explicit solution to the scalar radio-interferometer imaging equation on the sphere, is valid for non-coplanar interferometer configurations, and removes any special status for the conventional $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$3-term by formulating the problem spherically rather than cartesianly (Carozzi, 2015).

For transit interferometers, the relevant construction is the spherical harmonics $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$4-mode decomposition. After Fourier transforming in sidereal time, the visibility system factorizes into independent blocks,

$$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$5

so a very large inverse problem becomes a sequence of smaller problems, one for each $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$6. The maximum-likelihood estimator is

$$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$7

and the same framework yields the point spread function, transfer function, noise covariance matrix, and noise power spectrum $$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$8. The PAON-4 and Tianlai analyses further show that regular spacing produces more pronounced side lobes and more structure in the noise power spectrum than irregular spacing, particularly in the north-south direction (Zhang et al., 2016).

A different reconstruction problem arises when a telescope records broadband time-ordered data with a wavelength-dependent, asymmetric beam. In that setting, the full polychromatic sky is formally described by

$$B(\Omega_k)=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} b_{\ell m}Y_{\ell m}(\Omega_k),$$9

but the full position-by-wavelength inverse problem is underconstrained. The practical objective becomes the recovery of optimal linear combinations of wavelengths, or “colors,” defined by vectors $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$0. These are found as the eigenvectors of

$\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$1

with the largest eigenvalues corresponding to the best-constrained colors. When the scan is isotropic, the spherical harmonic basis makes the inverse noise covariance block diagonal in $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$2; simulations show that maps in multiple colors can be reconstructed accurately from both full-sky and partial-sky scans, and that at least two independent color maps are recoverable if the antenna pattern is sufficiently asymmetric and wavelength-dependent (Cantrall et al., 2022).

3. Harmonic-space fusion across instruments, frequencies, and lines

A central use of multi-wavelength spherical harmonics maps is the combination of surveys with complementary angular response. In the South Pole Telescope and Planck analysis, three maps of the millimeter-wave sky were formed by combining SPT-SZ data in bands centered at 95, 150, and 220 GHz with Planck temperature data in the 100, 143, and 217 GHz bands. The fusion is performed directly in spherical harmonic space: $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$3 with inverse-noise weights

$\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$4

This weighting places most of the large-angular-scale information in Planck and the smaller-scale information in SPT-SZ. After combination, a final Gaussian beam with $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$5 is applied and the map is reconstructed by an inverse spherical harmonic transform. Simulations show agreement of the output angular power spectra with the input at the $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$6 level over most scales, with slight excess at very high $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$7 above $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$8 (Chown et al., 2018).

For line-intensity mapping, multi-wavelength generality is built into the formalism itself. Because the mapping from observed frequency to comoving distance depends on the spectral line, separate spherical Fourier-Bessel decompositions can be constructed for different rest-frame transitions such as the $\tilde{v}_{\ell m}=4\pi(-i)^\ell b_{\ell m}.$9 line, carbon monoxide rotational lines, or singly ionized carbon fine-structure lines. This allows different lines to be placed on a common $k_0$0 basis and makes interloper contamination diagnostically useful rather than merely adversarial: if a contaminating line is mapped with the wrong frequency-to-distance relation, the inferred $k_0$1 develops a characteristic anisotropy in its $k_0$2- and $k_0$3-dependence (Liu et al., 2016).

The same curved-sky formalism also clarifies foreground structure. Spectrally smooth foregrounds are confined to low radial $k_0$4 in the idealized basis, while interferometric chromaticity spreads them into the foreground wedge in the $k_0$5 plane. This is the spherical counterpart of familiar flat-sky diagnostics, but without the need to partition a wide survey into locally planar patches (Liu et al., 2016).

4. Pixelization, multi-resolution storage, and map arithmetic

Harmonic representations are only one side of the problem; storage and exchange formats matter when different instruments resolve very different angular scales. HEALPix has become a standard in high-energy and gravitational-wave astronomy and is widely used for spherical data, but its homogeneous all-sky resolution is a limitation when small, highly localized regions require much finer sampling than the rest of the sky. The multi-resolution extension implemented in mhealpy addresses this by allowing only selected regions to be recursively subdivided, while the remainder of the sphere stays coarse (Martinez-Castellanos et al., 2021).

In this framework, the finest map order is $k_0$6, with $k_0$7, and each pixel is identified by an order $k_0$8 and a NESTED index $k_0$9. Several indexing schemes coexist. A Multi-Order List groups pixels by order; a Range Set represents pixels through index ranges at maximum resolution, with

$b_{\ell m}(k_0)$0

and NUNIQ packs $b_{\ell m}(k_0)$1 into a single unique integer for compact storage and serialization. The library supports efficient pixel querying, arithmetic operations between maps, adaptive mesh refinement, plotting, and serialization into FITS files compatible with IVOA MOC and NUNIQ conventions (Martinez-Castellanos et al., 2021).

Map arithmetic in a multi-resolution setting requires explicit semantics for how pixel values behave under refinement. mhealpy distinguishes density-like quantities, such as $b_{\ell m}(k_0)$2 or probability density, from histogram-like quantities, such as counts or probability. When a pixel at order $b_{\ell m}(k_0)$3 is promoted to order $b_{\ell m}(k_0)$4, it is split into $b_{\ell m}(k_0)$5 children; density-like values are copied to every child, whereas histogram-like values are divided by $b_{\ell m}(k_0)$6. Binary operations are carried out after mesh alignment, and by default the output mesh is as fine as necessary to avoid information loss (Martinez-Castellanos et al., 2021).

The scientific relevance is immediate in joint multi-wavelength and multi-messenger analyses. The paper’s examples include a map for the localization of SSS17a/AT 2017gfo in which the source region is represented at sub-arcsecond resolution with $b_{\ell m}(k_0)$7 while the rest of the sky remains coarse, and a GW170817 localization in which adaptive refinement makes pixel density track probability density. A notable limitation is explicit: spherical harmonic transforms are not presently implemented in mhealpy for multi-resolution maps, even though most other HEALPix-style operations are supported in a mesh-agnostic way (Martinez-Castellanos et al., 2021).

5. Algorithms, scaling laws, and throughput for many maps

As multi-wavelength analyses increasingly involve hundreds or thousands of sky maps, the spherical harmonic transform itself becomes a limiting operation. The HP2SPH algorithm addresses this on HEALPix grids by combining the double Fourier sphere method, a non-uniform fast Fourier transform, and Slevinsky’s fast spherical harmonic transform. After upsampling and shifting HEALPix rings onto a tensor-product grid, the method computes bivariate Fourier coefficients and then converts them to spherical harmonic coefficients. For a grid with $b_{\ell m}(k_0)$8 pixels, the stated complexity is $b_{\ell m}(k_0)$9 per map, with an initial setup cost of $T(\mathbf r)$0. This compares favorably with the $T(\mathbf r)$1 complexity of current HEALPix methods when many maps must be analyzed on the same grid (Drake et al., 2019).

The same work reports an accuracy gain rather than only a speed gain. Numerical experiments show better accuracy over the entire angular power spectrum of synthetic data than the current HEALPix methods, with a convergence rate at least two times higher. For real CMB data, the method produces results that are visually the same as the HEALPix methods, which indicates robustness in the presence of realistic noise rather than only idealized smooth inputs (Drake et al., 2019).

A complementary acceleration strategy targets large batches of maps directly. The fastSHT scheme reformulates the dominant part of the transform for many sky maps as a BLAS-3 matrix-matrix multiplication rather than repeated BLAS-2 matrix-vector operations. The matrix form is

$T(\mathbf r)$2

where $T(\mathbf r)$3 contains FFT-processed ring data for many maps and $T(\mathbf r)$4 contains the associated Legendre factors. This enables 2–10 times speedup on CPUs and up to 30 times speedup on GPUs relative to the standard Libsharp-HEALPix program. The method is explicitly optimized for large numbers of maps rather than single-map workloads, and GPU memory becomes a limiting factor at very high $T(\mathbf r)$5 (Tian et al., 2022).

These scaling results are particularly consequential for multi-frequency CMB analyses, component separation, null tests, and simulation campaigns, all of which require repeated transforms over large map ensembles. A plausible implication is that the practical unit of analysis is shifting from the single all-sky map to the map stack, with throughput determined as much by batching strategy as by asymptotic complexity.

6. Scientific uses, misconceptions, and cross-domain extensions

In astronomy, the principal scientific roles of multi-wavelength spherical harmonics maps are now well differentiated. In CMB and millimeter-wave cosmology they are used to combine instruments with complementary beam and noise properties (Chown et al., 2018). In radio interferometry they support exact or maximum-likelihood reconstruction on the sphere, including non-coplanar or transit geometries (Carozzi, 2015, Zhang et al., 2016). In intensity mapping they provide a curved-sky replacement for flat-sky $T(\mathbf r)$6 analyses and a natural basis for foreground and interloper diagnostics (Liu et al., 2016). In multi-messenger astronomy they motivate multi-resolution HEALPix products that can represent a gravitational-wave probability map and a sharply localized optical counterpart within a common spherical framework (Martinez-Castellanos et al., 2021).

Several recurring misunderstandings are resolved by the literature. A spherical harmonics map is not necessarily a uniformly pixelized all-sky image: multi-resolution HEALPix maps are explicitly designed to violate that assumption while preserving hierarchical spherical structure (Martinez-Castellanos et al., 2021). Conversely, a multi-resolution spherical map is not automatically a harmonic-domain object: mhealpy supports querying, plotting, and arithmetic, but spherical harmonic transforms are not presently implemented for such maps (Martinez-Castellanos et al., 2021). Likewise, the usual flat-sky decomposition into perpendicular and parallel Fourier modes is not the general definition of a spectral sky map; it is an approximation whose validity is limited to narrow fields, whereas the spherical Fourier-Bessel basis is constructed precisely to avoid that restriction (Liu et al., 2016). In interferometric imaging, the cartesian $T(\mathbf r)$7-term is not a fundamental physical obstruction on the full sphere but a feature of the planar approximation; in the spherical formalism it has no special significance (Carozzi, 2015).

Outside astronomy, closely related multi-channel spherical map constructions appear in illumination modeling and rendering. MixLight uses spherical harmonics and spherical Gaussians jointly, with spherical harmonics modeling low-frequency ambient light and spherical Gaussians modeling high-frequency light sources; the parameters are estimated per RGB channel, and experiments on the Web Dataset show that this parametric method has better generalization performance than non-parametric methods (Ji et al., 2024). Harmonics Virtual Lights project spherical-light contributions onto spherical harmonics and evaluate outgoing luminance in closed form, with $T(\mathbf r)$8 complexity for circular symmetric lobes and $T(\mathbf r)$9 in the general case, where $$\overline{T}_{\ell m}(k)\equiv \sqrt{\frac{2}{\pi}}\int d\Omega\int dr\,r^2\,j_\ell(kr)\,Y_{\ell m}^*(\hat{\mathbf r})\,T(\mathbf r),$$0 is the number of SH bands (Mézières et al., 2022). These results are not astronomical, but they demonstrate that multi-channel spherical harmonic maps are a general representational technology rather than a discipline-specific convention.

Taken together, these developments define a field in which “multi-wavelength spherical harmonics maps” are not a single file format or algorithm, but a family of spherical representations linking harmonic analysis, inverse problems, map fusion, adaptive pixelization, and high-throughput computation. The unifying idea is that angular information on the sphere and spectral diversity across bands should be treated jointly, with the choice of basis, reconstruction method, and storage model determined by survey geometry, instrumental response, and the scientific scales of interest.