Constraining primordial non-Gaussianity by combining next-generation galaxy and 21 cm intensity mapping surveys

Surveys of the matter distribution contain ‘fossil’ information on possible non-Gaussianity that is generated in the primordial Universe. This primordial signal survives only on the largest scales where cosmic variance is strongest. By combining different surveys in a multi-tracer approach, we can suppress the cosmic variance and significantly improve the precision on the level of primordial non-Gaussianity. We consider a combination of an optical galaxy survey, like the recently initiated DESI survey, together with a new and very different type of survey, a 21 cm intensity mapping survey, like the upcoming SKAO survey. A Fisher forecast of the precision on the local primordial non-Gaussianity parameter fNL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\textrm{NL}}$$\end{document}, shows that this multi-tracer combination, together with non-overlap single-tracer information, can deliver precision comparable to that from the CMB. Taking account of the largest systematic, i.e. foreground contamination in intensity mapping, we find that σ(fNL)∼4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma (f_{\textrm{NL}}) \sim 4$$\end{document}.


Introduction
Constraining the local primordial non-Gaussianity (PNG) parameter f NL presents an important challenge for next-generation large-scale structure surveys.Local-type PNG affects the power spectrum of dark matter tracers by inducing a scale-dependent correction to the tracer bias [1,2].This correction is strongly suppressed on small to medium scales, but on very large scales, it is ∝ f NL (H 2 0 /k 2 ).This means that ultra-large survey volumes are required for high-precision constraints.
The Planck survey of the cosmic microwave background (CMB) gives the current state-of-the-art constraint σ(f NL ) = 5.1 from the bispectrum [3].Future CMB surveys will lead to an improvement on the Planck precision, but not by much and not sufficient to access σ(f NL ) 1.This level of precision would enable us to rule out many single-field and multi-field inflationary scenarios [4][5][6].The hope is that this goal can be achieved by using the extra modes in 3-dimensional surveys of the large-scale structure (see e.g.[7][8][9][10][11][12]).However, extracting f NL from these surveys faces formidable difficulties: • Observational systematics on very large scales [13][14][15][16][17][18][19].For our simplified Fisher analysis, we neglect these systematics, except for the foreground contamination of 21 cm intensity mapping (see below).• Astrophysical complications entailed in the modelling of tracer clustering, in particular, the uncertainties in modelling the scale-dependent bias [20,21].
We use a simplified model since a full treatment requires simulations, beyond the scope of this paper.• A theoretical systematic arising from the neglect of lensing magnification and other relativistic lightcone effects on the power spectrum, which can bias the best-fit f NL [22][23][24][25][26][27][28][29][30][31][32][33][34][35].Although the bias on best-fit can be large, the error σ(f NL ) is not significantly affected and we therefore omit these lightcone effects.• Cosmic (or sample) variance that dominates ultra-large scale measurements.
We deal with this problem via a multi-tracer approach -which also helps to mitigate uncorrelated systematics.
Recent constraints on f NL have used the BOSS galaxy [36] and eBOSS quasar [37] samples, with the tightest constraint to date of σ(f NL ) = 21 from eBOSS [38].The much larger volumes of upcoming surveys should facilitate significant improvement on this precision.But if we rely on individual surveys using the power spectrum of single tracers, it turns out that σ(f NL ) 1 is not achievable -even if we neglect systematics and complexities in scale-dependent bias [11].The fundamental problem is cosmic variance.
In this paper, our aim is not to identify survey combinations that can achieve σ(f NL ) 1 -for this purpose, one typically requires the very high number densities of photometric surveys [48,52].Instead, our aim is to investigate the improvements over single-tracer constraints when using a new type of spectroscopic large-scale structure survey -neutral hydrogen (HI) intensity mapping of the 21 cm emission line -in combination with spectroscopic galaxy surveys.Such a pair of spectroscopic samples has very different clustering biases and systematics, which accentuates the advantages of a multi-tracer analysis.
We use a simple Fisher forecast on mock surveys, which are based on next-generation surveys that are starting to observe or close to starting.The two surveys we have in mind are: the Dark Energy Spectroscopic Instrument (DESI), with the Bright Galaxy Sample (BGS) and Emission Line Galaxy (ELG) samples [62,63], and the Square Kilometre Array Observatory (SKAO), with intensity mapping surveys in Bands 1 and 2 [57,64].Galaxy number count surveys are well established, but a cosmological HI intensity mapping survey has not yet been implemented.However, pilot surveys on the SKAO precursor telescope MeerKAT have already been used to: (a) measure the cross-power between MeerKAT and WiggleZ galaxy surveys [65], and (b) start developing a pipeline for the planned SKAO surveys [66,67].
We combine pairs of these surveys at low and high redshifts, using the Fisher single-tracer constraints in the non-overlap volume and the multi-tracer constraints in the overlap volume.The results show a reasonable improvement for the low-redshift pair, and a significant improvement for the high-redshift pair, compared to the standard single-tracer forecasts for σ(f NL ).The combination of all low-and high-redshift surveys improve on Planck, with σ(f NL ) ∼ 4.This constraint is based on avoiding foreground contamination of HI intensity mapping on very large scales.Foreground contamination of HI intensity mapping has a strong effect on the precision: if we neglect this foreground contamination, we find that σ(f NL ) ∼ 1.

Multi-tracer power spectra
At first order in perturbations, the density (or temperature) contrast of a tracer A is related to the matter density contrast δ by where b A is the (Gaussian) clustering bias, µ = k • n is the projection along the line-of-sight direction n, and f is the linear growth rate.We assume b A has a known redshift evolution β A , following [68], so that where the amplitude b A0 is a free parameter.We highlight the point that a multi-tracer approach reduces the impact of this assumption on σ(f NL ) [35].
The Fourier power spectra at tree-level are given by where the linear matter power spectrum (computed using CLASS [69]) is In the presence of local PNG, the bias acquires a scale-dependent correction [1,2,21]: Here T (k) is the matter transfer function (normalized to 1 on very large scales), D is the growth factor (normalized to 1 today) and g in is the initial growth suppression function, defined deep in the matter era.For ΛCDM [70] The non-Gaussian bias factor b Aφ is halo-dependent and should be determined with the aid of simulations [20,71,72].A simplified model reduces it to where the critical collapse density is given by δ c = 1.686 and p A are halodependent constants.In the simplest (universal) halo model, p A = 1 for all tracers A. But this model is not consistent with simulations [20,71,72].An improved (but still over-simplified) model allows the constant p A to vary with tracer.For galaxies chosen by stellar mass, simulations indicate that the rough approximation [71] p g ≈ 0.55 , (8) is an improvement, and we use this for the DESI-like samples.For HI intensity mapping, we use the approximation [72] p H ≈ 1.25 . ( The f NL terms in the power spectra are of order H 2 0 /k 2 on ultra-large scales, k k eq , where T ≈ 1. Local PNG therefore dominates the power on ultralarge scales -but these are also the scales where cosmic variance is largest.We deal with the cosmic variance via a multi-tracer analysis, which includes the information from all auto-and cross-power spectra of the two (or more) tracers (see section 4).
In both cases, we have a low-and a high-redshift survey.The sky and redshift coverage of the individual and overlapping surveys are given in Table 1.The one-parameter model (2) for the Gaussian clustering bias of the galaxy surveys follows [29] b The background comoving number density ng for the BGS and ELG samples is modelled following [29], which leads to the smoothed curves shown in Figure 1.
For HI IM, the background brightness temperature is modelled via the fit given in [66,76]: which is also shown in Figure 1.

Noise
In galaxy surveys, the shot noise (assumed to be Poissonian) is Then the total signal that we measure for the galaxy auto-power spectrum is where P gg is given by (4).
In IM surveys, there is shot noise, but also thermal (or instrumental) noise.The Poissonian shot noise (see [77] for non-Poissonian corrections) is derived in a halo-model framework and given by [74,78], where nH is the effective average comoving number density of HI, ρH is the average comoving HI density, N h is the halo mass function (average comoving halo number density per mass) and M H is the average HI mass in a halo of mass M .However, on the linear scales that we consider, the shot noise is much smaller than the thermal noise and can be safely neglected [74,78].The thermal noise in HI IM depends on the sky temperature in the radio band, the survey specifications and the array configuration (single-dish or interferometer).For the single-dish mode of SKAO-like surveys, the thermal noise power spectrum is [79][80][81]: T sys (z) TH (z) where ν 21 = 1420 MHz is the rest-frame frequency of the 21 cm emission, t tot is the total observing time, and the number of dishes is N d = 197 (with dish diameter D d = 15 m).The system temperature is modelled as [82]: where T d is the dish receiver temperature (see [73]).The total signal is then The noise for all surveys is shown in Figure 2.
In the case of the cross-power spectrum P gH , cross-shot noise arises if the halos that host galaxies and HI overlap.Following the arguments given in [73,83], we assume that the cross-shot noise may be neglected.The total cross-power signal is then PgH = P gH with P shot gH ≈ 0 .

Intensity mapping beam and foregrounds
HI IM surveys like the SKAO surveys have poor angular resolution, which results in power loss on small transverse scales, i.e. for large k ⊥ = (1 − µ 2 ) 1/2 k.This effect is typically modelled by a Gaussian beam factor [79]: (20) HI IM surveys are also contaminated by foregrounds that are much larger than the HI signal.Since these foregrounds are spectrally smooth, they can be separated from the non-smooth signal on small to medium scales.However, on very large radial scales, i.e. for small k = µk, the signal becomes smoother and therefore the separation fails.A comprehensive treatment includes simulations of foreground cleaning of the HI signal (e.g.[15,18]).For a simplified Fisher forecast we can instead use a foreground avoidance approach, by excising the regions of Fourier space where the foregrounds are significant.This means removing large radial scales, which can be modelled via an exponential suppression factor [75,84]: We choose a typically used value: Figure 3 illustrates the consequences of foreground avoidance for the monopoles of the HI power spectra, Note the consequence of foreground avoidance for the monopole of the HI power spectra P AH,0 , which are suppressed on large scales.It is clear that large-scale (small k) power in P AH,0 is lost and that the effect of local PNG on these large scales is suppressed by foreground avoidance.Note that f NL reduces the HI power on very large scales at z = 0.3, since b H < 1.25, which leads to b Hφ < 0 in (7).Also note that the beam effect on P AH,0 , which tends to suppress small scale (large k) modes, is negligible at the low redshift, but becomes apparent at higher redshift.
In summary, the effects of the radio telescope beam and radio foreground avoidance lead to the following modifications of the power spectra P AH : We note that these equations are consistent with the recent paper [85] on foreground cleaning via a transfer function, which shows that the effect of foregrounds is the same in the HI auto-power as in the cross-correlation power spectra with galaxies.However, we should also point out that the modelling of foreground removal is an ongoing project which has not yet reached maturity.

Fisher forecast
For a multi-tracer combination of two dark matter tracers g and H, the data vector of power spectra is We exclude the noise from the data vector, since the noise is independent of the cosmological and nuisance parameters.(The noise appears in the covariance, as given below.)The standard cosmological parameters are measured on medium to small scales and are effectively independent of the ultra-large scale parameter f NL .Fixing their values could bias the best-fit value of f NL but is unlikely to have any significant impact on σ(f NL ).We include in the Fisher analysis the cosmological parameters that directly affect the large-scale amplitude (σ 8,0 ) and shape (n s ) of the power spectrum.We therefore consider the following cosmological and nuisance parameters: Here we assume that the degeneracies between b A and σ 8,0 , and between TH and σ 8,0 , have been broken by other surveys that focus on medium to small scales.
The covariance for the multi-tracer power spectra is given by [86,87]: Here ∆k and ∆µ are the bin-widths for k and µ and k f is the fundamental mode, corresponding to the longest wavelength, which is determined by the comoving survey volume of the redshift bin centred at z: Then the multi-tracer Fisher matrix in a redshift bin is where ∂ α = ∂ / ∂ϑ α .We choose the bin-widths and k min following [88][89][90]: The smallest scale (largest k) is chosen to ensure that linear perturbations remain accurate, since we do not require information from nonlinear scales [91]: The galaxy and HI IM surveys do not have the same sky area and redshift ranges (see Table 1).The overlap redshift ranges are obvious; for the sky areas, we assume a nominal overlap area of 10,000 deg 2 .The multi-tracer applies in the overlap redshift range and overlap sky area.However, we can add the independent Fisher matrices obtained in the non-overlapping regions of the two individual surveys to the multi-tracer Fisher matrix in the overlapping region [73].In detail, the non-overlap region is • the non-overlap sky area of each survey, across the full redshift range of each survey; • the overlap sky area, across the non-overlap parts of the redshift ranges of each survey.
Then the full Fisher matrix (denoted by g ⊗ H) is Figure 4 shows the 1σ error contours computed from this Fisher matrix after marginalising over the bias parameters in (27).We use the fiducial values σ 8,0 = 0.8102, n s = 0.9665 and f NL = 0, keeping other cosmological parameters fixed to their Planck 2018 best-fit values [92].Table 2 lists the marginalised σ(f NL ) constraints, including the improvements (in parentheses) that follow when we ignore the HI IM foregrounds, i.e. when k fg = 0.The results for BGS, Band 1 and Band 2 are consistent with [35], which uses the angular power spectrum rather than the Fourier power spectrum used here.
Table 2 Cumulative marginalised error σ(f NL ) for: the single-tracer surveys; the low-z and high-z multi-tracer pairs -including single-tracer information from the non-overlap volumes as in (33); and the sum of the multi-tracer pairs.Values in parenthesis correspond to ignoring foregrounds, i.e., with k fg = 0.

Conclusion
We applied a simplified Fisher forecast analysis in order to gain insights into the improvements on σ(f NL ) from a new multi-tracer combination of largescale structure surveys that will be made possible by the upcoming HI intensity mapping surveys with the SKAO.For the spectroscopic galaxy samples, we used mock surveys based on DESI.This allowed for multi-tracer pairs at low and high redshifts, in order to see the level of improved precision from the higher volume at high redshifts.
We included all available Fisher information from these mock surveys, using the multi-tracer Fisher matrix on the overlap volume (sky area and redshift range) and the single-tracer Fisher matrices on the non-overlap volumes, as in (33).Our results for σ(f NL ) are summarised in Table 2, which gives the singletracer results for the 4 samples, the full multi-tracer results from (33) for the low-and high-redshift pairs, and finally the results from the Fisher sum of these pairs.The impact of foreground avoidance can be seen from the results (in parenthesis) when foregrounds are ignored, corresponding to k fg = 0 in (21).
As expected, the foreground avoidance severely degrades the constraining power of the HI intensity surveys.Without foregrounds, the high-redshift Band 1 survey would perform as well as the multi-tracer combination with the ELG survey: σ(f NL ) = 1.8.The reduced precision on f NL for single-tracer HI intensity constraints is consistent with the findings of [15], where foreground cleaning is simulated.However, using the multi-tracer with a spectroscopic galaxy survey helps to reduce the effect of foregrounds on σ(f NL ), as shown by Table 2.
The multi-tracer analysis is also expected to mitigate other systematics in the galaxy and HI intensity samples.For example, HI intensity mapping is also affected by radio frequency interference, receiver noise, calibration errors and polarisation leakage.Incorporating these and the galaxy systematics requires extensive simulations, beyond the scope of our paper.

Fig. 2
Fig. 2 Shot noise for galaxy surveys (left) and thermal noise for IM surveys (right).

Table 1
Sky area and redshift range of surveys.