Generalized Needlet ILC (GNILC)

• GNILC is a component-separation method that extends traditional ILC by leveraging localized, multi-dimensional needlet analysis and adaptive model-order selection.
• It employs needlet localization, local covariance whitening, and eigen-decomposition to dynamically separate foreground contamination from target signals in CMB and HI observations.
• GNILC has been effectively applied to Planck dust/CIB separation and HI intensity mapping, providing robust, scale-dependent signal recovery with minimized foreground leakage.

Generalized Needlet Internal Linear Combination (GNILC) is a component-separation method that extends needlet-based Internal Linear Combination from the case of a single fixed mixing vector to the case of a locally varying, multi-dimensional signal subspace. In the Planck formulation, it is a multi-dimensional ILC operating in a spherical needlet domain, with a data-driven, locally varying signal subspace dimension $m$ selected using model-order statistics such as the Akaike Information Criterion (AIC), after whitening the local data covariance with a prior nuisance covariance built from CIB, CMB, and noise angular power spectra (Collaboration et al., 2016). In 21-cm intensity mapping, GNILC preserves an HI subspace rather than a single spectral vector, using a prior HI covariance to identify and project out foreground-dominated directions locally in angular scale and sky position (Olivari et al., 2015). Subsequent work has adapted GNILC to BINGO and SKA-like surveys, extended it to the expanded GNILC (eGNILC) framework, used GNILC-inspired diagnostics in constrained polarization pipelines, and employed GNILC foreground reconstructions to build residual templates for cosmological-parameter debiasing (Mericia et al., 2022, Dai et al., 2024, Carones et al., 2024, Carones, 23 Oct 2025).

1. Method class and defining idea

GNILC belongs to the ILC family, but it is not the standard single-component ILC. Standard ILC minimizes variance subject to unit response to a known mixing vector, and NILC performs the same operation locally in needlet space. GNILC is “generalized” because it uses priors on nuisance components or, in HI applications, a prior HI covariance, to estimate locally in needlet space and across the sky the dimension and basis of the informative subspace where correlated astrophysical emission dominates over nuisances (Collaboration et al., 2016, Carones et al., 2024).

In the Planck high-frequency application, GNILC was developed because Galactic thermal dust emission and cosmic infrared background anisotropies have very similar modified-blackbody spectral shapes, making a purely spectral separation inefficient. The method therefore exploits spatial information, specifically well-measured angular power spectra of nuisance components, in addition to multi-frequency data, to disentangle thermal dust from CIB anisotropies (Collaboration et al., 2016). In HI intensity mapping, the motivation is different but structurally analogous: the cosmological 21-cm signal does not admit a single fixed spectral scaling vector across channels, so a naive ILC is not directly applicable; GNILC instead uses a prior on the HI power spectrum to whiten the data and infer, per needlet band and per sky location, the effective dimensionality of the foreground subspace, which is then projected out (Caro et al., 2 Sep 2025).

A common misconception is that GNILC is simply NILC with a different covariance estimate. The published formulations indicate a stronger distinction. Unlike standard NILC, which targets a fixed mixing vector, and unlike MILCA, which uses fixed dimensionality, GNILC is explicitly multi-dimensional and locally adaptive, with model-order selection varying over needlet scale and sky position (Collaboration et al., 2016). In HI work, the preserved object is an HI subspace, not a single channel-law response (Olivari et al., 2015).

2. Needlet localization and local covariance estimation

GNILC operates in a spherical needlet frame to localize the separation simultaneously in angular scale and sky position. In the Planck formulation, for a sky field with spherical harmonic coefficients $a_{\ell m}$, the needlet coefficient at scale $j$ and sampling location $\xi_{j,k}$ is

$$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$

where the needlet windows satisfy

$\sum_j h_j(\ell)^2=1$

to conserve power in synthesis (Collaboration et al., 2016).

In HI and polarization applications, analogous needlet-domain filtering is used. The polarized multifrequency data can be filtered by compact bandpass windows $b_\ell^j$ to produce needlet coefficient maps

$$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$

and local covariances are computed per needlet scale and spatial domain by Gaussian convolution or local averaging (Carones et al., 2024). In HI work, local sample covariances are estimated from needlet coefficient maps over domains whose size is chosen to balance locality against finite-sample bias (Olivari et al., 2015).

The local empirical covariance is central. In the Planck dust/CIB implementation, for needlet maps $x_a^j(p)$ and $x_b^j(p)$, the covariance around pixel $a_{\ell m}$0 is

$a_{\ell m}$1

with $a_{\ell m}$2 a local Gaussian-weighted domain (Collaboration et al., 2016). In BINGO and related HI analyses, the same logic is used with frequency-frequency covariances built in local needlet domains, with the locality controlled by the needlet windows and the ILC bias parameter $a_{\ell m}$3 (Mericia et al., 2022).

This needlet localization is not a cosmetic implementation choice. It allows GNILC to adapt both to sky regions where foreground complexity changes and to multipole ranges where the signal-to-foreground ratio changes. The Planck work reports that in practice $a_{\ell m}$4 decreases at high latitude and at small scales where dust signal-to-noise is low; in some high-latitude, small-scale domains, $a_{\ell m}$5, which induces a spatially varying effective resolution (Collaboration et al., 2016).

3. Whitening, eigenmodes, model order, and the multi-dimensional ILC filter

The defining GNILC step is local subspace estimation by whitening and eigendecomposition. In the Planck thermal-dust application, the nuisance covariance is

$a_{\ell m}$6

and the whitened covariance is

$a_{\ell m}$7

If $a_{\ell m}$8 matches the true nuisance covariance, then $a_{\ell m}$9, so eigenvalues near unity reflect nuisance while significant departures above unity indicate local dust signal-to-nuisance power (Collaboration et al., 2016).

In HI intensity mapping, the whitening is written in the complementary way, using a prior HI covariance:

$j$0

or, in the SKAO formulation,

$j$1

Modes with eigenvalues $j$2 are foreground-dominated relative to the HI prior, and GNILC selects the number of foreground components locally by minimizing an AIC criterion (Olivari et al., 2015, Caro et al., 2 Sep 2025).

The Planck AIC expression is

$j$3

where $j$4 are the eigenvalues of the whitened covariance attributed to the nuisance subspace and $j$5 is the number of samples in the local needlet domain (Collaboration et al., 2016). In the 2015 HI formulation, the corresponding criterion is

$j$6

The selected $j$7 determines the locally varying foreground rank, and the complementary subspace is treated as the preserved signal subspace (Olivari et al., 2015).

Once the subspace basis has been identified, GNILC applies a multi-dimensional ILC filter. In the Planck notation, mapping the estimated signal subspace back to the original domain gives

$j$8

and the estimator preserving that subspace while minimizing variance is

$j$9

This generalizes the standard single-component ILC weight

$\xi_{j,k}$0

In the HI notation, the equivalent multi-dimensional weight matrix is

$\xi_{j,k}$1

or, for the reconstructed channel-space signal,

$\xi_{j,k}$2

(Collaboration et al., 2016, Olivari et al., 2015).

These equations clarify another frequent misunderstanding: GNILC does not require a parametric spectral model for the signal subspace itself. What it requires is a prior covariance for whitening and dimension selection. The preserved basis is estimated locally from the data after whitening, rather than fixed globally in advance (Collaboration et al., 2016, Olivari et al., 2015).

4. Planck dust and cosmic infrared background separation

A major mature application of GNILC is the separation of Galactic thermal dust emission from cosmic infrared background anisotropies in Planck 2015 temperature maps. The method was implemented because previous Planck dust products treated dust and CIB together and suffered noticeable CIB leakage, especially at high Galactic latitude, biasing dust spectral parameters and increasing their dispersion (Collaboration et al., 2016).

Using the Planck PR2 temperature maps, GNILC produced improved all-sky maps of Planck thermal dust emission with reduced CIB contamination at 353, 545, and 857 GHz. By reducing CIB contamination, it provided more accurate estimates of the local dust temperature and dust spectral index over the sky with reduced dispersion, especially at high Galactic latitudes above $\xi_{j,k}$3. The reported averages are $\xi_{j,k}$4 and $\xi_{j,k}$5 over the whole sky, and $\xi_{j,k}$6 and $\xi_{j,k}$7 on $\xi_{j,k}$8 of the sky at high latitudes (Collaboration et al., 2016).

The same reconstruction also enabled CIB extraction. Subtracting the new CIB-removed thermal dust maps from the CMB-removed Planck maps gives access to the CIB anisotropies over $\xi_{j,k}$9 of the sky at Galactic latitudes $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$0 according to the abstract, and the detailed analysis reports GNILC CIB maps over approximately $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$1 overall. In GHIGLS fields, the GNILC CIB correlates strongly with the Planck 2013 CIB maps with Pearson $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$2 and slopes near unity, while the GNILC CIB at 353 GHz shows negligible correlation with H I, $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$3, over $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$4 of the sky, indicating minimal residual Galactic contamination (Collaboration et al., 2016).

The angular power-spectrum behavior is also distinctive. The dust angular power spectrum at 353 GHz becomes steeper, approximately $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$5, with high-$$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$6 amplitude reduced by a factor of approximately $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$7 at $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$8 relative to P13, consistent with CIB removal. Residual maps, defined as Planck minus CMB minus GNILC dust, have power consistent with the sum of Planck CIB best-fit and Planck noise spectra at 353 GHz; at 545 and 857 GHz a small residual CIB remains (Collaboration et al., 2016).

The delivered products include CIB-removed GNILC thermal dust maps at 353, 545, and 857 GHz, GNILC CIB maps at the same frequencies, dust optical depth $$\beta_{j,k}=\sum_{\ell,m} h_j(\ell)\,a_{\ell m}\,Y_{\ell m}(\xi_{j,k}),$$9, spectral index $\sum_j h_j(\ell)^2=1$0, and temperature $\sum_j h_j(\ell)^2=1$1 maps at $\sum_j h_j(\ell)^2=1$2 resolution, and effective beam maps. The effective local resolution is spatially varying because the largest needlet scale with $\sum_j h_j(\ell)^2=1$3 sets the local beam; the Planck analysis reports that at least $\sum_j h_j(\ell)^2=1$4 of the sky retains $\sum_j h_j(\ell)^2=1$5 resolution (Collaboration et al., 2016).

5. HI intensity mapping and 21-cm foreground removal

In HI intensity mapping, GNILC addresses the fact that the cosmological 21-cm signal is faint while Galactic synchrotron, free-free emission, radio sources, and related contaminants dominate the observed variance by orders of magnitude. The 2015 HI analysis studies GNILC for a general HI intensity mapping experiment and reports that, for simulated radio observations including HI emission, Galactic synchrotron, Galactic free-free, radio sources and $\sum_j h_j(\ell)^2=1$6 thermal noise, the method reconstructs the HI power spectrum for multipoles $\sum_j h_j(\ell)^2=1$7 with $\sum_j h_j(\ell)^2=1$8 accuracy on $\sum_j h_j(\ell)^2=1$9 of the sky for a redshift $b_\ell^j$0 (Olivari et al., 2015).

The same work emphasizes that GNILC is robust to the complexity of the foregrounds. In the detailed simulations, adding free-free emission, point sources, and a spatially varying synchrotron spectral index still yields approximately $b_\ell^j$1 average normalized absolute error over $b_\ell^j$2 in the 20-channel case. The method also remains robust when sinusoidal spectral distortions are added to the synchrotron spectrum, with average error remaining approximately $b_\ell^j$3–$b_\ell^j$4 for 20 channels in the tested cases (Olivari et al., 2015).

BINGO-specific analyses pushed this application further. In simulated BINGO data, GNILC combined with a power-spectrum debiasing procedure produced a near optimal reconstruction of the HI signal using an 80 bins configuration, resulting in a power spectrum reconstruction average error over all frequencies of $b_\ell^j$5. The same study reports robustness against different synchrotron emission models and shows that adding an extra channel with CBASS foreground information reduced the estimation error of the 21 cm signal (Mericia et al., 2022).

The debiasing step is explicit in the BINGO pipeline. The reconstructed noiseless HI angular power spectrum estimate is

$b_\ell^j$6

where $b_\ell^j$7 is the additive noise bias estimated by Monte Carlo projection through fixed GNILC weights, and $b_\ell^j$8 is the multiplicative transfer function estimated from projected pure-HI realizations (Mericia et al., 2022). This indicates that practical GNILC use in HI surveys may require explicit treatment of both multiplicative signal suppression and additive projected noise.

More recent SKA-like studies report that GNILC, PCA, Need-PCA, GMCA, and Need-GMCA have comparable results, recovering the HI signal within $b_\ell^j$9 accuracy across the frequency channels in the multipole range $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$0. The same analysis states that GNILC has significantly lower foreground leakage at $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$1 compared to PCA and GMCA, owing to its dynamic estimation of the number of removed foreground components across scales and locations, and that the cleaning methods are insensitive to polarization leakage in the tested setup with a $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$2 Galactic plane mask (Caro et al., 2 Sep 2025).

6. Extensions, diagnostics, comparisons, and limitations

Several direct extensions of GNILC have been proposed. The expanded GNILC framework, eGNILC, performs the Discrete Cosine Transform along the frequency axis, calculates the eGNILC bias to modify the criterion for determining the foreground degrees of freedom, and embeds Robust Principal Component Analysis in mixing-matrix computation to obtain a blind component-separation method. In that framework, the reconstructed signal covariance is biased as

$$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$3

so the eGNILC bias depends on the foreground degrees of freedom $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$4 and the effective domain size $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$5, but not the underlying 21-cm signal (Dai et al., 2024). The same work reports that with varying Airy-disk beam the power spectra of 21-cm can be effectively recovered at multipoles $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$6 for SKA-MID and $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$7 for BINGO, with no instrumental noise, with SKA-MID exhibiting $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$8 power loss and BINGO exhibiting $$\beta_\nu^j(\hat n)=\sum_{\ell,m}(a_{\ell m,\nu} b_\ell^j)Y_{\ell m}(\hat n),$$9 power loss (Dai et al., 2024).

A modified GNILC also appears in comparison work with Needlet Karhunen–Loève cleaning. That study introduces MGNILC, which incorporates an approximation of the foregrounds to improve performance. Tested on simulated maps including polarized foregrounds, GNILC, MGNILC, GMCA, ICA, and PCA are compared with NKL. NKL is reported to provide the best performance in both tests, with a factor of $x_a^j(p)$0 to $x_a^j(p)$1 improvement over GNILC at $x_a^j(p)$2 in the higher-redshift case and $x_a^j(p)$3 in the lower-redshift case, but none of the methods were found to recover the power spectrum satisfactorily at all BAO scales (Podczerwinski et al., 2023). This indicates that GNILC’s success depends on the observing regime and on the structure of the contamination being targeted.

GNILC has also become a diagnostic tool in CMB polarization work. In optimized constrained moment ILC, a GNILC-inspired diagnosis provides the number of locally detectable foreground degrees of freedom $x_a^j(p)$4 per needlet scale and sky region. The ocMILC pipeline then sets the number of moment-deprojection constraints to $x_a^j(p)$5 and optimizes which moments to deproject, which pivot parameters to use, and how strongly to enforce each constraint (Carones et al., 2024). In that role, GNILC is not merely a separator but a local dimension estimator for foreground complexity.

A further downstream use is cosmological-parameter debiasing. In LiteBIRD-like simulations, GNILC is used to reconstruct cleaned multifrequency foreground-emission maps; these maps are combined with the weights adopted for CMB reconstruction to estimate the spatial distribution of foreground residuals after component separation. The debiased residual-template spectrum is then added to the likelihood model with a free amplitude parameter, and the reported result is unbiased estimates of the tensor-to-scalar ratio $x_a^j(p)$6 regardless of its input value, the assumed foreground model, or the adopted masking strategy (Carones, 23 Oct 2025).

The limitations emphasized across the literature are consistent. GNILC relies on nuisance or signal priors for whitening, so inaccuracies in the prior covariance can affect whitening and model-order selection (Collaboration et al., 2016). In low-SNR domains, AIC selection can be sensitive, and in the Planck application $x_a^j(p)$7 in some high-latitude, small-scale domains, reducing effective resolution and yielding spatially varying beams (Collaboration et al., 2016). In HI applications, underestimating the HI prior can degrade recovery, while in eGNILC highly correlated adjacent channels can cause numerical issues in $x_a^j(p)$8 inversion and frequency-dependent beams can induce large-$x_a^j(p)$9 errors if untreated (Caro et al., 2 Sep 2025, Dai et al., 2024). These results suggest that GNILC is best regarded as a localized, prior-augmented, multi-dimensional ILC framework whose effectiveness depends on the fidelity of the covariance model, the stability of local covariance estimation, and the angular and spectral structure of the foregrounds.