Multi-tracing the primordial Universe with future surveys

. The fluctuations generated by Inflation are Gaussian in the simplest models, but may be non-Gaussian in more complex models. Since these fluctuations seed the large-scale structure, any primordial non-Gaussianity of local type will leave an imprint on the tracer power spectrum on ultra-large scales. In order to combat the problem of growing cosmic variance on these scales, we use a multi-tracer analysis that combines different tracers of the matter distribution to maximise any primordial signal. In order to illustrate the advantages that can be brought by multi-tracing future surveys, we consider two pairs of (galaxy survey, 21cm intensity mapping survey), one at high redshift ( 1 ≤ z ≤ 2 ) and one at very high redshift ( 2 ≤ z ≤ 5 ). The 21cm surveys are in interferometer mode, and are idealised versions of HIRAX and PUMA. We implement foreground avoidance filters and use non-trivial models of the interferometer thermal noise. The galaxy surveys are idealised versions of Euclid and MegaMapper. Via a simple Fisher forecast we illustrate the potential of the multi-tracer. Our results show a ∼ 20 − 70% improvement in precision on local primordial non-Gaussianity from the multi-tracer.


Introduction
Galaxy surveys in combination with cosmic microwave background (CMB) surveys and supernova Ia obervations have laid the foundations for precision cosmology.The next generation of surveys will deliver even higher precision by covering larger sky areas and deeper redshifts with better sensitivity.Each of these surveys is destined to make major advances in the precision and accuracy of the cosmological parameters of the standard LCDM model.At the same time, future surveys will look for signals beyond LCDM.One such signal arises if there is significant non-Gaussianity in the primordial perturbations that seed the CMB fluctuations and the large-scale structure.
Primordial non-Gaussianity (PNG) is a key probe of Inflation, which is currently the best framework that we have for the generation of seed fluctuations.The local type of PNG, described by the parameter f NL , if found to be nonzero, will rule out the simplest Inflation models [1].Furthermore, if 0 < |f NL | ≲ 1, then many other models can also be ruled out [2,3].In order to achieve this, a precision of σ(f NL ) < 1 is required.
The best current 1σ constraint comes from the Planck survey, via the CMB bispectrum [4]: f NL = −0.9± 5.1 (68% CL) . (1.1) The CMB power spectrum does not carry an imprint of local PNG, and the same is true for the matter power spectrum.However, the power spectrum of a biased tracer, such as galaxies, does contain a signal of local PNG on ultra-large scales -due to the phenomenon of scale-dependent bias [5][6][7].Constraints from future CMB surveys will not be able to reach σ(f NL ) < 1 due to cosmic variance.In the case of galaxy surveys, the current best constraints are σ(f NL ) ≳ 20 (e.g.[8,9]).Future galaxy surveys will improve on this, potentially reaching the CMB level of precision and beyond [10][11][12][13].Although galaxy power spectrum constraints on f NL will be a powerful independent check on the CMB results, single-tracer constraints are unlikely to achieve σ(f NL ) < 1, again because of cosmic variance, which is highest on precisely the scales where the local PNG signal is strongest.Since the local PNG signature in large-scale structure depends on the particular tracer, it is possible to improve the precision on f NL by combining the power spectra of different tracers.This multi-tracer approach can effectively remove cosmic variance [14][15][16][17].As a result, the multi-tracer of future galaxy surveys with intensity mapping of the 21cm emission of neutral hydrogen (HI IM surveys), can deliver significant improvements in the precision on f NL (e.g.[12,18,19]) and in some cases it can in principle reach σ(f NL ) < 1 (e.g.[20][21][22][23][24]).The multi-tracer can also mitigate some of the observational systematics in single-tracer surveys.
Radio-optical combinations have been intensively investigated (e.g.[12,[18][19][20][21][22][23]25]).They are well suited to the multi-tracer, since they combine surveys with very different clustering biases and systematics.Some of these investigations consider radio continuum surveys and photometric galaxy surveys.Here we follow [19,21,23,25] in considering pairs made of a spectroscopic galaxy survey and an HI intensity mapping survey.Previous papers have all used HI intensity mapping in single-dish mode.Here we use HI IM surveys in interferometer mode.Surveys in this mode have different limiting scales, foreground effects and thermal noise, compared to single-dish mode surveys.Single-dish mode surveys probe larger scales than those in interferometer mode and are thus expected to provide better precision on f NL .Nevertheless, it is interesting to see what is achievable in interferometermode surveys and whether it is worthwhile to combine them with galaxy surveys.In addition, an interferometer-mode survey on the proposed PUMA array [26] reaches higher redshifts than any proposed single-dish survey.Combining a PUMA-like survey with a galaxy survey like MegaMapper [27] enables us to investigate the highest (post-reionisation) redshift multi-tracer pair of an HI IM and a galaxy survey.
We use nominal simplified surveys since we perform Fisher forecasting which does not take into account most systematics.Our main aim is to demonstrate the advantage of the multi-tracer at high to very high redshifts with combined galaxy-HI probes, rather than to make forecasts for specific surveys.
The paper is structured as follows.In section 2, we provide a brief demonstration of multi-tracer power spectrum of two tracers.In section 3, we presented the characteristics of the galaxy and HI intensity mapping surveys we are using.We discuss the survey's instrumentation and brief noise requirements in section 3.1.The avoidance of foreground contamination for HI intensity mapping is presented in section 3.2.We also discuss the limited range of scales needed to obtain the cosmological information from both spectroscopic and radio surveys in section 3.3.The Fisher forecast is discussed in Section 4, along with how we combined the data from the two surveys.Section 5 provides a summary of our main conclusions.

Multi-tracer power spectra
At linear order in perturbations, the Fourier number density or brightness temperature contrast of tracer A is where δ is the matter density contrast, b A is the (Gaussian) clustering bias for the tracer A, µ = k • n is the projection along the line-of-sight direction n, and f = −d ln δ/d ln(1 + z) is the linear growth rate in LCDM.We use a single-parameter model for the biases [28], where we assume that the redshift evolution α A (z) is known.The power spectra for 2 tracers are given by δ so that where P is the matter power spectrum (computed with CLASS [29]).When B = A we use the shorthand P A ≡ P AA .
In the presence of local PNG, the bias acquires a scale-dependent correction: where 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 = (1 + z in )D in is the initial growth suppression factor, with z in deep in the matter era.For ΛCDM, g in = 1.274.The modelling of the non-Gaussian bias factor b Aϕ is in principle determined from halo simulations, but currently it remains uncertain [30][31][32][33][34][35].A simplified model is given by where the critical collapse density is given by δ c = 1.686 and p A is a tracer-dependent constant in the primordial non-Gaussian bias b Aϕ , arising from assembly bias.For the simplest case of a universal halo model, p A = 1 for all tracers A [5,6].
For our fiducial model, we use (2.6) with the following p A values (derived from simulations): for galaxies chosen by stellar mass [30], and for HI intensity mapping [31] p HI ≈ 1.25 . (2.8) In fact, the galaxy samples we use are not chosen by stellar mass but we are not aware of simulations to estimate p g values for these samples.Consequently we consider two further values: p g ≈ 1 and 1.5 , (2.9) and we find the effect on σ(f NL ) of these changes.
The Gaussian clustering biases are as follows: The background comoving number densities for Hα [42] and LBG [43] are taken as nHα (z) = 3.63 z −0.910 e 0.402 z − 4.14 For HI IM, the background brightness temperature is modelled as [44]:

Noise
For galaxy surveys, g = Hα or LBG, the shot noise is assumed to be Poissonian, so that the total observed auto-power is In HI IM surveys, HI = H or P, the dominant noise on linear scales is the thermal (or instrumental) noise [41,45], so that shot noise can safely be neglected and the observed HI IM power spectrum is where P therm HI depends on the model of the sky temperature in the radio band, the survey specifications and the survey mode (single-dish or interferometer).For the interferometer (IF) mode that we consider, the thermal noise power spectrum is [46][47][48][49]: T sys (z) THI (z) where k ⊥ = k 1 − µ 2 is the transverse mode, ν 21 = 1420 MHz is the rest-frame frequency of the λ 21 = 21 cm emission, t tot is the total observing time, r is the radial comoving distance, H is the conformal Hubble rate, and the system temperature is modelled as [48]: where T d is the dish receiver temperature, given in Table 2.The effective beam area is with dish diameter D d and the field of view is The physical baseline density of the array distribution is modelled following [48].For a given baseline of length L, we use the fitting formula for close-packed arrays given in [48] (see also [49,50]): where L s = N s D d .H (HIRAX-like) is a square close-packed array, while P (PUMA-like) is a hexagonal close-packed array with 50% fill factor.This means that half the dish sites are randomly empty, so the array is equivalent in size to that of twice the number of dish elements N d but with a quarter of the baseline density [50].Then it follows that N 2 s = N d for H and N 2 s = 2N d for P [48,51].The values of the fitting parameters a, b, c, e in (3.13) are given in Table 1, from [48] (their Appendix D, last equation).
The baseline densities Ln phys b for the H and P surveys are shown in Figure 2. Note that the thermal noise (3.9), which is ∝ 1/n phys b , tends to infinity at the maximum baseline D max .
Table 1.Specifications of HI IM surveys.H is HIRAX-like, with 256 dishes in phase 1 and 1024 in phase 2; P is PUMA-like, with 5k dishes in phase 1 and 32k in phase 2. N d is the number of dishes, D d is the dish diameter, T d is the dish receiver temperature, D max is the maximum baseline, t tot is the observing time.The parameters N s and a, • • • , e are given in (3.13).The cross-noise power spectrum between HI IM and galaxy surveys depends on the extent of overlapping halos that host the HI emitters and the galaxies.Typically it is neglected, as motivated in [52,53].Then the observed cross-power spectra are PgHI = P gHI . (3.14) The noise auto-power for galaxy surveys is shown in Figure 3.The HI IM thermal noise power spectra are displayed in Figure 4.

HI IM foreground avoidance
HI IM surveys are contaminated by foregrounds that are much larger than the HI signal.The cleaning of foregrounds from the HI signal relies on the fact that the foregrounds are nearly spectrally smooth (see e.g.[54,55] and references there in for recent advances).However, the cleaning is not perfect, and in some regions of Fourier space the signal loss is large.For simplified Fisher forecasts, we assume that cleaning is successful, except in the known regions of Fourier space where cleaning is expected to be least effective.Then we use foreground filters that avoid these regions.There are two regions of foreground avoidance for IF-mode surveys.factor In order to find the effect on f NL constraints of changing k fg , we consider a fiducial value that follows [56] and a less optimistic value: 2. A physical baseline in the IF array probes different angular scales k ⊥ at different frequencies (and hence different z).As a consequence, monochromatic emission from a foreground point source can contaminate the signal in a wedge-shaped region of Fourier space [57,58].This contamination occurs within N w primary beams of a pointing, where N w = 0 is the ideal case of no contamination, which occurs in principle if calibration of the interferometer is near-perfect [59].Currently this precision in calibration and stability is not achievable.In order to avoid the wedge region we therefore impose the filter [47,48,57,58,60] k ∥ > A w (z) k ⊥ where A w (z) = r(z)H(z) sin 0.61N w θ fov (z) . (3.17) We investigate the effect on f NL constraints of changing N w by using a fiducial value and a less optimistic value, following [50]: The fiducial value corresponds to foreground contamination from sources within 1 primary beam of the pointing, while the less optimistic value means that contamination extends to 3 primary beams.
Figure 5 illustrates D fg and A w .Foreground filters lead to the following modifications of the HI IM temperature contrasts: where Θ is the Heaviside step function.Note that P HI → D 2 fg ΘP HI and P gHI → D fg ΘP gHI .

Maximum and minimum scales
Since the local PNG signal in power spectra is on ultra-large scales, we only require linear perturbation theory.This means that we can exclude scales where linear perturbations break down, i.e. k < k max , where [61][62][63][64] Then we choose a very conservative value of k max,0 = 0.08 h Mpc −1 .The minimum wavenumber k min in each redshift bin is the fundamental mode k f , given by the comoving volume of the bin of width ∆z = 0.1 centred at z: The minimum and maximum scales are shown in Figure 6.

Fisher forecast
For a given redshift bin, the multi-tracer Fisher matrix for the combination of two dark matter tracers g and HI is where ∂ α = ∂ / ∂ϑ α and ϑ α are the parameters to be constrained.Here P is the data vector of the power spectra, P = P g , P gHI , P HI .(4.2) Note that the noise is excluded since it is independent of the cosmological and nuisance parameters.The noise enters the covariance [65][66][67][68]: Here ∆k and ∆µ are the bin-widths for k and µ, which we choose following [69,70]: In this paper our focus is on the improvements that the multi-tracer can deliver for pairs (galaxy, HI IM) of high and very high redshift surveys.We assume that most of the standard cosmological parameters have been measured by Planck [71].Since f NL affects the amplitude and shape of the power spectrum on ultra-large scales, we include in the Fisher analysis the cosmological parameters A s and n s that directly affect these properties of the power spectrum.Furthermore, surveys with the tracers A = g, HI are assumed each to have already established the redshift evolution α A (z) of the Gaussian clustering biases, leaving only the amplitude parameters b A0 in (2.2) to be constrained (see [28] for justification of this assumption).We also assume that the average HI brightness temperature THI (z) has been determined.In summary, we assume that the Planck survey and the single-tracer measurements have already determined the standard cosmological and astrophysical quantities.We incorporate the uncertainties in the bias amplitude parameters and in the amplitude and slope of the primordial power spectrum, since our focus is on the local primordial non-Gaussianity and on parameters that are strongly linked to it.The multi-tracer approach is known to significantly reduce the impact of uncertainties in the clustering biases (e.g.[14,15,20,21,72]).
These assumptions lead to optimistic forecasts for σ(f NL ) -but our aim is to find the relative improvement from the multi-tracer, rather than the precise values of σ(f NL ).Consequently, we consider the following cosmological and nuisance parameters: The fiducial values for the LCDM cosmological parameters A s , n s are taken from Planck [71] and we assume fNL = 0. Derivatives in (4.1) with respect to f NL , A s and b A0 are computed analytically, while the n s derivative is computed using the 5-point stencil method, with step-length 0.1 [73].
The Fisher forecast results, using the fiducial value of the non-Gaussian galaxy assembly bias parameter p g in (2.9) and the fiducial values of the foreground-avoidance paramaters k fg and N w in (3.16) and (3.18), are shown in Figure 7 and Figure 8, where we have marginalised over the bias nuisance parameters.These plots display the 1σ error contours for the multi-tracer pairs A ⊗ B = Hα ⊗ H (Euclid-like ⊗ HIRAX-like, high redshift) and LBG ⊗ P (MegaMapper-like ⊗ PUMA-like, very high redshift).The numerical values of σ(f NL ) are shown in Table 2.
Table 2. Marginalised 68% CL errors on f NL from galaxy surveys A = Hα (Euclid-like), LBG (MegaMapper-like) and HI IM inteferometer-mode surveys A = H (HIRAX-like) and P (PUMA-like).Both H and P have phases 1 and 2, with final number of dishes given in parentheses.Multi-tracer combinations are A⊗B = Hα⊗H (high redshift) and LBG ⊗ P (very high redshift).The final column shows the percentage improvement in precision from the multi-tracer for each single tracer.3 displays the changes in σ(f NL ) when using the alternative estimates for p g in (2.9) and the less optimistic values of k fg and N w in (3.16) and (3.18).Table 3.As in Table 2 but showing the effect on the multi-tracer σ(f NL ) when: (a) we change the non-Gaussian galaxy assembly bias parameter p g from the fiducial value p g = 0.55 to the alternative choices in (2.9) and fix p HI = 1.25 as in (2.8);(b) we enlarge the set of foreground avoidance parameters (k fg , N w ) from the fiducial values (0.005 hMpc −1 , 1) to include the less optimistic values (0.01 hMpc −1 , 3), and we consider all possible pairs (k fg , N w ).

Conclusion
We used Fisher forecasting to estimate the precision on local primordial non-Gaussianity that is achievable by combining future galaxy and HI intensity mapping surveys (taking foregrounds and interferometer thermal noise into account).We use a simplified model, incorporating uncertainties in the Gaussian clustering bias parameters and in the amplitude and slope of the primordial power spectrum, but fixing other parameters.Although this means that our σ(f NL ) forecasts are optimistic, it allows us to estimate the relative improvement in precision from the multi-tracer.The results of the Fisher analysis, shown in Figure 7 and Figure 8 and summarised in Table 2, confirm the effectiveness of using a multi-tracer combination of HI intensity mapping and spectroscopic galaxy samples, in order to enhance the precision in future measurements of local primordial non-Gaussianity.We have shown this for the first time in the case that the HI intensity mapping survey is performed in interferometer mode, as opposed to the single-dish mode used in previous forecasts.For the fiducial choice of non-Gaussian galaxy assembly bias parameter p g and of intensity mapping foreground parameters k fg , N w , the relative improvement from the multi-tracer over the best single-tracer constraint (from the galaxy samples) is ∼ 20 − 30% (see Table 2).
The multi-tracer also improves precision on n s and A s , as shown in Figure 7 and Figure 8.It is noticeable that the HI intensity mapping interferometer-mode surveys provide stronger constraints on n s , A s than the galaxy surveys, while the reverse is true for f NL .The reason is that the n s , A s signals are mainly dependent on intermediate scales, while the f NL signal is based on ultra-large scales.The HI interferometer-mode surveys provide better constraints on intermediate scales than the galaxy surveys.However, on ultra-large scales these surveys lose constraining power due to the foreground filters.
Our assumption on the tracer-dependent non-Gaussian galaxy assembly bias, i.e. that b Aϕ ∝ b A − p A as in (2.8) and (2.9), also affects galaxies and HI intensity mapping differently.
In the simplest universality model, p A = 1, there is no difference between galaxies and HI intensity in the form of the non-Gaussian galaxy assembly bias.In the fiducial case with p g < 1, it follows that b gϕ is larger than in the p = 1 model.Consequently, the constraints for galaxy surveys on f NL are improved compared to p = 1.By contrast, p HI > 1 means that b HIϕ is larger than in the p = 1 model, so that f NL constraints are degraded for HI intensity mapping surveys compared to the p = 1 model.
Single-tracer Fisher constraints for surveys similar to our mock surveys have been computed previously: e.g., [25] for Euclid; [40] for MegaMapper and PUMA; [50] for HIRAX and PUMA.These constraints are qualitatively compatible with ours -taking into account the following differences: • Previous work uses the universality model for b Aϕ , while we use eqrefe2.7,(2.8) and (2.9).
• Our mock surveys have differences in the redshift range, sky area and observing time (for HI intensity mapping) compared to the surveys considered in previous work.
• There are differences amongst the previous works, and with our paper, in the cosmological and nuisance parameters that are used and marginalised over.
To our knowledge, there are no previous works using the multi-tracer pairs that we analyse, so that comparisons of our multi-tracer results with previous results is not possible.The best single-tracer and multi-tracer f NL precision is delivered at the highest redshifts, by LBG (MegaMapper-like) and LBG ⊗ P (PUMA-like) (Figure 8).The LBG survey on its own already achieves σ(f NL ) < 1 (this is compatible with [40]), while P constraints (which are compatible with [50]) are weaker because of foreground contamination on ultra-large scales.When LBG and P are combined, we achieve an improvement in σ(f NL ) of ∼ 30% over the single-tracer LBG.
Our fiducial foreground avoidance parameters are fairly optimistic, but we expect that further advances in foreground treatments will emerge.In particular: • density reconstruction methods, such as forward model reconstruction, are being developed to retrieve more long wavelength radial modes lost to foregrounds (e.g [64,74,75]); • advances in the stability and calibration of interferometers will reduce the impact of the wedge and in principle will be able in future to reduce its impact to negligible.
In addition, there are theoretical uncertainties in the non-Gaussian galaxy assembly bias parameter p g .Our fiducial choice p g = 0.55 is based on galaxy samples selected by stellar mass, whereas the galaxy samples we use may have different p g .We therefore assessed the impact on σ(f NL ) of different p g choices [see (2.9)], as well as less optimistic intensity mapping foreground-avoidance parameters [see (3.16), (3.18)].The results are summarised in Table 3.It is apparent that the increase in p g leads to the largest increases in σ(f NL ).For example, if we consider the two best multi-tracer pairs, then Table 3 shows that with the foreground parameters fixed at their fiducial values, the increase in the errors when changing p g from 0.55 to 1.5 is 124% for Hα ⊗ H 1024 and 58% for LBG ⊗ P 32k.This makes sense, since increasing p g reduces b gϕ and hence reduces the scale-dependent bias.Note that the degradation is significantly weaker in the best-performing multi-tracer LBG ⊗ P 32k.For the foreground parameters, the increase in N w to 3 primary beams has very little impact on σ(f NL ), while increasing the radial mode suppression parameter k fg leads to a larger increase in σ(f NL ).In detail, increasing N w from 1 to 3 leads to increases in σ(f NL ) of 1% for Hα ⊗ H 1024 and 6% for LBG ⊗ P 32k.Changing k fg from 0.005 to 0.01h/Mpc leads to an increase in the errors of 9% and 7% respectively.
We conclude that the galaxy-HI multi-tracer constraints are not very sensitive to HI foreground avoidance parameters but are significantly degraded by increases in the non-Gaussian galaxy assembly bias parameter.Nevertheless, the improvements from the multi-tracer over the best single tracer is robust to the increases in p g , k fg and N w .The improvement ranges from 10 to 30%.
Fisher forecasts are of course highly simplified compared to observational and data reality.In practice there are numerous systematics to be modelled and there are major computational constraints on the simulations necessary to build a data pipeline.Nevertheless, a simplified Fisher analysis has shown that multi-tracing galaxy and HI intensity mapping interferometer-mode surveys at high to very high redshifts is worth exploring further.Furthermore, the multi-tracer will alleviate the impact of nuisance parameters and systematics in the galaxy and HI intensity mapping samples.
Finally, we note that our forecasts use a plane-parallel approximation in the Fourier power spectra.Future work will relax this assumption by including wide-angle effects (see e.g.[72,76,77]).

Figure 1 .
Figure 1.For the galaxy and HI IM surveys: clustering biases (left panel); comoving number densities (right panel, red, left-hand y-axis) and HI IM temperature (right panel, blue, right-hand y-axis).

Figure 3 .
Figure 3. Redshift evolution of galaxy shot noise.

Figure 5 .
Figure 5. Radial foreground damping function (left).Schematic to illustrate the foreground wedge; the wedge angle increases with redshift (right).

3 ]Figure 6 .
Figure 6.Minimum and maximum wavenumbers (left) and comoving volume (right) over the redshift range of the surveys.

Figure 7 .
Figure 7. Forecasted 1σ (68% CL) marginalised uncertainty contours for the local primordial non-Gaussianity parameter f NL , together with the cosmological parameters A s and n s , computed from the Hα (Euclid-like) galaxy survey and the H (HIRAX-like) intensity mapping survey.Top: 256 H dishes; bottom: 1024 H dishes.

Figure 8 .
Figure 8.As in Figure7but for the LBG galaxy survey and the P (PUMA-like) intensity mapping survey.Top: 5000 P dishes; bottom: 32,000 P dishes.