title: Power spectrum multipole expansion for HI intensity mapping experiments: unbiased parameter estimation
authors: Paula S. Soares Steven Cunnington Alkistis Pourtsidou Chris Blake
text: 
  
  firstpage–lastpage
  : Semantic API Migrations
      Michael F.P. O'Boyle
      December 30, 2023
  =========================
  
  
  
  
  
  We assess the performance of the multipole expansion formalism in the case of single-dish intensity mapping, including instrumental and foreground removal effects. This formalism is used to provide MCMC forecasts for a range of and cosmological parameters, including redshift space distortions and the Alcock-Paczynski effect. We first determine the range of validity of our power spectrum modelling by fitting to simulation data, concentrating on the monopole, quadrupole, and hexadecapole contributions. 
  We then show that foreground subtraction effects can lead to severe biases in the determination of cosmological parameters, in particular the parameters relating to the transverse BAO rescaling, the growth rate and the bias (α_⊥, T_ fσ_8, and T_ b_σ_8, respectively). We attempt to account for these biases by constructing a 2-parameter foreground modelling prescription, and find that our prescription  leads to unbiased parameter estimation at the expense of increasing the estimated uncertainties on cosmological parameters. In addition, we confirm that instrumental and foreground removal effects significantly impact the theoretical covariance matrix, and cause the covariance between different multipoles to become non-negligible. Finally, we show the effect of including higher-order multipoles in our analysis, and how these can be used to investigate the presence of instrumental and systematic effects in intensity mapping data.
  
  
  
  
  
  cosmology: large-scale structure of Universe – cosmology: theory – cosmology: observations – radio lines: general
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  The standard cosmological model, ΛCDM, describes a universe with zero spatial curvature, containing cold dark matter and dark energy in the form of a cosmological constant (Λ), which drives late time cosmic acceleration. It has 6 parameters, 5 of which have been measured to within 1% precision through observations of the Cosmic Microwave Background <cit.>. Large Scale Strucrure (LSS) surveys have also provided observations that are in very good agreement with ΛCDM <cit.>. LSS surveys in particular are able to probe whether general relativity is the correct description of gravity on cosmological scales by measuring the logarithmic growth rate of structure (f)  <cit.>. This parameter can be measured through the redshift space distortion (RSD) signature on the 2-point statistics of galaxy clustering <cit.>.
  
  Neutral hydrogen () Intensity Mapping (IM) is a novel technique that is able to efficiently and rapidly observe a very wide redshift range, including high redshifts, z>3, that are inaccessible by current and forthcoming optical galaxy surveys (see <cit.> for a review). In particular, IM treats the 21cm sky as a diffuse background and measures its intensity in large voxels, as opposed to detecting individual galaxies <cit.>. In the post-reionization Universe, neutral hydrogen resides inside galaxies where it is self-shielded from ionization; it can thus be used as a tracer of the underlying matter distribution.  Using IM, it is possible to map the 3D LSS of the Universe, and probe the underlying cosmology through the power spectrum.
  
  Current IM detections come from cross-correlating IM maps from the Green Bank Telescope (GBT) or the Parkes radio telescope with optical galaxy surveys, probing the clustering of neutral hydrogen at z<1 <cit.>. More specifically, the GBT has constrained the combination of the abundance (Ω_) and linear bias (b_) at z=0.8, Ω_ b_ r = [4.3 ± 1.1] × 10^-4, using cross-correlation with the WiggleZ optical galaxy survey, where r is the galaxy-hydrogen correlation coefficient <cit.>. A detection through auto-correlation is yet to be made due to residual systematic effects. However, since most of these systematic effects do not correlate with  optical galaxy surveys, they are mitigated in cross-correlation. 
  
  The Square Kilometre Array (SKA)[www.skatelescope.orgwww.skatelescope.org] will be a radio observatory able to reach unprecedented statistical precision on IM measurements (see e.g. <cit.>,  <cit.>), assuming that systematic effects are controlled or mitigated. In the case of IM, large galactic and extragalactic foregrounds dominate over the signal by several orders ot magnitude. However, in principle we can differentiate these dominant foregrounds from the signal since they are expected to be smooth in frequency <cit.>.
  
  In addition, it is not yet fully understood how the systematic effects present in an IM survey affect its noise properties, in particular the covariance matrix. The effect of foreground removal has been shown to mostly affect power on large scales. The telescope beam significantly damps power on scales smaller than its resolution, but it also affects larger scales <cit.>.  Analytically, both of these effects carry over to the theoretical covariance matrix in the case of the IM power spectrum <cit.>.
  
  In this paper, we build upon the work of <cit.> (hereafter C20) and <cit.>. <cit.> studied the modelling of the IM power spectrum including observational effects, and C20 extended this into a comprehensive simulations and data analysis pipeline[https://github.com/IntensityTools/MultipoleExpansiongithub.com/IntensityTools/MultipoleExpansion] for analysing the IM power spectrum multipoles taking into account instrumental and foreground removal effects. Here we extend on C20 to perform cosmological parameter estimation with the IM power spectrum, using simulations that include the relevant instrumental and foreground removal effects. In particular, we use Markov Chain Monte Carlo (MCMC) analyses to forecast uncertainties for a range of and cosmological parameters. We aim to realistically assess what a future SKA-like IM survey will be able to constrain, and in particular how foreground removal affects cosmological parameter estimation. We are interested in both precision and accuracy, i.e. we pay particular attention to the requirement of unbiased parameter estimation.
  
  The paper is structured as follows: In sec:model, we describe the observed IM power spectrum, including modelling of the telescope beam and foreground removal, decomposed into multipoles. In sec:sim, we describe our IM simulations. In sec:results, we test the range of validity of our model using the simulations, report our MCMC analysis results, look into how instrumental and systematic effects affect the noise covariance matrix, and investigate whether higher order multipoles can add useful information. We conclude in sec:conclusion. 
  
  Throughout this paper, we assume a flat ΛCDM cosmology consistent with the Planck15 analysis <cit.>, with Ω_M = 0.307, Ω_b = 0.048, Ω_Λ = 0.693,  σ_8 = 0.823, n_s = 0.96 and Hubble parameter h=0.678.
  
  
  
  § MODEL
  
  
  
  
   §.§ Redshift space distortions
  
  
  Redshift space distortions (RSD) introduce anisotropies in the observed power spectrum. In order to account for this, we consider the power spectrum as a function of the directional wave vector k⃗, which can be decomposed into its module k and the cosine of the angle θ between the wave vector and the line-of-sight (LoS) component μ = k⃗· k_∥≡cosθ. We model RSD by considering the Kaiser effect <cit.>, which is a large-scale effect dependent on the growth rate f, and the Fingers-of-God (FoG) effect <cit.>, which is a small-scale non-linear effect that depends on the velocity dispersion of the tracer objects (σ_v). The anisotropic power spectrum can be written as: 
  
      P_(k, μ) =  ( T_ b_ + T_ f μ^2)^2 P_M(k)/1+(k μσ_v / H_0)^2 + P_SN  ,
  
  
  where P_SN = T_^2 (1/n) is the shot noise, n is the number density of objects, P_M(k) is the underlying matter power spectrum, and T_ is the mean brightness temperature, modelled as <cit.>:
  
  
      T_(z) = 180Ω_(z)h(1+z)^2/H(z)/H_0 mK  .
  
  
  
  
   §.§ Alcock-Paczynski effect
  
  
  When we measure the power spectrum using intensity mapping, we first measure redshifts and then transform these into distances. In order to do this transformation, we must assume a cosmology. If the assumed cosmology does not match the real one, we get further anisotropies in the power spectrum measurements. This is known as the Alcock-Paczynski (AP) effect <cit.>. In the transverse and radial directions respectively, we model these anisotropies as (see e.g. <cit.>):
  
      α_⊥ = D_A(z)/D_A(z)^ f  ,
          α_∥ = H(z)^ f/H(z)  ,
  
  
  where throughout the paper the superscript or subscript `f' refers to the fiducial value, in this case our fiducial, assumed cosmology. Here we note that we follow the notation and method of <cit.> and do not include the degeneracy with the sound horizon at radiation drag (r_s) in these factors, as is done when performing a BAO-only analysis. This is because we are performing a full shape analysis. These factors distort the perpendicular and parallel to the LoS wave vectors as: 
  
      k_⊥ = k_⊥^ f/α_⊥  , 
           k_∥ = k_∥^ f/α_∥  .
  
  
  It is useful to define the factor F_ AP = α_∥ / α_⊥, which helps describe how k and μ become distorted, and how to recover the true underlying value from the fiducial value:
  
  
      k = k^ f/α_⊥[ 1 + (μ^ f)^2 (F_ AP^-2 - 1) ]^1/2  ,
          μ = μ^ f/F_ AP[ 1 + (μ^ f)^2 (F_ AP^-2 - 1) ]^-1/2  .
  
  
  The power spectrum can then be described in terms of this effect as:
  
  
      P_(k^ f, μ^ f) = α_∥^-1α_⊥^-2[ ( T_ b_ + T_ f μ^2)^2 P_M(k)/1+(k μσ_v / H_0)^2 + P_SN]   .
  
  
  
  
   §.§ Telescope beam smoothing effect
  
  
  The telescope beam introduces one of the main instrumental effects in the case of single-dish intensity mapping experiments. We can model this effect using a damping term dependent on the physical smoothing scale of the beam (see e.g. <cit.>).  Assuming the telescope beam can be modelled as a Gaussian, this is defined as R_beam = σ_θ r(z), where σ_θ =  θ_FWHM / ( 2√(2ln(2)) ), θ_FWHM is the full-width-half-maximum of the beam in radians, and r(z) is the comoving distance to a redshift z. The Fourier transform of the telescope beam damping term is:
  
  
      B_⊥(k,μ) = exp(-k^2 R_beam^2(1-μ^2)/2)   ,
  
  
  and the power spectrum becomes:
  
  
      P_(k^ f, μ^ f) =     B_⊥^2(k,μ)/α_∥ α_⊥^2
          ×[ ( T_ b_ + T_ f μ^2)^2 P_M(k)/1+(k μσ_v / H_0)^2 + P_SN]   .
  
  
  
  We should also note that for surveys that are limited in frequency resolution, a similar effect will occur on the small radial scales. In our case, and in general for single dish experiments, the frequency resolution is very good (much better than the angular resolution and radial non-linear dispersion effects) and does not cause any discernible effects in our power spectrum measurements, so we do not include it in our modelling. In cases where this might be relevant, a way to account for this smoothing effect is described in <cit.>. Given the resolution in the radial direction is set by the frequency channel bandwidth δ_ν, the smoothing effect can be modelled as: 
  
      B̃_∥(k,μ) = 
          sin(kμ s_∥/2)/kμ s_∥/2  ,
  
  where s_∥ = [c/H(z)](1+z)^2 (δ_ν/ν_21), with ν_21 being the rest emission frequency.
  
  
  
  
   §.§ Re-normalising by σ_8
  
  
  In this work, we calculate the underlying non-linear matter power spectrum given a fiducial, assumed cosmology using the  package [https://nbodykit.readthedocs.io/en/latest/https://nbodykit.readthedocs.io] <cit.>, which uses the CLASS Boltzmann solver <cit.>, and we choose the Halofit prescription <cit.>. It is useful to parametrise this template calculated matter power spectrum P_M(k) by σ_8, which is the RMS of the density fluctuations within a sphere of radius 8 h^-1Mpc (see e.g. <cit.> for a more detailed description):
  
  
      P_M,8(k) = P_M(k)/σ_8^2  .
  
  
  Including this, our final power spectrum model becomes:
  
  
      P_(k^ f, μ^ f) =     B_⊥^2(k,μ)/α_∥ α_⊥^2
          ×[ ( T_ b_σ_8 + T_ f σ_8 μ^2)^2 P_M,8(k)/1+(k μσ_v / H_0)^2 + P_SN]   .
  
  
  The set of parameter combinations that can be measured using this model is:
  
  
      θ⃗ = {α_∥,  α_⊥,  T_ fσ_8,  T_ b_σ_8,  σ_v,   P_SN}  .
  
  
  We note that, in comparison to optical galaxy surveys, we have an additional degeneracy coming from the mean brightness temperature T_, which is proportional to Ω_. Previous works using Fisher matrix forecasts (see e.g. <cit.>) assume this is a known quantity and keep it fixed, but here we choose to include it since Ω_ is quite poorly constrained <cit.>. We also note that, as suggested in <cit.>, this degeneracy can be broken by using information from the non-linear regime of structure formation and perturbation theory modelling, but this would require precise and well-calibrated interferometric observations. In this work we are assuming a survey in single-dish mode <cit.> and we are focusing on the beam and foreground removal effects (with the latter being a major issue for both single dishes and interferometers).
  
  
  
   §.§ Multipole expansion
  
  
  We can expand the anisotropic power spectrum P_(k, μ) in terms of Legendre polynomials as
  
  
      P_(k, μ)=∑_ℓ P_ℓ(k) ℒ_ℓ(μ)   ,
  
  
  where ℒ_ℓ(μ) is the ℓ^th Legendre polynomial:
  
      ℒ_0 = 1 ,   ℒ_2 = 3μ^2 -1/2 ,   ℒ_4 = 35μ^4 -30μ^2 + 3/8 , 
          ℒ_6 = 231μ^6 - 315μ^4 + 105μ^2 - 5/16 .
  
  
  Our full model, expanded into power spectrum multipoles, is then given by:
  
  
      P_ℓ(k^ f)=     2 ℓ+1/2 (α_∥^-1 α_⊥^-2) ∫_-1^1μ^ f ℒ_ℓ(μ^ f)  B_⊥^2(k,μ)  
          ×[ ( T_ b_σ_8 + T_ f σ_8 μ^2)^2 P_M,8(k)/1+(k μσ_v / H_0)^2 + P_SN]   .
  
  
  We consider the monopole (P_0), quadrupole (P_2), hexadecapole (P_4) and 64-pole (P_6) in our analysis. 
  
  
  
   §.§ Modelling the effect of foreground removal
  
  
  In order to model the effect of foreground removal on the power spectrum, we introduce a damping term inspired by the survey volume damping function (see e.g. <cit.>). This is given in Fourier space by:
  
  
      B_ vol(k, μ) =     ( 1 - exp{ -( k_⊥/k^ min_⊥)^2 }) 
          ×( 1 - exp{ -( k_∥/k^ min_∥)^2 })  ,
  
  
  which describes how we are not able to access modes smaller than k^ min_⊥ or k^ min_∥ in the perpendicular and parallel to the LoS directions. If we assume a survey box to have comoving distance dimensions given by [L_x, L_y, L_z], we have that the smallest (largest) modes (physical scales) accessible in the perpendicular and parallel to the LoS directions are: k^ min_⊥ = 2π/√(L_x^2 + L_y^2) and k^ min_∥ = 2π / L_z.
  
  We assume that the process of foreground removal similarly removes power from modes along the parallel and perpendicular to the LoS directions based on the foreground properties and survey geometry. In this case, we are particularly considering the effects of an Independent Component Analysis (ICA) foreground removal technique (see e.g. <cit.>). This component separation technique does not try to assume a specific form for the foreground contamination, but relies on the fact that the sources of the foregrounds are statistically independent and can be isolated from the cosmological signal. The foregrounds are mostly smooth in frequency (except in the presence of effects such as polarization leakage), while the cosmological signal is not smooth in frequency, since it traces the structure of matter in the Universe. This allows component separation techniques to separate the foregrounds from the desired underlying cosmological signal, and remove them. However, these techniques can confuse the signal with the foreground, usually in the largest scale limits of the particular box, where the signal looks smooth. This leads to a loss of signal, which we try to model using a damping function across radial and transverse modes.
  
  In this work we introduce a two-parameter damping prescription to model the effects of foreground removal. The two parameters, N_⊥ and N_∥, vary the scale of modes being damped by foreground removal. If N_⊥ and N_∥ equal zero, this corresponds to no damping. We expect the combined damping factors N_∥ k^ min_∥ to be greater than N_⊥ k^ min_⊥ since foreground removal mainly removes signal along the radial (LoS) direction. 
  
  Although we could have quoted the N_∥ k^ min_∥ and N_⊥ k^ min_⊥ values together as just two parameters, e.g. k^ FG_∥ and k^ FG_⊥, we choose instead to keep this form where we have the N_⊥ and N_∥ parameters present. The main reason for this is that N_⊥ and N_∥ are independent of the box dimensions (universal for similar conditions), while N_∥ k^ min_∥ and N_⊥ k^ min_⊥ depend on k^ min_∥ and k^ min_⊥ (specific to our case, will vary depending on the particular geometry of different simulations or surveys). We expect a user to retrieve similar N_⊥ and N_∥ to us if using a similar foreground removal method (e.g.  with N_ IC=4, <cit.>) regardless of the box dimensions, and thus find this the most relevant parameter to quote. 
  
  However, we note that we would expect N_⊥ and N_∥ to be larger for more aggressive foreground removal methods (with higher N_ IC), which are employed in real data to deal with more complicated foregrounds, noise and systematics (see e.g. <cit.>). In particular, real foregrounds might experience polarization leakage, an effect which hinders their spectral smoothness, and would require a more aggressive N_ IC choice to fully remove (see e.g. <cit.> for more insight on the effect of polarized foregrounds). In addition, the model is only valid for cases where the survey geometry is constant (i.e., L_x, L_y, L_z are not changing with redshift). Further consideration would be needed in order to successfully apply this model to a lightcone with realistic survey geometry.
  
  The damping term for modelling the effects of foreground removal is given in Fourier space by:
  
  
      B_ FG(k, μ) =     ( 1 - exp{ -( k/N_⊥ k^ min_⊥)^2 ( 1 - μ^2 ) }) 
          ×( 1 - exp{ -( k/N_∥ k^ min_∥)^2 μ^2 })   ,
  
  
  and the power spectrum model in the presence of foreground removal effects is:
  
  
      P_(k^ f, μ^ f) =     B_⊥^2(k,μ) B_ FG (k,μ)/α_∥ α_⊥^2
          ×[ ( T_ b_σ_8 + T_ f σ_8 μ^2)^2 P_M,8(k)/1+(k μσ_v / H_0)^2 + P_SN]  .
  
  
  When applying the multipole expansion formalism to this model, we obtain:
  
  
      P_ℓ(k^ f)=     2 ℓ+1/2 (α_∥^-1 α_⊥^-2) ∫_-1^1μ^ f ℒ_ℓ(μ^ f) B_⊥^2(k,μ) B_ FG(k,μ) 
          ×[ ( T_ b_σ_8 + T_ f σ_8 μ^2)^2 P_M,8(k)/1+(k μσ_v / H_0)^2 + P_SN]   .
  
  
  
  
  § SIMULATIONS
  
  
  
  The simulated data we use in this investigation are the same as in C20 and we refer the reader there for a more in depth introduction. For completeness, we provide a summary here of the cosmological signal simulations (sec:CosmoSig) and the foreground simulations (sec:ForegroundSims).
  
  
  
   §.§ Cosmological signal
  
  
  The source of our simulated cosmological data is from the MultiDark-Galaxies data <cit.> and the catalogue produced from the SAGE <cit.> semi-analytical model application. These galaxies were produced from the dark matter cosmological simulation MultiDark-Planck (MDPL2) <cit.>, which follows the evolution of 3840^3 particles in a cubical volume of 1  Gpc^3 h^-1 with mass resolution of 1.51×  10^9h^-1M_⊙ per dark matter particle. The cosmology adopted for this simulation is based on Planck15 cosmological parameters <cit.>, with {Ω_M,Ω_b, Ω_Λ, σ_8, n_s, h} = {0.307, 0.048, 0.693, 0.823, 0.96, 0.678}.
  
  As in our previous work (C20), we use the data from the z=0.82 redshift snapshot. At this redshift the box size with comoving distance dimensions L_x = L_y = L_z = 1000 Mpc h^-1 approximately corresponds to a sky area of 29×29 deg^2 with a redshift depth of Δ z = 0.5. Using Nearest Grid Point (NGP) assignment, we bin the catalogue of galaxies onto a grid with voxel dimensions N_x=N_y=N_z=225. We checked that using a higher resolution grid with N_x=N_y=N_z=512 made no discernible difference in our analysis. From the survey volume, we have that k_ min = 2π/V^1/3 = 0.006  h Mpc^-1, and use bins of width Δ k = 0.013  h Mpc^-1 to avoid correlations between bins.
  
  The SAGE catalogue we use has cold gas mass outputs for each galaxy from which we can compute a brightness temperature in each pixel of our map. However, since the simulation has a finite mass resolution, the lowest mass halos (≲ 10^10 h^-1M_⊙) which also contain will not be properly sampled (see C20 for further discussion). This is an important limitation to consider, since it affects the observed brightness, bias, and probability distribution of in our simulation. But the lowest mass haloes, albeit more abundant, contribute the least to the total brightness, and their exclusion has little effect on the bias <cit.>. To ensure a realistic global signal is present in the data we rescale the mean temperature such that the abundance is consistent with a value obtained in real data analyses at this redshift,  Ω_∼ 4.3 × 10^-4 <cit.>. This provides our simulated data with a realistically distributed signal with a mean value of T_ = 0.13 mK.
  
  
  
  We aim to emulate an upcoming SKA1-MID-like experiment <cit.> and we therefore include simulated instrumental effects from the radio telescope beam. Using the diameter of an SKA dish (D_max = 15 m) we can calculate the beam size of such an experiment from
  
  
      θ_FWHM = 1.22  λ_21/ D_max (1+z)   ,
  
  
  which for our case yields a beam of θ_FWHM=1.78deg, equivalent to R_beam = σ_θ r(z) = 26.16 Mpc h^-1 (see sec:beammodel). After combining our simulated temperature fluctuation field and foregrounds, we convolve the final simulation cube with the telescope beam described above.
  
  Using the  simulation without RSD or systematic effects, we measure the linear bias of the simulation to be  b_ (k) = √(P_ (k) / P_M,linear(k)), and by averaging this quantity over the large, linear k-scales we obtain b_ = 1.16 ± 0.04. We roughly measure the upper limit on the shot noise to be P_ SN = 2.5  mK^2  Mpc^3  h^-3. This is done by looking at our measurement of P_ (k) in a high resolution grid simulation of N_side=512, and seeing where the power spectrum `plateaus' at high k. At these scales the shot noise should be dominating over the faint cosmological signal, and this gives us an idea of what the upper limit of the shot noise is.
  
  
  
   §.§ Simulating the effect of foregrounds
  
  
  In order to generate realistic foregrounds we utilize the Global Sky Model (GSM) <cit.>, which extrapolates maps from real data at the desired frequency. For the redshift depth of our simulated data we can assume a frequency range of 673<ν<903 MHz and we therefore generate N_z=225 maps spanning this range. Therefore, the simulated foregrounds have a realistic evolving spectral index which is still smooth relative to the cosmological signal and will allow for successful component separation in the foreground clean. 
  
  In reality, foreground signals are likely to be more complex and include contributions from free-free emission, extra-galactic point-sources and suffer effects from polarization leakage. This often requires a more aggressive foreground clean than what is typically required on a simulation from the GSM alone. In order to add additional complexity to the simulated foregrounds, we also generate realizations of diffuse emission from a model power spectrum which aims to describe different foreground sources. We follow the details outlined in Table 1 of C20 to produce these contributions, which include models of extragalactic point sources and free-free emission. We then combine these realizations with the GSM outputs to complete the full sky foreground data.
  
  We then need to transform these full-sky maps into flat-sky data with the same dimensions as our cosmological simulation. To do this we define an angular coordinate for each pixel on the flat-sky map, which we match to a pixel in the [http://healpix.sourceforge.nethttps://healpix.sourceforge.io/] <cit.> map with the closest angular coordinate. While this approach is an approximation and may affect some angular coherence in the foreground maps, it will have no impact on the foreground as a contaminant to our data. We add these flattened foreground maps onto the cosmological maps to contaminate them and create the requirement for a foreground clean. 
  
  For the foreground cleaning, we use the blind foreground removal method Fast Independent Component Analysis () <cit.>, and refer the reader to <cit.> for a more detailed description. As discussed when introducing the foreground modelling in sec:FGmodel and in C20, this method removes the foregrounds by assuming that the raw, uncleaned data can be written as a linear equation, where the signal can be broken up into statistically independent components:
  
  
      x = As + ε = ∑_i=1 ^N_IC=ma_i s_i + ε  ,
  
  
  where m describes the number of independent components, x is the raw, uncleaned data, A is the mixing matrix which describes the amplitude of the independent components, s are the m independent components, and ε is the residual which includes noise and the cosmological signal. As an input, we choose N_IC = 4 in accordance with previous studies <cit.>. The independent components in this case are the foregrounds, and by appropriately identifying and removing these we are left with ε, which contains our cosmological signal.
  
  This process of foreground removal is imperfect, and tends to confuse the signal with the foreground at large scales, especially in small k_∥ modes where the cosmological signal also appears smooth in frequency. This leads to cosmological signal being removed, which affects the amplitude of the power spectrum (see e.g. <cit.> for further discussion). We show in our analysis that it is possible to account for this effect using a model with free parameters that we let vary.
  
  
  
   §.§ Instrumental Noise
  
  
  Instrumental noise is determined by the telescope configuration. For an SKA-like single-dish experiment, we assume the pixel noise is well represented by a Gaussian random field with spread given by:
  
  
      σ_pix = T_sys √(4π  f_sky/Ω_beam N_dish t_obs ν)  ,
  
  
  from which the noise power spectrum is then given by P_N = σ_pix^2 V_pix, where V_pix is the voxel volume given by:
  
  
      V_pix = Ω_beam∫^z + Δ z / 2_z - Δ z / 2 dz dV/dzdΩ  ,
  
  
  where Ω_ beam is defined on <ref>, and
  
  
      dV/dzdΩ = cr(z)^2/H(z)  .
  
  
  We assume SKA1-MID-like parameters for the noise (see <cit.>), shown in <ref>. The calculated noise power spectrum from these specifications is P_N = 4   mK^2 Mpc^3 h^-3. 
  
  
  
  We plot a comparison between our  z=0.82 cosmological signal in the absence of foregrounds, the foreground cleaned signal, the estimated shot noise and the instrumental noise, as seen in <ref>. We see that the instrumental noise dominates over the signal for k > 0.25 h Mpc^-1. We also note that pathfinder surveys for the SKA will have higher noise levels, but in this work we focus on the prospects of using IM for precision cosmology, hence why we choose to use SKA1-MID-like specifications.
  
  
  
  Combining our model, the  simulation power spectrum results, and the noise specifications, we can see how well our model agrees with the  data in <ref>. We also plot our foreground cleaned data (N_IC = 4) against the foreground model in <ref>, using guesses for the parameters N_⊥ and N_∥ found by eye (N_⊥=2, N_∥=2). For the fiducial model, we choose to use the estimated b_ = 1.16 and P_ SN = 2.5  mK^2 Mpc^3 h^-3, and we guess by eye the velocity dispersion parameter to be σ_v = 200  km/s. Note that when performing an MCMC analysis and checking for biased parameter results, we will not try to recover the `fiducial' values of the shot noise or velocity dispersion (since these are only rough estimates), but we will try to recover the estimated bias (b_) as discussed at the end of sec:CosmoSig, which is degenerate with T_ and σ_8. The full set of fiducial cosmological parameters is outlined on <ref>. These are used in the fiducial model and also in the covariance matrix calculations.
  
  
  
  
  
  
  
  
  
   §.§ Covariance matrix
  
  
  There are three main sources of error arising from the considered  IM experiment: sample variance, instrumental noise and shot noise. The covariance per k and μ bin (neglecting mode coupling), is:
  
  
      σ^2 (k, μ) = ( P_(k, μ) + P_N )^2/N_ modes(k, μ)  ,
  
  
  where N_ modes(k, μ) is the number of modes in each k and μ bin with widths Δ k and Δμ, respectively:
  
  
      N_ modes(k, μ) = k^2 Δ k Δμ/8π^2V_ sur  ,
  
  
  where V_ sur is the volume of the survey. For a survey scanning a sky area of Ω_tot this is given by:
  
  
      V_sur = Ω_tot∫^z_ max_z_ min dz dV/dzdΩ  .
  
  
  By neglecting mode coupling, we are assuming that the different k-bins are uncorrelated. We have assessed this assumption by comparing the multipole power spectrum errors obtained using a jackknife test to those obtained theoretically, assuming a diagonal covariance matrix, and found them to be consistent, meaning the diagonal covariance matrix is sufficient for our purposes (this assumption was also tested in C20 in the same way, and also found to be sufficient). We also note that we do not have a window function introducing correlations between different modes.
  
  The covariance matrix of the power spectrum multipoles is comprised of the sub-covariance matrices of each multipole, and those between different multipoles (i.e. the matrix is not diagonal, as it is essential to model the non-zero covariance between multipoles, see CovSection). The sub-covariance matrix for power spectrum multipoles ℓ and ℓ' is <cit.>:
  
  
      C_ℓℓ^'(k) = (2ℓ + 1)(2ℓ^' + 1)/2∫^1_-1 dμ σ^2 (k, μ) ℒ_ℓ(μ) ℒ_ℓ^'(μ)   .
  
  
  It follows from this that the total error on each multipole is given by <cit.>:
  
  
      σ_P_ℓ(k)      = √((2ℓ + 1)^2/2∫^1_-1 dμ σ^2 (k, μ) ℒ_ℓ^2(μ))
           = (2ℓ + 1)√(∫^1_0 dμ( P_(k, μ) + P_N )^2 ℒ_ℓ^2(μ)/N_ modes(k, μ))  .
  
  
  We can compute [S/N](k), the total signal-to-noise ratio per k bin, as:
  
  
      [ S/N](k) ^2  = Θ⃗^T (k) C^-1(k) Θ⃗ (k)   ,
  
  
  where Θ⃗(k) is a vector describing the power spectrum per k bin: 
  
  
      Θ⃗(k) = [P_0 (k), P_2 (k), P_4 (k)]
  
  
  Similarly, the log-likelihood is proportional to the χ^2 statistic:
  
  
      logℒ∝ -1/2χ^2  , 
          χ^2 = ΔΘ⃗^T C^-1ΔΘ⃗ ,
  
  
  where ΔΘ⃗ is the difference between our model prediction and the measurement from our simulation for all multipoles and k bins. For example, if we were considering N_ℓ = 3 multipoles and N_k = 20 k-bins in each multipole, our covariance matrix would have dimensions N_ℓ N_k × N_ℓ N_k = 60 × 60 and the vector ΔΘ⃗ would have length N_ℓ N_k = 60.
  
  We calculate the theoretical [S/N](k) using our model (<ref>) and the fiducial parameter values from the simulation (<ref>), as well as the assumed instrumental noise P_N = 4   mK^2 Mpc^3 h^-3. We plot the result for each combination of multipoles in <ref>. As expected, including higher order multipoles yields the best [S/N](k) results, with the quadrupole adding most of the additional information at both low and high k compared to the case where only the monopole is considered. We can see that the hexadecapole and 64-pole also add information, especially around k ∼ 0.15  h Mpc^-1.
  
  
  
  
  
  § RESULTS
  
  
  
  
   §.§ Model validation
  
  
  In this section, we first aim to test our model's range of validity, i.e. for which k range can we trust our model to return unbiased results for the fiducial cosmological parameters of our simulation? The parameters we vary are {α_∥, α_⊥, T_ fσ_8,T_ b_σ_8, σ_v, P_ SN}. We know the fiducial values of the following parameters: {α_∥, α_⊥, T_ fσ_8,T_ b_σ_8}, these are outlined in <ref>. For all of our MCMC analyses, we keep cosmological parameters in the covariance matrix fixed to the fiducial values. We impose an upper limit of P_ SN = 6  mK^2  Mpc^3  h^-3 and σ_v = 600  km/s on the shot noise and velocity dispersion parameter priors. All other priors are flat positivity priors. The MCMC analysis is performed using the publicly available  package [https://emcee.readthedocs.io/en/latest/https://emcee.readthedocs.io] <cit.>. We vary our model (<ref>) in the log-likelihood (<ref>) using 500 walkers and 2000 samples.
  
  In order to validate our model, we run an MCMC analysis for different k_ max limits, stopping when we see that our MCMC results are biased outside of 2σ. This test is meant to check that we are not going to k values that are too large, where our theory modelling breaks down and yields biased parameter estimates. We do this only up to k = 0.24  h Mpc^-1, as beyond this we find that the signal to noise per k bin drops to below 15 (see <ref>). In addition, we know that our beam starts to dominate at k_ beam = π/R_ beam = 0.12  h Mpc^-1, so we only consider the range up to 2k_ beam = 0.24  h Mpc^-1, beyond which we assume the beam entirely dominates over the cosmological signal. We tested going beyond this k_ max limit, but found no improvement on parameter uncertainties as expected.
  
  For the case of the galaxy power spectrum, it has been shown that the monopole and quadrupole contain most of the cosmological information <cit.>. The hexadecapole contains additional information, but spectroscopic galaxy surveys have found that it needs to be considered to a smaller k_ max than the monopole and quadrupole due to non-linear effects (see e.g. <cit.>, <cit.>). Following these studies, we test whether this is also the case in IM. First, we determine at which point the parameter estimation results from the MCMC become biased outside of 2σ for the monopole and quadrupole only and find that parameter estimates are unbiased up to k_ max = 0.24  h Mpc^-1. Then, using the k_ max determined for the monopole and quadrupole, we start adding the hexadecapole with different k_ max limits and check when results become biased, finding that we can only include it up to a restricted range of k_ max = 0.09  h Mpc^-1.
  
  We summarise our MCMC analysis results within the model's determined range of validity in <ref>, where we quote the parameter uncertainties (1σ) with and without the restricted hexadecapole. We show the MCMC results in <ref>. It is clear from these results that adding the hexadecapole in a restricted range decreases the error margins on the parameters. We plot our power spectrum model calculated with the best fit values from the MCMC analyses in <ref>.
  
  Our findings are consistent with the literature. <cit.> studied the precision of a generic line intensity mapping experiment, using a nearly identical model with the one here and synthetic data, and found that including the hexadecapole improves the precision of BAO scale measurements by 10-60% in a Fisher matrix analysis. In the case of spectroscopic optical galaxy clustering, including the hexadecapole in an MCMC analysis has proved beneficial in decreasing uncertainties, also requiring a restricted k-range for the hexadecapole in order to obtain unbiased results (see e.g. <cit.>). 
  
  With regard to the BAO scale parameters, our findings are qualitatively similar to what has been found in <cit.>. This study shows that while the large SKA telescope beam smears out the isotropic BAO peak signature at z > 1, it is possible to use the radial 21cm power spectrum to measure H(z) from the BAO peak at percent level precision for redshifts z < 2.5. Our measurement of the radial AP effect parameter, which is a function of H(z), instead relies on the overall shape of the 21cm power spectrum, which is susceptible to systematic effects. Nonetheless, we considered systematic effects due to the beam and foreground removal in our analysis, and found similar, unbiased sub-10% percent level constraints on the radial AP parameter related to the expansion rate. In addition, even in the presence of a large beam, we also find sub-1% percent level constraints on the transverse AP parameter related to the angular distance.
  
  
  
  
  
  
  
   §.§ Covariance Matrix
  
  
  Here we discuss our choice of using a non-diagonal covariance matrix (<ref>), meaning we consider covariance between different multipoles. We first describe the details of this non-diagonal covariance matrix. 
  
  The full covariance matrix is symmetric, and is composed of each diagonal sub-covariance matrix (diagonal since we neglect mode coupling, and since we do not have a survey window function introducing further correlations between modes). Each sub-covariance matrix has dimensions N_k × N_k where N_k is the number of k bins considered per multipole, so overall the full covariance matrix has dimensions N_ℓ N_k × N_ℓ N_k (where N_ℓ is the number of multipoles being considered). That is:
  
  
      C = 
          [ C_00 C_02 C_04;      C_22 C_24;           C_44 ]
          ,
  
  
  where each sub-covariance matrix is:
  
  
      C_ℓℓ^' = 
          [ C_ℓℓ^'(k_0)           0           0;           0 C_ℓℓ^'(k_1)           0;           0           0           ⋱ ].
  
  
  The full diagonal covariance matrix, in the case where we do not consider the covariance between different multipoles, also has dimensions N_ℓ N_k × N_ℓ N_k and is given by:
  
  
      C_ diag = 
          [ C_00    0    0;    0 C_22    0;    0    0 C_44 ]
          ,
  
  
  where each sub-covariance matrix is also given by <ref>, but only for the case of ℓ = ℓ^'. 
  
  For more detailed discussion of the covariance of galaxy power spectrum multipoles under the Gaussian assumption, and in particular the significance of the covariance between different multipoles, see e.g. <cit.> and <cit.>.
  
  
  
    §.§.§ Effect of the telescope beam
  
  
  In order to compare the diagonal and non-diagonal cases, we calculate the [S/N] per k bin for each case with and without a telescope beam damping term. Results for ℓ = 0,2 are given in <ref>. For the case where the telescope beam damping term is present, we can see that using a non-diagonal covariance matrix makes a difference at both low and high k. At low k, including covariance between multipoles seems to increase the [S/N], while for higher k it decreases it. For the case where we do not include the telescope beam, we can see that the [S/N] per k bin does not differ significantly between including or excluding off-diagonal terms.
  
  It is interesting to also compare how the telescope beam changes correlations between different multipoles. We can calculate the correlation matrix from the covariance matrix as:
  
  
      Corr_ℓℓ^'(k) = C_ℓℓ^'(k)/√(C_ℓℓ(k) C_ℓ^'ℓ^'(k))  .
  
  
  We plot the correlation matrix in the case of no telescope beam and compare it to the case of including a telescope beam with R_ beam = 26.2   Mpc h^-1 (the same beam used in our simulation) in <ref> (top-row). We can clearly see that the presence of the telescope beam increases the correlations in the off-diagonal terms, i.e. the beam makes the different multipoles more correlated. This is because the telescope beam damping term breaks the orthogonality of the multipoles. See CovAppendix for a more detailed discussion and derivation. 
  
  We demonstrate that these differences due to the telescope beam in the correlation matrix and in the [S/N] per k carry over to an MCMC analysis in the foreground-free case. We perform the MCMC analysis for the monopole and quadrupole at the determined k_ max=0.24 h Mpc^-1 for the non-diagonal and diagonal covariance matrix cases, and quote results for both cases on <ref>. From these results we can determine that the errors on our model parameters, as seen on <ref>, increase when ignoring the covariance between different multipoles (with the exception of the α_∥ parameter, where the error slightly decreases). In both cases, we obtain unbiased parameter estimates.
  
  In the case of considering only a diagonal covariance matrix, the percentage uncertainties on the α_⊥ parameter (<ref>) are approximately 2 times smaller than those on α_∥, which is consistent with results from optical galaxy surveys (see e.g. <cit.>). When considering the covariance between different multipoles, the fractional uncertainties on α_⊥ become approximately 5 times smaller than on α_∥. We attribute this to the telescope beam having a significant effect on the correlation between different multipoles, but leave further investigation of the effect of the telescope beam and non-diagonal covariance matrix on the AP parameters for future work.
  
  We conclude that, when in the presence of a telescope beam, using a non-diagonal covariance matrix is important for the following reasons:
  
      
    * It makes a difference in the [S/N] per k bin result, increasing [S/N] in the low k limit but decreasing it in the higher k limit;
      
    * The different multipoles are non-negligibly correlated due to the telescope beam;
      
    * It decreases the uncertainties in most cosmological parameter estimates obtained using MCMC.
  
  
  
  
  
  
  
  
   §.§ Effect of foreground removal
  
  
  Here we aim to assess the validity of our foreground model. For our simulation, we have k^ min_⊥ = 0.004  h Mpc^-1 and k^ min_∥ = 0.006  h Mpc^-1, making the foreground damping scales N_⊥ k^ min_⊥ = 0.009  h Mpc^-1 and N_∥ k^ min_∥ = 0.013  h Mpc^-1 respectively. We note that, although we find N_⊥=2, N_∥=2 to fit our data well by eye, this does not mean that the same amount of power is being damped on both the perpendicular and parallel to the LoS directions. Indeed, when looking at the damping scales, we can see that N_∥ k^ min_∥ > N_⊥ k^ min_⊥, meaning that more power is being damped in the parallel to the LoS direction.
  
  To motivate our foreground model further, we attempt to compare it to a measurement of the power spectrum decomposed into perpendicular and parallel modes, P(k_⊥,k_∥). We compare P(k_⊥,k_∥) in the foreground-free case to P(k_⊥,k_∥) in the foreground removed case by plotting the ratio of these, and compare it to our foreground model B_FG(k, μ) (<ref>) with N_⊥=2, N_∥=2 (<ref>). We also plot the difference between these, finding that they are in agreement and that differences are below 10% on all scales. As seen in <ref>, both our model and the data show more power being damped on small k_∥ modes, as expected. This comparison was also carried out in <cit.>, which found similar agreement with a similar model.
  
  
  
  
  We perform an MCMC analysis with the foreground subtracted data in four different cases, first only considering the monopole and quadrupole only up to k_ max = 0.24  h Mpc^-1 and later considering the inclusion of the hexadecapole up to the limit found in the foreground-free case k_ max = 0.09  h Mpc^-1. We check that up to these limits, results are unbiased as in the foreground-free case (except for case 1, where results are biased). For most cases, we are varying the parameters {α_∥, α_⊥, T_ fσ_8,T_ b_σ_8, σ_v, P_ SN} and use the same priors limits as in the foreground-free MCMC analysis case. In one of the cases (case 3), we vary two additional parameters from our foreground model, namely N_⊥, N_∥, bringing the full list of parameters we vary to {α_∥, α_⊥, T_ fσ_8,T_ b_σ_8, σ_v, P_ SN, N_⊥, N_∥}. For N_⊥, N_∥, we impose flat positivity priors.
  
  Case 1: The foreground-free model. First, we consider the foreground-free model (<ref>) to demonstrate how it yields biased parameter estimates (specifically, the parameters α_⊥, T_ f σ_8, and T_ b_σ_8 become biased outside of the 2σ limit). In this case we do not include the foreground model in the covariance matrix.
  
  Case 2: The fixed foreground model. Next, we consider our foreground model (<ref>) and keep the N_⊥, N_∥ parameters fixed to the best fit guesses found by eye (N_⊥=2, N_∥=2). In this case we include the foreground model in the covariance matrix.
  
  Case 3: The varied foreground model. Here we consider the foreground model (<ref>) but let N_⊥ and N_∥ be nuisance parameters that we vary. Here we also include the foreground model in the covariance matrix. We also compare with the case of not including the foreground model in the covariance matrix, which causes the covariance matrix values to be larger and consequently we find that this increases errors in the parameters but does not cause them to become biased. This is relevant to the case of a real data analysis, where we would not know the fiducial N_⊥ and N_∥ values in advance to fix in the covariance matrix, and would probably need to adopt this more conservative case. If end-to-end simulations were available that allowed for N_⊥ and N_∥ to be accurately determined for real data, the less conservative case could be adopted instead. Alternatively, an iterative process could also be employed with real data. We would start by assuming N_⊥ and N_∥ in the covariance matrix, run a parameter estimation, re-generate the covariance based on the best-fitting values, and re-run the parameter estimation until there is convergence.
  
  Case 4: The k_ min-cut model. Finally, we investigate what happens when we exclude the largest scales where foreground subtraction has the most impact. We impose a k_ min limit on the foreground subtracted data and try to recover cosmological parameters using the foreground-free model (<ref>), which we have seen would yield biased parameter results if considering the full k-range (case 1). Here we do not include the foreground model in the covariance matrix. We find that the limit k_ min = 0.05  h Mpc^-1 is sufficient to then recover unbiased parameter estimates, and we quote the uncertainties on these on <ref>. Note that including the hexadecapole with this cut, or a more restricted cut, yields biased results due to the considerable impact that foreground removal has on the hexadecapole. For all parameters, the uncertainty obtained with the k_ min-cut method is larger than the uncertainty obtained using any other method. Furthermore, the varied N_⊥, N_∥ method yields smaller uncertainties, and does not require a prior selection of a k_ min limit.
  
  Results from the MCMC analyses for cases 1 to 3 can be found in <ref>. We quote the different uncertainties on the parameters for all cases on <ref>.
  
  
  
  
  We now consider our varied foreground model in more depth for case 3. Still letting N_⊥ and N_∥ be nuisance parameters and including the foreground model in the covariance matrix, we compare the MCMC analysis results when excluding or including the hexadecapole at a restricted range of k_ max = 0.09  h Mpc^-1. Results can be found in <ref>, and 1σ uncertainties on <ref>. We can see that as in the foreground-free case, adding the hexadecapole at a restricted range still allows us to retrieve unbiased cosmological parameters with smaller uncertainties than without it. Including the hexadecapole also seems to make the posteriors more Gaussian-like, in particular for N_⊥ and N_∥. This motivates further the inclusion of the hexadecapole in parameter estimation analyses, particularly in foreground removed data. We plot our foreground model with best fit parameters from the case 3 MCMC analyses (with and without the hexadecapole) in <ref>.
  
  
  
  
  
  
  
    §.§.§ Covariance matrix
  
  
  We find that the effect of foreground removal significantly impacts the covariance matrix, and discuss this effect further keeping in mind that this is specific to our choice of modelling, simulations and survey specifications (which determine the instrumental noise level). For real data, one would need realistic end-to-end simulations specific to a given experiment in order to robustly include the effects of foregrounds in the covariance matrix. 
  
  We compare the theoretical covariance matrix with and without foreground removal effects included. As seen in <ref>, where we include the foreground removal model with N_⊥=2, N_∥=2 in the covariance matrix, this makes a significant difference in the correlation between the different multipoles' large scale modes.
  
  When performing the different MCMC analyses with the monopole and quadrupole for the foreground subtracted case, we showed in case 3 that including the foreground model in the covariance matrix decreases errors in the cosmological parameters of interest. However, it requires knowing the best fit N_⊥, N_∥ beforehand, which might be unlikely with real data. Nonetheless, this test shows that we are able to retrieve unbiased parameter estimates using our foreground model in either case of including or not including foreground removal effects in the covariance matrix, but with different resulting parameter uncertainties.
  
  As an additional indicator of how well our models fit the simulation measurements, we have also looked at the reduced χ^2 of our best-fits: χ^2_ red = χ^2/dof, where dof is the degrees of freedom found by subtracting the model parameters from the number of data points. It is useful to look at χ^2_ red because if it is much larger than 1, that usually indicates an incorrect model or underestimated errors, and if it is much smaller than 1 then we could be overestimating the errors/overfitting. We calculate the χ^2_ red for the monopole, quadrupole (up to k = 0.24  h Mpc^-1) and for the hexadecapole (up to k = 0.09  h Mpc^-1). For the foreground-free case, we find χ^2_ red≃ 0.6. For the foreground removed case, we find χ^2_ red≃ 1 when including the foreground removal effects in the covariance matrix. As expected, when trying to fit the foreground-free model to foreground removed data, we obtain a best-fit χ^2_ red≃ 2.6 (and heavily biased cosmological parameter estimates, see <ref>), confirming the need for an appropriate foreground model.
  
  Although the reduced χ^2 is a useful check and indicator that our model is appropriate for fitting our measurements, our main findings and model validation come from the MCMC analyses, which recovers the cosmological parameters within 2σ errors of the fiducial values.
  
  
  
   §.§ Higher order multipoles
  
  
  Here we investigate the effect of including higher order multipoles in our analysis (see e.g. <cit.> and <cit.> for examples of higher order multipoles being considered in galaxy and halo two-point correlation functions, respectively). The P_6 64-pole (or hexacontatetrapole) encompasses non-linear velocity information, vanishing in the case of considering only the linear Kaiser RSD effect <cit.>. In the case of our simulations (ignoring beam and foreground effects), where non-linear velocity effects are present (such as the FoG effect), the 64-pole is non-zero but still expected to be very small. However, we show in <ref> that the 64-pole is significantly affected by the telescope beam and foreground removal effects similarly to the other multipoles, meaning its signal is boosted due to these systematic and instrumental effects.
  
  
  
  Regarding the correlation matrix, we again find that the beam and foreground removal effects significantly affect the correlations between the 64-pole and other multipoles, as seen in <ref>.
  
  
  
  We test the effect of adding the 64-pole to our parameter estimation pipeline, first in the foreground-free case. We find that we can add the 64-pole up to the same restricted range as the hexadecapole (k_ max = 0.09  h Mpc^-1), and that it does improve results by decreasing the errors on our parameters while maintaining the estimates unbiased within 2σ (see <ref>).
  
  
  
  We also tested whether including the 64-pole in the foreground removed case would make a difference, and indeed it did. When we added the 64-pole in the restricted range in our analysis (Case 3, varied foreground model with foreground effects included in the covariance matrix), we obtained unbiased results for all cosmological parameters up to k_ max = 0.08  h Mpc^-1, a slightly more restricted range than we find for the hexadecapole. This is likely due to how the foreground removal effect suppresses our cosmological signal in the covariance matrix, thus decreasing the error budget, combined with the 64-pole being highly non-linear.
  
  Removing the foreground effect from the covariance matrix yields a much larger error budget, and we tried including the 64-pole in this case. We found that indeed we obtain unbiased results up to k_ max = 0.09  h Mpc^-1 in this case but with very large uncertainties on our parameters, as seen in <ref>.
  
  Our results show that in the absence of foregrounds, the 64-pole can improve constraints without biasing parameter estimates. In the foreground removed case, the 64-pole does not improve constraints, but the 64-pole could still be useful in analysing foreground cleaned IM data. This is because its underlying cosmological signal is quite weak, but it is highly sensitive to the effects of the beam and foreground removal, or other unidentified systematics. It could thus be used as a further check for any residual systematic effects that might be present in the data.
  
  
  
  § CONCLUSIONS
  
  
  The aim of this work was to perform a comprehensive cosmological parameter estimation with the IM power spectrum multipoles, and investigate the level of uncertainties future surveys like the SKA can realistically obtain, requiring unbiased estimates. We used modelling and simulations of IM that account for effects of the telescope beam and foreground removal, and performed MCMC analyses on these. We also showed how the beam and foreground removal effects impact the covariance matrix and higher order multipoles. We summarise our main findings and conclusions below:
  
  
      
    * In the absence of foregrounds, we are able to retrieve unbiased estimates for cosmological parameters using our model, with below 10% percent level uncertainties (and for the transverse AP parameter, below 1% percent level uncertainty). Including the hexadecapole in our analysis does not bias parameter estimates if we only consider it at a restricted range. Even at a restricted range, including the hexadecapole significantly decreases parameter uncertainties for all cases considered. In particular, when including the restricted hexadecapole, we are able to retrieve the growth rate parameter (T_ fσ_8) with 8.8% uncertainty.
      
    * In the presence of a telescope beam and foreground removal effects, it is crucial to include the modelling of these in the covariance matrix as it makes a significant difference. In particular, the covariance matrices between different multipoles become non-negligible as these effects change the correlations between multipoles.
      
    * If we do not account for the effects of foreground removal in the modelling, we obtain significantly biased parameter estimates (see also the very recent study by <cit.> for the case of primordial non-gaussianity measurements).
      
    * We therefore develop a 2-parameter foreground model to account for the removal of modes that occurs due to foreground cleaning. With no assumptions about the foreground removal process (i.e. by letting these parameters vary), we use this model to try and recover unbiased cosmological parameter estimates and succeed, finding that the two extra free parameters are enough to model the effects of foreground removal in our case.
      
    * We find that we are able to model the effects of foreground removal, and recover the growth rate parameter (T_ fσ_8) uncertainty to be 13.3%, slightly larger than in the foreground-free case. The other cosmological parameters also experience a slight increase in uncertainties, but they are not as significant.
      
    * We investigate the effect of including the 64-pole in our analysis. We find that for the foreground-free case, it improves parameter uncertainties without biasing them, but worsens constraints in the foreground removed case. However, we propose that the 64-pole could be a useful tool to investigate systematic effects in foreground cleaned intensity mapping data, since it is highly sensitive to these.
  
  
  The results in this paper are dependent on our choice of simulation and noise modelling. It would be interesting in future work to test the robustness of our modelling against more complex foreground simulations, for example including polarization leakage. Considering additional noise and systematic effects, such as (1/f) noise and RFI flagging in our simulations would also be worthwhile.
  
  Although we looked specifically at the single dish intensity mapping case in this paper, our findings could also be relevant for the interferometer case. The interferometer case would have much better angular resolution, but additional complications due to different instrumental systematic effects. However, foreground removal effects are analogous and equally important for both cases, and  in particular has been applied to interferometric data already (see e.g. <cit.>).
  
  To further investigate our results for the covariance matrix, future plans include calculating the covariance matrix using a suite of simulations, and obtaining a more robust estimate of how systematic effects impact the IM covariance matrix.
  
  We hope our findings can be useful for analysing intensity mapping data from the MeerKAT single-dish survey <cit.>, in particular by using multipole expansion and our modelling prescriptions for understanding systematic effects. Our formalism can also help the preparation of forthcoming observations by providing realistic forecasts.
  
  
  
  § ACKNOWLEDGEMENTS
  
  We are grateful to the anonymous reviewer for very useful comments and suggestions that improved the quality of this paper. We thank Seshadri Nadathur, Mario Santos and Marta Spinelli for useful discussions and feedback. PS is supported by the Science and Technology Facilities Council [grant number ST/P006760/1] through the DISCnet Centre for Doctoral Training. SC is supported by STFC grant ST/S000437/1. AP is a UK Research and Innovation Future Leaders Fellow, grant MR/S016066/1, and also acknowledges support by STFC grant ST/S000437/1. This research utilised Queen Mary’s Apocrita HPC facility, supported by QMUL Research-IT <http://doi.org/10.5281/zenodo.438045>. We acknowledge the use of open source software <cit.>. Some of the results in this paper have been derived using the  and  package. We thank New Mexico State University (USA) and Instituto de Astrofisica de Andalucia CSIC (Spain) for hosting the Skies & Universes site for cosmological simulation products. 
  
  
  
  § DATA AVAILABILITY
  
  
  The data underlying this article will be shared on reasonable request to the corresponding author.
  
  
  
  
  
  
  
  mnras
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § COVARIANCE MATRIX IN THE PRESENCE OF A TELESCOPE BEAM
  
  
  In order to better understand why the telescope beam is increasing correlations between different multipoles, we consider a toy power spectrum model with and without the telescope beam effect. 
  
  We begin with the case of no telescope beam. First, assume we have a simple, isotropic matter power spectrum (no RSD): P_(k) = T_^2   b_^2 P_ m(k). Let us also set P_N = 0, b_ = 1 and T_ = 1 to obtain P_(k) = P_ m(k). This yields  σ^2 (k, μ) = σ^2 (k) = P_ m^2(k) / N_ modes(k). The sub-covariance matrices become:
  
  
      C_ℓℓ^'(k) = (2ℓ + 1)(2ℓ^' + 1)/2P_ m^2(k)/N_ modes(k)∫^1_-1 dμ ℒ_ℓ(μ) ℒ_ℓ^'(μ)   .
  
  
  In the absence of RSD and any other anisotropic effect in the power spectrum, the quadrupole and hexadecapole are null. Using <ref> we can confirm that the off-diagonal covariance matrix terms that include these multipoles are also null, as expected: C_02(k) = C_04(k) = C_24(k)=0.
  
  Next we consider the same power spectrum (P_N = 0, b_ = 1 and T_ = 1) but with a telescope beam damping term, such that P_(k, μ) = P_ m(k) B^2_⊥(k,μ) = P_ m(k) exp(-k^2 R_beam^2(1-μ^2)). We assume R_beam = 1 Mpc h^-1 for simplicity. This yields a covariance per k and μ bin of: 
  
  
      σ^2 (k, μ) = P_ m^2(k)  exp(-2k^2 (1-μ^2))/N_ modes(k)  ,
  
  
  and sub-covariance matrices given by:
  
  
      C_ℓℓ^'(k) =     (2ℓ + 1)(2ℓ^' + 1)/2P_ m^2(k)/N_ modes(k)
          ×∫^1_-1 dμ  e^-2k^2 (1-μ^2)ℒ_ℓ(μ) ℒ_ℓ^'(μ)   ,
  
  
  yielding the following off-diagonal covariance matrices:
  
  
      C_02(k)      = 5/2P_ m^2(k)/N_ modes(k)∫^1_-1 dμ ℒ_0(μ) ℒ_2(μ) e^-2k^2 (1-μ^2)
           = P_ m^2(k)/N_ modes(k)· 60k - 5√(2π)e^-2k^2(4k^2 + 3) erfi(√(2)k) /32k^3  ,
    
  
  
      C_04(k)      = 9/2P_ m^2(k)/N_ modes(k)∫^1_-1 dμ ℒ_0(μ) ℒ_4(μ) e^-2k^2 (1-μ^2)
           = P_ m^2(k)/N_ modes(k)·9/512 k^5 [3√(2π)e^-2k^2 (16k^4 + 40k^2 + 35) 
          ×erfi(√(2)k) +20(4k^3 - 21k)]   ,
  
  
  
      C_24(k)      = 45/2P_ m^2(k)/N_ modes(k)∫^1_-1 dμ ℒ_2(μ) ℒ_4(μ) e^-2k^2 (1-μ^2)
           = P_ m^2(k)/N_ modes(k)·45/4096 k^7 [4(304k^5 - 600k^3 + 1575k) 
           - 3√(2π)e^-2k^2(64k^6 + 208k^4 + 500k^2 + 525)erfi(√(2)k)]  ,
  
  
  where erfi(x) is the imaginary error function. From the results above, we can see that even in the absence of RSD, when we include the effect of a telescope beam the power spectrum does not have a vanishing covariance between the different multipoles.
  
  
  
  
  
  
  
  
  
  
  