Projected two- and three-point statistics: Forecasts and mitigation of non-linear RSDs

The combination of two- and three-point clustering statistics of galaxies and the underlying matter distribution has the potential to break degeneracies between cosmological parameters and nuisance parameters and can lead to significantly tighter constraints on parameters describing the composition of the Universe and the dynamics of inflation. Here we investigate the relation between biases in the estimated parameters and inaccurate modelling of non-linear redshift-space distortions for the power spectrum and bispectrum of projected galaxy density fields and lensing convergence. Non-linear redshift-space distortions are one of the leading systematic uncertainties in galaxy clustering. Projections along the line of sight suppress radial modes and are thus allowing a trade-off between biases due to non-linear redshift-space distortions and statistical uncertainties. We investigate this bias-error trade-off for a CMASS-like survey with a varying number of redshift bins. Improved modelling of the non-linear redshift-space distortions allows the recovery of more radial information when controlling for biases. Not modelling non-linear redshift space distortions inflates error bars for almost all parameters by 20%. The information loss for the amplitude of local non-Gaussianities is smaller, since it is best constrained from large scales. In addition, we show empirically that one can recover more than 99% of the 3D power spectrum information if the depth of the tomographic bins is reduced to 10 $h^{-1}$Mpc.


INTRODUCTION
Our understanding of the Universe has been shaped by observations of the cosmic microwave background (CMB) and large-scale structure of the Universe (LSS) over the last thirty years (Planck Collaboration et al. 2018;Riess et al. 2016;Alam et al. 2017). Future CMB lensing experiments are expected to refine this picture (Ade et al. 2019;Abazajian et al. 2019). Upcoming galaxy surveys that trace the matter distribution of the LSS such as LSST (LSST Science Collaboration et al. 2009), SPHEREx (Doré et al. 2014), Euclid (Amendola et al. 2018) and DESI (DESI Collaboration et al. 2016) will contribute complementary information about the late time evolution of the Universe. Moreover, those surveys will achieve exquisitely small statistical errors due to the vast volumes they cover and high number density of tracers they resolve. Given their three-dimensional origin, upcoming LSS datasets are predicted to eventually contain more information about cosmological parameters than the CMB.
The early Universe's density distribution was very close to a Gaussian random field (Planck Collaboration et al. 2019) which is fully described by the two-point correlation function or its Fourier E-mail: ol248@cam.ac.uk transform, the power spectrum. The subsequent non-linear evolution changed the matter distribution which manifests itself in a modification of the power spectrum on small scales and non-vanishing higher order correlation functions. The bispectrum, which is the Fourier transform of the three-point-correlation function, is known to contain most of the non-linear information on mildly non-linear scales. In addition, it allows us to break degeneracies between bias and amplitude parameters (Scoccimarro 2000;Sefusatti et al. 2006).
CMB-lensing captures the integrated effect of matter onto CMB photons along their path from the surface of last scattering through the LSS to us. Accordingly, lensing spectra are described by projected spectra of the LSS. But there are also very valid reasons to study galaxy clustering statistics in projection. Firstly, tomographic surveys infer the redshift bins of objects and not their precise positions. Moreover, CMB lensing -galaxy clustering cross-correlations require a two-dimensional clustering analysis. Lastly, projections offer a way to suppress non-linear redshift-space distortions (RSDs).
RSDs are generated by galaxies' peculiar velocities parallel to the line of sight (LOS). The resulting Doppler redshift is degenerate with the redshift from the Hubble flow which is used to determine the radial positions. On large scales, the effect is well described perturbatively, but on smaller scales one has to resort to empirical models. Work has been done in this direction (Scoccimarro et al. 1999;Taruya et al. 2010;Gil-Marín et al. 2014;Slepian & Eisenstein 2017) but the fundamental issue of potentially biased estimates caused by inaccurate modelling remains.
In this work we are quantifying the parameter shifts due to inaccurate RSD modelling. In particular, we study how those biases depend on the chosen projection depth and RSD model used. Given the exquisitely small statistical errors of upcoming surveys, it is worthwhile to make estimators more robust in order to confidently leverage the small statistical uncertainties. Agarwal et al. (2020) recently investigated the parameter shifts arising from an incomplete or incorrect account of bias parameters and selection effects.
In addition, we forecast error bars to investigate the relation between statistical and systematic uncertainties. The constraining power of the galaxy bispectrum for future galaxy surveys has been studied in Yankelevich & Porciani (2019); Karagiannis et al. (2018); Agarwal et al. (2020) and power spectrum forecasts for CMB-lensing -galaxy-clustering cross-correlations were performed in Schmittfull & Seljak (2018).
While we employ the flat sky approximation in this work, there is a growing literature that studies angular (cross-)correlation functions (Schöneberg et al. 2018;Grasshorn Gebhardt & Jeong 2018;Simonović et al. 2018;Slepian 2018;Fang et al. 2020;Slepian et al. 2019;Campagne et al. 2017b;Tansella et al. 2018a,b;Campagne et al. 2017a;Moradinezhad Dizgah et al. 2020). The step from flat to curved sky is conceptually straightforward in our framework using the FFTlog-algorithm (Assassi et al. 2017), but we leave this for future work.
This paper is structured as follows: We first introduce the power spectrum and bispectrum for matter and galaxy (cross-)correlations in section 2. We then project those into observable 2D spectra in section 3. In section 4 we discuss the inference techniques used in our analysis. Our results are presented in section 5 and we conclude in section 6.

STATISTICS IN 3D
In this section we review the leading order power spectra and bispectra at late times. Starting with the matter predictions, we subsequently include galaxy biasing, RSDs and primordial non-Gaussianities (PNGs) into the model. We conclude this section with an overview of all matter-galaxy power spectra and bispectra.

Matter power spectrum and bispectrum
The linear matter power spectrum, P mm , can be efficiently computed using Boltzmann codes such as (Lewis et al. 2000;Challinor & Lewis 2011). 1 Using the matter power spectrum as input, the tree-level matter bispectrum is given by (Bernardeau et al. 2002) B mmm (k 1 , k 2 , k 3 ) = 2P mm (k 1 )P mm (k 2 )F 2 (k 1 , k 2 ) + 2 perm. , (1) where the second order gravitational kernel, F 2 , is with µ being the cosine of the angle between the two vectors k 1 and k 2 . The second order kernel for the velocity divergence is given by 1 https://camb.info/ G 2 will be relevant for the perturbative description of RSDs.

Biasing
Galaxies do not directly trace the underlying matter distribution, leading to a so called bias relation between galaxy and the matter distribution. At large scales the bias relation can be described perturbatively, for a review see Desjacques et al. (2018). Following McDonald & Roy (2009); Baldauf et al. (2012); Chan et al. (2012); Lazeyras & Schmidt (2018); Abidi & Baldauf (2018), we express the galaxy over-density at late times as a function of the three bias parameters (b 1 , b 2 , b s 2 ) together with some stochastic bias or shot-noise, , caused by the discrete nature of galaxies δ g (k) = b 1 δ(k) + (k) We truncated the expansion at second order, because we will be working with the tree-level power and bispectrum. The operator generating the Fourier representation of the square of the tidal tensor, S 2 , is given by Following Scoccimarro et al. (2001); Schmidt (2016); Desjacques et al. (2018), we model the shot-noise as Poissonian, which leads to the following non-zero spectra P = 1/n, wheren is the co-moving average number density of the tracers and all correlators between stochastic and matter densities fields are vanishing. Thus, the leading order galaxy power spectrum and bispectrum are given by and

Redshift space distortions
Galaxies' radial distances are measured via their redshifts. However, peculiar velocities parallel to the LOS give rise to a Doppler redshift that is degenerate with the cosmological redshift and does bias distance measurements. On large scales, RSDs are caused by the coherent infall of galaxies into gravitational potentials and lead to an enhancement of modes parallel to the LOS. This effect can be treated perturbatively and incorporated into the spectra by means of redshift kernels Z i that replace the gravitational kernels F i . The first two are given by Scoccimarro et al. (1999) where f = d ln D/d ln a refers to the logarithmic growth rate, the logarithmic derivative of the growth factor D. µ i is the cosine of the angle between the wave vector and the LOS, i.e. µ i = k i, /k i .
Here, k 2 i j = (k i + k j ) 2 and µ i j k i j = µ i k i + µ j k j . 2 On small scales, the high and incoherent velocities of galaxies within potential wells, give rise to the Finger-of-God (FoG) effect, i.e., structures appear elongated along the LOS. This corresponds to a damping of smallscale modes along the LOS. As it is intrinsically non-linear, a perturbative treatment is not possible and phenomenological modelling is needed (Scoccimarro et al. 1999). For this work we follow Peacock & Dodds (1994); Ballinger et al. (1996) and use a Gaussian damping prefactor for the perturbative predictions for the power spectrum and a similar form for the bispectrum We treat σ P and σ B as free parameters whose values are inspired by N-body simulations. Physically speaking, they correspond to the velocity dispersion of the galaxy velocity distribution. The general definition of (12) will shorten notation later in the paper.

Primordial non-Gaussianities
The assumption of Gaussian initial conditions can be tested by allowing small deviations and constraining their amplitude f NL . This is typically done by adding primordial non-Gaussianities (PNGs) of known shape. Those give rise to two types of new terms: Firstly, the added PNGs lead to a non-zero matter bispectrum at all times. Secondly, the bias expansion (4) has to be carried out both in the density and gravitational potential adding new bias terms (Dalal et al. 2008;McDonald 2008;Verde & Matarrese 2009;Giannantonio & Porciani 2010;Baldauf et al. 2011;Tellarini et al. 2015;Assassi et al. 2015). Since we use a fiducial cosmology without PNGs, only leading order (linear) terms in f NL are relevant for the Fisher forecasts and kept in the equations. Commonly the local, equilateral and orthogonal templates are used as proxies for more general PNGs. Those proxies are chosen to test particular aspects of inflation. The discovery of a non-zero PNG of the local shape would be a strong indicator for multifield inflation (Creminelli & Zaldarriaga 2004). A signal of the equilateral shape 2 When projecting with a symmetric window function, the third line of the Z 2 kernel integrates to zero since µ i j k i j = µ i k i + µ j k j = −k l, where 'l' labels the third vector.
arises in a wide range of non-vanilla inflationary dynamics, for instance with non-standard kinetic terms, (see Planck Collaboration et al. (2014) and references therein). Lastly, the orthogonal shape, probes derivative interactions in multifield inflation (Senatore et al. 2010). The mapping of the primordial gravitational potential φ to the late time density contrast δ is done using Poisson's equation together with the matter transfer function T normalized to unity on large scales where the Poisson factor is with the Hubble constant, H 0 , speed-of-light, c, growth factor D and matter density Ω m,0 . The linear matter power spectrum is obtained from the gravitational potential correlators at early times via A non-zero matter bispectrum at early times is proof of PNGs. The discussed proxies for PNGs are translated to late times via A complete bias expansion for isotropic and quadratic PNGs in Lagrangian space can be written using the field Ψ that captures the non-Gaussianities in the primordial potential (Assassi et al. 2015) In the squeezed limit, where the scale dependent bias terms are most relevant, the parameters (A, α) can be determined (Schmidt & Kamionkowski 2010). They are (1,0) for the local, (3,2) for the equilateral and (-3,1) for the orthogonal shape. Translating the expansion into Eulerian space, leads to the following additional terms in the bias expansion (Dalal et al. 2008;McDonald 2008;Verde & Matarrese 2009;Giannantonio & Porciani 2010;Baldauf et al. 2011;Assassi et al. 2015) The N 2 kernel originates from the linear displacement field that maps Lagrangian to Eulerian space (20) and is given by (Tellarini et al. 2015) Using the peak-background split, the two additional bias parameters, (b Ψ , b Ψδ ), can be computed in Lagrangian space in terms of the density bias parameters and the matter variance at the scale R of the tracers (Schmidt & Kamionkowski 2010) We choose for the smoothing scale R of the power spectrum the Lagrangian radius corresponding to tracers of mass M = 10 13 M .
Using the Sheth-Tormen mass function and translating Lagrangian bias parameters to Eulerian space leads to the following expression for the first order bias parameter (Schmidt & Kamionkowski 2010;Desjacques et al. 2011;Schmidt et al. 2013;Karagiannis et al. 2018 while the relevant second order parameter is given by (Giannantonio & Porciani 2010;Karagiannis et al. 2018) b The X serves as a placeholder for the local, equilateral and orthogonal templates. For local PNG, the scale dependent bias is caused by the Poisson factor needed to relate Ψ with the (observed) density contrast δ. For the equilateral and orthogonal PNGs this scale dependence is modified by the σ R,α terms, see (23), in the bias parameters.
The additional bias terms can be included into the Z kernels which become (Tellarini et al. 2015(Tellarini et al. , 2016) For equilateral PNGs, the scale dependent bias becomes constant on large scales and thus degenerate with b 1 . Moreover, its behavior on small scales is degenerate with derivative bias terms and probably unobservable (Assassi et al. 2015). Hence, we exclude it from the power spectrum and bispectrum forecasts. We model the additional stochastic term, Ψ , as Poissonian shot noise which leads to the following non-vanishing power spectrum (Schmidt 2016;Desjacques et al. 2018) All other correlators with the new stochastic term are either zero or higher order and thus discarded. This gives rise to three additional bias terms in the galaxy bispectrum

Observed cross-spectra
Putting all the ingredients together, we can express the matter-galaxy cross-power spectra and bispectra. When projecting, the matter fields will correspond to the lensing convergence that does not suffer from RSDs, accordingly, we use the F 2 kernels. The galaxy fields in contrast, are affected by biasing, RSDs and scale-dependent bias from PNGs and described by the Z i kernels. The leading order power spectra are given by where µ is the cosine of the angle of the wave vector with the LOS. The bispectra are given by B gmm =D B FoG (k 1, ) 2F 2 (k 1 , k 2 )Z 1 (k 1 )P(k 1 )P(k 2 ) + 2F 2 Z 1 (k 1 )P(k 1 )P(k 3 ) B mmm = 2F 2 P(k 1 )P(k 2 ) + 2 perm. + B prim (k 1 , k 2 , k 3 ).
For the sake of compactness, we omitted the arguments of the bispectra. The matter bispectrum depends only on the magnitude of the three wave vectors. RSDs break the statistical isotropy and introduce an explicit dependence on the projection of the wave vector of the galaxy fields on the line of sight.

PROJECTED STATISTICS
In this section we derive the projection integrals for the power spectra and bispectra defined in subsection 2.5. Throughout this paper we use the flat sky approximation and project along the LOS. To simplify calculations, we neglect the time evolution within tomographic bins. Fig. 1 gives an overview of the scenario we have in mind: The CMB lensing kernel determines how the matter between the surface of last scattering and the observer deflects photons from the CMB. In addition, we consider a galaxy survey at low redshift that allows to observe the late time universe directly and to crosscorrelate CMB lensing and galaxy clustering. Fig. 2 shows the composition of galaxy survey window functions in more detail. We keep the survey's volume and position constant and vary the number of tomographic bins within the survey. Beside the Gaussian window that is shown here, we also use a top-hat window since it allows the survey volume to be covered evenly.

Projections
Projections of the matter density contrast along the LOS, δ, with a window functionW, are straightforward in real space: This translates to Fourier space as For galaxy clustering, we assume either Gaussian or Top-hat profiles whose functional forms read in Fourier space Here, the bin size, l, is related to the survey size, d, via l = d/n. The comoving position of the j-th bin center is given by where χ c is the comoving distance to the center of the survey. The Top-hat window prevents volume effects when comparing 2D and 3D analysis or comparing 2D analysis with different numbers of bins (see Fig. 2), but its slow decay for large wave vectors poses numerical challenges. The exponential decay of the Gaussian window in Fourier space in contrast allows us to project the galaxy bispectrum efficiently. The lensing convergence describes the integrated effect of the matter fluctuations between the surface of last scattering surface at χ s and the observer. In an Einstein-deSitter Universe, the scale factor, a, cancels the time evolution of the linear density contrast. The window function for the linear convergence field can then be expressed via a spherical Bessel function of the first kind For numerical reasons, we split the projection integral (37) always into a time-dependent window function and a time-independent statistical field. This implies a different window function for each perturbative order of the convergence. In Appendix A we outline how to project higher orders of the lensing convergence in this framework.

Projected power spectrum
The kernel of projected density fields is uniquely specified by X i ∈ (δ, κ) and its positions χ i . Utilizing the symmetry of the window functions 3 makes cross-power spectra dependent on the two kernels involved and the distance between the centers of the two kernels. The projected power spectrum reads The apostrophe used after the correlator represents a suppressed factor of (2π) 2 δ D (k ⊥ + k ⊥ ), where δ D is the Dirac delta. With n tomographic bins, there are n(n + 1)/2 different (cross-)galaxy power spectra, n galaxy-CMB lensing spectra and one lensinglensing power spectrum.  . Linear power spectrum (blue curve) and its 2D projections with a Gaussian kernel of varying depth (orange to brown). The projections change both the overall amplitude and the power spectrum's shape at low k. Instead of the k ns scaling at low k, the projected power spectra become constant.
The impact of the projection depth on the linear galaxy autopower spectrum (i.e. χ 1 = χ 2 ) is illustrated in Fig. 3. When comparing projected power spectra to the linear power spectrum in 3D, power spectra differ by an overall factor that comes from the volume of the window function in Fourier space. Since projections only source power from smaller to larger scales, the projections lead to an enhancement the large scales where the 3D power spectrum is increasing, i.e., for wavenumbers below the wavenumber corresponding to matter-radiation equality and the peak of the 3D power spectrum. Moreover, the strength of the effect is increasing with decreasing projection depth i.e. wider projection kernels. In the large wave number limit, the integrand of the projection integral (42) becomes independent of the power spectrum and P 2D (k) ∝ P(k).
Since the projection of a 3D homogeneous, isotropic Gaussian random field yields a 2D homogeneous, isotropic Gaussian random field, the power spectrum's covariance can be computed as The fundamental wavenumber of a quadratic survey with area A is k f = 2π/A 1/2 and ∆k determines the bin size in k-space. The inclusion of cross-correlation between different tomographic bins is crucial to resolve all the information available in the survey.

Projected bispectrum
The projected cross-bispectrum correlating the density contrast of the tracers From a numerical perspective, the last expression for the bispectrum is advantageous because it consists of three FFTs -whose results can be cached for each (window, wave vector) combination -followed by a 1D integration. In contrast, the full 2D integration is significantly slower and does not allow to use caching. From now on, we will drop both the subscripts ⊥ and . Due to the angular dependence of the bispectrum, the 1D FFTs in the last part of equations (44) converge only for quickly decaying window functions. However, analytical progress can be made and we refer the interested reader to Appendix A for the details of the lensing projection integrals. Fig. 4 compares the different projected bispectra (44) for equilateral and squeezed configurations. We also show the 3D galaxy bispectrum and one sees that the projections have the strongest effect on largest scales.
Following Joachimi et al. (2009), we compute the bispectrum's covariance in two dimensions as

Theoretical uncertainties
Typically, one deals with theoretical uncertainties due to the perturbative nature of the solutions of the evolution equations by means of a hard (possibly time dependent) cut-off, k max , in the analysis.  the insight that the prediction's accuracy is lost gradually with increasing wave vectors and that the theoretical error can be estimated by the next-to-considered order of the perturbative expansion. We follow this proposal and model the theoretical uncertainties by a mean-zero Gaussian Process with Gaussian covariance function of correlation length r = r bao /2. We used half the correlation length than proposed by Baldauf et al. (2016) which is motivated by the observation that the correlation length should be the length scale on which the spectra are roughly constant and not the one on which we observe changes. Since we work in this paper at leading order, we use the parametrisation for the envelope of the one-loop power spectrum from Baldauf et al. (2016), as the scale of uncertainty where z eff is the survey's mean redshift. Gaussian Processes have the handy property that marginalizing over them modifies the regular convariance function as follows Relative importance of the one-loop correction to the linear power spectrum in 3D (blue curve) and for various projected power spectra (orange to brown). The ratio in 3D follows a power law, whereas the ratios of the projected power spectra become constant at low k. At large k, the ratios of the projected power spectra approach the 3D value from above. ties and with theoretical uncertainties of correlation lengths r BAO and r BAO /2. In the latter case, the error bars saturate earlier than in the former case.
Having outlined the formalism in 3D, we need to project the Gaussian Process. For the projection of the power spectrum covariance we refer to subsection 3.2 and recall that the projection of a Gaussian Process remains a Gaussian Process. The key observation for the projection of the theoretical covariance, C th , is that a Gaussian correlation function can be approximated by a finite sum where c min and c max are found empirically and ∆c = (c max − c min )/N. In practice, O(100) terms are sufficient to achieve subpercent precision. The separable representation has two advantages: instead of the 2D projection integral, one can perform 2N one dimensional integrations which is much faster in our setting. Moreover, the 1D integrals can be cached. This changes the scaling of the required integrations from quadratic to linear in (k max /∆k). Similarly, the scaling becomes linear in the number of redshift bins. The above approach to the bispectrum's theoretical uncertainties becomes computationally intractable when projecting. Due to the implicit wavevector ordering in the correlation function, the direct projection is a 4 dimensional integral. As we were not able to speed the computations sufficiently up, we do not use theoretical uncertainties for the projected bispectrum.
To compare the effect of theoretical uncertainties in 2D and 3D, we compare the projected envelope of the one-loop power spectrum with the projected (linear) power spectrum to their 3D counterpart in Fig. 6. While the relative deviation of the one-loop envelope from the linear power spectrum follows a power law in 3D (see (49)), the ratios in 2D deviate at small k and become constant. For large k, the 2D ratios converge towards the 3D ratio from above.

INFERENCE METHODS
In this section we first review Fisher forecasting and then outline our approach for computing parameter shifts due to inaccurate theoretical predictions.

Fisher forecasting
The Cramer-Rao bound provides a lower bound on the statistical error for any unbiased, linear estimator in terms of the inverse of the Fisher Information (matrix) where i, j label the parameters of interest, L is the likelihood function and all quantities are evaluated at the maximum-likelihood point. Assuming a Gaussian likelihood for the power spectrum and bispectrum, the Fisher Information can be calculated as (Tegmark et al. 1997) The theory vector µ contains the spectra of interest and C is the covariance. The derivatives in 3D are computed via finite differences (Smith et al. 2014). Using the product rule and (42), (44) allows us to compute the derivatives in 2D. The (un)marginalized error forecasts are then given by

Parameter estimation
We are interested in parameter shifts due to inaccurate (theoretical) modelling. In this scenario, we fit some theoretical model µ θ to the underlying ground truth µ true . Assuming a Gaussian distribution, this is done by choosing the parameters θ that maximize the following log-likelihood, Since we investigate small biases, we can linearize the theoretical model around the best fit parameters θ * as where In our case we use the ground truth parameters as best fit parameters. The likelihood of the linearized model has an explicit minimum where F is the Fisher Information and b is given by where for each component of b, one partial derivative is taken. (58) allows to compute the biases from using an inaccurate model.

Survey specifications
For our forecasts and bias estimations, we assume a moderate-sized galaxy survey like CMASS-like with n tomographic bins, each of depth l = d/n, where d is the survey depth along the line of sight. Fig. 2 illustrates the setting and the details of the survey are given in Table 1a.

RESULTS
In this section we first demonstrate empirically that one can recover the full 3D Fisher Information in projected surveys using a small enough projection depth. Next, we establish a consistent cut-off, k 2D max , as a function of the number of redshift bins that allows to compare 2D surveys with different projection depths without the need of theoretical uncertainties. We then use this cut-off to analyse the error bar-bias trade off as a function of the FoG model (11,12). We end the section with more optimistic forecasts.
Throughout this section, we use the finding that the crosscovariances between power spectra and bispectra are negligible at large scales (Song et al. 2015;Chan & Blot 2017;Yankelevich & Porciani 2019). When not stated otherwise, we use from section 5.2 on a CMASS-like survey as described in Fig. 1 and Table 1a.

3D-2D equivalence
We test the equivalence between three and two dimensional matter power spectrum analysis (without RSDs) empirically by performing both analyses and comparing the forecasted error bars. We control the sourcing of small-scale information to large scales with theoretical uncertainties. To this end, we use a cubic survey of side length 1000 h −1 Mpc . The corresponding fundamental frequency is roughly four times smaller than the correlation length of the theoretical uncertainties which ensures that the theoretical uncertainties are approximately constant over k-bin's with width of the fundamental frequency. The projections are done with Top-hat window functions of depths 1000/n, where n is the number of bins. Our findings are summarized in Fig. 7. We show the ratio between the error bars in the two dimensional setting and the values the three dimensional analysis as a function of the 2D cut-off scale, k 2D max . For a given number of bins, the error bars saturate in two dimensions, due to the projected theoretical uncertainties. In addition, those values converge to the 3D values as the number of bins increases.

Choosing the cut-off scale for projected spectra
As explained in section 3.4, it is unfeasible to directly implement theoretical uncertainties for the projected bispectrum. Thus we need another approach to control for theoretical uncertainties in the matter predictions. In this work we control these systematics by choosing a cut-off scale, k max , that ensures that all parameter shifts due to inaccurate matter modelling are below 20% of the corresponding error bars.
We estimate the parameter shifts by fitting a linear matter power spectrum to the (Smith et al. 2003) predictions as described in section 4.2. Both models are without RSDs. Fig. 8   0 Figure 7. The ratio of (unmarginalized) error bars in 3D and 2D, σ 3D x /σ 2D x , is shown as a function of k 2D max . We use the saturated value for the 3D uncertainty and refer to Fig. 5 where we studied the saturation in detail. For each parameter (row) and projection depth (column), the 2D-error bar saturates due to theoretical uncertainties. As the number of tomographic bins increases (from left to right), this saturated value approaches the 3D value.
illustrates the monotonic relation between cut-off scale and maximal relative biases for different projection depths. The cut-off is chosen to be the value where the maximal relative bias in the cosmological parameters is closest to 20%. Those values are marked in black in the figure and the numerical values are reported in Table 1b. We want to point out, that the precise values of those cut-offs are specific for the chosen survey specified in Table 1a.
Since the amount of imprecise, small scale information that gets sourced to larger scales increases with decreasing projection depth, we see that the cut-off decreases as the number of tomographic bins increases. The effective cut-off in 3D, where no sourcing happens, lies in between those extremes because of the different k dependence of the 2D and 3D covariance functions. The lensing cut-off is significantly larger for two reasons. Firstly, the kernel is very narrow in Fourier space and secondly, it peaks at early times, where non-linearities are small.
In Fig. 9 we compare the forecasted error bars from the galaxy power spectrum and bispectrum using the chosen cut-offs in two and three dimensions. As the number of tomographic bins increases, one gains information by resolving more of the modes parallel to the line of sight from cross-correlations between the tomographic bins but looses at the same time from the overall decreasing cut-off scale. For both the power spectrum (dashed line) and the bispectrum (dotted) we see an increase in information until ∼ 10 bins when the latter effects overtake and the information decreases again. This approach allows us to recover more than 80% of bias/amplitude parameters and more than 90% of cosmological parameters in a power spectrum analysis. In a pure bispectrum analysis, two thirds of the Fisher Information can be recovered compared to a 3D analysis.  Table 1a). The points with the black border are closest to 0.2 for each configuration and thus used as cut-offs in this paper. Their numerical value can be found in Table 1b. The 3D curve was determined from a cubic survey of side length l = (A · d) 1/3 with the same redshift and shot-noise as the 2D setting. Table 1. (a): Characterisation of the CMASS-like survey we use in this paper. The velocity dispersion parameters were obtained from fits against N -body simulations. (b): cut-offs for different projection depths that ensure all relative biases in the matter predictions are below 20% for a CMASS-like survey described on the right side. The 3D cut-off comes from a cubic survey with the same volume. The values correspond to the black points in Fig. 8.

Signal-to-noise
There are three effects that determine the signal-to-noise (SN) scaling with respect to the number of tomographic bins. 1) As one increases the number of redshift bins, the projection depth decreases and more radial signal is resolved. 2) k 2D max decreases, if at least one galaxy field is involved, and with smaller projection depth, the SN decreases too. 3) The galaxy selection function is changing when using a variable number of Gaussian profiles. Whereas the first two effects are relevant on all scales, the latter effect's size decreases with the number redshift bins, and is negligible from 4 bins onward. This justifies the use of Gaussian bins. Fig. 10 displays the SN for all individual power specta and bispectra in our fiducial cosmology including RSDs. Due to our conservative cut-offs, the galaxy clustering and lensing auto-power spectra have the most SN. The galaxy bispectrum's SN is strongly increasing with the number of bins and nearly reaches the lensing power spectrum SN at its maximum at 10 bins. The cross-spectra tend to be much smaller than those auto-correlations since the over- Signal-to-Noise δ g δ g δ g δ g δ g δ g κ δ g δ g κ κκ κκδ g κκκ Figure 10. Signal-to-noise scaling for different power specta and bispectra as a function of the number of tomographic bins. The dependence has three components: 1) With more tomographic bins, more radial information is resolved.
2) The cut-off, k 2D max , decreases with the number of tomographic bins which in turn decreases the SN. 3) The galaxy selection function is changing when using a variable number of Gaussian profiles. This effect is negligible from four bins onward. lap between the galaxy survey and the lensing kernel is small (see Fig. 1). The lensing bispectrum SN shows small fluctuations since the binning of its (constant) cut-off is changing to ensure it is consistent with the changing galaxy clustering binning.
Utilizing the covariance structure derived in section 3.2, there are analytical expressions for most of the power spectra's SN scalings and we report those in Appendix B. For the bispectrum, no analytical results exist. While for many power spectra the SN is independent of the chosen RSD model (see Appendix B) this prop-  Table 2. erty is approximately true for all the considered power specta and bispectra.

Error -bias tradeoff
In this subsection we investigate the error-bias trade-off in two scenarios. First we assume a ΛCDM cosmology and forecast the parameter uncertainties and shifts due to inaccurate FoG modelling. Next, we perform a similar analysis for the bias and PNG parameters for a fixed cosmology. In both cases, we assume a ground-truth Gaussian FoG model with velocity dispersions specified in Table 1a. We then study inaccurate models by changing the velocity dispersion parameters by −100% (i.e. no FoG modelling), −50%, −10%, 10% and 50%. Since we find that the biases are approximately independent of the sign of the shift in the velocity dispersions, we restrict ourselves to the first three cases and refer to them as model 1-3. The statistical uncertainties that we report together with biases are always the ones obtained from the inaccurate FoG model. In practice however, the uncertainties are nearly independent of the chosen FoG model. Moreover, we only show the results of the following five spectra combinations that we consider most interesting: galaxy power spectrum (δ g δ g ), all power spectra combined (PS), galaxy bispectrum (δ g δ g δ g ), galaxy power spectrum and bispectrum combined (gal) and all power spectra plus bispectra combined (total).

Cosmological forecasts
We fit the three above mentioned FoG models to the fiducial model and report the error bars and relative biases for all parameters. Fig. 11 shows the scaling of the relative biases and error bars for the model where the velocity dispersion is 50% smaller than its fiducial value (model 2). The five combinations can roughly be grouped into Table 2. This table summarises the optimal forecasted relative errors for three FoG models considered. For each spectra we give the number of tomographic bins such that the maximal relative bias is below 20% and the statistical uncertainties are minimal. On the left side, we report the marginalized relative 1σ uncertainties. On the right side, we show the corresponding relative biases. relative 1σ uncertainties [%] relative biases [%] δ g δ g PS δ g δ g δ g gal total δ g δ g PS δ g δ g δ g gal total Model 1:  three sets. The galaxy bispectrum has the least constraining power of the considered spectra and the derived error bars decrease very strongly with the number of tomographic bins until they reach a minimum around 10 bins (projection depth l 60 h −1 Mpc ). For smaller projection depths, the error bars increase again due to the decreasing k 2D max . The galaxy power spectrum and the combination of galaxy power spectrum and bispectrum scales more weakly with the number of resolved bins and have their minimum around 10 bins too. Finally, adding information from galaxy lensing to either the galaxy power spectrum alone or both power spectrum and bispectrum allows us to lower error bars once more. This is partly because CMB-lensing allows us to break the b 1 -σ 8 degeneracy. It does not help to constrain b s 2 better, as the trace-free part of the tidal tensor is rather uncorrelated to the trace of the tidal tensor that corresponds to the amplitude/bias parameters, which CMB-lensing can constrain well. The relative biases remain close to zero for all parameters, with the exception of σ 8 and b 1 , until a projection depth of l ∼ 100 h −1 Mpc , when they start growing quickly. Changing the FoG model, has a marginal impact on the error bars but shifts the curves of the relative biases left (right) if the FoG model becomes more (less) accurate.
In Table 2 we summarize our findings for the optimal error bars (left side) together with the relative biases (right side) conditioned on being smaller than 20%. We observe that the error bar difference across the spectra is significantly larger than the difference within the spectra across FoG models.
The decrease in Fisher Information for projection depths 60 h −1 Mpc implies that even with a velocity dispersion that is off by 10%, one is able to fully recover the available information in two dimensions. In contrast, a misspecification of 50% leads to a ∼ 10% increase in the error for cosmological parameters and more for bias parameters. The worst case scenario of not modelling the FoG effect at all, inflates the error bars by ∼ 20%.
The relative biases tend to be strongest in parameters that affect the amplitude strongly and tend to be positive, since the models considered underestimate the FoG damping.
The 2D contour plots of the best case Fisher Information matrix of model 2 (see Table 2) is displayed in Fig. 12. One sees that all spectra have approximately the same covariance structure and that CMB lensing helps to break the σ 8 -b 1 degeneracy. The positive correlation between b 1 and b 2 is explained as follows: The three bias/amplitude parameters σ 8 , b 1 and b 2 are all pairwise anticorrelated. However, the positive definiteness of the covariance matrix pushes the weakest among them, b 2 -b 1 , to a positive value in the joint analysis.

PNG forecasts
Assuming a known cosmology, we now forecast bias and f NL parameters for the local, equilateral and orthogonal shape. We perform a separate forecast for each template and FoG model and illustrate the dependence on the number of tomographic bins in Fig. 13.
The scale dependent bias in the galaxy power spectrum yields the best constraints for the local shape. Since the survey is roughly four times as wide as deep, the constraints from the power spectrum have no dependency on the projection depth. The error from the bispectrum analysis, in contrast, shows a strong dependence on the projection depths. For few bins, the error is very large but decreases quickly with decreasing projection depths. As with the cosmological parameters, the error bars become worse for more than 10 bins due to the decreasing cut-off. Since the equilateral shape does not lead to scale dependent bias, it can only be constrained from the bispectrum and the power spectra only contribute towards reducing the uncertainty in the bias parameters. The constraints for the orthogonal shape from the power spectrum and bispectrum are of similar order, since the scale dependent bias only scales as 1/k. The relative biases are close to zero for large projection depths and only start playing a role around 100 h −1 Mpc . We report the optimal forecasts with relative biases below 20% in Table 3. We observe that the error bar differences across spectra are larger than across FoG models given a spectra. For the local shape, the latter differences are basically zero, since only the largest scales are relevant. The constraints of the equilateral and orthogonal shape in contrast improve by ∼ 10% when modelling the FoG effect precisely.

Optimistic forecast
We believe an optimistic but still realistic scenario is given by an upper cut-off of k 2D max = 0.15 hMpc −1 and 6 redshift bins. This is in line with the choice of Karagiannis et al. (2018); Yankelevich & Porciani (2019). Agarwal et al. (2020) were able to work with larger cut-offs for the power spectrum by using separate values for the power spectrum and bispectrum and computing parameter shifts (58) relative to the next perturbative order instead of the fully nonlinear prediction. In Table 4 we summarize the error bars for this scenario in our fiducial cosmology with RSDs. Since, the cut-off for CMB lensing has previously been close to 0.15 hMpc −1 , there are no significant improvements there. In contrast, the error bars from the galaxy power spectrum shrink by a factor of 2 and for the bispectrum by a factor of 3 compared to the best case forecasts with the conservative cut-off from Table 1b. The combined error bars from an analysis with all spectra shrink by 30%. Adding a CMB prior leads to dramatic improvement for all parameters that can be constrained from CMB observations. This is expected because only the next generation galaxy surveys will contain a comparable information content to current CMB surveys. The CMB prior is based on Appendix A of Smith et al. (2014). Let us stress that Table 3. Optimal forecasted errors for a bias-f NL analysis with fixed cosmological parameters. We present separate forecasts for the three FoG models considered. For each spectra we give the number of tomographic bins such that the maximal relative bias is below 20% and the statistical uncertainties are minimal. On the left side, we report the marginalized 1σ uncertainties. On the right side, we show the corresponding relative biases.
1σ uncertainties relative biases [%] δ g δ g PS δ g δ g δ g gal total δ g δ g PS δ g δ g δ g gal total Model 1:  Table 4. Forecasts for the relative 1σ uncertainties of cosmological parameters in a CMASS-like survey with k max = 0.15 h −1 Mpc and 6 tomographic bins without (left side) and with (right) CMB prior.
relative 1σ uncertainties relative 1σ uncertainties without CMB prior [%] with CMB prior [%] δ g δ g PS δ g δ g δ g gal total δ g δ g PS δ g δ g δ g gal total CMASS is a galaxy sample from a moderate-sized survey volume and thus constraints from this sample should not be expected to be competitive with Planck constraints. Upcoming surveys will map regions that are 50-100 times larger than CMASS. This is expected to translate into seven to ten times tighter parameter constraints. Table 5 contains our PNG forecasts in the optimistic scenario. Since PNGs are best determined on large scales, adding small scale information decreases the constraints only by ∼ 20%. Since the large volume of future surveys translates into a smaller k f , constraints on local and orthogonal type non-Gaussinity are expected to improve by more than the volume related factor of seven to ten mentioned above. This is due to the scale dependent bias, which dominates on large scales and leads to additional survey volume dependence.

A simple signal compression approach
The computations can be sped up by removing the configurations from the analysis that contain the least information. Since most configurations contributing to the Fisher matrix are cross-correlations between distinct tomographic galaxy bins, we use the following simple implementation of the idea. Given a correlator with two or three galaxy fields, we only include configurations where the maximal distance between the two galaxy bins is less than some threshold. This reduces the scaling of bispectrum configurations from cubic to linear in the number of tomographic bins. Optimizing the threshold, we find that this speeds up computations, with and without RSDs, for 16 bins by a factor of more than 10 while still recovering more than 99% of the Fisher Information in all parameters. For a smaller total number of bins (deeper bins), the sped-ups are smaller. However, it is not possible to improve the bias-error trade-off for any of the considered FoG models. For further details we refer to Appendix C.

CONCLUSION
In this paper we quantified the statistical power of two-and threepoint correlators of projected density fields in constraining cosmological parameters and primordial non-Gaussianity. We investigated the trade-off between statistical errors and biases induced by imperfect modelling of non-linear redshift space distortions, in particular the Finger-of-God effect. We developed an efficient implementation of the projection integrals required to predict projected power spectra and bispectra.
Using a model for theoretical uncertainties, we have shown empirically that one can recover the full 3D Fisher Information in tomographic surveys with sufficiently small projection depths. Using a projection depth of l = 10 h −1 Mpc allows us to recover 99% of the 3D information. The full account for theoretical uncertainties is numerically not feasible for the 2D bispectrum due to the large number of four dimensional projection integrals that are needed for the non-sparse covariance matrix of the theoretical uncertainties. Instead, we control theoretical uncertainties by computing cut-offs that depend on the projection depth. Those cut-offs are chosen such that the maximal relative biases due to inaccuracies in the matter predictions are less than 20%. This approach allows to recover more than 80% of the information in bias/amplitude parameters and more than 90% of the information in cosmological parameters in a power spectrum analysis. In a bispectrum analysis, this approach allows us to recover 70% of the 3D Fisher Information.
Next, we studied the relation between FoG modelling and relative biases in the forecasted parameters for a CMASS-like survey. We found that the resulting biases are independent of whether one over-or underestimates the FoG damping. We found that not modelling the FoG effect inflates error bars by 20% when controlling for biases. A model that underestimates the velocity dispersion by 50% leads to an increase of 10% and a model whose velocity dispersion differs by 10% is able to recover the full information while maintaining all relative biases smaller than 20%.
We performed a similar analysis for PNGs of the local, equilateral and orthogonal shape. Here, the necessity to model the FoG effect depends crucially on whether or not one best constrains the PNGs from the scale dependent bias or from the template. Whereas in the former case, one can recover most of the information without modelling non-linear RSDs, one can improve the error bars by 10% and 20% for the template dominated orthogonal and equilateral shape respectively with an accurate FoG model.
In a more optimistic scenario where one can control systematic uncertainties up to k max = 0.15 h −1 Mpc , one can expect further improvements of more than 100% for spectra that contain galaxy clustering information. This translates into 30% improvements of combined (PS, total) error bars. Let us stress that future surveys would further improve constraints by 700-1000% due to the much larger volumes enabling LSS constraints as tight as those provided by Planck. With these surveys, the projected power spectrum and bispectrum provide a conservative yet powerful analysis toolkit.
Lastly, we studied the impact of dropping cross-bin galaxy correlations. We find that this is not a tool to improve the biaserror trade-off. However, it is possible to reduce the number of clustering bispectrum configurations considerably without losing much Fisher Information. In practice we were able to speed up the most numerically demanding configurations by a factor of more than 10 while losing less than 1% of Fisher Information in all parameters.
Throughout this paper, we made several simplifying assumptions that could be lifted in future work. For instance, one could use the FFTLog algorithm to go beyond the flat sky approximation used here (Assassi et al. 2017). One could also improve the theoretical modelling by taking more orders of the perturbative expansion into account in order to push k max higher and closer to the more optimistic value mentioned above (see section 5.5 and Tables 4, 5). Taruya A., Nishimichi T., Saito S., 2010, Phys. Rev. D, 82, 063522 Tegmark M., Taylor A. N., Heavens A. F., 1997 Tellarini M., Ross A. J., Tasinato G., Wands D., 2015, J. Cosmology Astropart. Phys., 2015, 004 Tellarini M., Ross A. J., Tasinato G., Wands D., 2016, J. Cosmology Astropart. Phys., 2016, 014 Verde L., Matarrese S., 2009, ApJ, 706, L91 Yankelevich V., Porciani C., 2019, MNRAS, 483, 2078 APPENDIX A: CMB LENSING WINDOW FUNCTION As outlined in section 2.5, we include the time dependency of the matter fields into the window functions when projecting (see (37)). While we ignore the time dependence in the thin tomographic bins of galaxy clustering, it cannot be neglected for CMB-lensing. This leads to different projection kernels for each perturbative order of the convergence field. In this work we are interested in the window functions at the first two orders, but the formalism outlined here is fully general.
The lensing window function in real space without the time evolution is given by where χ s is the comoving distance of the surface of last scattering and we centered the coordinate system at χ s /2. In an Einstein-deSitter Universe, the time evolution can be separated from the density contrast at all orders and the n th order perturbation comes with a factor D n /a.
Our strategy is to approximate the time evolution with polynomials since they allow us to analytically integrate the combined window function. We do this by fitting a sixth order polynomial centered around χ s /2 to the growth factor using the least square method and weights 1/D. The sixth order approximation leads to a relative error of less than 10 −4 for the second order. We refer to Fig. 1 to see the behavior of linear and quadratic lensing function.
The Fourier transform of the lensing window function can be expressed in terms of spherical Bessel functions, which is why we review them quickly in subsection A1 before we actually Fourier transform the polynomial real-space approximation to the window function in subsection A2. Finally, due to the slow decay of the window functions, one has to be careful with the one dimensional integral in the separation of the bispectrum (44). We outline in A3 how one can make sense of supposedly divergent one dimensional integrals and separate the bispectrum.

A2 Fourier transform of monomial times CMB lensing window
The Fourier transform of each monomial can be written as where we substituted q = χ s k/2. This equation is then solved in terms of spherical Bessel functions as follows: The window function in Fourier space is then given by a weighted sum of the relevant monomials.

A3 Integral identities for spherical Bessel functions
Having obtained an analytical expression for the lensing window function at all orders, we now discuss the projection integral of the bispectrum (44). Due to the slow decay of the lensing window function, the angular dependency leads to divergent 1D integrals when separating Terms that involve the power spectrum decay quickly enough, so that they can be solved via standard Fast Fourier Transform (FFT) methods. The other lensing terms either suffer from a very slow decay or even diverge. However, using the the Bessel function's integral representation in terms of a Legendre polynomial P l allows to compute Fourier transforms of spherical Bessel functions with polynomial coefficients for n 0 as follows 4 In summary, the relations presented here allow us to perform the Fourier transformations analytically that cannot be solved with the FFT method. This allows to separate the bispectrum projection integral and to use caching which makes it possible to include all clustering-lensing cross-spectra into the analysis.

APPENDIX B: POWER SPECTRUM SIGNAL-TO-NOISE
The analytical results for the power spectrum signal-to-noise (SN) with Gaussian covariance are summarized as follows: The total SN scales quadratically in k max /k f . Moreover, the total SN scales linearly with the number of included auto-spectra assuming all crossspectra are included in the signal too. These two features are independent of shot-noise, redshift-space distortions and only depend on the assumed Gaussian covariance structure.
In this Appendix we first show the scaling with k max /k f for a single auto-spectra in subsection B1 and then demonstrate the scaling with the number of bins in subsection B2 which implies the k max /k f scaling outlined before holds there too.

B1 Single spectra
Assuming a Gaussian covariance structure and using (43), we can compute the (k max /k f ) scaling of the SN of a single auto-spectrum 4 For n=0, that is simply the inverse Fourier transform and one recovers the real space window function.
Here, we start the binning in k-space at the fundamental frequency k f and use ∆k as k-binning. So the relation between the number of kbins and k max is given by n = k max −k f ∆k . Ignoring this discretization effect, we see that SN scales quadratically in k max /k f .
Cross-power spectra P XY in contrast have a different covariance structure Accordingly there is no analytic result and the SN is determined by the ratio P 2 XY /(P X P Y ).

B2 Combining spectra
As different auto-spectra are correlated, their individual signal-tonoises do not simply add up. The cross-covariance of two autospectra has the form (43) C X X,YY (k i , k j ) ∝ δ K i j P 2 XY (k).
Including the appropriate P XY into the analysis allows to remove the correlation between the auto-spectra. We demonstrate this explicitly for two and three different auto-spectra and the general case follows from induction. We start with two galaxy clustering bins. Due to homogeneity, we know that different scales are uncorrelated, so the result from the previous subsection applies, and we only have to investigate the effect of the two bins for one particular k. The scaling with k max can then be derived from above. The signal vector is then given by S = P 0 P 0 P 1 where P 0 is the auto-power spectrum and P 1 the cross-power spectrum. The covariance is given by where we dropped the mode counting prefactor k 2 f 2πk i ∆k . One sees that the P 2 1 term prohibits a simple adding of the SN's from bins one and two. However, adding P 1 decorrelates the two bins effectively and one obtains SN δ g δ g (2 bin) = SC −1 S = 1 = 2 SN δ g δ g (1 bin).
This effect relies not on the two galaxy clustering bins having the same autospectra. Next, we include CMB-lensing κκ = K and the required clustering-lensing cross-correlations ( κδ i = L i ) into the signal vector. S = P 0 P 0 P 1 L 0 L 1 K .
We have already seen the covariance of S = P 0 P 0 P 1 , so it remains to compute the covariance of S 2 = L 0 L 1 K and their cross-covariance. They are given by: In total, we obtain: and again, adding the cross-correlations allows adding the autosignal-to-noises directly SN tot = SC −1 S = 3 SN δ g δ g (1 bin).
This result illustrates the problem of small scale information being projected onto larger scales. Assuming a fixed survey volume, the total SN increases linearly with the number of bins despite the fact that k max is fixed. The infinite information comes from small scales that pollute the large scales increasingly as the kernel shrinks in real space.
In addition, it turns out that adding cross-spectra as signal to an analysis that also contains auto-spectra does not add any information beyond 'decorrelating' the included auto-spectra. I.e. adding the clustering-lensing cross-spectra does not improve the clustering or lensing SN, it only benefits a joint analysis where it breaks the correlation.

APPENDIX C: RESTRICTING MAXIMAL CORRELATION LENGTH
The number of cross-bispectra configuration is prohibitively large if one wants to perform an MCMC fit or estimate empirical covariance matrices from simulations. This is why several authors have studied compression techniques (Fergusson et al. 2012;Schmittfull et al. 2015;Byun et al. 2017;Gualdi et al. 2018). In this Appendix, we reduce the number of galaxy cross-spectra by introducing a maximal correlation length for all spectra that contain two or more galaxy fields ( δ g δ g , δ g δ g κ and δ g δ g δ g ). We define the maximum correlation length as mcl = d/n · (max(χ) − min(χ)) (C1) where d is the depth of the survey, n the total number of tomographic bins, and χ labels the galaxy bins with larger values assigned to bins that are further away from the observer. Restricting the maximal correlation length reduces the power spectrum scaling from n 2 to n · mcl and for the bispectrum n 3 to n · mcl 2 . We summarize the maximal correlation length needed in order to recover at least 99% of the Fisher information of the analysis with all cross-bin correlations in Table C1. We find that the maximum correlation length is independent of the FoG model used. Since all combinations of spectra show roughly the same maximal correlation length needed (see Table C1), we show in Figure C1 only the scaling for two spectra in full detail. One observes that restricting the maximal correlation lengths leads not better error-bias trade-offs than the ones we discussed in subsection 5.4.1. This is due to the fact that biases arise from inaccurate modelling of radial modes, which is precisely the type of information one excludes when decreasing the correlation length.  Figure C1. Error bars (solid lines) and fives times the biases (dashed lines) for joint clustering-lensing spectra as a function of the maximal correlation length (see (C1)) using FoG model 2. The intersection of dashed and solid lines, marks the point where the relative bias becomes 20%. The top two rows correspond to a power spectrum analysis and the bottom two rows show an analysis with all power spectra and bispectra. Each subplot shows the scaling for a different parameter.