Accurate Fourier-space statistics for line intensity mapping: Cartesian grid sampling without aliased power

Estimators for $n$-point clustering statistics in Fourier-space demand that modern surveys of large-scale structure be transformed to Cartesian coordinates to perform Fast Fourier Transforms (FFTs). In this work, we explore this transformation in the context of pixelised line intensity maps (LIM), highlighting potential biasing effects on power spectrum measurements. Current analyses often avoid a complete resampling of the data by approximating survey geometry as rectangular in Cartesian space, an increasingly inaccurate assumption for modern wide-sky surveys. Our simulations of a $20\,{\times}\,20\,\text{deg}^2$ 21cm LIM survey at $0.34\,{<}\,z\,{<}\,0.54$ show this assumption biases power spectrum measurements by ${>}\,20\%$ across all scales. We therefore present a more robust framework for regridding the voxel intensities onto a 3D FFT field by coordinate transforming large numbers of Monte-Carlo sampling particles. Whilst this unbiases power spectrum measurements on large scales, smaller-scale discrepancies remain, caused by structure smoothing and aliasing from separations unresolved by the grid. To correct these effects, we introduce modelling techniques, higher-order particle assignments, and interlaced FFT grids to suppress the aliased power. Using a Piecewise Cubic Spline (PCS) particle assignment and an interlaced FFT field, we achieve sub-percent accuracy up to 80% of the Nyquist frequency for our 21cm LIM simulations. We find a more subtle hierarchical improvement in results for higher-order assignment schemes, relative to the gains made for galaxy surveys, which we attribute to the extra complexity in LIM from additional discretising steps. Python code accompanying this paper is available at github.com/stevecunnington/gridimp.

Parameter inference from LIM will require clustering statistics such as -point correlation functions to which models can be fitted.Optimal estimators often transform observations into Fourier space which allows faster computation.Furthermore, in the linear regime, Fourier modes will evolve independently, thus simplifying modelling since perturbations of different length scales are less correlated in Fourier space (Peebles 1980).
★ E-mail: steven.cunnington@manchester.ac.ukWhilst it is theoretically possible to transform all sky voxel coordinates into Cartesian coordinates and perform a direct summation to obtain the Fourier modes of the fluctuation field, the number of voxels within forthcoming intensity maps will render the direct summation method computationally unfeasible.For this reason, discrete Fast Fourier Transform (FFT) algorithms Cooley & Tukey (1965) will be needed which require sampling the LIM data onto a regular grid with Cartesian coordinates.A similar process happens with galaxy surveys whose catalogues of galaxy positions in 3D are too numerous for direct summation and thus must also be assigned to a grid to run an FFT.
There has yet to be a dedicated study into the regridding of LIM data with sky coordinates (R.A., Dec. ) onto a regular 3D Cartesian grid with lengths  x ,  y ,  z in ℎ −1 Mpc.Previous studies have approximated the LIM sky footprint to be rectangular, thus assuming equal voxel size everywhere and assigning Cartesian lengths to the map edges (e.g.Dillon et al. 2013;Wolz et al. 2017;Ihle et al. 2022).This allows an FFT on the map directly without the need for resampling the voxel intensities onto a new grid.This can be sufficient for small sky surveys, but as future LIM observations begin to span wider angles, this approximation becomes increasingly inaccurate.One way to resample the data is with a Monte-Carlo integration technique whereby sampling particles are distributed across the map and assigned the intensity values of the voxels they fall within.The position of these particles are then transformed to Cartesian coordinates and distributed into cells on a grid which can be FFTed.The Monte-Carlo sampling technique has mainly been used in simulations for doing the reverse; sampling a Cartesian cosmological cube into spherical sky voxels (e.g.Alonso et al. 2014).In Blake (2019), a Monte-Carlo sampling was utilised in simulations for validating the modelling of observational effects for Hi intensity mapping, but a dedicated study into its performance as a regridding tool is still warranted.
A concern for any regridding technique is how it will distort the true field properties from the raw observations.For example, aliased contributions can be introduced by perturbations across small separations no longer resolved by the finite size of the new grid cells.This has been previously studied for radio interferometry where converting visibilities into image space requires sampling the visibility values onto a grid so an FFT can be performed to produce the image.Early experiments used straightforward nearest neighbour schemes for assigning the visibilities into cells on a Cartesian grid (Hogg et al. 1969), but this was shown to produce unwanted aliasing.Convolutions with more sophisticated window functions have since been explored (McMullin et al. 2007;Ye et al. 2020;van der Tol et al. 2018;Monnier et al. 2022;Barry & Chokshi 2022), where often the instrumental response of the interferometer needs to be incorporated into the convolution kernel.
For auto-correlator (e.g.single-dish or single-receiver) experiments such as SPHEREx, FAST, MeerKAT, SKAO, COMAP SPHEREx Collaboraton et al. (2014); Li & Pan (2016); Santos et al. (2017); SKA Cosmology SWG (2020); Cleary et al. (2021), LIM observations are quite different.Whilst interferometers make angular measurements directly in Fourier space, auto-correlator observations provide time-ordered data which is made into sky maps in real-(or configuration)-space, represented by 3D arrays of voxels, each with sky angular coordinates and frequencies corresponding to the redshifted spectral line.Resampling these voxel positions which carry intensities is akin to resolved galaxy surveys which are also required to assign weighted point-like galaxies with 3D sky coordinates onto a regularly spaced Cartesian grid for an FFT.
The effects of discretised gridding on the power spectra of large galaxy catalogues have been well studied (Jing 2005;Cui et al. 2008;Yang et al. 2009;Colombi et al. 2009;Jasche et al. 2009;Sefusatti et al. 2016).Aliasing in the power spectrum arises from a periodic summation of contributions at multiples of the sampling frequency of the grid.Techniques can be used such as higher-order interpolations beyond a zero-order Nearest Grid Point (NGP) mass assignment (Hockney & Eastwood 1988).These are effectively low-pass filters that suppress the contributions from the unresolved small-scale modes but have to be increasingly non-local to entirely remove aliasing.Any mass assignment choice requires some correction to the power spectrum since the measured power spectrum from an FFT grid will not be equivalent to the true power spectrum of the raw field, but instead, one that is convolved with a window function  2 mas ( ) whose shape is determined by the mass assignment scheme adopted (Jing 2005).Other techniques such as interlacing the FFT grids have also been shown to remove the odd summations of the sampling frequency multipoles and suppress aliasing (Sefusatti et al. 2016).
In this work, we use simulations to examine the performance of Monte-Carlo sampling techniques for regridding LIM data into Cartesian coordinates.We present a modelling pipeline that approximates many of the distorting effects from the mapping and then regridding of the LIM field.We also explore some of the techniques used by galaxy survey analyses, such as higher-order mass assignment and interlacing, to see if they can be applied to the Monte-Carlo sampling particles to mitigate the aliased power they introduce.
The paper structure sees Section 2 detail the simulations we adopt and the Monte-Carlo-style sampling process used to transform the LIM data onto the final Cartesian grid.We demonstrate the distortions to power spectrum measurements caused by this process if no corrections are applied.Section 3 focuses on improving the accuracy of the power spectrum measurements relative to a modelled expectation.This is done by extending the modelling to include various discretiastion effects as well as including higher-order interpolation schemes for the sampling particles and interlacing the FFT fields.We discuss and conclude the results in Section 4.
We adopt some consistent language and notation conventions throughout the paper; referring to 3D discretised field bins as voxels when in sky (RA., Dec., ) coordinates, and cells when in Cartesian (x, y, z) coordinates.We also sometimes distinguish between angular and radial binning of sky data with pixels and frequency channels, respectively.We use the Greek vector  to refer to an angular position in sky coordinates at a frequency , and the Latin vector  for a position in Cartesian space.We refer to the input LIM fluctuation field (in Cartesian space) by  0 (), the array of HEALPix maps at each frequency with (, ), and the final resampled Cartesian grid, on which the FFT is performed, as  G ().

LIM REGRIDDING WITH MONTE CARLO SAMPLING
Normalized temperature fluctuations,  LIM , from line emitters are a biased tracer of the underlying dark matter field  m , in the post-epoch of reionisation era.This is simply expressed with where  is the linear bias of the lines brightness temperature field.
A LIM experiment will observe brightness temperature fluctuations given by where  are the absolute temperatures of the field, generally not accessible by a LIM experiment, and  is the mean brightness temperature.Surveys record the continuous temperature fluctuations,  LIM , binning time-ordered data into voxels on a sky map (, ), using a scheme such as HEALPix (Górski et al. 2005) pixelisation for all frequencies within the survey's wavebands.A Fourier-space analysis of the temperature fluctuations, such as a power spectrum estimation, will be imprecise if done directly on the sky maps  since its geometry is not a regular spaced grid in Cartesian space.A resampling process onto a Cartesian grid  G () is therefore required to coordinate transform the collection of voxels with angular pointings  (R.A. and Dec.) and frequency  (MHz), related to redshift through  =  0 / − 1, where  0 is the rest frame frequency of the target spectral line.The regridded data in Cartesian space can then be discrete Fourier-transformed following where  are weights used to optimise the analysis.In this work, for simplicity, we use unity weighting but inverse thermal noise weighting is more applicable in real data (see Blake (2019) for a derivation of optimal weighting for LIM).An estimator for any Fourier-space -point clustering statistic can then be applied to δG such as a power spectrum estimation, given by where  cell is the volume of a single cell in the regridded field. .This provides the 3D voxel array in sky coordinates (R.A., Dec, ).Using a Monte-Carlo method with sampling particles, the sky voxel intensities can be regridded onto a 3D Cartesian grid (red boundaries) in comoving space ( x ,  y ,  z in ℎ −1 Mpc) which can be Fast Fourier transformed.The method we adopt (detailed in Section 2.1) creates pseudo-random particles by generating  p per input sky voxel ( p = 10 for this example), with one particle at the voxel centre and the other  p − 1 distributed uniformly randomly within the voxel.
In Figure 1 we sketch out the process of sky pixelisation in a real LIM data reduction pipeline, followed by a regridding into Cartesian cells (following a Monte-Carlo sampling technique which we detail shortly).This shows how the continuous emission of the spectral feature is discretised into sky pixels, which are represented by the bluedashed lines.Depending on the instrument and map-making choices, these can be quite broad and the 0.3 deg size pixels shown are representative of the MeerKAT 21cm IM pilot survey data (Wang et al. 2021).Once pixelised maps have been created for each frequency channel of the instrument, the stacked maps need to be resampled onto the Cartesian grid.Example grid boundaries are shown in Figure 1 by the red-solid lines, the size of which are free to be chosen, which we discuss shortly.The plot also shows example Monte-Carlo sampling particles.Each particle is assigned a sampling weight, equal to the intensity of its host sky map voxel.In this example, we use 10 for the number of sampling particles per map voxel ( p ).These particles can then be trivially transformed 1 and their intensities populated onto the red grid in Cartesian space.A Cartesian grid cell value is determined by the average intensity from all sampling particles that fall within the cell (at least for the nearest grid point assignment scheme).Figure 1 only shows the 2D collapsed projection in angular space but the discretisation into voxels and resampling onto a Cartesian grid is 3-dimensional.

Monte-Carlo sampling summary
Here we detail the exact steps we follow for performing a resampling of the pixelised LIM onto a Cartesian grid.The scatter points in Figure 1 demonstrate the pseudo-random Monte Carlo sampling technique we adopt in this work whereby a constant number ( p ) of particles are assigned to each voxel in the sky map.The exact sampling process is summarised by the following steps: •  p particles per input sky voxel are generated with one at the voxel centre and the remaining  p − 1 particles randomly distributed with uniform probability density within the voxel.• Each th particle's mass is assigned the full amplitude (  ) of the sky voxel it originated from (this is later normalised in the new Cartesian grid).• The (R.A., Dec., ) sky coordinates of the particles are transformed into (x, y, z in ℎ −1 Mpc) Cartesian coordinates.• At this point, an interpolation scheme must be chosen which we explore in detail later.For now, we assume a nearest grid point assignment where the full mass of every th particle is added to the Cartesian cell  G (  ) it falls within.• We normalise all Cartesian cells by dividing by the count of particles that fall within it.Therefore the amplitude of each Cartesian cell will be given by where  p provides the option for weighting each sampling particle (in this work we assume  p = 1 for all particles2 ).An alternative choice would be to normalise the masses (or intensities) assigned to each sampling particle i.e. assigning (  )/ p .However, we found more consistent performance with unnormalised particles and a Cartesian grid cell averaging.We provide an analytical sketch motivation for this choice in Appendix A.
This approach differs subtly from previous work that also utilised a Monte-Carlo-style integration for regridding.Blake (2019) used a fully random process whereby a large amount of particles are randomly distributed across the survey sky map.We found our choice of using a constant number of sampling particles per input map voxel yielded more consistent results across different scenarios.This is likely because the intensities from the map are being divided by consistent quantities i.e. the number of sampling particles generated for each voxel is independent of , the position on the map.This is not necessarily the case in a fully randomised sampling process since an extra element of stochasticity is introduced.The sampling count in each voxel would also vary with  if we opted to generate a constant number of sampling particles per output Cartesian cell, causing position-dependent averaging problems.
The choice of cell size in Cartesian space (red grid in Figure 1), is a balance between ensuring enough sampling particles fall within all cells across the grid but not making the cells too large and further limiting the Nyquist frequency in the analysis.One could set  p = 1 which is equivalent to treating the map voxel centres as a discrete set of points whose coordinates can be transformed to Cartesian space.However, if this is done on too fine a grid, then there will be empty cells within the footprint of the survey, which creates ringing artifacts and mode coupling in the Fourier transform which are difficult to model.Figure 2 illustrates the empty cell effect with a toy example; a 20 × 20 deg simulated survey spanning 900 <  <1038 MHz (details of our main simulations are explained in the following subsection).Choosing to simply transform the sky map voxels and bin them into the nearest FFT cells on a  G = 128 3 grid results in numerous empty cells as shown by the top row.Each panel shows the count of sampling particles (just central particles at the voxel centres for this first case) within the cells.Instead of just transforming the voxel centre coordinates, we can instead implement a more abundant sampling process where we fill each voxel with  p randomly positioned sampling particles.This is shown in the second row where we use  p = 5 sampling particle per map voxel.This improves the empty cell problem but does not completely resolve it.A further simple improvement can be made by lowering the resolution of the Cartesian grid as done in the final row, halving the resolution to a  G = 64 3 grid.This results in a completely covered field and will avoid introducing additional ringing in a Fourier transform.
The disadvantage in decreasing the Cartesian grid resolution, as in the final row of Figure 2, is it artificially reduces the cutoff frequency beyond which small-scale (large-) Fourier modes are inaccessible.The Nyquist frequency of the grid is defined as  Nyq = /, where  = / G is the cell width of the grid, and  is the grid length down one side in ℎ −1 Mpc.If the underlying fluctuations contain power at  >  Nyq , which we know to be the case with cosmological perturbations (even if damped by a low-aperture telescope), then the resulting gridding guarantees a loss of fidelity.Thus restricting this is essential for regridding.The lost power beyond the Nyquist frequency and the distortion it can cause through aliasing, are discussed in detail in the latter sections.To ensure a completely sampled Cartesian grid therefore, a balance is sought between limiting the number of Monte Carlo sampling particles per voxel, to keep computational efficiency under control, and avoiding a final Cartesian grid with too low a resolution and an artificially limited Nyquist frequency.
The general survey geometry in Cartesian space is displayed by Figure 2, which demonstrates how the constant sky area through frequency, which results in a larger spatial area being covered at larger distances (lower frequencies), provides the truncated pyramid footprint.This highlights the imprecision of the assumption that the survey footprint is approximately cubical in Cartesian space, even for a relatively small (∼20 2 deg 2 , 0.37 <  < 0.58) survey, which is already close to being surpassed by pathfinder LIM experiments.As sky area increases further, wide-angle effects which make the assumption that the line-of-sight is in the same parallel direction for every pixel, is no longer valid.This is covered in the literature (Castorina & White 2018; Blake et al. 2018) and becomes a non-negligible issue for very wide surveys, complicating anisotropic effects such as redshift space distortions, in addition to distorting the volume.In this work, we assume these issues are sub-dominant and avoid simulating anisotropies such as redshift space distortions for simplicity, since a robust simulation would require a projection of the velocity field along the line-of-sight in the frequency direction, which will not remain consistently parallel along any Cartesian axis.

Simulated intensity maps
To validate the Monte-Carlo regridding process introduced in the previous sub-section, we construct a suite of simulations covering two different scenarios.To provide statistically significant validation, we demand a large number of realisations of each simulation version.
For this reason, we use fast lognormal mocks as the raw input field.The exact cosmology and emission line properties are not overly important, since in this work we are only examining what effect the discretisation and regridding processes have.This can be evaluated by studying the power spectra of the regridded fields relative to the input truth, which will be the main focus for the rest of the paper.For completeness, however, a brief overview of the assumed input cosmology and lognormal generation process is provided in Appendix B. We ensure the input starting simulation,  0 , has a high resolution ( 0 = 5123 cubic voxels) relative to the subsequent sky maps and resampled Cartesian grid to mimic the reality of mapping a continuous field.
The suite of simulations used for this work are generated by following the below presciption: (i) The size of the simulated survey is chosen in terms of its span in R.A, Dec and frequency (or redshift).From this, we compute the size of the Cartesian grid in which the survey footprint is entirely enclosed.
(ii) A lognormal mock  0 is generated on the  0 = 512 3 input grid.Every cell in this input grid has the correct Cartesian coordinates relative to an observer at the origin to be consistent with the survey footprint.We align the mock such that the line-of-sight through the centre of the survey is parallel with the Cartesian z-axis.
(iii) We construct a sky map, or lightcone 3 , from the input mock, using HEALPix pixels at each frequency channel of the mock survey.The size of pixels and frequency channels are chosen to emulate realistic surveys or test certain scenarios.This involves an inverted Monte Carlo integration process, similar to that which we implement for the regridding to Cartesian space where all particles are transformed from Cartesian to sky coordinates and placed within the nearest (R.A., Dec, ) map voxel.We ensure each map voxel is well-sampled.Distributing sampling particles into sky coordinates to build a map is not dissimilar from real LIM map-making where calibrated time-ordered data stamps are collated into sky voxels.This can provide another source of aliasing which we discuss later.
(iv) From the sky maps, we execute the sampling process (outlined in Section 2.1) to provide the regridded field  G from which we can study the changes relative to the input field  0 .For the final regridding, a choice of resolution (number of FFT cells  G ) is made which we initially set to approximately double the size of the HEALPix sky voxels to ensure good sampling.Later in the paper, we revist the choice of  G , increasing it so the Nyquist frequency of the analysis is pushed to higher-.
Table 1 summarises the two simulation versions we adopt in this work.We begin with the Cubic voxels simulation version which we artificially design to have approximately cubic dimensions for the voxels and Cartesian cells in the final grid  G .This is to provide an isotropic Nyquist frequency, to simplify the evaluation of aliasing effects.In reality, this may not be the case for real LIM experiments which can have very different resolutions along the line-of-sight (through frequency), relative to their broad angular pixels.We explore the consequences of this anisotropic Nyquist frequency in the other Fine channel simulation version.
The choices made for the simulated surveys should not bear a huge influence on the generality of our results and conclusions.Nevertheless, we have chosen to approximately emulate a pathfinder 21cm intensity mapping experiment at radio wavelengths such as MeerKLASS (Santos et al. 2017).The nside = 256 resolution of the map pixels we adopt throughout is approximately consistent with the 0.3 deg MeerKLASS pixel size (Wang et al. 2021), approximately 1/3 of the beam size.We also do not expect the choice of using a HEALPix pixelisation scheme to have an impact on the generality of our results.The current calibration pipeline for MeerKAT pilot data (Wang et al. 2021) uses the Zenith Equal Area (ZEA) map projection method (Astropy Collaboration et al. 2013).Thus we anticipate imminent validation of different pixelisation schemes which preceded the same Monte Carlo sampling techniques we outline in this work.
The 20 × 20 deg survey size is arbitrarily chosen, but again is close to the current size of MeerKLASS pilot survey observations.This creates 7182 filled sky pixels covering an approximately square patch.To generate cubic voxels we chose the frequency range and number of channels such that they span approximately the same depth  z as the spatial span in  x and  y on the Cartesian grid.We can also implement a smoothing to the field perpendicular to the line-of-sight to emulate the beam profile of the telescope.For this, we use a Gaussian profile assuming the dish size for the Square Kilometre Array Observatory (SKA Cosmology SWG 2020), MeerKAT's successor, which has  dish = 15 m.This defines the FWHM of the Gaussian beam profile in radians from  FWHM ≈ / eff  dish , where  eff is the median frequency of the survey, hence we assume a frequency-independent (constant size) beam.
Figure 3 shows a realisation from the Cubic voxels simulation When defining the size of the input Cartesian grid given the chosen survey size, we add some zero-padding around the borders.This is to allow for higher-order interpolation schemes (investigated later) in the map-to-grid resampling step, where particle assignment can be non-local, hence voxels around the edge will generate sampling particles outside the original footprint.This is seen in the final  G panel where there are some empty pixels due to the zero-order NGP assignment currently adopted.For higher-order schemes, this space becomes more filled.
The Monte-Carlo resampling technique appears to be working byeye.Original structure from the input  0 can be seen in the final  G grid, albeit at a lower resolution.The edges of  G appear to show a greater loss in fidelity but this is merely a plotting effect caused by the averaging along the line-of-sight axis.Since the sky map  cuts a conical (or truncated pyramid) footprint from the input  0 , there will be more spatial coverage at higher-z.Closer to the edges of the z-averaged map in  G there will be fewer slices covering the larger spatial separations, hence a lower sampling of the input  0 and an apparent poor recovery of its structure.

Biased power spectra from mapping and regridding
As briefly discussed in the introduction, previous pathfinder LIM survey analysis often avoided resampling the sky maps altogether by assuming the maps are approximately rectangular in Cartesian space.Before assessing the performance of our resampled maps, we briefly comment on the accuracy of this assumption for power spectrum measurements, to serve as motivation.We emulated the simplified approach by assuming that the -stacked HEALPix maps from our 20 × 20 deg 2 Cubic voxels simulation formed a rectangle in Cartesian space, with equal-sized voxels throughout.We assigned comoving dimensions to the edges of this rectangle, calculated with  x =  c ( eff ) RA , where  RA = 20 deg (in radians). y is calculated similarly but for Declination.The line-of-sight length is given by  z =  c ( max ) −  c ( min ).We then Fourier transform the sky map with these dimensions and measure the power spectrum.We found this approach caused a >20% bias across all scales relative to the input model and was most pronounced at large scales (small-) reaching ∼30%.The discrepancy increases at larger scales due to the true conical geometry of the survey, which means the approximation that separations on the sky are equivalent to separations in Cartesian space, is inaccurate.Whilst this approximation has been sufficient for many surveys with small sky coverage and low signal-to-noise, future surveys will require a more accurate approach, which is the motivation behind our resampling techniques.
The regridding of the LIM simulations from sky coordinates onto a Cartesian grid appears to be working well by eye, as seen by Figure 3.However, the true test of performance will come from a statistical analysis of the field relative to a model with the same fiducial cosmology used to generate the input field  0 .The model we compare to therefore follows the input power spectrum discussed in Appendix B and is interpolated over the same grid of modes as δG ( ); where  m ( , ) is the matter power spectrum for the fiducial cosmology at redshift .We include the option of modelling the telescope's beam which damps perpendicular modes  ⊥ .For the frequencyindependent Gaussian beam, this damping is modulated by the beam size  beam ( eff ) =  c ( eff ) FWHM ( eff )/2 √ 2 ln 2. We examine the effects of the beam and discuss further potential effects from more realistic, complex beam patterns later in the paper.The ′ notation on  ′ mod in Equation 6 is to notify the different iterations of model improvements which will be useful later when extending this base model.
The *  G ( ) part of Equation 6represents the convolution of the model with the survey window function, and optional optimal weighting.This is necessary since mapping the input field onto the survey footprint leaves empty pixels in the final regridded field  G outside the survey footprint, demonstrated by the white space in the boarders of Figure 2. To correct for this in the model, we must convolve with a survey window function which we simply define as  G () = 1 for every non-zero cell on the grid covered by the survey footprint, and zero otherwise.More complex survey window functions can be implemented to account for non-uniform survey completion.Some apodization (or smoothing) of the footprint edges may help any discontinuous edge effects in the mask but is not something we investigate in detail in this work.Convolving the model with this binary window function, and normalising by it and any optional survey weighting can be expressed for a generic model power spectrum  mod with The model from Equation 6is averaged into the same spherical shell bins as the power estimation for the simulations, which we chose to be linear in .
Figure 4 shows the power spectrum estimation of the regridded field for the Cubic voxels simulation (solid line), initially with no beam i.e  beam = 0.For the simulation results, we average over the 100 realisations for each version.We found this was enough realisations to reach converged results for all tests in this paper.The shaded regions show the 1 range over the realisations to indicate the scatter from the simulations.Since the main focus in this work is optimising accuracy, we do not show these uncertainty regions in subsequent plots, since their size is driven by the survey size choice and not any analysis treatments we implement.The red dotted vertical lines (also in subsequent similar plots in the paper) represent the Nyquist frequency  Nyq and 0.5 Nyq of the Cartesian grid for the perpendicular x and y directions.For the Cubic voxels simulation, the z-direction Nyquist frequency is also equal to the x and y directions by design.For the Fine channels version, the z-direction has a higher Nyquist frequency  Nyq = 0.865 ℎ Mpc −1 , due to the increased resolution along the line-of-sight.
The model (shown by the black lines in the top-panel of Figure 4) is directly compared to in the bottom panel, where perfect agreement would see ()/ mod () = 1.There is poor agreement for the blue solid line which has no further treatment than what we have already discussed.We discuss the results from the Fine channels version (dashed lines) shortly.The modelling should be improved by recognising that in the construction of the  G () field, we performed an interpolation of the sampling particles using a NGP assignment scheme.As identified in Jing (2005), assigning a particle distribution that traces a continuous tracer field () onto a grid ( G ) is equivalent to a convolution of the continuous field with a mass assignment scheme window function  mas ().In Fourier-space this convolution becomes a multiplication, thus our LIM data resampled onto the grid  G can have the leading-order effects corrected by simply dividing the field by the Fourier transform of the window function Wmas ( ) (Sefusatti et al. 2016) For the NGP assignment scheme that we currently use, we are simply .We show the case where no compensation for the NGP interpolation is applied to the regridded field (blue lines) and the case where compensation is applied using Equation 8 and Equation 9(orange lines).We do not include the effects of the telescope beam in these results.Measured power spectra PG are averaged over 100 simulation realisations.The 1 scatter from the 100 simulations is shown by the shaded regions.
discretising particles into cells of size  = / G .This is modelled by a top-hat function which the Fourier transform of is a sinc function, sinc() = sin()/.For a 3D field, these window functions remain separable i.e.  () =  () () () (Hockney & Eastwood 1988) and the interpolation window function for NGP is given by The orange results in Figure 4 show the results from applying this compensation directly to the field defined by Equation 8. We see how this is providing some correction for the damping due to the NGP interpolation (orange solid line) and now power is enhancing as the Nyquist frequency is approached.This shows the impact of aliased power caused by small-scale separations between sampling particles not resolved by the cells in  G .Aliasing will still be present without the compensation in Equation 8 (blue lines), but the results are overly dominated by the smoothing from the particle assignment.The following section is devoted to how further improvements can be made both to mitigate this aliasing and to also include modelling for prior discretisation to the field in the map-making stage i.e. performing  0 () → (, ).
Before exploring further modelling and improved interpolation techniques, we consider the case where the LIM frequency channels have better resolution than the pixels, which will emulate a realistic situation for many LIM sureys.We use the Fine channels simulation version (see Table 1), where the number of frequency channels is more than trebled, increasing the maximum  ∥ accessible.We also increase the redshift range probed which is equivalent to increasing the depth of the survey in comoving space.We do all this in such a way to keep the central redshift  eff unchanged to aide result comparison.This departs from the simple case of having approximately cubic voxels and cells.Now there will be a different Nyquist frequency depending on the axis direction which will have interesting effects for aliasing.
The results for the Fine channel simulation are shown by the dashed lines in Figure 4.For the case where we compensate for the NGP interpolation (orange dashed line), there is less enhanced power compared to the Cubic voxels (orange solid line).This is because, there will be less aliasing along the line of sight when there is higher frequency resolution, hence aliased contributions to the spherically averaged modes will be suppressed.Analysing the 2D cylindricallyaveraged power spectrum decomposed into ( ⊥ ,  ∥ ) modes demonstrates this point.Figure 5 shows the 2D power for the cubic voxels (left panel) and fine channel case (right panel).We draw the same Nyquist frequency (red dashed line) from Figure 4 for reference.This is defined as the minimum Nyquist frequency which comes from the perpendicular dimensions, hence is the same for both panels.The grey area for the Cubic voxels is an area of -space inaccesible by the grid resolution, but can be accessed by the Fine channels simulation.The blue regions show the area of -space where aliasing dominates.For the cubic voxels, power enhances relatively isotropically as |  | = √︃  2 ∥ +  2 ⊥ approaches the Nyquist frequency.For the fine channels however, aliasing does not dominate along the parallel axis until much higher  ∥ .Therefore, aliasing in the spherically averaged modes around the reference scale (the minimum Nyquist frequency  Nyq ) will wash out as seen in the 1D power for the Cubic voxels in Figure 4. Aliasing remains present in regions of -space even with fine channel resolution and thus we still require an extension on the techniques we have used so far, which we explore in the following section.Before that however, we consider the impact from the telescope beam.

Why aliasing remains even with a large instrumental beam
To complete this section, we consider the impact from the telescope's instrumental beam.So far the results presented have been for the case with no perpendicular smoothing to the field caused by the beam.One might naively assume that the smoothing will completely suppress aliasing effects, however, this is not necessarily the case since smoothing from the beam comes before the sampling particle assignments.Therefore, any aliasing effects originating from particle assignment (or time stamp binning in map-making), will still inject aliased power, artificially enhancing the beam-damped power spectrum.Since fluctuations still exist beyond the Nyquist frequency, albeit damped by the beam, the power spectrum will still be distorted by aliasing.For this reason, we also expect aliasing to remain in the presence of any smoothing caused by instrumentation.Higherfrequency LIM experiments e.g.SPHEREx (SPHEREx Collaboraton et al. 2014), will have a less precise spectral resolution, which can be characterised by a radial damping, similar to the beam but in an orthogonal direction to it.Hence, the effects of aliasing from resampling the field are likely to be an issue for any LIM experiment aiming for precision cosmology.
Figure 6 shows the impact from the beam in the regridded power spectrum measurement.We smooth the input field  0 with a  FWHM = 1.15 deg beam, equivalently  beam = 9.6 ℎ −1 Mpc.This is not a perfect emulation of reality as the beam will actually act perpendicularly in the curved sky which will not be perfectly perpendicular in the Cartesian frame.We discuss the interesting impact from the beam further in Appendix C. Despite using an oversimplified Gaussian beam applied in the Cartesian frame, it provides a sufficient test to show that fields smoothed by a beam do not escape aliasing.Since we can trivially model the impact of this simple beam by including  beam in Equation 6, we continue to implement it in the rest of the paper for some extra completeness, but defer a more robust implementation of it to future work, as also discussed in Appendix C.

MITIGATING ALIASING AND MAPPING EFFECTS
In Section 2 we introduced the Monte-Carlo resampling which achieved a decent accuracy on large scales but introduced some effects at smaller scales.Even after correcting for the smoothing to the field (Equation 8) from the interpolation of sampling particles, a discrepancy between power spectrum measurement and model remained.This will be caused by aliased power introduced by the sampling particles as well as the additional discretisation step from the prior voxelisation into the sky maps.In this section we extend the modelling used in the previous section to approximate the effects from the sky mapping, with the aim of improving the accuracy of the final power spectrum estimation.Additionally, we explore some higher-order interpolation schemes for the Monte Carlo sampling particles beyond the simple nearest grid point (NGP) assignment scheme used by default so far, which we showed introduces aliased power.

Modelling the sky mapping
An exact model for the additional step of discretisation from mapping the "continuous" field  0 into (R.A., Dec. ) voxels on (, ), would require consideration of the curved sky and the spherical shell voxels that the temperature fluctuations are discretised into.However, following Blake (2019), by making a narrow-angle assumption, we can make the approximation that the HEALPix pixelisation and frequency channel binning are separable effects on  ⊥ and  ∥ respectively.
For the angular pixelisation, we utilise the HEALPix pixel window function4  pix (ℓ) which describes the damping to the harmonic angular power spectrum  ℓ →  ℓ  2 pix (ℓ), given by where is the spherical harmonic transform of a pixel in terms of the spherical harmonic functions  ℓ .For the final Cartesian grid  G (), we assume ℓ ∝  ⊥ and model the damping caused by the pixelisation as ( ) → ( ) 2 pix (ℓ eff ), where ℓ eff =  ⊥ / c ( eff ), and  c is the comoving distance to the central redshift of the survey.Thus the damping due to the angular pixelisation can be expressed with There is also a smoothing effect along the line-of-sight from the sky mapping due to the discrete frequency channels of the LIM survey.This acts as a top hat function along the frequency axis, which we approximate as being consistently aligned along  ∥ .For channels with depth Δ chan , this is modelled with a sinc function in Fourier space given by

Additional aliasing from sky mapping
During the creation of the HEALPix sky maps (, ), aliased frequencies will be inserted into the sky map, from the default NGP assignment in the map making.Since the map undergoes further regridding onto the Cartesian field, the original field in which these map-based effects arise is not accessible, thus we opt to apply a further correction in our modelling instead.There is scope to address this problem at the map-making stage with higher-order assignment of the time-ordered data into HEALPix (or an alternative scheme's) map pixels.Due to the extra requirement of simulating time-ordered data to test this, we leave this investigation to future work.To mitigate the bias we instead follow an approach closely aligned to that outlined in Jing (2005) (see also Blake et al. (2010) for an application with real data), whereby the power spectrum is summed over modes displaced by integers of 2 Nyq to correct for this discretisation at sky map level.This requires the grid onto which the model is interpolated be representative of the sky voxel sizes.Once again, this is not perfectly possible due to sky curvature but the pixel widths and channel depths can be approximated which is simply done by choosing the number of Fourier cells ( Gx ,  Gy ,  Gz ) to equal the number of map pixels and frequency channels.If the overall Fourier grid size is calculated in a way that is consistent with the Cartesian size of the map space, the FFT cells will be approximate to the voxel sizes.We can then perform a summation of displaced modes onto the model which should approximately suppress the map voxelisation effects.This has the added benefit of sampling onto a more resolved grid than before, where the cell size was approximately double the map voxel sizes.Providing we avoid the empty cell problem as presented by Figure 2, this will have the beneficial effect of not artificially decreasing the Nyquist frequency below that which is defined by the map voxel sizes.This keeps discrepant effects contained to smaller scales which can be disregarded in a final analysis.Any aliasing caused by the Monte Carlo sampling particles onto the final Cartesian grid can still be corrected for at the field level, so this approach partitions the effects with the sky mapping discretisation effects handled in the modelling, and Cartesian regridding effects addressed at field-level.We summarise the final modelling approach with the below two steps; (i) The first extension to the original  ′ mod base model (Equation 6) is to include damping from the map voxel discretisation, achieved with the smoothing functions in Equation 12 and Equation 13; As before, this model is convolved with the survey window functions following Equation 7. The results from this relative to the original model are shown in Figure 7 (orange line).
(ii) The final extension applied to the model, after convolution with the survey window function, is the correction to mitigate any bias from the NGP assignment into sky map voxels.For this we adopt a correction following (Jing 2005) which can equivalently be applied with the summation over modes; where  is a vector of integers, a summation over which displaces  by 2 Nyq in all dimensions.In principle, this can extend indefinitely but it is usually sufficient to just consider the leading order displacement.We explored this and found no improvement when going beyond a 3 3 grid of  = {−1, 0, 1}.In Equation 15, the model requires interpolation onto a grid whose cells are a close approximation to the sky voxel size, so that the displaced Nyquist frequencies and NGP window function are emulating the sky map effects.
As before, the plotted model power spectra are all averaged into the same spherical shell bins as the power estimation for the simulations.
Results for the full model, including the summation over modes,  6), ( ′′ mod Equation 14) and ( ′′′ mod Equation 15).Green shaded areas show the sub-5% and 1% accuracy regions.Measured power spectra PG are averaged over 100 simulation realisations.Blue line shows the limited model from Section 2. The orange line includes damping from the discretisation of the field into the map pixels and frequency channels.The red line further includes the summation over modes.are shown by the red line in Figure 7.There is no clear improvement by introducing the summation over modes, but this is still for the case where the zeroth order NGP assignment scheme is being used for the Monte Carlo sampling particles.Extending beyond this should mitigate the enhanced power shown by the red line which will be a result of residual aliasing, mostly from the Monte Carlo resampling.

Higher-order mass assignment schemes
We have already shown in Section 2.3 how the smoothing effects from the sampling particle assignments can be modelled by a convolution with the window function (Equation 9) in Fourier space, which is simply a product of separable top-hat functions along each dimension.As shown in Hockney & Eastwood (1988), iteratively convolving with further top-hat functions is an effective approach for delivering higher-order interpolations to suppress aliasing.These higher-order routines generate 3D cloud shapes and the portion of density from each shape that falls within a particular cell determines the value assigned to that cell from the initial mass introduced.An easier way to illustrate the interpolations is with 1D assignment function shapes.The first-order function, cloud-in-cell (CIC) is piecewise linear, as apposed to NGP which is piecewise constant.Beyond this, the triangular-shaped cloud (TSC) scheme becomes a continuous first derivative, with increasingly higher-order derivatives beyond this.Following Sefusatti et al. (2016), we test up to the third-order interpolation scheme now commonly referred to as the piecewise cubic spline (PCS).
We express the 1D assignment function shapes for the four schemes we consider below as a function of , the normalised particle distance from a cell-centre  G  , i.e.  = ( G  − )/, where once again  = / G is defined as the cell size; • Nearest Grid Point (NGP): This describes the practical interpolation of the sampling particles onto the regridded field  G . Figure 8 (top panel) displays these assignment functions in configuration space as a function of , showing the increasing non-locality of the higher order functions, hence computational demand increases with , since more neighboring cell assignments require calculation for each particle.In this work, we utilise the public python package pmesh5 , also used for other cosmological analysis codes such as nbodykit6 (Hand et al. 2018).The window functions in Fourier space for each of these schemes is then simply the sinc functions raised to the power  which depends on the choice of mass assignment scheme; where  = {1, 2, 3, 4} for NGP, CIC, TSC and PCS respectively.We display these window functions in the bottom panel of Figure 8 which shows the desired effect of making the assignment more local in Fourier-space and should suppress aliasing contributions.The price to pay for making the assignment more local in Fourier space is the configuration-space assignment is increasingly non-local (as shown by the top panel), thus increasing the computational cost.Even more non-local interpolations can be considered e.g. using Daubechies wavelets (Cui et al. 2008).Results in Hand et al. (2018) showed that despite their improved performance at smaller scales, some larger scale discrepancy is introduced, hence we do not investigate these schemes in this work.
In Figure 9, we show the accuracy of the higher-order interpolation scheme results (dashed lines), relative to the full model case  ′′′ mod Equation 15.We discuss these results in more detail soon, but first, we introduce the additional field-level treatment of interlacing.

Interlacing
Sefusatti et al. ( 2016) investigated using interlaced FFT grids as a means for suppressing aliasing, as initially introduced in Hockney & Eastwood (1988).Interlacing involves interpolating onto a second FFT grid with the cell positions all shifted by half a cell size (/2) in all directions, so the discrete Fourier transform of this second field will be given by The interlacing technique is inspired by the "butterfly" operation which the FFT is based upon.As analytically shown in Sefusatti et al. (2016), the interlaced field will cancel out the contributions from the odd integer phases of the aliased modes which should include the leading order contributions.
Figure 9 presents the higher-order interpolation scheme results, both with and without interlaced fields.Relative to the full model case  ′′′ mod (Equation 15).This shows how the enhanced interpolations are generally suppressing the aliasing and allowing a sub-percent agreement with the model power.
There is a less hierarchical order of results in Figure 9 than expected, and as seen in equivalent investigations for galaxy catalogues  15) for the higher-order mass assignment schemes, both with and without interlacing.The green shaded area shows the sub-1% accuracy region.Measured power spectra PG are averaged over 100 simulation realisations.
(e.g.Jing 2005; Sefusatti et al. 2016;Hand et al. 2018).In these studies, where a straightforward gridding of point-like galaxies is performed, the higher-order schemes deliver a clear improvement in accuracy.For our results however, we appear to have diminishing returns and it even appears that the TSC scheme performs better than the higher-order PCS at high- (although this is not the case at larger linear scales as we will show).We suspect that the lack of conclusive successive improvement with increasingly more sophisticated schemes is due to the complication we face from the additional discretisation step caused by the mapping into sky voxels.This introduces a NGP interpolation of the continuous field, which we then apply a further interpolation to in the Monte-Carlo resampling to the Cartesian grid.Whilst our approximate models for these effects are performing well, they will not untangle the complicated mix of particle assignments that we believe make the final higher-order interpolations less efficient at suppressing residual aliasing.We were able to test this hypothesis in a Cartesian-only investigation where the sky mapping step was done onto another Cartesian field, thus avoiding curved-sky complications.This toy case reproduced similar results where some small residual aliasing remains from the extra mapping step and there were still no clear hierarchical improvements as higher-order interpolation schemes were implemented.We leave further investigation to future work into whether an enhanced mapmaking process can be introduced that avoids this consequence and suppresses aliased power further.
The results from Figure 9 sufficiently demonstrate how sub-percent accuracy for regridded LIM data is achievable far up to the Nyquist frequency of the Cartesian grid.Despite a better performance of TCS for high-, on medium to large scales (0.07 <  < 0.25 ℎ Mpc −1 ), which are of most interest for cosmological experiments since perturbations are more linear, the best performance comes from the interlaced PCS scheme as one would expect.On these scales, there is slight evidence for the hierarchical ordering of results.This is shown by Table 2 where we show the average bias across these scales for each scheme.We see a near consistent gradual decrease in the bias for higher-order schemes with interlacing.We chose this scale range to avoid the largest scales where cosmic variance will start to dom- inate and potential masking/modelling effects may be less accurate, and also the smallest scales where the residual aliasing is having the most impact and would be discarded in real analysis.The Table 2 results thus motivate the use of higher-order interpolations with interlaced fields, for high precision cosmology with LIM data.We also tested these processes on the Cubic voxels simulation where aliasing is more prevalent at lower- and found equally conclusive results concerning the need for damping and aliasing treatments.

DISCUSSION & CONCLUSION
Line intensity maps provide matter tracing fluctuations across a 3D space in sky coordinates, where the intensity in each voxel corresponds to a position in (R.A., Dec., ).Analysing these fluctuations in Fourier space, e.g. with a 2-point power spectrum estimation, requires the sky map to be resampled onto a regularly spaced Cartesian grid to allow a computationally viable FFT to be performed.
The voxel coordinates can be trivially transformed to Cartesian coordinates and their intensities binned into the new Cartesian cells.However, if the choice of resolution for the new Cartesian grid is too high, empty cells will appear in the field and create unwanted ringing effects in the Fourier transform.Too low a resolution artificially reduces the Nyquist frequency, limiting the scale range in which the analysis can be performed.
In this work, we showed how a Monte-Carlo-style integration technique, whereby a number of sampling particles are assigned to each map voxel and then transformed to Cartesian space, is an effective way to resample LIM data.However, we showed how small-scale separations between the sampling particles can inject unwanted aliasing into a power spectrum estimation.Inspired by previous work (Jing 2005;Cui et al. 2008;Sefusatti et al. 2016;Hand et al. 2018;Blake 2019) which addresses a similar effect but mostly from gridding of point-like galaxies, we presented corrections that can be applied both at field-level and in the modelling to unbias the effects from the regridding to Cartesian space.We then explored higher-order mass assignment schemes and interlaced fields as techniques to suppress the aliasing introduced from the resampling particles.We found that these can be effective solutions and allow sub-percent accuracy relative to the modelled expectation.We found there was a less conclusive hierarchy of results for the schemes compared to other results from galaxy survey simulations (e.g.Jing 2005; Sefusatti et al. 2016), but attributed this to the more complex process LIM data goes through, containing an additional interpolation step in the map-making where the continuous data is discretised into sky voxels.The complex mix of this initial interpolation, followed by the further interpolation at the Cartesian regridding stage is responsible for a less efficient im-provement in accuracy from higher-order schemes.However, our results from Figure 9 and Table 2 show that the PCS scheme with an interlaced field delivered the most accurate results, well within sub-percent agreement for the scales of interest.
Whilst the process of resampling to a Cartesian grid can be avoided by using a harmonic-space  ℓ power spectrum or a configurationspace correlation function, these estimators come with additional challenges and are less utilised in current LIM analyses.For example, a Fourier-space power spectrum analysis is adopted in Masui et al. (2013); Anderson et al. (2018); Keenan et al. (2022); Wolz et al. (2022); Ihle et al. (2022); Cunnington et al. (2022).In all these works some transformation to Cartesian space is required.In many cases, this is simplified by assuming the survey footprint is approximately rectangular in Cartesian space and assigning Cartesian equal dimensions to all the sky voxels in the map.However, as we showed, this assumption becomes increasingly inaccurate for wider sky surveys, already causing > 20% bias for our 20 × 20 deg 2 simulated survey.Experiments like MeerKLASS (Santos et al. 2017) are now gathering data for a > 4,000 deg 2 radio intensity mapping survey, thus currently adopted approximations will no longer be suitable.The highly adaptable pipeline we have presented in this work provides the alternative framework for accurate analysis with future LIM.
The techniques we have outlined are all compatible with the other data reduction techniques necessary in LIM analysis.For example, foreground cleaning is a major part of an unbiased LIM pipeline and one we have not discussed in this work.However, blind foreground removal techniques (Alonso et al. 2015;Carucci et al. 2020;Cunnington et al. 2021a;Spinelli et al. 2021) can be trivially performed on the original sky map data where the continuum foreground remains coherent along the frequency (line-of-sight) direction, for optimal component separation.Once foreground cleaning is performed, the regridding steps we have outlined in this work can then be performed with no modification necessary.Furthermore, the current adopted technique of using a transfer function to correct for signal loss in the foreground clean (Cunnington et al. 2023) can still be implemented.The only change needed here would involve emulating the same regridding step that was performed on the data, on the injected mocks used for the transfer function computation.Whilst we have not explicitly checked these assumptions, we see no obvious reason that our regridding and modelling techniques will cause unforeseen problems to LIM foreground cleaning.We mention this to motivate an explicit check with future work.
The field-level corrections we presented will be directly applicable to higher-order -point estimators such as the bispectrum (Randrianjanahary et al. 2023), and also one-point statistics (Bernal 2023).Furthermore, much of the modelling of discretistaion effects can be tweaked slightly to apply to high-order estimators, as similarly shown in Cunnington et al. (2021b) for modelling observational effects in the bispectrum.It would also be interesting to see if the higher-order assignment schemes deliver a clearer improvement in performance for these statistics which are more sensitive to non-linear and non-Gaussian fluctuations.
In addition to validating the claim that foreground removal presents no additional challenges for the regridding process (or vice-versa), we intend to further expand our investigation in future work.This work used simple lognormal mocks cut to a survey footprint.For a more robust test of the pipeline, we aim to extend these simulations to include a redshift evolving lightcone (as in Sato-Polito et al. 2022) and to also assign luminosities following a well-motivated luminosity function (as in Niemeyer et al. 2023).These more realistic simulations will also introduce non-linear effects, boosting power on small scales, thus potentially enhancing aliased power.Furthermore, as discussed in Appendix C, anisotropic observational effects such as a non-Gaussian, frequency-dependent beam will present additional challenges for LIM analysis, not just for regridding and aliasing mitigation techniques.We have also not included contributions to the field from thermal noise, residual foregrounds or systematics.It will also be necessary to include the case where LIM data are cross-correlated with galaxy surveys.These surveys should exhibit different aliasing contributions and hence validating that we can simultaneously handle their effects will be essential for future precision cosmology which will benefit from LIM-galaxy cross-correlations.
The Python code which performs the LIM simulations, regridding, higher-order interpolations, power spectra estimation and modelling discussed in this paper, is made available at github.com/stevecunnington/gridimp.Method (i) should perfectly emulate reality providing that a curved sky beam kernel can be formulated, which is no trivial task.Method (ii) is analytically simpler but would require more computational demand since the simulation pipeline will now require higher nside maps and the additional resampling to the final resolution.The idea behind the latter method is that the discretisation onto a high nside map will produce negligible discretisation so applying the beam here will not artificially suppress the effects we wish to investigate.Introducing an enhanced beam model beyond a Gaussian profile will also increase the realism of the simulation, as explored in previous work (e.g.Matshawule et al. 2021;Spinelli et al. 2021).We leave the pursuit of these challenges to future work.
With the caveat that our application of the beam is a somewhat simplified process, we found some interesting effects when analysing the distortions and aliasing in the anisotropic 2D ( ⊥ ,  ∥ ) power spectra.Figure C1 shows another version of Figure 5, comparing the model agreement between the Cubic voxels and Fine channels simulation, but now with the beam applied.We we see a similar conclusion that the discrepancies between measurement and model are pushed to higher  ∥ for the Fine channels simulation, beyond the lowest Nyquist frequency (red dashed line) determined from the perpendicular directions.However, at the higher  ∥ resolved by the grid, we now see lots of under-estimated power relative to the model (red regions) which we did not see in the resolved modes without the beam (Figure 5).The power spectrum modelling we use should include a perfect representation of the beam that we applied to the Cartesian field in the simulations, therefore it is difficult to explain the source of the negatively biased modes.The fact that this negative discrepancy is most pronounced for the high- ∥ modes in the Fine channels version is also interesting since the beam is damping high- ⊥ modes.Evidence was presented in the Villaescusa-Navarro et al. (2017) to show the damping effect from the beam can impact parallel scales.This was in the context of parallel baryon acoustic oscillations in the presence of a large beam, but it could be related to what we see.Further investigation with a more realistic beam is warranted in future work to fully understand this effect.For now, it is sufficient to see that when we include the full modelling treatments onto a more resolved Cartesian grid (see Figure C2, following steps in Section 3) we see excellent agreement for the area of -space ( ≲ 0.2 ℎ Mpc −1 ) of interest for the 21cm LIM survey we have simulated.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure1.Toy illustration of LIM pixelisation and Cartesian gridding processes necessary for analysis in Fourier-space.The background image represents a continuous fluctuation field for the spectral line (field is mean-centred and mean-normalised).LIM surveys will discretise their calibrated observations of the continuous field into broad sky pixels (blue dashed boundaries in insert).This provides the 3D voxel array in sky coordinates (R.A., Dec, ).Using a Monte-Carlo method with sampling particles, the sky voxel intensities can be regridded onto a 3D Cartesian grid (red boundaries) in comoving space ( x ,  y ,  z in ℎ −1 Mpc) which can be Fast Fourier transformed.The method we adopt (detailed in Section 2.1) creates pseudo-random particles by generating  p per input sky voxel ( p = 10 for this example), with one particle at the voxel centre and the other  p − 1 distributed uniformly randomly within the voxel.

Figure 2 .
Figure 2. Sampling coverage on the final Cartesian grid for a toy 20 × 20 deg simulated 21cm LIM survey spanning 900 <  <1038 MHz (0.37 <  < 0.58).The colour scale show the counts of sampling particles falling in the grid cells.Every panel shows a slice halfway through the array along the z,y and x directions going left to right, as indicated by the red directional arrows.Top-left text in each row gives details for each scenario. G is the number of cells on the side of the cubic Cartesian grid. p is the number of sampling particles per pixel of the input map being regridded.

Figure 3 .
Figure 3. Mapping and regridding steps for one realisation of our simulated sky surveys covering 10 • < R.A. <30 • and 10 • < Dec. <30 • .The red flowchart arrows and text indicate which steps in the simulation process (outlined near the beginning of Section 2.2) are being performed.Left-panel shows the starting simulation; a lognormal realisation generated from a power spectrum model.The size of the input field is determined so that the grid encloses the survey footprint, with some border padding to allow for non-local particle interpolation.Middle-panel shows the full map in sky coordinates produced by sampling  0 into HEALPix pixels and frequency channels with chosen resolutions nside and .The middle-panel inset zooms in on the survey footprint showing the 20 × 20 deg 2 field, also in sky coordinates.Right-panel shows the final Cartesian grid resampled from the sky voxels with a chosen interpolation scheme for the sampling particles (we used nearest grid point assignment for this demonstration).

Figure 4 .
Figure 4. Accuracy of the measured power spectrum of the regridded field relative to the model (black lines) for the two simulation versions (Cubic voxels and Fine channels).We show the case where no compensation for the NGP interpolation is applied to the regridded field (blue lines) and the case where compensation is applied using Equation 8 and Equation9(orange lines).We do not include the effects of the telescope beam in these results.Measured power spectra PG are averaged over 100 simulation realisations.The 1 scatter from the 100 simulations is shown by the shaded regions.

Figure 5 .
Figure 5. Cylindrical ( ⊥ ,  ∥ ) power spectrum comparison to model for the Cubic voxels (left panel) and Fine channel (right panel) simulations.Blue regions represent overly enhanced power relative to the true model driven by aliasing.The red dashed line shows the same minimum Nyquist frequency as in Figure 4 calculated from the cell resolution along one of the dimensions in the perpendicular directions.Colour-bar range is fixed to avoid saturation.

Figure 6 .
Figure 6.Accuracy of the measured power spectrum relative to the model for the case with and without smoothing effects from the telescope beam.Measured power spectra PG are averaged over 100 simulation realisations.

Figure 7 .
Figure 7. Percent accuracy of the measured power spectrum relative to increasingly sophisticated models ( ′ mod Equation6), ( ′′ mod Equation14) and ( ′′′ mod Equation15).Green shaded areas show the sub-5% and 1% accuracy regions.Measured power spectra PG are averaged over 100 simulation realisations.Blue line shows the limited model from Section 2. The orange line includes damping from the discretisation of the field into the map pixels and frequency channels.The red line further includes the summation over modes.

Figure 8 .
Figure 8. Window functions in 1D for the different mass assignment schemes considered in this work.Top panel shows the different interpolations in configuration space used to assign the particle intensities into cells as a function of the normalised distance from the cell .The bottom panel shows the Fourier transform of the window functions.

Figure 9 .
Figure 9. Percent accuracy of the measured power spectrum relative to full model  ′′′ mod (Equation15) for the higher-order mass assignment schemes, both with and without interlacing.The green shaded area shows the sub-1% accuracy region.Measured power spectra PG are averaged over 100 simulation realisations.

Figure C1 .
Figure C1.Cylindrical ( ⊥ ,  ∥ ) power spectrum comparison to model for two simulation versions.Same as Figure 5 but with smoothing effects from the telecope beam.The red dashed line shows the minimum Nyquist frequency calculated from the cell resolution along one of the dimensions in the perpendicular directions.Colour-bar range is fixed to avoid saturation.

Figure C2 .
Figure C2.Cylindrical ( ⊥ ,  ∥ ) power spectrum for the Cubic voxels simulation (with the beam) relative to the final full model  ′′′ mod (Equation15) including the treatments for discretisation and aliasing from the sky map voxels.This shows the 2D version of Figure9where PCS and interlacing has been used in the regridding of the data.The red dashed line shows the minimum Nyquist frequency calculated from the cell resolution along one of the dimensions in the perpendicular directions.Colour-bar range is fixed to avoid saturation.

Table 1 .
Specifications for simulated data.Both simulation versions have 100 different realisations which we average results over.We chose to emulate a generic Hi 21cm LIM survey centred at 1 GHz, corresponding to a central effective redshift of  eff = 0.435.A 3D grid size is calculated to fully enclose the survey footprint which defines the size of the input lognormal simulations.The final Cartesian grid has these same dimensions ( x ,  y ,  z ) but much broader cells ( 0x ,  0y ,  0z →  Gx ,  Gy ,  Gz ) initially chosen to approximately double the size of the HEALPix voxels.
* Note that in some cases, which we clearly identify, we do not include the beam, equivalently setting  FWHM = 0.

Table 2 .
Average percentage bias relative to model in results from Figure9for scales 0.07 <  < 0.25 ℎ Mpc −1 .For completeness, we include the NGP results not plotted in Figure9.Percentages are the average of the absolute bias i.e. 100 × | PG / mod − 1|, to avoid positive and negative biases cancelling.