SIMPLE: Simple Intensity Map Producer for Line Emission

We present the Simple Intensity Map Producer for Line Emission (SIMPLE), a public code to quickly simulate mock line-intensity maps, and an analytical framework to model intensity maps including observational effects. SIMPLE can be applied to any spectral line sourced by galaxies. The SIMPLE code is based on lognormal mock catalogs of galaxies including positions and velocities and assigns luminosities following the luminosity function. After applying a selection function to distinguish between detected and undetected galaxies, the code generates an intensity map, which can be modified with anisotropic smoothing, noise, a mask, and sky subtraction, and calculates the power spectrum multipoles. We show that the intensity auto-power spectrum and the galaxy-intensity cross-power spectrum agree well with the analytical estimates in real space. We derive and show that the sky subtraction suppresses the intensity auto-power spectrum and the cross-power spectrum on scales larger than the size of an individual observation. As an example application, we make forecasts for the sensitivity of an intensity mapping experiment similar to the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX) to the cross-power spectrum of Ly$\alpha$-emitting galaxies and the Ly$\alpha$ intensity. We predict that HETDEX will measure the galaxy-intensity cross-power spectrum with a high signal-to-noise ratio on scales of $0.04\, h\,\mathrm{Mpc}^{-1}


INTRODUCTION
Line-intensity mapping (LIM) is a promising survey strategy for constraining cosmological parameters and studying astrophysics of galaxies and intergalactic gas (e.g., Kovetz et al. 2017;Bernal & Kovetz 2022).Instead of relying on detections of individual galaxies to trace the large-scale structure of matter, LIM measures the integrated line emission from all galaxies and the intergalactic medium (IGM) as a biased tracer of the matter distribution.It collects light from all emitters, including those that are too faint for individual detection at high redshift, and requires lower resolution and shorter integration times than traditional galaxy surveys.
Several LIM surveys are in operation (Santos et al. 2016;Keating et al. 2016Keating et al. , 2020;;DeBoer et al. 2017; maja@mpa-garching.mpg.deCONCERTO Collaboration et al. 2020;Gebhardt et al. 2021;Cleary et al. 2022) or in preparation (Doré et al. 2014;Vieira et al. 2020;Sun et al. 2021;Switzer et al. 2021;CCAT-Prime Collaboration et al. 2022).They target lines ranging from the 21 cm line in radio frequencies to the ultraviolet Lyman-α (Lyα) line, which trace atomic or molecular gas.Although the 21 cm line is emitted by neutral atomic hydrogen in the IGM and in neutral pockets within the interstellar medium (ISM) of galaxies, other lines targeted by LIM experiments are associated with star formation in the ISM (e.g., Bernal & Kovetz 2022).Neglecting scattering and diffuse emission of the Lyα line, which can illuminate the circumgalactic medium (CGM) and IGM (see, e.g., Byrohl et al. 2021;Byrohl & Nelson 2023), these line intensities are therefore sourced within galaxies and influenced by their astrophysical properties such as their star-formation rate.
Mock intensity maps are necessary to model the signal beyond the capabilities of analytical models, estimate statistical covariance, and explore observational effects such as foregrounds, sky subtraction, and survey footprints on summary statistics (e.g., Cunnington et al. 2023).Modeling techniques of line-intensity maps have to compromise between astrophysical complexity, volume, and computational feasibility.On the one hand, cosmological hydrodynamical simulations can be postprocessed using the astrophysical properties of galaxies and the gas to infer the line emission (e.g., Moriwaki et al. 2019;Silva et al. 2021;Kannan et al. 2022;Byrohl & Nelson 2023;Liang et al. 2023).Halo catalogs from N-body simulations can also be combined with galaxy evolution models to predict line luminosities and produce intensity maps (e.g., Lidz et al. 2011;Gong et al. 2012;Li et al. 2016;Chung et al. 2019;Spina et al. 2021;Béthermin et al. 2022;Moradinezhad Dizgah et al. 2022;Sato-Polito et al. 2022).Although these approaches provide accurate small-scale clustering and include astrophysical dependences of the line intensity, it is not feasible to produce enough realizations for covariance estimation or parameter inference.To speed up the calculation, one can model the underlying dark matter density field through Lagrangian perturbation theory or masspeak patch and apply various post-processing steps to model the astrophysical dependence of line intensities (Mesinger et al. 2011;Silva et al. 2013Silva et al. , 2015;;Mesinger et al. 2016;Heneka et al. 2017;Heneka & Mesinger 2020;Mas-Ribas et al. 2023;Chung et al. 2022;Roy et al. 2023).The line intensity can also be modeled by multiplying the total matter density in a fast lognormal simulation by a bias (Alonso et al. 2014;Rubiola et al. 2022).These approaches are fast and enable simulating many realizations of large volumes, but do not account for the shot noise contribution to the power spectrum from the discreteness of the line-emitting sources.Finally, Obuljen et al. (2022) use a field-level forward-modeling approach to simulate intensity maps based on effective field theory.
To include shot noise, one can generate a galaxy catalog from a lognormal galaxy number density field via Poisson sampling.Lognormal simulations take advantage of the roughly lognormal probability density function (PDF) of matter and galaxy density distributions, measured both in N-body simulations (e.g., Kayo et al. 2001;Shin et al. 2017) and in galaxy surveys (e.g., Clerkin et al. 2017).The lognormality of the PDF of the density contrast δ implies that the logarithmic transformation field, ln (1 + δ), is a Gaussian random field whose statistics are defined entirely by its two-point correlation function or power spectrum.Agrawal et al. (2017) present a public lognormal mock generator for galaxy catalogs, including self-consistent velocities of the galaxies.The velocities enable redshift-space distortion (RSD) modeling beyond the linear Kaiser model (Kaiser 1987), which is critical for any large-scale structure measurement in redshift space such as LIM.
The lognormal galaxy catalog generator of Agrawal et al. (2017) has been extended to generate weak lensing fields (Makiya et al. 2021).In this paper, we extend it to quickly and self-consistently generate intensity maps and galaxy catalogs, which we call the Simple Intensity Map Producer for Line Emission (Simple).1Given a luminosity function for any emission line, we assign line luminosities to galaxies in a lognormal galaxy simulation, apply a selection function to obtain a galaxy catalog, and calculate the intensity on a grid.We can add noise, smooth the intensity map, model the sky subtraction, and apply a mask before calculating the galaxy and LIM auto and cross-power spectra.The simplicity of this approach enables us to quickly generate many mock intensity maps, e.g., to estimate the covariance matrices of the LIM power spectra.
As an example of its capabilities, we use the Simple framework to forecast the sensitivity of the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX; Gebhardt et al. 2021) to the cross-correlation of Lyαemitting galaxies (LAEs) and the Lyα intensity.An analytical forecast for the intensity auto-power spectrum in HETDEX was presented in Fonseca et al. (2017).
This paper is structured as follows.Section 2 derives the power spectrum formalism.In Section 3 we describe the steps of the Simple code to generate line-intensity mocks.Section 4 explains an example mock setup for HETDEX that is used in the rest of the paper.We validate the Simple code by comparing the power spectrum multipoles with theoretical predictions in Section 5. Section 6 presents the forecast for the LAE-Lyα intensity cross-correlation of HETDEX.Section 7 discusses the limitations of the Simple framework.We conclude in Section 8.
We use the following Fourier convention: where the tilde denotes quantities in Fourier space.We refer to real space in contrast to redshift space, and configuration space in contrast to Fourier space.Throughout this paper, we assume a flat Λ cold dark matter (ΛCDM) cosmology with H 0 = 67.66km s −1 Mpc −1 , Ω b,0 h 2 = 0.022, Ω m,0 h 2 = 0.142, ln 10 10 A s = 3.094, and n s = 0.9645.

POWER SPECTRUM MODELING
Consider a fluctuation δA(x) = A(x)− A(x) of a field A, such as the intensity A(x) = I(x) and the normalized galaxy number density A(x) = n(x)/ n(x) = 1 + δ g (x).Here •(x) denotes the mean field at location x over many realizations, for example the mean intensity as a function of redshift or the galaxy selection function as a function of position.The dimensionless correlation function of fields A and B in configuration space is defined as (2) The corresponding power spectrum with the dimension of volume is the Fourier transform of the correlation function, 2.1.Galaxy and intensity auto-and cross-power spectra Following the standard approach of Peebles (1980), we model the galaxy number N (x) as a Poisson point process in infinitesimal volume elements δV , so that the occupation number in each cell is N i ∈ {0, 1}.The expectation value of N in one cell is n [1 + δ g (x)] δV , where n is the mean number density of galaxies in the entire volume and δ g is the galaxy overdensity.Let each galaxy have a line luminosity L i , which is sampled from a luminosity function dn/dL.We require that the integral dL dn dL =: n converge so that φ(L) := dn dL n−1 is a PDF, for example, by setting a minimum luminosity.
The specific intensity in a cell with luminosity L(x) is or where c is the speed of light, H(z) is the Hubble expansion rate, X I is a redshift-dependent conversion factor, and λ 0 and ν 0 are the rest-frame wavelength and frequency of the line, respectively.For simplicity, we will refer to the specific intensity as 'intensity' with the symbol I.
Given a function f (A) of a continuous random variable A with PDF φ(A), we can calculate its expectation value as f (A) = dA f (A)φ(A).As the PDF of IδV at position x is given by φ (IδV ) = φ(L)X I (x), the first and second moments of IδV are given by where N = nδV .Following the integration approach of Feldman et al. (1994), we find where δ D is the Dirac delta function.The second term in the previous expression assumes Poisson shot noise.
Similarly, the cross-correlation with the galaxy density contrast of detected galaxies, δ g (x) = n(x)/ n(x) − 1, is given by where the (second) shot noise term only contains galaxies that contribute to the galaxy catalog and the intensity map, i.e., only detected galaxies with non-zero luminosity in the target line, denoted in the expression as G g ∩ G I .We define the weighted Fourier transform as δI(k) = d 3 x w I (x)δI(x)e ik•x with a dimensionless weight w I (x), which can represent a survey footprint. 2We define the estimator for the intensity power spectrum where V box is the volume of a cuboid enclosing the survey used to compute the Fourier transform, so that PII (k) has the dimension of volume times intensity squared.We find that 2 Blake (2019) follow a similar approach but use weights in units of (intensity) −1 .
where P II is the power spectrum defined in Eq. ( 3), and the window function is defined as W I (k) = d 3 x e ik•x w I (x) I(x) .Similary, we define the estimator for the cross-power spectrum as with the dimension of volume times intensity, where '*' is the complex conjugate operator, and where the galaxy window function is W g (k) = d 3 x e ik•x w g (x) with a dimensionless weight, w g (x).Finally, we define the galaxy power spectrum estimator as where

Smoothing and noise
Limited observational resolution can be modeled by smoothing the intensity map.Suppose that the intensity is smoothed with a smoothing kernel D(x), i.e., Ĩs (k) = Ĩ(k) D(k).Then the factors in the power spectrum estimators change as Examples for smoothing kernels are a Gaussian or tophat smoothing in the line of sight (LOS) direction, mimicking the line-spread function or binning of intensity into frequency channels and a Gaussian smoothing in the angular direction, mimicking the beam smoothing or point-spread function where s and s ⊥ define the smoothing lengths parallel and perpendicular to the LOS, respectively.Similarly, k = µk and k ⊥ = 1 − µ 2 k are the components of the wave number parallel and perpendicular to the LOS, respectively, and µ is the cosine of the angle between the LOS and the wavevector k.Uncorrelated noise ∆I n with zero mean can be added to the intensity to model instrumental and sky noise.Therefore, I s,n (x) = I s (x) + ∆I n (x), and the noise is characterized by its variance ∆I 2 n (x) = σ 2 I (x), where σ I is the standard deviation of the noise in each voxel.This adds a term to the power spectrum estimator: The second term still holds for larger than infinitesimal voxel volumes δV .

Modeling the sky subtraction
The sky subtraction is a common step in the reduction pipeline of ground-based optical data to remove zodiacal light, aurora, and airglow (e.g., Wyse & Gilmore 1992).The sky spectrum consists of a continuum component and emission lines.With the assumption that the sky spectrum is homogeneous on scales of the focal plane size, one can estimate the sky foreground by taking an average of the spectra in 'empty' areas on the detector, i.e., those that do not contain an object above a fiducial detection limit.The sky foreground spectrum is typically subtracted from the data.
Sky subtraction also removes roughly the average intensity per redshift slice of the cosmological signal of interest within the scale used for the sky foreground estimation, for example, the focal plane radius (see, e.g., the discussion in Lujan Niemeyer et al. 2022).Therefore, the sky subtraction decorrelates fluctuations on larger scales.This does not affect the galaxy clustering of detected galaxies because the detection of galaxies is not influenced by the zero point on larger scales.However, it does affect the intensity auto-power and cross-power spectra even when all galaxies are detected because the intensity zero point is changed.While optical groundbased LIM experiments such as HETDEX suffer from a loss of large-scale power due to sky subtraction, groundbased LIM experiments in other wavelengths such as radio and submillimeter may also lose large-scale power due to similar issues.
We can model the effect of the sky subtraction on the intensity map by calculating the contribution of the line intensity to the estimated sky foreground.We smooth the intensity map with a two-dimensional spherical tophat kernel in the plane perpendicular to the LOS with the size of the area used for the sky spectrum estimation.Specifically, the estimated contribution to the sky foreground, i.e., the top-hat-smoothed intensity map, is where Dsky (k) is the Fourier transform of the twodimensional spherical top-hat kernel given by where a and b denote the directions perpendicular to the LOS and J 1 is the Bessel function of the first kind and of first order and s f is the radius of the area used for estimating the sky.
As Dsky (k) is real-valued, we obtain This shows that the power spectrum is suppressed on scales larger than s f .

GENERATING MOCK INTENSITY MAPS
This section describes the framework of the public code Simple to generate mock intensity maps.In a nutshell, the code follows these steps: 1. Generate a galaxy catalog in real and redshift space using the lognormal code of Agrawal et al. (2017).
2. Randomly assign line luminosities to galaxies following an input luminosity function.
3. Assign redshift and flux to each galaxy.Apply the input selection function to distinguish between detected and undetected galaxies.
4. Paint an intensity map and a galaxy density map using detected, undetected, or all galaxies in each map, and optionally smooth the intensity map.
5. Optionally, generate and add an intensity noise map, apply a mask or weights, and model the sky subtraction.
6. Calculate the auto-and cross-power spectra.
The main input parameters to Simple are cosmological parameters, a luminosity function, a linear galaxy bias for all galaxies, the central redshift of the box, and the rest-frame wavelength or frequency of the target line.If no tabulated matter power spectrum is provided as input, the cosmological parameters are used to generate the matter power spectrum using the Eisenstein & Hu fitting function (Eisenstein & Hu 1998).They are also used to calculate the luminosity and angular diameter distances in later calculations.
Along with the luminosity function, one has to specify the minimum luminosity to obtain a finite number of galaxies.This defines the number of galaxies to simulate with the lognormal galaxy catalog generator of Agrawal et al. (2017), which produces their positions, velocities, and redshift-space positions.This procedure assumes a flat sky and a single redshift.We assign a luminosity to each galaxy by randomly drawing from the luminosity function.Unless otherwise specified, we define the first axis of the simulation box as the LOS and assign the redshift to each galaxy according to its distance from the observer inferred from the position in this simulation axis.We use this redshift to convert luminosities to fluxes.Alternatively, a single redshift can be assigned to all galaxies.
We then apply a selection function based on an input flux limit above which a galaxy is detected, or on the target galaxy number density as a function of redshift.This produces a galaxy catalog with detection flags.We then calculate an intensity map and a galaxy number density map using the nearest-grid-point (NGP) assignment scheme.One can specify which galaxies contribute to the intensity and galaxy maps, i.e., detected, undetected, or all galaxies.
From the luminosities we calculate the intensity map, either in terms of a specific intensity or brightness temperature where ρ L = Ng i=0 L i /δV is the total emissivity in each voxel, i.e., the sum of the luminosities of all galaxies in that voxel divided by the voxel volume, and k B is the Boltzmann constant.
If specified in the input, the intensity map is smoothed with a Gaussian kernel perpendicular to the LOS, imitating the beam smoothing.Along the LOS, one can apply Gaussian or top-hat smoothing, imitating a line-spread function or binning in redshift/wavelength/frequency channels, respectively (see Section 2.2).One can add random Gaussian noise with the provided standard deviation per voxel σ I , and subtract the estimated sky foreground (see Section 2.3) from 2.9 1.5 Table 1.Summary of the differences between the low-z and high-z HETDEX mocks.
the intensity map.If a mask is specified in the input, the galaxy number density and intensity maps are multiplied by the mask.The mask is equivalent to the weights introduced in Section 2.1.Finally, one can calculate the summary statistics, i.e. the auto-power spectra of galaxies and the intensity and the cross-power spectrum.We compute the intensity and galaxy power spectra and the cross power spectrum using the estimators defined in Section 2.1, where V box is the volume of the simulation box.We use the fast Fourier transform (FFT) to calculate δI and δg , keeping only the independent modes, and calculate the corresponding k and µ values of the cells.For the quadrupole, we multiply the mesh δA(k) δB * (k) by the second Legendre polynomial evaluated at the µ values of the mesh.For each k bin, we collect cells whose k values fall into the bin and calculate their mean.We calculate the mean k value of each bin in the same way.

EXAMPLE MOCK SETUP: HETDEX
In this section, we introduce the setup for a HETDEXlike Lyα LIM experiment.HETDEX is primarily a cosmological galaxy survey that aims to map approximately one million LAEs through their Lyα emission line at z ∈ [1.88, 3.52] (Gebhardt et al. 2021).HETDEX uses the Visible Integral-field Replicable Unit Spectrograph (VIRUS; Hill et al. 2021) on the Hobby-Eberly Telescope (HET), which consists of 78 integral-field unit spectrographs (IFUs), each of which contains 448 fibers that are 1.5 in diameter, with a spectral resolution of 5.6 Å.
HETDEX observes a total area of 540 deg 2 without target pre-selection, expecting 460, 000 IFU observations.Because each IFU spans 51 × 51 in the sky, this amounts to an effective fill factor of the survey of f survey 0.17.The layout of the IFUs across the 18 diameter focal plane leaves an IFU-sized gap between adjacent IFUs and a hole in the center of the focal plane, which is used for other instruments (see Figure 1).Three six-minute exposures per HETDEX observation fill in the gaps between fibers, but not between IFUs, so that each individual HETDEX observation has a fill factor of 1/4.6.The HETDEX data reduction pipeline offers a sky subtraction mode that estimates the sky spectrum from all IFUs simultaneously, i.e., roughly from a discontinuous circular area 9 in radius.
We divided the redshift range of HETDEX into two parts with a similar LOS distance of 622 (624) h −1 Mpc for the low-z (high-z) sample.The low-z (high-z) sample covers z ∈ [1.88, 2.57] (z ∈ [2.57, 3.52]) and is centered around z = 2.22 (z = 3.04).At the mean redshift, the survey area of 540 deg 2 translates into a comoving area of 2.38 h −2 Gpc 2 (3.23 h −2 Gpc 2 ).We generate 7000 (9000) mocks of cubic volumes with side length 622 h −1 Mpc (624 h −1 Mpc).By averaging the power spectra over 7 (9), cubic boxes, we effectively get 1000 power spectra of a box that is 1.14 (1.09) times the size of the low-z (high-z) volume.In the last step, we multiply the covariance matrix of each redshift slice by the respective factor to correct for this oversampling.The largest accessible scale with this box size is k min = 0.01 hMpc −1 .
We set the cosmological parameters to the fiducial cosmology (see Section 1) and set the linear galaxy bias to b = 1.5.We simulate the Lyα intensity at the restframe wavelength of 1215.67 Å.In the following, we use the specific intensity I λ .For both redshift sections we adopt the Lyα luminosity function of the EWgt60 sample of Konno et al. (2016) for galaxies at z = 2.2, which is given by a Schechter function with L * = 4.87 × 10 42 erg s −1 , φ * = 3.37 × 10 −4 Mpc −3 , and α = −1.8.We set the minimum luminosity to 4 × 10 40 erg s −1 , based on Figure 4 of Gronke et al. (2015), and do not set the maximum luminosity.Chiang et al. (2013) investigated the impact of sparse sampling for galaxy redshift surveys, in particular HET-DEX, and found that the voxel size used for the power spectrum calculation has to be at least twice as large as the separation between IFUs to obtain an unbiased power spectrum.For this reason, we grid the simulation box with a resolution of 2 h −1 Mpc, i.e.N mesh = 311 (312), which corresponds to a Nyquist frequency of k Ny = 1.57hMpc −1 .Each voxel in our mock encompasses roughly eight HETDEX voxels, i.e., the length perpendicular to the LOS roughly corresponds to two IFU side lengths (102 ) and the LOS length roughly corresponds to two spectral bins (4 Å), which is smaller than the spectral resolution.This means that in our HETDEX-like experiment, one has to average the intensity of the fibers in each of the larger voxels.
Each cubic box at low-z (high-z) contains roughly 14.2 (14.4) million galaxies, both detected and undetected.We use a wavelength-dependent flux limit, enforcing that there are on average 2.5 detected galaxies per IFU, summed over the entire HETDEX LOS range (see Gebhardt et al. 2021).This corresponds to a mean galaxy number density of ng 2.0 × 10 −3 h 3 Mpc −3 .We rescaled the measured flux noise standard deviation from Figure 18 in Hill et al. (2021) to obtain a wavelength-dependent flux limit that satisfies this condition.
Figure 2 shows the resulting mean ng (z) in one realization of low-z and high-z HETDEX mocks.We compare this to the distribution of LAEs from the HETDEX public source catalog I (Mentuch Cooper et al. 2023), where we selected LAEs with the 'lae' flag in the 'source type' column.The shapes agree well.A decrease in the detected LAEs at z 2.7 is caused by a mask applied at the center of 50% of the HETDEX detectors and an increase in night-sky emission.The real number density in the high-z volume can differ if the luminosity function at z 3 is different from that at z 2. However, because of the good agreement with the detected galaxies in HETDEX, we continue to use the luminosity function at z = 2.2.
The high angular resolution of HETDEX allows one to mask individual detected galaxies without masking entire IFU-sized voxels.Therefore, we focus on the intensity map of only undetected sources.This intensity map does not share galaxies with the catalog of detected galaxies, so there is no contribution of shot noise to the cross-power spectrum between galaxies and intensity. 3 We do not apply beam smoothing perpendicular to the LOS because the point-spread function of VIRUS is smaller than the size of an IFU/voxel.We smooth the intensity map along the LOS with a Gaussian kernel with standard deviation σ = 2.38 Å(FWHM = 5.6 Å) to imitate the line-spread function, i.e., spectral resolution, of VIRUS (Hill et al. 2021).This corresponds to σ = 1.76 h −1 Mpc (1.27 h −1 Mpc) for the low-z (high-z) volumes.For the sky subtraction, we set the focal plane radius to 9 , which corresponds to a distance scale of s f = 10.0 h −1 Mpc (11.6 h −1 Mpc) for the low-z (high-z) part.
To add HETDEX-like noise, we transform the measured wavelength-dependent 5σ sky flux noise per resolution element per fiber in VIRUS (see Fig. 18 in Hill et al. 2021), 5σ F , into a specific intensity noise per fiber, σ I λ , by dividing by 5 × π(0.75 ) 2 × 5.6 Å .Then we divide by √ 3 × 448 to account for the averaging over fibers within an IFU in three dithers.Our voxels are shorter along the LOS than one spectral resolution element of VIRUS (5.6 Å; see Hill et al. 2021).Because the 3 Combining the catalog of detected galaxies with the full intensity map of all sources adds information with respect to the autopower spectrum of detected galaxies alone.This is true even if the intensity map is dominated by bright, detected galaxies (i.e., if the line luminosity function flattens at the faint end) because intensity fluctuations have a different (luminosity-weighted) bias from number-count fluctuations.We leave the study of this case for future work.
noise below this scale is correlated, the factor of 5.6 Å ∆λ to obtain the noise in the voxel with LOS length ∆λ is an overestimate of the correlated noise.We therefore do not apply this factor.Because the area of a voxel encompasses four IFU areas, we divide this noise level by N IFU (x), where N IFU (x) is the number of IFUs observed in the voxel at position x.In summary, we convert the 5-σ flux noise per resolution element 5σ F to the intensity noise at position x by calculating (29) This results in, on average, σ I λ = 2.9 (1.5) × 10 −20 erg s −1 cm −2 arcsec −2 Å−1 in each voxel in the low-z (high-z) boxes.
To obtain a mask on the coarse grid, we first generate a HETDEX-like mask with double resolution, so that each cell corresponds roughly to one IFU.We generate a mask of VIRUS-like tiles of ones on and zeros in between IFUs and keep the mask constant along the LOS.If we filled the entire area with observations, we would have a fill factor of f obs = 0.23, similar to the focal plane fill factor of VIRUS f VIRUS obs = 1/4.6= 0.22.To match the effective fill factor of f survey 0.17 due to sparser observations, we randomly remove 26% of the individual observations after applying the VIRUS-like mask.Then we downsample this mask to the same resolution as our simulated maps by averaging eight adjacent cells, which is equivalent to NGP assignment.The result is shown in Figure 1.We generate 7 (9) separate masks for the low-z (high-z) boxes and then average the power spectra of the boxes with different masks.
We calculate the power spectrum monopoles and quadrupoles in linearly spaced bins from k min = 0.04 hMpc −1 to k max = 1 hMpc −1 with ∆k = 0.04 hMpc −1 .We summarize the differences between the low-z and high-z HETDEX mocks in Table 1.

VALIDATION
In this section, we show that the results of our simulations in real space agree with the expected power spectra given in Section 2. We perform this comparison in real space because the lognormal algorithm precisely reproduces the input power spectra in real space, while the redshift-space power spectra deviate from the expectation in linear approximation (see Agrawal et al. 2017).
The setup of the validation mocks is almost identical to the low-z HETDEX-like mocks described in Section 4.However, we reduce the intensity noise by a factor of 30 and the flux limit by a factor of 5 to reduce the shot noise for validation, and do not apply sky subtraction.We apply one of the HETDEX-like masks to all intensity maps.We calculate the validation power spectra in real space and average over 1000 mocks with side length 622 h −1 Mpc.All other input parameters are the same.We calculate the intensity-intensity, galaxygalaxy, and galaxy-intensity power spectrum monopoles and quadrupoles, where undetected galaxies contribute to the intensity map.
To obtain the model, we evaluate the input power spectrum that includes the galaxy bias on a mesh with the same k and µ values as obtained from FFT of the mock maps.We also evaluate the damping functions from the intensity smoothing on this mesh.We calculate the window function by multiplying the mask by the mean expected intensity or the galaxy number density expected from the flux limit and the luminosity function as a function of redshift.We use FFT for the convolution of the power spectrum with the window function, add the shot noise, and multiply the result by the damping functions.Then we add the intensity noise power spectrum to the intensity auto-power spectrum mesh.We calculate the power spectrum multipoles of these results in the same way as for the mock (see Section 3).
Figure 3 shows the power spectrum monopole and the quadrupole in real space measured from the validation mocks.The shot noise (intensity noise) power spectrum is subtracted from the galaxy (intensity) auto-power spectrum monopole.The shot noise is not subtracted from the intensity auto-power spectrum because it contains information about the luminosity function.The quadrupole is affected by the anisotropic mask, the selection function, and the smoothing.It also shows the analytical prediction presented in Section 2, as well as the relative residuals between the power spectra of the mock and the model.The measured real-space power spectra from the mock agree with the model at all k modes.We also tested the validation in the cases where all or only detected galaxies contribute to the intensity map and find excellent agreement.

HETDEX FORECAST
In this section, we forecast the sensitivity of HETDEX to LIM power spectra.We use the mock setup described in Section 4 in redshift space and calculate the galaxy and intensity auto-power spectra, as well as their crosspower spectrum with and without sky subtraction.As explained earlier, only undetected sources contribute to the intensity map.
The upper panels of Figure 4 show the monopole and quadrupole power spectra in redshift space.The bottom panels show the cumulative signal-to-noise ratios (SNR) calculated as and the quadrupole (right) power spectra in real space.The shot noise (intensity noise) power spectrum is subtracted from the galaxy (intensity) auto-power spectrum monopole.The bottom panels show the relative residuals between the mock and model power spectra.The units are uI = I λ 2.8 × 10 −22 erg s −1 cm −2 arcsec −2 Å−1 and ug = 1.The shaded areas show the 1σ error of the mean given by the standard deviation of the different realizations divided by the square root of the number of realizations.The gray area shows deviations within 1%.
where k N denotes the maximum wavenumber considered.Here, ΘAB, AB, is the mean of AB, is the estimator of the power spectrum monopole ( = 0) or quadrupole ( = 2) calculated from the ith realization of maps A, B ∈ {I, δ g }.
Here, Q is the constant shot noise of the galaxy auto-power spectrum monopole, the noise power spectrum for the intensity auto-power spectrum monopole, and zero otherwise.The noise power spectrum of the sky-subtracted intensity auto-power spectrum was calculated by generating 10 pure noise mocks, performing the sky subtraction, calculating their monopole and quadrupole power spectra, and averaging them over the 10 realizations.For simplicity, we will leave out the subscript AB, from now on.We do not subtract the shot noise from the intensity auto-power spectrum because it contains information about the luminosity function.
−1 is the inverse of the covariance matrix whose elements are defined as Only elements of the covariance matrix up to the maximum k N are considered for matrix inversion to obtain As explained in Section 4, we multiply the covariance matrix by the factor f V = 1.14 (1.09) at low-z (high-z) to correct for slightly larger simulated volumes than the HETDEX volume.
We predict that an ideal HETDEX-like experiment can detect the galaxy-intensity cross-power spectrum monopole at 80σ (122σ) in the low redshift (high redshift) part of the survey despite the significant loss of large-scale power from sky subtraction.The quadrupole of the galaxy-intensity cross-power spectrum should be detected at 12σ (17σ).Although the intensity autopower spectrum monopole is dominated by intensity noise, the SNR is 11 (31) after removing the noise power spectrum at low redshift (high redshift).The intensity auto-power spectrum quadrupole should be detected at 4σ (9σ) after removing the noise power spectrum quadrupole., galaxy (orange), and sky-subtracted intensity (purple) auto-power spectra, and the intensity-galaxy cross power spectra without (green) and with (red) sky subtraction in the high-z volume.The quadrupole of the sky-subtracted intensity auto-power spectrum is not shown because of its large amplitude.The shaded areas show the square root of the diagonal elements of the corresponding covariance matrices.The units are uI = uI ss = I λ 4.6 × 10 −22 erg s −1 cm −2 arcsec −2 Å−1 and ug = 1.The dotted lines show the intensity noise and the galaxy shot noise power spectra.The bottom panels show the cumulative SNR up to a given k bin after subtracting the shot noise from the galaxy auto-power spectrum and subtracting the intensity noise power spectrum from the intensity auto-power spectra.The solid (dashed) lines correspond to the high-z (low-z) HETDEX-like mocks.
Figure 5 shows the correlation matrices of the data vectors Φ = (P 0 (k 0 ), P 0 (k 1 ), . . ., P 0 (k N ), P 2 (k 0 ), . . ., P 2 (k N )) (33) for the different power spectra that include the shot noise in the high-z HETDEX mock.The correlation matrix R of a vector Φ is given by R Φ mn = C Φ mn / C Φ mm C Φ nn , where the covariance matrix is given in Eq. ( 32).
The galaxy auto-power spectrum monopole has a large off-diagonal correlation, which increases with increasing k.This is mainly due to the mode coupling introduced by the convolution with the complicated window function.The galaxy auto-power spectrum monopole and quadrupole have a systematic low-level cross-correlation visible as stripes at constant k of the quadrupole.This may result from the imperfect integration of the µ values, in which the Legendre polynomials are no longer orthogonal (see Appendix D.2.3 of Agrawal et al. 2017).
The correlation matrices of the intensity auto-power spectrum and the cross-power spectrum without sky subtraction are dominated by diagonal elements.They have an off-diagonal correlation between the monopole and quadrupole power spectra at the same k, which is slightly higher than the noise around most off-diagonal elements.Although this monopole-quadrupole correlation is positive without sky subtraction, it is larger and negative with sky subtraction.The off-diagonal correlation for the intensity auto-power spectrum monopole is negligible because the covariance is dominated by the uncorrelated intensity noise.

DISCUSSION
The Simple framework is a simple and fast LIM simulation scheme with two main approximations.First, the galaxy distribution is modeled by a lognormal random field following the input power spectrum that inculdes a linear bias (Agrawal et al. 2017).The galaxies are obtained using Poisson sampling and the galaxy velocities follow the linear continuity equation.By construction, it is only accurate on the scales where the input power spectrum is accurate, which, for example, is only true on linear scales for the Eisenstein & Hu fitting function ( Eisenstein & Hu 1998).We may also underestimate the small-scale power spectrum because we do not model the 1-halo term, non-linear bias, or assembly bias.It is possible to improve the density and velocity distribution of mocks at the expense of decreasing the speed of the code.
The second approximation of the Simple framework is that the luminosity of a galaxy is randomly sampled from a luminosity function.Therefore, it does not depend on galaxy properties such as the star-formation rate or on its environment.This procedure misses any nontrivial connection between a specific galaxy population and their clustering, limiting any astrophysical analysis to the luminosity function.For instance, Simple can be used to study whether a given summary stastistic is an unbiased estimator of the luminosity function and its precision, but external computations must be done to connect this to astrophysical properties of interest.Alternatives like empirical determination of the astrophysical properties, semi-analytic models of galaxy evolution, or hydrodynamical simulations defeat the purpose of Simple, given their computational cost.

Limitations specific to the Lyα line
Our approach only models the emission from galaxies that is contained in the luminosity function.This does not capture all relevant physical emission processes of the Lyα line, as it neglects the recombination of ionized hydrogen and the collisional excitation of neutral hydrogen in the CGM and IGM (e.g., Dijkstra 2019).Photons originating from a galaxy can also be observed far away from its source due to scattering in the CGM and IGM (Byrohl et al. 2021).Byrohl & Nelson (2023) find that while photons produced through diffuse emission in the IGM are negligible in the global Lyα luminosity budget, photons that originate from the ISM or CGM of galaxies and scatter in the IGM contribute substantially to the Lyα luminosity budget.One can test the significance of the scattered contributions for LIM by modifying the luminosity function to include emission from the CGM calibrated on observations or by smoothing a small fraction of the intensity map to imitate scattering on larger scales than the voxel size.One can also add an intensity component that is directly proportional to the matter distribution to mimic diffuse gas emission.
Scattering of Lyα photons in the ISM strongly affects the escape fraction of Lyα photons and the emission peak wavelength (e.g., Hashimoto et al. 2013;Blaizot et al. 2023).The luminosity function contains the observed Lyα luminosity of detected galaxies, which consists only of photons that escaped from the ISM.However, the Simple framework does not account for the non-cosmological redshift of the Lyα line due to ra-diative transfer within the ISM, which may cause an anisotropic effect similar to the Fingers-of-God effect (Byrohl et al. 2019).It also does not account for radiative transfer in other environments, especially absorption, which may cause an anisotropic effect similar to the RSD, but with the opposite sign (e.g., Zheng et al. 2011;Behrens et al. 2018).These two effects can be modeled in an extension of the Simple code by shifting the Lyα line and calculating the optical depth along the LOS as a function of density and velocity, which can be obtained from a lognormal simulation.

Limitations of the HETDEX forecast
This HETDEX forecast is meant as an order-ofmagnitude forecast for the detectability of the crosspower spectrum of sky-subtracted intensity and detected galaxies.There are several possible improvements for a more accurate forecast in addition to the above Lyα modeling improvements: one can account for the evolution of the luminosity function; include galactic extinction; include interloper contamination of the intensity map and masking of bright interlopers; model the mask, selection function, and non-Gaussian noise more accurately using the HETDEX data; account for false positive galaxy detections; and model the intensity data reduction more accurately.We leave these improvements for future work.

SUMMARY AND CONCLUSIONS
We have presented the publicly available Simple code to quickly generate intensity maps that include observational effects such as noise, anisotropic smoothing, sky subtraction, and masking.It is based on a lognormal simulation of galaxies and random assignment of luminosities to these galaxies.Although this approach does not contain dependence of the line luminosity on galaxy properties or of the galaxy properties on the cosmological environment, it provides a fast and versatile way to generate mock intensity maps.These can be used to study survey systematics and calculate covariance matrices of power spectra.
We also derived an analytical model and showed that sky subtraction suppresses the power spectrum on scales larger than the focal plane size of the telescope.We validated the SIMPLE code by showing that the power spectra in real space agree precisely with those of the analytical model.This is the advantage of a lognormal mock, where the output power spectrum is designed to agree with the input power spectrum in real space (Agrawal et al. 2017).
As an application, we generated mock intensity and galaxy number density maps for a HETDEX-like LIM survey in redshift space that included a realistic mask, selection function, and intensity noise.We calculated the cross-power spectra of detected galaxies and the intensity of undetected galaxies after subtracting the sky spectrum and the respective covariance matrix in two redshift regimes within the survey.We predict that HETDEX will detect the cross-power spectrum monopole and quadrupole with high significance.
In summary, Simple has been designed to provide fast but reliable realizations that allow for changes in cosmology, luminosity function, and observational specifications with little effort.For example, one can numerically estimate covariance matrices for different setups.Because lognormal realizations overestimate the off-diagonal elements of the covariance matrix due to strong non-Gaussianity, Simple returns a conservative estimate especially for the galaxy auto-power spectrum (Blot et al. 2019).Therefore, this code is an excellent complement to slower, more physical simulations with the potential to guide the analysis of LIM surveys.

Figure 1 .
Figure 1.On-sky cut-out of a HETDEX-like mask applied to the HETDEX and validation mocks.The right panel is a zoom-in of a small area in the left panel.Each roughly hexagonal element with a hole in the center corresponds to one HETDEX observation with the VIRUS focal plane layout.The green squares in the right panel show the individual IFUs in the higher-resolution mask that was used to generate the downsampled, gray mask.The holes between observations show randomly removed individual observations to match the target fill-factor of HETDEX.

Figure 2 .
Figure 2. Mean detected galaxy number density as a function of redshift in one low-z and one high-z HETDEX-like mock (black).The dashed line shows the transition from the low-z to the high-z mock.The red histogram shows the distribution of detected LAEs in HETDEX with SNR > 5.5 from the HETDEX public source catalog I (Mentuch Cooper et al. 2023).

Figure 3 .
Figure3.Comparison of the mock power spectra (solid lines) with the analytical model (dashed lines) of the monopole (left) and the quadrupole (right) power spectra in real space.The shot noise (intensity noise) power spectrum is subtracted from the galaxy (intensity) auto-power spectrum monopole.The bottom panels show the relative residuals between the mock and model power spectra.The units are uI = I λ 2.8 × 10 −22 erg s −1 cm −2 arcsec −2 Å−1 and ug = 1.The shaded areas show the 1σ error of the mean given by the standard deviation of the different realizations divided by the square root of the number of realizations.The gray area shows deviations within 1%.

Figure 4 .
Figure4.Power spectra of the HETDEX-like mocks in redshift space.The top left (right) panel shows the monopole (quadrupole) of the intensity (blue), galaxy (orange), and sky-subtracted intensity (purple) auto-power spectra, and the intensity-galaxy cross power spectra without (green) and with (red) sky subtraction in the high-z volume.The quadrupole of the sky-subtracted intensity auto-power spectrum is not shown because of its large amplitude.The shaded areas show the square root of the diagonal elements of the corresponding covariance matrices.The units are uI = uI ss = I λ 4.6 × 10 −22 erg s −1 cm −2 arcsec −2 Å−1 and ug = 1.The dotted lines show the intensity noise and the galaxy shot noise power spectra.The bottom panels show the cumulative SNR up to a given k bin after subtracting the shot noise from the galaxy auto-power spectrum and subtracting the intensity noise power spectrum from the intensity auto-power spectra.The solid (dashed) lines correspond to the high-z (low-z) HETDEX-like mocks.