\section{Introduction}
Line intensity mapping (LIM) has emerged as a promising tool to study the three-dimensional large-scale structures by mapping a particular spectral line emission, and infers the line-of-sight distance of the emission sources from the frequency-redshift relation. LIM measures the aggregate emission of all sources to constrain the bulk properties of the galaxies; whereas in traditional galaxy surveys, only the brighter sources can be individually detected. This relatively low spatial resolution and point-source sensitivity requirement in LIM enables the use of small apertures to efficiently scan a large survey volume out to high redshifts. 

Several spectral lines have been proposed for LIM survey. The 21 cm hyperfine emission from neutral hydrogen  \citep{1990MNRAS.247..510S,1997ApJ...475..429M,2008PhRvL.100i1303C,2008MNRAS.383..606W}, the CO rotational lines \citep{2008A&A...489..489R,2010JCAP...11..016V,2011ApJ...730L..30C,2011ApJ...741...70L, 2011ApJ...728L..46G,2014MNRAS.443.3506B, 2013ApJ...768...15P, 2015JCAP...11..028M,2015ApJ...814..140K,2016MNRAS.457L.127B,2016ApJ...817..169L,2016ApJ...830...34K,2017MNRAS.464.1948F,2017MNRAS.468..741B,2019ApJ...872..186C}, the [\ion{C}{ii}] 157.7 $\mu$m fine-structure line \citep{2012ApJ...745...49G,2014ApJ...793..116U,2015ApJ...806..209S,2015MNRAS.450.3829Y,2017MNRAS.464.1948F}, and the Lyman-$\alpha$ emission line \citep{2013ApJ...763..132S,2014ApJ...785...72G,2014ApJ...786..111P, 2016MNRAS.455..725C,2016MNRAS.457.3541C,2017MNRAS.464.1948F,2018MNRAS.481.1320C} are amongst the most studied lines in the LIM regime.

One of the main challenges in LIM is the astrophysical foreground contaminations, including the continuum emission and line interlopers. Although the continuum foregrounds are usually a few orders of magnitude brighter than the lines (a situation more severe for 21cm than for other lines), their smooth spectral feature can be used to distinguish from the line signals. This has been extensively studied in the context of 21 cm LIM \citep[e.g.,][]{2006PhR...433..181F,2006ApJ...648..767M,2009ApJ...695..183B,2012MNRAS.419.3491L,2012ApJ...756..165P,2012MNRAS.423.2518C,2015ApJ...815...51S}. The line interlopers, which originate from sources residing in different redshifts emitting spectral lines in the same observed frequency channel, is another pressing issue for LIM experiments. The two most studied line de-confusion techniques, source masking and cross-correlation, typically rely on external data sets that trace the same cosmic volume: the masking technique makes use of a galaxy survey catalog to identify and remove bright interloper sources \citep{2015MNRAS.452.3408B,2015MNRAS.450.3829Y,2015ApJ...806..209S,2018ApJ...856..107S}, whereas cross-correlation of an LIM survey with an external (or internal)  data set can help extract signals of interest \citep{2009ApJ...690..252L,2010JCAP...11..016V,2012ApJ...745...49G, 2014ApJ...785...72G,2015aska.confE...4C,2015ApJ...806..209S,2010Natur.466..463C,2013ApJ...763L..20M, 2016MNRAS.457.3541C,2019ApJ...872..186C}.  In addition to these two methods, \citet{2016ApJ...825..143L} and \citet{2016ApJ...832..165C} use the anisotropy of the interloper power spectrum arising from projection to the target line redshift to separate the lines. \citet{2020ApJ...894..152G} distinguish the lines from the same projection effect but using the multipole power spectrum. \citet{2014arXiv1403.3727D} propose to use angular fluctuations of the light to reconstruct the 3D source luminosity density.

Most of the existing line de-confusion methods (e.g. cross-correlation and power spectrum anisotropy) only extract the two-point statistics  (power spectrum or correlation function) but lose the phase information of individual line maps, which are valuable for cosmological parameter constraints and systematics control in the data. With individual line maps, one can extract information beyond two-point statistics in the non-Gaussian intensity maps, especially ones from the epoch of reionization (EoR). For example, \citet{2017MNRAS.467.2996B} and \citet{2019ApJ...871...75I} show that the one-point statistics of the intensity field can help constrain the luminosity function model. In addition, individual line maps can be used directly as density tracers for various cross-correlation, multi-tracer analysis, de-lensing of the cosmic microwave background (CMB), and perform consistency tests on different spatial regions with different foreground properties.

In this work, we develop a technique to extract individual line intensity maps from an LIM data set with blended interlopers. Using the fact that when multiple spectral lines emitted by a source are observable in an LIM survey, the redshift of the source can be pinned down by fitting to a set of spectral templates that are unique at each redshift. Without any external tracers or spectroscopic follow-up observations, individual line maps can be directly derived. For demonstration, we apply our technique to simulated data of an LIM survey targeting the EoR [\ion{C}{ii}] line with multiple low-redshift CO interlopers. In this case, the intensity field of the low-$z$ CO lines ($0.5\lesssim z\lesssim 1.5$) can be reconstructed since they can be detected in multiple spectral channels.

\citet{2015ApJ...806..234K} first explored the map-space line de-confusion using the multi-line wavelength information in the context of a pencil-beam spectroscopic survey. In this work, we explore the technique in the LIM regime that has a much lower sensitivity and the spectral resolution. In this regime, our template-fitting-based technique can obtain the signal-to-noise ratio (S/N) of the desired signals by 
using the data from multiple frequency channels.

The recent work by \citet{2020MNRAS.496L..54M} demonstrated the feasibility of LIM phase-space de-confusion with deep learning. They show that, in the absence of noise or foreground components, their algorithm can reconstruct the individual line maps that are mixed in the LIM data set. Their training data generation relies on the assumption of the signal clustering and the line luminosity model, whereas in this work, we develop a line de-confusion technique that only makes use of the spectral feature of the lines, which is more robust against the model uncertainty and the noise.

This paper is organized as follows. First we introduce the model and the survey parameters we used to generate the mock LIM survey data in Sec.~\ref{S: model}. Then, Sec.~\ref{S:methods} describes our line de-blending technique. Sec.~\ref{S:results} presents the results on the fiducial setup. In Sec.~\ref{S:discussion}, we present the performance of the technique with more practical considerations, and discuss its applications and extensions. The conclusion are given in Sec.~\ref{S:conclusion}. Throughout this paper, we consider a flat $\Lambda$CDM cosmology with $n_s=0.97$, $\sigma_8=0.82$, $\Omega_m=0.26$, $\Omega_b=0.049$, $\Omega_\Lambda=0.69$, and $h=0.68$, consistent with the measurement from Planck \citep{2016A&A...594A..13P}.

\section{Mock Light Cone Construction}\label{S: model}
For each spatial pixel in an LIM survey, there are a set of spectral channel measurements. Hereafter, the term ``light cone''  refers to the collection of the spectral measurements in a single pixel. The line de-confusion method introduced in this work is performed on a pixel-by-pixel basis, which only utilizes the spectral information in an individual line of sight (light cone), without taking into account the spatial clustering information, which we leave for future work.

We test our de-confusion technique on simulated light cones. Since the clustering information is not relevant to our technique, we generate light cones that are based solely on the spectral line luminosity function models and not on the clustering properties. Thus, all light cones are independent from one another, and we also ignore the line-of-sight clustering in this work. This allows for both the light cone construction and de-confusion by parallelization to speed up without affecting the quantification of performance.

As a demonstration of the technique, we assume an LIM experiment targeting the redshifted [\ion{C}{ii}] fine-structure emission from the EoR: the LIM data set contains multiple low-$z$ CO rotational transitions from low redshifts as interlopers. We note that the technique can be readily applied to any line-confusion problem at other wavelengths.

\subsection{Line Signal Models}\label{S:model_sig}
We model the line emissions of the redshifted [\ion{C}{ii}] emission and five low-redshift CO $J$-transitions: \{CO(2--1), CO(3--2), CO(4--3), CO(5--4), CO(6-5)\}. In reality, in the sub-mm spectral range of interest (generally in the $\sim$ 200-300 GHz range), there are higher CO $J$ lines (that are fainter), Galactic and extragalactic dust continuum emissions (the cosmic infrared background), CMB radiation, and atmospheric emissions that can all contribute to the measurements. Since none of them will produce strong spectral features that impact the performance of our technique, we will not include them in the light cones. Instead, in Sec.~\ref{S:bg}, we will demonstrate that the continuum foreground mitigation in the data analysis process has negligible a impact on our technique performance.

We use the [\ion{C}{ii}] and CO luminosity function models provided by \citet{2016MNRAS.461...93P}, in which a semi-analytic model including the effect of radiative transfer was used to estimate the CO and [\ion{C}{ii}] luminosity functions as constrained by current observations. Here we adopt their fitted Schechter luminosity functions to construct our light cones.

\subsection{Survey Parameters}\label{S:survey_params}
We consider a mock experiment that has similar survey parameters as the two ongoing EoR [\ion{C}{ii}] LIM experiments, TIME \citep{2014SPIE.9153E..1WC} and CONCERTO \citep{2018A&A...609A.130L}. 

The mock survey covers 200--305 GHz with 70 evenly spaced spectral channels ($\delta\nu = 1.5$ GHz), and the $\Omega_{\rm pix} = 0.43^2$ arcmin$^2$ pixel size. We assume the instrument noise is white and has a Gaussian distribution, with four different per-pixel noise levels of standard deviation $\sigma_n=$ \{$10^3$, $5\times10^3$, $10^4$, $5\times10^4$\} Jy sr$^{-1}$. These values are comparable to the range of expected instrument noise in TIME and CONCERTO. 

\subsection{Light Cone Generation}\label{S:model_lightcone}
\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_lc.pdf}
\caption{\label{F:gen_lc}The steps of constructing a light cone. Left: $N(z)$, the effective number of $\ell_*$ sources in each redshift bin in this light cone. Middle: the signals from all six spectral lines and the sources in the 70 spectral channels. Right: the mock observed light cone (blue) consists of the signals (middle panel) and the $\sigma_n = 10$ kJy Gaussian noise. The black line is the signal component (same as the middle panel) for reference.}
\end{center}
\end{figure*}

Based on the assumed luminosity function models and survey parameters, we populate the light cones with sources drawn from random realizations of the Schechter function model.  We first define a fiducial Schechter function \footnote{For the fiducial Schechter function, we choose the $\phi_*$ and $\alpha$ values of CO(1--0) for $z\leq5$ and [\ion{C}{ii}] for $z>5$
in \cite{2016MNRAS.461...93P}. We interpolate or extrapolate the Schechter function parameters to the desired redshift from the values of $z=[0,1,2,3,4,6]$ given by \cite{2016MNRAS.461...93P}.
} for all of the lines as a function of redshift, 
\begin{equation}\label{E:lfunc}
\Phi(\ell)=\phi_\ast\, (\mathcal{L}/\mathcal{L}_*)^\alpha\, e^{-(\mathcal{L}/\mathcal{L}_*)}.
\end{equation}
We discretize the luminosity and redshift into bins of $\Delta (\mathcal{L}/\mathcal{L}_*)$ in 100 luminosity bins in log-space from $\mathcal{L}/\mathcal{L}_*=10^{-3}$ to $10$, and 2000 redshift bins of $\Delta z=5\times10^{-4}$ in linear space from $z=0$ to $z=10$.  
The expectation value of the source counts $N$ within each $\Delta \mathcal{L}$ and $\Delta z$ bin is given by
\begin{equation}
\langle N\rangle = \Phi(\ell) [\Delta (\mathcal{L}/\mathcal{L}_*)] [\Omega_{\rm pix} D_A^2(z) \frac{d\chi}{dz}\Delta z],
\end{equation}
where the last square bracket is the co-moving volume of the voxel defined by spatial pixel $\Omega_{\rm pix}$ and redshift bin $\Delta z$, $D_A$ is the co-moving angular diameter distance, and $\chi$ is the co-moving distance. 

For each $(\mathcal{L}, z)$ bin, we assign its source counts to a Poisson random number with expectation value $\langle N\rangle$, and then integrate along $\mathcal{L}$ to get the total luminosity in each redshift bin $\mathcal{L}_{\rm tot}(z)$. We define $N(z) \equiv \mathcal{L}_{\rm tot} / \mathcal{L}_*$ as the ``effective number counts'' per redshift bin. The left panel of Fig.~\ref{F:gen_lc}  shows the $N(z)$ of one example light cone.

Next, we assign the line luminosity signals in each redshift bin to $N(z)\,\ell_*^{\rm line}$, where $\ell_*^{\rm line}$ is the $\ell_*$ value of the line in the \citet{2016MNRAS.461...93P} model. We  then project the line signals to their corresponding spectral channels to make the light cones. The middle panel of Fig.~\ref{F:gen_lc} shows an example light cone spectrum from the $N(z)$ shown in the left panel.

Finally, we add a Gaussian random fluctuation with a an rms $\sigma_n$ value to each channel to account for the instrumental noise. The right panel of Fig.~\ref{F:gen_lc} shows the same light cone with a $\sigma_n = 10$ kJy Gaussian noise.

This light cone construction procedure assumes that all of the spectral lines have the same Schechter function parameters ($\phi_*$ and $\alpha$) and thus the same luminosity function shape, and the CO spectral line energy distribution (SLED) is also fixed. That is, by construction, all of the sources have the same line luminosity ratio sets by the relative value of $\ell_*^{\rm line}$. In the main parts of this work, we use this fixed SLED model, and we test the impact of adding the SLED variation on the performance of our technique in Sec.~\ref{S:model_uncertainty}.

We further point out that, even though the luminosity function has been sampled to the faint end ($10^{-3}\ell_*$), the light cone signals are still dominated by a few bright peaks. This indicates that the emission field can be well described by the few bright sources in the data. In other words, to extract the line emission field from a single line, one only needs to determine the redshift and luminosity of those bright sources. This is the main concept of our de-confusion technique, detailed in Sec.~\ref{S:methods}.

\section{Methods}\label{S:methods}
\subsection{Formalism}\label{S:methods_matrix}

The intensity of a light cone in frequency channel $\nu_i$ can be expressed as the linear combination of signals from all $N_z$ redshift bins and the noise $n_i$,
\begin{equation}\label{E:elements}
I(\nu_i) = \sum_{j = 1}^{N_z} \widetilde{A}_{ij}\,N(z_j) + n_i,
\end{equation}
where $N(z_j)$ is the effective number of $\ell_*$ sources in redshift $z_j$ (Sec.~\ref{S:model_lightcone}), and $\widetilde{A}_{ij}$ converts $N(z_j)$ to the observed intensity in channel $\nu_i$. 

$\widetilde{A}$ is an $N_{\rm ch}\times N_{\rm z}$ matrix, where we have $N_{\rm ch}=70$ spectral channels, and we use $N_{\rm z}=2000$ ($\Delta z=5\times 10^{-4}$) redshift bins. In principle, we can set $\Delta z$ to infinitesimally small values;  however, in practice the redshift resolution is limited by the instrument spectral resolution, as the SLEDs of nearby galaxies can be highly degenerate when their respective spectral lines fall in the same set of observed channels. Setting $N_{\rm z}=2000$ in fact gives a much finer redshift resolution compared to the spectral channel width. We will use this fine redshift resolution as a starting point, and discuss the strategy to reduce redshift bins and therefore remove redundant information at the end of Sec.~\ref{S:methods_matrix}.

Most of the elements in $\widetilde{A}_{ij}$ are zeros except for which the sources at $z_j$ emit a spectral line at the observed frequency $\nu_i$. In this case,
\begin{equation}
\widetilde{A}_{ij} \equiv I_*^{\rm line}(z_j) = \ell_*^{\rm line}(z_j)\frac{1}{4\pi D_L^2(z_j)\delta\nu_i\Omega_{\rm pix}},
\end{equation}
where $\ell_*^{\rm line}$ is the line luminosity of an $\ell_*$ source in the model, $D_L$ is the luminosity distance, and $I_*^{\rm line}$ is defined as the observed line intensity of an $\ell_*^{\rm line}(z_j)$ source at $z_j$. Fig.~\ref{F:Istr} shows the $I_*^{\rm line}$ in our assumed model\footnote{The glitch in CO(5--4) at $z\sim 1$ is due to a feature in the Schechter luminosity function model in \cite{2016MNRAS.461...93P}. In their CO(5--4) luminosity function fit at $z=1$, a slightly higher $L_*$ and and a lower $\alpha$ value are derived that cause the slight apparent  discontinuity in redshift, even though the luminosity function does not show an abrupt change at $z=1$.}. Note that Eq.~\ref{E:elements} assumes that all of the sources at the same redshift have the same SLED, and we will first build our technique based on this assumption. In reality, the SLED varies across galaxy type. In Sec.~\ref{S:model_uncertainty}, we will show that our method also works in the realistic level of SLED variation.

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{Istr.pdf}
\caption{\label{F:Istr}Intensity of the CO lines $I_*$ from the sources of characteristic luminosity $\ell_*$.}
\end{center}
\end{figure}

Eq.~\ref{E:elements} can be written in the matrix form,
\begin{equation}\label{E:matrixAraw}
\mathbf{I} = \widetilde{\mathbf{A}}\mathbf{N} + \mathbf{n},
\end{equation}
where $\mathbf{I}$ and $\mathbf{n}$ are $N_{\rm ch}$-element column vectors, $\mathbf{N}$ is an $N_z$-element column vector, and $\widetilde{\mathbf{A}}$ is an $N_{\rm ch}\times N_{\rm z}$ matrix. The top panel of Fig.~\ref{F:A} shows the $\widetilde{\mathbf{A}}$ matrix in our model. The $\widetilde{\mathbf{A}}$ matrix is mostly zeros, and the six curves from left to right  are the six spectral lines, CO(2--1), CO(3--2), CO(4--3), CO(5--4), CO(6--5), and [\ion{C}{ii}].

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{A_raw.pdf}
\includegraphics[width=\linewidth]{A_zbins.pdf}
\caption{\label{F:A} Top: $\widetilde{\mathbf{A}}$ matrix with our model on 70 frequency channels and 2000 redshift bins. $\widetilde{\mathbf{A}}$ are zeros (white) for the majority of the elements, and the six curves from left to right correspond to the six spectral lines, CO(2--1), CO(3--2), CO(4--3), CO(5--4), CO(6--5), and [\ion{C}{ii}]. The color scale indicates line intensities in the fiducial model. The gray shaded regions are the redshifts at which sources can be observed in multiple lines, and thus can be reconstructed in our technique. Bottom Left: $\mathbf{A}$ matrix with $N_{\rm ch}\times N_{\rm z} = 70 \times 265$ size, which is the reduced and normalized $\widetilde{\mathbf{A}}$. The color scale represents the intensities with the fiducial model SLED normalized within each column, i.e. for all column $j$, $\sum_{i=1}^{N_{\rm ch}}\mathbf{A}^2_{ij}=1$.,  Bottom Right: the redshift of the 195 multi-line redshift bins in $\mathbf{A}$, which are the redshift bins that have multiple CO lines observable in our mock survey. The colors label the pairs of detectable CO lines. The redshifts not covered by multiple CO lines cannot be reconstructed with our technique ($0.15\leq z\leq0.51$; $0.72\leq z\leq0.89$; $ z\geq1.87$).}
\end{center}
\end{figure*}

\begin{deluxetable*}{ccccccc}
\tablenum{1}
\tablecaption{\label{T:bands}Frequencies and Redshifts of the Six Defined Broad Bands}
\tablewidth{0pt}
\tablehead{
\colhead{Name} & \colhead{Line} & \colhead{$\nu$ bin index} & \colhead{$\nu$ (GHz)} &  \colhead{$\langle z \rangle (z_{\rm min} - z_{\rm max})$}
}
\startdata
J3 high    &CO(3--2)   & 51--69 (19 bins)   & 200.75--227.75   & 0.61 (0.51--0.72)\\ 
J4 low     &CO(4--3)   & 0--24 (25 bins)  & 268.25--304.25     & 0.61 (0.51--0.72)\\ 
J4 high    &CO(4--3)   & 41--69 (29 bins)   & 200.75--242.75   & 1.09 (0.89--1.30)\\ 
J5 low     &CO(5--4)   & 0--36 (37 bins)  & 250.25--304.25    & 1.09 (0.89--1.30)\\ 
J5 high    &CO(5--4)   & 37--69 (33 bins)   & 200.75--248.75   & 1.59 (1.30--1.87)\\ 
J6 low     &CO(6--5)   & 0--42 (43 bins)  & 241.25--304.25     & 1.56 (1.26--1.87)\\ 
\enddata
\tablecomments{The CO(2--1) and CO(3--2) overlapping redshifts of $0.13<z<0.5$ can also be reconstructed, but they are only covered by four frequency channels, which makes it difficult to quantify the reconstruction performance with sufficient statistical power. Therefore, we ignore these redshifts in our analysis.}
\end{deluxetable*}

The goal of our line de-confusion technique is to solve for the source count vector $\mathbf{N}$ in a given observed light cone data $\mathbf{I}$ and a model $\widetilde{\mathbf{A}}$ matrix. With the $\mathbf{N}$ solution, the intensity map of individual spectral lines $\mathbf{I}_{\rm line}$ can be reconstructed by $\mathbf{I}_{\rm line}= \widetilde{\mathbf{A}}_{\rm line}\mathbf{N}$, where $ \widetilde{\mathbf{A}}_{\rm line}$ is $\widetilde{\mathbf{A}}$   when only the target spectral line signals are turned on.

We further rewrite Eq.~\ref{E:elements} by normalization and reducing the nuisance information. First, we normalize the columns in $\widetilde{\mathbf{A}}$ and move the normalization factor $\mathbf{I}^{\rm norm}_j$ to their $\mathbf{N}$ element:
\begin{equation}\label{E:matrixA}
\mathbf{I} = \mathbf{A}\widetilde{\mathbf{N}} + \mathbf{n},
\end{equation}
where 
\begin{equation}
\begin{split}
\mathbf{A}_{ij} &=  \widetilde{\mathbf{A}}_{ij}\,/\, \mathbf{I}^{\rm norm}_j,\\
\widetilde{\mathbf{N}}_j &= \mathbf{N}_j\, \mathbf{I}^{\rm norm}_j,\\
\sum_{i=1}^{N_{\rm ch}}&\mathbf{A}^2_{ij}=1.
\end{split}
\end{equation}

Next, we will reduce the nuisance or redundant elements in Eq.~\ref{E:matrixA}. The columns of $\mathbf{A}$ are the basis spanning the observed data space. In constructing $\widetilde{\mathbf{A}}$, we simply design the columns to be equally spaced redshift bins as shown in the top panel of Fig.~\ref{F:A}. However, this natural basis is highly degenerate. To remove the nuisance information, we first discard the redshift bins that are zero vectors in $\widetilde{\mathbf{A}}$. The sources in these redshifts do not emit lines in observable frequencies, and thus no information can be used to constrain their $N(z)$.  Second, for the redshift bins containing only one CO line, their normalized columns in $\mathbf{A}$ are identical to other columns having a [\ion{C}{ii}] signal in the same frequency channel. In other words, given an observed data $\mathbf{I}$, we cannot distinguish the origin of the source with a single line emission. Therefore, we combine these identical columns into a single column.  In conclusion, we keep the columns in $\mathbf{A}$ with redshift bins that can be observed with multiple spectral lines, plus an identity matrix for the redshift bins that only have a single detectable line.

For all of the redshifts that can be observed in multiple lines, we design the size of the redshift bins based on the following two competing considerations. On the one hand, the redshift bins have to be small enough to faithfully represent the emitting source distribution. 
On the other hand, finer redshift bins give larger $\widetilde{\mathbf{N}}$ size,  and therefore more unknown parameters to be solved. The information can be compressed by combining some neighboring redshift bins, which are highly degenerate, since they have signals in the same channels with similar amplitudes. Therefore, we design the redshift bins using the following procedures: we (1) generate $\mathbf{A}$ with fine redshift bins ($\Delta z = 5\times 10 ^{-4}$ from $z = 0$ to 10)\footnote{Note that with these fine bins, the same frequency bin can map to multiple redshift bins instead of an one-to-one mapping.}, (2) keep the columns with multiple lines, (3) then identify the group of neighboring columns that have signals in the same sets of channels, and (4) keep the medium bin and discard the others. With this process, we get 195 non-degenerate columns. Hereafter, ``multi-line redshift bins'' refers to these 195 redshift bins that can be detected in multiple channels. Finally, we append the $N_{\rm ch}$-sized identity matrix (that account for  the single-line redshift bins) to these 195 columns to generate $\mathbf{A}$. The bottom left panel in Fig.~\ref{F:A} shows matrix $\mathbf{A}$ which has size $N_{\rm ch} \times N_z = 70\times 265$ (195 multi-line redshifts plus 70 columns of identity matrix). The bottom right panel shows the redshift of the 195 multi-line bins, and their color labels are the pairs of detectable CO lines. We also define six broad bands from these pairs of lines by binning groups of channels. Table~\ref{T:bands} lists the definition of the broad bands. Note that the CO(2--1) and CO(3--2) overlapping redshifts of $0.13<z<0.5$ can also be reconstructed, but they are only covered by four frequency channels, which makes it difficult to quantify the reconstruction performance with sufficient statistical power. Therefore, we ignore these redshifts in our analysis.

Our line de-confusion technique can only solve $\widetilde{\mathbf{N}}$ in the 195 multi-line redshift bins. The last 70 elements in $\widetilde{\mathbf{N}}$ that correspond to the identity matrix in $\mathbf{A}$ are nuisance parameters, since they represent degenerate ``single-line'' signals from different redshifts -- that is to say, we cannot reconstruct the signals in the single-line redshifts, which are the regions not covered by the lines in The bottom left panel of Fig.~\ref{F:A} ($0.15\leq z\leq0.51$; $0.72\leq z\leq0.89$; $ z\geq1.87$). In the following analysis, we will only focus on the reconstruction of $\widetilde{\mathbf{N}}$ in the 195 multi-line redshift bins.

\subsection{Sparse Approximation}\label{S:methods_sparse}
The key step in our de-confusion technique is to solve for $\widetilde{\mathbf{N}}$ in Eq.~\ref{E:matrixA}, given the observed spectrum $\mathbf{I}$ and model $\mathbf{A}$. This type of linear system has been extensively studied in the context of CMB map making, in which $\mathbf{I}$, $\mathbf{n}$, and $\widetilde{\mathbf{N}}$ in Eq.~\ref{E:matrixA} can be analogized to the time-ordered data, time-stream noise, and the pointing matrix, respectively. However, contrary to the map-making problem, our system is an ill-posed problem as there are more unknown variables ($N_z = 265$) than the input data points ($N_{\rm ch} = 70$). Thus the standard map-making algorithm \citep[e.g.,][]{2002PhRvD..65b2003S} cannot be applied.

In Eq.~\ref{E:matrixA}, the columns of $\mathbf{A}$ form a basis for $\mathbf{I}$, and the solution $\widetilde{\mathbf{N}}$ is the  linear combination coefficient. The columns of $\mathbf{A}$ form an over-complete basis, since the $\mathbf{A}$ matrix is only of rank $N_{\rm ch}$, and thus the solution $\widetilde{\mathbf{N}}$ is not unique. Indeed, for any given observed data $\mathbf{I}$, there are infinite $\widetilde{\mathbf{N}}$ that can perfectly fit the input.

Nevertheless, Eq.~\ref{E:matrixA} can be solved with the ``sparse'' condition, which means the preferred solution of $\widetilde{\mathbf{N}}$  is the one with a small number of nonzero elements. With this constraint, we can solve Eq.~\ref{E:matrixA} with the following well-defined optimization problem:
\begin{equation}\label{E:sparse_eqn}
\underset{\left \| \widetilde{\mathbf{N}} \right \|_0}{\mathrm{argmin}}\frac{1}{N_{\rm ch}}\left \| \mathbf{I} - \mathbf{A}\widetilde{\mathbf{N}} \right \|^2_2 < \epsilon^2,
\end{equation}
where the $\ell_0$-norm $\left \| \cdot  \right \|_0$ is the number of nonzero elements, and $\epsilon$ sets the threshold of error tolerance of the fit.

This type of problem is known as ``sparse approximation,'' which has been extensively studied in the context of signal processing and compressive sensing \citep{Candes:2006:RUP:2263435.2272020,Donoho:2006:CS:2263438.2272089}. The sparse approximation algorithms solve the sparse representation of the signal in a ``dictionary'' that is composed of a set of ``atoms,'' and represent the signal in the data in terms of the linear combination of a few atoms in the dictionary. In Eq.~\ref{E:sparse_eqn}, the dictionary is the matrix $\mathbf{A}$, and the atoms are the column vectors of $\mathbf{A}$.

The sparse approximation can only be applied if $k \equiv\left \| \widetilde{\mathbf{N}} \right \|_0 \ll N_{\rm ch}$. Note that the sparsity of the problem is quantify by $k/N_{\rm ch}$ but not $k/N_z$, since $k/N_z$ can always be designed to be arbitrarily small by choosing a large basis (fine redshift bins). However, the degree of freedom in the solution is restricted by the input size $N_{\rm ch}$, and thus the $k\ll N_{\rm ch}$ condition prohibits the algorithm from using more parameters than the input degree of freedom to over-fit the data.

In general, LIM light cones are not sparse, since there is always a large number of faint sources in the typical luminosity function (e.g. Schechter function), so all of the elements in $\widetilde{\mathbf{N}}$ are nonzero. However, as mentioned in Sec.~\ref{S:model_lightcone}, the light cone signals are dominated by only a few bright sources, and the intensity field can be well described by them. Consequently, the parameter $k$ in our problem can be quantified by the ``effective'' number of these bright sources per voxel that contribute most of the emission. Following  \cite{2019ApJ...877...86C}, we define the effective number $N_{\rm eff}$ as:
\begin{equation}
N_{\rm eff}(z) \equiv \frac{\left( V_{\rm vox} \, \int d\ell \, \Phi(\ell,z) \, \ell \right)^2}{V_{\rm vox} \, \int d\ell \, \Phi(\ell,z) \, \ell^2},
\end{equation}
where $V_{\rm vox}$ is the voxel size of the redshift bin. Note that $\Phi$ is the number of sources per luminosity per voxel, and therefore $N_{\rm eff}$ is dimensionless and is  
proportional to the voxel size. $N_{\rm eff}$ can be interpreted as the reciprocal of the effective shot noise in LIM, which is analogous to the $1/N$ shot noise in a galaxy power spectrum. If the luminosity function $\Phi$ follows the Schechter function form, then $N_{\rm eff}$ is (approximately) the number of sources brighter than $\ell_*$, which contributes the majority of the emission. 

We can estimate $k$ by the cumulative $N_{\rm eff}(z)$ along the line of sight per light cone. Fig.~\ref{F:Neff} shows the cumulative $N_{\rm eff}$ in our model. While $N_{\rm eff}\sim 100$ from $z=0$ to 10, the only relevant range is $z\lesssim 2.5$, where CO lines fall in the observed frequency range (see Fig.~\ref{F:A}, top panel). Above $z\sim2.5$, only the high-$z$ [\ion{C}{ii}] lines can be observed, but they are much fainter then the CO signals and the assumed noise level, so they can be treated as background fluctuations. Therefore, for $z\lesssim2.5$, we find $N_{\rm eff}\sim 10\ll N_{\rm ch}=70$, so the sparse condition is qualified in our problem.

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{Neff.pdf}
\caption{\label{F:Neff} The redshift cumulated $N_{\rm eff}$ in our model (black dashed line). The red segments mark the multi-line redshift ranges that have multiple observable CO lines. Thus the de-confusion technique can be applied.}
\end{center}
\end{figure}

\subsection{The Matching Pursuit Algorithm}\label{S:methods_MP}
We use the matching pursuit (MP) algorithm first introduced by \citet{1993ITSP...41.3397M} to solve for Eq.~\ref{E:sparse_eqn}. The MP algorithm iteratively selects an ``atom'' in the ``dictionary'' to project out part of the signals in the data, and keep track of the current solution of the signal and residual for the next step until the stopping criteria is met. In our case, the column vectors in $\mathbf{A}$ are the atoms that form the dictionary space. A detailed description of the MP algorithm is in Appendix~\ref{A:MPA}.

In each step of the MP algorithm, the selected atom is the one that has the maximum inner product with the residual. The S/N of the signals in each step is the ratio of that maximum inner product to the instrument noise level $\sigma_n$ (see Appendix~\ref{A:proof_u_var} for the proof). Therefore, if we set the stopping criteria to be the maximum inner product smaller than $m$ times of $\sigma_n$, then this is an $m$-$\sigma$ detection threshold on the signals (e.g., $m=5$ for a 5$\sigma$ detection).  Note that the detection threshold here is based on the combined information in multiple spectral channels projected onto the dictionary space.

The choice of detection threshold ``$m$'' is a trade-off between the purity and completeness of the source extraction. Higher ``$m$'' values give a higher purity map, whereas lower ``$m$'' values pick out fainter sources at the cost of increased false detections from noise. The optimal value of ``$m$'' depends on the instrument sensitivity and the purpose of the reconstruction map. For example, to reconstruct the line luminosity function, one might use a higher threshold to reduce the false detections at the faint end; whereas to constrain the large-scale structure, a lower threshold is preferred to reduce the shot noise in the power spectrum. An analytical  formalism to determine the optimal threshold is to make use of the Fisher information framework \citep{2019ApJ...877...86C}, where one calculates the expectation value of the desired observable (e.g. power spectrum) as a function of threshold for a given signal model and noise level, and estimates the threshold that optimizes the Fisher information. Alternatively, one can simply perform test simulations with different thresholds to determine the optimal value. Both approaches can provide guidance on choosing the optimal threshold for the problem at hand, and a detailed investigation is beyond the scope of this paper.

\section{Results}\label{S:results}
We present line de-confusion results in the simple case where  mock light cones and the template ($\mathbf{A}$) are both generated from the same signal model (Sec.~\ref{S:model_lightcone}). We demonstrate that in this scenario, our technique is capable of extracting low-$z$ CO signals in the presence of realistic instrumental noise. We discuss the robustness of the performance against uncertainties in the signal model and contamination from astrophysical foregrounds, and extend the application to spectral lines in different wavelengths in Sec.~\ref{S:discussion}.

We quantify the reconstruction performance by computing two statistics on the true and reconstructed data: (1) the Pearson correlation coefficient (Sec.~\ref{S:Pearson}) and (2) the voxel intensity distribution (VID; Sec.~\ref{S:VID}). The former quantifies the phase-space information, whereas the latter captures the one-point statistics that describes the distribution of voxel intensities.

Finally, we present results with a variety of instrument noise levels $\sigma_n$ and the reconstruction threshold $m$. For each test, we use 2500 mock light cones to calculate the correlation coefficient and VID, and estimate errors with 100 noise realizations.

\subsection{Visualization of Example Results}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_im.pdf}
\caption{\label{F:gen_im} Visualization of the phase-space reconstruction with $\sigma_n=10^4$ Jy sr$^{-1}$. We place the 400 toy model light cones into a $20\times 20$ pixel map, and show the reconstruction results on one of the spectrum bins (274 GHz). The top left panel is the true signal intensity map from all of the spectral lines (CO(3--2), CO(4--3), CO(5--4), CO(6--5), [\ion{C}{ii}]). The top right panel is the observed intensity map including line signals (top left panel) and noise. The bottom three panels show the true input and the 10$\sigma$, 5$\sigma$, and 1$\sigma$ reconstructed (right) emission field of the three spectral lines in the multi-line regime where the signal can be reconstructed.}
\end{center}
\end{figure}

Fig.~\ref{F:gen_im} visualizes the reconstruction results on one of the spectrum bins (274 GHz) that contains five lines (CO(3--2), CO(4--3), CO(5--4), CO(6--5), [\ion{C}{ii}]) in the observable band. We place the 400 light cones into a $20\times 20$ pixel map, and show the true and the reconstructed intensity fields. Note that the input signals do not exhibit the spatial clustering because each light cone is generated independently. The top left panel shows the total signal intensity from all five lines, and the top right panel shows the observed intensity including the total line signals and a Gaussian instrument noise with $\sigma_n=10^4$ Jy sr$^{-1}$. In this channel, three of the lines (CO(4--3), CO(5--4), CO(6--5)) are in the multi-line regime, so they can be reconstructed with our algorithm. The three bottom panels compare the true input to the reconstructed intensity maps for these lines with 10$\sigma$, 5$\sigma$, and 1$\sigma$ reconstruction threshold as the MP algorithm stopping criteria.

The choice of threshold is a trade-off between the completeness and the purity in the reconstructed map. As shown in Fig.~\ref{F:gen_im}, 
in the high threshold case (10$\sigma$), the MP algorithm only extracts a few bright sources that are above the reconstruction threshold. As the threshold decreases, more sources are being reconstructed, at the cost of increased false positive detection from the noise fluctuations or the interlopers. This example provides a visual depiction of the reconstruction algorithm. To further quantify the reconstruction performance, we consider two summary statistics in the following sections.

\subsection{Pearson Correlation Coefficient}\label{S:Pearson}
We quantify the reconstruction performance by the Pearson correlation coefficient between the true and the reconstructed maps in each channel. The Pearson correlation coefficient is defined by
\begin{equation}\label{E:pearson}
r = \frac{\sum\limits_{i=1}^{N_{\rm lc}}\left ( I_{\rm true}^i - \left \langle I_{\rm true}  \right \rangle \right )\left ( I_{\rm rec}^i - \left \langle I_{\rm rec}  \right \rangle \right )}{\sqrt{\sum\limits_{i=1}^{N_{\rm lc}}\left ( I_{\rm true}^i - \left \langle I_{\rm true}  \right \rangle \right )^2}\sqrt{\sum\limits_{i=1}^{N_{\rm lc}}\left ( I_{\rm rec}^i - \left \langle I_{\rm rec}  \right \rangle \right )^2}},
\end{equation}
where $N_{\rm lc}=2500$ is the number of light cones, $I_{\rm true}^i$ and $I_{\rm rec}^i$ are the true and the reconstructed line intensity maps, respectively, at the $i$th light cone. Fig.~\ref{F:gen_r} shows the results of the correlation coefficient with a $\sigma_n = 10^4$ Jy sr$^{-1}$ noise level and a 5$\sigma$ reconstruction threshold. Our reconstructed map achieves $\sim 80\%$ correlation with the true input map at $z\lesssim 1.5$.

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r.pdf}
\caption{\label{F:gen_r} Pearson correlation coefficient $r$ between the true and the reconstructed maps on 2500 light cones with $\sigma_n=10^4$ Jy sr$^{-1}$ 5$\sigma$ reconstruction. The bands are the 1$\sigma$ scatter of 100 noise realizations with the sample line signal. The gray bands are the $r$ value with the white noise map for reference.}
\end{center}
\end{figure}


\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_bands.pdf}
\caption{\label{F:gen_r_bands} Average correlation coefficient $r$ in each broad band (Table~\ref{T:bands}) 
 within a range of noise level $\sigma_n$ and the reconstruction threshold $m$. The error bars are the rms of 100 noise realizations.}
\end{center}
\end{figure*}


Fig.~\ref{F:gen_r_bands} shows the correlation coefficient $r$ within a range of noise level $\sigma_n$ and the reconstruction threshold $m$. For simplicity, we show the average $r$ value of each broad band defined in Table~\ref{T:bands}. The key findings are summarized below.

\begin{itemize}
\item At a fixed $m$ value, $r$ decreases as $\sigma_n$ increases. This is because under the same purity (same $m$), the detection threshold $m\sigma_n$ is higher for higher $\sigma_n$, and thus fewer sources have been reconstructed.
\item The six defined bands correspond to three pairs of lines from different redshift bands (Table~\ref{T:bands}). The lines within each pair are strongly correlated because they are the signals from the same sources, and thus they are reconstructed in the same MP iteration. The pairs of lines in the same redshift are reconstructed in the same MP iteration, and as a result, they are highly correlated. 
\item Because of the purity and completeness trade-off, the maximum correlation $r$ happens at the intermediate threshold $m$ (except for the lowest-noise $\sigma_n=10^3$ Jy sr$^{-1}$ case, discussed in the next enumerated point). 
\item Correlation coefficient $r$ has very low dependency on $\sigma_{\rm th}$ in the $\sigma_n=10^3$ Jy sr$^{-1}$ case. This can be understood by comparing the noise level $\sigma_n$ to the quantity $I_*$, the intensity of the $\ell_*$ source in the Schechter function (Fig.~\ref{F:Istr}). At the redshift range in which we perform the reconstruction ($0.5\lesssim z\lesssim 1.9$), $10^4\lesssim I_*\lesssim 10^5$ Jy sr$^{-1}$, which indicates that $m\sigma_n < I_*$ for all of the $m$ values considered in Fig.~\ref{F:gen_r_bands} ($m=1\sim10$). In the Schechter luminosity function, the sources $\gtrsim \ell_*$ contribute the majority of the information in the intensity field \citep{2019ApJ...877...86C}, and thus the correlation coefficient $r$ is not sensitive to the change in reconstruction threshold if $m\sigma_n\ll\ell_*$. 
\end{itemize}


\subsection{VID}\label{S:VID}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_VID.pdf}
\caption{\label{F:gen_VID} VID of the true input (black) and the reconstructed (red) maps in the six broad bands with noise level $\sigma_n=10^4$ Jy sr$^{-1}$ and reconstruction threshold $m=4$.  The error bars are the rms of 100 noise realizations. The gray curves are the VID of the total observed map, which includes signals from all of the lines and noise. For reference, the blue and green dashed lines mark the noise level $\sigma_n$ and $I_*$, respectively.}
\end{center}
\end{figure*}

The Pearson correlation coefficient traces the phase-space variations between the true and reconstructed maps, but it cannot distinguish a systematic constant offset, i.e., if the reconstructed line signals are systematically lower or higher than the true input. Therefore, we check the consistency of the reconstructed and true input maps using the one-point statistics, VID. Note that we do not directly compare the mean intensity of each map since we only reconstruct the bright sources that are above the threshold and neglect all faint sources in the reconstructed map, so the mean intensity is not expected be faithfully recovered.

The VID of an LIM map contains information beyond the power spectrum, and is valuable for LIM targeting a late-time universe where the large-scale structure is highly non-Gaussian and cannot be fully described by two-point statistics. For example, \citet{2017MNRAS.467.2996B} showed that the VID can constrain the luminosity function model parameters, and \citet{2019ApJ...871...75I} demonstrated that a joint analysis of the power spectrum and VID improved the constraining power on the source luminosity function. 

Fig.~\ref{F:gen_VID} compares the VID of the true and reconstructed maps in the six broad bands with noise level $\sigma_n=10^4$ Jy  sr$^{-1}$ and reconstruction threshold $m=4$. The gray lines compose the VID of the total observed map, which includes signals from all of the lines and noise. The black lines are the VID of the target line maps, and the red data points are the VID of the reconstructed target line map.While the VID of the total observed map is one to two orders of magnitude above that of the target line signal, our reconstruction technique can faithfully recover the VID of the signal to slightly below the $I_*$-scale, the characteristic source luminosity in the Schechter function.

In this realistic survey setup with the assumed signal model, we show that our method can successfully reconstruct the VID of the CO signal down to $\sim \ell_*$ scales. This provides a strong constraint on the CO luminosity function, as both are 1D statistics of the intensity field and are closely related. The CO luminosity function at various redshifts provides valuable insight on the formation and evolution of galaxies  across cosmic time. Specifically, CO is a tracer of H$_2$ gas in the interstellar medium, and can therefore be used to study the evolution of the molecular gas content and their distribution as a function of time \citep{2014ApJ...782...79W,2019ApJ...882..138D,2019ApJ...872....7R}. We note though that the expected S/N on the luminosity function depends on the assumed model.

\section{Discussion}\label{S:discussion}
\subsection{Model Uncertainty}\label{S:model_uncertainty}
\begin{figure*}[h!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_bands_var.pdf}
\caption{\label{F:gen_r_bands_var} Comparing $r$ of the no SLED variation (red), 20\% variation (cyan), 50\% variation (yellow), and 100\% variation (purple) with the 4$\sigma$ reconstruction threshold. The values are the average of $r$ within the channels of the band, and the error bars are the rms of 100 noise realizations of all of the spectral bins in each band.}
\includegraphics[width=\linewidth]{gen_VID_var.pdf}
\caption{\label{F:gen_VID_var} VID of the no SLED variation (red), 20\% variation (cyan), 50\% variation (yellow), and 100\% variation (purple), and the true (black) maps in the six broad bands with noise level $\sigma_n=10^4$ Jy sr$^{-1}$ and reconstruction threshold $m=4$.  The error bars are the r.m.s. of 100 noise realizations. The gray curves are the VID of the total observed map, which includes signals from all of the lines and noise. For reference, the blue and green dashed lines mark the noise level $\sigma_n$ and $I_*$.}
\end{center}
\end{figure*}

For the results presented in Sec.~\ref{S:results}, the light cone signals and the dictionary template $\mathbf{A}$ are both generated from the same assumed signal model (Sec.\ref{S: model}). However, in reality, the variation in SLED across galaxies will affect the reconstruction performance. To test how the SLED uncertainties affect the reconstruction, we apply three different SLED model variations and bias levels (at $20\%$, $50\%$, and $100\%$) to the mock data, and run the reconstruction with the same dictionary template $\mathbf{A}$. We detail the definition of variation and bias below. \citet{2015A&A...577A..46D} measured  multiple CO lines of ULIRGs at $z\sim 1.5$, and estimated a $\sim 20\%$ variation on the CO SLED ratio for their sample. Therefore, an assumed $50\%$ or $100\%$ variation can be more extreme than realistic variations.

\subsubsection{SLED model variation}

\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_bands_bias.pdf}
\caption{\label{F:gen_r_bands_bias} Comparison $r$ of the no SLED bias (red), 20\% bias (cyan), 50\% bias (yellow), and 100\% bias (purple) with a 4$\sigma$ reconstruction threshold. The values are the average of $r$ within the channels of the band, and the error bars are the rms of 100 noise realizations of all of the spectral bins in each band.}
\includegraphics[width=\linewidth]{gen_VID_bias.pdf}
\caption{\label{F:gen_VID_bias} VID of no SLED bias (red), 20\% bias (cyan), 50\% bias (yellow), 100\% bias (purple), and the true maps (solid lines) in the six broad band with noise level $\sigma_n=10^4$ Jy /sr and reconstruction threshold $m=4$.  The error bars are the r.m.s. of 100 noise realizations. The grey curves are the VID of the total observed map that includes signals from all the lines and noise. For reference, the blue and green dashed lines mark the noise level $\sigma_n$ and $I_*$.}
\end{center}
\end{figure*}

First, we test the case with SLED model variations. For each line of each source in the mock light cone, we assign a line luminosity $L^{\rm line}$:
\begin{equation}\label{E:SLED_var}
L^{\rm line} = L^{\rm line}_{\rm fid} (1+\delta_L),
\end{equation}
where $L^{\rm line}_{fid}$ is the fiducial luminosity from our model, and $\delta_L$ is a zero-mean Gaussian random variable with the standard deviations of $0.2$, $0.5$, and $1.0$ for the $20\%$, $50\%$, and $100\%$ variation cases, respectively.

Fig.~\ref{F:gen_r_bands_var} shows the reconstructed correlation coefficients $r$ with a 4$\sigma$ reconstruction threshold ($m=4$) with different SLED variations. We see that introducing $20\%$ SLED variations in the model has a very mild impact on the reconstructed correlation; whereas with 50\% (100\%) variation, the correlation coefficient drops by about 10\% ($\sim 50\%$). Fig.~\ref{F:gen_VID_var} compares the VID. We see that the SLED fluctuation has a negligible impact on the VID reconstruction in the 20\%, 50\%, and 100\% variation cases. Hence, we conclude that our technique is robust against realistic level of CO SLED fluctuations.

\subsubsection{Model offset}

In addition to SLED variation, we also test whether our model gives biased estimates of the average SLED by assigning line luminosity $L^{\rm line}$ as
\begin{equation}
L^{\rm line} = L^{\rm line}_{\rm fid} (1+b_L),
\end{equation}
where $b_L$ is a constant offset that we assign to the model lines, and it is applied to all of the sources in the mock light cones. For the 20\% bias level, we apply $b_L=(+0.1,-0.1,+0.1,-0.1)$ for the four CO lines CO(3--2), CO(4--3), CO(5--4), and CO(6--5) to ensure the SLED ratio between neighboring lines is $~20\%$, since our algorithm is only sensitive to the line ratio between two neighboring CO lines.

Fig.~\ref{F:gen_r_bands_bias} and Fig.~\ref{F:gen_VID_bias} show the reconstructed correlation coefficients $r$ and VID with a 4$\sigma$ with three SLED bias levels. Similar to Fig.~\ref{F:gen_r_bands_var}, the correlation only drops significantly when the bias is tuned to 100\%. Therefore, our technique is also robust against a realistic level of potential CO SLED bias.


\subsection{Application to LIM in Other Wavelengths}
The technique developed in this work is not restricted to the [\ion{C}{ii}] and CO lines' blending problem. It can in principle be applied to a range of LIM experimental setups. As a demonstration, we apply our method to reconstruct near-infrared lines in an SPHEREx-like survey.

SPHEREx is an ongoing NASA MIDEX mission to conduct an all-sky near-infrared spectro-imaging survey \citep{2014arXiv1412.4872D}.\footnote{\url{http://spherex.caltech.edu}} SPHEREx will carry out the first all-sky spectral survey at wavelengths between 0.75 and 5 $\mu$m with 96 spectral channels and a ${6{''}.2}$ pixel size. Ly$\alpha$ (121.6 nm), H$\alpha$ (656.3 nm), H$\beta$ (486.1 nm), [\ion{O}{ii}] (372.7 nm), and [\ion{O}{iii}] (500.7 nm) are the five prominent lines detectable by SPHEREx across a range of redshifts especially in the LIM regime.

The line signal model is described in Appendix \ref{A:SPHEREx}. We generate a near-infrared LIM mock data with a ${6{''}.2}\times {6{''}.2}$ pixel size and a  $5\sigma$ point-source sensitivity of $m_{\rm AB}=22$ (similar depth as the SPHEREx deep fields), and run our de-confusion algorithm on the mock light cones. Fig.~\ref{F:gen_r_spherex} shows the results of the correlation coefficients between the true and a 3$\sigma$-threshold reconstructed intensity maps. The reconstructed map achieves $\sim 80\%$ correlation with the true input map at $z\lesssim 3$, and decreases toward higher redshifts, as the $I_*$ of the lines approaches the noise level $\sigma_n$. At $z\sim 5$, the brightest line, H$\alpha$, is redshifted into SPHEREx bands with a high spectral resolution of $R\sim 130$ and suffers less signal dilation, resulting in an increase of S/N on a single source detection, and thus $r$ slightly rebounds at this redshift.

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_spherex.pdf}
\caption{\label{F:gen_r_spherex} Pearson correlation coefficient $r$ between the true and the reconstructed maps on 2500 light cones for SPHEREx-like mock data. The bands are the rms. of the value in 100 noise realizations with the sample line signal. The gray bands are the correlation coefficient with the uncorrelated white noise map for reference.}
\end{center}
\end{figure}

To account for uncertainties in modeling the SLED, we apply a realistic level of line luminosity variations from \citet{2006ApJ...642..775M}. We apply 10\%, 50\%, and 100\% SLED variations (Eq.\ref{E:SLED_var}) to the line luminosity ratio of H$_\alpha$/H$\beta$, [\ion{O}{ii}]/H$\alpha$, and [\ion{O}{iii}]/[\ion{O}{ii}], respectively. Fig.~\ref{F:gen_r_spherex_var} shows the correlation coefficients between the true and a 3$\sigma$-threshold reconstructed intensity maps. Comparing with the fixed SLED case (Fig.~\ref{F:gen_r_spherex}), only the $H\beta$ line shows a significant decrease in the performance. The intensity map of the three brighter lines (H$\alpha$, [\ion{O}{ii}], and [\ion{O}{iii}]) can still be extracted with $\gtrsim 70\%$ correlation compared to the true input.

We conclude that our algorithm can reasonably well reconstruct the phase-space LIM signal in a SPHEREx-like experiment, given the expected variation of (redshifted) optical line ratios. The technique can be generalized to  different LIM experimental applications, and the reconstructions are fairly robust against uncertainties in the SLED modeling.

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{SPHEREx_r_var.pdf}
\caption{\label{F:gen_r_spherex_var} Pearson correlation coefficient $r$ between the true and the reconstructed maps on 2500 light cones for SPHEREx-like mock data with a realistic level of SLED variation. The bands are the rms of the value in 100 noise realizations with the sample line signal. The gray bands are the correlation coefficient with the uncorrelated white noise map for reference.}
\end{center}
\end{figure}

\subsection{Foreground Subtraction}\label{S:bg}
\begin{figure*}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_bands_bg.pdf}
\caption{\label{F:gen_r_bands_bg}Comparing $r$ of the reconstructed map with the no background subtraction (red), subtracting the mean of the whole data cube (cyan), and subtracting a first-order polynomial in the spectral direction for each light cone (yellow), using a 4$\sigma$ reconstruction threshold. The values are the average of $r$ within the channels of the band, and the error bars are the rms of 100 noise realizations of all of the spectral bins in each band.}
\includegraphics[width=\linewidth]{gen_VID_bg.pdf}
\caption{\label{F:gen_VID_bg} VID of no background subtraction (red),subtracting the mean of the whole data cube (cyan), subtracting a first order polynomial in spectral direction for each light cone (yellow), and the true (black) maps in the six broad band with noise level $\sigma_n=10^4$ Jy /sr and reconstruction threshold $m=4$.  The error bars are the r.m.s. of 100 noise realizations. The grey curves are the VID of the total observed map that includes signals from all the lines and noise. For reference, the blue and green dashed lines mark the noise level $\sigma_n$ and $I_*$.}
\end{center}
\end{figure*}

In addition to line interlopers, LIM data are subjected to strong continuum foreground from various sources. For the frequency range considered in this work ($\sim$ 200-300 GHz), the dominant foregrounds are the atmospheric emission, dust continuum and the CMB; whereas for LIM in the near infrared, e.g. SPHEREx \citep{2014arXiv1412.4872D}, the zodiacal light and the galaxy stellar continuum are the dominant  continuum foregrounds. Even though these  foregrounds are brighter than the sought-after line signals, their spectral responses are expected to be smooth and are distinct from the spectral line features, so that the continuum foregrounds can be separated and mitigated, for example by a smooth function fit such as a low order polynomial\footnote{Some foreground components are also spatially smooth (e.g. zodiacal light) that can be filtered in the spatial domain as well.}. Here we test how the foreground mitigation process affects our line reconstruction results.

We consider two cases of foreground mitigation. First we emulate the foreground removal process in the presence of an approximately constant foreground in both the spatial and spectral dimensions, for example the zodiacal light. In this case, before running the reconstruction, we subtract the mean value of the whole data cube, i.e. the mean intensity in $N_{\rm lc}\times N_{\nu}$ voxels. The second case is to emulate the continuum subtraction process of the galaxy stellar or dust continuum, which are expected to have smooth spectra but different in each light cone, since each light cone contain different galaxies with different continuum spectrum. We fit and subtract a first order polynomial function to the spectrum of each light cone before running the reconstruction.

The results of a 4-$\sigma$-threshold reconstruction with different noise $\sigma_n$ level are shown in Fig.~\ref{F:gen_r_bands_bg} and Fig.~\ref{F:gen_VID_bg}. We can see compared to the no background subtraction case, the $r$ value is even higher in these two tests. This is because our reconstruction only extract the bright lines, and the fainter lines act as a background for the MP algorithm. The signals from the fainter lines introduce not only fluctuations but also a bias in the data, since the line signals are always positive unlike the zero-mean noise. The reconstruction performance is improved after background subtraction because this bias level is also removed during this process. For the VID results, we see that there is no significant difference compared to the no background subtraction case. This is again due to the fact that the background level is much fainter than the brightness of the sources being extracted with our algorithm, so the background subtraction have no impact on the reconstruction. In conclusion, the background subtraction in the LIM data reduction pipeline will not affect our line reconstruction technique.

\subsection{Prior with External Catalogs}\label{S:ext_cats}
\begin{figure*}[ht!]
\includegraphics[width=\linewidth]{gen_r_bands_prior.pdf}
\caption{\label{F:gen_r_bands_prior}Comparing $r$ of the reconstructed map with no external catalog prior (red), and utilizing the external catalog prior to fit the catalog sources before running the MP reconstruction (cyan), using 4-$\sigma$ reconstruction threshold. The values are the average of $r$ within the channels of the band, and the error bars are the r.m.s. of 100 noise realizations of all the spectral bins in each band.}
\end{figure*}

Our analysis uses the LIM data itself without invoking any external information. In practice, initial LIM survey fields are designed to in part overlap with existing photometric or spectroscopic galaxy surveys, and thus there will be information provided by external galaxy catalogs to aid the line de-confusion problem. Ignoring redshift uncertainties of the external catalogs, one simple approach to incorporate the external information is to force the MP algorithm to first select the redshift bins that contain galaxies from the catalog. After  iterating through the catalog sources, we then continue the normal MP procedure until hitting the stopping criteria. 

To test the effect of including prior knowledge from external catalogs, we generate a mock catalog by selecting sources with CO(5--4) flux greater than $150\, L_\odot$ Mpc$^{-2}$ ($6\times 10^{-17}$ W m$^{-2}$) in each light cone. The source density in the catalog is $\sim 1.2$ per light cone (integrated along line of sight) for a $0.43^2$ arcmin$^2$ pixel solid angle. The flux cut corresponds to an $L_*$ galaxy at $z\sim 2$. According to \citet{2012ApJ...752..113H}, such a $L_*$ galaxy has an absolute magnitude $M_{\rm AB}\sim -23$ in the optical, which gives an apparent magnitude of $m_{\rm AB}\sim 21.8$, approximately the depth of the assumed optical catalog.

With this mock external catalog, we identify the redshift bins containing the catalog sources, regardless of the noise level and threshold value. After projecting out these components, we run the MP algorithm on the residual data as per usual until we hit the stopping criteria. The results of a 4-$\sigma$-threshold reconstruction with different noise levels are shown in Fig.~\ref{F:gen_r_bands_prior}. The reconstruction shows improved results  for all noise levels. However, the huge improvement in the highest-noise case cannot be interpreted as a successful reconstruction of source intensities. At this high noise level, the signals are well buried under the noise, and when we fit the data with the catalog source redshift templates, the extracted components are dominated by noise rather than signal amplitude. Thus the improved correlation is merely due to the position information imposed by the external catalogs.

\subsection{Comparing with the Limit of No Interlopers}
\begin{figure}[htbp!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_txtn.pdf}
\caption{\label{F:gen_r_txtn} Pearson correlation coefficient $r$ between the true line maps and the same maps adding the $\sigma_n=10^{-4}$ Jy sr$^{-1}$ noise on 2500 light cones. The bands are the rms of the value in 100 noise realizations with the sample line signal. The gray bands are the $r$ value shown with the uncorrelated white noise map for reference.}
\end{center}
\end{figure}

In the case of no interloper lines (i.e. only single line emission in the data), the best estimator of the single-line intensity map is the observed map (regardless of foregrounds). To compare our reconstruction performance with this limiting case, we calculate the correlation coefficient between the input single-line map with the same single-line map plus the instrument noise.

Fig.~\ref{F:gen_r_txtn} shows the results using $\sigma_n=10^4$ Jy sr$^{-1}$.  Comparing to Fig.~\ref{F:gen_r}, which has the same noise level, the single-line-plus-noise case has a lower correlation than the reconstruction in Fig.~\ref{F:gen_r}. Especially for the fainter lines (e.g., CO(3--2) at $z\sim 0.6$ and CO(4--3) at $z\sim1.1$), our reconstruction map has a much better correlation. This can be explained by the fact that in our algorithm, the sources are detected in the template space rather than in a single voxel. That is, if a source can be observed in two frequency channels, we extract the source by projecting the signals in these two channels to the template space, which is effectively combining the information from both channels. Consequently, we are able to achieve a better S/N on the fainter lines because of the greater sensitivity of their brighter counterpart.

\subsection{Improving Cross-correlation Uncertainty}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{SNRxcorr.pdf}
\caption{\label{F:SNRxcorr} Comparing the S/N on the cross-power spectrum of one CO line and the external tracer, using the total observed map ($S/N^{\rm tot}$) and the reconstructed CO map ($S/N^{\rm r}$). We consider two different values of $\delta r_{\rm co,co^{r}}$ and $r_{\rm co,tot}$ that covers the range of realistic parameter values in our model. The value of $r_{\rm co,co^{r}}\sim 0.7-0.9$ according to the results in Sec~\ref{S:results}, and therefore using the reconstructed map instead of the total observed map in cross-correlation can reduce the uncertainty (i.e. $SNR^{\rm r} > SNR^{\rm tot}$).}
\end{center}
\end{figure}

\begin{figure}[ht!]
\begin{center}
\includegraphics[width=\linewidth]{gen_r_all.pdf}
\caption{\label{F:gen_r_all} Pearson correlation coefficient $r$ between the true input and the total observed maps on 2500 light cones with $\sigma_n=10^4$ Jy sr$^{-1}$. The bands are the rms of the value in 100 noise realizations with the sample line signal. The gray bands are the $r$ value shown with the uncorrelated white noise map for reference.}
\end{center}
\end{figure}

Another useful application of the reconstruction technique is for a more precise measurement of cross-correlation between LIM and other tracers. Cross-correlation analysis not only serves as a validation of a cosmological signal in LIM, since the cross-correlation is less susceptible to foreground contamination and other systematic effects, but also provides valuable astrophysical and cosmological information. 

While cross-correlating the total observed LIM map with an external tracer gives an unbiased estimator of cross spectrum between the target line and the external tracer, the presence of continuum foregrounds and interlopers increases the error of this measurement. Here we present a simple argument that for cross-correlation analysis, using a reconstructed map instead of the total observed map can effectively reduce the error bars on the cross-power spectrum. 

The cross-power spectrum errors between two fields $\delta P_{12}$ in a single mode are given by
\begin{equation}\label{E:xerr}
\delta P_{1,2}^2 = \frac{1}{2}\left ( P_{1,2}^2 + \delta P_1\delta P_2 \right ),
\end{equation}
where $P_{1,2}$ is the cross-power spectrum, and $\delta P_1$ and $\delta P_2$ are the errors on the auto power spectrum in the two fields. For a single $k$ mode, $\delta P_1=P_1$, $\delta P_2=P_2$, where $P_1$ and $P_2$ are the total power spectra (including signals and noise) measured in two fields. Eq.~\ref{E:xerr} can be expressed in terms of the cross-correlation coefficient: $r_{1,2}=P_{1,2}/\sqrt{P_1P_2}$,
\begin{equation}\label{E:xerr2}
\delta P_{1,2}^2 = \frac{1}{2}P_{1,2}^2\left (1+\frac{1}{r_{1,2}^2}\right ),
\end{equation}

Say we have an external galaxy sample that traces one of the target CO lines, then we can write the total observed LIM data as the combination of target CO line ($I_{\rm CO}$), other interloper lines ($I_{\rm interlopers}$), and the noise ($I_n$),
\begin{equation}
I_{\rm tot} = I_{\rm CO} +I_{\rm interlopers} +I_n.
\end{equation}
The expectation value of the cross spectrum between the observed total map and galaxy is the same as the cross spectrum with only the target CO, since the other components are not correlated with the large-scale structure at the same redshift, so $\left \langle P_{\rm g,tot} \right \rangle=\left \langle P_{\rm g,co} \right \rangle$. 

For simplicity, we assume the galaxies are perfectly correlated with the target CO line field on the scale of interest, so $r_{\rm g,co}=1$. This implies the galaxy field and CO field always have the same correlation $r$ with any given field $x$, $r_{\rm g,x}=r_{\rm co,x}$.

From Eq.~\ref{E:xerr2}, the error on the galaxy-CO cross spectrum measured by cross-correlating the galaxy field with the total observed LIM data is
\begin{equation}
\begin{split}
\delta P_{\rm g,co}^{\rm tot} &= P_{\rm g,tot}\sqrt{\frac{1}{2}\left(1+\frac{1}{r_{\rm g,tot}^2}\right)}\\
&= P_{\rm g,co}\sqrt{\frac{1}{2}\left(1+\frac{1}{r_{\rm co,tot}^2}\right)}
\end{split}
\end{equation}

On the other hand, if we cross-correlate the galaxies with the reconstructed CO map (${\rm CO^r}$), the error is
\begin{equation}
\delta P_{\rm g,co{^r}} = P_{\rm g,co{^r}}\sqrt{\frac{1}{2}\left(1+\frac{1}{r_{\rm g,co{^r}}^2}\right)}.
\end{equation}
However, $P_{\rm g,co{^r}}$ is a biased estimator of $P_{\rm g,co}$ because of the error in the reconstruction. If we assume that the reconstructed map roughly preserves the same power as the true map, $P_{\rm co^{r}}\approx P_{\rm co}$, then we can write\footnote{$P_{\rm g,co{^r}}=r_{\rm g,co{^r}}\sqrt{P_gP_{\rm co^{r}}}\approx r_{\rm g,co{^r}}\sqrt{P_gP_{\rm co}}=r_{\rm g,co{^r}}P_{\rm g-co}$, where the last equality uses the assumption $r_{\rm g,co}=1$.} $P_{\rm g,co{^r}}=r_{\rm g,co{^r}}P_{\rm g,co}=r_{\rm co,co{^r}}P_{\rm g,co}$. Therefore, we have to de-bias cross spectrum $P_{\rm g,co{^r}}$ by factor $1/r_{\rm co,co{^r}}$, $P_{\rm g,co}=P_{\rm g,co{^r}}/r_{\rm co,co{^r}}$. The value of $r_{\rm co,co{^r}}$ cannot be directly inferred from the data, so we have to estimate it by simulating the possible range of signals; this introduces an extra error term to the $\delta r_{\rm co,co{^r}}$ due to the uncertainty in $r_{\rm co,co{^r}}$,
\begin{equation}
\begin{split}
\delta P_{\rm g,co}^{\rm r} &= P_{\rm g,co}\sqrt{\left ( \frac{\delta P _{\rm g,co^{r}}}{P_{\rm g,co^{r}}} \right )^2 + \left ( \frac{\delta r _{\rm co,co^{r}}}{r_{\rm co,co^{r}}} \right )^2}\\
&=P_{\rm g,co}\sqrt{\frac{1}{2}\left ( 1+\frac{1}{r^2_{\rm co,co^{r}}} \right ) + \left ( \frac{\delta r _{\rm co,co^{r}}}{r_{\rm co,co^{r}}} \right )^2}.
\end{split}
\end{equation}

In summary, the S/N on the galaxy-CO cross spectrum $P_{\rm g,co}$ using the total observed map and the reconstructed map is
\begin{equation}
\begin{split}
SNR^{\rm tot}=\frac{P_{\rm g,co}}{\delta P_{\rm g,co}^{\rm tot}}&=\frac{1}{\sqrt{\frac{1}{2}\left(1+\frac{1}{r_{\rm co,tot}^2}\right)}},\\
SNR^{\rm r}=\frac{P_{\rm g,co}}{\delta P_{\rm g,co}^{\rm r}}&=\frac{1}{\sqrt{\frac{1}{2}\left ( 1+\frac{1}{r^2_{\rm co,co^{r}}} \right ) + \left ( \frac{\delta r _{\rm co,co^{r}}}{r_{\rm co,co^{r}}} \right )^2}}.
\end{split}
\end{equation}

Note that the S/N in both cases converges to unity when $r$ is unity and the de-bias error is zero, which is the limit of sample variance\footnote{The power spectrum cross-correlation coefficient can be derived from the Pearson correlation (Eq.~\ref{E:pearson}) with a weighting on pixels. Since our pixels are generated and reconstructed independently of each other, the Pearson correlation coefficient here is an unbiased estimator of the power spectrum correlation coefficient. Thus here we will use the value of Pearson correlation we derived for the power spectrum correlation coefficient.}. 

Fig.~\ref{F:SNRxcorr} shows the $S/N^{\rm tot}$ and $S/N^{\rm r}$ with two different values of $\delta r_{\rm co,co^{r}}$ and $r_{\rm co,tot}$ as a function of $r_{\rm co,co^{r}}$. According to the calculations in Sec.~\ref{S:model_uncertainty}, the value of $r_{\rm co,co^{r}}$ ranges from $\sim 0.7$ to $\sim 0.9$ (except for the most noisy case, $\sigma_n=5\times 10^4$ Jy sr$^{-1}$). To estimate the realistic $r_{\rm co,tot}$ value, we calculate the correlation of the input CO line maps and the observed map with $\sigma_n=10^4$ Jy sr$^{-1}$, the same noise level as in Fig.~\ref{F:gen_r}. The results are shown in Fig.~\ref{F:gen_r_all}. We find that $r_{\rm co,tot}$ are around 0.2-0.5 in this case. With this range of parameters, Fig.~\ref{F:SNRxcorr} indicates that $S/N^{\rm r}$ is better than $S/N^{\rm tot}$, which means that in our model, using the reconstructed map instead of the observed map in cross-correlation can reduce the uncertainty. 

\subsection{Estimating the complete interloper population}
Here we discuss a potential extension of the technique, which is capable of estimating the {\it total} power of an interloper population.  The method uses an incompletely reconstructed interloper sample and an incomplete, external  tracer of the interloper density field. Once the interloper contribution is fully quantified, the high-redshift signal of interest in an LIM dataset can be estimated without bias. As an example, we write the observed intensity of a [\ion{C}{ii}] LIM dataset as 
\begin{equation}
    I^{\rm obs} = I_{\rm CO} + I_{\rm C\,II} + \delta_n^{\rm obs},
\end{equation}
where $I_{\rm CO}$ is the total  CO interloper intensity, $I_{\rm C\,II}$ is the \ion{C}{ii} signal, and $\delta_n^{\rm obs}$ is the instrumental noise. For simplicity, here we only consider the contribution of one CO rotational line as the foreground.

Given a reconstruction threshold, we reconstruct the bright CO emissions using our technique as:
\begin{equation}
    I^{\rm rec}_{\rm CO} = \alpha I_{\rm CO_b} + \delta_n^{\rm rec}, 
\end{equation}
where $I^{\rm rec}_{\rm CO}$ is the reconstructed CO intensity and $I_{\rm CO_b}$ is the intensity of bright CO sources in the reconstruction, which is a subset of the total CO population. There are two sources of error in the reconstruction: a multiplicative term $\alpha$ proportional to the bright CO intensity, where $\alpha$ can be greater or smaller than unity, and an additive term $\delta_n^{\rm rec}$ describing random (or misidentified) fluctuations about the true CO intensity, which is uncorrelated with the CO field. 

The auto power spectrum of the reconstructed CO map is
\begin{equation}
\begin{split}
    \left \langle  I_{\rm CO}^{\rm rec}I_{\rm CO}^{\rm rec} \right \rangle =& \bar{\alpha}^2 \left \langle  I_{\rm CO_b}I_{\rm CO_b} \right \rangle  + \left \langle \delta_n^{\rm rec}\delta_n^{\rm rec} \right \rangle \\
    =& \bar{\alpha}^2  \left \langle  I_{\rm CO_b} \right \rangle^2 b_{\rm CO_b}^2 P(k) + N^{\rm rec},
\end{split}
\end{equation}
where, in the linear regime that we consider here, $b_{\rm CO_b}$ is the cosmological clustering bias of the bright CO sources, $P(k)$ is the matter density field, and $N^{\rm rec}$ is the auto power spectrum of the $\delta_n^{\rm rec}$ term as a noise bias.

Using an external galaxy sample $g$ in the same redshift range as the CO interlopers, we cross-correlate $g$ with the observed and reconstructed maps, respectively, and consider only linear clustering scales:
\begin{equation}
    \left \langle g I_{\rm obs} \right \rangle  = r \left \langle  I_{\rm CO} \right \rangle b_{\rm CO} b_{\rm g} P(k),
\end{equation}
\begin{equation}
    \left \langle g I_{\rm CO}^{\rm rec} \right \rangle = \left \langle g \alpha I_{\rm CO_b} \right \rangle = \bar{\alpha} r_b \left \langle  I_{\rm CO_b} \right \rangle b_{\rm CO_b} b_{\rm g} P(k), 
\end{equation}
where $\langle \alpha \rangle = \bar{\alpha}$, $b_g$ is the bias of the galaxy tracer, $b_{\rm CO_b}$ is the bias of the bright CO population, and $r$ ($r_b$) is the astrophysical stochastic cross-correlation parameter between the galaxy and CO (CO$_b$) populations.  

Finally, we construct an estimator of the full CO power spectrum as
\begin{equation}\label{E:Pco_est}
\begin{split}
   &\left \langle  \widehat{I_{\rm CO} I_{\rm CO}} \right \rangle = \left \langle  I_{\rm CO}^{\rm rec}I_{\rm CO}^{\rm rec} \right \rangle\left (  \frac{\left \langle g I_{\rm CO}^{\rm obs} \right \rangle}{\left \langle g I_{\rm CO}^{\rm rec} \right \rangle} \right )^2 \\
   &= \left (\bar{\alpha}^2 \left \langle  I_{\rm CO_b} \right \rangle^2 b_{\rm CO_b}^2 P(k) + N^{\rm rec}  \right ) \left ( \frac{r\left \langle I_{\rm CO}\right \rangle b_{CO}}{\bar{\alpha} r_b\left \langle I_{\rm CO_b}\right \rangle b_{CO_b} } \right )^2 \\
   &= \left ( \left \langle  I_{\rm CO_b} \right \rangle^2 b_{\rm CO_b}^2 P(k) + \frac{N^{\rm rec}}{\bar{\alpha}^2}  \right ) \left ( \frac{r\left \langle I_{\rm CO}\right \rangle b_{CO}}{ r_b\left \langle I_{\rm CO_b}\right \rangle b_{CO_b} } \right )^2 \\
\end{split}
\end{equation}

If we can pick a reconstruction threshold such that the reconstruction of a subset of the CO population (bright CO sources) includes the majority of CO emitters, then $r \approx r_b$. Furthermore, if the reconstruction noise is negligible, $N^{\rm rec}/\bar{\alpha}^2 \rightarrow 0$, then this estimator becomes an  unbiased estimator of the full CO power spectrum on large scales:
\begin{equation}
    \left \langle  \widehat{I_{\rm CO} I_{\rm CO}} \right \rangle =
   \left \langle  I_{\rm CO}^{\rm rec}I_{\rm CO}^{\rm rec} \right \rangle\left (  \frac{\left \langle g I_{\rm CO}^{\rm obs} \right \rangle}{\left \langle g I_{\rm CO}^{\rm rec} \right \rangle} \right )^2 = \left \langle  I_{\rm CO} \right \rangle^2 b_{\rm CO}^2 P(k).
\end{equation}
In principle, this argument holds regardless of the luminosity limit of the CO population used for the reconstruction {\it and} regardless of the magnitude limit of the galaxy sample used for the cross-correlation estimate, as long as the galaxy and CO samples overlap spatially. This is potentially a powerful approach to access the entire interloper CO population without the need to identify  the faint, undetected source contributions. Subtracting this CO estimate then provides an unbiased estimate of the high-redshift \ion{C}{ii} power spectrum in the same LIM dataset, which is highly desirable. 

With an external galaxy catalog along, one can also estimate the total CO power on  large scales where $r\approx 1$, with the estimator: $\left \langle g I_{\rm CO}^{\rm obs} \right \rangle^2 / \left \langle g g \right \rangle$. However, the advantage of the  Eq.~\ref{E:Pco_est} estimator is that it only requires $r \approx r_b$, which can be valid on smaller scales where $r < 1$, as long as the reconstruction has sufficient quality in terms of the noise level purity
and completeness. Therefore, with the reconstructed map, one can potentially extract the total CO power to smaller scales.

We note that a high reconstruction threshold $m$, is required in order to achieve a (nearly) noiseless reconstruction ($N^{\rm rec}/\bar{\alpha}^2 \rightarrow 0$). On the other hand, $r \approx r_b$ becomes invalid if the threshold is too high such that majority of the fainter sources is missed in the reconstruction. In addition, a high threshold tends to boost the shot noise in the reconstructed  power spectrum, as well as uncertainties in the CO power spectrum estimator (Eq.~\ref{E:Pco_est}). There is clearly a trade-off between the fidelity of the reconstructed signals and the uncertainty to estimate the desired signals. A detailed simulation is necessary to determine the optimal threshold to minimize the effect of bias from $N_{\rm rec}$ and the variance from the reconstruction shot noise, and evaluate the performance of the CO power spectrum estimator. We leave this investigation to future work.

\subsection{Comparing Foreground Cleaning Capability with Masking}
\begin{figure*}[ht!]
\includegraphics[width=\linewidth]{res_mask.pdf}
\caption{\label{F:res_mask} S/N on [\ion{C}{ii}] shot-noise power spectrum after CO line foreground removal with masking (black dotted) and sparse reconstruction (red) with different instrument noise level $\sigma_n$ and reconstruction threshold $m$. The orange dashed line marks the $m=4$ threshold for reference. We also show the S/N before cleaning (black dashed), and the limiting case of the no CO signals (green) for comparison.}
\end{figure*}

Our algorithm can serve as a foreground mitigation method for \ion{C}{ii} LIM measurement by identifying the bright CO foreground signals. In this section, we quantify this "foreground cleaning" performance using our algorithm and compare it with the masking method, where a ``cleaned map'' is obtained by masking out voxels that contain bright CO sources identified with an external source catalog. 

In the following, we compare two cases: (1) in the masking case, the ``cleaned map'' is the observed map (including all of the lines and instrument noise) masked using external CO catalogs; and (2) the ``cleaned map'' derived from our algorithm is the observed map subtracting the CO reconstructed map. 

The external catalog considered here is the same as the one described in Sec.~\ref{S:ext_cats}: a flux cut on CO(5--4) at the level of $150\, L_\odot$ Mpc$^{-2}$ ($6\times 10^{-17}$ W m$^{-2}$), which gives $\sim 1.2$ galaxies per light cone ($0.43^2$ arcmin$^2$ pixel solid angle) and corresponds to a $m_{\rm AB}\sim 21.8$ threshold in the optical band. Note that this masking threshold is comparable to the ``case A'' masking in \citet{2018ApJ...856..107S}, although here we consider a simpler model that ignores the scatter in the line luminosity model. For this masking scenario, the ``cleaned map'' is the observed map minus any voxels that contain the sources in the external catalog.

In our algorithm, the CO sources and their spectra are reconstructed and removed from the data iteratively, and the residual can be regarded as a ``cleaned map'' that is free of bright CO sources. However, we are only capable of cleaning the multi-line redshift bins in our algorithm. The signal identified in the single-line redshift bins in our algorithm is the combination of the remaining CO and \ion{C}{ii} signals. If we remove all of the reconstructed single-line signals, we will over-subtract \ion{C}{ii} in the cleaned map. Therefore, for the single-line redshift bins ($0.15\leq z\leq0.51$; $0.72\leq z\leq0.89$; $ z\geq1.87$), we clean the data by masking  voxels that contain external catalog sources.

As a figure of merit, we calculate the S/N of the respective \ion{C}{ii} shot-noise power spectra (since there is no clustering signal in our mock light cones) while including the residual CO as part of the noise contribution. 

The error on the shot-noise power spectrum $\delta P_{\rm sh}$ in a map is
\begin{equation}
\delta P_{\rm sh} = \frac{P_{\rm tot}}{\sqrt{N_{\rm mode}}}=\frac{V_{\rm vox}\sigma^2_{\rm tot}}{\sqrt{N_{\rm mode}}},
\end{equation}
where $P_{\rm tot}$ is the total power spectrum of the cleaned map on the shot-noise scales, which is proportional to the total voxel variance in the map, $\sigma^2_{\rm tot}$. $N_{\rm mode}$ is the number of k-space modes used to measure the shot noise. $N_{\rm mode}$ is usually of the order of the total number of voxels, so we choose $N_{\rm mode}=6000$, similar to the number of voxels in TIME. The shot-noise power of the \ion{C}{ii} signal is given by
\begin{equation}
 P_{\rm CII}=V_{\rm vox}\sigma^2_{\rm CII},
\end{equation}
where $\sigma^2_{\rm CII}$ is the voxel variance of the \ion{C}{ii} signal map. The S/N of the \ion{C}{ii} shot-noise power spectrum is then:
\begin{equation}\label{E:CII_SNR}
 \frac{P_{\rm CII}}{\delta P_{\rm sh}}=\frac{\sigma^2_{\rm CII}}{\sigma^2_{\rm tot}}\sqrt{N_{\rm mode}}.
\end{equation}

We test the foreground cleaning performance with the same set of 2500 mock light cones described in Sec.~\ref{S:results}, and calculate the \ion{C}{ii} shot-noise S/N using Eq.~\ref{E:CII_SNR}. The $\sigma^2_{\rm CII}$ is the variance of the input \ion{C}{ii} map, and $\sigma^2_{\rm tot}$ is the variance of the cleaned map. For simplicity, $\sigma^2_{\rm CII}$ and $\sigma^2_{\rm tot}$ are the variance from all of the frequency channels.

In Fig.~\ref{F:res_mask}, the black dotted lines and the red lines show the \ion{C}{ii} shot-noise S/N with the cleaned map obtained from masking and from our algorithm, respectively. For reference, the black dashed lines are the \ion{C}{ii} shot-noise S/N of the observed map before cleaning, and thus the map includes the contribution from all of the lines and the instrument noise.

According to Fig.~\ref{F:gen_r_bands}, the optimal threshold that gives the maximum $r$ value is $m\sim4$. In Fig.~\ref{F:res_mask}, $m=4$ is marked with an orange dashed line, and we see that in the highest-noise case for this threshold, masking (black dotted line) using an external (deep) catalog  performs slightly better than reconstruction (red line) because it is difficult for the MP algorithm to extract the signals from noisy data directly. For realistic noise levels (between $\sigma_n=5\times10^3$ and $10^4$  Jy sr$^{-1}$), the reconstruction outperforms masking.

We also compare this result with the limiting case where there are only \ion{C}{ii} and instrumental noise in the data (green line), i.e. $\sigma^2_{\rm tot} = \sigma^2_{\rm CII} + \sigma^2_n$. For a  small reconstruction threshold $m$, the reconstructed S/N is better than this limit, which indicates that the reconstruction over-fits and misidentifies noise fluctuation as a signal and removes them from the cleaned data. We see that for the two realistic noise levels at $m\sim4$, the reconstructed S/N is lower than this limit, and thus indicates that overfitting is not an issue at this threshold. Also note that in the highest-noise case, the masking SNR is close to the noise-plus-\ion{C}{ii} limit (green line), for the following reason: since the external catalog goes much deeper than the noise level, noise fluctuation dominates over the line signals after masking, and thus variance in the masked map is close to the noise variance.

Finally, we point out that the reconstructed \ion{C}{ii} shot-noise S/N (red line) converges to a constant instead of increasing with smaller threshold $m$ values. This is because of the fact that in the low threshold limit, the reconstruction residual in the cleaned map is subdominant compared to the (masked) single-line redshift bin signals being added back to the residual, and thus the cleaned map S/N does not depend on the threshold value. 

To sum up, for a realistic noise level, our reconstruction performs better than masking in terms of foreground cleaning capability, given our signal model and the external catalog considered in this work. We note that this conclusion depends on the line luminosity function model and the depth of the external catalog for masking, and we leave a more detailed analysis to future work.

\subsection{Technique Extensions}
In this section, we outline some directions for extending the current framework to further improve upon the line reconstruction in future work. 

\subsubsection{Template Generalization}
Currently we use a single SLED template for all redshift bins in the reconstruction; this can be easily generalized to incorporate multiple spectral templates to account for the redshift evolution and SLED variation of the signals. In addition, the  extension can help differentiate the emission from different types of galaxies that have different SLEDs. For instance, if we have two different SLED models for early- and late-type galaxies, respectively, we can incorporate them by having two columns in $\mathbf{A}$ for every redshift bin such that the reconstruction can infer not only the redshift and luminosity of the sources, but also their galaxy type from the SLED templates. 

\subsubsection{Alternative Sparse Approximation Algorithm}
The MP algorithm adopted in this work optimizes the $\ell_0$ norm in Eq.~\ref{E:sparse_eqn}, which is the direct sum of the nonzero elements in $\widetilde{\mathbf{N}}$. We can improve the algorithm by including prior information on the expected value of each element in $\widetilde{\mathbf{N}}$, which depends on the voxel size and source luminosity function. For instance, instead of using the MP algorithm, one can obtain the sparse solution by solving the following  $\ell_1$-norm regularization equation:
\begin{equation}\label{E:L1_reg}
\underset{\widetilde{\mathbf{N}} }{\mathrm{argmin}}\left \| \mathbf{I} - \mathbf{A}\widetilde{\mathbf{N}} \right \|^2_2+\lambda\left \| \mathbf{w} \cdot\widetilde{\mathbf{N}} \right \|_1,
\end{equation}
where the parameter $\lambda$ determines the regularization strength for preventing overfitting, which has a similar effect as the stopping criteria in the MP algorithm; the prior information on the number density of the sources in each redshift bin can be encoded in the weight vector $\mathbf{w}$ in this expression.

\subsubsection{Clustering Information}
In this work, we only perform the pixel-by-pixel line de-confusion using the information in the spectral correlation due to the multiple lines emitted from the same source to reconstruct the signals. The clustering information of the galaxies, which is neglected in this work, could provide additional information on the emission field. We can generalize this framework by incorporating the clustering information from the known galaxy two-point correlation, and perform the reconstruction on an ensemble of pixels to simultaneously fit for the spectral correlation and clustering. For example, if we have a theoretical model for the line-of-sight two-point correlation function of $\mathbf{N}(z)$, $\xi_{\rm th}(\mathbf{N}(z))$, we can add another $\ell_2$-norm regularization term to Eq.~\ref{E:L1_reg},
\begin{equation}
    \lambda_{\rm clus}\left \| \xi(\widetilde{\mathbf{N}}(z)) - \xi_{\rm th}(\mathbf{N}(z))  \right \|_2^2.
\end{equation}
This will enforce the algorithm to give higher priority to solutions close to the theoretical correlation function. Similarly, with an external catalog that traces the same large-scale structure, one can also constrain the algorithm with cross-correlation:
\begin{equation}
    \lambda_{\rm clus}^{x}\left \| \xi^x(\widetilde{\mathbf{N}}(z), \delta_{\rm ext}) - \xi^x_{\rm th}(\mathbf{N}(z), \delta_{\rm ext})  \right \|_2^2,
\end{equation}
where $\delta_{\rm ext}$ is the density field of the external tracer, e.g. galaxy samples, and $\xi^x$ is the cross-correlation between $\mathbf{N}(z)$ and $\delta_{\rm ext}$. We leave further investigation that makes use of clustering information to future work.

\section{Conclusion}\label{S:conclusion}
We develop a spectral line de-confusion technique for LIM experiments, where different spectral lines emitted by sources at different redshifts can be observed in the same frequency channel and then confused. Unlike most of the previously proposed methods that decompose the line signals in the power spectrum space, we perform a phase-space de-confusion that reconstructs  the individual line intensity maps, if multiple spectral lines of a redshifted source population are observable. The reconstructed line intensity maps are direct data products of an LIM experiment, and can be used to trace the underlying density field for various science applications.

Our method is based on the information that multiple spectral lines emitted by redshifted sources are mapped onto distinct observed frequencies, which give deterministic features in the observed spectrum that can be fitted by a template. With a set of spectral template models and assuming the sparse approximation, we fit the LIM data iteratively with the MP algorithm.

As an example, we consider an LIM survey with similar survey parameters as the ongoing EoR [\ion{C}{ii}] experiments, TIME, and CONCERTO. The intervening CO line intensity maps at $0.5\lesssim z\lesssim1.5$ can be extracted with our technique, since multiple CO rotational transitions are observable. We demonstrate that with the assumed signal model and realistic noise level, our reconstructed CO maps reach $\sim 80\%$ spatial correlation with the true maps. In addition, in our assumed signal model and realistic survey setup, the VID of individual lines can be correctly extracted with a high S/N ratio down to the $\ell_*$-scale. The CO luminosity function derived from the VID measurement can provide information on galaxy formation and evolution as traced by the CO distribution across cosmic time. The reconstruction performance is robust against a realistic level of line ratio uncertainties and continuum foreground mitigation process.

In addition to probing  the large-scale luminosity and density fields, the reconstructed line intensity maps can also be used for a variety of applications. As a demonstration, we show that using the reconstructed map instead of the original LIM dataset can effectively reduce uncertainties in cross-correlation measurements, and improve the performance of interloper masking to reveal the high-redshift line emissions. Furthermore, given that the reconstructed intensity map, even if incomplete, traces the matter density on large scales, we construct an estimator capable of estimating the total interloper power. The estimator invokes the cross-correlation of the reconstructed map with an external density tracer such as galaxies in the linear clustering regime. This approach has the potential to fully specify the interloper and high-redshift source populations and warrants future investigation.

While we mainly discuss the application for an EoR [\ion{C}{ii}] LIM experiment in this paper, this technique is not restricted to this setup. We demonstrate that our technique can successfully extract redshifted optical line signals from a SPHEREx-like experiment in the near-infrared. The technique is a general framework that can be readily applied to mitigate line de-confusion problems in LIM experiments and enhance the science returns.

\acknowledgments
We are grateful for the helpful discussions with Chun-Lin Liu, Patrick Breysse, Matthieu B{\'e}thermin, Emmanuel Schaan, Matt Orr, the TIME collaboration, the Caltech ObsCos group, and the participants of the conference ``Lines in the Large Scale Structure.'' We would like to thank the anonymous referee for
valuable comments that improved the manuscript. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

