Constraining Primordial non-Gaussianity with Bispectrum and Power Spectum from Upcoming Optical and Radio Surveys

We forecast constraints on primordial non-Gaussianity (PNG) and bias parameters from measurements of galaxy power spectrum and bispectrum in future radio continuum and optical surveys. In the galaxy bispectrum, we consider a comprehensive list of effects, including the bias expansion for non-Gaussian initial conditions up to second order, redshift space distortions, redshift uncertainties and theoretical errors. These effects are all combined in a single PNG forecast for the first time. Moreover, we improve the bispectrum modelling over previous forecasts, by accounting for trispectrum contributions. All effects have an impact on final predicted bounds, which varies with the type of survey. We find that the bispectrum can lead to improvements up to a factor $\sim 5$ over bounds based on the power spectrum alone, leading to significantly better constraints for local-type PNG, with respect to current limits from \textit{Planck}. Future radio and photometric surveys could obtain a measurement error of $\sigma(f_{\mathrm{NL}}^{\mathrm{loc}}) \approx 0.2$. In the case of equilateral PNG, galaxy bispectrum can improve upon present bounds only if significant improvements in the redshift determinations of future, large volume, photometric or radio surveys could be achieved. For orthogonal non-Gaussianity, expected constraints are generally comparable to current ones.


INTRODUCTION
So far, cosmological analyses of Large Scale Structure (LSS) surveys have relied nearly exclusively on matter and galaxy power spectrum estimation.It is however wellknown that important extra-information can be extracted via higher-order correlation functions, such as the matter and galaxy bispectrum (three-point correlation function of Fourier modes), which allow both probing the non-linear regime of structure growth and setting constraints on primordial non-Gaussianity (PNG) (see e.g.Bernardeau et al. 2002, and references therein).The latter is in particular an important and general prediction of inflationary theories.It is a direct product of non-linear interactions during the infla-E-mail: dionysios.karagiannis@pd.infn.it† Marie Sk lodowska-Curie fellow tionary or reheating stage.The bispectrum of the primordial curvature perturbation field, arising from such interactions, is characterized by a dimensionless, amplitude parameter, fNL, and by a shape function F (k1, k2, k3).While fNL defines the strength of the PNG signal, the shape describes the functional dependence of the bispectrum on different Fourier space triangles.Both of them are strongly model dependent and provide significant information on the physical mechanisms at work during inflation.Three very important shapes, encompassing a large amount of scenarios, are the so called local shape (NG signal peaking on "squeezed triangles", i.e. k1 k2 ∼ k3), equilateral shape (k1 ∼ k2 ∼ k3) and foldedshape (k1 ∼ k2 ∼ k3/2).
Currently, the tightest experimental fNL bounds, including a large number of different shapes, come from Planck Cosmic Microwave Background (CMB) bispectrum measurements (Planck Collaboration et al. 2016b).Bispec-trum measurements of LSS data have been already obtained (Scoccimarro et al. 2001b;Feldman et al. 2001;Verde et al. 2002;Marín et al. 2013;Gil-Marín et al. 2014, 2017), but the current level of sensitivity is not sufficient to generate useful PNG bounds (current LSS power spectrum constraints on local fNL are more interesting (Padmanabhan et al. 2007;Slosar et al. 2008;Xia et al. 2010bXia et al. , 2011;;Nikoloudakis et al. 2013;Agarwal et al. 2014;Karagiannis et al. 2014;Leistedt et al. 2014), albeit still not competitive with the CMB (Verde et al. 2000a)).On the other hand, bispectrum estimates of fNL with future LSS data do have in principle great potential to improve over CMB bounds, at least for specific shapes.This is because 3D LSS surveys, covering large volumes and probing a wide range of scales, have access to a much larger amount of modes, with respect to 2D CMB maps.However, LSS measurements will also be very challenging, due to late-time non-linearities, expected to produce much larger non-Gaussian signatures than the primordial component.These contributions therefore need to be understood and subtracted with exquisite accuracy.
The issue of theoretical modelling of non-linear effects and of higher order LSS correlators has indeed been long debated in the literature (see e.g.Bernardeau et al. 2002, for a review) and the interest in producing accurate and realistic LSS primordial bispectrum forecasts has been steadily increasing in recent years.Important contributions in this direction include the work of Scoccimarro et al. (2004) and Sefusatti & Komatsu (2007) -where the bispectrum of galaxies is used for the first time to forecast the constraining power of LSS surveys on measuring the amplitude of PNG -and the study of Song et al. (2015), where information from power spectrum and bispectrum of galaxies is combined -also including redshift space distortion effects -in order to constrain growth parameters and galaxy bias terms.Additional contributions were then made by Tellarini et al. (2016), who took into account the second order tidal bias term (McDonald & Roy 2009;Baldauf et al. 2012;Chan et al. 2012), as well as the bivariate bias expansion (Giannantonio & Porciani 2010) in the redshift space galaxy bispectrum, in order to constrain the amplitude of local PNG.Finally, the authors of Baldauf et al. (2016) pointed out the importance of including uncertainties in the theoretical modelling of the signal (theoretical errors) and properly propagating them into the final error bar estimates.
Many more details on these issues -including a more detailed description of improvements and refinements in bias and redshift space distortion modelling -are provided in sections 2.2 and 2.3.2.Here we point out that many of the analysis ingredients mentioned above were considered separately and independently in previous forecasts, with different works considering the importance of specific new terms, without accounting however for all of them at once (for example, theoretical errors are studied in detail in the real space treatment of Baldauf et al. (2016), whereas redshift space distortions are accounted for in detail in Tellarini et al. (2016), without including theoretical errors).Here, for the first time, we consistently include all these terms and produce as complete and realistic as possible PNG forecasts, in terms of fNL parameters, combining power spectrum and bispectrum constraints.One advantage of using the bispectrum is that it opens the possibility to explore the full range of primordial shapes, including the equilateral and orthog-onal ones.These shapes are very little explored in previous PNG LSS studies.In addition, we also include -for the first time in an actual forecast -trispectrum contributions to fNL arising from the bias expansion in the galaxy bispectrum, which were pointed out as potentially important (Sefusatti 2009;Jeong & Komatsu 2009).We note that such term could play a significant role in constraining the signal from non-local shapes.
Another important issue in a LSS PNG analysis is of course that of establishing which survey design and which statistical probe provide the best fNL constraints.Clearly, adding modes by going to smaller scales does in principle improve sensitivity.The obvious caveat is that such approach requires non-linear scales, where the non-PNG contribution gets very large and hard to model.Moreover, late time non-linearities couple different modes.This unavoidably produces a saturation of the available information.To estimate in detail this effect, a full calculation of the bispectrum covariance is needed in the evaluation of the signal-to-noise ratio for fNL.All these issues are still open, and they are currently under a significant amount of scrutiny in the literature (Heavens et al. 1998;Crocce & Scoccimarro 2006b,a;Pietroni 2008;Bernardeau et al. 2008;Wagner et al. 2010;Crocce et al. 2012;Baumann et al. 2012;Carrasco et al. 2012;Gil-Marín et al. 2012, 2014;Lazanu et al. 2016Lazanu et al. , 2017;;Lazanu 2017).
In this work we consider the alternative approach: we only investigate large and quasi-linear scales; we look at galaxy clustering statistics at high redshift, where nonlinearities become important at much smaller scales.Forthcoming radio continuum surveys seem ideally suited to this purpose, as they will probe wide areas of the sky over a large redshift range.While the power spectrum of radio continuum has been already considered in the literature (Xia et al. 2010a(Xia et al. ,c, 2011;;Raccanelli et al. 2012;Camera et al. 2015a;Raccanelli et al. 2017), the bispectrum of radio continuum surveys has not been studied so far.We devote particular attention to it in this work.The drawback with radio continuum sources is the lack of a direct determination of their redshifts.Our analysis therefore considers the possibility to extract redshift information via clustering-based estimation methods (Ménard et al. 2013).We follow the implementation for forthcoming radio surveys developed in Kovetz et al. (2017).For comparison, we also present forecasts for optical surveys, both spectroscopic and photometric.
The paper is structured as follows: In Sec.2.1, we provide a brief review of the theory behind matter clustering and PNG; in Sec.2.2 we discuss in more detail the relationship between matter and galaxy statistics; in Sec.2.3 we present the expressions for galaxy power spectra and bispectra considering all the relevant effects; in Sec. 3 we describe the survey specifications that we have used for our analysis; in Sec.4.1 we describe the Fisher matrix forecast method; in Sec.4.2 we explain how to account for theoretical errors in the forecast.We summarize our theoretical framework in Sec. 5.The main results are then shown in Sec. 6.Finally, we conclude in Sec. 7. In the Appendix we present a summary of standard perturbation theory (Appendix A), as well as a small review on the results of the MPTbreeze formalism (Appendix B).In addition, we present the derivation of the bias parameters from the adopted peak-background split method (Appendix C).Finally in Appendix D, we present the kernels used for the redshift space galaxy models.

Matter power spectrum and bispectrum
In linear theory, the matter power spectrum can be related to the primordial gravitational potential power spectrum, via with Here k ≡ |k| is the magnitude of the wavevector k and D(z) is the linear growth factor, originating from the growing mode of the linear fluid equations, where it is normalized to unity today (D(0) = 1).T (k) is the matter transfer function normalized to unity on large scales k → 0 and c is the speed of light.In this work, matter non-linearities will be accounted for in the framework of Standard Perturbation Theory (SPT) or Renormalized Perturbation Theory (RPT).Useful results in SPT are reviewed in Appendix A. To be consistent with the SPT approach we will only use the linear power spectrum and the tree-level bispectrum in SPT, since we will restrict our analysis to only mildly non-linear scales (here we consider kmax(z) = 0.1/D(z) h/Mpc).There is nevertheless additional information encoded in small, nonlinear scales, but modelling this regime is challenging, and will not be considered here.
Based on the inflationary model considered, primordial gravitational potential bispectra are produced in different shapes.In this work we consider the three most widely used shapes, the local (Salopek & Bond 1990;Gangui et al. 1994;Verde et al. 2000a;Komatsu & Spergel 2001), equilateral (Creminelli et al. 2006) and orthogonal (Senatore et al. 2010) types of PNG.
The three primordial bispectrum shapes, defining PNG of the local, equilateral and orthogonal types, and are defined respectively as B loc φ (k1, k2, k3) = 2f loc NL [P φ (k1)P φ (k2) + 2 perms] , (4) B eq φ (k1, k2, k3) = 6f equil NL {−[P φ (k1)P φ (k2) + 2 perms] − 2[P φ (k1)P φ (k2)P φ (k3)] 2/3 + [P Large local PNG can be produced in multi-field scenarios, while sizeable equilateral and orthogonal bispectra can be obtained, for example, in single-field scenarios with nonstandard kinetic terms and higher derivative interactions.The folded shape, mentioned in the introduction, can arise in models with excited initial states and is obtained as a linear combination of equilateral and orthogonal shapes.In fact the orthogonal shape, while difficult to visualise, is by construction orthogonal to the equilateral ones.For more details see e.g.Planck Collaboration et al. (2016b), and references therein.

Bias
The bias of tracers i.e. the relation between the statistics of observed objects and the underlying distribution of matter, can be understood as the combination of two components: the bias of dark matter halos and how the selected tracers populate the halos.Here we first present our adopted modelling for the bias of halos and then in Sec.2.2.3 we use a halo occupation distribution model to connect this to the bias of selected tracers (see e.g.Desjacques et al. 2016, for a review).

Halo bias
In order to use the tree-level galaxy bispectrum to forecast the amplitude of PNG coming from future LSS surveys, knowledge of the bias expansion up to second order is needed.A complete set of the bias terms was presented in Assassi et al. (2014), Senatore (2015) and Mirbabayi et al. (2015), where the notion of local-in-matter bias relation (Coles 1993;Fry & Gaztanaga 1993) is generalized to a local bias perturbative expansion over operators O representing all possible local gravitational observables along the fluid trajectory.Such operators are built from powers of the tidal tensor ∂i∂jΦ which can be decomposed into the trace ∇ 2 Φ ∝ δ and the traceless part sij, i.e. the tidal field (Catelan et al. 2000;McDonald & Roy 2009;Chan et al. 2012;Baldauf et al. 2012).
Higher-derivative terms must also be taken into account in the general local bias expansion (Desjacques et al. 2016).They are constructed from powers of spatial derivatives over each operator O (i.e. higher than two derivatives on Φ) and incorporate non-gravitational contributions in the halo formation process (e.g.gas heating, feedback processes).These terms introduce a limiting spatial scale, into the bias expansion, R * , at which they become important.This fundamental cut-off indicates the scale where a perturbative description of bias breaks down (see e.g.Desjacques et al. 2016, for a discussion).In the case of dark matter halos it is proportional to the Lagrangian radius R(M ) = (3M/4πρ m ) 1/3 .Including all the above, the most general expansion up to second order in the Eulerian framework with Gaussian initial conditions is (Coles 1993;Fry & Gaztanaga 1993;Fry 1996;Catelan et al. 1998Catelan et al. , 2000;;McDonald & Roy 2009;Elia et al. 2011;Chan et al. 2012;Baldauf et al. 2012;Assassi et al. 2014;Senatore 2015;Mirbabayi et al. 2015) where τ is the conformal time and x are the spatial comoving coordinates in the Eulerian frame.In addition s 2 = sijs ij is the simplest scalar that can be formed from the tidal field, ε E is the leading stochastic bias contribution (Dekel & Lahav 1999;Taruya & Soda 1999;Matsubara 1999) and ε E δ is the stochastic counterpart of the linear bias.The secondorder tidal field bias term, following the convention of Baldauf et al. (2012), is given by b Note that this expression is not expected to be exact, since it assumes a formation time of τ * → 0, as well as a conserved evolution (see e.g.Desjacques et al. 2016, and references therein for a detailed discussion).
Although the subscript R is missing from the fields in the bias expansion, they are all assumed smoothed at a radius R > R * .In the above expansion we have ignored all higher-derivative terms (∼ O(R 2 * )) up to second order, since we consider sufficiently large scales (k 1/R * ) wherein their contribution is suppressed by (R * k) 2 .Further, on large scales we can also ignore the velocity bias (see e.g.Desjacques et al. 2016).
The bias coefficients of the general bias expansion can be defined on large scales in the context of the peakbackground split (Kaiser 1984;Bardeen et al. 1986;Cole & Kaiser 1989).Here the mass function proposed by Sheth & Tormen (1999) (ST hereafter) will be used to derive the Eulerian local-in-matter bias terms as a function of the peak height (Scoccimarro et al. 2001a).Details on the peakbackground split argument can be found in Appendix C.

Non-Gaussian effects
In the presence of PNG of the local type, longwavelength fluctuations of the primordial gravitational potential modulate the initial condition perturbations of small scales, due to the induced coupling between the two.This leads to to a scale-dependent bias correction on large scales (Dalal et al. 2008;Slosar et al. 2008;Matarrese & Verde 2008;Afshordi & Tolley 2008;Verde & Matarrese 2009;Desjacques & Seljak 2010;Schmidt & Kamionkowski 2010;Desjacques et al. 2011b;Schmidt et al. 2013).Additional operators involving the primordial Bardeen potential Φ (without any derivatives) must be added in the general bias expansion in order to model the scale-dependent corrections (in the same spirit as in the work of (McDonald 2008;Giannantonio & Porciani 2010;Baldauf et al. 2011)).In the case of an arbitrary quadratic PNG, the full set of operators was derived in Assassi et al. (2015), where up to second-order in perturbations and linear in fNL the bias expansion in the Eulerian frame becomes where q are the spatial coordinates in the Lagrangian frame and ε E Ψ is the stochastic counterpart of the field Ψ.The lat-ter is a non-local transformation of the primordial Bardeen potential given by, where α can take real values, which depend on the PNG type.This field originates from the generalization of the local ansatz in the case of a general quadratic PNG (e.g.Schmidt & Kamionkowski 2010;Scoccimarro et al. 2012;Schmid & Hui 2013;Assassi et al. 2015).
By generalizing the peak-background split argument (see Appendix C), the values of the additional non-Gaussian bias parameters can be calculated.The leading bias term in the squeezed limit for the case of a general nonlocal quadratic non-Gaussianity is (Desjacques et al. 2011b;Schmidt et al. 2013) ) where δc is the spherical collapse threshold in an Einsteinde Sitter Universe and the superscript X in fNL denotes one of the three non-Gaussian types considered here.We defined )dk, and different primordial shapes correspond to different choices of parameters: α = 2, A = 3 for the equilateral shape and α = 1, A = −3 for the orthogonal configuration, while the usual local case is described by A = 1 and α = 0 (see e.g.Dalal et al. 2008;Slosar et al. 2008;Giannantonio & Porciani 2010, and Appendix C for more details).The higher-order bias parameter result is given in the Eulerian framework by For the local case the above formula was derived in Giannantonio & Porciani (2010).Here we extend it to the equilateral and orthogonal cases [Eq.( 11)], by using the peak-background split argument (see Appendix C for details of the full derivation).In the case of equilateral PNG, the scale-dependent bias term approaches a constant value on large scales and thus becomes degenerate with the linear bias parameter, excluding the possibility of constraining f equil NL through this non-Gaussian correction.It must be noticed that the transfer function inside M (k) in Eq. ( 20) also introduces a scale dependence, due to PNG, in the bias at small scales.However, for large k (i.e.k ∼ 1/R * ), another strong degeneracy arises, this time between the PNG bias contribution and higher-order derivative bias terms (Assassi et al. 2015).Due to the presence of degeneracies on both large and small scales, the power spectrum does not essentially bear any constraining power on equilateral PNG.Note that, in our analysis, we never consider scales (kmax < 1/R * ≈ 0.25 h/Mpc at redshift z = 0) where higher-order derivative terms produce significant contributions.
The presence of PNG will affect the halo number density and will therefore introduce a small scale-independent correction in the Gaussian bias parameters.In Desjacques et al. (2009) and Sefusatti et al. (2012) expressions for these offsets 1 , up to the quadratic local-in-matter Eulerian bias parameter, were derived, where ν(M, z) = δc/σR(M, z) is the height of the peak (see Appendix C for a discussion), σR(M, z) is the variance of the smoothed density field on mass scale M and RNG is the non-Gaussian correction to the mass function, as derived in LoVerde et al. (2008).

Galaxy bias
In order to connect the halo bias formalism, presented in the previous two sections, with the bias of galaxies, one needs a method for describing the way galaxies populate the dark matter halos.One way to achieve this is to adopt a halo model (Cooray & Sheth 2002), where a Halo Occupation Distribution (HOD) function is used to provide the mean number density of galaxies ( N M ) per halo of a given mass M .The HOD model will be used in our Fisher matrix analysis only for predicting the fiducial values of the galaxy bias parameters, in the case of the radio continuum surveys.For optical surveys we will not use a HOD.We will instead use specific bias values, taken from the literature describing relevant galaxy populations (see Sec. 3.2 for details).The galaxy bias coefficients can be obtained from a weighted average of the halo bias over the range of host halo masses corresponding to the desired galaxy type as where n h (M, z) is the mean number of halos above mass M , given by the halo mass function, and N (M, z) is the mean number of galaxies per dark matter halo of mass M , given by the HOD model.Finally Mmin is the minimum mass a halo must have to contain a central galaxy and ng(z) is the mean galaxy density given by The HOD used here is a simple three-parameter model proposed by Tinker et al. (2005) and immediately adopted in the literature (e.g.Conroy et al. 2006;Sefusatti & Komatsu 2007).
) where M1 is the mass required for a halo to contain a second satellite galaxy.In Conroy et al. (2006) 0.76 log 10 (M1) + 2.3, by fitting the HOD free parameters on N -body simulations at different redshifts and number densities.In addition, it has been shown that the ratio log 10 (M1/Mmin) = 1.1 is almost redshift and mean number density independent.We will use both relationships in order to simplify the HOD model, leaving Mmin as the only free parameter.Finally, we set Mmin so that ng matches the expected mean galaxy number density of the survey in each redshift bin.

Real space
The halo bias expansion discussed in the previous section can be generalized for any kind of dark matter tracer.In the case of galaxies, the bias expansion in the Eulerian framework can be expressed up to second order in Fourier space, for non-Gaussian initial conditions, by Eqs. ( 7) and (8) as where the kernels S2(k1, k2) and N2(k1, k2) are defined as: The kernel S2 arises from the Fourier transform of the tidal field scalar s 2 (McDonald & Roy 2009;Baldauf et al. 2012), while the kernel N2(k1, k2) originates from the Fourier transform of the linear displacement field in the mapping between the Eulerian and Lagrangian frames.For non-Gaussian initial conditions the displacement field is coupled to the primordial gravitational potential, introducing additional terms in the tree-level galaxy bispectrum (Giannantonio & Porciani 2010;Baldauf et al. 2011).Henceforth, we drop the superscript E from the bias parameters, since we consider galaxy statistics at the time of observation.Before presenting the galaxy power spectrum and bispectrum model, we should note that we only consider up to linear terms in fNL, since we assume an fNL = 0 fiducial cosmology.In addition we would like to stress that for the bias expansion an implicit smoothing scale R is assumed for the density field and thus WR(k) factors are not retained in the expressions; this is allowed since we are working on scales much larger than the smoothing radius (kmax 1/R).Following Heavens et al. (1998), we can evolve the non-linear density field up to the desired order, smooth it with a filter and then apply the general bias expansion.The galaxy The shape of the galaxy bispectrum for Gaussian initial conditions, together with its highest contributing terms.In each panel the bispectrum, normalized to its absolute maximum value, is plotted as a function of k 2 /k 1 and k 3 /k 1 for all configurations in the case of fixed k 1 .We consider the following three values: k 1 = 0.01, 0.05, 0.1 h/Mpc, where the triangle sides satisfy k 3 k 2 k 1 .In the first row the Gaussian galaxy bispectrum B (G) without the trispectrum contribution is plotted [Eq.( 21)], in the second we plot the tree-level matter bispectrum (B G ) as predicted by SPT, in the third the quadratic bias term indicated as P 2 is plotted (i.e.P L m (k 1 )P L m (k 2 )+2 perm) and finally in the last row we plot the tidal bias term contribution (i.e.B S 2 = S 2 (k 1 , k 2 )P L m (k 1 )P L m (k 2 ) + 2 perm).Note that for the latter the absolute value is plotted, where a white dotted line shows the separation between the positive (left side) and negative (right side) values.For a detailed explanation, see the main text (Sec.2.3.1).
power spectrum can be written up to tree level contributions (Dalal et al. 2008;Slosar et al. 2008;Matarrese & Verde 2008;Verde & Matarrese 2009), where P L m is the linear matter power spectrum [Eq.( 1)], For the galaxy bispectrum, by keeping terms up to O(δ 4 ), and for Gaussian initial conditions we obtain where BG is the tree-level gravity-induced matter bispectrum (see e.g.Verde et al. 1998Verde et al. , 2000b;;Bernardeau et al. 2002, and references therein) and T δ is the trispectrum of the non-linear matter overdensity (i.e.δ(k)).Here Bε is the bispectrum of the leading stochastic field (i.e.ε(k)), Pεε δ is the cross power spectrum between ε and the next-to-leading order stochastic field (i.e.ε δ ).The full result including non-Gaussian terms up to linear order in fNL is then given by where δD( i ki)T δ G δ = δG(k1)δ(k2)δ(k3)δ(k4) , with δG being the Gaussian part of the density field (i.e. the Gaussian part of the primordial curvature perturbation, linearly propagated via Poisson equation), and Pεε Ψ is the cross power spectrum between ε and εΨ.The fiducial values of the stochastic terms in Eqs. ( 19), ( 21) and ( 22), are taken to be those predicted by Poisson statistics (as k 1/R * ), and are given by (Schmidt 2016;Desjacques et al. 2016): The non-Gaussian stochastic bias term in Eq. ( 22) can also be written as b1Pεε Ψ = bΨPεε δ .The Gaussian part of the galaxy bispectrum [Eq.( 21)], without the trispectrum contribution, is plotted in Fig. 1 for all the triangle configurations generated after keeping k1 fixed and k1 k2 k3.Three different values are chosen for k1, i.e. k1 = 0.01 h/Mpc, k1 = 0.05 h/Mpc and k1 = 0.1 h/Mpc.In addition, the highest contributing terms to B (G) ggg , are also shown in Fig. 1.These include the treelevel gravity-induced matter bispectrum (BG), the quadratic bias term (P L m (k1, z)P L m (k2, z)+2 perm), denoted in the plot as P 2 , and the tidal bias term (S2(k1, k2)P L m (k1)P L m (k2) + 2 perm), denoted as BS 2 .In order to show the shape dependence on the triangle configurations, the amplitude of each term is normalized to the maximum value in each panel.
In the second row of Fig. 1 the tree-level matter bispectrum signal is plotted, as derived in SPT.We can see that this term peaks mainly at the elongated (k1 = k2 + k3) and folded (k1 = 2k2 = 2k3) configurations, while for the squeezed triangles (k1 k2 k3) its contribution vanishes.This is due to presence of the second order SPT kernel F (s) 2 (ki, kj) [Eq.(A10)] in BG, which disappears in the squeezed limit and has a maximum at the folded/elongated triangles (see e.g.Sefusatti & Komatsu 2007;Sefusatti 2009;Jeong & Komatsu 2009, for a discussion).As we approach large scales (first column in Fig. 1) we observe that the maximum signal of the matter bispectrum is now at the equilateral triangles (k1 = k2 = k3).This is due to the fact that in this regime the matter power spectrum increases as a function of k and therefore we can get an excess in the signal of BG when all sides of the triangle are equally large, i.e. equilateral configurations (e.g.Jeong & Komatsu 2009).The quadratic bias term (shown in the third row of Fig. 1) follows a similar behaviour, where the only difference from BG is the absence of the F (s) 2 kernel.This leads to an enhancement at the squeezed limit and a suppression for the folded/elongated triangles.
The tree-level galaxy bispectrum contribution, proportional to the tidal bias (i.e.BS 2 ), peaks on elongated configurations, as we can see from the bottom row of Fig. 1.Note that, due to the presence of the S2 kernel, BS 2 can have an amplitude with a negative sign.In order to avoid the saturation of the colour maps in Fig. 1, we show the absolute value of BS 2 and we use a white dotted line to separate negative and positive BS 2 regions.For most configurations, BS 2 is negative on small scales (right column of Fig. 1), while the occurrence of positive values increases on large scales.Note here that, for all equilateral and for most isosceles triangles the tidal bispectrum term is negative, independently of the scale.This behaviour can be explained by the nature of the S2(ki, kj) kernel, which takes its maximum positive value (for simplicity ki = k1 and kj = k2) when k1 = ak2, where a > 1 (i.e.elongated and folded triangles), and its maximum negative value for k1 = k2 (i.e.equilateral triangles).In the folded limit we have, , and since the matter power spectrum increases towards large scales, the peak of the signal moves towards this configuration.On the other hand, for isosceles configurations (k1 > k2 = k3) the resulting tidal term can be positive or negative depending on the relative size of k1 with respect to the other sides of the triangle.
Let us now discuss in greater detail the trispectrum terms generated by the bias expansion in the galaxy bispectrum.The importance of the non-linear bias term in Eq. ( 21) was recognised in the work of Sefusatti (2009) and Jeong & Komatsu (2009) for increasing the sensitivity of galaxy bispectrum to the non-Gaussian initial conditions.Following the notation of Sefusatti (2009), each individual trispectrum term is expressed as (24) Here the upper indices (i.e.i, j, l, m) denote the expanding order of the non-linear density field.The primordial matter trispectrum (i.e.T1111 = T (4) ) for the local case depends on f 2 NL and gNL and can therefore be neglected in the Fisher analysis performed here.However, we should note that such term has a dominant large scale behaviour for the squeezed configurations and is larger than the non-Gaussian correction to the galaxy power spectrum [Eq. (19)].The part of the trispectrum generated by the non-linear gravitational coupling (i.e.T (6) δ ) exhibits no fNL dependence up to tree-level and therefore it can be ignored in the linear regime considered here.The important contribution for PNG constraints is a coupling term between a non-zero primordial bispectrum and the tree-level gravitational contribution, given by where kij = ki + kj and kij = |kij|.The above equation has a linear dependence on f X NL .On large scales it generates a signal that dominates over non-linear terms, for essentially all triangle configurations in the case of local non-Gaussianity (Sefusatti 2009;Jeong & Komatsu 2009).Therefore, the constraints on f X NL can be significantly improved, as we will show in Sec. 6.When terms proportional to the primordial field Ψ are also considered in the bias expansion, the corresponding trispectrum corrections in Eq.  1 for the Gaussian part).The panels display the amplitude of the galaxy bispectrum, normalized to the respective maximum value (note that this implies that a direct comparison of the color scale between different panels is meaningless).The non-linear evolution of the matter field is treated here with the MPTbreeze perturbation theory scheme (see Appendix B).In the middle panel the trispectrum loop quadratic bias correction (the b 2 trispectrum term in Eq. ( 21), i.e.
perm), is plotted.Finally, in the bottom panel the tidal bias term trispectrum correction (i.e.
All the terms plotted here peak at the squeezed limit for this type of PNG.
(22) (i.e.T δ G δ ) exhibits one occurrence of δG and therefore will be missing a permutation in Eq. ( 25).However, since we only consider T (5) in the tree level matter trispectrum, all the terms in Eq. ( 22) with T δ G δ will be O(f2 NL ) and hence they can be ignored.The only remaining O(fNL) trispectrum contribution is the one coming from the T δ term of Eq. ( 21).Its amplitude is shown in the colour maps of Fig. 2, 3 and 4, for the three PNG types considered here.Moreover, the non-Gaussian part of the galaxy bispectrum [Eq.( 22)] is also shown.
The colouring in the plots shows the shape of the non-Gaussian terms of the galaxy bispectrum, as a function of k3/k1 and k2/k1, for three different fixed values of k1 = 0.01, 0.05, 0.1 h/Mpc.In the local case, both the primordial bispectrum and the trispectrum correction peak in the squeezed limit.A small difference between the two is originating from the presence of the tidal kernel S2.The PNG contribution can be easily disentangled from the nonprimordial part of the galaxy bispectrum, as already pointed out earlier.On larger scales a small increase in the signal is also observed for all configurations, due to the behaviour of the matter power spectrum and of the F (s) 2 kernel in this regime (see Sec. 2.3.1 for a discussion).
In the equilateral case, the non-Gaussian scale dependence is easily observed.The peak of the signal moves towards large scales, for equilateral configurations.The PNG contribution of the galaxy bispectrum, considered here, includes only trispectrum corrections and BI (see Sec. 2.2.2), which explains why the panels of Fig. 3 display a similar behaviour.Note that the observed scale dependence is stronger and more localized, in the equilateral configurations, with respect to the non-primordial bispectrum signal (see first row in Fig. 1), in the large scale regime.This scale dependence can in principle provide a unique signature for measuring f equil NL .The same general scale dependent behaviour is observed also for orthogonal models and it can also improve the fNL constraints also this case.
The trispectrum integrals present an ultraviolet divergence, which is automatically cured by adding the smoothing filter with a finite value of R. Nevertheless, this introduces a dependence on the smoothing scale in the integration of the trispectrum for the three shapes we consider here 2 .This makes the results rely upon a non-fundamental quantity, which is unsatisfactory.On large scales this dependence on the smoothing radius goes like 1/σ 2 R , as was also noted in Jeong & Komatsu (2009).In order to cancel it, σ 2 R is included explicitly in front of the trispectrum integral and later on is reabsorbed by the bias parameters.The "new"  2 but for the equilateral type of PNG.For this type of PNG the galaxy bispectrum is taken to be that of Eq. ( 21) (i.e.B N G is here the sum of the trispectrum bias corrections), since any term proportional to the field Ψ in Eq. ( 22) (introduced to model the scale dependent bias corrections) is excluded as we discuss in Sec.2.2.2.MNRAS 000, 000-000 (0000) bias coefficients so obtained are then considered free parameters in the Fisher matrix analysis.
An alternative approach would be to use a perturbation theory that applies a renormalized technique, like the renormalized perturbation theory (RPT) (Crocce & Scoccimarro 2006b,a), time renormalized group model (Pietroni 2008) and renormalization of bias (McDonald 2006;Schmidt et al. 2013;Assassi et al. 2014;Senatore 2015;Mirbabayi et al. 2015).Regardless of the approach taken, the final result for the statistics of galaxies must be the same, therefore here we will use the MPTbreeze formalism (Bernardeau et al. 2008;Crocce et al. 2012) which simplifies greatly the computational effort of RPT (see Appendix B).The reasoning behind this choice is the exponential cut-off that is generated within this formalism and removes the UV divergence of the trispectrum integral by suppressing the small scales contribution.For the physical motivation and the details of the MPTbreeze formalism we refer the reader to (Bernardeau et al. 2008;Crocce et al. 2012).We note that, a reduction in the signal originating from intermediate scales is expected in this approach, due to the drop of the matter power spectrum and bispectrum beyond these scales (k > 0.15 h/Mpc at z = 0).This is shown in Lazanu et al. (2016), Lazanu et al. (2017) and Lazanu (2017), where a comparison between different perturbation theories is performed using simulations.The resulting power spectrum and bispectrum in the case of MPTbreeze are the same as those defined in Eqs. ( 19) and ( 22) respectively, but multiplied by a suppression factor that accounts for non-linear coupling [Eq.( B15)] (see Appendix B for a quick review).

Redshift Space
In order to model the statistics of galaxies in redshift space, the effect of redshift space distortions (RSD) (Kaiser 1987;Hamilton 1998) and the Finger-of-God (FOG) must be taken into account (Jackson 1972).The first can be treated perturbatively (Verde et al. 1998;Scoccimarro et al. 1999), where the kernel formalism of SPT (see Appendix A) can be generalised to include the redshift distortions and the bias terms [Eq.( 17)] and the second is described phenomenologically by a damping factor D P FOG and D B FOG for the power spectrum and bispectrum respectively.Since our analysis is restricted to large scales, we only require up to the second order redshift kernels in order to derive the linear power spectrum and tree-level bispectrum in redshift space.For the general, non-local, PNG the power spectrum and bispectrum are given by where the general non-Gaussian redshift kernels up to second order are given explicitly in Appendix D. They are reported in a slightly different notation compared to Tellarini et al. ( 2016), and we also neglect O(f 2 NL ) contributions.The term T (5) in redshift space, after neglecting all O(f 2 NL ) contributions, derived by using standard PT formalism, the bias expansion of Eq. ( 17) and RSD up to second order, is given by where f is the linear growth rate.Here 2 (k1, k2) are the Gaussian parts of the redshift kernels Z1 and Z2 respectively [Eqs.(D3) and (D4)], while for Z G,b 2 we also exclude the two SPT kernel contributions.The terms denoted FP i and GP i are the ith permutation of Eq. ( 25), where the letter F and G represent the SPT kernel used in the expression at hand (e.g.GP 4 = δ (1) A redshift space model similar to the one in Eq. ( 27) was used in Tellarini et al. (2016) to study expected constraints on the PNG amplitude for a list of future LSS surveys.However, in this reference, only the local case without trispectrum contributions and stochastic bias terms is considered.Besides these new terms, redshift uncertainties are also included in our redshift model, as discussed in the next paragraph.
The FOG term models the damping effect of the clustering power induced by the FOG effect on linear scales.Peacock & Dodds (1994) and Ballinger et al. (1996) find that a good model for the power spectrum as an exponential function and an analogous function can be written for the bispectrum Here we consider the fiducial values of σP = σB = συ(z), where συ is the usual linear, one dimensional velocity dispersion.Besides the FOG effect, the redshift uncertainty of galaxy surveys must be also taken into account.The redshift error, σz, can be translated into a position uncertainty along the line of sight.The treatment of this effect is the same as in the case of FOG (see e.g.Seo & Eisenstein 2003), where the only difference is the fiducial value of the σ parameter, which will be σr = cσz(z)/H(z).Considering both effects gives the final form of the damping factors in Eqs. ( 26) and ( 27), with the σ parameters given by σ 2 v = σ 2 υ + σ 2 r .These multiplicative factors introduce a suppression of the signal for all scales with a large component along the line of sight, affecting mostly small scales.In other words, only k-modes satisfying kµσv 1 are not dominated by noise and can contribute to the power spectrum and bispectrum measurements.
The redshift space bispectrum is characterized by five variables, three of them defining the triangle shape (e.g. the magnitude of the three wavenumbers, k1, k2, k3) and the remaining two the orientation of the triangle with respect to the direction of the line of sight ẑ, which we consider to be the polar angle µ1 and the azimuthal angle φ.The bispectrum will now be B s g (k1, k2, k3) = B s g (k1, k2, k3, µ1, φ).Taking the average over angles (i.e. the monopole term in the Legendre expansion) of Eq. ( 27), in a similar fashion as it was done in Kaiser (1987) for the power spectrum, one can obtain (Sefusatti et al. 2006;Gil-Marín et al. 2012) where The terms P r g and B r ggg are the real space galaxy power spectrum and bispectrum, given by Eqs. ( 19) and ( 22) respectively.The redshift space bispectrum presented above, as described in Sefusatti et al. (2006), is derived after an additional average over θ12 = acos( k1 k2) and by dropping the dependence on the second-order PT velocity kernel [Eq.(A14)] and the FOG effect.This is a good approximation on large scales since these two effects partially cancel out.
In the cases of local and orthogonal PNG the galaxy power spectrum and bispectrum is described by the full model of Eqs. ( 19), ( 22) and ( 26), ( 27) for real and redshift space respectively.In Eq. ( 27) we will keep only the O(fNL) terms, while the full form was written down for completeness.For the equilateral case, as we discussed before, the scale dependent bias contribution is degenerate and will not be used.Nevertheless we will use the Gaussian power spectrum, excluding its signal contribution from constraining fNL, and the bispectrum without the terms proportional to the non-local primordial field Ψ.In this case, the trispectrum bias contribution can compensate for the missing large scale signal and improve the expected constraints on the non-Gaussian amplitude, as we will show in Sec. 6.

GALAXY SURVEYS
In this Section we describe the specifications of the radio continuum and optical galaxy surveys we used in this work.

Future radio surveys
We focus on future radio surveys to investigate whether the high-redshift, full sky nature of those datasets will provide better constraints on non-Gaussianity parameters, despite the lower precision in redshift information.
We forecast results for two types of radio continuum surveys, assuming either a 10µJy or a a 1µJy flux limit.In both cases we consider 30, 000 deg 2 .
Given that we want to forecast the advantages of having both full sky and high-z data, we use radio continuum surveys, that are however plagued by the fact that the sources' redshifts are in principle not known.This would allow only  .The galaxy bias, defined from the weighted average over the halo one, as a function of redshift for the two radio surveys considered here, a survey flux limited at 10µJy and one with a 1µJy flux limit.The details of the surveys are listed in Table 1.
the computation of an angular projected power spectrum and bispectrum, in the standard case.
However, recently some techniques have been developed in order to provide such surveys with statistical redshift information; here we follow the Clustering-Based Redshift (CBR) technique, developed in Schneider et al. (2006), Newman (2008), Ménard et al. (2013) and studied for some cosmological applications, including the SKA in Kovetz et al. (2017).A disadvantage of this method is that the provided redshift information can have a large uncertainty.In fact it is safe to assume that the redshift error σz(z) is equal to the width of the bin considered.The resulting predicted redshift distribution and redshift errors for the radio surveys used here are presented in Table 1.
The galaxy bias for the two radio galaxy surveys is given by the modelling described in Sec.2.2.3, while the results for the linear and quadratic bias are plotted in Fig. 5

Future optical surveys
In the forthcoming years a plethora of large scale structure surveys will provide maps of various galaxy types over large volumes.The large number of modes and high redshifts probed by these surveys will give the possibility to produce tight constraints on the amplitude of PNG.In this section we present the forecasts for spectroscopic and a photometric experiments, covering a variety of redshift ranges and volumes.
For the spectroscopic survey we consider the redshift range of 0.7 z 2 over 15, 000 deg 2 , while we model the redshift distribution as in Table 2 from the one in Orsi et al. (2010).
In Table 2 we show the main specifications of the survey used in this work, which is not too dissimilar from an Euclidclass or DESI-class survey (Laureijs et al. 2011;Font-Ribera et al. 2014).The adopted size of redshift bins is ∆z = 0.1, while the spectroscopic redshift error is given by σz(z) ≈ 0.001(1 + z).
The galaxy sample is assumed to involve a single tracer, whose linear bias is given by b1(z) = 0.76/D(z) (Font-Ribera et al. 2014).The higher order Eulerian bias coefficients are given from the halo bias predictions (see Appendix C), where ν is determined by b1 The results for the linear and quadratic cases can be seen in Fig. 7.
In the photometric survey case, we take inspiration from the future LSST survey (LSST Science Collaboration et al. 2009), and we use a redshift distribution as in Fig. 6, where we divide our sample in eight redshift bins covering the range 0 z 3 over 18, 000 deg 2 .The distribution of the true redshifts of galaxies inside each photometric bin i is given by (Ma et al. 2006;Zhan 2006) with z ph being the photometric redshift, and where the integration limits define the extent of bin i.Following the work  33), is also shown (dashed, dotted, etc. lines).
of Ma et al. (2006), we model the photometric redshift error at each redshift by a Gaussian, given by: where the photometric redshift bias z bias and the photometric rms error σz are functions of the true redshift.The latter is chosen here to be σz(z) = 0.05(1+z).The fiducial value of the redshift bias is chosen to be z bias (z) = 0. Note that the above probability must be normalised to ∞ 0 dz ph P (z ph |z) in order to ensure the positivity of the photometric redshifts.
The overall galaxy true redshift distribution n(z) = d 2 N/dzdΩ is given by the functional form, n(z) ∝ z α exp[−(z/z0) β ] ( Wittman et al. 2000), with α = 2, β = 1 and z0 = 0.3125.The normalization of the overall galaxy redshift distribution is fixed by requiring that the total number of galaxies per steradian (i.e.ntot = ∞ 0 dzn(z)) to be equal to the cumulative galaxy counts of the survey.The overall normalized redshift distribution n(z) is plotted together with the true redshift distribution in each photometric bin in Fig. 6.The final number densities for each photometric redshift bin are listed in Table 3.The fiducial value for the linear bias is chosen to be b1(z) = 1/D(z) (LSST Science Collaboration et al. 2009;Font-Ribera et al. 2014;Passaglia et al. 2017), where the derivation of the higherorder bias parameters follows the prescription described before.The final results for the linear and quadratic biases are plotted in Fig. 7.We use a surface number density of ntot = 40 gal/arcmin 2 .These specifications are chosen to mimic an LSST-like, representative future photometric survey (see e.g.LSST Science Collaboration et al. 2009;Zhan & Tyson 2017).

Fisher Matrix
In order to forecast the precision in measurements of the non-Gaussianity amplitude fNL and bias parameters, we make use of the Fisher matrix formalism.The Fisher information defines the minimum error on an unknown parameter given an observable.In the case of the galaxy power  spectrum, the Fisher matrix is given by while for the bispectrum we have (Scoccimarro et al. 2004) where p α,β are the unknown parameters of interest and the derivatives are evaluated at the fiducial value of the parameter vector.We model power spectrum and bispectrum as discussed in Sec.2.3.The sum over the triangles, which is denoted in Eq. ( 36) as T , is given by applying the triangle condition on the k modes, while the ordering considered is k1 k2 k3.The wavenumbers are binned with a bin size of ∆k, taken here to be the fundamental frequency of the survey k f = 2π/V 1/3 , between a minimum value kmin = k f (largest scales probed by survey) and kmax (smallest scales considered).The angular bin sizes are taken here to be ∆µ1 = 0.1 and ∆φ = π/25.
The parameter vector considered here per bin consists of the non-Gaussian amplitude, three bias parameters, stochas-tic power spectrum and bispectrum contributions, the linear growth rate and the velocity dispersion, The stochastic bias contributions (i.e.Pε, Pεε δ and Bε) are considered here as nuisance parameters and they are marginalised over to acquire the Fisher sub-matrix for the parameters of interest.Cosmological parameters are fixed throughout this work.Since we can determine them with high accuracy, using other probes (e.g.CMB, BAO etc.) we do not expect a significant impact of propagating their errors on fNL, which is instead constrained via power spectrum scale dependence on very large scales, or via the bispectrum.This is actually verified for the CMB primordial bispectrum in Liguori & Riotto (2008), and for that case it is also shown that the level of degeneracy between cosmological parameters and fNL is small 3 .The final forecasts for the non-Gaussian amplitude are derived after marginalising over the remaining parameters and summing the contribution from all the available redshift bins.We aim for a complete and conservative analysis, therefore we will restrict the analysis to the linear regime, and exclude non-linear scales.For the redshift evolution of kmax(z) we consider, kmax(z) = 0.1/D(z) h/Mpc where it slowly varies with redshift while it stays inside the linear and semi-linear regime, ensuring the validity of the bias expansion as well as SPT itself.This limit in the magnitude of the k mode holds also for the bispectrum, where the sides of the triangle satisfy kmin k1, k2, k3 kmax.
In our Fisher matrix analysis, only the diagonal part of the covariance matrix (i.e.∆P 2 and ∆B 2 ) is taken into consideration, neglecting all the cross-correlations between different triangles (bispectrum) and k-bins (power spectrum).We adopt a Gaussian approximation for the variance terms: where s123 = 6, 2, 1 for equilateral, isosceles and scalene triangles respectively.The volume of the fundamental shell in Fourier space is V f = k 3 f .In addition Ptot(k, z) = P s g (k, z) + 1/ng, where the stochastic contribution is excluded from P s g and the remaining term accounts for the shot noise.Note that for the Fisher matrix in redshift space, the normalization for the range of the two angles, i.e.Nµ = µmax − µmin and N φ = φmax − φmin, must be applied.The full covariance of the two estimators is outlined in Sefusatti et al. (2006), where the off-diagonal elements are defined by higher than three point correlators.Although the above results are for redshift space statistics, the reduction to real space is straightforward.
Recently, the authors of Chan & Blot (2017) used dark matter N -body simulations, including four halo samples with different number densities, in order to study the full covariance of the power spectrum and bispectrum estimators.They focused on extracting an integrated signal-tonoise ratio for all bins/triangles, checking how this gets degraded when off-diagonal covariance elements and non-Gaussian contributions to the variance are accounted for.While this does not include a specific study of the degradation of error bars for PNG or other cosmological parameters, their results can provide useful guidelines to assess the validity of our diagonal covariance approximation and the error on the parameter forecast we introduce by employing it.
For the dark matter power spectrum, Chan & Blot (2017) show that the correlation coefficient between different modes never exceeds ∼ 15% at z = 0, up to kj = 0.1 h/Mpc.For z = 1, the correlation reaches at most ∼ 20%, up to kj ∼ 0.15 h/Mpc.In all cases, correlations decrease away from the diagonal.On the scale range considered here, the non-Gaussian corrections to the diagonal part of the covariance is negligible at z = 0, as well as for z = 1.For halos, these corrections can be up to ∼ 10% for the same scale range, in the case of the small density halo samples, with the exact value depending on the redshift.Furthermore, the results agree with the covariance model predicted by PT up to k = 1 h/Mpc for the abundant halo sample.Therefore we conclude that the exclusion of the off-diagonal part in the galaxy power spectrum covariance will introduce an error of up to 10 − 15%, for the scale range and redshifts considered here, depending on the number density of the sample.In addition, the non-Gaussian corrections to the variance are negligible for the high density samples and scale range considered here.
For the dark matter bispectrum, it is shown that the non-Gaussian corrections to the diagonal Gaussian part (for equilateral configurations) is ∼ 8% at z = 0 and k = 0.1 h/Mpc and that PT predictions agree with the numerical results up to k ∼ 0.15 h/Mpc.For higher redshifts (i.e.z = 0.5, 1) the corrections are at the level of a few percent, up to k ∼ 0.16 h/Mpc, while the PT predictions are in good agreement with the results up to k ∼ 0.2 h/Mpc and k ∼ 0.3 h/Mpc for z = 0.5 and z = 1 respectively.In addition, it is shown that the correlation coefficient, used to test couplings between different triangles, is consistent with zero for the large scales and for the redshift slices considered.
When, instead of dark matter, we consider halos, triangle couplings and non-Gaussian corrections strongly depend on the redshift and density of the sample.Corrections to the diagonal Gaussian part are negligible at z = 0 and k = 0.1 h/Mpc in the case of a high density sample (n ∼ 10 −3 [Mpc/h] −3 ).On the other hand, for a low density case (n 2 × 10 −5 [Mpc/h] −3 ) the correction is found to be between few percent and ∼ 10% at low redshift, and increasing up to ∼ 90% at high redshifts.
The samples considered in this work have a high number density for the majority of the redshift bins (except only for some high redshift slices, where non-linearities are anyway less important), ng 10 4 (h/Mpc) 3 (see Table 1 and Appendix 3).Therefore we do not expect the exclusion of the non-Gaussian part in the covariance to overestimate much the S/N ratio and to have a large impact on the final PNG forecasts.
As explicitly shown in Chan & Blot (2017), on the large scales considered here the full non-Gaussian contribution can be well approximated by including perturbative corrections to the power spectrum appearing in the bispectrum variance expression, obtaining where P NL g (k3) is given by Eq. ( 19) after replacing the linear matter power spectrum with the non-linear one, as predicted by the HALOFIT algorithm (Smith et al. 2003;Takahashi et al. 2012).The reason that HALOFIT, instead of the PT one-loop matter power spectrum, is used, lies in the fact that the latter leads to overestimating the actual variance in the weakly non-linear regime.We will use this expression later on in our analysis, in order to estimate in more detail the effect of neglecting non-Gaussian corrections in our forecasts.
Another aspect to consider is the covariance between the power spectrum and the bispectrum (PB) when the two are used jointly.This was provided, in the Gaussian case, in Sefusatti et al. (2006) and used in Song et al. (2015).In Chan & Blot (2017), a comparison between the S/N ratios coming from N -body simulations is made in order to test the effect of the P B cross-covariance.In the case of dark matter, they show that, in the redshift range 0 z 1, on linear scales, for abundant samples as those considered here, P B corrections to the overall signal-to-noise amount at most to ∼ 10%.This justifies neglecting the cross-covariance in our forecasts.The combined Fisher matrix will therefore be just obtained as

Theoretical Errors
The final ingredient of our analyses is the consideration of theoretical errors in the Fisher matrix formalism.They quantify the uncertainties on the modelling of the matter perturbations and the bias expansion in the statistics of galaxies.In the majority of parameter forecasting analyses they are neglected, while a perfect knowledge of the theoretical model is considered.
Here we will follow and extend the treatment of Baldauf et al. (2016), where theoretical errors e are defined as the difference between the true theory and the fiducial theoretical prediction.This formalism considers the true theory to be the model which takes into account at least one more perturbative order than the fiducial one.These errors are bounded by an envelope E and their variation as a function of wavenumbers is characterized by ∆k [Eq.( 42)].The value of the correlation length ∆k is taken to be that of the smallest coherence length of the total power spectrum, that is the scale of the Baryon Acoustic Oscillations (BAO), ∆k = ∆BAO = 0.05 h/Mpc (see Baldauf et al. 2016, for an extensive discussion).The error covariance matrix is written as: where i, j are the indices of the different momentum configurations (i.e.number of bins and triangles for the power spectrum and bispectrum respectively).The correlation coefficient ρij accounts for the correlations between the momentum configurations, considered to follow a Gaussian distribution, and given by For a diagonal error covariance (i.e.ρii = 1 and ρij = 0 for i = j) and a fixed ∆k, the envelope E(k) would be independent of the bin size, contrary to the statistical errors.This means that, for uncorrelated bins, choosing a smaller bin size will reduce the effect of the theoretical errors.The presence of an off-diagonal ρij ensures that this does not happen and the relative impact of errors is independent from the size of the k bins.After marginalising over the theoretical errors e, the final covariance that will be used in the Fisher matrix analysis becomes just the sum of the variance of the power spectrum and bispectrum estimators (Eqs.( 38) and (39) respectively) with the theoretical error covariance [Eq.( 41)].
One of the goals of this work is to test the effect of theoretical errors on the expected parameter constraints coming from high redshift LSS surveys.The Universe is more linear at high redshifts and hence, for the scales considered in this analysis, we would not expect to see a significant impact on the constraints solely from the theoretical uncertainties attributed to PT.In the formalism proposed by Baldauf et al. (2016), the envelope is fitted up to two-loops in matter perturbations for both power spectrum and bispectrum while for the bias expansion they consider only the linear bias.Here, we proceed in extending their approach to include the theoretical uncertainties coming from the localin-matter bias terms (i.e.b1, b2, b3, etc.) that appear up to the 1-loop expression of the galaxy power spectrum and bispectrum.This set of terms has been shown to provide an accurate description by comparing with simulations and galaxy catalogues (see e.g.Scoccimarro et al. 2001b;Feldman et al. 2001;Verde et al. 2002;Marín et al. 2013;Gil-Marín et al. 2014).Note that the inclusion of all the bias terms at each order (e.g.including tidal terms and other operators, see Desjacques et al. (2016) for a review) is potentially important and we will consider their contribution in the theoretical error formalism in the near future (see Sec. 6.1.1 for a discussion).
In order to quantify the theoretical uncertainties, we fit the galaxy power spectrum and bispectrum envelope after including the local-in-matter bias terms up to 1-loop, in addition to the 1-loop matter expressions originating from the description of PT, while assuming Gaussian initial con-ditions.More precisely, for the galaxy power spectrum we take into account all the terms originating from the bias expansion, that have a dependence on b1, b2 and b3, up to 1-loop (i.e.Pg,11, P I g,22 and P I g,13 of Sefusatti ( 2009)).For the galaxy bispectrum, the bias terms considered are all those with a dependence on b1, b2, b3 and b4, up to 1loop, while we exclude those with a dependence on PNG initial conditions (see e.g.Sefusatti 2009, for details).For the matter expansion, we consider up to 1-loop terms for both power spectrum and bispectrum in the cases of SPT (see e.g.Bernardeau et al. 2002, for a review) and MPTbreeze (see Appendix B for details).The envelopes for the power spectrum and bispectrum in the the SPT and MPTbreeze schemes are given by MPTbreeze.
(43) and where k = (k1 + k2 + k3)/3.Note that the fitted coefficients in the above envelopes exhibit a small dependence on the fiducial cosmology, which is negligible for our purposes.The numerical values presented correspond to the cosmological parameters considered in this work.The SPT envelopes are presented here for completeness, since, throughout this work the MPTbreeze description of the matter perturbations will be used.
A slightly different methodology to quantify the effect of theoretical errors was developed in Audren et al. (2013).The procedure followed is similar to the one used here.The main difference is that they define theoretical uncertainties by fitting a correction function to the HALOFIT, instead of the desired loop order given by PT.Their technique is applied only to the power spectrum and the relative error to the non-linear P is added to the diagonal part of the covariance matrix.The level of correction is quantified by the precision of HALOFIT, where only the linear bias is taken into account.In order to quantify the difference between the two approaches, we tested this methodology for the power spectrum, since the fitting function is provided only for this case, while a generalisation to the halo bispectrum is not straightforward.To perform the comparison, we adjusted the envelope function of Eq. ( 43) by removing the higher order bias contributions, as well as by using only the diagonal part of the error covariance [Eq.( 41)].In the case of local PNG the difference on the final fNL constraint turns out to be ∼ 2%, while for orthogonal PNG is ∼ 8%.Therefore the two approaches produce consistent results, whenever a comparison is possible.We prefer to use the approach of Baldauf et al. (2016) because in this way we take also into account higher-order bias corrections, which are important for the redshift range we consider.
At low redshifts using only the diagonal part of the error covariance underestimates the effect of theoretical uncertainties, but not at higher redshift, where non-linear corrections are less important.Given that most of the volume and thus the statistical power is at high redshift we expect the difference between the two methodologies to be at the level of a few percent, in the case where we consider only the linear bias terms for the galaxy power spectrum.

SUMMARY OF THE ANALYSIS METHOD
Here we summarize the modelling used for the galaxy power spectrum and bispectrum, as well as all the steps followed in the next section in order to test the various effects that have an effect on the fNL forecasts.
(i) For the matter power spectrum and bispectrum, the MPTbreeze tree-level description is used throughout.In order to cross-validate our results, we have also considered the scenario of using SPT.The forecast results are in agreement between SPT and MPTbreeze (see Sec. 6.1.1).(ii) Non-Gaussian initial conditions are assumed and a bias expansion up to second order is used [Eq.( 17)].(iii) The forecasts are performed by following the Fisher matrix formalism (see Sec. 4.1), where a diagonal covariance is used for both power spectrum and bispectrum as described in Eqs. ( 38) and ( 39) respectively.The analysis is restricted to scales k kmax = 0.1/D(z) h/Mpc.(iv) We start by considering an idealised scenario, in which we adopt the approximate monopole statistics defined by Eqs. ( 31) and ( 32), for the galaxy power spectrum and bispectrum respectively.Here, the use of the term "approximate" stems from the fact that, in Eq. ( 32), we neglect a number of terms appearing in the full bispectrum monopole expression, as shown e.g. in Tellarini et al. (2016).As explicitly checked, the terms we neglect have a small impact on the large scales included in our analysis, making our bispectrum formulae accurate enough for this preliminary step.See Sefusatti et al. (2006) for a further discussion on this point.The idealisation level is of course much greater for the photometric/radio survey case as in reality the redshift errors are very large.Note that, for this warm up exercise, we are also assuming perfect determination of redshifts, neglecting theoretical errors and also excluding the trispectrum term shown in Eq. ( 22).Despite its lack of realism, it is quite useful to start with this simplified case, for several reasons.Firstly, a nontrivial amount of previous literature includes only real space forecasts, or angle averaged statistics in redshift space, without theoretical errors.Therefore, explicitly including this initial step facilitates comparisons.Secondly, by proceeding in this way we are able to better isolate the relative impact of different effects, such as theoretical and redshift errors, when we add them in subsequent steps.Thirdly, as explained in more detail in the following, off-diagonal terms in the theoretical error covariance are often neglected in our full redshift space analysis, due to computational limitations.On the other hand, off-diagonal terms can be in most cases fully included at this initial stage.This provides important guidelines to assess the accuracy of the diagonal covariance approximation, used later when needed (we anticipate here that we find this diagonal approximation to be quite good).Finally, by comparison of our preliminary monopole forecasts with the final, full redshift space analysis, we are able to assess the amount of PNG information lost by averaging over angles (i.e. the monopole).We will refer to this step of the analysis as to our "idealised" case, in the following.(v) The second step is to add the theoretical error covariance given in Eq. ( 41) to the variance of Fisher matrix, in order to account for the uncertainties in the theoretical model (see Sec. 4.2 for a discussion).Here we account only for the exclusion of 1-loop corrections in the matter perturbations for both power spectrum and bispectrum.For the bias expansion we only quantify the error in excluding the 1-loop contributions that are related to the local-in-matter bias coefficients, i.e. b1, b2, b3, etc. (vi) The effect of the theoretical errors on forecasts is shown as a function of redshift in Fig. 8 and 9 for the two radio continuum surveys considered here.The effect on the forecasts coming from the summed signal over all redshift bins is shown in Tables 4 and 9 for the radio and optical surveys respectively.(vii) The third step is to move to redshift space and include the full RSD treatment up to second order.The galaxy power spectrum and bispectrum model in redshift space are given by Eqs. ( 26) and ( 27) respectively.Note that the trispectrum term in Eq. ( 27) is still excluded for now.Only the diagonal part of the theoretical error covariance is used in the redshift space models.As also mentioned just above, we argue in Sec.6.1.1 that the effect of the off-diagonal part on the final fNL forecasts is small.Note, finally, that we are still not including redshift errors.Therefore our forecast up to this point are still unrealistic for optical photometric and for radio surveys, for which these errors are important.The goal, up to here, remains that of assessing the impact of each ingredient added at each new step of the analysis.Moreover, forecasts with no redshift errors included provide upper limits to the performance of realistic radio continuum surveys, eventually approachable by improving current techniques for redshift determination.(viii) In addition to the RSD effect, we consider redshift uncertainties, which are modelled like Eqs. ( 29) and (30) for the power spectrum and bispectrum respectively (see Sec. 2.3.2 for a discussion).The effect of RSD, theoretical errors and redshift uncertainties to the fNL forecasts is shown in Tables 5, 6 and 11 for the radio and optical surveys respectively.(ix) Finally, we take into account the trispectrum term in the galaxy bispectrum for both the idealised and the full RSD models (see Eqs. ( 22) and ( 27)).The effect of this trispectrum correction to the final PNG forecasts is shown, for the idealised case, in Tables 7 and 10 for the radio and optical surveys respectively.For the fi- Radio continuum, 1 Jy Figure 8.The effect of theoretical errors on the bias parameters and on the three f N L parameters per redshift bin (i.e. two times the value of σz in Table 1) for the 1 µJy radio continuum survey.The dashed lines represent our "idealised model", using a simplified monopole statistic for the bispectrum (see item iv in Sec.5), without theoretical or redshift errors (i.e. the starting step of our analysis, as explained in the main text); for the power spectrum we used the linear predictions with PNG [Eq.( 31)], while for the bispectrum of galaxies we used Eq. ( 32).The non-linear evolution is treated within the formalism of MPTbreeze.The solid lines represent the same model, but including theoretical errors, as described in Sec.4.2.Forecasts for the power spectrum are plotted in blue, for the bispectrum in green and for both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.The dotted grey line indicates the best constraints on the PNG amplitude, as given by Planck Collaboration et al. (2016b).
nal "full" model (i.e.RSD+theoretical errors+redshift errors+trispectrum), the fNL forecast are shown in Tables 8 and 11 for the radio continuum and optical surveys respectively.(x) We summarize our final forecast results on the amplitude of PNG coming from future LSS surveys in Table 12.

RESULTS
In this section we present the results of our forecasts for radio continuum and optical surveys, obtained with the procedure summarized in Sec. 5.For the local and orthogonal PNG, we consider the power spectrum and bispectrum model that was described in detail in Secs.2.3.1 and 2.3.2, for real and redshift spaces respectively [by Eqs. ( 19), ( 22), ( 26) and ( 27)].In the case of the equilateral type of non-Gaussianity, we use the same expressions, but without the corrections from the primordial local gravitational potential Ψ due to degeneracies with the bias parameters (see Sec. 2.2.2).The final results are derived after marginalising over the nuisance stochastic bias parameters (see Sec. 4.1).The linear matter power spectrum is computed with CAMB (Lewis et al. 2000), while the cosmological parameters are those determined by Planck Collaboration et al. (2016a): h = 0.6774, Ωch 2 = 0.1188, Ω b h 2 = 0.0223, ns = 0.9667, ∆ 2 ζ = 2.142 × 10 −9 , τ = 0.066.

Radio surveys
In this section we present the results of our forecasts for the parameter vector [Eq.( 37)], considering the two cases described earlier, i.e. the 10µJy and the 1µJy flux limit high redshift surveys.The modelling used is described in the previous paragraphs.Complications in the modelling are added in several steps, while the intermediate results are also presented.The reader interested in the results for the "full" case can skip to Sec. 6.1.3.

Theoretical errors effect
Our main goal at this stage will be quantifying the effect of theoretical errors, by testing their impact on our initial, idealised forecasts.Since we are not yet including redshift errors and other important effects in the analysis, we immediately warn the reader that our fNL forecasts are still unrealistic in this section.The purpose here is only the relative comparison between the two cases, with and without theoretical errors, in order to assess the impact of the latter.Due to the demanding computational effort needed to perform the inversion of the bispectrum covariance matrix (a 10 6 × 10 6 size matrix for the final redshift bins), which is no longer diagonal when the theoretical errors are included [Eq. ( 41)], we use the full covariance matrix only for the three lowest redshift bins (that's the largest computationally affordable amount; we choose the lowest three redshifts because the effect of off-diagonal terms is largest there), while for the rest we consider only the diagonal contribution.This means that cross-correlations between modes are excluded.This can in principle affect the impact of the theoretical errors.However, we performed tests to check the effects of these off-diagonal components, up to the redshift bin allowed by the computational resources available, and we observed that the effect of the off-diagonal terms becomes actually negligible at high-z.This is reasonable: the Universe is more linear at higher redshifts, therefore the loop corrections, up to the scales we consider, are expected to be suppressed.Our approximations work therefore very well.Let us note that, in the case of the power spectrum, we always use the full covariance matrix, since no computational issues arise in this case.Let us also note here again that the theoretical modelling was performed using MPTbreeze, which has an embedded cut-off function at high-k.This means that higher order contributions, as well as the theoretical error effect will be suppressed on small scales.In order to check whether this has a significant effect, we have performed a similar analysis using SPT.The results were consistent with those presented here, throughout the range of scales chosen for our analysis.
The comparison is shown in Fig. 8 and 9 and the forecasts for fNL are quantitatively reported in Table 4.In addition, the ratio of the marginal 1σ error is presented for the two radio continuum surveys in Fig. 10 and 11, for the cases where theoretical errors are taken into account and when they are omitted.
Note that the effect of theoretical errors doesn't get smaller at high redshifts in all cases, as seen in Fig. 10 and 11, and in particular this is evident for the power spectrum case.The argument of a more linear Universe at high redshifts, which leads to a reduction in the contribution of the higher loop corrections, holds for the case of matter and it is not generally true for the biased tracers.The reason is that at large redshifts, and in the scale range considered here, galaxies become more biased, while the matter loop corrections become less important.The final effect of theoretical errors is determined by adding the monotonic decreasing contribution of matter higher-order terms and the increasing contribution coming from higher-order coefficients in the bias expansion (this also explains the spikes observed in Fig. 10 and 11).Therefore, depending on the redshift evolution of the bias parameters, the depth of the survey and its volume, a different behaviour of the impact of the theoretical errors can be observed, when bias loop corrections are taken into account in the envelope fitting (see Sec. 4.2 for details).This shows the importance of extending the formalism of Baldauf et al. (2016), to include theoretical errors attributed to the bias expansion, as we did.
It would seem at this stage that the combined information from the power spectrum and bispectrum of galaxies can provide very tight forecasted constraints on the amplitude of local PNG.For the 1 µJy radio continuum survey, an error of σ(f loc NL ) < 0.13 can be achieved after using the total signal from all redshift bins.In the orthogonal case, this survey can provide σ(f orth NL ) ∼ 3, while for equilateral PNG, already at this stage the expected constraints are weaker, Radio continuum, 1 Jy Figure 10.Same as Fig. 8, where here the ratio of the expected 1σ errors between the case with theoretical errors (denoted in the plots with the upper index TE) and without is shown.The spikes observed in the expected forecasts, obtained when including the effect of theoretical errors can be attributed to the trade-off between the contributions coming from higher-order terms in the matter and bias expansions (see main text in Sec.6.1.1 for a discussion).9, where here the ratio of the 1σ errors between the case with theoretical errors (denoted in the plots with the upper index TE) and without is shown.Table 4. Forecasts for the non-Gaussian amplitude f NL for the three shapes (local, equilateral and orthogonal) over all redshift bins, from the two surveys considered here (the 1 µJy -left and the 10 µJy -right, radio continuum surveys), in the idealised analysis, no theoretical errors case and the case including theoretical errors.The expected constraints were derived from the power spectrum (P), bispectrum (B) and by combining the two (P + B).
compared to the other shapes, as we can see from Table 4. Let us however stress again that these results will deteriorate considerably when considering realistic redshift space measurements and especially when accounting for errors in the determination of the redshift of radio sources (see next subsection).As specified since the beginning of this section, at this stage we do not focus yet on the absolute value of the expected constraint, but on the assessment of the effect of theoretical errors.
Forecasts on bias parameters are shown in the left panels of Fig. 8 and 9 for each of the surveys considered here.Again, rather than quantitative assessments, which will be refined in the redshift space section later on, it is useful at this level to focus on qualitative behaviours.We see, as expected, that the power spectrum provides the main signal for constraining b1 (blue lines in Fig. 8 and 9), while the contribution from the bispectrum (green lines) is minimal as we can also see in the results coming from the combination of the two (red lines).For the quadratic bias parameter, as well as the tidal bias, the expected constraints are weaker.For both these bias coefficients, the constraining signal originates solely from the bispectrum, since these terms appear in loop corrections of the linear galaxy power spectrum, which we do not consider here, and therefore this reduction in the constraining power is justified.In the combined power spectrum and bispectrum case, an improvement is observed in the statistical error of b2 and b s 2 , due to the tight constraint on the linear bias provided by the power spectrum.In addition, the presence of the tidal bias term breaks the degeneracy between the linear and quadratic terms, improving the predicted errors on b2.For both bias parameters, radio continuum surveys (mainly the high redshift bins) contribute enough signal in order to achieve a few percent precision measurements.However the introduction of theoretical errors deteriorates the expected constraints to a 10 − 20% precision at high redshift, mainly due to the uncertainty introduced by the exclusion of higher bias terms in the power spectrum and bispectrum of galaxies.

Redshift Space Distortions and redshift
uncertainties.
Both the effect of RSD and of uncertainties in the determination of redshifts must be fully taken into account in any realistic galaxy forecast, since galaxies are observed in this coordinate system.In the second step of our analysis, we start by comparing the idealised model (see Sec. 2.3.2) with the results derived from including the full RSD, up to second order (i.e. for the galaxy power spectrum and bispectrum we use Eqs.( 26) and ( 27) respectively).At this stage, we are still neglecting redshift errors.We also consider separately the effect of adding RSD and theoretical errors.The goal for this section is to specifically check the impact of including RSD and of going beyond angle averaged statistics, via a relative comparison with previous results.Actual, realistic forecasts, including redshift errors, will finally be presented in the next subsection.The final forecasts use the signal coming from the power spectrum, bispectrum and their combination.Results for bias and PNG parameters, as a function of redshift, are displayed in Fig. 12 and 13.The expected constraints on fNL, coming from the summed signal over all redshift bins, are reported in Table 5 and 6 (column marked by "RSD (no z-errors)").
The difference between the model including RSD only and the preliminary, idealised analysis is evident for the bias forecasts, especially the bispectrum derived ones.As we can see in Fig. 12 and 13, the effect of RSD alone is instead overall small for PNG, with the bispectrum again being affected the most.
As we can see in the two first columns of Tables 5 and  6, the forecasts on fNL become tighter for all three PNG cases when the full RSD treatment is applied (column denoted as "RSD (no z-error)") with respect to the results of the idealised case (column denoted as "Idealised").This improvement is expected, since the full RSD model contains additional signal compared to the approximating monopole [Eqs.( 31) and ( 32)], used in the idealised case.This indicates that the supplementary information contained in the higher bispectrum multipoles is important and it should be taken into account by future LSS surveys in order to provide tight constraints on the PNG amplitude.Radio continuum, 1 Jy Figure 12.Expected constraints on the linear and on the two quadratic bias parameters from the tree-level bias expansion and on the amplitudes of the PNG parameter f NL for each redshift bin in the case of the 1 µJy radio continuum survey when including RSD and FOG effects.The dashed lines represent the idealised case, with theoretical errors included, as described in Sec.6.1.1.The dashed-dotted lines represent the model which accounts for RSD and FOG effects, assuming perfect knowledge of redshifts, while the solid lines take also into account the redshift uncertainties.The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.The expressions used for the modelling of the power spectrum and bispectrum are given in Eqs. ( 26) and ( 27) respectively, while for the latter the trispectrum term is excluded.These results consider only the diagonal part of theoretical error covariance [Eq.( 41)] for the bispectrum, while for the power spectrum the full error covariance is used.MPTbreeze has been used to derive the matter two-and three-point correlators.
As a further step, in addition to RSD, we include theoretical errors in order to test the combined effect on the fNL expected constraints for the local, equilateral and orthogonal shapes.A comparison of the 1σ error bars on the non-Gaussian amplitude obtained respectively considering the idealised models, the RSD model and the RSD + theoretical errors model, without accounting for redshift errors, is displayed in Table 5 and 6.Here, the off-diagonal part of the theoretical error covariance matrix is ignored in the case of the bispectrum, due to the significant computational cost of inverting this matrix; off-diagonal elements should now be included for each orientation of the triangles [see Eq. ( 36)].Excluding the off-diagonal terms underestimates the mode coupling and the overall effect of theoretical errors.In order to quantify this effect, we compare forecasts with and without off-diagonal terms, for the computationally affordable idealised case.Since RSD corrections are not considered in the theoretical error treatment, the effect of the off-diagonal terms in the expected constraints can be transferred to the present forecasts.The marginal error provided by the bispectrum for the PNG amplitude, in the local case, increases in this test by ∼ 1%, if we include the full theoretical error covariance matrix for the initial redshift bins (see Sec. 6.1.1),instead of only the diagonal part.For the equilateral and orthogonal cases, we observe ∼ 10% and ∼ 5% en-hancements respectively, after marginalising over the other free parameters, where the degradation of the bias parameters, due to the inclusion of theoretical errors, is propagated on fNL.Note that these quantitative results are valid only for the radio surveys we consider here, since the effect of theoretical errors depends on redshift as well as the fiducial values of the bias parameters.With these caveats in mind, we see that the effect of theoretical errors follows a similar pattern to the one already discussed in the previous section, as expected.In the case of local non-Gaussianity, the impact of theoretical errors is small, while it is very large for the equilateral shape.
Having considered RSD, we now include redshift uncertainties, modelled like the FOG in Eqs. ( 29) and (30).Note that the power spectrum and bispectrum, including the stochastic bias terms, which at large scales resembles the Poisson shot noise, is multiplied by the damping factor, containing the redshift uncertainties, as shown in Eqs. ( 26) and ( 27)] respectively.
The forecasted constraints on PNG amplitudes and bias, after including the redshift errors, are presented as functions of redshift in Fig. 12 and 13.The final 1σ error bars on fNL are displayed in Table 5 and 6 under the column named "RSD+z-errors+Theoretical errors".In the case of local PNG, the effect of redshift uncertainties is small and  Table 5. Forecast results for the amplitude of PNG (f NL ) in the case of three different primordial shapes (local, equilateral and orthogonal), when considering respectively the idealised model, the model taking into account RSD and FOG effects and the redshift space model with theoretical errors, for the 1 µJy radio continuum survey.The redshift uncertainties are considered separately and added on top of all the previous effects under the column denoted "RSD+z-errors+Theoretical errors".The expected constraints on the bias parameters, up to quadratic order [Eq.( 17)], are also presented.The results are derived after marginalising over the unknown parameters.Note that the observed improvement in the P+B(equil) results, with respect to B(equil), originate from the tighter forecasts on b 1 provided by P(equil).
the expected constraints from both radio surveys are still tighter than those originating from Planck , for more or less the whole redshift range.The final error on f loc NL is degraded by a factor of ∼ 1.5 for a 1 µJy survey with respect to the case with a perfect redshift determination.This shows that future radio surveys can tightly constrain PNG of the local type, even with large redshift uncertainties, arising from the statistical nature of redshift determination.For the 1 µJy case, an error of σ(f loc NL ) ∼ 0.18 can be achieved from the combined P + B signal, i.e. an expected ∼ 30 times improvement with respect to the constraints of Planck .In the case of orthogonal PNG the degradation due to redshift uncertainties is much larger, as seen in Fig. 12 and 13.The final power spectrum + bispectrum forecasts for the orthogonal PNG amplitude, adding contributions from all redshift bins, show a ∼ 7 times deterioration, compared to the case with no z-errors, and is slightly worse than the current Planck constraint.
The degradation is even larger for PNG of the equilateral type.In this case, expected constraints provided by both radio surveys are far weaker than those coming from Planck , for the whole redshift range (dashed-dotted lines in Fig. 12 and 13).The predictions on the f equil NL parameter coming from all redshifts and the combined P and B signal get degraded by a factor ∼ 29, as seen in Table 5 and 6.The shape-dependence of this error degradation is due to the functional form of the damping factor, which takes a minimum value -and hence has the maximum effect -for the equilateral configurations, giving DFOG = exp[−3(kµσv) 2 ].The impact of the damping factor is instead minimized for squeezed triangles, on intermediate scales, and for folded ones, on large scales.
Therefore we can conclude that radio continuum samples, in combination with clustering-based redshift estimation, can provide tight constraints for the local PNG amplitude, with an important contribution from bispectrum measurements.In order to achieve tight constraints for other shapes, though, the precision in the determination of redshift should at least match the one achieved by photometric surveys (it is currently estimated to be about an order of magnitude worse).Theoretical errors, as seen in Table 5 and 6, are less relevant for the fNL predictions coming from surveys with large redshift uncertainties, since the effect of the latter overshadows completely the impact of the first.
The power spectrum redshift space forecasts presented here (i.e."RSD+z-errors+Theoretical errors" columns in Table 5 and 6) for the two radio continuum surveys are consistent with those presented in Raccanelli et al. (2017), after taking into account the different flux limits used in the latter.
The relative errors on bias parameters, shown in Fig. 12 and 13, coming from joint power spectrum and bispectrum estimation, increase by a factor of two when redshift errors are included.For the linear bias, the increase originates mainly from the deterioration of the power spectrum results.For the two quadratic bias terms, the relative errors become larger than unity, since the signal comes only from the bispectrum which is more affected by redshift errors.
Since redshift errors are the main limiting factor in our radio continuum forecasts, it is natural to consider alternative radio datasets in which such issue does not arise.To this purpose, it is interesting to look at HI spectroscopic surveys.The precise redshift information obtained in this case, together with the large volume of the survey, can in principle be very promising for PNG studies.The most interesting aspect is that very tight forecasts could be obtained not only for local-type, but also for equilateral-type NG, by adding the bispectrum and the trispectrum term to previous PNG power spectrum studies, using HI surveys, such as Camera et al. (2015b).A full bispectrum forecast for the HI SKA galaxy sample requires a dedicated study, and goes beyond the scope of this paper.It is the object of a forthcoming work.

The effect of the trispectrum term
In our final step, we include the trispectrum contribution in the galaxy bispectrum [Eq.( 22)], together with theoretical errors.This term will be considered for both the idealised and the full RSD models.The trispectrum term, discussed here, is present only for non-Gaussian initial conditions and exhibits a scale dependence, coming from the primordial bispectrum [Eq.( 25)].It can then be used to tighten the constraints on PNG amplitudes.Interestingly, this is limited not only to the local case, but it also includes the equilateral shape.The importance of the trispectrum contribution was originally investigated by Sefusatti (2009) and Jeong & Komatsu (2009) for local non-Gaussianity, but the term is included in an actual forecast here for the first time.The scale dependence induced in the galaxy bispectrum by the quadratic bias trispectrum correction is analogous to the scale dependent bias of the power spectrum.It can be calculated exactly for any type of PNG, without the need of using any kind of squeezed limit approximation, and 1/M (k) ∝ k −2 .Due to this, the degeneracy of f equil NL with the linear bias on large scales can now be broken (see the discussion in Sec.2.2.2).
The results for the PNG amplitude, after including the trispectrum contribution, as well as theoretical errors, are presented in Fig. 14, where we compare them to the expected constraints coming from the idealised, preliminary model.An improvement is observed for all PNG amplitudes and especially the expected constraints for the equilateral shape show an impressive enhancement.We have checked that the theoretical errors do not affect significantly σ(f equil NL ) in this case, since the scale dependent signal contribution compensates for the theoretical uncertainties.The bias forecasted constraints are unaffected, since the trispectrum correction disappears for chosen fNL = 0 fiducial value.The final predicted errors on the PNG amplitudes for the three cases considered here are presented in Table 7, after including the trispectrum term from the bias expansion.We see that, in  19) for the power spectrum and Eq. ( 22) for the bispectrum].The full theoretical errors covariance is considered here.The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.
principle, the inclusion of trispectrum corrections allows for significant improvements for all shapes, including equilateral.
One concern is that the two trispectrum corrections, i.e.T M P T 1112 and S2T M P T 1112 , depend both on fNL and on bias parameters b1, b2 and bS 2 .Therefore, they could generate degeneracies between PNG and bias terms, which is indeed the case.However, we explicitly checked that, for large volume surveys such as the 1 µJy radio continuum survey considered here, where very large scales are included, such degeneracies are broken.
The extension of this formalism to redshift space, by using the derived expression in Eq. ( 27) with the trispectrum correction [Eq.( 28)], is performed here for the two radio continuum surveys.This model will be considered as our "full" treatment case.The resulting forecasts of the PNG amplitude are presented for each redshift bin in Fig. 15.Forecasted constraints for the integrated signal over the whole redshift range are presented in Table 8.The modelling for both correlators is described in detail in Sec.6.1.2(i.e.Eqs. ( 26) and ( 27) for power spectrum and bispectrum respectively), where the effects of redshift errors and FOG are taken into account.Note that, only the diagonal part of theoretical error covariance is used in these redshift space forecasts (see Sec. 6.1.2for a discussion).
As seen in Fig. 15 and Table 8, the improvement provided by the trispectrum term in the forecasts of f loc NL is negligible, while for the orthogonal PNG type the contribution is minimal.The constraining power of the trispectrum term is reduced when the RSD effect is taken into account, contrary to the idealised, preliminary analysis (Table 7).This can be mainly attributed to the presence of redshift errors.On the other hand, the scale dependence provided by the trispectrum term still significantly improves the final fore-  26) and ( 27) respectively, while for the latter the trispectrum term is excluded.The solid lines represent the same model when the trispectrum correction [Eq.( 28)] is taken into account.Redshift errors and FOG effects are considered for both cases, while we consider only the diagonal part of theoretical error covariance [Eq.( 41)].The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.
casted constraints for the equilateral PNG type, in combination with clustering-based redshift estimation.More precisely, the forecasts coming from the redshift space bispectrum are tighter by a factor of ∼ 2.5, while for the combined signal from the two correlators the improvement can reach up to a factor of ∼ 2.
The "RSD+trispectrum" results, presented in Table 8, correspond to our "full" treatment and they are considered as the final forecasts on fNL for the three PNG types provided by the radio continuum samples.These forecasts show a significant expected improvement with respect to current f loc NL Planck constraints, while for orthogonal PNG the 1σ forecast results are on the same level as Planck.For the equilateral type, the forecasts show worse expected constraining power with respect to what already achieved by Planck.Nevertheless, we demonstrated the importance of the trispectrum bias term in improving the expected constraints provided by large volume LSS surveys, even in the case of large redshift uncertainties.If large volume surveys (either optical or radio) will allow for accurate redshift measurements, at some stage, not only local models will be measured with high sensitivity, but also all other shapes.This will be further discussed in the optical survey forecasts to follow (Sec.6.2).
6.1.4Non-Gaussian corrections to the bispectrum variance.
In this work, as discussed in Sec.4.1, only the diagonal part of the bispectrum covariance is used in the Fisher matrix formalism.In addition, for the variance we use the  31) and ( 32)], while the full effect of the theoretical errors is taken into account here.The modelling used for the power spectrum is given by Eq. ( 19) and for the bispectrum by Eq. ( 22), where for the latter the trispectrum bias term is taken into account.Table 8.Forecast 1σ results for the local, equilateral and orthogonal PNG in redshift space.These results consider FOG effects and redshift errors are considered, while only the diagonal part of theoretical error covariance is taken into account here.The modelling used is for the power spectrum Eq. ( 26) and for the bispectrum is given by Eq. ( 27), where for the latter the trispectrum bias term [Eq.( 28)] is taken into account.The RSD model without the trispectrum contribution (i.e."RSD") is shown for comparison purposes.
predictions of PT up to tree-level.Therefore, at this point, we would like to test the effect of excluding higher-order corrections.Using Eq. ( 40), we test the effect on the expected fNL constraints coming from the 1 µJy radio continuum, and in particular those originating from the real-space model [Eq.( 22)] as well as the redshift space bispectrum [Eq.( 27)], following the procedure outlined at the end of section Sec.4.1.For local PNG, expected constraints in real space deteriorate by ∼ 39%, while in redshift space the deterioration increases to ∼ 67%.For equilateral PNG, the effect seems to be smaller, a degradation of ∼ 12% and ∼ 17% is observed for the real and redshift space case, respectively.This should be kept in mind when quoting final results.A full Fisher matrix analysis, including full covariances, rather than the simplified estimates provided here, will be object of a future study.

Optical surveys
In this section we show forecasts for a spectroscopic and a photometric survey (see Sec. 3.2 for details on the specifications used).These will be ultimately compared with those originating from the two radio continuum surveys (see Sec. 3.1), presented in the previous section.Additionally, these surveys will allow us to test the full effect of theoretical uncertainties, mainly because their smaller volumes reduces significantly the computational time.We will in particular include off-diagonal theoretical error terms in our preliminary, idealised forecasts, up to the highest redshift bins (for nearly all scenarios; few exceptions will be pointed out case by case).
The optical survey analysis will broadly follow the same scheme as adopted in the previous section, based on adding realistic features and higher order corrections step by step, on top of the initial idealised model, in order to check separately their impact.

Idealised treatment
As before, we start with the idealised case and show the effect of adding theoretical errors.The full error covariance is here used only up to the 6th bin for the photometric case.For the remaining bins, only the diagonal contribution is considered, since for the high redshift bins the off-diagonal Table 9. Forecast 1σ results for the three PNG type considered here, originating from a spectroscopic and a photometric survey.The idealised model is as usual specified by Eqs. ( 31) and ( 32), and we show the full effect of including theoretical errors.
terms have a minimum effect on the final PNG amplitude forecasts, as discussed extensively in Sec.6.1.1.The full theoretical error covariance will instead always be used over the whole redshift range for the power spectrum, as done before.
Besides discussing theoretical errors, we also show the effect of the inclusion of trispectrum corrections.The marginalised fNL forecasts in both cases are shown in Fig. 16.The effect of the full theoretical error covariance is presented in Table 9 for the two optical surveys.As done also for the radio continuum surveys, the practicality of studying the idealised case first lies in the fact that a quantification of the effect of the full theoretical error covariance is needed, in order to propagate it to the RSD treatment (next section), where the computational effort is large.It is evident that the behaviour follows the same pattern as in the case of radio continuum.A degradation in the final PNG forecasts is observed ranging between 25 − 50% depending on the PNG type, as well as the size and redshift range of the survey (see Table 9), where for the equilateral case the effect is maximum.Excluding the off-diagonal elements of the theoretical error covariance introduces an almost negligible error for local PNG for all optical surveys, while for the equilateral and orthogonal forecasts an underestimation of the theoretical uncertainties by 9% and 4% is observed respectively.The forecasts for all thee PNG types coming from the photometric survey are improved with respect to the spectroscopic one.
Adding trispectrum corrections produces some improvements.For local PNG, low redshift bins produce an order of magnitude improvement, while high redshifts have a much smaller effect.The same behaviour is also observed for the orthogonal shape, but in this case the improvement is much smaller (∼ 2.5 times).For equilateral PNG, a small improvement is observed only for the very low redshift bins, while for the large redshifts the expected constraints become even worse.The lack of improvement observed in the small volume spectroscopic survey is due to degeneracies between fNL and bias parameters, which are enhanced by the presence of the trispectrum term.As seen in Fig. 17, f loc NL does not show any degeneracies with the three bias terms (i.e.b1, b2 and b s 2 ), considered here as free parameters, while f orth NL shows some degeneracy mainly with b2.On the other hand, f equil NL shows strong degeneracy with all three and hence this explains the observed degradation in results (see Fig. 16 and Table 10) after the inclusion of the trispectrum term.Such degeneracies can be broken if large enough scales are avail-  19) for the power spectrum and Eq. ( 22) for the bispectrum).The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross-covariance.
able in the survey (see Fig. 14 for the radio continuum cases).For example, an overall gain in the forecasts provided by the photometric survey can be observed, since they probe a larger part of the sky and hence grant access to larger scale modes.This fact will partially break the degeneracies of fNL with the bias parameters, providing the observed improvement over the results of the idealised model.More specifically, for the local-type PNG, the trispectrum corrections reduce the 1σ expected constraints on f loc NL originating from the photometric survey beyond unity for the high redshift bins.For the equilateral PNG case, a compelling improvement is observed for both surveys.
Here we would like to point out the importance of the bispectrum for the PNG amplitude forecasts coming from LSS surveys.As we can see in Fig. 16, the gain in the final fNL expected constraints of the local type is small when the signal from the bispectrum is added on top of the one generated by the power spectrum.However, this changes dramatically for the equilateral and orthogonal types, where the bispectrum is the main source of signal.
We present the marginalised idealised forecasts for the PNG amplitude in Table 10.We see that in principle the trispectrum corrections can lead to significant improvements for large volume photometric surveys.In total analogy with the radio continuum case, the introduction of redshift errors at the expected level for next generation surveys will however wash away the enhancement.We explicitly show this in the next section.
All the previous results, assume a Gaussian and diagonal bispectrum covariance matrix.By using Eq. ( 40), proposed in Chan & Blot (2017), we can effectively resum the contributions to the diagonal coming from higher order terms.Doing so will degrade the results presented in Table 10.In the case of the spectroscopic survey a ∼ 35% and ∼ 16% increase is observed for the 1σ forecast errors in the case of local and equilateral PNG type respectively.For the photometric survey the deterioration of the forecasted constraints on local and equilateral PNG is ∼ 29% and ∼ 10% respectively.The same trend seen for the radio continuum surveys is also observed here, i.e. the local expected constraints are highly affected by the higher order corrections in the bispectrum variance, while a minimum effect is observed for the equilateral PNG case.

Full redshift space treatment
Our next step is to move to redshift space, considering RSD, redshift errors and the FOG smearing effect.As usual, for the sake of comparison, we will use the bispectrum RSD model with and without the trispectrum bias term given in Eq. ( 28).Regarding theoretical errors, only the diagonal part of the error covariance [Eq.( 41)] will be used in redshift space, due to the high computational cost.The results are presented as a function of redshift bins in Fig. 18 for the spectroscopic and photometric surveys.These surveys can produce improvements on current Planck bounds for the local shape, but not for the other two, mostly due to theoretical errors.
Adding trispectrum corrections has a minimal effect for spectroscopic surveys, because such corrections need very large scales to be effective and the probed volumes are here too small to make a difference.
The situation changes a bit for the photometric survey.This survey probes a much larger volume, therefore trispectrum corrections can in principle play a useful role here, as shown explicitly in the previous section.At the same time, however, we are now affected by large photometric redshift errors.Due to the latter, all PNG predicted constraints are now degraded (Fig. 16) with respect to the idealised forecasts, where redshift uncertainties are not considered.The full redshift space forecasts, coming from the integrated signal over the entire redshift range, are presented in Table 11 for the galaxy power spectrum, bispectrum and their combination.Both cases, with and without the trispectrum term, are shown.To better quantify the effect of redshift uncertainties on our final fNL predicted constraints, we present also the results of the RSD model, taking σz(z) = 0.The local PNG case is affected the least by redshift uncertainties and therefore the observed improvements in the f loc NL pre-  31) and ( 32)], after considering in the galaxy bispectrum the trispectrum bias term shown in Eq. ( 22).
dicted constraints for the photometric survey are justified.Its constraining power is significantly better than that of Planck, for bins with z 1. Adding the trispectrum correction has a minimal impact on the bispectrum forecasts for the orthogonal and local types, while it can have an overall improvement in the equilateral results of the large-volume surveys.Still, large redshift errors make the trispectrum contribution far from enough to improve over current bounds.As we show in Sec.6.2.1, the effect of neglecting offdiagonal terms from the theoretical error covariance is expected to be small for surveys with wide redshift range.In order to quantify this, we performed a full error covariance analysis for the spectroscopic survey.The survey volume size shrinks in this case the computational effort tremendously, compared to radio continuum cases, and makes a full numerical analysis feasible.The deterioration level in the final PNG expected constraints is consistent with what reported for radio continuum.More specifically, we observe a < 1% and ∼ 7% degradation for local and equilateral PNG cases respectively.
The "RSD+trispectrum" results in Table 11 are considered as our final forecasts on the three PNG types originating from the optical surveys.Comparing these results with those available in the literature, we see that our power spectrum spectroscopic results, after excluding the theoretical errors, are consistent with the forecasts in Giannantonio et al. (2012) for the local and orthogonal PNG.On the other hand our forecasts on the local type, coming from the spectroscopic survey, are slightly worse than those reported in Tellarini et al. (2016), where the difference originates from the fact that we take into account the presence of theoretical errors, as well as the effect of redshift uncertainties.The  26) and ( 27)].Three different versions of the model are examined here: without considering the redshift uncertainties ("RSD (no z-errors)+Theoretical errors"), including the redshift errors ("RSD+z-errors+Theoretical errors") and adding to the latter the trispectrum contribution ("RSD+z-errors+Theoretical errors+trispectrum"). Regarding the theoretical errors, only the diagonal part of the covariance is used for all versions.26) and ( 27) for the power spectrum and bispectrum respectively, without the presence of the trispectrum term [Eq.( 28)] in the galaxy bispectrum.Solid lines represent the same model but with the trispectrum correction included.The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.The FOG effect, as well as the redshift error, is taken into account in both models.Only the diagonal part of the theoretical errors covariance is used here.
differences are slightly higher for redshift space bispectrum forecasts, owing to the higher impact of theoretical errors and to the inclusion of redshift errors (small, but not completely negligible also for spectroscopic surveys).
As done for the radio analysis, we also tested the effect of neglecting non-Gaussian corrections in the bispectrum variance, by using Eq. ( 40).In the case of the spectroscopic survey, a degradation of ∼ 42% and ∼ 20% is observed in the forecasts for the local and equilateral PNG types respec-tively.For the photometric probe we get a deterioration of ∼ 32% and ∼ 10% for the final predicted constraints on local and equilateral PNG respectively.This is not overall a very large effect, in the context of a Fisher forecast.However, it does make it clear again that such NG contributions to the bispectrum variance are not negligible, and future analyses will have to include them to achieve high accuracy.

CONCLUSIONS
In this paper we have investigated possible future constraints on the amplitude of the non-Gaussian parameter fNL for three types of PNG shapes -local, equilateral and orthogonal -and on galaxy bias parameters, through galaxy power spectrum and bispectrum measurements on large scales using a Fisher matrix approach.We thoroughly accounted for a large number of effects in modelling the gravitational non-Gaussian contributions, including a full second order treatment of bias and RSD, going beyond real space angle averaged fNL forecasts presented in other works.We carefully investigated the propagation of theoretical uncertainties, following the approach introduced in Baldauf et al. ( 2016), but extending it to bias loop-corrections.All these effects were to a various extent included in previous literature, but for the first time we accounted for all of them at once in a single forecast analysis.The cross-correlation between power spectrum and bispectrum was ignored in this work, and the bispectrum covariance approximated as diagonal.However, for the large scales considered here, we discussed how this should have a small impact, based on recent results in the literature (Chan & Blot 2017).We also employed, as standardly done, a Gaussian approximation for the bispectrum variance.We presented an explicit estimate of the effect of ignoring non-Gaussian contributions to the variance, by considering leading non-Gaussian corrections.We found that, in the worst case scenario, such corrections can degrade constraints by ∼ 50%.A more detailed study of the bispectrum covariance, including PNG corrections, will be included in a future work.Likewise, it will be important Table 12.Summary of 1σ limits for the three PNG types considered, from radio continuum and optical surveys derived from combining the power spectrum and bispectrum and accounting for RSD, the trispectrum term and theoretical errors.See text for more details.
to account in the future for the effect of relativistic corrections on the bispectrum, especially for large volume surveys (see e.g.Di Dio et al. 2016;Raccanelli et al. 2016;Umeh et al. 2017;Bertacca et al. 2018).
In addition to the previous ingredients, we improved the modelling of the galaxy bispectrum by considering a complete second order bias expansion, which includes for the first time the trispectrum term [Eq.( 17)].We only consider the zeroth order (tree-level) expansions in the matter fields, because we are only interested in the large-scale contributions.For dark matter, we have used the MPTbreeze perturbation theory, based on Renormalised Perturbation Theory, which provides a natural cut-off in the non-linear regime, such that ultraviolet divergences are automatically removed.The final bispectrum model used for our forecast is presented in Eq. ( 22).
Our results are based on radio continuum (with 10 µJy and 1 µJy flux limits) and optical (spectroscopic and photometric) surveys.In the case of radio surveys, we assume a clustering-based estimate of redshifts.
We have summarised our main results in Table 12, where we have reported the forecasts derived by combining power spectrum and bispectrum for the three non-Gaussian shapes, together with the Planck temperature and polarisation constraints (Planck Collaboration et al. 2016b), as a comparison.
For the local shape, LSS measurements can provide important improvements over current CMB bounds, even using only the power spectrum, but the bispectrum will allow doing much better, with a further improving factor ∼ 5.For equilateral PNG, the bispectrum is by far the most useful statistics, but it is plagued by the fact that theoretical errors peak in the equilateral limit.Therefore, improving upon current forecasts by including smaller scales and more modes, at low redshift, will require exquisitely accurate modelling of late-time non-linearities.On the other hand, large future radio continuum surveys will provide access to much larger volumes and higher redshifts.
The trade-off for these surveys is however represented by large errors in the determination of redshifts, which become dominant and now overshadow the effect of theoretical errors.In a very idealised case where galaxy redshifts could be accurately known for all objects, even in presence of significant theoretical errors, large improvements (e.g. up to a factor ∼ 5 for radio) with respect to Planck equilateral constraints are obtained.This constraining power mostly comes from very large scales and trispectrum contributions, which display a ∼ k −2 scale-dependence in the equilateral case.Such contributions therefore deserve further attention.In a more realistic case, redshift errors affect these forecasts significantly.The gain obtained by including the trispectrum term is actually still large when z-errors are included, but it becomes insufficient to improve over current equilat-eral Planck bounds.The final forecasts are in the first two columns of Table 12.
We note that the issue above cannot be circumvented simply by looking at the trispectrum term in fortchoming optical spectroscopic surveys, because they have too small a volume to make such term significant.We argued that looking instead at HI spectroscopic surveys is another interesting approach, which we will pursue in a forthcoming work.The other possible way forward would be finding strategies for better determination of redshifts in either future photometric or radio continuum surveys.
With this formalism, one can calculate the n-point correlation function in RPT for an arbitrary number of loops, but the actual computations are difficult because they involve solving numerically a set of integro-differential equations.Nevertheless, this method provides a well-defined perturbative expansion in the non-linear regime, which is not the case in SPT.

APPENDIX C: DERIVATION OF BIAS PARAMETERS WITH THE PEAK-BACKGROUND SPLIT METHOD
The peak-background split approach relates the bias parameters to the mean tracer abundance, where the density fluctuation field is decomposed into a low-amplitude longwavelength linear fluctuation and a noisy short wavelength one, δ h (x) = δ h,l (x) + δ h,s (x).The short-wavelength fluctuations ride on top of the large-wavelength ones, where we treat the matter perturbations as a superposition of small and large modes separated by a cut-off wavenumber.The short wavelength modes are the sources of the dark matter halos, while the long-wavelength modes increase or decrease the background density in large patches of the sky.In this picture the peaks of the small modes that are located over the peaks of long-wavelength modes will be clustered more than the average ones, as well as being the first to collapse forming galaxy clusters.This explains why galaxy clusters are more clustered than the galaxies themselves.
The presence of the long-wavelength linear fluctuation from the Gaussian case is that now the long and short wavelength modes are coupled to each other.For the local case, after splitting the Gaussian part of the primordial gravitational potential into long and short wavelength fluctuations, ΦG = φ l + φs and substituting into the local expansion of the non-Gaussian potential field Φ we get Φ = φ l + fNLφ 2 l + (1 + 2fNLφ l )φs + fNLφs. (C16) The most important term here is the coupling term, (1 + 2fNLφ l )φs, between long and short modes, since the long wavelength linear fluctuations will introduce a scale dependence rescaling in the amplitude of the short modes.In the case of a general non-local non-Gaussianity this rescaling of long-wavelength part of the initial density fluctuations can be parametrized through where the is an infinitesimal parameter, which becomes = fNLφ l for local non-Gaussianity.The modulation in the primordial large wavelength density mode will affect the variance of the small scale modes introducing additional dependences in the number of collapsed objects.The short wavelength variance will transform at the lowest order to The Jacobian ,J ≡ dlnν dlnM , will also be transformed for the non-local case into (Desjacques et al. 2011b,a) The peak-background split argument of Eq. (C6) can be easily generalized to where the average halo number density n h is taken after substituting in Eq. (C1), Eq. (C18) and Eq.(C19).The leading non-Gaussian bias bΨ can be derived from the above relation as a special case for N = 0.The local non-Gaussian result can be then calculated from it for α = 0, where for a mass function with a universal form we get which yields the well-known formula for the scale-dependent bias, after being plugged into Eq.( 20).
The first higher order non-Gaussian bias parameter b Ψδ can be easily derived through Eq. (C20) for N = 1 as Following the same steps as for the local-in-matter bias parameters, we relate the Lagrangian bias to the desired Eulerian one throughout Eq.(C11) as (Giannantonio & Porciani 2010;Baldauf et al. 2011) Combining the above equations with Eq. (C22) and the Eulerian results for the local-in-matter bias parameters we can derive Eq. ( 11).The non-Gaussian mass function considered here is the LoVerde mass function, where the fractional correction RNG will introduce a scale-independent offset in the bias parameters originating from the partial derivative in Eq. (C20).These terms will depend on fNL due to the presence of a non-zero skewness S3 in the Edgeworth expansion of the LV mass function.These corrections will be: and
Figure1.The shape of the galaxy bispectrum for Gaussian initial conditions, together with its highest contributing terms.In each panel the bispectrum, normalized to its absolute maximum value, is plotted as a function of k 2 /k 1 and k 3 /k 1 for all configurations in the case of fixed k 1 .We consider the following three values: k 1 = 0.01, 0.05, 0.1 h/Mpc, where the triangle sides satisfy k 3 k 2 k 1 .In the first row the Gaussian galaxy bispectrum B(G) without the trispectrum contribution is plotted [Eq.(21)], in the second we plot the tree-level matter bispectrum (B G ) as predicted by SPT, in the third the quadratic bias term indicated as P 2 is plotted (i.e.P L m (k 1 )P L m (k 2 )+2 perm) and finally in the last row we plot the tidal bias term contribution (i.e.B S 2 = S 2 (k 1 , k 2 )P L m (k 1 )P L m (k 2 ) + 2 perm).Note that for the latter the absolute value is plotted, where a white dotted line shows the separation between the positive (left side) and negative (right side) values.For a detailed explanation, see the main text (Sec.2.3.1).

Figure 3 .
Figure3.Same as Fig.2but for the equilateral type of PNG.For this type of PNG the galaxy bispectrum is taken to be that of Eq. (21) (i.e.B N G is here the sum of the trispectrum bias corrections), since any term proportional to the field Ψ in Eq. (22) (introduced to model the scale dependent bias corrections) is excluded as we discuss in Sec.2.2.2.

Figure 4 .
Figure 4. Same as Fig. 2 but for the orthogonal type of PNG.In this case the equations used are the same as for the local PNG.
Figure5.The galaxy bias, defined from the weighted average over the halo one, as a function of redshift for the two radio surveys considered here, a survey flux limited at 10µJy and one with a 1µJy flux limit.The details of the surveys are listed in Table1.

Figure 6 .
Figure6.The overall true redshift distribution n(z) (solid black line) is plotted over the redshift range considered for the photometric survey.The redshift distribution for each tomographic bin, given by Eq. (33), is also shown (dashed, dotted, etc. lines).

Figure 7 .
Figure 7.The linear and quadratic galaxy bias for the optical surveys considered here.The details for the derivation are described in Sec.3.2 for each individual survey.

Figure 9 .
Figure 9. Same as in Fig. 8, but for the 10 µJy radio continuum survey.

Figure 11 .
Figure11.Same as Fig.9, where here the ratio of the 1σ errors between the case with theoretical errors (denoted in the plots with the upper index TE) and without is shown.

Figure 13 .
Figure 13.Same as Fig. 12, but for the 10 µJy flux limited radio continuum survey.

Figure 14 .
Figure14.Fisher forecasts on f NL for local, equilateral and orthogonal shapes in the case of a 1 µJy survey (left) and a 10 µJy on (right).The dashed lines represent the starting, idealised model, while the solid lines represent the same model when the trispectrum correction is considered [i.e.Eq. (19) for the power spectrum and Eq.(22) for the bispectrum].The full theoretical errors covariance is considered here.The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.

Figure 15 .
Figure15.Fisher forecasts on f NL for local, equilateral and orthogonal shapes in the case of a 1 µJy survey (left) and a 10 µJy on (right).The dashed lines represent the model for the power spectrum and bispectrum in redshift space, given by Eqs.(26) and (27) respectively, while for the latter the trispectrum term is excluded.The solid lines represent the same model when the trispectrum correction [Eq.(28)] is taken into account.Redshift errors and FOG effects are considered for both cases, while we consider only the diagonal part of theoretical error covariance [Eq.(41)].The expected constraints from the power spectrum are plotted in blue, from the bispectrum in green and from both power spectrum and bispectrum in red, without taking into account their cross term in the covariance.

Figure 17 .
Figure17.The joint 1σ forecasts on the three galaxy bias parameters (i.e.b 1 , b 2 and b s 2 ) and PNG from the z = 1.8 redshift bin of the spectroscopic survey.In each panel the two-dimensional joint forecasts are shown for all the combinations between the parameters, after marginalising over the remaining.The model used is the approximating monopole for both power spectrum and bispectrum [Eqs.(31) and (32)], after considering in the galaxy bispectrum the trispectrum bias term shown in Eq. (22).

Figure C1 .
Figure C1.The Eulerian halo bias parameters up to fourth order redshift z = 0 as a function of the halo mass.

Table 1 .
The basic specifications for the two radio surveys considered here for each redshift bin.The shell volume is in units of (Gpc/h) 3 and the mean number density in 10 −4 (h/Mpc) 3 .The redshift distributions are obtained with the CBR technique (see Sec. 3.1 for details).We take the redshift uncertainty of sources as the same value of half the width of the bin it resides in.

Table 2 .
as a function of redshift.The basic specifications for a spectroscopic survey for each redshift bin.The shell volume is in units of (Gpc/h) 3 and the mean number density in 10 −4 (h/Mpc) 3 .

Table 3 .
The basic specifications for a photometric survey for each redshift bin.The derivation of these values is described in Sec.3.2.The shell volume is in units of (Gpc/h) 3 and the mean number density in 10 −3 (h/Mpc) 3 .

Table 6 .
Same as Table5, but for the 10 µJy flux limited radio continuum survey.

Table 7 .
Forecast 1σ results for the local, equilateral and orthogonal PNG.The calculations are performed by using the approximating monopole [Eqs.( The idealised model is shown for comparison purposes.

Table 10 .
Forecast 1σ results for the three PNG type considered here, originating from a spectroscopic and a photometric survey.Left column: idealised case with + theoretical errors [Eqs.(31) and (32)].Right column: idealised + theoretical errors + trispectrum corrections .The full effect of the theoretical errors is taken into account here, as described at the beginning of Sec.6.2.

Table 11 .
Forecast marginalised 1σ results for the three PNG type considered here, originating from a spectroscopic and a photometric survey.The full RSD model for both power spectrum and bispectrum is considered [Eqs.(