LIMpy: A Semianalytic Approach to Simulating Multiline Intensity Maps at Millimeter Wavelengths

Mapping of multiple lines such as the fine-structure emission from [C ii] (157.7 μm), [O iii] (52 and 88.4 μm), and rotational emission lines from CO are of particular interest for upcoming line intensity mapping (LIM) experiments at millimeter wavelengths, due to their brightness features. Several upcoming experiments aim to cover a broad range of scientific goals, from detecting signatures of the epoch of reionization to the physics of star formation and its role in galaxy evolution. In this paper, we develop a semianalytic approach to modeling line strengths as functions of the star formation rate (SFR) or infrared luminosity based on observations of local and high-z galaxies. This package, LIMpy (Line Intensity Mapping in Python), estimates the intensity and power spectra of [C ii], [O iii], and CO rotational transition lines up to the J levels (1–0) to (13–12) based both on analytic formalism and on simulations. We develop a relation among halo mass, SFR, and multiline intensities that permits us to construct a generic formula for the evolution of several line strengths up to z ∼ 10. We implement a variety of star formation models and multiline luminosity relations to estimate the astrophysical uncertainties on the intensity power spectrum of these lines. As a demonstration, we predict the signal-to-noise ratio of [C ii] detection for an EoR-Spec-like instrument on the Fred Young Submillimeter Telescope. Furthermore, the ability to use any halo catalog allows the LIMpy code to be easily integrated into existing simulation pipelines, providing a flexible tool to study intensity mapping in the context of complex galaxy formation physics.


INTRODUCTION
Observations of redshifted line emissions from atomic and molecular gas in galaxies and intergalactic medium trace the underlying dark matter density fluctuations.Several factors influence the strength of various spectral lines, including the star formation history (SFH), metallicity, and the host halo mass of galaxies.At high redshifts, z 6, it is an arduous task to resolve each individual galaxy in a survey field to understand the physics behind galaxy formation, evolution, and their connection to the intergalactic medium.Multi-line intensity mapping (MLIM) encapsulates integrated emissions from both luminous and faint sources, providing rich information about galaxy clustering, star formation rate density (SFRD), and galaxy luminosity functions (Visbal & Loeb 2010;Visbal et al. 2011;Kovetz et al. 2017;Bernal & Kovetz 2022).The detection of different atomic and molecular line intensities at a particular observational frequency probes the Universe at different epochs (or redshifts); thus, it provides a unique opportunity to construct a three-dimensional (3D) map of the Universe by the measurement of several line emission using a couple of observational frequencies.
Detecting the power spectra of fine structure lines, such as [CII] and [OIII], at high redshift (z 6) has the potential to reveal the sources of reionization and their clustering properties (Dumitru et al. 2019;Padmanabhan 2019;Padmanabhan et al. 2022;Karoumpis et al. 2022;Sun et al. 2022).Furthermore, mapping the Universe using various rotational transitions of CO (J-level transition) can probe the formation of structures at high redshifts, offering insights into the process of reionization and the star formation history of the first-generation galaxies (Kovetz et al. 2017;Breysse et al. 2022).Employing a tomographic approach to exploring the Universe enables the measurement of key quantities, including the growth factor of structures, the Hubble constant, and the equation of the state of dark energy (Kovetz et al. 2017;Karkare & Bird 2018;Bernal et al. 2019b;Silva et al. 2021).A joint analysis of all lines could prove valuable for constraining the inflationary paradigm by limiting f N L (Moradinezhad Dizgah & Keating 2019;Bernal et al. 2019a;Chen & Pullen 2022).Detecting the 21 cm−[CII] cross-power spectrum and cross-bispectrum signals can aid in mitigating the low-redshift contamination of 21 cm data, and incorporating 21 cm observations may enhance constraints on astrophysical parameters (Beane & Lidz 2018;Dumitru et al. 2019;Schaan & White 2021).
Several observational efforts have been made to detect the 21 cm line emission to study the cosmic dawn, the epoch of reionization (EoR), and late-time structure formation.Intensity mapping of other lines, such as fine-structure emission from carbon [CII] line (157.7 µm), doubly ionized oxygen [OIII] (88.4 µm), and rotational emission lines from CO are of particular interest of upcoming LIM experiments (Suginohara et al. 1998;Righi et al. 2008;Lidz et al. 2011;Carilli 2011;Fonseca et al. 2017;Gong et al. 2017;Kovetz et al. 2017;Chung et al. 2020;Padmanabhan 2018Padmanabhan , 2019;;Dumitru et al. 2019;Chung et al. 2019;Kannan et al. 2022;Murmu et al. 2021;Karoumpis et al. 2022).Despite the several advantages offered by the MLIM technique, there are key challenges to detecting a particular line emission across a broad redshift range in the presence of foreground and instrumental noise.As many of the lines emitted from other sources can be redshifted to the same observational frequency channel, it will create line confusion by adding extra emission to the particular line emission we aim to detect (Lidz & Taylor 2016;Cheng et al. 2016).These lines are called 'interlopers', and their contamination presents an obstacle to the detection of a particular line coming from the sources at a certain redshift.Moreover, the uncertainty of the star formation history (SFH) in galaxies and its relation with the mass of the host halos arises from the lack of observational data, particularly at high redshift when the reionization process occurred.However, the uncertainties in star formation and their relation with the host halos can be well explored by the different high-resolution simulations, such as UniverseMachine (Behroozi et al. 2019), IllustrisTNG (Pillepich et al. 2018), and Emerge (Moster et al. 2018).
Multiple experiments like FYST1 (Aravena et al. 2021), SPHEREx2 (Doré et al. 2018), TIME (Crites et al. 2014), CONCERTO3 (Ade et al. 2020a), COMAP4 (Cleary et al. 2021), EXCLAIM (Ade et al. 2020b), aim to cover a broad range of scientific goals, from the detection of EoR signatures to the formation of stars in galaxies.To explore the synergies among these experiments requires modelling and simulating the desired signal over a broad redshift range.In this work, we develop a package, LIMpy5 , to model and simulate several line emissions up to z 10.We implement a range of models for the star formation histories and the relations between multi-line luminosity and these star formation histories, based on analytic expressions and simulations.Collecting many models in one place allows us to explore the astrophysical uncertainties of the amplitude and shape of the signals as well as the level of contamination from interlopers.Additionally, we determine the power spectrum of line intensities both analytically through the halo model and from simulations.We adopt a map-making approach by utilizing halo catalogues generated from N-body simulations, which enables us to forecast the detectability of power spectra for future line intensity mapping (LIM) experiments.In our analysis, we incorporate the effects of beam convolution and develop realistic simulations of instrumental noise, comparing the results with those of the analytic halo models.By combining the intensity signals of various lines, LIMpy offers valuable insights into the astrophysical processes governing these emissions and their potential detectability in future experiments.Ultimately, this approach paves the way for a deeper understanding of the underlying physics and the potential impact of interlopers on the observed signals.
The LIMpy package can simulate multi-line intensity maps relatively quickly, which is useful for interpreting the signals once the observation of a particular line is made.This package also comes with several analysis techniques for calculating the three-dimensional isotropic power spectrum (P 3D (k)), the anisotropic power spectrum (P 3D k , k ⊥ ) in 3D, and the angular power spectrum in 2D (C ), so that line intensity maps can be analyzed in different ways to extract the maximum encoded information.Simulations of multi-line intensity maps across a broad redshift range are helpful in performing cross-correlations between two line intensity maps at the same redshift, as they probe the same sources and underlying dark matter density fluctuations.Furthermore, scanning the Universe at different redshifts with the MLIM technique is not only a promising probe but also carries an opportunity to perform cross-correlations with galaxy surveys and CMB secondary anisotropies, e.g., CMB weak lensing, thermal and kinetic Sunyaev-Zel'dovich (tSZ and kSZ, respectively) effects (Sato-Polito et al. 2021;Schaan & White 2021;Chung 2022).
This paper is structured as follows: in Section 2, we provide an overview of the theoretical framework for line intensities, discussing their connection to the SFR of galaxies.In Section 3, we introduce the halo-model formalism used to calculate the power spectrum of line intensities at specific redshifts and present the results obtained through this approach.In Section 4, we showcase the simulation results by describing the steps taken to generate intensity maps, and we present the detectability of the CII 158 signal in Section 5. Finally, in Section 6, we summarize our findings and conclusions.By exploring the theoretical and simulation-based aspects of line intensity mapping, we provide insights into astrophysical modelling uncertainties and the potential for future observational efforts in this area.
Throughout this study, we assume a flat ΛCDM universe with cosmological parameters as defined by the Planck TT, TE, EE+lowE+lensing results (Planck Collaboration 2018).In the rest of this paper, we denote atomic line emission by writing together the line name and its wavelength in micrometre, e.g.CII 158.For molecular line emission from CO, we denote the lines with the upper rotational transition level to the lower level, e.g., CO (1-0).We followed the same naming convention in LIMpy, and line names can be passed to calculate the necessary quantities.

THEORY OF LINE INTENSITY MAPPING
The rest-frame frequency of a particular line emission, ν rest , at redshift z em will be observed at present by an instrument with an observational frequency ν obs , such that ν obs = ν rest /(1 + z em ).An instrument can be designed to probe a bright line from a broad redshift range to understand the different physical processes at that time by selecting several frequency channels.For instance, the Epoch of Reionization Spectrometer (EoR-Spec) on FYST with the observational frequency from 220-410 GHz is set up to detect the CII 158 line emission across the broad redshift range of 7.6 to 3.6 (Aravena et al. 2021).In Figure 1, we show the redshift evolution of a few bright-line emissions that fall mainly in the FYST's EoR-Spec frequency coverage, 220 -480 GHz.All the lines that intersect the horizontal line representing a frequency channel will carry the information from the redshifts corresponding to the intersection.If EoR-Spec on FYST aims to detect the CII 158 lines from z ∼ 7.3 using ν obs = 220GHz, all other lines that cross the 220 GHz frequency line, such as all CO transitions from CO (2-1) to OIII 88,and OI 145,etc., will also come from different redshifts into the same frequency channel.In this case, one has to clean the signals of other lines to detect the desired line emission.
The detection of CII 158 and OIII 88 line emissions at z 6 will play a crucial role in understanding the epoch of reionization.In contrast, detecting higher J-ladder transition will provide us with information about structure formation and galaxy evolution during the postreionization era.To quantify their relative contributions to the total observed signal, we calculate the power spectra of the signals using both the halo model approach and N-body simulations.We show a few faint lines, such as OI 145, OI 63, OIII 52, which could act as interlopers because of their redshift overlaps with FYST's EoR-Spec frequencies.Simple modelling of these lines is also important to understand the foreground contamination of CII 158 and OIII 88.The FYST's EoR-Spec survey will scan the sky with a frequency range of 220-410 GHz at a spectral resolution of R ∼ 100.This corresponds to the redshift coverage of z CII158 ∼ 3.6−7.6for CII 158 line and z OIII88 ∼ 7.3 − 14.5.We show the redshifts for the line emission corresponding to a few central frequencies of FYST's EoR-Spec, such as 220 GHz, 280 GHz, 350 GHz, and 410 GHz.With the LIMpy code, the intensities and power spectra of any selected lines can be generated at any redshift between z ∼ 0-10.However, in this paper, we show the results only for the redshifts mentioned in Table 2.
We describe the workflow of the LIMpy package in Figure 2. The main ingredients are fed to the code as input to initialize it.In Section 2.1 we describe the in-build models of star formation histories of galaxies that can be passed to the code by mentioning the model name according to the documentation.The default cosmological parameters are based on the Planck 2018 paper.Users can modify these parameters by modifying the input file.There are several models that convert the SFR to the different line luminosities, and these models can be passed by changing the inputs.Once the basic cosmological and astrophysical parameters are initialized, then the code can calculate the power spectrum either based on the halo-model approach or simulations.Next, we calculate the power spectrum and forecast the signal-to-noise ratio for a particular experiment, providing the configurations of a telescope for the white noise calculation.In principle, different noise sources, such as atmospheric, foreground contamination, instrumental white noise, etc., can be passed in the code altogether to calculate the signalto-noise ratio.The final goal is to make forecasts for parameters based on particular observations or mock data.Either of Fisher forecasts or MCMC algorithms can be applied for parameter estimation using the LIMpy modules.
In the following subsections, we review and summarize the basic properties of star formation histories and their relationship with atomic and molecular line luminosities.These models are based on several assumptions, and changing those assumptions will typically lead to a change in results.The detailed analysis and interpretation of SFR based on galaxy formation models are out of the scope of this paper.Our main goal is to quantify the SFR as a function of halo mass M halo and z so that we can calculate the necessary quantities to estimate the power spectra of line emissions over a broad redshift range.

Empirical models of SFR
One of the most complex problems in the field of modern astrophysics is how stars form in galaxies and their role in the process of galaxy evolution.The SFR across cosmic time and its relation with halo mass are key to understanding galaxies' morphology, chemical and physical properties.Several simulation suites of galaxy formation incorporate complicated astrophysical processes in galaxies and are capable of shedding light on the SFR−M halo relation across cosmic time (e.g., Crain et al. 2015;Springel et al. 2018;Henden et al. 2018;Behroozi et al. 2019).Multi-wavelength observations of galaxies in the UV by HST and Far-IR observations by the Herschel telescope reconstructed the cosmic star formation density out to redshift z 10 (Madau & Dickinson 2014).Due to the lack of observational data at high redshift (z 4), statistical errors on the SFRD increase significantly.We aim to reconstruct the SFR empirically from several models and simulations that vary the mass of host halos across a wide range of redshifts.We incorporate five SFR models that can be used to produce line intensity maps, namely Behroozi19 from UniverseMachine (Behroozi et al. 2019), Tng300 and Tng100 (Springel et al. 2018;Pillepich et al. 2018) from IllustrisTNG, as well as fitting functions such as Silva15 (Silva et al. 2015) and Fonseca16 (Fonseca et al. 2017).
In Silva15 (Silva et al. 2015), the average SFR is extracted from the post-processed simulated galaxy catalogue (De Lucia & Blaizot 2007;Guo et al. 2011), in which the minimum halo mass is set to 10 8 M /h.The SFR − M halo scaling relation is applicable over a broad redshift range, from z = 0 to 20.This empirical relation can be approximated to z 20 for studying the high redshift line intensities, particularly for OIII 88 and OIII 52.In the Silva15 SFR model, SFR is parameterized as a function of two power-law terms of halo mass, whereas the Fonseca16 model uses the SFR based on the same simulated catalogue as the three power law exponent of halo masses.The parameters for the SFR function according to the Fonseca16 model are given for the redshift range 0-10, and we keep the SFR fixed for the redshift range 10 − 20 as the same as for the SFR at redshift z = 10.Various physical processes, such as Several built-in star formation models and multi-line luminosity models are implemented as inputs to the code.Based on these input choices, the package will calculate the power spectrum relying on either the halo model approach or painting the line luminosities on an externally provided halo catalogue.The package can make line intensity maps, and if the specification of an experiment is provided, it can calculate the signal-to-noise ratio.Furthermore, LIMpy can be used for parameter estimation based on Markov Chain Monte Carlo (MCMC) or Fisher matrix methods.These methods incorporate observational data to infer the parameters that best describe the underlying astrophysical processes.
galaxy mergers, the effect of the environment, feedback, etc., are involved in the SFR across a broad range of redshifts, which is very complex to model.Hence, we adopt fitting SFR functions such as Fonseca et al. (2017) and Silva et al. (2015).For comparison, we use the output of SFR for each halo from IllustrisTNG simulations done in box size for L = 100 cMpc and 300 cMpc (hereafter TNG100 and TNG300, respectively) (Springel et al. 2018;Pillepich et al. 2018).
We adopt the output of the UniverseMachine simulations6 to infer the SFR across the redshift z = 0 − 10 ( Behroozi et al. 2019).The empirical methods for tracking down the SFR of each halo across the redshift range are constrained by the observations such as galaxy UV luminosity functions, observed stellar mass functions, quenched fractions, etc.Furthermore, to evaluate the best-fit function from the scattered SFRs from the TNG100 and TNG300 simulations, we take 100 bins in halo mass from 10 10 M /h to 10 15 M /h and take the median value of all the SFRs that fall into each mass bin.We fix the minimum mass of the line emission sources to M min = 10 10 M /h throughout this paper.We use these star formation histories to understand the multi-line luminosities and their astrophysical uncertainties due to the scatter of SFR for a fixed halo mass.
In Figure 3, we show how the SFR varies with the masses of host dark matter halos as predicted by different empirical models.The ratio between the maximum and minimum SFR for a fixed halo mass varies at different redshifts due to the complex physical process of star formation.For M halo = 10 11 M /h, this ratio becomes 250, 96, and 11 at redshifts z ≈ 1, 4, and 6, respectively.This figure provides insights into the evolution of star formation in halos of varying masses and at different cosmic epochs, shedding light on the complex interplay between dark matter, gas, and other astrophysical processes that shape the growth and evolution of galaxies over time.The exact nature of the SFR − M halo relationship is complex and depends on several factors, including the efficiency of gas cooling, the ability of gas to collapse into small-scale structures, and feedback processes such as supernova explosions that can regulate star formation (Conroy & Wechsler 2009).Understanding the  2017).This plot captures the complex picture of star formation history in halos, as all dark matter halos with the same mass do not form the same amount of stars.The uncertainty in the SFR can propagate to the luminosity of the various emission lines.Careful consideration of the uncertainties in the SFR is necessary when interpreting intensity mapping observations and making predictions about the underlying astrophysical processes that drive the formation and evolution of galaxies over cosmic time.
underlying physical mechanisms that govern the star formation process in galaxies and their dependence on halo mass and redshift is crucial for developing a comprehensive picture of galaxy formation and evolution (Sun et al. 2022).For simplicity and optimization purposes, we use the average SFR − M halo relation to estimate the multi-line luminosities, ignoring the dependencies of other astrophysical parameters related to the complex star formation history in haloes.

SFR -L line relation
A crucial question that arises in the development of Line Intensity Mapping (LIM) models is the identification of the key factors that trace the observed multi-line luminosities in galaxies.It is assumed that multi-line luminosities trace the star formation histories, and SFR can be converted to line luminosities using a power-law relation.In the previous subsection, we modelled how the SFR depends on the mass of halos, and we discussed how the multi-line luminosities are related to the SFR so that for a given halo mass, we can estimate the multiline luminosities.We incorporated several L line − M halo relations in LIMpy to study the modelling uncertainties in the intensity maps.
The luminosity of these lines scales with the SFR as L line = R line × SFR for Visbal10 model, where R line is the conversion factor (Visbal & Loeb 2010), which does not evolve with the redshifts.Assuming all galaxies have the same R line , their values are given by 6 × 10 6 and 2.3 × 10 6 in L /(M /yr) unit for CII 158 and OIII 88 lines, respectively.The values of R line for J-ladder transitions of CO molecules are given in the Visbal & Loeb 2010, see Table 1.In the Silva15-m1 (and m2, m3, and m4) model, the luminosity of CII 158 is modelled as a power law relation; log L CII 158 = α + β log(SF R).The four sets of models are given by the different values of α and β that we specify in LIMpy by the name Silva15-m1, Silva15-m2, Silva15-m3, and Silva15-m4 (Silva et al. 2015).For the Fonseca16 model, the luminosity of CII 158, OIII 88, OI 145, OI 63, and OIII 52 lines can be expressed as the same power law relation with SFR, but the coefficients are different than Silva15 model as mentioned in Fonseca et al. (2017).For the Silva15 and Fonseca16 models, the coefficients α and β do not change with the redshift, but the multi-line luminosity of lines varies with the halo mass only because of the evolution of SFR with redshift.The redshift evolution of the coefficients is captured with a modified version of the power in the Lagache18 model (Lagache et al. 2018).In addition, we also model the scaling relation from low redshift observations of the molecular line from the ALMA-Alpine experiment that can be written in a We use several scaling relations to model CO line emissions in our study.One such relation is based on the spectra obtained from the Herschel SPIRE Fourier Transform Spectrometer (Kamenetzky et al. 2016).The luminosity of CO molecular emission, L CO , is found to depend on the FIR luminosity of galaxy samples L FIR .In order to calculate the luminosity of all lines for a given SFR, we convert it into L FIR , where L FIR = 1.1 × 10 10 × SF R (Carilli 2011).Additionally, we incorporate another model for CO molecular transitions based on ALMA observations (Greve et al. 2014).This model allows us to estimate the luminosity of full rotational transitions of CO molecules using data from Herschel SPIRE-FTS and ground-based telescopes Greve et al. 2014, see Table 3.The full J-ladder rotational transitions of CO are essential for understanding the relative contributions in the interlopers which contaminate the desired signal that we aim to detect.Using multiple models allows us to estimate the CO line emissions, which helps us to account for the uncertainties associated with these relations.
Figure 4 shows the evolution of different line luminosities for the Silva15 star formation model.The plot shows the ratio between the maximum and minimum luminosities of the CII 158 lines, the CO (7-6) line, and the OIII 88 line as a function of redshift for a fixed minimum halo mass of M min = 10 10 M /h.At a redshift of z ∼ 3.8, the ratio between the maximum and minimum CII 158 line luminosities is approximately 45.However, for the same minimum halo mass, the ratio for CO (7-6) luminosity becomes 110 at z ∼ 2.66 and 310 at z ∼ 0.96, highlighting the increasing spread in luminosity ratios as the redshift decreases.The ratio for the CO (7-6) line decreases to 1.5 and 6, respectively, if the minimum halo mass is set to 10 11 M /h, suggesting that the choice of minimum halo mass can significantly impact the luminosity ratios.Finally, for the OIII 88 line luminosity with M min = 10 10 M /h, the ratio between the maximum and minimum luminosity is 6 and 8 at redshifts z= 14.5 and 7.3, respectively, indicating that this line is less sensitive to changes in redshift than the other lines considered in the plot.

ANALYTIC MODEL
In this section, we employ a halo model formalism to compute the power spectrum of multi-line intensities.The intensity of lines emitted at z em can be expressed .The power spectra of the CII 158 line at redshifts 7.6, 5.8, 4.4, and 3.6.The dashed lines in each panel show the contribution to the signal from the clustering term alone, while the solid lines represent the total signal, including both the clustering and shot noise terms.The solid red lines in each panel are the mean line of all the models, representing the average signal predicted by the different theoretical models that we consider here.The shaded grey region represents the scales where the dominant term is shot noise, while the yellow-shaded area roughly corresponds to the scales where the clustering term is larger than the shot noise term.The shape and amplitude of the CII 158 have important implications for studying ISM physics and structure formation at high redshift, as they provide insights into the large-scale structure of matter in the Universe and the properties of the sources that generate the CII 158 signal. as (1) In this equation, c represents the speed of light in a vacuum, and H(z em ) denotes the Hubble parameter at the redshift of line emission.The halo mass function is represented by dn/dM .Throughout our study, we utilize the Tinker halo mass function for our calculations (Tinker et al. 2008).Here, M min refers to the minimum mass of the halos contributing to the intensity maps, while M max signifies the upper mass threshold of the sources.
The power spectrum of fluctuations of the line intensity is the summation of 1-halo and 2-halo terms that can be written as (2) In the above equation, P m (k, z) is the matter power spectrum and P shot is the shot noise term.We calculate the matter power spectrum using CAMB (Lewis & Challinor 2011) under the linear approximation.The bias of the line emission, b line , is proportional to the bias of line emitting sources, which can be written as Here, b h is the bias of dark matter halos.We use Colossus7 package to calculate the bias of dark matter halos and halo mass function (Diemer 2018).The bias of lines accounts for the clustering properties of the power spectrum, which has significant contributions at large scales.
Finally, the shot noise term of the power spectrum is proportional to the line luminosity function, which is given by The shot noise term has the same contributions at all scales.
Figure 5 presents the power spectrum of CII 158 for various models, highlighting the clustering and total signal comprising both the clustering and shot noise terms.Despite using the same SFR as input, the L line − SF R relation, the dispersion on the amplitude of power spectra at these redshifts due to the modelling differences exceed one order of magnitude at these redshifts.Some models are based on the L line − M halo relation, while others convert SFR to CII 158 line luminosities.For the latter case, we first calculate the SFR for different halo masses and then determine the line luminosities based on the L line − SF R relation.Similarly, we can generate the power spectra of OIII 88 and molecular lines from CO (1-0) to CO (13-12) for available models.For an example case, we show the power spectra of OIII 88 and CO (7-6) at redshifts corresponding to FYST's EoR-Spec in Appendix B. In this way, we can assess the contribution of these lines to the total signal at a particular frequency channel and determine their detectability in the presence of interlopers.
In Equations ( 3) and ( 4), we made the assumption that there is only one star-forming source in each halo.However, this is not the case for high-mass halos.To account for multiple star-forming sources in halos, we employ the halo occupation distribution (HOD) model.The mean occupation functions for central and satellite galaxies in a halo of mass M h are given by (Zheng et al. 2005): In the above equations, N cen (M h ) and N sat (M h ) represent the average number density of central and satellite galaxies, respectively.M th denotes the threshold halo mass required to host a central galaxy, while M cut represents the minimum mass necessary for hosting satellite galaxies.σ logM is the width of the transition in the step-like error function, α g refers to the power law exponent, and M 1 is the mass normalization factor.The HOD-model parameters are given as log M th = 10 8 , σ logM = 0.287, log M cut = 12.95, log M 1 = 13.62, and α = 0.98 (Zheng et al. 2005).By incorporating the HOD model, we can more accurately account for the distribution of star-forming sources in halos, particularly in high mass halos.The inclusion of the HOD model provides a more comprehensive picture of the power spectrum of multi-line intensities.Figure 6 presents the percentage difference in the power spectra of CII 158 lines resulting from the inclusion of the HOD model in our calculations.The HOD model accounts for the line emissions both from central and satellite galaxies, and its inclusion can significantly affect the 1-halo term of power spectra.The plot shows the increase in the power spectra due to the HOD model for different redshifts and scales, represented by the percentage difference compared to the power spectra without HOD.At a scale of k ∼ 5 h Mpc −1 , the power spectra of CII 158 lines increase by 73%, 25%, 2%, and 0.1% at redshifts 3.6, 4.4, 5.8, and 7.6, respectively, highlighting the significant impact of the HOD model on the power spectra of CII 158 lines.

SIMULATED MAPS OF MLIM
In addition to theoretical modelling, LIMpy also performs the simulation of multi-line intensity maps at several redshifts.We utilize the semi-numerical cosmological simulation, 21cmFAST8 , to generate the dark matter halo catalogue.We execute the simulation on a box with a length, L= 800 cMpc (≈ 544 cMpc/h) (Mesinger et al. 2011).The initial conditions for generating the perturbations in density are set at z = 300.The density field evolves over cosmic time following linear perturbation theory.Next, we generate snapshots of several redshifts for different lines corresponding to the FYST's EoR-Spec frequency channels.We set the minimum mass of the halos to M min = 10 10 M /h.The simulation setup uses a number of grids along the box length, N grid = 1024, meaning that the total number of dark matter particles is N 3 grid .This enables us to step down to a length of 0.53 cMpc/h.After obtaining the halo field or resolved catalogues from simulations, we assign a specific line luminosity to those halos based on an SFR model and line luminosity model.
The primary advantage of line intensity mapping simulations is that we do not need to resolve individual sources.Consequently, we use a low-resolution simulation with a smaller N grid .This approach reduces the simulation time.We save all the intensity grids at different redshifts to calculate the power spectrum.
The 3D line intensity power spectrum in a simulation box is expressed by: Here, V box represents the total volume of the simulation box, and Ĩ is the Fourier transform of the intensity grid.We then perform the Fourier transform of the intensity grid using the NumPy FFT module.The intensity of each cell is calculated as (Dumitru et al. 2019): In principle, CO transitions with J up ≥ 4 at low redshifts will act as interlopers for all four of FYST's EoR-Spec frequency channels.However, we only show the CO (7-6) transition as an example case.We project the intensity grid of the length of 35 Mpc, 16 Mpc, 60 Mpc for CII 158, CO (7-6), and OIII 88 lines, which roughly correspond to the frequency resolution of the EoR-Spec on FYST at the central frequency, 280 GHz.
Once we generate the intensity grid with a reasonable value of N grid to achieve the M min for line emitting sources, we need to incorporate the effect of frequency resolution (δν obs ) as an experiment cannot resolve the sources along the redshift axis for the smaller value than the δν obs .Therefore, the effective number of grids along the redshift axis will depend on the frequency resolution of an experiment.With a resolution of N grid = 1024, we select halos above the mass M min 10 10 M /h.The 21cmFAST code does not explicitly resolve each halo in the simulation box, but it generates a halo field seminumerically that can be accurately compared with the output from N-body simulations (Mesinger et al. 2011;Mas-Ribas et al. 2022).This process saves memory usage and makes it run faster.If δν obs is 2.8 GHz around the central observational frequency, ν obs = 220 GHz for an experiment to probe CII 158 line emission, the corresponding redshift resolution is δz ≈ 0.07.This δz corresponds to the box length along the redshift axis .We display simulated maps of CII 158 (first row), OIII 88 (second row), CO (7-6) (third row) line intensities at redshifts corresponding to the central frequencies of the FYST's EoR-Spec experiment.The simulation boxes are generated from 21cmFAST package for a 544 cMpc/h box, and we keep fixed the same initial condition for all the simulations at different redshifts.The columns correspond to the specific observational frequency, as indicated by the column titles.The fourth row of the plot compares the dimensionless power spectra of these lines based on the Tng300 star formation model and the Visbal10 line luminosity model (Visbal & Loeb 2010).For visualization purposes, Gaussian beam convolutions were applied to the maps with full-width half-maximum (FWHM) beam sizes of 58, 45, 37, and 32 arc seconds from left to right.
(we define it to be along the z axis of the cartesian coordinate system) δL z ≈ 44 cMpc/h and δN grid,z ≈ 90.Therefore, the experiment with this configuration cannot resolve the sources along the redshift axis that fall between z = 5.80 and 5.87.In this case, the total number of cells of the same simulation box becomes 1024 × 1024 × 11 as the grid points along the z-axis reduces to N grid,z = 1024/δN grid,z ≈ 11.We take the average intensity of all the intensity grids within this frequency resolution.We note if there is no mention of ν obs or the length corresponding to ν obs exceeds the length of the simulation box along the z-direction, the code does not apply the effect of the frequency resolution and calculates the power spectrum based on the grid points that were used to generate the halo catalogue.
The code presented in this paper is versatile and can accommodate any type of halo catalogue that is provided as input.Our code requires the halo catalogue to contain two essential pieces of information: the halo mass in units of M /h, and the halo positions (x, y, z) in Cartesian coordinates specified in units of Mpc/h.By accepting any halo catalogue as an input, the code offers the flexibility to utilize halo catalogues generated from full N-body simulations and perform detailed astrophysical analyses.This functionality is especially useful in studying the properties and evolution of halos in largescale structure simulations, as well as in exploring the connection between halo properties and other astrophysical observables.
Next, we apply beam convolution techniques to recreate the actual observation.Although the actual beam of an experiment can have a complex pattern, for simplicity we assume the beam can be approximated to be Gaussian.Then the beam pattern is characterized by θ FWHM , the beam size in arcminute unit at the full width at half maximum (FWHM), and the standard deviation of the beam is σ beam = θ FWHM / (8 log 2).We used the Astropy 9 package to perform the beam convolution on the simulated line intensity maps.The beam convolution does not change the power spectrum at large scales but reduces the power significantly at small scales, above the scales that correspond to the beam size.
The shape and amplitude of power spectra of different lines at various redshifts provide valuable information about the properties of the intergalactic medium and galaxy populations.In Figure 7, we present the simulated intensity maps of CII 158, CO (7-6), and OIII 88 line emissions from the halos at several redshifts corresponding to FYST's EoR-Spec central observational 9 https://www.astropy.org/frequencies.We project the intensity grid of the length of ≈ 1.3 cMpc/h for CII 158, CO (7-6), and OIII 88 lines.We show the three-dimensional power spectra of intensity maps without performing the beam convolution to show both clustering and shot noise terms.However, for the visualization, we convolve the intensity maps with the Gaussian beam, and the FWHM values are varied according to the EoR-Spec on FYST ν obs (Aravena et al. 2021).For all intensity maps, we consider M min = 10 10 M /h the redshifts mentioned in Figure 6, except the OIII 88 intensity map at z ∼ 14.5.Since there are no high-mass halos 10 10 M /h present at such a high redshift, we show the halos whose masses are larger than 10 9 M /h for that case.
At ν obs ∼ 220 GHz, the power spectrum of CO (7-6) lines is ∼ 350 times larger than CII 158 power spectrum at k ∼ 0.1 h/Mpc, but the ratio becomes 40 to 60 times higher at the shot noise dominated scales, k 1 h/Mpc.At this frequency corresponding to z ∼ 14.5, the OIII 88 signal is negligible as there are very few lineemitting sources.This comparison shows it is impossible to detect OIII 88 from such high redshift by the ongoing and planned MLIM experiments.However, at ν obs ∼ 410 GHz, the CII 158 signal becomes larger than CO (7-6) by a factor of ∼ 6 at k = 0.1 h/cMpc and ∼ 2 at k = 1 h/Mpc.Furthermore, for the same redshift at z ∼ 7.4, the OIII 88 power spectrum is approximately 4.5 and 1.7 times larger than the CII 158 power spectrum at k ∼ 0.1 h/Mpc and 1 h/Mpc, respectively.Therefore, by using two frequency bands, 220 and 280 GHz, we could detect OIII 88 and CII lines and perform crosscorrelation studies as they come from the same sources.

DETECTABILITY
In this section, we forecast the detectability of CII 158 and OIII 88 by considering the specifications of an EoR-Spec experiment.The signal-to-noise ratio is proportional to the number of observed modes present in the survey volume.To determine the number of modes between the wave number k and k+∆k, we use the following equation: where k i is the central wave number in the bin width, ∆k i .The survey volume of an experiment is given by (Gong et al. 2017;Dumitru et al. 2019)  In the above equation, λ line symbolizes the rest frame wavelength of the line emission, S A is the effective survey area of an experiment, and B ν is the frequency bandwidth.For an EoR-Spec-like experiment on FYST, we consider B ν = 40 GHz for all frequency channels.
To calculate the covariance matrix for the detection of P line (k), we employ the following formula: Here, P N represents the noise power spectrum.The source of noise can be a combination of white noise, atmospheric noise, and interloper contribution.However, in this paper, we only take the white noise contribution into account to forecast the signal-to-noise ratio (SNR) for the EoR-Spec experiment on FYST (Aravena et al. 2021).
The signal-to-noise ratio, (S/N ), is given by: By employing these equations, we can effectively forecast the detectability of the CII 158 and OIII 88 spectral lines, taking into account the white noise as the sole source of noise in our analysis.This allows us to estimate the signal-to-noise ratio for the EoR-Spec experiment on FYST and assess the overall feasibility of detecting these spectral lines.Figure 8 presents the detectability status of the CII 158 power spectrum at different redshifts corresponding to the frequency coverage of the EoR-Spec experiment on FYST.The plot compares the power spectra of CII 158 between the halo model and simulation, demonstrating the potential of EoR-Spec to detect the CII 158 signal at different redshifts.The halo mass function used to generate the halo catalogue is the Sheth-Torman mass function, while we use the Tinker mass function for the Halo model approach.The figure also displays the error bars for CII 158 lines forecasted based on the halo model approach for the different frequency bands of the EoR-Spec.The results of our analysis indicate that the EoR-Spec-like experiment will be able to detect the CII 158 signal at more than 350 σ for ten bins in the range of k min ∼ 0.1 cMpc/h to k max ∼ 10 cMpc/h at z ∼ 5.8.At higher redshifts, the signal-to-noise ratio is expected to be lower, with values of 26, 373, and 295 at z ∼ 7.4, 4.6, and 3.4, respectively.
Table 5 summarizes the signal-to-noise ratio (SNR) for the detection of CII,158 at various redshifts corresponding to the FYST's EoR-Spec frequency coverage.Given the considerable uncertainty in the amplitude of the power spectrum, we present SNR forecasts for three different scenarios: optimistic, moderate, and pessimistic, representing the largest, median, and weakest expected signals, respectively.The results demonstrate that an EoR-Spec FYST-like experiment has the potential to detect the CII 158 signal with high significance, offering a valuable opportunity to constrain theoretical models of galaxy formation and evolution.By probing the complex interplay between various astrophysical processes, these observations could provide crucial insights into the underlying physical mechanisms driving the growth and evolution of galaxies.

CONCLUSION
In the study of galaxy formation and line intensity mapping, uncertainties in astrophysical modelling of line intensity signals and variations in star formation histories within dark matter halos are crucial topics.In our research, we have developed a package that brings together various models to enable comparison and facilitate the elimination of models when making MLIM observations.The LIMpy package is a semi-numerical code that allows for the modelling of CII 158, OIII 88, and different CO J-ladder transitions in a single framework.We have implemented several star formation models, including those inferred from analytic prescriptions and state-of-the-art simulations such as IllustrisTNG and UniverseMachine, as well as empirical relations from the abundance matching approach.We assume that the SFR serves as a proxy for line luminosities and have included various SFR−L line scaling relations based on several best-fit models.Our primary objective in this study is to provide a tool for investigating the astrophysical and cosmological information derived from line intensity mapping while taking into account the modelling uncertainties inherent in such an approach.By integrating various models and scaling relations in a single framework, our approach can help interpret the observed MLIM signal.
The LIMpy package not only allows for the modelling of line intensity maps but also enables the exploration of cosmological parameters and their effects on these maps.Halo catalogues generated from various simulations, such as N-body or cosmological hydro simulations, can be inputted into the code to investigate the impact of these parameters on the line emissions.The code efficiently paints halos with line emissions based on the SFR and line luminosity models, making it an ideal tool for performing MCMC analyses on simulations to be constrained by MLIM observations.Furthermore, the simulated maps produced by LIMpy can be used to determine optimal statistics for analyzing observed MLIM data.For instance, future studies may employ voxel intensity distribution (VID) to analyze the simulated multi-line intensity maps outputted by the package (Ihle et al. 2019;Breysse 2022;Sato-Polito et al. 2021).Overall, LIMpy provides a versatile and powerful tool for studying both astrophysical and cosmological parameters at the map level and can facilitate current and future MLIM observations.
The CII 158 line emission, which exhibits a large scatter, raises critical questions about interpreting the MLIM signal at these redshifts during observations.The uncertainties in the astrophysical parameters pose a challenge in accurately constraining the properties of the intergalactic medium and galaxy populations responsible for the observed signals.Thus, low-noise MLIM observations are crucial to obtaining a robust analysis of the data and constraining the astrophysical parameters.This highlights the need for improving the modelling of the CII 158 power spectrum and other line intensities to maximize the scientific returns of MLIM experiments.
The high signal-to-noise ratio attainable through an EoR-Spec-like experiment for the CII 158 signal makes it an ideal tool for probing the parameters associated with the reionization process, such as the ionized bubble size and mean free path of photons.These observations could be used to reconstruct the luminosity function of galaxies and provide insight into the ionizing sources.Moreover, the frequency overlap of the EoR-Spec allows for the inter-line cross-correlation of CII 158 (220 GHz) with OIII 88 (410 GHz) to obtain a snapshot of the Universe at z ∼ 7.4.However, this analysis is subject to potential issues such as interloper contamination and the effect of beam convolution, which we neglected in our forecasts.Consequently, while our forecasts for detecting CII 158 are optimistic, further work is required to account for these factors and obtain a more accurate estimation of the detectability of the CII 158 signal.
However, the rotational emissions from CO molecules at low redshifts dominate over the CII 158 emissions from high redshifts, presenting a significant obstacle to studying the reionization epoch using CII 158.Additionally, extragalactic foregrounds, such as broadband emission from the cosmic infrared background (CIB), contribute to the contamination.To minimize contamination, linear combinations of maps reconstructed using different channels can be used.This approach can significantly reduce the bias for CII 158 detections.However, to obtain more realistic predictions, future work will explore bias due to interlopers, instrumental and atmospheric noise, and will build estimators to estimate the signal in their presence.

Figure 1 .
Figure 1.Redshift evolution of different atomic and molecular lines of interest in the redshift-frequency space.The four horizontal lines show the central frequencies of EoR-Spec on FYST, νcen, and the corresponding shaded regions represent the frequency bandwidths, ∆ν.Dashed lines show the rotational transitions of CO from the Jup to the (Jup − 1) level.

Figure 2 .
Figure2.The schematic flowchart of the LIMpy package.Several built-in star formation models and multi-line luminosity models are implemented as inputs to the code.Based on these input choices, the package will calculate the power spectrum relying on either the halo model approach or painting the line luminosities on an externally provided halo catalogue.The package can make line intensity maps, and if the specification of an experiment is provided, it can calculate the signal-to-noise ratio.Furthermore, LIMpy can be used for parameter estimation based on Markov Chain Monte Carlo (MCMC) or Fisher matrix methods.These methods incorporate observational data to infer the parameters that best describe the underlying astrophysical processes.

Figure 3 .
Figure3.We illustrate the assumed SFR models as a function of halo mass, three different redshifts: z ∼ 1 (left), z ∼ 4 (middle), and z ∼ 6 (right).The scatter points represent the star formation histories of individual halos in the TNG300 simulations(Springel et al. 2018), while the pink and purple solid lines depict the best-fit curves based on the TNG100 and TNG300 simulations.For comparison, we show the interpolated SFR fromBehroozi et al. (2019), and analytic models of SFR taken fromSilva et al. (2015) andFonseca et al. (2017).This plot captures the complex picture of star formation history in halos, as all dark matter halos with the same mass do not form the same amount of stars.The uncertainty in the SFR can propagate to the luminosity of the various emission lines.Careful consideration of the uncertainties in the SFR is necessary when interpreting intensity mapping observations and making predictions about the underlying astrophysical processes that drive the formation and evolution of galaxies over cosmic time.

Figure 4 .
Figure 4. Redshift evolution of the CII 158, CO76, and OIII 88 luminosities based on the models mentioned in the legends.Solid and dashed lines are the line luminosities for the halo mass 10 10 M /h and 10 11 M /h, respectively.The vertical-shaded regions are the redshift coverage of these lines corresponding to the frequency bandwidths of the EoR-Spec on FYST.similar form.The mean values of α and β are given in Table 2 of Schaerer et al. (2020).For modelling OIII 88 lines, we include the L OIII88 − SF R scaling relation defined in the code by Visbal10 (Visbal & Loeb 2010), Delooze14 (De Looze et al. 2014), Fonseca16 (Fonseca et al. 2017), Gong17 (Gong et al. 2017), Harikane20 (Harikane et al. 2020), and Kannan21 (Kannan et al. 2022).The SFR in the far-infrared is modelled in terms of the CII 158, OI 63, and OIII 88 line emissions from the Herschel Dwarf Galaxy Survey, and the scaling relations are obtained from De Looze et al. 2014, see Table 2.They find that OI 63 and OIII 88 trace the SFR better than the CII 158 lines and the disperse in the relation between SFR and L line for OIII 88 and OI 63 is improved by a factor of ∼ 2 compared with the CII 158 lines.In addition to that, we adopt the scaling relations between SFR and luminosity of OIII 88 lines based on the observations by ALMA at z ∼ 6 − 9 (Harikane et al. 2020).They find the ratio L OIII 88 /L CII 158 can be 10 times higher than this ratio at z ∼ 0, suggesting a strong redshift evolution of line luminosities.While we implement this scaling relation in LIMpy, the OIII 88 line luminosities can change across redshift due to the change of SFR.Furthermore, we use the scaling relation of OIII 88 that is derived from the observed luminosity function and SFRD at z 5 Gong et al. 2017, see Section 2.We use several scaling relations to model CO line emissions in our study.One such relation is based on the spectra obtained from the Herschel SPIRE Fourier Transform Spectrometer(Kamenetzky et al. 2016).The luminosity of CO molecular emission, L CO , is found to depend on the FIR luminosity of galaxy samples L FIR .In order to calculate the luminosity of all lines for a given SFR, we convert it into L FIR , where L FIR = 1.1 × 10 10 × SF R(Carilli 2011).Additionally, we incorporate another

Figure 6 .
Figure 6.Percentage difference in the CII 158 power spectrum due to the use of a HOD model.For this example, we calculate the power spectrum based on the Silva15 star formation model and Fonseca16 line luminosity model of CII 158.It shows that the inclusion of the HOD model is more important for the power spectrum of CII 158 at low redshift.
Figure7.We display simulated maps of CII 158 (first row), OIII 88 (second row), CO (7-6) (third row) line intensities at redshifts corresponding to the central frequencies of the FYST's EoR-Spec experiment.The simulation boxes are generated from 21cmFAST package for a 544 cMpc/h box, and we keep fixed the same initial condition for all the simulations at different redshifts.The columns correspond to the specific observational frequency, as indicated by the column titles.The fourth row of the plot compares the dimensionless power spectra of these lines based on the Tng300 star formation model and the Visbal10 line luminosity model(Visbal & Loeb 2010).For visualization purposes, Gaussian beam convolutions were applied to the maps with full-width half-maximum (FWHM) beam sizes of 58, 45, 37, and 32 arc seconds from left to right.

Figure 8 .
Figure 8. compare the power spectra of CII 158 lines at four different redshifts between the halo model and simulations.The figure also shows the error bars for CII 158 lines forecasted based on the halo model approach for the different frequency bands of the EoR-Spec on the FYST experiment.
to thank Anthony Challinor, Steve Choi, Andrea Lapi, Dominik Reichers, and Gordon Stacey for helpful discussions.AR is partially supported by the CCAT-prime collaboration.DV acknowledges the REU program through the CCAPS at Cornell University under the NSF award NST/AST-1950324.NB acknowledges support from NSF grant AST-1910021 and NASA grants 21-ADAP21-0114 and 21-ATP21-0129.AvE acknowledges support from NASA grants 22-ADAP22-0149 and 22-ADAP22-0150.

Table 2 .
We quote the forecasted cumulative signal-to-noise ratio of CII 158 lines at several redshifts for an FYST-like MLIM experiment with taking into account the effect of foregrounds.