Radio-optical synergies at high redshift to constrain primordial non-Gaussianity

We apply the multi-tracer technique to test the possibility of improved constraints on the amplitude of local primordial non-Gaussianity, $f_{\mathrm{NL}}$, in the cosmic large-scale structure. A precise measurement of $f_{\mathrm{NL}}$ is difficult because the effects of non-Gaussianity mostly arise on the largest scales, which are heavily affected by the low statistical sampling commonly referred to as cosmic variance. The multi-tracer approach suppresses cosmic variance and we implement it by combining the information from next-generation galaxy surveys in the optical/near-infrared band and neutral hydrogen (HI) intensity mapping surveys in the radio band. High-redshift surveys enhance the precision on $f_{\mathrm{NL}}$, due to the larger available volume, and HI intensity mapping surveys can naturally reach high redshifts. In order to extend the redshift coverage of a galaxy survey, we consider different emission-line galaxy populations, focusing on the H$\alpha$ line at low redshift and on oxygen lines at higher redshift. By doing so, we cover a wide redshift range $1 \lesssim z \lesssim 4$. To assess the capability of our approach, we implement a synthetic-data analysis by means of Markov chain Monte Carlo sampling of the (cosmological+nuisance) parameter posterior, to evaluate the constraints on $f_{\mathrm{NL}}$ obtained in different survey configurations. We find significant improvements from the multi-tracer technique: the full data set leads to a precision of $\sigma(f_{\mathrm{NL}})<1$.


Introduction
The ΛCDM model can successfully explain a wide range of cosmological observations, but still leaves some issues open.In addition to dark matter and dark energy, a phase of exponential expansion -cosmological inflation -is required to produce the fluctuations in the density field that seed the large-scale structure of the Universe.A key probe of inflation is primordial non-Gaussianity [1][2][3].
The simplest inflationary models generate primordial fluctuations that follow a Gaussian distribution, but many scenarios of inflation predict departures from Gaussianity, which may be quantified by the local primordial non-Gaussianity parameter, f NL [3,4].The effects of non-Gaussianity can be probed in various ways, the most prominent being a measurement of a non-vanishing primordial bispectrum [5][6][7][8][9].However, the power spectrum of clustering of biased tracers of the underlying matter density field also exhibits a peculiar scale-dependence on the largest scales [10,11].A measurement of f NL ̸ = 0 from the power spectrum would allow us to rule out entire classes of inflationary models.
Currently, the tightest constraints on f NL come from the bispectrum of the cosmic microwave background [5], but observations of the large-scale structure with next-generation cosmological surveys will attain competitive constraints on f NL [6,[12][13][14][15][16][17][18][19][20][21].In the optical/nearinfrared band, experiments like the European Space Agency's Euclid satellite, 1 the Nancy Grace Roman Space Telescope2 (Roman hereafter), the Legacy Survey of Space and Time at the Vera C. Rubin Observatory, 3 and the Dark Energy Spectroscopic Instrument DESI, 4 will perform large galaxy surveys with high sensitivity and up to high redshifts.At much longer wavelengths, the SKA Observatory5 (SKAO) will cover a large area of the southern sky using the neutral hydrogen (Hi) intensity mapping technique to map the unresolved emission in the 21-cm line from radio galaxies residing in the large-scale structure.Optical and radio galaxy surveys trace two (mainly) independent tracers of the underlying matter distribution and they are therefore complementary.This allows us to apply the multi-tracer method [22][23][24][25] to use jointly the information from both the galaxies and the Hi distribution: by correlating independent biased tracers of the large-scale structure the cosmic variance error is strongly mitigated.
The paper is organised as follows: in Section 2 we present the model of the power spectrum in auto-and cross-correlation, including the effects of non-Gaussianity on the power spectrum and we present the formalism of the multi-tracer technique.In Section 3, we describe the mock data sets employed in the analysis.In Section 4, we evaluate the signal-to-noise ratio to asses the significance of the detection of primordial non-Gaussianity.The analysis is detailed in Section 5, and the results in Section 6. Conclusions are drawn in Section 7.

Power spectrum and multi-tracer technique
The large-scale structure of the Universe can be studied by measuring its statistical properties.Two-point statistics have been the major target of cosmological observational campaigns, including next-generation surveys such as those mentioned above.In Fourier space, the twopoint statistic relevant for inhomogeneities in the matter distribution is the power spectrum of the perturbations of the matter over-density.Since the matter distribution cannot be observed directly, to sample it we use (biased) tracers, such as galaxies or Hi.In the linear regime, the power spectrum of two tracers can be written as a function of the wavevector k and the redshift z as where A, B are tracer labels and µ is the cosine of the angle between the wavevector k and the line-of-sight direction (which we take as opposite to that of incoming photons).In the expression above, b A is the linear bias of tracer A, f = −d ln D/d ln(1 + z) is the growth rate, where D(z) is the linear growth factor, and the term f µ 2 describes the effects of redshift-space distortions (RSD) in the linear regime.Finally, P lin ∝ D 2 (z) T 2 (k) P R (k) is the linear matter power spectrum, with T (k) the transfer function (normalised such that T = 1 for k → 0) and P R (k) the dimensionless power spectrum of the primordial curvature perturbation.
When the two tracers are the same, A = B, we talk about auto-correlation power spectrum, while if A ̸ = B, it is a cross-correlation.Moreover, analysing auto-correlations by themselves will be referred to as a single-tracer analysis, whereas considering both auto-and cross-correlations at the same time leads to a multi-tracer analysis.
The effect of primordial non-Gaussianity of the local type is to introduce a scale dependence in the linear bias [3,4,10,18,19,[26][27][28]: (2.2) Here, b Aϕ is the primordial non-Gaussian bias factor, where ϕ is the potential deep in the matter era.This bias, like the Guassian bias b A , depends on the properties of tracer Abut the relation between b A and b Aϕ is still poorly understood.Simulations indicate that this relation is sensitive to the assembly history and the selection criteria of the tracer [18,19,27,28].Since b Aϕ is completely degenerate with f NL , we cannot constrain f NL without a theoretical model for b Aϕ .Typically, a strong theoretical prior has been imposed by assuming universality of the halo mass function.In this case, where δ c ≃ 1.686 is the critical density contrast for collapse in the linear regime.Nearly all observational constraints from galaxy surveys on f NL to date rely on this universality assumption -with a few exceptions (e.g.[29,30]).Similarly, most forecast constraints are based on assuming the universality relation.Further progress in uncovering the relation between b A and b Aϕ is needed in order to break the degeneracy without invoking strong theoretical priors.Clearly the best-fit values of f NL will be biased by using incorrect models of b Aϕ , and f NL errors will be increased by uncertainties on b Aϕ in a given model.However, we highlight the point that a multi-tracer analysis is significantly more robust to these problems than a single-tracer analysis [20].Here we assume that the degeneracy is broken and that the constraints on f NL are not very sensitive to the detailed form of the relation between b Aϕ and b A .Then we use Eq.(2.3) as the simplest proxy for this relation.
The bias correction for non-Gaussianity is proportional to T (k) −1 k −2 , which implies that the effects become manifest at ultra-large scales, T → 1, where the cosmic variance is larger due to the decrease of observable independent modes.To overcome this limit, we apply the multi-tracer technique, which relies on the fact that the correlation of independent biased tracers of the same underlying matter distribution allows us to eliminate the error of the cosmic variance [22][23][24].This technique consists in analysing simultaneously the auto-and cross-correlation power spectra of two tracers in a single, joint data vector.Namely, in our case, P = P g g , P g Hi , P Hi Hi . (2.4) The full covariance matrix associated to the multi-tracer power spectra is given by [18,31] where PAB = P AB + P noise AB δ K AB , P noise AB is the (statistical or instrumental) noise, and δ K AB is the Kronecker delta.The term N k (z) represents the number of independent modes within an observed volume V (z) at given z and k, i.e.
with ∆k(z) and ∆µ the width of the bins in k = |k| and µ, respectively.(Note that, in principle, the k-binning may depend on redshift.)For a survey covering a sky fraction f sky , the comoving volume observed at redshift z is where χ(z) is the comoving distance up to redshift z.
On the ultra-large scales where the signal of local primordial non-Gaussianity is strongest, general relativistic 'projection' effects -i.e.effects from observing on the past lightcone -can also become non-negligible [32][33][34][35][36][37].The observed overdensity of tracer A is then of the form where n is the line-of-sight direction, v is the peculiar velocity of matter (assumed equal to that of tracer A).The second term is the standard redshift-space distortion of overdensity.The third term is the lensing distortion of overdensity, where s A is the magnification bias of tracer A and κ is the lensing convergence.The last term, δ dp A , contains Doppler and potential distortions of overdensity (including Sachs-Wolfe, integrated Sachs-Wolfe and time delay terms).A natural way to include all these effects in the 2-point correlations is via the angular power spectra (which also naturally include the full wide-angle effects).This allows for a theoretically correct analysis of local primordial non-Gaussianity in clustering, as performed in [38] and then generalised by e.g.[12,15,39].
The question is: what is the consequence for f NL of neglecting the general relativistic projection effects?This was first addressed in the case of lensing convergence in photometric surveys by [40] (see e.g.[41] for recent work on this case).Then all general relativistic effects were included (in spectroscopic surveys) by [42] (for recent work, see e.g.[43], which also uses a multi-tracer analysis, and [44]).It turns out that the estimate of the best-fitting value of f NL can be strongly biased by the neglect of relativistic projection effects -but the error on f NL is not significantly affected [43].Since we are performing forecasts of σ(f NL ) and not measurements of f NL , it is reasonable for us to neglect the relativistic projection effects and use the Fourier power spectra given by Eq. (2.1).

Data sets
The first step for building a synthetic data set both for galaxies and Hi is defining the observed cosmic volume, i.e. the depth and the sky coverage of the reference experiments.The observed volume, and therefore the redshift, define in turn the scales that can be probed by the surveys.The relation between the redshift and the observed scale arises from the fact that the observed sky volume is a redshift dependent quantity.Assuming that we are observing a sky cube of volume V , the minimum wavenumber, k min , corresponding to the largest scale observed, can be expressed as The interval in k is then different for each observed redshift and it is defined by this minimum value, decreasing with the redshift, while the upper limit is set to be the largest k before non-linearities need to be taken into account [14,45,46]  where n s is the slope of the primordial curvature power spectrum.This choice, resulting in 0.085 < k max Mpc < 0.16, is made in order to avoid non-linear scales where theoretical understanding is still poor.Moreover, since the primordial non-Gaussianity signal is strongest on the largest scales, we prefer to be conservative and not include even mildly non-linear scales.
In our analysis, we adopt 12 bins in the range z ∈ [0.85, 4.0], 10 linear bins in k, and 10 equi-spaced bins in µ spanning all the possible values between µ = −1 and µ = 1.The redshift bins have a variable width ∆z ∼ 0.2, the width being chosen to follow the transition between galaxy types [47], and the same division was adopted for the Hi sample in order to match the binning between radio and optical tracers.Finally, the binning in k is done using 10 bins in order to ensure that the condition ∆k ≥ k min is satisfied.

Galaxy samples
As a benchmark, we consider an emission-line galaxy survey.We assume a target sample at low redshift of Hα galaxies between z ≃ 0.9 and z ≃ 1.7.At higher redshift, oxygen lines, i.e. the Oiii line and the Oii doublet, will be used to identify galaxies, considering also contributions from the overlapping Hβ line.This is in line with the capabilities of present and near-future cosmological experiments like Euclid and Roman.Note, however, that for the former, oxygen-line galaxies will in fact be interlopers, which must be identified correctly to reach the completeness and purity required for the target sample.But once they have been identified, they can be used as additional samples of galaxies to extend the redshift range [47; see also 48,49].Specifically, we consider a sample of Hβ+Oiii galaxies between z ≃ 1.9 and z ≃ 2.6 and a sample of Oii galaxies in the range 2.8 ≲ z ≲ 4: the transition between different galaxy types is represented in the left panel of Fig. 1.For the sky coverage, we adopt 15 000 deg 2 .The impact of a different sky coverage is addressed in Appendix A.
Regarding the properties of the various galaxy samples, we parametrise their (linear) bias as a function of the galaxy type, the redshift, and the flux limit of the survey.We   tracer (z) HI Hα galaxies Hβ+OIII galaxies OII galaxies Linear bias of Hi and of different galaxy types.The galaxy bias is computed varying the flux limit F c : dotted and dash-dotted lines represent respectively higher and lower flux limits compared to the reference value F c = 2 × 10 −16 erg s −1 cm −2 (solid).Hα-line galaxies are characterised by a bias closer to unity and do not present a relevant dependence on the flux limit.All other galaxy types considered in our data set at higher redshift have a higher bias and a stronger dependence both on the redshift and on the flux limit of the survey.The green dashed line corresponds to the Hi bias parametrised as in Eq. (3.9).
adopted the 5-parameter model described in [50], which provides the fitting formula Here x is related to the flux limit of the survey through x = log 10 L min = log 10 4 π D 2 L (z) F c , where D L is the luminosity distance and cgs units are used.The five parameters are set according to the galaxy type: their values are presented in Table 1.The dependence of the galaxy bias on the redshift is displayed in Fig. 2 for different galaxy types and flux limits.We also investigate other possible parametrisations in Appendix A.
The expected number density of galaxies per unit volume in each redshift bin can be calculated by integrating the galaxy luminosity function Φ(L, z), i.e.
where the minimum luminosity is determined by the flux limit of the survey, while the maximum luminosity can be safely set to infinity.Table 2. Values of the parameters that characterise the luminosity function of the Hβ+Oiii and the Oii galaxy samples.Units are Mpc −3 for the number densities Φ * and erg s −1 for the luminosities L * .
order to capture the differences between galaxy types, following [47].For the Hα sample we adopted the broken power law luminosity function described in [51], where Φ * is fixed to its value at z = 0 while the reference luminosity L * is given by The values of the parameters that characterise the Hα luminosity function are α = −1.587, For the Hβ+Oiii and the Oii samples, we used a Schechter luminosity function, where the values of α, Φ * (z i ) and L * (z i ) are given in Table 2 for each galaxy type, according to the estimates of [52].
Further comments and analysis on the observed number density of interlopers are presented in Appendix A.
For a galaxy survey, the noise in the power spectrum corresponds to a shot-noise term given by the inverse of the observed galaxy number density in a given redshift bin, namely where zi is the mean redshift of the ith redshift bin and ni is the number density of galaxies per unit volume in that bin.The impact of shot noise is the more severe the higher the redshift, as a consequence of the flux limit of a spectroscopic survey: galaxies at high z are fainter and difficult to detect, and this results in a lower observed number density.We consider a flux limit of F c = 2 × 10 −16 erg s −1 cm −2 as reference and use the observed number galaxies calculated as in [47].Then we also let the flux limit vary from F c = 10 −16 erg s −1 cm −2 to F c = 3 × 10 −16 erg s −1 cm −2 to evaluate the impact of the observed galaxy number density on the forecast.Note that we choose the aforementioned values to include the capabilities of present and upcoming surveys, such as Euclid [53] and Roman [54].

Hi intensity mapping
The intensity mapping survey is considered to have the properties of the corresponding SKAO survey in the mid band proposed in the SKA Cosmology Red Book 2018 [55].Such a survey will cover a sky fraction f sky = 0.48 in single-dish mode, covering the whole redshift range considered for the analysis.We model the Hi bias according to the following parametrisation [56], b Hi (z) = 0.842 1 + 0.823 z − 0.0546 z 2 . (3.9) The dominant noise term in the case of the intensity mapping is due to thermal noise of the instrument, which we model as (3.10) In the expression above, t tot is the total observing time, N d is the number of dishes, while the Hi temperature and the system temperature are modelled respectively as . (3.12) For T Hi , we use the parametrisation given in [5,57,58], which captures the dependence of the average Hi brightness temperature on the comoving Hi density fraction, Ω Hi , whose full expression is given in Eq. (2.1) of [59].For T d , which is the dish receiver temperature, we instead follow [14].It is worth noticing that there is a difference of orders of magnitude between T Hi and T sys , meaning that the cosmological Hi signal is not the dominant contribution.Note also that the thermal noise is scale-independent and increases with redshift.
In principle, a shot noise term should be present as well, as the continuous Hi emission is fundamentally a collection of the emission of all unresolved Hi-line galaxies.It can be derived in a halo-model framework as [15,56,60] where ρ Hi is the average comoving Hi density, n h is the halo mass function, and M Hi the average mass of neutral hydrogen within a halo of mass M .However, on linear scales, and especially at high redshift, the amplitude of shot noise is much lower that that of thermal noise.Since this is the regime we are mainly interested in, we neglect the shot noise term in our analysis.Intensity mapping observations also suffer from a loss of signal at small scales in the perpendicular direction due to the low angular resolution.In an idealistic scenario, this can be modelled with a Gaussian beam perpendicular to the line-of-sight direction.It results into a damping factor affecting perpendicular modes k where the size of the beam θ b of a dish of diameter D d is given, in radians, by [64] In the case of the intensity mapping the foreground contamination must be taken into account [65][66][67][68].Foregrounds mostly affect the largest radial scales, where they cannot be separated from the cosmological signal.We decide to treat this effect with a foreground avoidance approach [60,[69][70][71], considering two different methods.The first one consists in applying an exponential suppression factor to the power spectrum in order to remove the radial modes k ∥ = µ k smaller than a critical scale k ∥fg .This can be modelled as The second approach is more conservative, i.e. all the scales corresponding to [70] k < N fg k min (3.17) are excluded from the analysis, where k min is the minimum wavenumber in the given redshift bin and N fg is a factor to be chosen.As a consequence, the survey's effective volume and the sky fraction must be rescaled accordingly by a factor N 3 fg .Note that this method not only limits the radial k ∥ , but also the transverse k ⊥ .
Both methods result in a loss of signal at large scales, which is problematic for constraining parameters such as f NL .In the following we will consider the exponential suppression as reference for the analysis, and then we will also compare it to the results that can be obtained from the more drastic approach.The resulting Hi power spectrum, once the damping factors are taken into account, reads or, in the conservative approach of Eq. (3.17), where Θ is the Heaviside step function.

Cross-correlation between data sets
The cross-correlation of two independent tracers of the matter distribution leads to the crosspower spectrum, in which the bias with respect to the matter power spectrum is a combination of the linear biases of both the tracers -in our case, galaxies and Hi.Concerning the largest scale achievable in this case, we choose to be conservative and to consider as reference the one of the experiment, which can only access the smallest of the two volumes.Since the observations of spectroscopic galaxies and Hi are independent, the galaxy shot noise and the thermal noise in the Hi power spectrum are uncorrelated; the variance of the cross power spectrum is therefore just given by the cosmic variance term, while the noise term vanishes, viz.
In fact, there is in principle a noise term for the cross-correlation, due to the overlap mass range between the dark matter haloes hosting Hi galaxies and Hα galaxies.However, it has been shown that this term is different from zero only at low redshifts (up to redshift z ∼ 0.5), and in general it is very small [15].Thus, Eq. (3.20) is a reasonable assumption for our analysis.As a further test of the robustness of our analysis in this respect, we re-run it using ] Galaxies: shot noise HI: thermal noise A comparison between the noise power spectrum for the galaxy survey and the Hi intensity mapping survey.Shot noise for the galaxy power spectrum is represented in blue and is computed for flux limits F c = 10 −16 erg s −1 cm −2 (dash-dotted), F c = 2 × 10 −16 erg s −1 cm −2 (solid) and F c = 3 × 10 −16 erg s −1 cm −2 (dotted).Thermal noise, in green, is computed for an Hi survey with the specifics given in Section 3.2.The shot noise is larger than the thermal noise at all redshifts, even when lowering the flux limit to F c = 10 −16 erg s −1 cm −2 .
the Fisher matrix formalism and considering a smaller value of k max , which would mimic the effect of having small-scale power washed out by the additional cross-correlation noise term.We find that, by fixing k max = 0.05 Mpc −1 for all the redshift bins, there is only a 2% worsening of the marginalised uncertainty on f NL with respect to the complete case.
Finally, the damping term at large k ⊥ due to the beam, D b , enters linearly in the crosspower spectrum and, in case Hi foreground are treated with the large scale correction given by Eq. (3.16), this needs to be taken into account also in the cross-correlation via

Detection significance
First of all, we consider what the overall detection significance of the signal will be.We estimate it via the signal-to-noise ratio (SNR), which in the multi-tracer scenario reads where the sum runs over the corresponding k-, µ-, and z-bins, and 'overlap' means that the multi-tracer data vector (including the cross-correlation) is computed considering only the overlapping sky area.It corresponds to the one observed by the experiment with the smallest sky coverage, which in our case is the emission-line galaxy survey, with f sky ≃ 0.36.Instead,  .Per-bin SNR for galaxies (blue), Hi (green), and their multi-tracer (bordeaux).Left panel: All noise terms are included.The SNR of galaxies increases at intermediate redshifts, corresponding to Hβ+Oiii galaxies, which are an advantageous combination of high linear bias and a still high observed galaxy number density.The SNR of the Hi intensity mapping survey is more stable.For the multi-tracer, SNR is driven by the Hi term.Right panel: The cosmic variance limited case (zero noise).Breaks in the monotonic trend are related to the spikes in k min due to the transition between different tracers (i.e. the impossibility to reach larger scales).
for the auto-correlations of galaxies and Hi and their cross-correlation alone, the expression above reduces to Since we are considering two surveys with non-perfectly matching sky area, in order not to throw away any data, we include in the full multi-tracer SNR also the contribution of the auto-correlation of the single tracers from the non-overlapping regions [72], namely This operation lets us retain all the information from the survey with the largest sky coverage, that would be otherwise lost, and it is very advantageous for observation at high redshift and large scales.
Recalling the expression of the covariance matrix associated to the power spectrum, Eq. (2.5), we can appreciate how the SNR depends on the combination of the effects of the noise term and the number of independent modes, which are different in the case of a galaxy redshift survey or an intensity mapping observation.The shot noise of the galaxy sample and the thermal noise of the Hi survey are both scale independent and they increase with redshift.As it can be seen in Fig. 3, the shot noise increases much more rapidly than the thermal noise.
Then, the number of independent modes depends on both the probed scale and the redshift; it increases with redshift, whereas it drops the larger the scale.It is also linearly related to the observed sky fraction: when the impact of foregrounds on the Hi power spectrum are modelled with the damping factor D fg , we have f sky,g < f sky,Hi ; viceversa, if the scale cut Left panel: SNR for the galaxy auto-correlation at different flux limits, F c = 2 × 10 −16 erg s −1 cm −2 (solid), F c = 3 × 10 −16 erg s −1 cm −2 (dotted), F c = 10 −16 erg s −1 cm −2 (dashdotted).In the last case, since the observed galaxy number density is higher at all redshifts and for all galaxy types, the increase of SNR for the Hβ+Oiii sample compared to the other samples is milder.Right panel: SNR for the Hi auto-correlation, computed for different foreground avoidance methods -the damping factor in Eq. (3.16) (solid), which is not too dissimilar to the ideal case without foregrounds (dash-dotted line); a scale cut with N fg = 2 (dotted line) leading to significantly lower SNR.
approach is used, we have f sky,g > f sky,Hi due to the rescaling of the effective volume and observed area.For the Hi, the beam damping and possibly the foreground damping play a role as well, since they result in a loss in the power of the signal.Tracer bias is another relevant element, since the amplitude of the primordial non-Gaussianity correction depends on the term b A − 1.This means that if the bias is close to unity, we expect the constraining power on f NL to be weaker.But we emphasise that this feature is model dependent, as it comes from the ansatz of Eq. (2.3).
For the benchmark configuration, all these factors together lead to a total SNR that is larger for the Hi signal than for galaxies, mainly because of the different sky coverage and the shot noise having more impact than the thermal noise, as depicted in the left panel of Fig. 4. The cosmic variance limited case is shown in the right panel of Fig. 4, since it helps to understand the impact of the value of the bias on the SNR and the effects of including the noise terms.As expected, the SNR increases almost monotonically with the redshift, thanks to the inclusion of larger and larger scales and the increasing value of the biases, especially for the different galaxy samples at higher redshift.When comparing this result with the cases where noise terms are included, we notice that there is a general reduction of the detection significance, especially at high redshift, where the noise is more intense and galaxies are more affected (the shot noise is much higher than the thermal noise).
Furthermore, it is interesting to analyse the responses of the galaxy samples, whose biases range across a wider interval because of the transition between galaxy types.Hβ+Oiii galaxies are the most robust sample, thanks to the combination of a high bias and an intermediate noise.By contrast, the Hα sample, despite its higher number density, has a bias closer to unity, and the SNR of the Oii sample is strongly suppressed by the shot noise.For these reasons, although the multi-tracer SNR in the cosmic variance limited case receives contributions from galaxies and Hi, in the complete case most of the information comes from the Hi intensity mapping.
Focusing on galaxies, we also compared the galaxy SNR obtained with different flux limits.The point of interest is the possible balance between variations in the shot noise and in the galaxy bias: for larger values of F c , the observed galaxy number density decreases but the galaxy bias is higher, and vice versa.The shot noise is the major source of variations.When considering F c = 10 −16 erg s −1 cm −2 , the SNR for the galaxy survey is almost doubled and becomes comparable with that of the Hi intensity mapping survey.By contrast, for F c = 3 × 10 −16 erg s −1 cm −2 the galaxy SNR drops, with SNR < 50 for z > 2.5.This behaviour corroborates our previous considerations regarding the impact of the noise term, in particular at z ≥ 2.5.The Hα sample partly deviates from this trend and we see a partial mitigation of the effects of the noise.In this case the galaxy number density responds the least to a change in flux limit, and the SNR is much more sensitive to variation in the galaxy bias, since it is close to unity.Overall, the possible gain due to a more moderate shot noise in the case with F c = 10 −16 erg s −1 cm −2 is reduced by having a lower bias, while the worsening expected for a less sensitive survey is alleviated by a higher bias.
On the other hand, concerning Hi intensity mapping, foregrounds are crucial, as can be seen when moving to the more aggressive foreground removal approach and cutting out the largest scales, or to the ideal situation without foregrounds.In the former case, there is a worsening by a factor of ∼ 3, while in the latter case the SNR remains almost unchanged.These results, summarised in Fig. 5, are shown in the full multi-tracer SNR, which is larger than the single-tracer SNRs in all cases.

Analysis
After having assessed the detectability of the signal, the aim of the subsequent analysis is to test to which extent the method is able to detect primordial non-Gaussianity imprints in the tracer power spectrum, estimating the uncertainty on the estimation of f NL .In particular, we focus on comparing the constraining power of the multi-tracer technique with respect to the auto-correlation of a single tracer.We consider different survey configurations, to test their performance and to have a deeper understanding on the impact of different effects, such as the foreground contamination of the 21-cm signal, or the flux limit of a galaxy survey.We analyse the following cases: • Benchmark scenario: Galaxy survey with F c = 2 × 10 −16 erg s −1 cm −2 and impact of foreground modelled with the damping factor D fg .
-Galaxy survey with F c = 1×10 −16 erg s −1 cm −2 and impact of foreground modelled with the damping factor D fg .
• Pessimistic scenarios: -Galaxy survey with F c = 2 × 10 −16 erg s −1 cm −2 and and impact of foreground modelled with the scale cut using N fg = 2.
-Galaxy survey with F c = 3×10 −16 erg s −1 cm −2 and impact of foreground modelled with the damping factor D fg .
The analysis is performed by sampling the likelihood with Markov chain Monte Carlo algorithms using the publicly available package emcee [73] 6 .In analogy with the expression of the signal-to noise-ratio Eq. ( 4.3), we build a Gaussian likelihood function including the multi-tracer combination together with non-overlap single-tracer information, namely ln L(d|θ) MT tot = ln L(d|θ) MT overlap + ln L(d|θ) AA no-overlap + ln L(d|θ) BB no-overlap , ( where d is the theoretical data vector and θ an array containing the free parameters of the model, whose best estimate can be obtained by maximising the likelihood.We consider as free parameters, beside f NL , the primordial spectral index n s and the biases of the tracers, meaning that the terms in Eq. (5.1) providing the single-tracer information do not directly contribute to the reconstruction of the bias of the other tracer; they indeed allow to get globally tighter constraints on the parameters because they do bring relevant information on the other parameters.
We adopt uniform priors on all parameters and assume as fiducial values the Planck 2018 cosmology [74], and the expected values for the biases [50,56].We also make sure that the priors are large enough not to bias our results on the recovered values and the posterior distribution.More quantitatively, our priors are several tens of times broader than the final constraints; only in one case being narrower, but still more than three times the size of the final marginalised errors.We decide to keep the bias and the spectral index as free parameters because they are clearly degenerate, albeit to a different extent, with f NL , while we expect the other cosmological parameters to be mostly constrained on smaller scales, or by independent measurements.Note that the amplitude of the primordial fluctuations, A s is also degenerate with f NL , and its inclusion in the parameter set may lead to a loosening of the constraints on PNG, if we assume the universality of the halo mass function.However, since cosmic microwave background measurements put tight constraints on A s from the power spectrum, whereas their constraining power on f NL is complementary, coming from the bispectrum, we can safely fix A s to the values measured e.g. by Planck.
Hence, our parameter set is θ = {b g,i , b Hi,i , n s , f NL }, with the further subscript 'i' denoting that the nuisance parameters for the bias(es) correspond to the value of the relevant quantity in each redshift bin.Since the biases of the tracers are redshift dependent, for each z-bin they correspond to one free parameter in the case of the auto-correlation power spectrum and two free parameters in the case of the cross-correlation power spectrum.This implies a rapidly increasing number of free parameters and a higher computational cost if one wants to exploit the information from more bins jointly.For this reason we first consider two redshift bins at a time (in total four free parameters for the auto-power spectrum, and six for the cross-power spectrum).This will be presented in Section 6.1.
Then, we also consider bins grouped in order to divide emission-line galaxies of the spectroscopic survey: this allows to improve the amount of information.Moreover the division of the data set according to the transition to different emission-line galaxy types, allow us to evaluate the impact of the galaxy and Hi biases on the uncertainty on f NL .Having assumed the universality relation for the assembly bias b ϕ , the constraining power depends on the amplitude of b A − 1.In particular, we expect a more pronounced dependence in the case of the auto-correlation of a single tracer, its bias entering quadratically the power spectrum, whereas the product of the biases of the two tracers in the cross-correlation would mitigate the effects of a sharp variation in the bias of one between galaxies and Hi.This analysis will be discussed in Section 6.2.
For the benchmark scenario we also calculate the uncertainty on f NL considering the whole redshift range.

Results
First of all, we present the results of the benchmark configuration analysed spanning over the whole redshift range: the posterior distribution is shown in Fig. 6.The recovered parameters are consistent with the fiducial values whether the analysis is done with the auto-correlation of a single tracer or with the multi-tracer technique.The uncertainty on the primordial non-Gaussianity parameter can be then evaluated by marginalising over all the other parameters.The evaluation of all the redshift range leads to the following constraints on f NL : • Galaxies auto-correlation: f NL = 0.0 ± 2.8; • Hi auto-correlation: f NL = 0.0 ± 2.3; • Multi-tracer technique: f NL = −0.01 ± 0.76.
Hi constraints are slightly tighter than those provided by galaxies because the thermal noise is lower than the shot noise and also because, in this setup, the Hi intensity mapping survey would cover a larger cosmic volume and this results to be advantageous even in presence of foregrounds, if their effect is modelled as a damping factor in the 21-cm signal.
The multi-tracer technique allows us to get an uncertainty on f NL that is more than halved with respect to those coming from the simple auto-correlation, despite the bigger amount of free bias parameters (a total of 26 free parameters for the multi-tracer technique, instead of the 14 needed for the auto-correlations alone).Moreover, reaching σ(f NL ) ≃ 1 is crucial: many inflationary models predict f NL ≃ 1, making σ(f NL ) ≃ 1 a relevant threshold to get to in order to precisely constrain this parameter [75,76].

Redshift dependence
By analysing the redshift bins two by two, it is possible to evaluate how the constraining power of this analysis varies as a function of the redshift.The primordial non-Gaussianity parameter fiducial value is always efficiently recovered with an uncertainty that decreases going to higher redshifts, as a consequence of the fact that larger volumes can be accessed at high z.This is shown in the left panel of Fig. 7.
Focusing on the auto-correlation of the galaxy sample, it can be noticed that f NL is not well constrained in the low-z bins, both because very large scales cannot be accessed at those redshifts, and because the bias term (b g − 1) is very small for the galaxies at those redshifts.This, in turn, lowers the effects of f NL of the power spectrum.On the other hand, at high redshift, in the bins corresponding to the Oii sample, the uncertainty on f NL saturates and, after reaching a lower limit, it slightly increases; this is due to the fact that the shot noise becomes too large and it starts to dominate the observed power spectrum.Thus, the advantage of accessing large scales is lost.A test was done considering a configuration in which the galaxy power spectrum is assumed to be cosmic variance limited, thus without shot noise: the constraints obtained on f NL improve at every redshift but in particular at high redshifts.These effect are not present in the auto-correlation of the Hi signal: the higher SNR and the larger sky coverage (despite the foregrounds) ensure a better constraining power for the Hi auto-correlation.Contour plot resulting from the MCMC analysis of the full data set, marginalising over the bias parameters in order to focus on the cosmological ones: blue, green, and bordeaux contours correspond respectively to the posterior distribution of the galaxy auto-correlations, of the Hi auto-correlation, and of the multi-tracer technique (this color code will be used in the next plots as well).The multi-tracer technique provides tighter constraints, and, thanks to the combination of independent biased tracers, the degeneracy between f NL and n s is strongly mitigated.
The multi-tracer technique results in be the best performing method to constrain primordial non-Gaussianity.It is not affected by the low SNR at high z of the power spectrum of galaxies, nor by its limited scale range at low z, nor by the loss of information due to foregrounds contaminating the Hi signal.The improvement on the constraints is up to 30% with respect to the best performing tracer in the corresponding bins.The gain is lower when one of the two tracers is individually performing significantly better: in this case that tracer gives the dominant contribution in the multi-tracer technique.
The constraining power on the others free parameters of the multivariate analysis is also taken into account.The primordial spectral index and the bias parameters are all recovered within a confidence level of 68% in both the auto-correlation of the single tracers and the multi-tracer technique, even in the redshift bins where galaxies cannot constrain f NL .When applying the multi-tracer technique, the uncertainty associated to these parameters is reduced as well, mostly as a consequence of reduced degeneracies between parameters.Left: 68% error bars on f NL from the analysis of two redshift bins at a time.Horizontal bars represent the total width of two bins considered.Right: Constraints from the analysis with redshift bins grouped in order to follow the division in galaxy types.For comparison, the constraint from the cosmic microwave background [5] is shown by the dashed grey line.

Tracer dependence
The analysis of more redshift bins at a time leads to a global improvement of the results, as displayed in the right panel of Fig. 7: the constraints on f NL are reduced up to a factor ∼ 2 thanks to the improvement of the amount of information, and there is a global gain in the performance also concerning the other free parameters.As in the analysis performed per redshift bins, f NL is recovered with an uncertainty decreasing with increasing redshift.
In the present analysis, the differences between the various tracers are highlighted, in particular concerning the bias term.The Hα galaxies still suffer from the scale limit imposed by the small surveyable volume at low redshift, as well as from the small bias.As a consequence, the constraints from their auto-correlation are less competitive.The best performing sample is the Hβ+Oiii sample, thanks to a combination of three advantageous conditions: first of all the effects of a primordial non-Gaussianity are enhanced by a large bias; secondly, it is still possible to observe a significant number of galaxies; and, lastly, the scale range reachable at those redshifts is extended thanks to much larger volume.For the Oii sample only two of those conditions are met: the large bias and the progressive enlargement of the scale range, however this sample is affected by its low density, which prevents the constraints to improve with respect to the previous sample.
Since in the case of the Hi intensity mapping, there is no transition between different samples, the considerations on the effects of the bias are less relevant because the trend of b Hi is smooth and moreover b Hi significantly deviates from unity even at low z.Thus the results are mainly driven by the accessible scale range, which is wider at high z, so that f NL is recovered with decreasing uncertainty at high redshift.
All such considerations remain valid when the multi-tracer technique is applied, but the results are more solid against possible degrading effects such as the galaxy shot noise or the small bias term.The gain with respect to the simple auto-correlation is comparable to what can be obtained with the analysis with only two redshift bins: to get a larger improvement more information should be added analysing all the redshift range available.

Flux dependence
By varying the flux limit of the galaxy survey, we can appreciate the impact of the shot noise on the final result.
The variation of the observed number density with flux cut is more evident for the oxygen-line samples, namely the Hβ+Oiii sample and the Oii sample; and especially for this last case, since it is the one at the highest redshift, the detectability of the sources is therefore very responsive to any change in the specifics of the survey.When varying the flux limit, the galaxy bias needs to be evaluated according to it.
The results are shown in Fig. 8.In general, they do not change significantly at low z, where constraining power is limited by having a very small term (b g − 1) and a small volume.The main differences arise at high redshift both in the analysis with two or four redshift bins, mostly because of the variation of the shot noise.
Considering the optimistic scenario of a survey with F c = 1 × 10 −16 erg s −1 cm −2 , all the parameters are well constrained and focusing on f NL the uncertainty decreases with respect to the benchmark configuration presented in Sections 6.1 and 6.2.The galaxy auto-correlation leads to constraints which are up to a factor of three tighter and the multi-tracer technique as well takes advantage from this improvement.Notably, the uncertainty on f NL does not saturate anymore at the redshifts corresponding to the Oii sample, which is therefore able to provide stronger constraints.
Moving to the more stringent scenario with F c = 3×10 −16 erg s −1 cm −2 , the performance of the analysis resents the higher shot noise, which for the Oii sample is always above the cosmological signal.This implies that the posterior distribution for f NL is prior dominated for z ≳ 2.7 in the case of the galaxy auto-correlation, and it is driven by the Hi signal in the case of the multi-tracer approach.The degradation of the constraints on f NL is less severe at intermediate z, corresponding to the Hβ+Oiii galaxies: despite the higher value of F c , this sample is still dense enough and the higher bias consents to precisely recover the primordial non-Gaussianity parameter.
For the SNR, the combination of the variation of the bias and of the shot noise deserves deeper investigation.It can be noticed that the balancing between the two competing terms becomes more effective for the Hα sample.In particular, there is full compensation in the first redshift bins.Finally, it is worth highlighting that, even if the flux limit has a great impact on the performance of the galaxy auto-correlation, the multi-tracer approach is less sensitive to the variation of the flux limit of the survey.

Foregrounds
Foregrounds are one of the main challenges in the detection of the Hi cosmological signal and the evaluation of their impact is a key point when constraining cosmological parameters.The consequences of the loss of signal due to foregrounds can be estimated considering the ideal case in which the Hi signal is foreground-free.In this case, with the analysis of Hi in auto-correlation the cosmological parameters and the bias parameters are recovered to their fiducial values with higher accuracy because there is no damping in the signal at large scales; the mean improvement is around 15% with respect to our standard scenario, with the biggest gain at low redshift.As a consequence of the larger amount of information especially at large scales, the multi-tracer technique provides tighter constraints on the primordial non-Gaussianity parameter as well.
On the other side, when adopting a more strict approach to avoid foregrounds, the loss of signal can be critical.If one needs to exclude all the modes below a given threshold, as  (left panels) and for a pessimistic scenario with F c = 3 × 10 −16 erg s −1 cm −2 (right panels).Top panels refer to the analysis discussed in Section 6.1 (as in Fig. 7 left).Bottom panels refer to the analysis discussed in Section 6.2 (as in Fig. 7 right).Foregrounds in the Hi intensity mapping are treated with the benchmark approach with the damping factor described in Eq. (3.16), so that the uncertainty on f NL from the Hi auto-correlation is the same as in Fig. 7.
explained in the last part of Section 3.2, the analysis suffers from a lack of information on those scales which are more sensitive to the presence of primordial non-Gaussianity.The worsening of the results is significant at all z and in the case of the auto-correlation of Hi the uncertainty on f NL is at least doubled with respect to the reference configuration.Thanks to the combination of galaxies and Hi in the multi-tracer technique the degradation of the constraints is mitigated.

Conclusion
In the next years cosmological surveys aiming to map the large-scale structure of the Universe will shed light on the characteristics of the primordial Universe, and thanks to the development of the Hi intensity mapping technique, the complementarity of two observables -galaxies in the optical band and Hi in the radio band -will be a powerful tool to exploit.  .Impact of foreground avoidance on the uncertainty on f NL .Left panels show forecasts without foregrounds, while right panels include a scale cut in the Hi intensity mapping.Top panels refer to the analysis discussed in Section 6.1 (as in Fig. 7 left), while bottom panels refer to the analysis discussed in Section 6.2 (as in Fig. 7 right).The flux limit of the galaxy survey is kept constant at F c = 3 × 10 −16 erg s −1 cm −2 , so the results for galaxy auto-correlation are the same as in Fig. 7.
In this context, we presented a forecast on the performance in measurements of the primordial non-Gaussianity parameter, f NL , done in a Bayesian framework, to assess the gain on the constraining power that can be reached when applying the multi-tracer technique.The scale dependent bias induced by the presence of a primordial non-Gaussianity leads to a change in the slope of the power spectrum of a matter tracer, mostly evident at large scales, where a multi-tracer approach is expected to grant the best performance thanks to the possibility of beating the cosmic variance.
Starting from a fiducial cosmology generated with CAMB, we simulated a synthetic data set for an emission-line galaxy survey, for which we also took into account interloper samples at high z in order to extend the redshift coverage.Concerning the Hi intensity maps, we simulated the signal including instrumental effects, such as the beam damping, and the loss of signal at large scales due to astrophysical foregrounds.
By means of MCMC methods, we evaluated the constraining power of the analysis in terms of the accuracy on the reconstruction of our free parameters (f NL , n s , and the biases of the tracers), focusing on the primordial non-Gaussianity parameter.The full multi-tracer likelihood was build including not only the data from the overlapping region between the galaxy survey and the Hi intensity mapping survey, but also what the individual tracers provide from non-overlapping regions, allowing us to retain all the available information, especially at large scales.
We found that this method is unbiased and that all the parameters are recovered to their fiducial value both with the single-tracer and the multi-tracer approach.However, when performing the analysis on the total data set, the multi-tracer technique provides an improvement of a factor greater than 2 on the uncertainty associated to the reconstructed parameters and especially on f NL , showing that it will be advantageous to go for a multitracer approach and therefore design surveys in different bands with maximised overlapping coverage both in redshift and in the observed sky patch.
When considering the results with respect to the redshift, it is evident that the constraining power improves when going to high redshift, when larger cosmic volumes are accessible, as long as the noise term does not exceed the cosmological signal.This is the case of the last redshift bins of the galaxy sample, where the observed number density of Oii interlopers is insufficient to provide tight constraints on f NL .This does not mean however that this sample must be excluded, because the multi-tracer technique results to be quite robust against the increase of the noise of one of the samples considered.Beside the effect of the scale range, we assessed the impact of the amplitude of the bias, which can be mostly appreciated from the analysis of the galaxy auto-correlation: the Hα sample, despite being the main sample of the modelled galaxy survey, preforms the worst even compared with the interloper samples precisely because of the values of the bias.This justifies again the need of exploiting also other samples beside the target one.
A caveat must be added in the discussion of the bias, since we used the universality relation for b ϕ .Recent literature demonstrates that this might not be the correct approach, and that further investigation of the assembly bias is needed [18,20,28,77].Identifying a better model for b ϕ would allow us to achieve more accurate results on primordial non-Gaussianity.However, the multi-tracer technique provides some mitigation here, since it is more robust to inaccurate models of the assembly bias [20].
We evaluated the impact of the shot noise associated to the galaxy power spectrum by varying the flux limit of the survey, finding that the main differences arise at high redshift and in the analysis of the galaxies in auto-correlation, while the multi-tracer technique is more stable against variation in the observed galaxy number density.Concerning the Hi intensity mapping data, foregrounds were taken into account with different approaches in order to estimate how much the constraints on f NL are affected by how the foreground avoidance is modelled.We found that there is a high variability in the results depending on how severe we model the loss of signal to be.The development of techniques to clean the Hi intensity mapping data from foregrounds without erasing the cosmological signal is therefore crucial: blind and non-parametric methods are typically employed [62,68,78,79] and it will be important to explore in detail the best settings, eventually in combination with techniques to reconstruct the Hi signal.
Globally, the constraints we found from our forecast on the presence of a primordial non-Gaussianity are tight, reaching their best value thanks to the multi-tracer technique.These are promising results, suggesting that next-generation survey both in the optical and the radio bands will allow to probe inflationary models and the dynamics of the early Universe with high accuracy.In order to do this, we should take advantage of the multi-tracer technique to overcome the cosmic variance.6, the widening of the galaxy auto-correlation contours is the most evident difference, whereas there is only a slight difference in the multi-tracer contours.The multi-tracer is confirmed to be more constraining than both single tracers in autocorrelation and also to break the degeneracy between f NL and n s .
while keeping fixed the Hα number density, i.e., As depicted in the left panel of Fig. 12, this leads to a gradual broadening of the constraints on f NL obtained with the galaxy auto-correlation, up to a ∼ 50% worsening in the case that only 25% of interlopers is detected.On the other hand, the multi-tracer uncertainty on f NL is at most only 10% larger than in the benchmark scenario.This confirms the stability of the multi-tracer technique even against a degradation of the quality of the data sets, thanks to the combination of more than one tracer -also in agreement with the results in Section 6.3.The strength of the features due to primordial non-Gaussianity is related to the value of the bias through the b A − 1 term entering the expression of the scale-dependent bias.This means that any inaccuracy in the model of the bias could be reflected in the constraints on f NL .Among the different galaxy samples and Hi, the interloper samples are again the ones whose modeling is the most difficult to describe, so we explored a different parametrisation of the bias of interloper galaxies at high redshift.We adopted the parametrisation of [48] for the intermediate sample (and for consistency also for the Hα sample) and, following [80], we Solid lines correspond to the fitting formula of [50], for a flux limit F c = 2 × 10 −16 erg s −1 cm −2 .Ddashed lines correspond to the models in [48] (Hα and Oiii galaxies) and in [80] (Oii galaxies).used a fixed bias of b g = 2.5 for the Oii sample, as shown in Fig. 11.The main goal of this analysis is to compare the results from the bias model of [50], which shows a steep evolution in redshift, and from other models, where the galaxy bias assumes smaller values even at high z.Hence, even a fixed value for the Oii galaxies can at first approximation be used to highlight this difference.We also emphasise that the fitting formula in [48] is obtained considering a flux limit of F c = 10 −16 erg s −1 cm −2 , which is lower than the one used in this analysis.Since a lower flux limit also implies a lower cumulative bias, these values are underestimated and therefore we expect to find pessimistic constraints from the artificially decreased term b A − 1.
Moreover, since there is a large uncertainty on the properties of Oii galaxies, we also performed the analysis excluding these galaxies.In this case we discarded the four corresponding redshift bins not only for the galaxy auto-correlation but also from the Hi data set and in the multi-tracer analysis.This scenario should provide the most conservative results, since the multi-tracer likelihood also lacks the Hi single-tracer information from the non-overlapping redshift interval.
The results of this forecast are shown in the right panel of Fig. 12 and are consistent with the previous tests: the galaxy auto-correlation is more prone to variations of the results while the multi-tracer is more stable.It is also interesting to note that a lower galaxy bias in the multi-tracer analysis impacts much less than cutting the redshift range at z = 2.6, since a cut in redshift affects the full data vector and not only one of the two tracers.In this regard, the different response to this cut in the single-tracer and multi-tracer cases is interesting.The Hi auto-correlation constraints respond the most, while the galaxy auto-correlation is less affected.This is in line with the overall results of the analysis, where we find the strongest constraining power for Hi at the highest redshift, where galaxies are already limited by the shot noise (see Fig. 7).The loss of signal from this redshift range is much more damaging Left panel: 68% error bars on f NL for f int from 0.25 (a quarter of the observable interlopers detected) to 1 (all detected).The blue line (galaxy auto-correlation) is much steeper the bordeaux line (multi-tracer): a depletion of the interloper samples downgrades the galaxy single-tracer results, but does not significantly affect the multi-tracer results.The green line (Hi autocorrelation) is shown as reference.Right panel: 68% error bars on f NL obtained with different bias models -considering the whole redshift range (solid lines) or excluding the last redshift bins (dashed lines).The galaxy auto-correlation is more sensitive to variation of the value of the bias, whereas the multi-tracer constraints show only a few percent change.Note that removing the redshifts above z = 2.6 has a different impact on the Hi and galaxy single-tracer and on the multi-tracer analysis, but the overall trends are preserved.
for the Hi single-tracer analysis.The impact on the multi-tracer method is, as expected, a balance between the one on the two auto-correlations.Furthermore, the variation of the galaxy bias has almost the same effects using the full z-range or using only the redshifts z < 2. 6.
In all the considered scenarios, the multi-tracer technique is confirmed to be the most powerful in constraining f NL .

B Fisher matrix formalism
In analogy to the expression of the signal-to-noise ratio and of the likelihood defined in Sections 4 and 5, the Fisher matrix is a

Figure 4
Figure 4. Per-bin SNR for galaxies (blue), Hi (green), and their multi-tracer (bordeaux).Left panel: All noise terms are included.The SNR of galaxies increases at intermediate redshifts, corresponding to Hβ+Oiii galaxies, which are an advantageous combination of high linear bias and a still high observed galaxy number density.The SNR of the Hi intensity mapping survey is more stable.For the multi-tracer, SNR is driven by the Hi term.Right panel: The cosmic variance limited case (zero noise).Breaks in the monotonic trend are related to the spikes in k min due to the transition between different tracers (i.e. the impossibility to reach larger scales).

1 HI: N fg = 2 Figure 5 .
Figure 5.Left panel: SNR for the galaxy auto-correlation at different flux limits, F c = 2 × 10 −16 erg s −1 cm −2 (solid), F c = 3 × 10 −16 erg s −1 cm −2 (dotted), F c = 10 −16 erg s −1 cm −2 (dashdotted).In the last case, since the observed galaxy number density is higher at all redshifts and for all galaxy types, the increase of SNR for the Hβ+Oiii sample compared to the other samples is milder.Right panel: SNR for the Hi auto-correlation, computed for different foreground avoidance methods -the damping factor in Eq. (3.16) (solid), which is not too dissimilar to the ideal case without foregrounds (dash-dotted line); a scale cut with N fg = 2 (dotted line) leading to significantly lower SNR.
Figure 6.Contour plot resulting from the MCMC analysis of the full data set, marginalising over the bias parameters in order to focus on the cosmological ones: blue, green, and bordeaux contours correspond respectively to the posterior distribution of the galaxy auto-correlations, of the Hi auto-correlation, and of the multi-tracer technique (this color code will be used in the next plots as well).The multi-tracer technique provides tighter constraints, and, thanks to the combination of independent biased tracers, the degeneracy between f NL and n s is strongly mitigated.

Figure 7 .
Figure 7. Left: 68% error bars on f NL from the analysis of two redshift bins at a time.Horizontal bars represent the total width of two bins considered.Right: Constraints from the analysis with redshift bins grouped in order to follow the division in galaxy types.For comparison, the constraint from the cosmic microwave background[5] is shown by the dashed grey line.

Figure 8 .
Figure 8. Marginalised uncertainty on f NL for an optimistic scenario with F c = 10 −16 erg s −1 cm −2 (left panels) and for a pessimistic scenario with F c = 3 × 10 −16 erg s −1 cm −2 (right panels).Top panels refer to the analysis discussed in Section 6.1 (as in Fig.7 left).Bottom panels refer to the analysis discussed in Section 6.2 (as in Fig.7 right).Foregrounds in the Hi intensity mapping are treated with the benchmark approach with the damping factor described in Eq.(3.16), so that the uncertainty on f NL from the Hi auto-correlation is the same as in Fig.7.

Figure 9
Figure 9. Impact of foreground avoidance on the uncertainty on f NL .Left panels show forecasts without foregrounds, while right panels include a scale cut in the Hi intensity mapping.Top panels refer to the analysis discussed in Section 6.1 (as in Fig.7left), while bottom panels refer to the analysis discussed in Section 6.2 (as in Fig.7 right).The flux limit of the galaxy survey is kept constant at F c = 3 × 10 −16 erg s −1 cm −2 , so the results for galaxy auto-correlation are the same as in Fig.7.

Figure 10 .
Figure 10.Contour plot resulting from the MCMC analysis of the Roman-SKAO data set, marginalising over the bias parameters.If compared with Fig.6, the widening of the galaxy auto-correlation contours is the most evident difference, whereas there is only a slight difference in the multi-tracer contours.The multi-tracer is confirmed to be more constraining than both single tracers in autocorrelation and also to break the degeneracy between f NL and n s .

Figure 11 .
Figure 11.Galaxy bias computed with various parametrisations for the different galaxy types.Solid lines correspond to the fitting formula of[50], for a flux limit F c = 2 × 10 −16 erg s −1 cm −2 .Ddashed lines correspond to the models in[48] (Hα and Oiii galaxies) and in[80] (Oii galaxies).
Figure 12.Left panel: 68% error bars on f NL for f int from 0.25 (a quarter of the observable interlopers detected) to 1 (all detected).The blue line (galaxy auto-correlation) is much steeper the bordeaux line (multi-tracer): a depletion of the interloper samples downgrades the galaxy single-tracer results, but does not significantly affect the multi-tracer results.The green line (Hi autocorrelation) is shown as reference.Right panel: 68% error bars on f NL obtained with different bias models -considering the whole redshift range (solid lines) or excluding the last redshift bins (dashed lines).The galaxy auto-correlation is more sensitive to variation of the value of the bias, whereas the multi-tracer constraints show only a few percent change.Note that removing the redshifts above z = 2.6 has a different impact on the Hi and galaxy single-tracer and on the multi-tracer analysis, but the overall trends are preserved.

F∂
MT  overlap αβ = F MT overlap αβ + F AA no−overlap αβ + F BB no−overlap αβ .(B.1)Here derivatives of the power spectrum are computed with respect to the parameter vectorθ α = {f NL , n s , b g (z), b Hi (z)}.The multi-tracer overlap contribution and the single-tracer non-overlap contribution areF MT overlap αβ = k,µ,z ∂ α P T Cov −1 (P , P ) ∂ β P , α P AB ∂ β P AB [∆P AB ] 2 .(B.3)The marginalised uncertainty on each free parameter of the model can be then evaluated as σ (θ α ) = F

Table 1 .
Values of the parameters of the fitting formula used for the galaxy bias.