<figure><img src="image_1.png" /></figure>

Bayesian Multi-line Intensity Mapping

Yun-Ting Cheng¹,² ID, Kailai Wang³ ID, Benjamin D. Wandelt⁴,⁵ ID, Tzu-Ching Chang¹,² ID, and Olivier Doré¹,² ID

¹California Institute of Technology, 1200 E. California Boulevard, Pasadena, CA 91125, USA; ycheng3@caltech.edu

²Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

³Department of Physics, Cornell University, 616 Thurston Avenue, Ithaca, NY 14853, USA

⁴Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

⁵Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

Received 2024 March 28; revised 2024 June 8; accepted 2024 June 11; published 2024 August 14

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.

*Unified Astronomy Thesaurus concepts: Large-scale structure of the universe (902); Cosmic background radiation (317); Cosmology (343)*

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 electro-
magnetic 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;

<figure><img src="image_2.png" /></figure>

distribution of this work must maintain attribution to the author(s) and the title
of the work, journal citation and DOI.

Original content from this work may be used under the terms
of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further
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,\ell'\nu\nu'}$. 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,\ell'\nu\nu'}$ with all combinations of observed frequency channels $\nu$ and $\nu'$. 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,\ell'\nu\nu'}$.

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,\ell'\nu\nu'}$) 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$\alpha$, [O III], H$\beta$, and [O II]—within the SPHEREx spectral coverage. We apply our algorithm to a simulated observed data covariance ($C_{\ell,\ell'\nu\nu'}$) 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 $\Lambda$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\ell'\nu\nu'}$, 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 $\nu$ and angular position $\hat{n}$ is given by the emission from continuum, spectral lines, foregrounds, and the noise:

$$ \nu I_{\nu}(\nu, \hat{n}) = \nu I_{\nu}^{\text{cont}}(\nu, \hat{n}) + \nu I_{\nu}^{\text{line}}(\nu, \hat{n}) + \nu I_{\nu}^{\text{FG}}(\nu, \hat{n}) + \nu I_{\nu}^{\text{n}}(\nu, \hat{n}). \quad (1) $$

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).<fnref n="6" /> 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_{\text{line}}$ number of spectral lines and the noise fluctuations:

$$ \nu I_{\nu}(\nu, \hat{n}) = \nu I_{\nu}^{\text{line}}(\nu, \hat{n}) + \nu I_{\nu}^{\text{n}}(\nu, \hat{n}) \\ = \sum_{i=1}^{N_{\text{line}}} \nu I_{\nu}^{i}(\nu, \hat{n}) + \nu I_{\nu}^{n}(\nu, \hat{n}). \quad (2) $$

<footnote n="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.</footnote>
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'_{\text{rf}}$, which gives the line intensity in the following form (see Appendix A.1 for detailed derivations):

$$ \nu I_{\nu}^{i}(\nu, \hat{n}) = \frac{c(1 + z_{i\nu})}{H(z_{i\nu})} M_{0,i}(\chi_{i\nu}, \hat{n}) A_{0}(\chi_{i\nu}), \quad (3) $$

where c is the speed of light, $H(z)$ is the Hubble parameter, and $z_{i\nu} = \nu'_{\text{rf}}/\nu - 1$ is the redshift of the ith line at the observed frequency$^7$ $\nu$. We define $\chi_{i\nu}$ as the comoving distance at redshift$^8$ $z_{i\nu}$, 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(\chi)$:

$$ M_{0,i}(\chi, \hat{n}) = b_i(\chi) M_{0,i}(\chi) [1 + \delta_m(\chi, \hat{n})], \quad (4) $$

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

$$ A_0(\chi) = D_A^2(\chi)/4\pi D_L^2(\chi), \quad (5) $$

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

We model the noise at the frequency band $\nu$ as a zero-mean Gaussian fluctuation with the variance $\sigma_n^2(\nu)$:

$$ \nu I_{\nu}^{n}(\nu, \hat{n}) \sim \mathcal{N}(0, \sigma_{n}^{2}(\nu)). \quad (6) $$

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,\ell'\nu'\nu}$'s, of all combinations of the frequency channels $\nu$ and $\nu'$. 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,\ell'\nu'\nu}$'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:

$$ C_{\ell} = C_{\ell}^{\text{line}} + C_{\ell}^{n} = \sum_{i=1}^{N_{\text{line}}} \sum_{i'=1}^{N_{\text{line}}} C_{\ell,ii'} + C_{\ell}^{n}. \quad (7) $$

Here, the boldface $C_{\ell}$ denotes the angular-power-spectrum matrix of size $N_{\nu} \times N_{\nu'}$, where $N_{\nu'}$ is the total number of spectral channels, and each element of $C_{\ell}$ is given by $C_{\ell,\ell'\nu'}$.

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

$$ \begin{aligned} C_{\ell,\nu\nu'} &= \sum_{i=1}^{N_{\text{line}}} \sum_{i'=1}^{N_{\text{line}}} C_{\ell,\ell'\nu'\nu'} \\ &= \sum_{i=1}^{N_{\text{line}}} \sum_{i'=1}^{N_{\text{line}}} \int d\chi W_{i\nu}(\chi) \int d\chi' W_{i'\nu'}(\chi') \\ &\qquad \cdot \int \frac{dk}{k} \frac{2}{\pi} k^3 P(k) D(\chi) D(\chi') j_{\ell}(k\chi) j_{\ell}(k\chi'), \end{aligned} \quad (8) $$

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

$$ W_{i\nu}(\chi) = \begin{cases} \frac{\nu}{\Delta\nu} b_i(\chi) M_{0,i}(\chi) A_0(\chi) & \text{if } \chi_{i\nu}^{\text{min}} < \chi < \chi_{i\nu}^{\text{max}} \\ 0 & \text{otherwise,} \end{cases} \quad (9) $$

where $\chi_{i\nu}^{\text{min/max}} = \chi(z_{i\nu}^{\text{min/max}})$ and $z_{i\nu}^{\text{min/max}} = (\nu'_{\text{rf}}/\nu^{\text{max/min}}) - 1$. $z_{i\nu}^{\text{min/max}}$ and $\chi_{i\nu}^{\text{min/max}}$ denote the corresponding redshift and comoving distance that can be probed by the ith line in the frequency channel $\nu$ spanning the frequency range $\nu_{\text{min}} < \nu < \nu_{\text{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

$$ C_{\ell,\ell'\nu\nu'} = \int \frac{d\chi}{\chi^2} W_{i\ell}(\chi) W_{i'\ell'}(\chi) D^2(\chi) P\left(\frac{\ell + \frac{1}{2}}{\chi}\right). \quad (10) $$

We define

$$ M_i(\chi) = b_i(\chi) M_{0,i}(\chi). \quad (11) $$

Hereafter, the term “bias-weighted luminosity” refers to $M_i(\chi)$ and, similarly, “bias-weighted intensity” refers to the quantity $b_i(\chi)\nu I_{\nu}^{i}(\chi)$. We also define

$$ A(\chi) = D(\chi) A_{0,i}(\chi). \quad (12) $$

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

$$ \begin{aligned} C_{\ell,\ell'\nu\nu'} = \frac{\nu}{\Delta\nu} \frac{\nu'}{\Delta\nu'} \frac{\Delta\chi_{\ell\ell'\nu\nu'}}{\chi_{\ell\ell'\nu\nu'}^2} A^2(\chi_{\ell\ell'\nu\nu'}) \\ \cdot M_i(\chi_{\ell\ell'\nu\nu'}) M_{i'}(\chi_{\ell'\ell\nu\nu'}) P\left(\frac{\ell + \frac{1}{2}}{\chi_{\ell\ell'\nu\nu'}}\right), \end{aligned} \quad (13) $$

where $\chi_{\ell\ell'\nu\nu'}$ and $\Delta\chi_{\ell\ell'\nu\nu'}$ denote the center and the width of the overlapping comoving distance between $\chi_{i\nu}^{\text{min}} < \chi < \chi_{i\nu}^{\text{max}}$ and $\chi_{i'\nu'}^{\text{min}} < \chi < \chi_{i'\nu'}^{\text{max}}$, respectively. For the autospectra of a line

<footnote>$^7$Throughout this manuscript, $\nu$ is referred to as the observed frequency, and the rest-frame frequency is denoted by $\nu_{\text{rf}}$.</footnote>

<footnote>$^8$The redshift $z$ and the comoving distance $\chi$ will be used interchangeably to describe the LOS distance.</footnote>
(i = i'), if there is no overlap between the filter profile, which is the case we consider in this work, the $C_{\ell,\nu\nu'\ell'\nu'}^{\text{line}}$s are only nonzero for the same spectral channel (i.e., $\nu = \nu'$). For two different lines (i $\neq$ i'), the $C_{\ell,\nu\nu'\ell'\nu'}^{\text{s}}$s are only nonzero when the two channels probe the two emission lines i and i' from the same redshift. Figure 3 presents an example angular-power-spectrum matrix $C_{\ell}^{\text{line}}$, where the line signal model is detailed in Section 4.

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

$$C_{\ell,\nu\nu'}^{n} = \sigma_{n}^{2}(\nu)\Omega_{\text{pix}}\delta_{\nu,\nu'}^{K}, \quad (14)$$

where $\delta^{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 $\ell$, $C_{\ell}^{d}$, is a random sample from a Wishart distribution with a scale matrix given by $C_{\ell}$ (Equation (7)) and the degree of freedom $n_{\ell}$, where $n_{\ell}$ is the number of $\ell$ 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\nu}$. 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$\beta$ 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}(\chi)$ and $M_{0,i}(\chi)$ 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}(\chi)M_{0,i}(\chi)$ from the data rather than the two terms individually. More detailed investigations considering the RSD effect are left for future work.

<figure><img src="image_3.png" /><figcaption>Figure 1. Top: SPHEREx spectral resolution of each channel. Bottom: SPHEREx surface brightness sensitivity per spectral channel in a 6"$^2$ 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 $\mu$m).</figcaption></figure>

## 3. Survey Setup

We demonstrate our technique with a survey setup similar to the deep-field survey of the SPHEREx mission<fnref n="9" /> (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 $\mu$m, with a pixel size of 6"$^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 $\mu$m (51 spectral channels), $R=35$ for band 4 between 2.42 and 3.82 $\mu$m (17 spectral channels), $R=110$ for band 5 between 3.82 and 4.42 $\mu$m (17 spectral channels), and $R=130$ for band 6 between 4.42 and 5.00 $\mu$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$^2$ each, with the noise rms ~50 times lower than its all-sky survey (Figure 1).

Here, we consider a survey setup similar to SPHEREx deep fields totaling 200 deg$^2$ (100 deg$^2$ in each ecliptic pole) in this study, corresponding to a sky fraction of $f_{\text{sky}} = 0.48\%$.<fnref n="10" /> We use the first four bands of SPHEREx spanning 0.75–3.82 $\mu$m and assume nonoverlapping top-hat filters equally spaced in

<footnote n="9">⁹ http://spherex.caltech.edu</footnote>

<footnote n="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.</footnote>
<table><caption>**Table 1**
Spectral Lines Modeled in This Work</caption><thead><tr><th>Line</th><th>λ<sub>rl</sub></th><th>SPHEREx Coverage</th><th>r<sub>i</sub></th><th>A<sub>i</sub></th></tr></thead><tbody><tr><td>Hα</td><td>0.6563</td><td>0.14 &lt; z &lt; 4.82(6.62)</td><td>1.27</td><td>1.0</td></tr><tr><td>[O III]</td><td>0.5007</td><td>0.50 &lt; z &lt; 6.63(8.99)</td><td>1.32</td><td>1.32</td></tr><tr><td>Hβ</td><td>0.4861</td><td>0.54 &lt; z &lt; 6.68(9.29)</td><td>0.44</td><td>1.38</td></tr><tr><td>[O II]</td><td>0.3727</td><td>1.01 &lt; z &lt; 9.25(12.4)</td><td>0.71</td><td>0.62</td></tr></tbody></table>

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

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''^2$ and the surface brightness sensitivity in each channel given by the public products.<fnref n="11" /> 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.

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 \gtrsim 4$.

We choose the multipole mode range of $50 < \ell < 350$ in our analysis. The minimum $\ell$ mode ($\ell_{min} = 50$) corresponds to SPHEREx's field of view of $3.5^\circ$ (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 $\ell$ mode ($\ell_{max} = 350$) is made to restrict

our analysis to linear clustering scale (see Appendix B for details).

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

$$n_{\ell}^{\alpha} = f_{\text{sky}} [(\ell_{\text{max}}^{\alpha})^2 - (\ell_{\text{min}}^{\alpha})^2], \quad (15)$$

where $f_{sky} = 0.48\%$ is the fraction of sky area in the 200 deg$^2$ 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:

$$\left[ \frac{L_i}{\text{erg s}^{-1}} \right] = r_i \left[ \frac{\text{SFR}}{M_{\odot} \text{yr}^{-1}} \right], \quad (16)$$

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):

$$\frac{\text{SFRD}(z)}{M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}} = 0.015 \frac{(1+z)^{2.7}}{1 + [(1+z)/2.9]^{5.6}}. \quad (17)$$

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

$$M_{0,i}(z) = \frac{dL_i(z)}{dV} = r_i \cdot \text{SFRD}(z). \quad (18)$$

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:

$$b_i(z) = \frac{\int dM \frac{dn}{dM}(M, z) b_h(M, z)}{\int dM \frac{dn}{dM}(M, z)}, \quad (19)$$

where $\frac{dn}{dM}$ is the halo mass function (Sheth & Tormen 1999) and $b_h$ is the halo bias (Sheth et al. 2001). Our bias prescription

<footnote n="11">¹¹ Data downloaded from: https://github.com/SPHEREx/Public-products/blob/master/Surface_Brightness_v28_base_cbe.txt.</footnote>
<figure><img src="image_4.png" /><figcaption>**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 $\sigma_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.</figcaption></figure>

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 \sim 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.

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 $C_{\ell}^{\text{line}}$ from our model at the lowest-$\ell$ bin centered at

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top: the line power spectrum matrix $C_{\ell}^{\text{line}}$ (Equation (23)) in the lowest-$\ell$ mode centered at $\ell = 91$. The signals from different pairs of lines are labeled in the figure. Bottom: the line 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).</figcaption></figure>

$\ell = 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.

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,\ell'}^{\text{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}_{l,\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)$}:

$$M_i(z) = \sum_{j=1}^{N_m} c_{ij} \hat{M}_j(z), \qquad (20)$$

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:

$$C_{l,\nu\nu'ii'jj'} = \sum_{j=1}^{N_m} \sum_{j'=1}^{N_m} c_{ij} c_{i'j'} \hat{C}_{l,\nu\nu'ii'jj'}, \qquad (21)$$

where

$$\hat{C}_{l,\nu\nu'ii'jj'} = \frac{\nu}{\Delta\nu} \frac{\nu'}{\Delta\nu'} \frac{\Delta\chi_{ii'\nu\nu'}}{\chi_{ii'\nu\nu'}^2} A^2(\chi_{ii'\nu\nu'}) \cdot \hat{M}_j(\chi_{ii'\nu\nu'}) \hat{M}_{j'}(\chi_{ii'\nu\nu'}) P\left(\frac{\ell + \frac{1}{2}}{\chi_{ii'\nu\nu'}}\right). \qquad (22)$$

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

$$C_{l}^{\text{line}} = \sum_{i=1}^{N_{\text{line}}} \sum_{i'=1}^{N_{\text{line}}} \sum_{j=1}^{N_{\text{m}}} \sum_{j'=1}^{N_{\text{m}}} c_{ij} c_{i'j'} \hat{C}_{l,ii'jj'}. \qquad (23)$$

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:

$$\frac{\partial C_l}{\partial c_{ij}} = \frac{\partial C_l^{\text{line}}}{\partial c_{ij}} = 2 \sum_{i'=1}^{N_{\text{line}}} \sum_{j=1}^{N_{\text{m}}} c_{i'j'} \hat{C}_{l,ii'jj'}. \qquad (24)$$

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

<figure><img src="image_6.png" /><figcaption>**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, \dots, 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$.</figcaption></figure>

### 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_j$s:

$$\hat{M}_j(z) = \begin{cases} z - z_j & \text{if } z_j < z \\ 0 & \text{otherwise.} \end{cases} \qquad (25)$$

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.

The linear combination of this basis set spans the piecewise linear functions anchored at $z = z_j + 1$ (i.e., $z = 0, 1, \dots, 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

$$m_{ij} = M_i(z_j + 1), \qquad (26)$$

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, \dots, 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 $\Theta$ from the data power spectra {$C_{\ell}^d$} in $N_{\ell}$ multipole bins using a Bayesian framework. The posterior probability distribution $p(\Theta|\{C_{\ell}^d\})$ is given by

$$p(\Theta|\{C_{\ell}^d\}) \propto \mathcal{L}(\{C_{\ell}^d\}|\Theta)\pi(\Theta), \qquad (27)$$

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

Here, our parameter set $\Theta$ 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 {$\theta_{ij}$}, where $\theta_{ij} = \log m_{ij}$ (see Appendix C for details). Our objective is to determine the maximum a posteriori solution for {$\theta_{ij}$}. By doing so, the positivity constraint on $m_{ij}$s will be automatically satisfied. We set flat

priors on {$\theta_{ij}$}, which effectively give the logarithmic priors on {$m_{ij}$}.<fnref n="12" />

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,\ell'}/\text{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 $\ell$ bin:

$$\begin{aligned} \log \mathcal{L}(\{C_{\ell}^d\}|\Theta) &= -\frac{1}{2}\sum_{\ell} n_{\ell} \log \mathcal{N}(C_{\ell}^d, C_{\ell}(\Theta)) \\ &= -\frac{1}{2}\sum_{\ell} n_{\ell}[\text{Tr}(C_{\ell}^d C_{\ell}^{-1}(\Theta)) \\ &\quad + \log |C_{\ell}(\Theta)| + N_{\ell} \log(2\pi)], \end{aligned} \qquad (28)$$

where $n_{\ell}$ is the number of modes in each $\ell$ bin (Equation (15)), and $C_{\ell}(\Theta)$ is the modeled power spectrum given the parameter set $\Theta$.

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,\ell'}/\text{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 {$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 $\Theta_{\max}$ that maximizes the log-likelihood, given the angular power spectra from the data, {$C_{\ell}^d$}:

$$\Theta_{\max} = \max_{\Theta} \log \mathcal{L}(\{C_{\ell}^d\}|\Theta). \qquad (29)$$

The Newton-Raphson algorithm iteratively updates the current parameter set from $\Theta_{\ell}$ to $\Theta_{\ell+1}$ through the equation:

$$\Theta_{\ell+1} = \Theta_{\ell} - \eta H^{-1}g, \qquad (30)$$

where $g = \nabla \log \mathcal{L}$ and $H = \nabla(\log \mathcal{L})\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 $\eta$ serves as the learning rate, determining the step size of the update. In each iteration of the Newton-Raphson update, we initiate with $\eta = 2$ and

<footnote n="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.</footnote>
subsequently verify whether the new proposed parameter set yields a higher log-likelihood value. If not, we reduce the value of $\eta$ by half until the condition is satisfied.

The Newton–Raphson optimization is performed with the logarithmically transformed parameter set $\{\theta_{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_{\ell}$ on parameters, thus expediting the optimization process. This approximation approaches the exact expression when $\Theta_{\ell}$ is in proximity to $\Theta_{\text{max}}$, and we have confirmed the successful convergence of our algorithm using this approximation.

### 5.3.2. Fisher Matrix

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

$$ 
\begin{aligned}
F_{\alpha\beta} &= -\left\langle \frac{\partial^2 \log \mathcal{L}}{\partial \theta_{\alpha} \partial \theta_{\beta}} \right\rangle \\
&= \frac{1}{2} \sum_{\ell} n_{\ell} \operatorname{Tr} \left( \mathbf{C}_{\ell}^{-1} \frac{\partial \mathbf{C}_{\ell}}{\partial \theta_{\alpha}} \mathbf{C}_{\ell}^{-1} \frac{\partial \mathbf{C}_{\ell}}{\partial \theta_{\beta}} \right),
\end{aligned}
\quad (31) $$

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

$$ \operatorname{Cov}(\theta_{\alpha}, \theta_{\beta}) = F_{\alpha\beta}^{-1}. \quad (32) $$

The Fisher matrix calculation requires the derivative of $C_{\ell}$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'(z')$ for the two given lines $i$ and $i'$ at redshift $z$ and $z'$ by the expression

$$ \operatorname{Cov}[M_i(z), M_i'(z')] = \sum_{j=1}^{N_{\text{m}}} \sum_{j'=1}^{N_{\text{n}}} \operatorname{Cov}[c_{ij}, c_{i'j'}] \hat{M}_j(z) \hat{M}_{j'}(z'). \quad (33) $$

## 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 $\{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, \dots, 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 $\{C_{\ell}^d\}$ in eight $\ell$ bins, and we fit for 28 parameters (seven parameters for each of the four lines) in total.

With the input spectra $\{C_{\ell}^d\}$, we first use the Newton–Raphson method to find the parameter set $\Theta_{\text{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 ($\{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\sigma$ 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.

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 $\gtrsim 10\sigma$ level at $z<2$. For the brightest line, $H\alpha$, the $10\sigma$ sensitivity can be achieved out to $z \sim 3$, with the signal-to-noise ratio (S/N) peaking at $\sim 100\sigma$ 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.

We can also estimate the constraining power on the SFRD using our constraints on the bias-weighted intensity. Here, we propagate our $H\alpha$ constraint to the SFRD, assuming that the bias-weighted intensity of $H\alpha$ 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\alpha$ 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 \lesssim 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.

## 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, $(\sigma_n^{\text{SPHEREx-deep}})^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><img src="image_7.png" /><figcaption>**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 {$C_l^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\sigma$ and $2\sigma$ constraints derived from the Fisher matrix.</figcaption></figure>

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$\alpha$, 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_{\text{sky}}^{\text{SPHEREx-all}} \approx 75\%$, we get a factor of ~150 more sky coverage than the deep field ($f_{\text{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
<figure><img src="image_8.png" /><figcaption>**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\sigma$ 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) \cdot \nu_l(z)$ for each line (colored lines). The solid, dashed, and dotted black lines mark the $1\sigma$, $3\sigma$, and $10\sigma$ sensitivity levels, respectively. Bottom four panels: $1\sigma$ constraints for each line relative to the truth input ($H\alpha$, [O III], $H\beta$, and [O II], from top to bottom, respectively).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The $1\sigma$ SFRD constraints from our fiducial case. We assume that the SFRD is proportional to the bias-weighted intensity of the $H\alpha$ line, and thus the sensitivity on the SFRD is propagated from our $H\alpha$ 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).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The S/N on the bias-weighted intensity at $z=1$ (top), $z=2$ (middle), and $z=3$ (bottom) for each line is shown as a function of the noise level relative to the fiducial SPHEREx deep-field noise, $\sigma_n^{\text{SPHEREx-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 ($\sigma_n = 0$). The solid, dashed, and dotted gray horizontal lines mark the $1\sigma$, $3\sigma$, and $10\sigma$ sensitivity levels, respectively. The solid and dotted gray vertical lines mark $\sigma_n^2/(\sigma_n^{\text{SPHEREx-deep}})^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.</figcaption></figure>

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
<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>

example, detector crosstalk. We set the eigenvalues of this correlated noise matrix to be similar to our fiducial case.

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.<fnref n="13" /> 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 ($\ell < 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-$\ell$ 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

<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

shown in Figure 10. While there are ~30 times more modes between $350 < \ell < 2000$ compared to our fiducial setup of $50 < \ell < 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-$\ell$ modes. Therefore, there is only a factor of ~1.6 gain in the line sensitivity by the inclusion of the higher-$\ell$ 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.

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.

### 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 $i$th interloper signal will be projected from the redshift $z_i$ to the target line redshift $z_t$, where $z_i = \lambda_i(1 + z_t)/\lambda_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

<footnote n="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.</footnote>
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

power spectrum anisotropic (Cheng et al. 2016; Lidz & Taylor 2016; Gong et al. 2020), which introduces nonzero $\ell > 0$ modes when expanding the 3D power spectrum with Legendre polynomials (Bernal et al. 2019; 2021; Gong et al. 2020).

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

$$P_{\ell}(k, z_{\ell}) = \sum_{i=1}^{N_{\text{line}}} P_{\ell}^{i}(k, z_{\ell}). \qquad (34)$$

The $P_{\ell}$ from the $\ell$th line is given by

$$P_{\ell}^{i}(k, z_{\ell}) = \frac{1}{(q_{\perp}^{i})^{2} q_{\parallel}^{i}} \frac{2\ell + 1}{2} \cdot \int_{-1}^{1} d\mu P^{i}(k_{i}, \mu_{i}, z_{i}) \mathcal{L}_{\ell}(\mu), \qquad (35)$$

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 $P^{i}(k_{i}, \mu_{i}, z_{i})$ is the intrinsic power spectrum of the $\ell$th line at redshift $z_{i}$ and the corresponding Fourier mode $k_{i}$ and cosine angle $\mu_{i}$. For the intrinsic power spectrum, $P^{i}(k_{i}, \mu_{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_{\ell})$ for the three lowest-$\ell$ modes.<fnref n="14" /> 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$\alpha$), 1.61 ([O III]), 1.68 (H$\beta$).

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 $\{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.

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 $\sim 1\sigma$ range at $z < 4$. Beyond this redshift, however, we obtain biased results for the faintest line, [O II], at $z \sim 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 $\{C_{\ell}^{d}\}$ with doubled redshift resolution, i.e., setting the redshift anchoring point spacing to $\Delta z = 0.5$ instead of the fiducial case of $\Delta 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.

### 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,\ell'}$. 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

<footnote n="14">¹⁴ Since the line power spectrum $P^{i}(k_{i}, \mu_{i}, z_{i})$ is an even function of the LOS cosine angle $\mu$, the odd multipole modes vanish.</footnote>
<figure><img src="image_14.png" /><figcaption>**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\sigma$ constraint of the bias-weighted intensity ($b_ℓ(z)U_{ℓ}(z)$) from our fiducial case (Section 6), which gives an S/N of 3.9 at $z = 2.5$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**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, \dots, 6$ (points on dashed lines).</figcaption></figure>

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_{ℓ,νν'}$ 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\alpha$ and $H\beta$ luminosity is fixed by atomic physics ($L_{H\beta}/L_{H\alpha} = 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).
<figure><img src="image_16.png" /></figure>

**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\sigma$ constraints from our inference. The input mock data contain sample variance fluctuations. Second panel: the S/N on $b(z) \cdot v_I(z)$ for each line (colored lines). The solid, dashed, and dotted black lines mark the $1\sigma$, $3\sigma$, and $10\sigma$ sensitivity levels, respectively. Bottom four panels: $1\sigma$ constraints for each line relative to the truth input (H$\alpha$, [O III], H$\beta$, and [O II], from top to bottom, respectively).

### 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.

<figure><img src="image_17.png" /></figure>

**Figure 15.** Reconstruction results of the same input signal as in Figure 14, with finer redshift anchoring points of $\Delta z = 0.5$ instead of the fiducial case of $\Delta z = 1$ (Figure 14).

### 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, allow-
ing 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 like-
lihood on the $C_{\ell,\ell'}/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,\ell'}$, 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 suscep-
tible 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\beta$ 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,\ell'}$), 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,\ell'}/S$).

We employ Bayesian analysis to infer the bias-weighted
intensity of each line from the data covariance $C_{\ell,\ell'}/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 $\gtrsim 10\sigma$ level at $z < 2$. For the brightest line, Hα, the $10\sigma$ sensitivity can be achieved out to $z \sim 3$, with S/N peaks at $\sim 100\sigma$ 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 $\nu$ is given by

$$ \nu I_{\nu}^{i}(\nu, \hat{n}) = \int d\chi \frac{d[\nu_{rf} L_{\nu}^{i}(\nu_{rf}, \chi, \hat{n})]}{dV} \frac{D_{A}^{2}(\chi)}{4\pi D_{L}^{2}(\chi)}. \quad (A1) $$

Here, $D_L$ is the luminosity distance, $D_A$ is the comoving angular diameter distance (equal to the comoving distance in a flat Universe), $\nu_{rf} = \nu(1+z)$ is the rest-frame frequency corresponding to the observed frequency $\nu$ at redshift $z$, and $L_{\nu}^{i}(\nu_{rf}, \chi, \hat{n}) = dL^{i}(\nu_{rf}, \chi, \hat{n})/d\nu_{rf}$ is the line profile (specific luminosity density) of the ith line at the angular position $\hat{n}$ and rest-frame frequency $\nu_{rf}$. Note that $\nu_{rf} L_{\nu}^{i}(\nu_{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_{rf}, \chi, \hat{n})$ as a Dirac delta function $\delta^D$ at the rest-frame line frequency $\nu_{rf}$. We define the comoving line luminosity density:

$$ M_{0,i}(\chi, \hat{n}) = \frac{dL_i(\chi, \hat{n})}{dV} = \frac{d}{dV} \int d\nu_{rf} L_{\nu}^{i}(\nu_{rf}, \chi, \hat{n}). \quad (A2) $$

With this, we get

$$ \begin{aligned} & \frac{d[\nu_{rf} L_{\nu}^{i}(\nu_{rf}, \chi, \hat{n})]}{dV} \\ &= M_{0,i}(\chi, \hat{n}) \nu_{rf} \delta^{D}(\nu_{rf} - \nu_{rf}^{i}) \\ &= M_{0,i}(\chi, \hat{n}) \nu \delta^{D}\left(\nu - \frac{\nu_{rf}^{i}}{1+z}\right) \\ &= M_{0,i}(\chi, \hat{n}) \nu \frac{d\chi}{d\nu} \Big|_{z_{i\nu}} \delta^{D}(\chi - \chi_{i\nu}) \\ &= M_{0,i}(\chi, \hat{n}) \frac{c(1+z_{i\nu})}{H(z_{i\nu})} \delta^{D}(\chi - \chi_{i\nu}), \end{aligned} \quad (A3) $$

where $z_{i\nu} = \nu_{rf}^i/\nu - 1$ is the redshift of the ith line at the observed frequency $\nu$, and $\chi_{i\nu}$ 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

$$ \nu I_{\nu}^{i}(\nu, \hat{n}) = \frac{c(1 + z_{i\nu})}{H(z_{i\nu})} M_{0,i}(\chi_{i\nu}, \hat{n}) A_{0}(\chi_{i\nu}), \quad (A4) $$

where we define

$$ A_0(\chi) = \frac{D_A^2(\chi)}{4\pi D_L^2(\chi)}. \quad (A5) $$

### A.2. Window Function

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

$$ \langle \nu I_{\nu}^{i}(\nu) \rangle_{R} = \frac{\int d\nu' \nu' I_{\nu}(\nu') R(\nu')}{\int d\nu' R(\nu')}. \quad (A6) $$

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

$$ R(\nu) = \frac{1}{\Delta\nu} \mathcal{I}(\nu; \nu^{\min}, \nu^{\max}), \quad (A7) $$
where $\Delta\nu = \nu^{\text{max}} - \nu^{\text{min}}$ is the filter width, and the indicator function $\mathcal{I}(\nu; \nu^{\text{min}}, \nu^{\text{max}})$ is defined by

$$ \mathcal{I}(\nu; \nu^{\text{min}}, \nu^{\text{max}}) = \begin{cases} 1 & \text{if } \nu^{\text{min}} < \nu < \nu^{\text{max}} \\ 0 & \text{otherwise.} \end{cases} \qquad (A8) $$

The window function of the power spectrum at the channel centered at $\nu$, denoted as $W_{it}(\chi)$ in Equation (10), relates to the line intensity field by the following expression:

$$ \langle \nu I_{\nu}^{i}(\nu) \rangle_{R} = \int d\chi \, W_{0,i\nu}(\chi), \qquad (A9) $$

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

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

$$ \begin{aligned} \langle \nu I_{\nu}^{i}(\nu) \rangle_{R} &= \frac{1}{\Delta\nu} \int d\nu' \nu' I_{\nu}^{i}(\nu') \mathcal{I}(\nu'; \nu^{\text{min}}, \nu^{\text{max}}) \\ &= \frac{1}{\Delta\nu} \int d\chi \, \nu' \frac{c(1+z)}{H(z)} M_{0,i}(\chi) A_{0}(\chi) \\ & \qquad \cdot \mathcal{I}(\chi; \chi_{i\nu}^{\text{min}}, \chi_{i\nu}^{\text{max}}) \frac{d\nu'}{d\chi} \Big|_{i\nu'} \\ &= \int d\chi \, \frac{\nu'}{\Delta\nu} M_{0,i}(\chi) A_{0}(\chi) \mathcal{I}(\chi; \chi_{i\nu}^{\text{min}}, \chi_{i\nu}^{\text{max}}), \quad (A10) \end{aligned} $$

where $\chi_{i\nu}^{\text{min/max}} = \chi(z_{i\nu}^{\text{min/max}})$ and $z_{i\nu}^{\text{min/max}} = (\nu_{\text{rf}}^{i}/\nu^{\text{min/max}}) - 1$. Therefore, the window function of the $i$th line at the channel $\nu$ is

$$ \begin{aligned} W_{i\nu}(\chi) &= \frac{\nu}{\Delta\nu} b_i(\chi) M_{0,i}(\chi) A_0(\chi) \mathcal{I}(\chi; \chi_{i\nu}^{\text{min}}, \chi_{i\nu}^{\text{max}}) \\ &= \begin{cases} \frac{\nu}{\Delta\nu} b_i(\chi) M_{0,i}(\chi) A_0(\chi) & \text{if } \chi_{i\nu}^{\text{min}} < \chi < \chi_{i\nu}^{\text{max}} \\ 0 & \text{otherwise.} \end{cases} \end{aligned} \qquad (A11) $$

### A.3. Poisson Noise

The Poisson noise power from lines $i$ and $i'$ in channels $\nu$ and $\nu'$ is given by

$$ \begin{aligned} C_{\ell, \nu \nu' i i'}^{P} &= \int d\chi \, D_{A}^{2}(\chi) \int dM \frac{dn}{dM} \\ & \quad \cdot \frac{\nu L_{\nu}^{i}(M, \chi) \nu' L_{\nu'}^{i'}(M, \chi)}{4\pi D_{L}^{2}(\chi) 4\pi D_{L}^{2}(\chi)} \\ &= \frac{\nu}{\Delta\nu} \frac{\nu'}{\Delta\nu'} \int d\chi \, D_{A}^{2}(\chi) \left[ \frac{1}{4\pi D_{L}^{2}(\chi)} \right]^{2} \\ & \quad \cdot \int dM \frac{dn}{dM} L^{i}(M, \chi) L^{i'}(M, \chi), \end{aligned} \qquad (A12) $$

where $dn/dM$ is the halo mass function (Sheth & Tormen 1999). Here, we assume the line luminosities $L^i$ are functions of the halo mass $M$ and ignore scatters in this relation. The shot noise power depends on the details of the $L^i-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^i-M$, which allows us to relate the

<figure><img src="image_18.png" /><figcaption>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.</figcaption></figure>

last integral in Equation (A12) to the luminosity density $M_{0,i}$ by

$$ \begin{aligned} & \int dM \frac{dn}{dM} L_{\nu}^{i}(M, \chi) L_{\nu'}^{i'}(M, \chi) \\ &= M_{0,i}(\chi) M_{0,i'}(\chi) \frac{\int dM \frac{dm}{dM} M^2}{\left(\int dM \frac{dm}{dM} M\right)^2}. \end{aligned} \qquad (A13) $$

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.

## Appendix B
Redshift and Multipole Ranges

Our choice of $z_{\text{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\sigma$ point-source sensitivity of $m \sim 21.6$ per channel (at $\lambda < 3.82$ µm) in the deep fields$^{15}$ (the blue solid line in the top panel of Figure 17).

<footnote>$^{15}$ https://github.com/SPHEREx/Public-products/blob/master/Point_Source_Sensitivity_v28_base_cbe.txt</footnote>
<figure><img src="image_19.png" /><figcaption>**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.</figcaption></figure>

Considering that the SPHEREx photometric redshift fitting can achieve high accuracy for galaxies with $\gtrsim 3\sigma$ per channel sensitivity, this corresponds to a masking depth of $\sim 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 \sim 0.7$, there is $\lesssim 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.

The choice of the maximum-$\ell$ 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 \gtrsim 0.2h\,\text{Mpc}^{-1}$ (Figure 18), and thus we set $\ell_{max}$ to correspond to a transverse comoving maximum-k mode of $0.2h\,\text{Mpc}^{-1}$ at $z_{min} = 0.7$ ($k_{max} \sim \ell_{max}/\chi(z_{min}) \sim 0.2$).

<figure><img src="image_20.png" /><figcaption>**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 $\ell$, 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 $\gtrsim 10\%$ at $k \sim 0.2$, corresponding to $\ell \sim 350$. Therefore, we choose $\ell_{max} = 350$ for our analysis.</figcaption></figure>

## Appendix C: Parameter Transformation

In this section, we present the relationship between the three sets of parameters: {$c_{ij}$}, {$m_{ij}$}, and {$\theta_{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

$$m_{ij} = M_i(z_j + 1) = \sum_{j=1}^{N_m} c_{ij} \hat{M}_j(z_j + 1), \quad (C1)$$

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 $J_i^{cm}$:

$$c_i = J_i^{cm} m_i, \qquad (C2)$$

and

$$m_i = (J_i^{cm})^{-1} c_i, \qquad (C3)$$

where

$$J_i^{cm} = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ -2 & 1 & 0 & \dots & 0 \\ 1 & -2 & 1 & \dots & 0 \\ 0 & 1 & -2 & \dots & 1 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & \dots & 0 & 1 & -2 & 1 \end{bmatrix}, \qquad (C4)$$

and

$$(J_i^{cm})^{-1} = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ 2 & 1 & 0 & \dots & 0 \\ 3 & 2 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \dots & 3 & 2 & 1 & 0 \\ \dots & 4 & 3 & 2 & 1 \end{bmatrix}. \qquad (C5)$$

The parameter set $\theta_{ij}$ is defined by the logarithm of $m_{ij}$:

$$\theta_{ij} = \log(m_{ij}), \qquad (C6)$$

and thus the Jacobian $J_i^{\theta m}$ is a diagonal matrix, with the diagonal elements given by $m_i$:

$$J_i^{m\theta} = \mathrm{Diag}(m_i). \qquad (C7)$$

It is important to note that this transformation is nonlinear, and unlike $J_i^{cm}$, the Jacobian $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:

$$\nabla^{\beta} \log \mathcal{L} = J^{\alpha\beta} \nabla^{\alpha} \log \mathcal{L}, \qquad (C8)$$

and

$$F^{\beta} = (J^{\alpha\beta})^{T} F^{\alpha} J^{\alpha\beta}, \qquad (C9)$$

where $\alpha$ and $\beta$ represent any two of the parameter sets ($\{c_{ij}\}$, $\{m_{ij}\}$, and $\{\theta_{ij}\}$). Note that the Jacobian matrices follow the chain rule: $J^{\alpha\beta} = J^{\alpha\gamma}J^{\gamma\beta}$. Therefore, for example, we can obtain $J^{c\theta}$ by $J^{c\theta} = J^{cm}J^{m\theta}$.

In our algorithm, we initially compute the likelihood derivative and Fisher matrix in $\{c_{ij}\}$, then transform them to $\{\theta_{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 $\ell$th multipole of the 3D power spectrum from interloper line $i$ projected to the target line $t$ is given by (Equation (35)):

$$P_{\ell}^{i}(k, z_{\ell}) = \frac{1}{(q_{\perp}^{i})^{2} q_{\parallel}^{i}} \frac{2\ell + 1}{2} \cdot \int_{-1}^{1} d\mu P^{i}(k_{i}, \mu_{i}, z_{i}) \mathcal{L}_{\ell}(\mu), \qquad (D1)$$

where $z_i = \lambda_i(1 + z_t)/\lambda_i - 1$ is the redshift of the $i$th 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

$$\begin{aligned} q_{\perp}^{i} &= \frac{D_{A}(z_{i})}{D_{A}(z_{t})}, \\ q_{\parallel}^{i} &= \frac{(1 + z_{i})/D_{H}(z_{i})}{(1 + z_{t})D_{H}(z_{t})}, \end{aligned} \qquad (D2)$$

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

$$P^i(k_i, \mu_i, z_i) = \left(1 + \frac{f_i}{b_i} \mu_i^2\right)^2 W(k_i, \mu_i, z) \cdot b_i^2(z_i) (\nu_{\nu,i}(z_i))^2 D^2(z_i) P(k_i). \qquad (D3)$$

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)\nu_{\nu,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 \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 $\mu_i$, respectively, are given by

$$\begin{aligned} k_i &= \frac{k}{q_{\perp}^i} [1 + \mu^2 ((F_{\mathrm{proj}}^i)^{-2} - 1)]^{1/2}, \\ \mu_i &= \frac{\mu}{F_{\mathrm{proj}}^i} [1 + \mu^2 ((F_{\mathrm{proj}}^i)^{-2} - 1)]^{-1/2}, \end{aligned} \qquad (D4)$$

where $F_{\mathrm{proj}}^i = q_{\parallel}^i/q_{\perp}^i$.

The window function $W(k, \mu, z)$ is given by

$$W(k, \mu, z) = \exp\{-k^2[\sigma_{\perp}^2(1 - \mu^2) + \sigma_{\parallel}^2\mu^2]\}, \qquad (D5)$$

where

$$\begin{aligned} \sigma_{\parallel} &= \frac{c(1 + z)}{H(z)R}, \\ \sigma_{\perp} &= D_{A}(z)\sigma_{\mathrm{beam}}, \end{aligned} \qquad (D6)$$

and $R$ is the spectral resolution, $D_A$ is the comoving angular diameter distance, and $\sigma_{\mathrm{beam}}$ is the beam size that we take for the SPHEREx pixel size of $6''/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.

## ORCID iDs

Yun-Ting Cheng https://orcid.org/0000-0002-5437-0504

Kailai Wang https://orcid.org/0009-0000-4756-4223

Benjamin D. Wandelt https://orcid.org/0000-0002-5854-8269

Tzu-Ching Chang https://orcid.org/0000-0001-5929-4187

Olivier Doré https://orcid.org/0000-0001-7432-2932

## References

Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, PhRvD, 105, 023520

Ade, P. A. R., Anderson, C. J., Barrentine, E. M., et al. 2020, JLTP, 199, 1027

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bartlett, D. J., Kammerer, L., Kronberger, G., et al. 2024a, A&A, 686, A209

Bartlett, D. J., Wandelt, B. D., Zennaro, M., Ferreira, P. G., & Desmond, H. 2024b, A&A, 686, A150

Beane, A., & Lidz, A. 2018, ApJ, 867, 26

Beane, A., Villaescusa-Navarro, F., & Lidz, A. 2019, ApJ, 874, 133

Bernal, J. L., Breysse, P. C., Gil-Marín, H., & Kovetz, E. D. 2019, PhRvD, 100, 123522

Bernal, J. L., Caputo, A., & Kamionkowski, M. 2021, PhRvD, 103, 063523

Bernal, J. L., & Kovetz, E. D. 2022, A&ARv, 30, 5

Béthermin, M., Gkogkou, A., Van Cuyck, M., et al. 2022, A&A, 667, A156

Breysse, P. C. 2022, arXiv:2209.01223

Breysse, P. C., Chung, D. T., Cleary, K. A., et al. 2022, ApJ, 933, 188

CCAT-Prime Collaboration, Aravena, M., Austermann, J. E., et al. 2023, ApJS, 264, 7

Chang, T.-C., Pen, U.-L., Bandura, K., & Peterson, J. B. 2010, Natur, 466, 463

Cheng, Y.-T., & Chang, T.-C. 2022, ApJ, 925, 136

Cheng, Y.-T., Chang, T.-C., Bock, J., Bradford, C. M., & Cooray, A. 2016, ApJ, 832, 165

Cheng, Y.-T., Chang, T.-C., & Bock, J. J. 2020, ApJ, 901, 142

Cheng, Y.-T., Wandelt, B. D., Chang, T.-C., & Doré, O. 2023, ApJ, 944, 151

Chung, D. T., Bangari, I., Breysse, P. C., et al. 2023, MNRAS, 520, 5305

Chung, D. T., Viero, M. P., Church, S. E., et al. 2019, ApJ, 872, 186

Cleary, K. A., Borowska, J., Breysse, P. C., et al. 2022, ApJ, 933, 182

CONCERTO Collaboration, Ade, P., Aravena, M., et al. 2020, A&A, 642, A60

Creque-Sarbinowski, C., & Kamionkowski, M. 2018, PhRvD, 98, 063524

Crites, A. T., Bock, J. J., Bradford, C. M., et al. 2014, Proc. SPIE, 9153, 91531W

Croft, R. A. C., Miralda-Escudé, J., Zheng, Z., et al. 2016, MNRAS, 457, 3541

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

Cunnington, S., & Wolz, L. 2023, arXiv:2312.0728

de Oliveira-Costa, A., Tegmark, M., Gaensler, B. M., et al. 2008, MNRAS, 388, 247

DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv:1611.00036

Diemer, B. 2018, ApJS, 239, 35

Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv:1412.4872

Doré, O., Werner, M. W., Ashby, M., et al. 2016, arXiv:1606.07039

Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018, arXiv:1805.05489

Elvin-Poole, J., Crocce, M., Ross, A. J., et al. 2018, PhRvD, 98, 042006

Feder, R. M., Masters, D. C., Lee, B., et al. 2023, arXiv:2312.04636

Feng, C., Cooray, A., Bock, J., et al. 2019, ApJ, 875, 86

Fonseca, J., Silva, M. B., Santos, M. G., & Cooray, A. 2017, MNRAS, 464, 1948

Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, PhR, 433, 181

Gebhardt, K., Mentuch Cooper, E., Ciardullo, R., et al. 2021, ApJ, 923, 217

Gong, Y., Chen, X., & Cooray, A. 2020, ApJ, 894, 152

Gong, Y., Cooray, A., Silva, M., et al. 2012, ApJ, 745, 49

Gong, Y., Cooray, A., Silva, M. B., et al. 2017, ApJ, 835, 273

Hamimeche, S., & Lewis, A. 2008, PhRvD, 77, 103013

Helgason, K., Ricotti, M., & Kashlinsky, A. 2012, ApJ, 752, 113

Hill, G. J., Gebhardt, K., Komatsu, E., et al. 2008, in ASP Conf. Ser. 399, Panoramic Views of Galaxy Formation and Evolution, ed. T. Kodama, T. Yamada, & K. Aoki (San Francisco, CA: ASP), 115

Hinton, S. R. 2016, JOSS, 1, 00045

Huterer, D., Cunha, C. E., & Fang, W. 2013, MNRAS, 432, 2945

Ihle, H. T., Chung, D., Stein, G., et al. 2019, ApJ, 871, 75

Kaiser, N. 1987, MNRAS, 227, 1

Karkare, K. S., Anderson, A. J., Barry, P. S., et al. 2022, JLTP, 209, 758

Keating, G. K., Bower, G. C., Marrone, D. P., et al. 2015, ApJ, 814, 140

Keating, G. K., Marrone, D. P., Bower, G. C., & Keenan, R. P. 2020, ApJ, 901, 141

Keating, G. K., Marrone, D. P., Bower, G. C., et al. 2016, ApJ, 830, 34

Kennicutt, R. C. J. 1998, ARA&A, 36, 189

Korngut, P. M., Bock, J. J., Akeson, R., et al. 2018, Proc. SPIE, 10698, 106981U

Kovetz, E. D., Viero, M. P., Lidz, A., et al. 2017, arXiv:1709.09066

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Li, T. Y., Wechsler, R. H., Devaraj, K., & Church, S. E. 2016, ApJ, 817, 169

Lidz, A., & Taylor, J. 2016, ApJ, 825, 143

Limber, D. N. 1953, ApJ, 117, 134

LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv:0912.0201

Ly, C., Malkan, M. A., Kashikawa, N., et al. 2007, ApJ, 657, 738

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJL, 763, L20

McBride, L., & Liu, A. 2023, arXiv:2308.00749

Millea, M., & Bouchet, F. 2018, A&A, 617, A96

Moradinezhad Dizgah, A., Nikakhtar, F., Keating, G. K., & Castorina, E. 2022, JCAP, 2022, 026

Morales, M. F., & Wyithe, J. S. B. 2010, ARA&A, 48, 127

Moriwaki, K., Filippova, N., Shirasaki, M., & Yoshida, N. 2020, MNRAS, 496, L54

Moriwaki, K., & Yoshida, N. 2021, ApJL, 923, L7

Moullet, A., Kataria, T., Lis, D., et al. 2023, arXiv:2310.20572

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Dulles, VA: University Science Books)

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Pritchard, J. R., & Loeb, A. 2012, RPPh, 75, 086901

Pullen, A. R., Chang, T.-C., Doré, O., & Lidz, A. 2013, ApJ, 768, 15

Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778

Roy, A., & Battaglia, N. 2024, ApJ, 969, 2

Sato-Polito, G., Bernal, J. L., Kovetz, E. D., & Kamionkowski, M. 2020, PhRvD, 102, 043519

Schaan, E., & White, M. 2021, JCAP, 2021, 068

Serra, P., Doré, O., & Lagache, G. 2016, ApJ, 833, 153

Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Silva, M., Santos, M. G., Cooray, A., & Gong, Y. 2015, ApJ, 806, 209

Silva, M. B., Santos, M. G., Gong, Y., Cooray, A., & Bock, J. 2013, ApJ, 763, 132

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Sun, G., Chang, T. C., Uzgil, B. D., et al. 2021, ApJ, 915, 33

Sun, G., Hensley, B. S., Chang, T.-C., Doré, O., & Serra, P. 2019, ApJ, 887, 142

Sun, G., Moncelsi, L., Viero, M. P., et al. 2018, ApJ, 856, 107

Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, MNRAS, 434, L46

Tegmark, M., Eisenstein, D. J., Strauss, M. A., et al. 2006, PhRvD, 74, 123507

Van Cuyck, M., Ponthieu, N., Lagache, G., et al. 2023, A&A, 676, A62

Vieira, J., Aguirre, J., Bradford, C. M., et al. 2020, arXiv:2009.14340

Visbal, E., & Loeb, A. 2010, JCAP, 2010, 016

Visbal, E., & McQuinn, M. 2023, ApJ, 956, 84

Visbal, E., Trac, H., & Loeb, A. 2011, JCAP, 2011, 010

Wandelt, B. D. 2013, Astrostatistical Challenges for the New Astronomy (New York: Springer), 1013

Yang, S., Pullen, A. R., & Switzer, E. R. 2019, MNRAS, 489, L53

Yue, B., Ferrara, A., Pallottini, A., Gallerani, S., & Vallini, L. 2015, MNRAS, 450, 3829

Zhang, L., Bunn, E. F., Karakci, A., et al. 2016, ApJS, 222, 3

Zhou, M., Tan, J., & Mao, Y. 2021, ApJ, 909, 51