A publishing partnership

The following article is Open access

Bayesian Multi-line Intensity Mapping

Yun-Ting Cheng, Kailai Wang, Benjamin D. Wandelt, Tzu-Ching Chang, and Olivier Doré

Published 2024 August 14 • © 2024. The Author(s). Published by the American Astronomical Society.
The Astrophysical Journal, Volume 971, Number 2 Citation Yun-Ting Cheng et al 2024 ApJ 971 159 DOI 10.3847/1538-4357/ad57b9

Download Article PDF

DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

Article metrics

273 Total downloads
Video abstract views

Dates

Received 2024 March 28
Revised 2024 June 8
Accepted 2024 June 11
Published 2024 August 14

Unified Astronomy Thesaurus concepts

Large-scale structure of the universe; Cosmic background radiation; Cosmology

Journal RSS

Create or edit your corridor alerts

What are corridors?

Abstract

Line intensity mapping (LIM) has emerged as a promising tool for probing the 3D large-scale structure through the aggregate emission of spectral lines. The presence of interloper lines poses a crucial challenge in extracting the signal from the target line in LIM. In this work, we introduce a novel method for LIM analysis that simultaneously extracts line signals from multiple spectral lines, utilizing the covariance of native LIM data elements defined in the spectral–angular space. We leverage correlated information from different lines to perform joint inference on all lines simultaneously, employing a Bayesian analysis framework. We present the formalism, demonstrate our technique with a mock survey setup resembling the SPHEREx deep-field observation, and consider four spectral lines within the SPHEREx spectral coverage in the near-infrared: Hα, [O iii], Hβ, and [O ii]. We demonstrate that our method can extract the power spectrum of all four lines at the ≳10σ level at z < 2. For the brightest line, Hα, the 10σ sensitivity can be achieved out to z ∼ 3. Our technique offers a flexible framework for LIM analysis, enabling simultaneous inference of signals from multiple line emissions while accommodating diverse modeling constraints and parameterizations.

Export citation and abstract BibTeX RIS

Previous article in issue

Next article in issue

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Line intensity mapping (LIM; for reviews, see Kovetz et al. 2017; Bernal & Kovetz 2022) is an emerging technique for studying the large-scale structure (LSS) of the Universe. By mapping the emission from specific spectral lines and determining their redshifts from observed frequencies, LIM traces the three-dimensional (3D) LSS using cumulative emissions from all sources. It serves as a promising method for bridging the gap in LSS probes between the recombination era explored by the cosmic microwave background and the lower-redshift Universe (z ≲ 3) accessible with current and upcoming galaxy surveys—e.g., the Sloan Digital Sky Survey (Tegmark et al. 2006), the Dark Energy Survey (Elvin-Poole et al. 2018; Abbott et al. 2022), DESI (DESI Collaboration et al. 2016), Euclid (Laureijs et al. 2011), the Rubin Observatory (LSST Science Collaboration et al. 2009), SPHEREx (Doré et al. 2014, 2016, 2018), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015). Additionally, LIM provides crucial constraints on the collective properties of the interstellar medium (ISM) of galaxies across cosmic time through the measurement of aggregate line emission.

The field of LIM, initially pioneered by 21 cm cosmology with a primary focus on probing the Dark Ages and the Epoch of Reionization (EOR; Furlanetto et al. 2006; Morales & Wyithe 2010; Pritchard & Loeb 2012), has since expanded to include atomic or molecular lines across a broad electromagnetic spectrum. In the submillimeter wavelengths, several LIM experiments targeting [C ii] and/or CO rotational ladders have reported preliminary detections or provided upper-limit constraints, including COPSS (Keating et al. 2015, 2016), mmIME (Keating et al. 2020), and COMAP (Breysse et al. 2022; Cleary et al. 2022), with more measurements anticipated from ongoing and upcoming experiments like FYST (CCAT-Prime Collaboration et al. 2023), CONCERTO (CONCERTO Collaboration et al. 2020), TIME (Crites et al. 2014; Sun et al. 2021), SPT-SLIM (Karkare et al. 2022), EXCLAIM (Ade et al. 2020), and the Terahertz Intensity Mapper (Vieira et al. 2020). In the optical and near-infrared regimes, the ongoing HETDEX experiment (Hill et al. 2008; Gebhardt et al. 2021) is conducting Lyα LIM at z ∼ 2–4. The upcoming near-infrared all-sky spectral survey SPHEREx will explore emissions from multiple lines, including Hα, [O iii], Hβ, and [O ii], among others (Feder et al. 2023). Last, the proposed next-generation far-infrared observatories, such as PRIMA (Moullet et al. 2023), are poised to conduct LIM across various far-infrared lines ([Ne ii], H₂, [S iii], and [Si ii], etc.).

The spectral coverage of most of the LIM experiments, except for 21 cm LIM, are accessible to multiple spectral lines from different redshift ranges. While analyzing signals from multiple lines has the potential to unveil more ISM physics, as different lines often trace different ISM environments and/or dust properties, it also poses a significant data analysis challenge, since the emissions from different lines are mixed in the LIM data set.

Extensive studies have explored different strategies for tackling the challenge of line confusion in LIM. One approach involves masking voxels—fundamental 3D LIM data elements defined by pixels and spectral channels—that contain bright interlopers, using external galaxy survey catalogs to mitigate the interloper line signal (Silva et al. 2015; Yue et al. 2015; Sun et al. 2018; Béthermin et al. 2022; Van Cuyck et al. 2023). External catalogs can also be used in cross-correlation with LIM data, enabling the extraction of signals from individual lines (Pullen et al. 2013; Silva et al. 2013, 2015; Croft et al. 2016; Chung et al. 2019; Yang et al. 2019; Visbal & McQuinn 2023). Furthermore, with the LIM data set, different line signals can also be distinguished by the anisotropy of the interloper power spectrum arising from projection to the comoving frame of the target line (Cheng et al. 2016; Lidz & Taylor 2016; Gong et al. 2020).

Another approach for mitigating line confusion in LIM involves leveraging the correlation between pairs of observed frequencies containing different lines from the same redshift. This correlation arises because these frequencies trace the same underlying LSS, while the interlopers in each channel of the pair are uncorrelated, originating from distinct line-of-sight (LOS) distances. Cheng et al. (2020) utilize this information to extract an intensity map from individual lines in LIM using a spectral-template-fitting technique. In addition, previous studies have explored the possibility of probing cross-spectra between two different lines within the same or different LIM data sets (Visbal & Loeb 2010; Visbal et al. 2011; Gong et al. 2012; Serra et al. 2016; Roy & Battaglia 2024). Furthermore, some proposed methods use combinations of multiple cross-line power spectra to estimate the autospectrum of individual lines (Beane et al. 2019; Schaan & White 2021; McBride & Liu 2023). The cross-bispectrum of two lines can also achieve similar results in the same spirit (Beane & Lidz 2018).

In this study, we present a novel method for simultaneously extracting the line signal from multiple spectral lines in a LIM data set, building upon the technique introduced in (Cheng et al. 2023; hereafter, C23). In C23, we developed an inference framework to reconstruct the LSS clustering, the spectral energy distribution (SED), and the LOS distribution of emitting sources from the cross-frequency angular power spectra, ${C}_{{\ell },\nu \nu ^{\prime} }$ s. C23 demonstrated the effectiveness of this technique by applying it to multiple broadband photometric maps. C23 found that sharp features in the SED significantly enhance parameter constraints, since they break the degeneracy between spectral and redshift information in 2D frequency maps. Given this finding, LIM emerges as a promising application of the technique, as line emissions are inherently sharp features in the SED. Therefore, this study aims to extend the method developed in C23 to LIM to reconstruct the redshift evolution of line emissions. In LIM, the 3D spatial distribution of line emission, which traces the underlying LSS, is encoded in 3D spectral–angular space. On large scales, all two-point-level information in a LIM data set is contained in the data covariance, or equivalently, in harmonic space, the cross-frequency angular power spectra ${C}_{{\ell },\nu \nu ^{\prime} }$ s with all combinations of observed frequency channels ν and $\nu ^{\prime}$ . By assuming only the homogeneity and isotropy of cosmological density fluctuations, along with knowledge of the rest-frame frequencies of spectral lines, we employ a Bayesian approach to extract the bias-weighted line intensity as a function of redshift for all spectral lines in LIM from the data covariance, ${C}_{{\ell },\nu \nu ^{\prime} }$ .

In contrast to many aforementioned methods for separating spectral lines in LIM, our approach avoids the common step of projecting the LIM data into 3D comoving space at an assumed central redshift. Instead, we directly model the covariance in the native data space, i.e., in the spectral–angular space, where the data covariance ( ${C}_{{\ell },\nu \nu ^{\prime} }$ ) naturally incorporates anisotropic projection, line–line correlations, and the redshift evolution of line intensity and cosmological fluctuations within our formalism.

To demonstrate our method, we implement it in a simulated survey resembling the SPHEREx deep-field configuration with its expected sensitivity, and consider four spectral lines—Hα, [O iii], Hβ, and [O ii]—within the SPHEREx spectral coverage. We apply our algorithm to a simulated observed data covariance ( ${C}_{{\ell },\nu \nu ^{\prime} }$ ) and assess the uncertainties with our inference of the input line signals.

This paper is organized as follows. Section 2 provides a detailed description of the LIM signal and power spectra. Our assumed survey setup and the models for line signals are outlined in Sections 3 and 4, respectively. Section 5 introduces our inference algorithm. The results of applying our technique to mock observed data are presented in Section 6, followed by additional discussions in Section 7. We highlight the unique advantages of our method in Section 8 and provide discussions on comparisons with relevant previous works in Section 9. Future prospects for extending this work are discussed in Section 10, and the conclusion is provided in Section 11. Throughout this work, we assume a flat ΛCDM cosmology consistent with the measurements from Planck (Planck Collaboration et al. 2020).

2. Power Spectrum Modeling

In this section, we present the formalism for the intensity field from line emission in LIM in Section 2.1 and for the auto- and cross-frequency angular power spectra, ${C}_{{\ell }}^{\nu \nu ^{\prime} }$ , in Section 2.2. Only the main expressions are presented here; more detailed derivations are provided in Appendix A.

2.1. Intensity Field

The intensity field at an observed frequency ν and angular position $\hat{n}$ is given by the emission from continuum, spectral lines, foregrounds, and the noise:

The continuum and the foreground usually have a smooth spectrum, which can be effectively mitigated by methods that filter out the smooth spectral component in the data (see Section 7.6 for discussions).⁶ The line signal might be altered during the high-pass-filtering process in foreground cleaning, and we assume this potential bias has already been corrected in this study. However, we emphasize that any line signal transfer function induced by the foreground filtering must be carefully characterized in practice.

Therefore, in this work, we assume that the LIM data set only contains the line emission from the N_line number of spectral lines and the noise fluctuations:

In typical LIM experiments, the spectral resolution is not sufficient to resolve the intrinsic line profile from sources, and thus we model the line profile as a Dirac delta function at the rest-frame line frequency ${\nu }_{\mathrm{rf}}^{i}$ , which gives the line intensity in the following form (see Appendix A.1 for detailed derivations):

where c is the speed of light, H(z) is the Hubble parameter, and ${z}_{i\nu }={\nu }_{\mathrm{rf}}^{i}/\nu -1$ is the redshift of the ith line at the observed frequency⁷ ν. We define χ_{i
ν} as the comoving distance at redshift⁸ z_{i
ν}, and ${M}_{0,i}(\chi ,\hat{n})={{dL}}_{i}(\chi ,\hat{n})/{dV}$ is the comoving line luminosity density.

On large scales, ${M}_{0,i}(\chi ,\hat{n})$ follows the underlying matter density field ${\delta }_{m}(\chi ,\hat{n})$ with a scale-independent luminosity-weighted bias factor b_i(χ):

where M_0,i(χ) is the mean luminosity density averaged over angular positions $\hat{n}$ , and

where D_A and D_L are the comoving angular diameter distance and luminosity distance, respectively.

We model the noise at the frequency band ν as a zero-mean Gaussian fluctuation with the variance ${\sigma }_{n}^{2}(\nu )$ :

We assume the noise in each frequency channel is independent, and thus there is no cross-channel noise covariance.

2.2. Angular Power Spectrum

We describe the covariance of the LIM data set in terms of the auto- and cross-angular power spectra, ${C}_{{\ell },\nu \nu ^{\prime} }$ s, of all combinations of the frequency channels ν and $\nu ^{\prime}$ . On large scales, ignoring the redshift space distortion (RSD) effect, the line emission field is isotropic, and its fluctuations can be fully described by a Gaussian probability distribution. Therefore, the ${C}_{{\ell },\nu \nu ^{\prime} }$ s capture the full two-point information from the LIM data set on large scales (Wandelt 2013). In this work, we focus only on two-point statistics. However, we note that the LIM field on small scales is highly non-Gaussian, and the power spectrum alone is insufficient to capture the full information. This has motivated previous studies to include one-point statistics to exploit more information from LIM data (Ihle et al. 2019; Breysse 2022; Chung et al. 2023).

The total power spectrum is the sum of the contribution from emission lines and noise:

Here, the boldface C _ℓ denotes the angular-power-spectrum matrix of size N_ν × N_ν, where N_ν is the total number of spectral channels, and each element of C _ℓ is given by ${C}_{{\ell },\nu \nu ^{\prime} }$ .

In this work, we only consider information from large (linear) scales, ignoring fluctuations from nonlinear clustering and Poisson noise. We verify that, for the scales considered in this work, the Poisson noise power is negligible (Appendix A.3). Therefore, the angular power spectrum from emission lines can be expressed as

where D(χ) is the linear growth factor, and P(k) is the linear matter power spectrum at the present time. W_{i
ν}(χ) is the window function of the frequency channel ν for the emission from line i. Here, we assume the observing filter profile is a narrow top-hat function centered at the observed frequency ν, with a width Δν that spans the frequency range ${\nu }_{\min }\lt \nu \lt {\nu }_{\max }$ . With this assumption, the window function can be expressed as (see Appendix A.2 for detailed derivations)

where ${\chi }_{i\nu }^{\min /\max }=\chi ({z}_{i\nu }^{\min /\max })$ and ${z}_{i\nu }^{\min /\max }=({\nu }_{\mathrm{rf}}^{i}/{\nu }^{\max /\min })-1$ . ${z}_{i\nu }^{\min /\max }$ and ${\chi }_{i\nu }^{\min /\max }$ denote the corresponding redshift and comoving distance that can be probed by the ith line in the frequency channel ν spanning the frequency range ${\nu }_{\min }\lt \nu \lt {\nu }_{\max }$ .

As most current and upcoming LIM experiments target a relatively small field size (typically around degree scales), we adopt the Limber approximation (Limber 1953), which is valid for small survey sizes (e.g., Huterer et al. 2013). The angular power spectrum can thus be simplified as

We define

Hereafter, the term "bias-weighted luminosity" refers to M_i(χ) and, similarly, "bias-weighted intensity" refers to the quantity ${b}_{i}(\chi )\nu {I}_{\nu }^{i}(\chi )$ . We also define

Then, the angular power spectrum from the lines can be expressed as

where ${\chi }_{{ii}^{\prime} \nu \nu ^{\prime} }$ and ${\rm{\Delta }}{\chi }_{{ii}^{\prime} \nu \nu ^{\prime} }$ denote the center and the width of the overlapping comoving distance between ${\chi }_{i\nu }^{\min }\lt \chi \lt {\chi }_{i\nu }^{\max }$ and ${\chi }_{i^{\prime} \nu ^{\prime} }^{\min }\lt \chi \lt {\chi }_{i^{\prime} \nu ^{\prime} }^{\max }$ , respectively. For the autospectra of a line ( $i=i^{\prime}$ ), if there is no overlap between the filter profile, which is the case we consider in this work, the ${C}_{{\ell },\nu \nu ^{\prime} {ii}^{\prime} }^{\mathrm{line}}$ s are only nonzero for the same spectral channel (i.e., $\nu =\nu ^{\prime}$ ). For two different lines ( $i\ne i^{\prime}$ ), the ${C}_{{\ell },\nu \nu ^{\prime} {ii}^{\prime} }$ s are only nonzero when the two channels probe the two emission lines i and $i^{\prime}$ from the same redshift. Figure 3 presents an example angular-power-spectrum matrix ${{\boldsymbol{C}}}_{\ell }^{\mathrm{line}}$ , where the line signal model is detailed in Section 4.

Under the assumption that the noise is white noise without cross-channel correlation, ${{\boldsymbol{C}}}_{\ell }^{n}$ is an ℓ-independent diagonal matrix:

where δ^K is the Kronecker delta.

In practice, the data usually exhibit correlated noise from the instrument and/or foreground residuals. In Section 7.2, we present the results of applying our algorithm in the presence of such correlated noise.

Finally, there are stochastic fluctuations in the real data, such that the observed power spectrum of the data at a multipole bin ℓ, ${{\boldsymbol{C}}}_{\ell }^{d}$ , is a random sample from a Wishart distribution with a scale matrix given by C _ℓ (Equation (7)) and the degree of freedom n_ℓ, where n_ℓ is the number of ℓ modes in the binned angular power spectrum, which depends on the bin width and the survey angular size, detailed in Section 3.

2.2.1. Caveats of our Power Spectrum Model

We employ the Limber approximation in this work. However, we note that there are a few caveats associated with this simplification. First, the accuracy of the Limber approximation depends on the comoving width of the window function W_{i
ν}. For narrow widths, the LOS modes may contribute to a non-negligible level. Furthermore, if there are lines with close rest-frame frequencies, such as the [O iii] and Hβ lines in the case we considered (see Section 3), additional correlation will occur between the two lines due to the correlation between close LOS distances, even if the two lines do not fall in the same spectral channel. While this LOS correlation adds additional information for the inference, this is being ignored in our current implementation using the Limber approximation. All these effects can be properly taken into account by using the exact expression in Equation (8) instead of relying on the Limber approximation, albeit at the cost of computing triple integrations. While these are important considerations in practice, for the purpose of demonstrating our technique, we defer more detailed investigations to future work.

Here, we ignore the RSD effect in our model. The RSD effect would introduce additional terms to the angular power spectrum in Equation (8), helping to break the degeneracy between b_i(χ) and M_0,i(χ) in the window function. By using the Limber approximation in Equation (10), our formalism only accounts for the transverse modes, which are not impacted by the RSD effect. Therefore, without the inclusion of RSD, our focus is on constraining the quantity b_i(χ)M_0,i(χ) from the data rather than the two terms individually. More detailed investigations considering the RSD effect are left for future work.

3. Survey Setup

We demonstrate our technique with a survey setup similar to the deep-field survey of the SPHEREx mission⁹ (Doré et al. 2014, 2016, 2018). SPHEREx is the next NASA Medium Class Explorer mission, scheduled to launch in early 2025. SPHEREx will carry out the first all-sky near-infrared spectro-imaging survey from 0.75 to 5 μm, with a pixel size of 6 farcs 2, through four consecutive surveys over the nominal 2 yr mission. SPHEREx consists of six H2RG detector arrays spanning six broad bands in the near-infrared, with low-resolution spectroscopy conducted by linear variable filters (Korngut et al. 2018). Each band contains 17 channels with different spectral resolutions: R = 41 for bands 1–3 at wavelengths between 0.75 and 2.42 μm (51 spectral channels), R = 35 for band 4 between 2.42 and 3.82 μm (17 spectral channels), R = 110 for band 5 between 3.82 and 4.42 μm (17 spectral channels), and R = 130 for band 6 between 4.42 and 5.00 μm (17 spectral channels). SPHEREx will scan the north and south ecliptic poles with a much higher cadence, due to its scanning strategy. Consequently, SPHEREx will produce two deep-field mosaic maps of ∼100 deg² each, with the noise rms ∼50 times lower than its all-sky survey (Figure 1).

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Top: SPHEREx spectral resolution of each channel. Bottom: SPHEREx surface brightness sensitivity per spectral channel in a 62 sky pixel (blue) and the corresponding instrument noise power spectrum in the deep fields (Equation (14); orange). The solid points and dashed lines represent the SPHEREx all-sky and deep-field sensitivity, respectively. We use the 96-channel configuration from the SPHEREx public products. The gray shaded region marks the 64 channels considered in this work (0.75–3.82 μm).
Download figure:
Standard image High-resolution image

farcs — **Figure 1.** Top: SPHEREx spectral resolution of each channel. Bottom: SPHEREx surface brightness sensitivity per spectral channel in a 62 sky pixel (blue) and the corresponding instrument noise power spectrum in the deep fields (Equation (14); orange). The solid points and dashed lines represent the SPHEREx all-sky and deep-field sensitivity, respectively. We use the 96-channel configuration from the SPHEREx public products. The gray shaded region marks the 64 channels considered in this work (0.75–3.82 μm).
Download figure:
Standard image High-resolution image

Here, we consider a survey setup similar to SPHEREx deep fields totaling 200 deg² (100 deg² in each ecliptic pole) in this study, corresponding to a sky fraction of f_sky = 0.48%.¹⁰ We use the first four bands of SPHEREx spanning 0.75–3.82 μm and assume nonoverlapping top-hat filters equally spaced in logarithmic frequency. The last two bands are not included, since they only probe the very high redshift emission from the four lines we consider in the SPHEREx spectral coverage (see Table 1). We consider an angular resolution of 6 farcs 2 and the surface brightness sensitivity in each channel given by the public products.¹¹ The SPHEREx public products are based on the previous design, which has 16 instead of 17 spectral channels in each band, and thus we also use the same configuration of 16 channels per band in this work, which gives 64 channels (four bands) in total.

Table 1. Spectral Lines Modeled in This Work

Line	λ_rf	SPHEREx Coverage	r_i	A_i
Hα	0.6563	0.14 < z < 4.82(6.62)	1.27	1.0
[O iii]	0.5007	0.50 < z < 6.63(8.99)	1.32	1.32
Hβ	0.4861	0.54 < z < 6.68(9.29)	0.44	1.38
[O ii]	0.3727	1.01 < z < 9.25(12.4)	0.71	0.62

Note. λ_rf—rest-frame wavelength (μm); the maximum redshifts, with and without parentheses, correspond to the first four bands and the full SPHEREx coverage with λ = 3.82 μm and 5 μm, respectively; r_i—line luminosity and the SFR ratio L_i/SFR (10⁴¹ erg s⁻¹ ${M}_{\odot }^{-1}$ yr); and A_i—dust extinction factor (mag).

Download table as: ASCII Typeset image

The top panel of Figure 1 shows the SPHEREx spectral resolution as a function of wavelength (spectral channel), and the bottom panel shows the expected noise rms per spectral channel per pixel in SPHEREx and the corresponding noise power spectrum given by Equation (14).

We consider the following four lines within the SPHEREx spectral coverage: Hα, [O iii], Hβ, and [O ii]. Table 1 summarizes their rest-frame wavelengths and the redshift ranges that SPHEREx can probe. For the purpose of demonstrating our technique, we only focus on these four lines in this work. We note that there are more spectral lines within the SPHEREx spectral range—such as Lyα, Paschen-α, [N ii], and [S ii]—and the line flux of some of them may be comparable to the four lines considered here (Feder et al. 2023). Therefore, in practice, the analysis for SPHEREx should account for all the prominent lines for a more realistic modeling.

We consider the line emission from sources within the redshift range 0.7 < z < 6. Removing detected local point sources at lower redshift helps improve the sensitivity in probing the diffuse line emission from fainter sources. We estimate that below our chosen redshift lower limit ${z}_{\min }=0.7$ , we can reliably detect and constrain the redshifts of the majority of galaxies with SPHEREx, thus allowing them to be masked before calculating the power spectrum (see Appendix B for details). The choice of the maximum redshift ${z}_{\max }$ does not affect our results; as shown in Figure 6, we have no sensitivity on the line signal at z ≳ 4.

We choose the multipole mode range of 50 < ℓ < 350 in our analysis. The minimum ℓ mode ( ${{\ell }}_{\min }=50$ ) corresponds to SPHEREx's field of view of 3 fdg 5 (on the smaller side), as fluctuations larger than the field size will be partially suppressed due to zodiacal light filtering in processing individual exposures and we do not model this effect. The choice of the maximum ℓ mode ( ${{\ell }}_{\max }=350$ ) is made to restrict our analysis to linear clustering scale (see Appendix B for details).

We use eight ℓ bins within the range of 50 < ℓ < 350. The bins are selected to contain approximately the same number of modes in each bin. For a given ℓ bin spanning ${\ell }\in [{{\ell }}_{\min }^{\alpha },{{\ell }}_{\max }^{\alpha })$ , the number of multipole modes ${n}_{{\ell }}^{\alpha }$ is given by

where f_sky = 0.48% is the fraction of sky area in the 200 deg² field. In reality, some pixels with bright foreground contamination will be masked, resulting in a reduction of the effective number of modes. We ignore this effect in our analysis.

4. Line Signal Modeling

Our model for line emission follows the prescription from Gong et al. (2017), which is built on an empirical star formation rate (SFR) and the line luminosity relation. We use the SFR density (SFRD) constraints to model the luminosity density for each line.

We assume a linear relation between the line luminosity L_i and the SFR:

and we use the L_i–SFR relations from Kennicutt (1998) and Ly et al. (2007) for Hα, [O ii], and [O iii]. For Hβ, we assume a fixed line ratio of L_Hβ/L_Hα = 0.35 (Osterbrock & Ferland 2006), which has been validated to have good agreement with simulations and observations by Gong et al. (2017). The line luminosity–SFR ratio (r_i) for each line is summarized in Table 1. Following Gong et al. (2017), we adopt the same dust extinction factors, also listed in Table 1.

Despite this model adopting a simple linear scaling for the L_i–SFR relation, Gong et al. (2017) have validated that the resulting line intensity is in agreement with another model based on simulations as well as the observational constraints from integrating the observed line luminosity functions. We also note that there are scatters in the L_i–SFR relation in reality, which will boost the power spectrum amplitude (Sun et al. 2019). For simplicity, we ignore the effect of this scatter in this work.

For the SFRD, we use the analytical fitting formula from Madau & Dickinson (2014):

The luminosity density of the line is then given by a linear scaling of the SFRD:

We model the luminosity-weighted bias b_i of the lines with the halo-mass-weighted bias, under the assumption that the line luminosity is proportional to the halo mass:

where $\tfrac{{dn}}{{dM}}$ is the halo mass function (Sheth & Tormen 1999) and b_h is the halo bias (Sheth et al. 2001). Our bias prescription assumes the linear relation between the line luminosity and the halo mass. We assume the same b_i(z) for all spectral lines. Note that although here we build models for M_0,i(z) and b_i(z) separately, our method can only constrain M_i(z), the product of these two terms.

Figure 2 shows our model of the bias-weighted luminosity density M_i(z) and the bias-weighted intensity ${b}_{i}(z)\nu {I}_{\nu }^{i}(z)$ of each line (see Equation (3) for the conversion from luminosity density to intensity). Since our line modeling is a linear scaling of the SFRD, M_i(z) follows the same redshift dependence as the SFRD, which exhibits a peak at z ∼ 2. The M_i(z) from different lines are linearly proportional to each other in our model, since we assume the same line bias and a linear scaling from SFRD for all lines.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Top: our model of the bias-weighted luminosity density M_i(z) as a function of redshift. Middle: the bias-weighted line intensity as a function of redshift. The gray shaded regions in the top and middle panels denote the range of redshift considered in this work (0.7 < z < 6). Bottom: the bias-weighted line intensity as a function of observed wavelength. The solid parts in all three panels denote the redshift/wavelength ranges that can be probed by the spectral range considered in this work (0.75–3.82 μm). The SPHEREx noise level σ_n is also shown at the bottom panel for comparison (black). We note that while noise fluctuations overwhelm the line intensity, the large-scale clustering power of the lines is not suppressed by the noise power, as shown in Figure 3.
Download figure:
Standard image High-resolution image

While this simple model is not able to capture the complex line emission mechanisms in reality, it is sufficient for our purpose of demonstrating the algorithm in this study. Furthermore, as detailed in Section 5.1, we introduce a flexible parameterization to fit for any redshift dependence of the M_i(z) function, which is not restricted to a certain functional form of M_i(z).

The top panel of Figure 3 presents the line power spectrum matrix ${{\boldsymbol{C}}}_{\ell }^{\mathrm{line}}$ from our model at the lowest-ℓ bin centered at ℓ = 91. The diagonal power is from the autocorrelation of each line at the same redshift/frequency. The "broadened" band along the diagonal line is the cross-power of the closely paired [O iii] and Hβ lines. Other off-diagonal correlations arise from different combinations of cross-power between pairs of lines, as marked in Figure 3. Coincidentally, the cross-power of Hα and the [O iii]–Hβ pair falls at overlapping frequency channels with [O ii] and the [O iii]–Hβ pair.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Top: the line power spectrum matrix ${{\bf{C}}}_{{\ell }}^{\mathrm{line}}$ (Equation (23)) in the lowest-ℓ mode centered at ℓ = 91. The signals from different pairs of lines are labeled in the figure. Bottom: the auto power spectrum of each line (colored lines) compared to the SPHEREx noise power spectrum (black dashed line). The short-wavelength cutoff corresponds to the minimum redshift ${z}_{\min }=0.7$ , where we assume no line signal below this redshift, as galaxies can effectively be masked given the SPHEREx deep-field depth (Appendix B).
Download figure:
Standard image High-resolution image

**Figure 3.** Top: the line power spectrum matrix ${{\bf{C}}}_{{\ell }}^{\mathrm{line}}$ (Equation (23)) in the lowest-ℓ mode centered at ℓ = 91. The signals from different pairs of lines are labeled in the figure. Bottom: the auto power spectrum of each line (colored lines) compared to the SPHEREx noise power spectrum (black dashed line). The short-wavelength cutoff corresponds to the minimum redshift ${z}_{\min }=0.7$ , where we assume no line signal below this redshift, as galaxies can effectively be masked given the SPHEREx deep-field depth (Appendix B).
Download figure:
Standard image High-resolution image

The bottom panel of Figure 3 shows the autospectrum of each line compared to the SPHEREx noise power spectrum. Although the SPHEREx noise overwhelms the mean intensity of all the line signals (Figure 2), the line emission traces the underlying density field, exhibiting large-scale clustering that boosts the line power spectrum above the noise fluctuations on large scales (Figure 3). Moreover, we emphasize that in the power spectrum space, from which we extract information, the noise is present only in the diagonal elements. The unique off-diagonal features in the frequency–frequency correlation space are unaffected by noise fluctuations.

5. Algorithm

This section describes our algorithm for constraining the line signals from the LIM data. Our method infers the bias-weighted luminosity density (M_i(z)) for each line from the auto- and cross-frequency power spectra ( ${C}_{{\ell },\nu \nu ^{\prime} }$ s). We first introduce our parameterization for M_i(z) (Section 5.1.1), then we describe our Bayesian inference framework (Section 5.2) and the algorithm for inferring the parameter constraints (Section 5.3).

5.1. Parameterization

Our goal is to infer the function M_i(z) for each line from the LIM data. This is achieved by first parameterizing M_i(z) and then fitting the defined parameters to constrain the M_i(z) functions. M_i(z) can be characterized with a small number of parameters, as M_i(z) is expected to be a smooth function with redshift. One approach is to use a parametric functional form to define a smooth curve for M_i(z). However, here we instead choose a series of linear basis models for M_i(z). As detailed in Section 5.1.1, this parameterization allows us to precompute the basis power spectra, ${\hat{C}}_{{\ell },\nu \nu }$ s, to significantly reduce the computational cost during the inference stage. To ensure flexibility in capturing any possible redshift dependence of the signal, we employ a basis function set that forms a piecewise linear function for M_i(z) (Section 5.1.2). The piecewise linear function can approximate any continuous functions, and thus we are not restricted to any prior assumptions about the shape of M_i(z) functions.

5.1.1. Linear Basis Decomposition

To parameterize M_i(z), we decompose it with a linear combination of basis functions $\{{\hat{M}}_{j}(z)\}$ :

and we fit for the coefficients {c_ij} to constrain M_i(z) for each line. With this decomposition, we can also express the power spectrum in terms of the linear combinations of the basis functions:

where

The total power spectrum from all the spectral lines can then be written as

In the parameter inference stage (Section 5.3), we also need the derivatives of the power spectrum with respect to the parameters. This can also be expressed as the linear combination of the basis components:

With this linear basis decomposition, the basis power spectra ${\hat{{\boldsymbol{C}}}}_{{\ell },{ii}^{\prime} {jj}^{\prime} }$ can be precomputed to greatly speed up the inference process.

5.1.2. Basis Functions

We use a piecewise linear function to describe the bias-weighted luminosity density M_i(z). The piecewise linear function can flexibly approximate any continuous function with a sufficiently fine segmentation, and it does not depend on any underlying assumption about the shape of the function being fitted. Since M_i(z) follows the global redshift dependence of the large-scale bias and the SFRD, it is expected to be a smooth function of redshift. Thus, only a few segments of the piecewise linear function are sufficient to approximate the M_i(z) functions.

The piecewise linear function can be expressed in terms of linear combinations of a series of rectified linear unit (ReLU) functions. Thus, we define our basis functions $\{{\hat{M}}_{j}(z)\}$ as ReLU functions with anchoring redshifts z_js:

We choose {z_j} = {−1, 0,..., 5}, which gives a total number of basis functions N_m = 7. The basis functions ${\hat{M}}_{j}(z)$ are shown in the top panel of Figure 4.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Top: the seven ReLU basis functions $\{{\hat{M}}_{j}(z)\}$ for decomposing the bias-weighted luminosity density (M_i(z)) for each line. Bottom: our fiducial input model M_i(z) (the blue solid line with dots). This is a piecewise linear function with anchoring redshifts at z = 0, 1,..., 6 (blue dots) that fit to the modeled M_i(z) (yellow curve; the same as the top panel of Figure 2). This piecewise linear function can be produced by the linear combinations of the seven basis functions $\{{\hat{M}}_{j}(z)\}$ shown in the top panel. Here, we show the Hα line as an example. The fiducial model for other lines is set with the same process. The gray shaded region denotes our redshift range of 0.7 < z < 6.
Download figure:
Standard image High-resolution image

The linear combination of this basis set spans the piecewise linear functions anchored at z = z_j + 1 (i.e., z = 0, 1,..., 6). This provides the flexibility to approximate any underlying M_i(z) functions, and we can also easily increase the accuracy at any specific redshift range by adding more ReLU basis functions with z_j around the desired redshifts. The optimal number of redshift anchoring points and their positions depend on the line signals and noise of particular surveys. Therefore, we leave this further investigation to future work.

While {c_ij} is the native parameter set in our formalism, the underlying line signal is best described by {m_ij}, defined as

which is the M_i value at the anchoring redshifts z = z_j + 1. There is a simple linear transformation between {m_ij} and {c_ij}, detailed in Appendix C. This transformation relation enables us not only to convert values between {m_ij} and {c_ij}, but also to propagate the parameter constraints, which are determined by the Jacobian of this transformation.

5.1.3. Fiducial Parameters

Our fiducial input parameters {c_ij} are set by matching their corresponding {m_ij} to the modeled M_i(z) (the top panel of Figure 2) at the seven anchoring redshifts (z = 0, 1,..., 6), as demonstrated in the bottom panel of Figure 4.

In reality, our piecewise linear model serves as an approximate representation of the true signal. Nonetheless, we employ this piecewise linear approximation as the fiducial input, providing a set of ground-truth parameters for evaluating our algorithm's performance.

We emphasize that while our fiducial input assumes the same shape for M_i(z) across all four lines, our algorithm fits each line separately. This allows us to reconstruct the redshift evolution of M_i(z) for each line independently. Further demonstration of this capability is provided in Section 7.5, where different M_i(z) functions are used to generate the mock signal for each line, and we validate that our algorithm can robustly extract the inputs, even when the input M_i(z)s are smooth curves that are not able to be perfectly described by our piecewise linear parameterization. The example in Section 7.5 also assumes very different shapes of M_i(z) for each line, to demonstrate the algorithm's robustness against model variations.

5.2. Bayesian Framework

Our parameter inference method follows the framework presented in C23. We constrain the parameter set Θ from the data power spectra $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ in N_ℓ multipole bins using a Bayesian framework. The posterior probability distribution $p\left({\boldsymbol{\Theta }}| \{{{\boldsymbol{C}}}_{{\ell }}^{d}\}\right)$ is given by

where ${ \mathcal L }$ and π are the likelihood and prior, respectively.

Here, our parameter set Θ consists of the coefficients of the basis components for M_i(z) for each line, i.e., {c_ij}. For the prior, we only enforce the positivity condition on M_i(z), which is equivalent to requiring all m_ij values to be positive in our piecewise linear model. Here, m_ij represents the M_i(z) function at the anchoring redshifts in our piecewise linear model, thus enforcing the positivity of m_ij guarantees that M_i(z) is positive across all redshifts. We implement this positivity constraint through a logarithmic transformation on m_ij, defining a new set of parameters {θ_ij}, where ${\theta }_{{ij}}=\mathrm{log}{m}_{{ij}}$ (see Appendix C for details). Our objective is to determine the maximum a posteriori solution for {θ_ij}. By doing so, the positivity constraint on m_ijs will be automatically satisfied. We set flat priors on {θ_ij}, which effectively give the logarithmic priors on {m_ij}.¹²

As we consider the two-point information in this work, the likelihood of the full LIM data set, i.e., the voxel intensities, can be described as a Gaussian distribution, and the cross angular power spectrum ${C}_{{\ell },\nu \nu ^{\prime} }$ s represent the covariance matrices of the Gaussian likelihood on the voxel intensity maps in the spherical harmonic space. As each multipole mode is independent, the log-likelihood function is the sum of normal distributions ${ \mathcal N }$ for each ℓ bin:

where n_ℓ is the number of modes in each ℓ bin (Equation (15)), and ${{\boldsymbol{C}}}_{{\ell }}\left({\boldsymbol{\Theta }}\right)$ is the modeled power spectrum given the parameter set Θ.

We note that our likelihood models the voxel intensity, the native data product from LIM, as a Gaussian distribution with covariances given by the ${C}_{{\ell },\nu \nu ^{\prime} }$ s. This captures the full two-point information in the field, which is a lossless representation on large scales, where the underlying signal is expected to be fully characterized by a Gaussian distribution.

5.3. Parameter Inference

With the observed angular power spectra from the LIM data $\{{{\boldsymbol{C}}}_{\ell }^{d}\}$ , we conduct parameter inference within the Bayesian framework outlined in Section 5.2. The inference process is similar to our prior work in C23, involving two steps: first, we employ the Newton–Raphson method to identify the parameter set corresponding to the maximum likelihood; then, we estimate the parameter constraints using the Fisher matrix derived from the maximum likelihood found with the Newton–Raphson method.

5.3.1. Newton–Raphson Method

The Newton–Raphson method is an iterative approach to find the maximum/minimum of a function. In this context, our objective is to find the parameter set ${{\rm{\Theta }}}_{\max }$ that maximizes the log-likelihood, given the angular power spectra from the data, $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ :

The Newton–Raphson algorithm iteratively updates the current parameter set from Θ_t to Θ_t+1 through the equation:

where ${\bf{g}}={\rm{\nabla }}\,\mathrm{log}\,{ \mathcal L }$ and ${\boldsymbol{H}}={\rm{\nabla }}(\mathrm{log}\,{ \mathcal L }){{\rm{\nabla }}}^{T}$ represent the gradient and Hessian matrix of the log-likelihood, respectively. For the complete expressions and detailed derivations, see Appendix E of C23. The parameter η serves as the learning rate, determining the step size of the update. In each iteration of the Newton–Raphson update, we initiate with η = 2 and subsequently verify whether the new proposed parameter set yields a higher log-likelihood value. If not, we reduce the value of η by half until the condition is satisfied.

The Newton–Raphson optimization is performed with the logarithmically transformed parameter set {θ_ij}. In this parameter space, the positivity condition for the prior is automatically satisfied, eliminating the need for additional prior constraints, such as a hard boundary for the disallowed parameter space.

Similar to the algorithm in C23, we utilize an approximated Hessian provided by Equation E33 in C23 to avoid computing the second derivatives of C_ℓ on parameters, thus expediting the optimization process. This approximation approaches the exact expression when Θ_t is in proximity to ${{\boldsymbol{\Theta }}}_{\max }$ , and we have confirmed the successful convergence of our algorithm using this approximation.

5.3.2. Fisher Matrix

After determining ${{\boldsymbol{\Theta }}}_{\max }$ using the Newton–Raphson method, we estimate the parameter constraints with the Fisher matrix at ${{\boldsymbol{\Theta }}}_{\max }$ . The Fisher matrix is given by

and the inverse of the Fisher matrix gives the covariance of the parameters:

The Fisher matrix calculation requires the derivative of C _ℓs. The derivatives can also be expressed as linear combinations of the basis functions given by Equation (24).

In Section 6, we also quantify the constraints on M_i(z) at any given redshift. From Equation (20), we can obtain the covariance of M_i(z) and ${M}_{i}^{\prime} (z^{\prime} )$ for the two given lines i and $i^{\prime}$ at redshift z and $z^{\prime}$ by the expression

6. Results

Here, we present the results of inference using our fiducial model with the SPHEREx deep-field noise level. The mock observed power spectra $\{{{\boldsymbol{C}}}_{\ell }^{d}\}$ are produced by the fiducial parameter set described in Section 5.1.3. That is, we set M_i(z)s to piecewise linear functions anchored at redshifts z = 0, 1,...,6, and the value of M_i(z) for each line is fixed to an analytical model described in Section 4. For each of the four lines, we fit for seven linear coefficients {c_ij} that define the M_i(z) function. In summary, our data consist of the 64 × 64 (64 spectral channels) symmetric matrices of $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ in eight ℓ bins, and we fit for 28 parameters (seven parameters for each of the four lines) in total.

With the input spectra $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ , we first use the Newton–Raphson method to find the parameter set ${{\boldsymbol{\Theta }}}_{\max }$ that maximizes the log-likelihood function (Section 5.3.1). The Newton–Raphson optimization in our case can efficiently converge within a few tens of steps. Then, we calculate the covariance on parameters using the Fisher matrix (Section 5.3.2).

Figure 5 displays the inference results on the fiducial model. We run the inference on cases with and without adding sample variance fluctuations to the data ( $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ ), respectively. The case with sample variance fluctuations (orange contours) represents a realistic scenario, and the results show that our algorithm gives parameter constraints within about a 1σ level of the truth, as expected. As a sanity check, we also run the case without sample variance fluctuations in the input data (blue contours). In this case, the likelihood function peaks at exactly the truth input values, verifying that our Newton–Raphson method can successfully locate the maximum a posteriori.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Parameter constraints for our fiducial case. We display only the c_ij constraints for the two basis components that are sensitive to the lowest-redshift line emission accessible by our survey setup for each line. The blue/orange colors represent the mock observed data power spectra $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ without/with sample variance fluctuations. The black dashed lines indicate the truth model input, while the blue/orange dots represent the maximum likelihood found by our Newton–Raphson algorithm. The contours show the 1σ and 2σ constraints derived from the Fisher matrix.
Download figure:
Standard image High-resolution image

**Figure 5.** Parameter constraints for our fiducial case. We display only the c_ij constraints for the two basis components that are sensitive to the lowest-redshift line emission accessible by our survey setup for each line. The blue/orange colors represent the mock observed data power spectra $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ without/with sample variance fluctuations. The black dashed lines indicate the truth model input, while the blue/orange dots represent the maximum likelihood found by our Newton–Raphson algorithm. The contours show the 1σ and 2σ constraints derived from the Fisher matrix.
Download figure:
Standard image High-resolution image

Next, we propagate our parameter constraints into the bias-weighted intensity for each line. We use Equation (33) to first derive the constraints on M_i(z), and we convert the M_i(z) to intensity with Equation (3). The results are shown in Figure 6. With our fiducial setup similar to the SPHEREx deep-field sensitivity, our algorithm can extract the bias-weighted intensity for all four lines at the ≳10σ level at z < 2. For the brightest line, Hα, the 10σ sensitivity can be achieved out to z ∼ 3, with the signal-to-noise ratio (S/N) peaking at ∼100σ around z = 1.5. While in reality, the sensitivity depends on both the underlying signal model and the systematic uncertainties not accounted for in our forecast, our results suggest a promising prospect for simultaneously detecting LIM signals from multiple lines with SPHEREx, by leveraging the information encoded in the correlation between lines in the spectral–angular space.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Top: our constraints on the bias-weighted intensity for each line from the fiducial case. The solid/dashed lines denote the input model within the redshift ranges accessible/inaccessible by the survey spectral coverage. We consider the line signal only from 0.7 < z < 6 (see Section 3). The colored shaded regions mark the 1σ constraints from our inference with the case that contains sample variance in the input mock data (the orange case in Figure 5). Second panel: the S/N on the bias-weighted b(z) · ν I_ν(z) for each line (colored lines). The solid, dashed, and dotted black lines mark the 1σ, 3σ, and 10σ sensitivity levels, respectively. Bottom four panels: 1σ constraints for each line relative to the truth input (Hα, [O iii], Hβ, and [O ii], from top to bottom, respectively).
Download figure:
Standard image High-resolution image

We can also estimate the constraining power on the SFRD using our constraints on the bias-weighted intensity. Here, we propagate our Hα constraint to the SFRD, assuming that the bias-weighted intensity of Hα is proportional to the SFRD and ignoring uncertainties in this conversion factor. In other words, we assume the same S/N for the bias-weighted intensity of Hα and the SFRD. The results are shown in Figure 7. We see that with the SPHEREx survey setup and applying our algorithm for the LIM analysis, we can derive a competitive constraining power on the SFRD at z ≲ 3. While our algorithm simultaneously infers multiple lines that trace the star formation history, allowing for a potentially tighter constraint on the SFRD through the joint information from all lines, this analysis will require additional modeling of the correlation between lines, which is beyond the scope of this work.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The 1σ SFRD constraints from our fiducial case. We assume that the SFRD is proportional to the bias-weighted intensity of the Hα line, and thus the sensitivity on the SFRD is propagated from our Hα constraints, as shown in Figure 6 (red shaded region). The red curve denotes our underlying model for the SFRD from Madau & Dickinson (2014). For comparison, we also display observational constraints from Cucciati et al. (2012; green) and Reddy & Steidel (2009; purple).
Download figure:
Standard image High-resolution image

7. Discussion

7.1. Dependence on the Noise Level

Our fiducial case, presented in Section 6, considers the SPHEREx deep-field sensitivity. Here, we explore how the model constraints depend on the noise level. We use the same model for the line signal as in the fiducial case, apply different scaling to the fiducial noise variance in the SPHEREx deep field, ${\left({\sigma }_{n}^{\mathrm{SPHEREx}-\mathrm{deep}}\right)}^{2}$ , across all frequency channels, and use the Fisher matrix to derive parameter constraints as a function of the noise level.

The results are displayed in Figure 8. We observe that the sensitivity increases as the noise level decreases and approaches the noiseless limit (plus symbols) for bright lines at lower redshifts, where the line power is significantly stronger than the noise.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** The S/N on the bias-weighted intensity at z = 1 (top), 2 (middle), and 3 (bottom) for each line is shown as a function of the noise level relative to the fiducial SPHEREx deep-field noise, ${\sigma }_{n}^{\mathrm{SPHEREx}-\mathrm{deep}}$ . The absence of [O ii] at z = 1 is because it falls outside the range that can be probed by our survey. The plus markers denote the noiseless case (σ_n = 0). The solid, dashed, and dotted gray horizontal lines mark the 1σ, 3σ, and 10σ sensitivity levels, respectively. The solid and dotted gray vertical lines mark ${\sigma }_{n}^{2}/{\left({\sigma }_{n}^{\mathrm{SPHEREx}-\mathrm{deep}}\right)}^{2}\,=\,1$ and 0.5, respectively. The gray dashed vertical line marks the SPHEREx all-sky survey sensitivity. The colored crosses scale the all-sky S/N by the square root of the sky coverage ratio to represent the constraining power in the all-sky survey.
Download figure:
Standard image High-resolution image

**Figure 8.** The S/N on the bias-weighted intensity at z = 1 (top), 2 (middle), and 3 (bottom) for each line is shown as a function of the noise level relative to the fiducial SPHEREx deep-field noise, ${\sigma }_{n}^{\mathrm{SPHEREx}-\mathrm{deep}}$ . The absence of [O ii] at z = 1 is because it falls outside the range that can be probed by our survey. The plus markers denote the noiseless case (σ_n = 0). The solid, dashed, and dotted gray horizontal lines mark the 1σ, 3σ, and 10σ sensitivity levels, respectively. The solid and dotted gray vertical lines mark ${\sigma }_{n}^{2}/{\left({\sigma }_{n}^{\mathrm{SPHEREx}-\mathrm{deep}}\right)}^{2}\,=\,1$ and 0.5, respectively. The gray dashed vertical line marks the SPHEREx all-sky survey sensitivity. The colored crosses scale the all-sky S/N by the square root of the sky coverage ratio to represent the constraining power in the all-sky survey.
Download figure:
Standard image High-resolution image

Figure 8 also provides insights into the sensitivity of detecting the line signal if SPHEREx extends beyond its nominal 2 yr survey. For instance, if SPHEREx extends its mission lifetime to 4 yr, the noise variance will integrate down by a factor of 2, as indicated by the vertical dotted line in Figure 8. In this case, we find about a factor of 2 sensitivity improvement for detecting the line signals at z = 3, whereas at z = 1, the improvement is less evident, especially for brighter lines such as Hα, as the signals are already not in the noise-dominated regime, even in the case of the fiducial sensitivity of the nominal 2 yr mission.

Similarly, we can also estimate the sensitivity if we apply our algorithm to the SPHEREx all-sky survey instead of in deep fields. The SPHEREx all-sky noise variance is about 50 times higher than the deep field (Figure 1), which is denoted by the dashed vertical lines in Figure 8. While the curves in Figure 8 indicate low sensitivity at this noise level, we note that we will have access to a much larger sky coverage in the all-sky survey. Fisher information on the parameters is proportional to the number of available modes and thus the sky coverage. Assuming SPHEREx all-sky coverage of ${f}_{\mathrm{sky}}^{\mathrm{SPHEREx}-\mathrm{all}}\approx 75 \%$ , we get a factor of ∼150 more sky coverage than the deep field (f_sky = 0.48%), and thus a ∼12 times boost of the S/N, indicated as the colored crosses in Figure 8. We find that at z = 1, the all-sky and deep-field sensitivity is similar, while for higher redshift, at z = 3, the deep field will perform better than the all-sky survey, as the higher-redshift line signals are fainter and thus more susceptible to noise fluctuations.

We emphasize that we cannot achieve infinite sensitivity to the line signal, even in the absence of instrument noise (the plus symbols in Figure 8). This fundamental sensitivity limit is due to the nature of the line confusion in the data. In the spectral–angular space, emission from different lines is mixed together, and hence acts as "line noise" for each other.

7.2. Presence of Correlated Noise

In this work, we assume that instrument noise is uncorrelated across channels, making the noise power spectra a diagonal matrix. In reality, the data usually contain correlated noise from the intrinsic instrumental noise, foreground residuals, and the continuum-filtering process.

To assess how correlated noise might affect the reconstruction results, we create a test case with correlated noise, as shown in the top panel of Figure 9. We assume continuous correlated noise that decays with channel separation to emulate typical foreground residuals. Additionally, we add a few off-diagonal streaks that manifest similar features to line correlation signals. These kinds of noise features can be caused by, for example, detector crosstalk. We set the eigenvalues of this correlated noise matrix to be similar to our fiducial case.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Top: an example of a correlated noise matrix. Bottom: comparison of the S/N on M_i(z) for each line with uncorrelated (solid) and correlated (dashed) noise. We assume the same signal model and survey setup as the fiducial case, using the noise matrix shown in the top panel for the correlated noise and a diagonal matrix with the same eigenvalues for the uncorrelated noise.
Download figure:
Standard image High-resolution image

We then quantify the S/N on M_i(z) for each line with a Fisher forecast (Figure 9, bottom panel). We compare the case of using the correlated noise shown in the top panel with the case of uncorrelated noise with the same eigenvalues. We find the parameter constraints to be almost identical in both cases.¹³ This indicates that our algorithm is not strongly affected by the presence of correlated noise. This is because the correlated line signal has a certain (and deterministic) pattern in the cross-power-spectrum space, so only correlated noise that manifests the exact pattern of the spectral correlation of the signal will induce significant degeneracy in our inference.

7.3. Information from Small Scales

Our fiducial setup limits the analysis to large scales (ℓ < 350), to model the signal from only the linear regime. This discards a huge amount of smaller-scale modes accessible by SPHEREx. To assess the information content from these higher-ℓ modes, we calculate the S/N on line reconstruction with a setup that extends the maximum multipole mode to ${{\ell }}_{\max }=2000$ , while keeping other assumptions the same as the fiducial model. The Poisson noise becomes non-negligible on smaller scales; therefore, in this calculation, we include a model of Poisson noise (Appendix A.3) for both cases. The results are shown in Figure 10. While there are ∼30 times more modes between 350 < ℓ < 2000 compared to our fiducial setup of 50 < ℓ < 350, the higher-multipole modes correspond to higher-k modes in the matter power spectrum P(k) with a lower power and thus are more susceptible to the (white) instrument noise, and thus they do not contain as much information as the lower-ℓ modes. Therefore, there is only a factor of ∼1.6 gain in the line sensitivity by the inclusion of the higher-ℓ modes. We also note that in reality, to extract information from small scales, one must include the effect of nonlinear clustering in the model, which we have ignored in this calculation.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Top: the S/N of line reconstruction in the fiducial setup (solid) and the case with ${{\ell }}_{\max }$ extending from 350 to 2000 (dashed). Bottom: the ratio of the inferred line sensitivity of the extended ${{\ell }}_{\max }$ case and the fiducial case.
Download figure:
Standard image High-resolution image

**Figure 10.** Top: the S/N of line reconstruction in the fiducial setup (solid) and the case with ${{\ell }}_{\max }$ extending from 350 to 2000 (dashed). Bottom: the ratio of the inferred line sensitivity of the extended ${{\ell }}_{\max }$ case and the fiducial case.
Download figure:
Standard image High-resolution image

Similarly, we can estimate the information loss if higher-multipole modes have to be discarded due to, for example, higher nonlinear bias of the line emission field compared to our current model. Figure 11 compares the S/N when reducing ${{\ell }}_{\max }$ from 350 to 200. While this corresponds to discarding 70% of multipole modes, the sensitivity on M_i(z) is only reduced by ∼30%, for the same reason as above: the large-scale modes are less noisy than the small-scale modes.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Top: the S/N of line reconstruction in the fiducial setup (solid) and the case with ${{\ell }}_{\max }$ reducing from 350 to 200 (dashed). Bottom: the ratio of the inferred line sensitivity of the extended ${{\ell }}_{\max }$ case and the fiducial case.
Download figure:
Standard image High-resolution image

7.4. Capability of Interloper Separation

Our method jointly constrains signals from multiple spectral lines in LIM data sets. Here, we highlight the capability of our method in extracting faint line signals from multiple interlopers with orders of magnitude stronger power. As the most commonly used summary statistic in LIM is the 3D power spectrum P(k), we present our model constraints in this representation. Specifically, we show the first few multipole moments of P(k), which capture the anisotropy of the line power spectrum. This anisotropy is due to the incorrect projection of interloper lines to the target line redshift. In LIM analysis, the ith interloper signal will be projected from the redshift z_i to the target line redshift z_t, where z_i = λ_t(1 + z_t)/λ_i − 1, making the interloper and the target line fall in the same observed frequency. The different projection in the transverse and the LOS directions makes the interloper power spectrum anisotropic (Cheng et al. 2016; Lidz & Taylor 2016; Gong et al. 2020), which introduces nonzero ℓ > 0 modes when expanding the 3D power spectrum with Legendre polynomials (Bernal et al. 2019; 2021; Gong et al. 2020).

The ℓth multipole of the total LIM power spectrum is the sum of the contribution from all lines:

The P_ℓ from the ith line is given by

where ${{ \mathcal L }}_{{\ell }}$ is the Legendre polynomial, ${q}_{\perp }^{i}$ and ${q}_{\parallel }^{i}$ are the projection factors in the transverse and the LOS directions, respectively, and ^Pi(k_i, μ_i, z_i) is the intrinsic power spectrum of the ith line at redshift z_i and the corresponding Fourier mode k_i and cosine angle μ_i. For the intrinsic power spectrum, ^Pi(k_i, μ_i, z_i), we also incorporate the RSD effect (Kaiser 1987) and the window functions due to the finite resolution in the LOS and transverse directions. Both effects introduce additional sources of anisotropy to the observed power spectrum. The full expression of these quantities is presented in Appendix D.

Figure 12 presents the power spectrum multipoles ${P}_{{\ell }}^{i}(k,{z}_{t})$ for the three lowest-ℓ modes.¹⁴ Here, we choose [O ii] as the target line and present the results at the target line redshift z_t = 2.5, corresponding to the observed wavelength at 1.3 μm, and the interlopers are from redshifts z_i = 0.99 (Hα), 1.61 ([O iii]), 1.68 (Hβ).

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Top: the 3D power spectrum multipoles P_ℓ(k) of the lines projected to the [O iii] frame at z = 2.5. We show the first three multipole modes at ℓ = 0 (left), ℓ = 2 (middle), and ℓ = 4 (right). Bottom: the ratio of P_ℓ(k) from each line to the target line [O iii]. The shaded region on the [O iii] line (green) denotes the 1σ constraint of the bias-weighted intensity (b_i(z)ν I_ν(z)) from our fiducial case (Section 6), which gives an S/N of 3.9 at z = 2.5.
Download figure:
Standard image High-resolution image

From Figure 12, we see that the [O ii] power spectrum is overwhelmed by interloper power by more than 2 orders of magnitude in all three multipole modes. Nevertheless, our method makes use of the information on the frequency correlations between lines and can successfully extract the bias-weighted intensity of [O ii] at z_t = 2.5 with an S/N of 3.9.

7.5. Robustness against Model Misspecification

All calculations up to this point have utilized the input line signal from our fiducial model (Section 5.1.3), constructed from our piecewise linear parameterization. While the piecewise linear function offers great flexibility to approximate any continuous function, with a limited number of anchoring points, it cannot perfectly capture realistic bias-weighted luminosity density functions, which are expected to be smooth curves. To validate that our algorithm can faithfully reconstruct signals with a different underlying input, we perform the following test. We generate the input mock data $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ using Gaussian functions for the bias-weighted luminosity densities M_i(z) for each line, as illustrated in Figure 13. In addition to testing the robustness of the reconstruction under our piecewise linear approximation, we also intentionally assign very different shapes of M_i(z) for each line compared to the fiducial case, to assess whether our algorithm can still reconstruct the input accurately.

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** The biased-weighted luminosity density for the four spectral lines in our test case for model misspecification (solid lines). The dashed lines are our fiducial case for comparison. The fiducial case assumes a piecewise linear function with anchoring redshifts at z = 0,1,..., 6 (points on dashed lines).
Download figure:
Standard image High-resolution image

Then, we apply our algorithm to infer the line signal from the data using our parameterization, i.e., approximating the input M_i(z) with the piecewise linear function. The results are shown in Figure 14. We find that even if our model cannot perfectly describe the true signal, our inference can still unbiasedly reconstruct the signals within a ∼1σ range at z < 4. Beyond this redshift, however, we obtain biased results for the faintest line, [O ii], at z ∼ 5. We checked that this type of bias depends on the specific realization of the sample variance. To further investigate this feature, we ran the inference on the same data power spectra $\{{{\boldsymbol{C}}}_{{\ell }}^{d}\}$ with doubled redshift resolution, i.e., setting the redshift anchoring point spacing to Δz = 0.5 instead of the fiducial case of Δz = 1. The results are shown in Figure 15. We find a higher best-fit likelihood value, and the inference on M_i(z) is less biased than the fiducial case. This indicates that, in practice, it is essential to optimize the trade-off between the complexity and flexibility of our model in determining the number of redshift parameters. We leave more detailed analysis of this optimization to future works.

Figure 14. Refer to the following caption and surrounding text. — **Figure 14.** Top: our constraints on the bias-weighted intensity for each line, where the line signal in the input mock data is shown in Figure 13, and we reconstruct the signal with our piecewise linear parameterization. The solid/dashed lines denote the input model within the redshift ranges accessible/inaccessible by the survey spectral coverage. The colored shaded regions mark the 1σ constraints from our inference. The input mock data contain sample variance fluctuations. Second panel: the S/N on b(z) · ν I_ν(z) for each line (colored lines). The solid, dashed, and dotted black lines mark the 1σ, 3σ, and 10σ sensitivity levels, respectively. Bottom four panels: 1σ constraints for each line relative to the truth input (Hα, [O iii], Hβ, and [O ii], from top to bottom, respectively).
Download figure:
Standard image High-resolution image

Figure 15. Refer to the following caption and surrounding text. — **Figure 15.** Reconstruction results of the same input signal as in Figure 14, with finer redshift anchoring points of Δz = 0.5 instead of the fiducial case of Δz = 1 (Figure 14).
Download figure:
Standard image High-resolution image

7.6. Implementation with Continuum Foregrounds

Our model ignores continuum foregrounds that will be present in LIM data in practice. The continuum includes both Galactic and extragalactic emission. The continuum usually has a smooth spectrum, which we assume is already being filtered out in the map space before computing the cross-power spectra ${C}_{{\ell },\nu \nu ^{\prime} }$ . There are extensive studies on continuum foreground removal strategies for LIM—for example, principal component analysis (de Oliveira-Costa et al. 2008; Chang et al. 2010; Masui et al. 2013; Switzer et al. 2013; Van Cuyck et al. 2023) and the asymmetric re-weighted penalized least-squares (Van Cuyck et al. 2023) and semi-blind component separation techniques based on independent component analysis (Zhang et al. 2016).

Part of the extragalactic continuum is emitted by the same galaxies that emit the line signal. Thus, a joint framework that simultaneously fits for lines and continuum can provide additional constraints on both galaxy physics and the underlying LSS. We defer this investigation to future studies.

8. Advantages of Our Method

8.1. Multiline Inference across Redshifts

Our method infers the bias-weighted line intensity as a function of redshift from multiple lines. This differs from many previously proposed LIM analysis techniques that operate in the 3D power spectrum (P(k)) space, approximating the line signal from different frequencies as originating from the same redshift (detailed in Section 9.1). The exceptions are some works that investigate the "antisymmetric cross power spectrum" between the CO and H i lines during the EOR (Sato-Polito et al. 2020; Zhou et al. 2021). These studies show that the cross-correlation of the two lines at slightly offset LOS distances is not symmetric under exchange, yielding a nonzero antisymmetric cross-spectrum that contains additional information to constrain the EOR parameters. We emphasize that while we adopt the Limber approximation, which ignores any LOS correlation in this work, the full expression of the ${C}_{{\ell },\nu \nu ^{\prime} }$ matrices will incorporate this information, as the redshift evolution model has been encoded in our formalism.

Additionally, some line deconfusion techniques treat interlopers as nuisance emission and seek strategies to mitigate interloper contamination in order to detect the target line. Here, our analysis is not restricted to extracting only one target line, allowing us to jointly reconstruct the signal from multiple spectral lines.

The intensity from multiple lines can provide a wealth of information. For example, many atomic or molecular lines serve as good tracers for the star formation history, and combining constraints from multiple lines can enhance the understanding of the global SFRD evolution. The ratio between the line signals also offers a valuable astrophysical census. For instance, given that the intrinsic ratio between Hα and Hβ luminosity is fixed by atomic physics (L_Hβ/L_Hα = 0.35), deviations from this ratio probe the dust attenuation law. However, extracting these astrophysical constraints requires breaking the degeneracy between bias and intensity, either by using ancillary information or by jointly modeling the bias and intensity based on more detailed simulations.

8.2. Straightforward Implementation

Our model uses the auto-/cross-spectra between spectral channels, which is the covariance of the native form of the LIM data product. Many instrumental effects and foregrounds are naturally being described in this spectral–angular space (e.g., the filter transmission profile, noise correlations, and atmospheric emission). Although, in this work, we disregard these components in our analysis, our technique provides a more compatible framework for incorporating these effects in reality.

In contrast, many other LIM analysis methods convert the data into the comoving space of the target line and compute the 3D power spectrum P(k). This conversion relies on the assumed cosmological model and, furthermore, systematics may arise from errors in interpolation and projection (Cunnington & Wolz 2023).

8.3. Flexibility

Our flexible framework can accommodate any form of the signal model. In this study, we choose a piecewise linear model to fit the bias-weighted luminosity density M_i(z), providing good flexibility in approximating various functional forms. However, any parameterization for the signal can be implemented within our framework, as long as we can express the likelihood's dependence on the parameters. Additionally, incorporating any prior assumptions is straightforward, either through designing the parameterization of M_i(z) or encoding them into the Bayesian prior. While we only apply a positivity prior on M_i(z) in this work, one can introduce other prior information, such as a prior on the line ratio between certain pairs of lines.

Moreover, despite fixing the cosmological model in this analysis, one can also simultaneously fit for cosmological information and the line signal. This approach was demonstrated in C23, where we jointly fit for the matter power spectrum P(k) and the spectral and redshift dependence of the emission.

8.4. Generalizability

While we demonstrate signal reconstruction using cross-channel correlations within the same LIM data set, our framework can be extended to include correlations with other data sets. For instance, we can cross-correlate with other LIM surveys and photometric/spectroscopic galaxy catalogs to derive joint constraints from multiple probes.

9. Comparison with Other LIM Analysis Methods

In this section, we summarize the main differences between our technique and a few other LIM analysis methods in the literature.

9.1. 3D Cross-power Spectrum

Several analyses have investigated the detectability of cross-correlation between different lines, either within the same LIM data set or with different experiments that probe lines within the same cosmic volume (Visbal & Loeb 2010; Visbal et al. 2011; Gong et al. 2012; Serra et al. 2016; Roy & Battaglia 2024). Multiple line–line cross-spectra can also be used to reconstruct the autospectra of the target line (Beane et al. 2019; Schaan & White 2021; McBride & Liu 2023). These analyses consider cross-correlation on the 3D power spectrum P(k), requiring the projection of LIM data into comoving space and assuming the projected line signal at a fixed redshift. In contrast, our method does not involve projection before computing the cross-power spectrum, allowing us to model the redshift evolution of the line signal instead of assuming the emissions are from a single redshift.

The 3D power spectrum provides a simpler basis for extracting information from the LOS modes, which are useful for breaking the bias and intensity degeneracy from the RSD effect. The LOS correlations of the underlying density fluctuations are ignored in this work, due to the assumption of a top-hat window function and the Limber approximation. Without these simplifications, our formalism could also extract the LOS correlation. However, modeling LOS modes in the angular correlation space requires a double-Bessel-function integration (Equation (8)), whereas in the 3D power spectrum space, the signal is simply a Fourier transform of the field. Therefore, our frequency angular correlation and the 3D power spectrum are suitable for different analysis purposes and will be complementary to each other in practice.

9.1.1. Angular Power Spectrum Covariance

Feng et al. (2019) also employ a Bayesian framework to analyze the auto-/cross-frequency power spectrum in LIM in the context of component separation for the cosmic near-infrared background. However, they approximate the likelihood on the ${C}_{{\ell },\nu \nu ^{\prime} }$ s as a Gaussian distribution, which is only valid in the limit of high S/N and with a large number of modes (Hamimeche & Lewis 2008). In contrast, in our analysis, our data vector is the native LIM data product—the voxel intensity map (in the spherical harmonic space), with a covariance matrix given by ${C}_{{\ell },\nu \nu ^{\prime} }$ , which can fully capture the information from the data at the two-point level.

Furthermore, the analysis in Feng et al. (2019) fits the observed power spectra with a set of amplitudes associated with predefined theoretical templates, making it more susceptible to model misspecification. In contrast, our framework employs flexible parameterization, allowing it to accommodate signals that might not be well described by existing models.

9.1.2. Power Spectrum Anisotropy

The anisotropy of the 3D power spectrum P(k) of interloper lines, upon projection to the target line redshift (Section 7.4), has been proposed as a strategy for line separation in LIM studies (Cheng et al. 2016; Lidz & Taylor 2016; Gong et al. 2020). However, this technique encounters difficulties with interloper lines that have rest-frame frequencies very close to that of the target line (such as Hβ and [O iii]), due to the minimal effect of projection. In contrast, our method, which utilizes the complete information available in the spectral–angular space, is capable of successfully extracting the signals from all lines, even in cases where there is a closely located interloper.

9.1.3. Pixel-space Spectral Template Fitting

Cheng et al. (2020) introduce a technique to extract the intensity map from individual lines in LIM, also using cross-frequency information. By fitting the spectrum in each LIM pixel to a large dictionary of spectra that encode which frequency channel contains line emission for a given redshift, they can successfully infer the redshift of the line emitters that are brighter than the noise level and thus reconstruct the line intensity map from those sources.

Our method, while not capable of reconstructing the individual line signal in the map space, makes use of the emission from all sources, rather than only detectable sources. Therefore, it is not restricted to the relatively low-noise and low-confusion regime, as in Cheng et al. (2020). Furthermore, our angular power spectrum utilizes both the spectral correlations from multiple lines and the spatial correlation of the underlying cosmological field, whereas the pixel-by-pixel fitting in Cheng et al. (2020) does not involve spatial clustering information.

9.1.4. Machine Learning

Machine learning (ML) has shown promise in decomposing emissions from different lines in LIM (Moriwaki et al. 2020; Moriwaki & Yoshida 2021), leveraging information beyond the two-point statistics that is not captured in the power spectrum. However, the effectiveness of ML approaches hinges on the availability of a comprehensive training set that covers the full range of potential signals and systematic variations. This requirement is challenging for many current LIM experiments, given our limited understanding of the underlying signal models, foreground, and instrument systematics.

10. Future Work

Here, we outline directions toward better realism and broader applications of our technique for future studies.

Our previous work of C23 (continuum) establishes the inference algorithm for constraining continuum emission with broadband photometry from the frequency–frequency cross-correlation, while this study applies a similar framework for the line emission. A joint inference with both continuum and line emission will exploit the information from the full SEDs of galaxies.

Furthermore, while we only perform inference on cross-frequency correlations ( ${C}_{{\ell },\nu \nu ^{\prime} }$ ), our framework can also be extended to incorporate cross-correlation with other tracers, such as photometric or spectroscopic galaxies (Cheng & Chang 2022).

Finally, Our analysis framework can also serve as a tool to search for unknown spectral features that trace the LSS (Y.-T. Cheng et al. 2024, in preparation). This is of great interest in searching for dark matter candidates decaying into photons that may manifest as unexpected lines in the LIM data (Creque-Sarbinowski & Kamionkowski 2018; Bernal et al. 2021).

11. Conclusion

In this work, we introduce a novel technique for analyzing LIM data with multiple line signals. While the presence of multiple line signals in LIM poses a challenging analysis issue, known as interloper contamination, our method leverages the correlated information from multiple lines to perform joint inference on all lines simultaneously. This is enabled by the correlated signal from lines originating from the same redshift, manifesting as unique off-diagonal signals in the covariance of the LIM data ( ${C}_{{\ell },\nu \nu ^{\prime} }$ s).

We employ Bayesian analysis to infer the bias-weighted intensity of each line from the data covariance ${C}_{{\ell },\nu \nu ^{\prime} }$ s. Without relying on any external data set, and only making use of assumptions of the signal homogeneity and isotropy, as well as the positivity condition on the bias-weighted intensity of all lines, our method enables the full exploitation of the information in the data on large scales, where the line emission field is Gaussian and can be fully characterized by two-point statistics.

We apply our method to mock LIM power spectra generated from a survey setup similar to the SPHEREx deep-field observation, considering four lines within the SPHEREx spectral coverage: Hα, [O iii], Hβ, and [O ii]. We demonstrate that our algorithm can constrain the bias-weighted intensity of all four lines at the ≳10σ level at z < 2. For the brightest line, Hα, the 10σ sensitivity can be achieved out to z ∼ 3, with S/N peaks at ∼100σ around z = 1.5. We also show that our method is robust against model misspecification.

This work lays the foundation for broader applications in analyzing LIM data with various spectral features, which is timely, as many LIM experiments are expected to come online in the near future.

Acknowledgments

We would like to thank the anonymous referee for valuable comments that improved the manuscript. We are grateful to Asantha Cooray, Richard Feder, Adam Lidz, Jordan Mirocha, and Mike Zemcov for constructive discussions and comments on a draft manuscript, as well as Ari Cukierman, Brandon Hensley, the SPHEREx science team, and the participants in the workshop "Present and Future of Line Intensity Mapping," held at the Max Planck Institute for Astrophysics, Munich, for helpful discussions regarding this work. Y.-T.C. acknowledges the Balzan Cosmological Studies Travel Grant and the hospitality of Institut d'Astrophysique de Paris, where part of this work was conducted. K.W. acknowledges the support of the JPL SURF program. Y.-T.C. acknowledges support by NASA ROSES grant 18-2ADAP18-0192. B.D.W. acknowledges support by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche; the INFINITY NEXT project grant under the DIM ORIGINES Equipements 2023 program of the Île-de-France region; and the Simons Collaboration on "Learning the Universe." The Flatiron Institute is supported by the Simons Foundation. T.-C.C. acknowledges support by NASA ROSES grant 21-ADAP21-0122. Part of this work was done at Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We acknowledge support from the SPHEREx project under a contract from the NASA/GODDARD Space Flight Center to the California Institute of Technology.

Software: astropy (Astropy Collaboration et al. 2013, 2018), COLOSSUS (Diemer 2018), ChainConsumer (Hinton 2016).

Appendix A: Line Intensity and Power Spectrum Derivations

Here, we present a detailed derivation of the line intensity field and the window function of the power spectrum, as discussed in Sections 2.1 and 2.2.

A.1. Line Intensity Field

The intensity from the ith line at the angular position $\hat{n}$ and observed frequency ν is given by

Here, D_L is the luminosity distance, D_A is the comoving angular diameter distance (equal to the comoving distance in a flat Universe), ν_rf = ν(1 + z) is the rest-frame frequency corresponding to the observed frequency ν at redshift z, and ${L}_{\nu }^{i}({\nu }_{\mathrm{rf}},\chi ,\hat{n})={{dL}}^{i}({\nu }_{\mathrm{rf}},\chi ,\hat{n})/d{\nu }_{\mathrm{rf}}$ is the line profile (specific luminosity density) of the ith line at the angular position $\hat{n}$ and rest-frame frequency ν_rf. Note that ${\nu }_{\mathrm{rf}}{L}_{\nu }^{i}({\nu }_{\mathrm{rf}},\chi ,\hat{n})\,=\nu {L}_{\nu }^{i}(\nu ,\chi ,\hat{n})$ .

The spectral resolution of typical LIM experiments cannot resolve the intrinsic line profile from sources. Therefore, we approximate ${L}_{\nu }^{i}({\nu }_{\mathrm{rf}},\chi ,\hat{n})$ as a Dirac delta function δ^D at the rest-frame line frequency ${\nu }_{\mathrm{rf}}^{i}$ . We define the comoving line luminosity density:

With this, we get

where ${z}_{i\nu }={\nu }_{\mathrm{rf}}^{i}/\nu -1$ is the redshift of the ith line at the observed frequency ν, and χ_{i
ν} is its corresponding comoving distance. Here, c is the speed of light, and H(z) is the Hubble parameter.

Inserting this into Equation (A1), we obtain

where we define

A.2. Window Function

The intensity field of the ith line in an observed filter with the filter profile function R(ν) is given by

In this study, we consider the channel centered at a given observed frequency ν to be a top-hat function spanning ${\nu }^{\min }\lt \nu \lt {\nu }^{\max }$ . Thus, the filter profile function is

where ${\rm{\Delta }}\nu ={\nu }^{\max }-{\nu }^{\min }$ is the filter width, and the indicator function ${ \mathcal I }(\nu ;{\nu }^{\min },{\nu }^{\max })$ is defined by

The window function of the power spectrum at the channel centered at ν, denoted as W_{i
ν}(χ) in Equation (10), relates to the line intensity field by the following expression:

where we define W_{0,i
ν}(χ) = W_{i
ν}(χ)/b_i(χ). To obtain the window function for the power spectrum W_{i
ν}(χ), we first derive the expression of the intensity field in terms of the χ integration.

Inserting Equations (A4) and (A6) into Equation (A7), we get

where ${\chi }_{i\nu }^{\min /\max }=\chi ({z}_{i\nu }^{\min /\max })$ and ${z}_{i\nu }^{\min /\max }=({\nu }_{\mathrm{rf}}^{i}/{\nu }^{\min /\max })\,-1$ . Therefore, the window function of the ith line at the channel ν is

A.3. Poisson Noise

The Poisson noise power from lines i and $i^{\prime}$ in channels ν and $\nu ^{\prime}$ is given by

where dn/dM is the halo mass function (Sheth & Tormen 1999). Here, we assume the line luminosities Lⁱ are functions of the halo mass M and ignore scatters in this relation. The shot noise power depends on the details of the Lⁱ–M relation. Many previous studies have modeled this relation, ranging from scaling relations to semi-analytical models (e.g., Li et al. 2016; Fonseca et al. 2017; Moradinezhad Dizgah et al. 2022). Here, for simplicity, we assume a linear relation for Lⁱ–M, which allows us to relate the last integral in Equation (A12) to the luminosity density M_0,i by

Figure 16 compares the clustering and Poisson noise of the line power spectrum, as well as the SPHEREx noise level, in our fiducial case at our lowest- and highest-multipole bins. For our chosen scales and redshift range, the Poisson noise is much lower than the clustering signal and the noise. We also checked that the inclusion of Poisson noise has negligible effects on inference, and thus we ignore Poisson noise in our analysis.

Figure 16. Refer to the following caption and surrounding text. — **Figure 16.** The clustering (black) and Poisson noise (brown) auto power spectrum from all four lines. The colored dotted lines denote the clustering terms of individual lines. The gray dashed line shows the SPHEREx noise power spectrum for reference. We show the power spectrum in the lowest- (top) and highest- (bottom) multipole modes in our fiducial setup.
Download figure:
Standard image High-resolution image

Appendix B: Redshift and Multipole Ranges

Our choice of ${z}_{\min }=0.7$ in Section 3 is based on the expected point-source sensitivity depth of SPHEREx. We estimate that below this redshift, we can reliably detect and constrain the redshift of the majority of galaxies with SPHEREx, and thus mask them to improve the sensitivity on diffuse line emission from fainter sources.

SPHEREx is expected to achieve a 5σ point-source sensitivity of m ∼ 21.6 per channel (at λ < 3.82 μm) in the deep fields¹⁵ (the blue solid line in the top panel of Figure 17). Considering that the SPHEREx photometric redshift fitting can achieve high accuracy for galaxies with ≳3σ per channel sensitivity, this corresponds to a masking depth of ∼22.7 at 2 μm (the orange solid line in the top panel of Figure 17). Using a model of the galaxy luminosity function across redshift and wavelength from Helgason et al. (2012), we estimate that below z ∼ 0.7, there is ≲10% of the integrated galaxy emission from sources below the masking depth (the bottom panel of Figure 17), and thus we ignore any galaxy emission at z < 0.7 in our model.

Figure 17. Refer to the following caption and surrounding text. — **Figure 17.** Top: the 5σ SPHEREx point-source sensitivity per channel on all-sky (blue dashed) and in deep fields (blue solid) in AB magnitudes. The orange line denotes the 3σ sensitivity per channel in deep fields, which is the depth we assumed for the point-source masking limit. Bottom: the fraction of total galaxy intensity below m = 22.7 at 2 μm as a function of redshift.
Download figure:
Standard image High-resolution image

The choice of the maximum-ℓ mode ( ${{\ell }}_{\max }=350$ ) is made to restrict our analysis to linear clustering scales. At ${z}_{\min }=0.7$ , the effects of nonlinear clustering enter at k ≳ 0.2h Mpc⁻¹ (Figure 18), and thus we set ${{\ell }}_{\max }$ to correspond to a transverse comoving maximum-k mode of 0.2h Mpc⁻¹ at ${z}_{\min }\,=0.7$ ( ${k}_{\max }\sim {{\ell }}_{\max }/\chi ({z}_{\min })\sim 0.2$ ).

Figure 18. Refer to the following caption and surrounding text. — **Figure 18.** Top: linear (black dashed) and nonlinear (blue) matter power spectrum at z = 0.7. The nonlinear power spectrum is obtained with the `syren-halofit` package (Bartlett et al. 2024b, 2024a). Bottom: the fractional difference between the linear and nonlinear power spectrum. The angular multipole mode ℓ, corresponding to the transverse k mode, is marked at the top axis. At z = 0.7, the nonlinear clustering power deviates from the linear power by ≳10% at k ∼ 0.2, corresponding to ℓ ∼ 350. Therefore, we choose ${{\ell }}_{\max }=350$ for our analysis.
Download figure:
Standard image High-resolution image

Appendix C: Parameter Transformation

In this section, we present the relationship between the three sets of parameters: {c_ij}, {m_ij}, and {θ_ij}. We will first describe the transformation between these parameters and then present the formalism for expressing the likelihood gradient and the Fisher matrix on one set of parameters in terms of another set of parameters.

The relationship between {c_ij} and {m_ij} is

where {z_j} = {−1, 0,...,5} and $\{{\hat{M}}_{j}(z)\}$ are the ReLU functions defined in Equation (25). Representing c_ij and m_ij as the N_m-sized vectors c _i and m _i, respectively, these two vectors follow a linear transformation relation defined by the Jacobian matrix ${{\boldsymbol{J}}}_{i}^{{cm}}$ :

and

where

and

The parameter set θ_ij is defined by the logarithm of m_ij:

and thus the Jacobian ${{\boldsymbol{J}}}_{i}^{\theta m}$ is a diagonal matrix, with the diagonal elements given by m _i:

It is important to note that this transformation is nonlinear, and unlike ${{\boldsymbol{J}}}_{i}^{{cm}}$ , the Jacobian ${{\boldsymbol{J}}}_{i}^{\theta m}$ depends on the specific parameter values.

During parameter inference, we concatenate parameters for all lines (is) into a single vector. In this case, the corresponding Jacobian matrix J is a block-diagonal matrix, with each block representing the Jacobian for each line ( J _i).

With the Jacobian matrix between the parameter sets, we can transform the likelihood derivatives and Fisher matrices between different parameter sets using the following relations:

and

where α and β represent any two of the parameter sets ({c_ij}, {m_ij}, and {θ_ij}). Note that the Jacobian matrices follow the chain rule: J ^{α
β} = J ^{α
γ} J ^{γ
β}. Therefore, for example, we can obtain J ^{c
θ} by J ^{c
θ} = J ^cm J ^{m
θ}.

In our algorithm, we initially compute the likelihood derivative and Fisher matrix in {c_ij}, then transform them to {θ_ij} during the Newton–Raphson optimization. Subsequently, we perform another transformation to {m_ij} to quantify the model constraints.

Appendix D: 3D Power Spectrum Multipoles

Here, we provide a detailed expression for the interloper power spectrum multipoles, following the prescription from Bernal et al. (2021).

The ℓth multipole of the 3D power spectrum from interloper line i projected to the target line t is given by (Equation (35)):

where z_i = λ_t(1 + z_t)/λ_i − 1 is the redshift of the ith line that contaminates the target line signal at the same observed frequency. ${{ \mathcal L }}_{{\ell }}$ is the Legendre polynomial, and ${q}_{\perp }^{i}$ and ${q}_{\parallel }^{i}$ are the transverse and LOS projecting factors from the interloper redshift z_i to the target line redshift z_t, respectively, given by

where D_A is the comoving angular diameter distance, and H(z) is the Hubble parameter. The intrinsic line power spectrum is given by

Here, we only consider the clustering power on large scales, where the power spectrum is proportional to the square of the bias-weighted intensity b_i(z)ν I_ν,i(z_i). D(z) is the linear growth rate, and P(k) is the linear power spectrum at the present day. The first term introduces intrinsic anisotropy from the RSD effect (Kaiser 1987), where ${f}_{i}\approx {{\rm{\Omega }}}_{m}^{0.55}({z}_{i})$ is the linear growth rate at z_i. The corresponding Fourier mode and cosine angle of the interloper, k_i and μ_i, respectively, are given by

where ${F}_{\mathrm{proj}}^{i}={q}_{\parallel }^{i}/{q}_{\perp }^{i}$ .

The window function W(k, μ, z) is given by

where

and R is the spectral resolution, D_A is the comoving angular diameter distance, and σ_beam is the beam size that we take for the SPHEREx pixel size of 6 farcs 2.

We note that a more comprehensive 3D power spectrum modeling requires accounting for the redshift evolution of signals along the LOS, which is beyond the scope of this study.

Footnotes

6
We note that some foregrounds, like Galactic and atmospheric foregrounds, may contain spectral line emission or absorption. Thus, in practice, additional data-processing steps may be necessary to remove these foreground line features before implementing our analysis.
7
Throughout this manuscript, ν is referred to as the observed frequency, and the rest-frame frequency is denoted by ν_rf.
8
The redshift z and the comoving distance χ will be used interchangeably to describe the LOS distance.
9
http://spherex.caltech.edu
10
SPHEREx will also provide all-sky coverage with shallower depth, which can potentially be used for LIM analysis. However, the lower redundancy over all sky will make reliable all-sky maps harder to construct. Therefore, we focus on the deep field in this study.
11
Data downloaded from: https://github.com/SPHEREx/Public-products/blob/master/Surface_Brightness_v28_base_cbe.txt.
12
In practice, one can incorporate external information, such as line luminosity function constraints, into priors. The choice of different priors may have a non-negligible impact on inference (Millea & Bouchet 2018). We leave these considerations for future work.
13
We note that marginalized constraints on individual parameters depend nontrivially on the noise covariance, and therefore the S/N for certain parameters in the presence of correlated noise could be higher than in the noncorrelated case.
14
Since the line power spectrum ^Pi(k_i, μ_i, z_i) is an even function of the LOS cosine angle μ, the odd multipole modes vanish.
15
https://github.com/SPHEREx/Public-products/blob/master/Point_Source_Sensitivity_v28_base_cbe.txt

Please wait… references are loading.

Bayesian Multi-line Intensity Mapping

Article metrics

Share this article

Dates

Abstract

1. Introduction

2. Power Spectrum Modeling

2.1. Intensity Field

2.2. Angular Power Spectrum

2.2.1. Caveats of our Power Spectrum Model

3. Survey Setup

4. Line Signal Modeling

5. Algorithm

5.1. Parameterization

5.1.1. Linear Basis Decomposition

5.1.2. Basis Functions

5.1.3. Fiducial Parameters

5.2. Bayesian Framework

5.3. Parameter Inference

5.3.1. Newton–Raphson Method

5.3.2. Fisher Matrix

6. Results

7. Discussion

7.1. Dependence on the Noise Level

7.2. Presence of Correlated Noise

7.3. Information from Small Scales

7.4. Capability of Interloper Separation

7.5. Robustness against Model Misspecification

7.6. Implementation with Continuum Foregrounds

8. Advantages of Our Method

8.1. Multiline Inference across Redshifts

8.2. Straightforward Implementation

8.3. Flexibility

8.4. Generalizability

9. Comparison with Other LIM Analysis Methods

9.1. 3D Cross-power Spectrum

9.1.1. Angular Power Spectrum Covariance

9.1.2. Power Spectrum Anisotropy

9.1.3. Pixel-space Spectral Template Fitting

9.1.4. Machine Learning

10. Future Work

11. Conclusion

Acknowledgments

Appendix A: Line Intensity and Power Spectrum Derivations

A.1. Line Intensity Field

A.2. Window Function

A.3. Poisson Noise

Appendix B: Redshift and Multipole Ranges

Appendix C: Parameter Transformation

Appendix D: 3D Power Spectrum Multipoles

Footnotes