Assessing non-linear models for galaxy clustering I: unbiased growth forecasts from multipole expansion

TL;DR

This paper evaluates the TNS model's effectiveness in using mildly non-linear scales of galaxy clustering data to accurately forecast the growth rate of cosmic structure, emphasizing the importance of multipole expansion.

Contribution

It provides a detailed assessment of the TNS model's validity range and demonstrates how multipole expansion improves growth rate constraints in galaxy surveys.

Findings

Including the hexadecapole improves constraints when scales are restricted.

Unbiased growth rate estimates require a limited range of scales for the hexadecapole.

Multipole expansion formalism effectively quantifies constraints and biases.

Abstract

We assess the performance of the Taruya, Nishimichi and Saito (TNS) model for the halo redshift space power spectrum, focusing on utilising mildly non-linear scales to constrain the growth rate of structure f. Using simulations with volume and number density typical of forthcoming Stage IV galaxy surveys, we determine ranges of validity for the model at redshifts z = 0.5 and z = 1. We proceed to perform a Bayesian MCMC analysis utilising the monopole, quadrupole, and hexadecapole spectra, followed by an exploratory Fisher matrix analysis. As previously noted in other forecasts as well as in real data analyses, we find that including the hexadecapole can significantly improve the constraints. However, a restricted range of scales is required for the hexadecapole in order for the growth parameter estimation to remain unbiased, limiting the improvement. We consistently quantify these effects by employing the multipole expansion formalism in both our Fisher and MCMC forecasts.

Figures

Tables

Table 1. Table 1: Number of bins, k max ​ [ h / Mpc ] subscript 𝑘 max delimited-[] ℎ Mpc k_{\rm max}[h/{\rm Mpc}] , and fiducial parameters found by a least χ 2 superscript 𝜒 2 \chi^{2} fit to the PICOLA data.

z	$N_{\mathrm{bins}}$	$k_{\max }$	$b_1$	$b_2$	$N$	${\sigma }_v$
$0.5$	$36$	$0.227$	$1.506$	$0.091$	$-272$	$8.99$
$1$	$44$	$0.276$	$1.897$	$-0.318$	$504$	$8.09$

Table 2. Table 2: Marginalised percent errors (at 1 ​ σ 1 𝜎 1\sigma ) on f 𝑓 f from the MCMC and Fisher analyses at z = 0.5 𝑧 0.5 z=0.5 and z = 1 𝑧 1 z=1 . We utilise the monopole and quadrupole up to the k max subscript 𝑘 max k_{\rm max} given in Table 1 , and the hexadecapole up to a conservative (restricted) k max , 4 subscript 𝑘 max 4 k_{\rm max,4} , as detailed in the main text. For comparison we show the results for the full P ​ ( k , μ ) 𝑃 𝑘 𝜇 P(k,\mu) Fisher matrix case, and warn that this falsely small error corresponds to biased f 𝑓 f estimates.

Analysis	$z=0.5$	$z=1$
MCMC: $P_0+P_2$	$3.6\%$	$3.0\%$
Fisher: $P_0+P_2$	$3.9\%$	$3.1\%$
MCMC: $P_0+P_2+{P_4|}_{\mathrm{restricted}}$	$3.2\%$	$2.6\%$
Fisher: $P_0+P_2+{P_4|}_{\mathrm{restricted}}$	$3.8\%$	$2.9\%$
Fisher: Full $P(k,\mu )$, biased	$2.4\%$	$1.4\%$

Equations

$$P_{(\mathrm{TNS})}^{\mathrm{S}}(k,\mu )=$$

$$\displaystyle\qquad\qquad\qquad\qquad+b_{1}^{3}A(k,\mu)+b_{1}^{4}B(k,\mu)+b_{1}^{2}C(k,\mu)\Big{]},$$

$$D_{\mathrm{FoG}}^{\mathrm{Lor}}(k^2{\mu }^2{\sigma }_v^2)=\frac{1}{1+(k^2{\mu }^2{\sigma }_v^2)/2},$$

$$P_{\ell }^{\mathrm{S}}(k)=\frac{2\ell +1}{2}{\int }_{-1}^{1}d\mu P^{\mathrm{S}}(k,\mu ){\mathcal{P}}_{\ell }(\mu ),$$

$${\chi }_{\mathrm{red}}^2(k_{\max })=\frac{1}{N_{\mathrm{dof}}}\sum _{k=k_{\min }}^{k_{\max }}\sum _{\ell ,{\ell }^{\prime }=0,2}\left[P_{\ell ,\mathrm{data}}^{\mathrm{S}}(k)-P_{\ell ,\mathrm{model}}^{\mathrm{S}}(k)\right]\text{mbox}{Cov}_{\ell ,{\ell }^{\prime }}^{-1}(k)\left[P_{{\ell }^{\prime },\mathrm{data}}^{\mathrm{S}}(k)-P_{{\ell }^{\prime },\mathrm{model}}^{\mathrm{S}}(k)\right],$$

$$F_{ij}=\frac{V_{\mathrm{s}}}{4{\pi }^2}\sum _{\ell ,{\ell }^{\prime }}{\int }_{k_{\min }}^{k_{\max }}dk{k}^2\frac{\partial P_{\ell }^{\mathrm{S}}(k)}{\partial p_i}{\left[{\widetilde{\mathrm{Cov}}}^{\ell {\ell }^{\prime }}(k)\right]}^{-1}\frac{\partial P_{\ell }^{\mathrm{S}}(k)}{\partial p_j},$$

$${\widetilde{\mathrm{Cov}}}^{\ell {\ell }^{\prime }}(k)=\frac{(2\ell +1)\left(2{\ell }^{\prime }+1\right)}{2}\times {\int }_{-1}^{1}d\mu {\mathcal{P}}_{\ell }^{\mathrm{S}}(\mu ){\mathcal{P}}_{{\ell }^{\prime }}^{\mathrm{S}}(\mu ){\left[P(k,\mu )+\frac{1}{\overline{n}}\right]}^2.$$

$$F_{ij}^{(2\mathrm{D})}=\frac{V_{\mathrm{s}}}{4{\pi }^2}{\int }_{-1}^{1}d\mu {\int }_{k_{\min }}^{k_{\max }}dk{k}^2\frac{\partial P(k,\mu {)}^{\mathrm{S}}}{\partial p_i}{\left\{P(k,\mu )+\frac{1}{\overline{n}}\right\}}^{-2}\frac{\partial P(k,\mu {)}^{\mathrm{S}}}{\partial p_j}.$$

$$\Delta p_i=\sqrt{(F^{-1}{)}_{ii}}.$$

$$r(p_i,p_j)=\frac{(F^{-1}{)}_{ij}}{\sqrt{(F^{-1}{)}_{ii}(F^{-1}{)}_{jj}}}.$$

Full text

1 Introduction

Forthcoming large-scale structure (LSS) surveys such as EUCLID111www.euclid-ec.org (Laureijs et al., 2011), WFIRST222https://wfirst.gsfc.nasa.gov/, and the Dark Energy Spectroscopic Instrument (DESI)333www.desi.lbl.gov, are promising to deliver exquisite cosmological measurements and test the laws of gravity with unprecedented precision. The standard cosmological model, $\Lambda$CDM, has been shown to provide an excellent fit to the data from a suite of CMB and LSS surveys (Ade et al., 2016; Anderson et al., 2013; Song et al., 2015; Beutler et al., 2017), assuming that General Relativity is the correct description of gravity across all scales. A plethora of modified gravity models have been proposed, motivated by the possibility of explaining the late time accelerated expansion of the Universe without the need of a cosmological constant (see Clifton et al., 2012, for a comprehensive review). Future LSS measurements with Stage IV spectroscopic galaxy surveys are aiming to measure the logarithmic growth rate of structure $f$, which strongly depends on cosmology and gravity (Guzzo et al., 2008). This will be achieved by probing the redshift space distortion (RSD) signature with the galaxy power spectrum or correlation function (Blake et al., 2011; Reid et al., 2012; Macaulay et al., 2013; Beutler et al., 2014; Gil-Marín et al., 2016a; Simpson et al., 2016).

The volume and number of galaxies probed by Stage IV surveys are very large, meaning that these surveys will not be limited by statistical uncertainties, but rather by our ability to deal with systematic effects and theoretical uncertainties. The latter include non-linear effects on the redshift space galaxy clustering. In the case of RSD measurements, these have to be understood and modelled properly so that we can confidently utilise the data from mildly non-linear scales to improve the constraints on $f$.

A way to attack this challenge is by using perturbation theory based models (Bernardeau et al., 2002; Kaiser, 1987; Scoccimarro, 2004) that can be easily extended to include non-standard theories of gravity and dark energy (see Bose and Koyama, 2016, 2017; Bose et al., 2018a, b, for example). These can be combined with phenomenological ingredients to model non-linear physics (Taruya et al., 2010; Senatore and Zaldarriaga, 2014; de la Bella et al., 2017). Additionally, a model for the galaxy bias is required to relate the dark matter and galaxy distributions .

An important factor, which is the main focus of this paper, is our ability to use these prescriptions to model non-linear structure formation in an unbiased manner, and how this affects the growth of structure parameter estimation with Stage IV surveys.

For this purpose, we consider the commonly used TNS model (Taruya et al., 2010). This model can reproduce the broadband power spectrum including RSD from simulations at linear and mildly non-linear scales (Nishimichi and Taruya, 2011; Taruya et al., 2013; Ishikawa et al., 2014; Zheng and Song, 2016; Gil-Marín et al., 2016a, b; Bose et al., 2017; Bose and Koyama, 2016). The model is combined with the generalised bias prescription of McDonald and Roy (2009). This is very similar to what has been used in part of the BOSS survey data analysis (Beutler et al., 2017). We use a set of COLA simulations (Tassev et al., 2013; Howlett et al., 2015; Valogiannis and Bean, 2017; Winther et al., 2017) to determine a range of validity for the models at redshifts $z=0.5$ and $1$. We then perform a Bayesian MCMC analysis with the goal of assessing the constraining power of Stage IV-like surveys while also confirming the validity of the models across the range of scales we determined previously. Our focus is getting unbiased constraints on $f$ using information from the first three multipoles: the monopole ($P_{0}$), quadrupole ($P_{2}$), and hexadecapole ($P_{4}$) of the power spectrum. We then move on to an exploratory Fisher matrix forecast analysis using the full anisotropic power spectrum, $P(k,\mu)$, in addition to the multipole expansion formalism. Our final estimates for the unbiased measurements on $f$ are determined consistently from the MCMC and Fisher multipole expansion forecasts analyses.

This paper is organised as follows: In section 2 we present the TNS-based biased tracer RSD model we assemble, and the fits to simulations to determine the fiducial nuisance parameters and model range of validity. In section 3 we perform a Bayesian MCMC analysis and present results, followed by a Fisher matrix analysis in section 4. We then present a comparison between Fisher and MCMC using the multipole expansion formalism. In section 5 we summarise our findings and conclude. Our Fisher matrix codes have been made publicly available at https://github.com/Alkistis/GC-Fish-nonlinear.

2 TNS model description and comparison to simulations

2.1 Model description

The TNS model (Taruya et al., 2010) is based on standard Eulerian perturbation theory (SPT), which assumes that the background space-time is homogeneous and isotropic, and that we work within the Newtonian regime at mildly non-linear scales. These scales are far within the horizon but with $\delta,\theta\ll 1$, where $\delta$ and $\theta$ are the density and velocity perturbations respectively. We also assume that gravity is described by general relativity. The explicit model is given by

[TABLE]

where $f$ is the logarithmic growth rate, $\mu$ is the cosine of the angle between $k$ and the line of sight and $P_{g}$ are the 1-loop galaxy power spectra with the bias model of McDonald and Roy (2009) implicitly included 444We assume the local Lagrangian assumption and hence only first and second order bias are considered: $b_{1}$ and $b_{2}$ as well as a stochasticity term $N$., while the $A$, $B$ and $C$ terms are non-linear perturbative corrections arising from the transformation to redshift space. The terms in brackets are all constructed within SPT555We refer the reader to Beutler et al. (2017) and Taruya et al. (2010) for explicit expressions for the galaxy spectra and correction terms., while the prefactor, $D_{\rm FoG}$, is phenomenological. Here we choose a Lorentzian form (Gil-Marin et al., 2012; Percival and White, 2009; Taruya et al., 2010; Peacock and Dodds, 1994)

[TABLE]

where $\sigma_{v}$ is a redshift-dependent free parameter and represents the velocity dispersion of perturbations at cluster scales. We again refer the reader to Beutler et al. (2017); Taruya et al. (2010) for the formulas for the perturbative components of the model, along with the explicit dependency on the independent free bias parameters $\{b_{1},b_{2},N\}$. As was done in Beutler et al. (2017), we ignore all other bias terms under the local Lagrangian assumption (Sheth et al., 2013; Chan et al., 2012; Saito et al., 2014; Baldauf et al., 2012).

We emphasise again that this model is very similar that used in the recent BOSS analysis (Beutler et al., 2017). We choose it to make use of its thorough validation with simulations and mock catalogues. It is also worth noting that it has shown robustness when considering alternative theories of gravity (for example Bose and Koyama, 2016). However, we stress that there are key differences in our model. We choose a Lorentzian rather than a Gaussian damping factor in Equation 2 666This functional form was shown to give a significantly better fit to COLA multipoles in Bose et al. (2019a)., we include the $C(k,\mu)$ term in Equation 1 and use SPT rather than RegPT (Taruya et al., 2012), the latter being the biggest difference. SPT is known to suffer from divergences in the loop expansion at low redshifts (see Carlson et al., 2009, for example), which the re-normalisation scheme of the RegPT approach addresses. Despite this, SPT clearly does well for $z\geq 1$ (Carlson et al., 2009; Osato et al., 2018). Our work finds that it also works well for $z=0.5$ given our derived goodness of fit to the data. This was also suggested in other works (Carlson et al., 2009; Vlah et al., 2015; Bose and Koyama, 2017). Further, in Bose et al. (2019a) the authors also show that SPT provides a better fit to the redshift space halo multipoles at both $z=0.5$ and $z=1$ than when using the RegPT prescription. Given these findings, it would be very interesting to repeat the BOSS data analysis with this variant of the TNS model.

The full set of nuisance parameters in this model is therefore $\{\sigma_{v},b_{1},b_{2},N\}$. In our MCMC and Fisher matrix analyses we will vary these nuisance parameters along with the cosmological parameter of interest, the growth of structure $f$. Our chosen set of parameters is restricted. Perhaps most importantly, we do not include variations of the background “shape” parameters or the Alcock-Paczynski effect. There are two reasons for this. Firstly, our goal is to demonstrate the trade-off between our constraining power on $f$ and the bias on its estimation as a function of the $k$-ranges used from the monopole, quadrupole, and hexadecapole spectra. We then wish to provide a consistent comparison between Fisher matrix and MCMC forecasts. This can be achieved without an extensive set of parameters. Secondly, varying extra parameters in the MCMC is computationally expensive, since all the model components have to be calculated for every sample. Using a restricted set mitigates this problem. However, optimising our code for speed is under development and we hope to present an analysis with all parameters of interest in future work. We also note that while it is customary to present constraints on $(f\sigma_{8})(z)$, and indeed the BOSS analysis uses this parametrisation, other analyses have opted to present constraints on $f$ alone, e.g. Blake et al. (2011). The reason for this choice in a real data analysis is to test if the Planck best-fit model also predicts the observed growth of structure by the galaxy survey.

Finally, we note that great progress has been made since the TNS model was first proposed, in particular in relation to the relationship between bias and redshift space distortions (see Desjacques et al. (2018a) for a review). Despite this, the TNS model with this minimal bias model remains one of the simplest and most economical phenomenological models. For example, a fully consistent treatment of the halo power spectrum in redshift space discussed requires many more free parameters (Desjacques et al., 2018b).

2.2 Comparison to Simulations

This section is dedicated to determining a rough range of validity for Equation 1 as well as fiducial values for the nuisance parameters. To do this, we make use of a set of four Parallel COmoving Lagrangian Acceleration (PICOLA) simulations (Howlett et al., 2015; Winther et al., 2017) of box length $1024\,\mbox{Mpc}/h$ with $1024^{3}$ dark matter particles and a starting redshift $z_{\rm ini}=49$. These are all run within the same $\Lambda$CDM cosmology taken from WMAP9 (Hinshaw et al., 2013): $\Omega_{m}=0.281$, $\Omega_{b}=0.046$, $h=0.697$, $n_{s}=0.971$ and $\sigma_{8}(z=0)=0.844$.

Halo catalogues from these simulations are constructed using the Friends-of-Friends algorithm with a linking length of 0.2-times the mean particle separation. The halo spectra are measured at redshifts of $z=0.5$ and $z=1$ and use all halos above a mass of $M_{\rm min}=4\times 10^{12}\,M_{\odot}$. This corresponds to a number density of $n_{\rm h}=1\times 10^{-3}\,h^{3}/\mbox{Mpc}^{3}$ which is similar to that estimated for Stage IV surveys galaxy number density around the redshifts considered (Amendola et al., 2018).

As is commonly done in real data analyses, the power spectrum is decomposed into its multipoles. The PICOLA multipoles are measured using the distant-observer approximation 777That is, we assume the observer is located at a distance much greater then the box size ($r\gg 1024\,\mbox{Mpc}/h$), so we treat all the lines of sight as parallel to the chosen Cartesian axes of the simulation box. Next, we use an appropriate velocity component ($v_{x},v_{y}$ or $v_{z}$) to disturb the position of a matter particle. and averaged over three line-of-sight directions. We further average over the four PICOLA simulations. We note here that (PI)COLA is an approximate method and should be validated against full N-body to ensure any comparisons to the measurements are meaningful. In Izard et al. (2016) comparisons of COLA with full N-body measurements are made, and sufficient agreement is found in the halo monopole and quadrupole. Furthermore, in Bose et al. (2019b), comparisons of COLA with full N-body are made for the halo monopole, quadrupole and hexadecapole using the same COLA code and simulation specifications used in this work. The authors find the COLA approach to be sufficiently accurate at the scales and redshifts considered here.

On the theoretical side, the multipoles are expressed as

[TABLE]

where $\mathcal{P}_{\ell}(\mu)$ denote the Legendre polynomials and $P^{\rm S}(k,\mu)$ is given by Equation 1. For our fitting analysis, we utilise only the monopole ($\ell=0$) and quadrupole ($\ell=2$) since they include most of the clustering information. The hexadecapole is later considered in section 3, where we perform an MCMC analysis on the PICOLA data.

Our maximum scale in $k$ (the model’s range of validity) will be denoted as $k_{\rm max}$. To determine this we follow the following procedure:

We fix all cosmological parameters including the growth rate $f$ and perform a least $\chi^{2}$ fit to the PICOLA data by varying the model nuisance parameters888This is done through a brute force search of parameter space using the code MG-Copter (Bose and Koyama, 2016).. In our $\chi^{2}$ analysis we use all data bins from $k_{\rm min}=0.006\,h/{\rm Mpc}$ to $k_{\rm max}$, with $k_{\rm max}$ being varied in the range $0.125\,h/{\rm Mpc}\leq k_{\rm max}\leq 0.3\,h/{\rm Mpc}$. 2. 2.

We take the $95\%$ ($2\sigma$) confidence intervals ($2\Delta\chi^{2}_{\rm red}$) on a $\chi^{2}$ distribution with $N_{\rm dof}$ degrees of freedom. Since $N_{\rm dof}$ is large in our analysis the errors are approximately symmetric. 3. 3.

The criterion we use to calculate the final $k_{\rm max}$ for the rest of the analysis is the maximum k-value which gives $\left[\chi^{2}_{\rm red}(k_{\rm max})~{}-~{}2\Delta\chi^{2}_{\rm red}(k_{\rm max})\right]~{}\leq~{}1$.

Roughly, this gives an indication of which range of scales the model is reliable at without biasing cosmological parameter estimates at the required $2\sigma$ level999This procedure is validated in section 3.. The reduced $\chi^{2}$ statistic is given by

[TABLE]

where $\mbox{Cov}_{\ell,\ell^{\prime}}$ is the Gaussian covariance matrix between the different multipoles and $k_{\rm min}=0.006\,h/{\rm Mpc}$. The number of degrees of freedom $N_{\rm dof}$ is given by $N_{\rm dof}=2\times N_{\rm bins}-N_{\rm params}$, where $N_{\rm bins}$ is the number of $k-$bins used in the summation and $N_{\rm params}$ is the number of free parameters in the theoretical model. Here, $N_{\rm params}=4$ for the TNS model of Equation 1. Finally, the bin-width we use is $\Delta k=0.006\,h/{\rm Mpc}$.

We use linear theory for the covariance matrix between the multipoles (see Appendix C of Taruya et al., 2010, for details). This has been shown to reproduce N-body results up to $k\leq 0.3\,h/\mbox{Mpc}$ at $z=1$. A linear covariance also seems to work well at $z=0.5$ up to $k\leq 0.2\,h/\mbox{Mpc}$ at $z=0.5$, as shown very recently in Sugiyama et al. (2019). We note that for analysing real data from Stage IV surveys validated analytical approximations will likely be used alongside numerical covariances constructed using mocks. In the covariance matrix we assume a number density of $n_{h}=1\times 10^{-3}\,h^{3}/\mbox{Mpc}^{3}$ and a survey volume101010This volume corresponds to a survey with sky area $A_{\rm sky}\simeq 16,000\,(14,000)\,\mbox{deg}^{2}$ and bin width $\Delta z=0.2\,(0.1)$ around $z=0.5\,(1)$. of $V_{\rm s}=4\,\mbox{Gpc}^{3}/h^{3}$.

In Figure 1 we show the minimized $\chi^{2}_{\rm red}(k_{\rm max})$ for $z=0.5$ and $z=1$ for the TNS model, with the associated $2\sigma$ error bars. We determine $k_{\rm max}=0.227\,h/{\rm Mpc}$ and $k_{\rm max}=0.276\,h/{\rm Mpc}$ at $z=0.5$ and $z=1$ respectively. The larger $k_{\rm max}$ at $z=1$ is expected due to less non-linear structure formation at higher redshifts. We summarise the best fit nuisance parameters and details of the fit in Table 1. We also plot the best fit TNS multipoles against the PICOLA data in Figure 2. In the bottom panels of Figure 2 we show the residuals, that is the difference in theoretical prediction to simulation measurement divided by the errors coming from the covariance matrix. Linear theory (Kaiser, 1987) is also shown in green as a reference, where we use $b_{1}$ as measured from the simulations and the fiducial value of $f$.

3 MCMC analysis

In this section we perform Bayesian MCMC analyses at redshifts $z=0.5$ and $1$. We model our log-likelihood using Equation 4, and vary the growth rate $f$ and nuisance parameters of the TNS model outlined in section 2.

Our approach has two purposes. The first is to check how biased the $f$ estimates can be at our derived $k_{\rm max}$ chosen in Table 1. The second is to provide estimates for the $f$ constraints using our version of the TNS model and Stage IV-like specifications, and assess the improvement when adding the hexadecapole.

Our results at $z=0.5$ are shown in the top panel of Figure 3. We first utilise the monopole and quadrupole spectra ($P_{0}$,$P_{2}$) only at the range of scales determined in section 2. The Figure shows the TNS model’s recovery of $f$ at $k_{\rm max}=0.227\,h/$Mpc (red contours). We can see that the fiducial value is safely recovered with a $2\sigma$ criterion (same as in section 2), and we have also checked that using a higher $k_{\rm max}$ the estimates become more biased.

We then add the hexadecapole, $P_{4}$. We find that the estimates of the growth rate $f$ is biased if we take the hexadecapole up to the same $k_{\rm max}$ found in Table 1. That is because the TNS model is not flexible enough to account for the hexadecapole at this range of scales. This has also been seen in the BOSS data analysis in Beutler et al. (2017). However, motivated by Figure 2, we can consider the hexadecapole up to a more conservative value, $k_{\rm max,4}=0.129\,h/$Mpc, without biasing the $f$ estimate (blue contours). This is again similar to what has been done in the BOSS data analysis in Beutler et al. (2017). Our MCMC estimates for the fractional error on $f(z=0.5)$ with this process are $3.6\%$ using the monopole and quadrupole, and $3.2\%$ when adding the hexadecapole with the restricted range of scales. We will refer to this case as $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ throughout the paper.

Our results at $z=1$ are shown in the bottom panel of Figure 3. We follow a similar procedure as before, with $k_{\rm max}=0.276\,h/$Mpc for the monopole and quadrupole, and $k_{\rm max,4}=0.05\,h/$Mpc for the hexadecapole in order for the $f$ estimate to remain unbiased at our required $2\sigma$ level. We have checked that a higher $k_{\rm max,4}$ biases estimates of $f$. This $k_{\rm max,4}$ is lower than the one used at $z=0.5$ which may seem counterintuitive. We can explain this in terms of the model’s flexibility. At $k_{\rm max}=0.276\,h/{\rm Mpc}$ the model is already being severely tested and therefore it cannot account for $P_{4}$ to any higher $k_{\rm max,4}$ without becoming biased. Our MCMC estimates for the marginalised fractional error on $f(z=1)$ with this process are $3\%$ using the monopole and quadrupole, and $2.6\%$ when adding the hexadecapole with the restricted range of scales, $P_{0}+P_{2}+P_{4}|_{\rm restricted}$.

Before moving on to the Fisher matrix analysis, it is useful to perform another set of MCMC analyses to get a sense of the tradeoff between bias and constraining power in $f$ depending on $k_{\rm max}$. At redshift $z=0.5$, in addition to the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ presented above we perform two additional MCMC analyses using the full $P(k,\mu)$ (equivalent to $P_{0}+P_{2}+P_{4}$ with equal $k_{\rm max}$ for monopole, quadrupole, and hexadecapole). We first set $k_{\rm max}=0.227\,h/{\rm Mpc}$, followed by a more conservative $k_{\rm max}=0.129\,h/{\rm Mpc}$. The results are shown in Figure 4. For the full $P(k,\mu)$ with $k_{\rm max}=0.227\,h/{\rm Mpc}$ (red dashed line) we find $f=0.520\pm 0.013$, a heavily biased estimate compared to the true value of $f=0.733$ in the PICOLA simulations (thick black dashed-dot line). For the more conservative full $P(k,\mu)$ with $k_{\rm max}=0.129\,h/{\rm Mpc}$ (green dotted line) we see that the true $f$ is within the $1\sigma$ contour, $f=0.698\pm 0.037$. For the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case, the estimate is unbiased at the $2\sigma$ level, $f=0.782\pm 0.025$.

4 Fisher matrix analysis

In this section we are going to perform an exploratory Fisher matrix analysis to forecast constraints on $f$ using the TNS model. Since we wish to compare the Fisher results with the ones from our MCMC analysis, we will use the Fisher matrix formalism written in terms of multipoles. This allows us to choose different ranges for the monopole, quadrupole, and hexadecapole spectra, to mimic the procedure followed in our MCMC analysis in the previous section. We will also perform a comparison between this method with the most commonly used full anisotropic power spectrum method. This allows us to quickly assess how the requirement for unbiased $f$ estimates limits the improvement with respect to the case where the hexadecapole is added assuming the same $k_{\rm max}$ as in the monopole and quadrupole analysis.

Fisher matrix formalism using multipole expansion:

Here we summarise the multipole expansion formalism for Fisher forecasts. We refer the reader to Taruya et al. (2011) for a comprehensive description. The authors of Taruya et al. (2011) use the TNS model equiped with a linear deterministic bias, noting that this assumption is idealistic and a more accurate prescription is essential for realistic forecasts. They then perform a Fisher matrix analysis to investigate the relative contributions of the different multipoles ($P_{0},P_{2}$ and $P_{4}$) if taken at the same $k_{\rm max}$, and compare with the full $P(k,\mu)$ formalism. Hence, an important distinction between the forecasts performed in this work and that of Taruya et al. (2011) is (i) our implementation of a non-linear bias prescription that corresponds to a more accurate theoretical template that could readily be applied to the data (ii) the use of different $k_{\rm max}$ for the different multipoles to ensure unbiased forecasts and (iii) an MCMC analysis alongside the Fisher matrix analysis to consolidate our findings.

In terms of multipoles, the Fisher matrix for a set of parameters $\{p\}$ is given by (Fisher, 1935; Tegmark, 1997; Seo and Eisenstein, 2007; Taruya et al., 2011)

[TABLE]

with $\widetilde{\mathrm{Cov}}^{\ell\ell^{\prime}}$ being the reduced covariance matrix:

[TABLE]

Note that in our analysis, we promote $k_{\rm max}\rightarrow k_{{\rm max},\ell}$ in Equation 5, in order to be able to study the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case. Also, since we want to mimic our assumptions in the MCMC analysis, we use linear covariance, which means that $P(k,\mu)$ in the brackets of Equation 6 is given by the linear formula (Kaiser, 1987).

Fisher matrix formalism using the full (2D) anisotropic power spectrum:

Considering the full power spectrum signal in redshift space, the Fisher matrix becomes (Tegmark, 1997; Seo and Eisenstein, 2007)

[TABLE]

Using the Fisher matrix formalism, the forecasted errors on parameter $p_{i}$, marginalised over all other parameters, are given by the square root of the diagonal of the inverse of the Fisher matrix as

[TABLE]

We are now ready to present our Fisher matrix analysis results. In Figure 5 we compare

the full $P(k,\mu)$ Fisher matrix results111111We have checked that this is equivalent to using the multipole expansion formalism taking $P_{0}$, $P_{2}$, and $P_{4}$ up to the same $k_{\rm max}$ (Taruya et al., 2011). This is expected since higher order multipoles carry little information (Hand et al., 2017). with the results using only $P_{0}$ and $P_{2}$ up to the $k_{\rm max}$ determined in Table 1, and the ones found when adding the hexadecapole with the restricted range of scales in order for the $f$ estimates to remain unbiased ($P_{0}+P_{2}+P_{4}|_{\rm restricted}$).

In both redshifts, $z=0.5$ and $z=1$, it is clear that the full $P(k,\mu)$ treatment gives much better constraints on $f$. More specifically, using the full $P(k,\mu)$ formalism we find a $2.4\%$ constraint on $f$ compared to $3.8\%$ in the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case at $z=0.5$. At $z=1$, we find $1.4\%$ compared to $2.9\%$ in the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case, a factor of $\sim 2$ difference.

4.1 Comparison between Fisher matrix and MCMC analysis

We will now present the comparison between Fisher matrix and MCMC results, in the unbiased, $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case. For this purpose we will focus on comparing the $(f,\sigma_{v})$ contours, since $\sigma_{v}$ is the most correlated nuisance parameter with the cosmological parameter of interest $f$. We include the full range of comparison plots in Appendix A.

The results are shown in Figure 6. At $z=0.5$, we find remarkable agreement between Fisher and MCMC, demonstrating that the multipole expansion formalism method is both reliable and robust. At $z=1$, we have plotted two Fisher ellipses: the dotted one, denoted by $\mathrm{Fisher|_{deg}}$, corresponds to the Fisher matrix results presented in Figure 5. We see that the area of the $(f,\sigma_{v})$ ellipse is larger than the MCMC one.

Inspecting Figure 5 at $z=1$, we see that there is a near-perfect degeneracy between the $(b_{2},N)$ parameters121212A contributing factor to this could be the loop integral involving $b^{2}_{2}$, that asymptotes to a constant at low $k$ (Desjacques et al., 2018a).. We can quantify this using the correlation coefficient $r$ given by

[TABLE]

This characterises parameter degeneracies: $r=0$ means $p_{i}$ and $p_{j}$ are uncorrelated, while $r=\pm 1$ means they are completely (anti)correlated. In the case we are concerned with at $z=1$, we find $r(b_{2},N)=-0.999$. This signals a possible instability in the Fisher matrix that might be responsible for the disagreement with the MCMC. To investigate this, we impose a conservative $\sim 50\%$ prior on the $N$ parameter and rerun the Fisher matrix analysis. This breaks the degeneracy just enough to mitigate the instability, and gives the excellent agreement shown with the solid black line. However, we emphasise that despite this instability, the marginalised constraint on $f$ is in both cases (with and without the prior correction) in very good agreement with the MCMC. This is because it is the $(f,\sigma_{v})$ degeneracy that affects mostly the final, marginalised error on $f$. The careful reader will notice that there are some differences between the Fisher and MCMC results, especially regarding the improvement when including the hexadecapole at $z=0.5$. This seems to give only marginal gains in the constraining power when considering the Fisher matrix, disagreeing with the improvement shown by the MCMC forecast. Given that the MCMC contours are not perfect Gaussian ellipses we would not expect perfect agreement with the Fisher ellipses. The disagreement is worse at $z=0.5$ and the MCMC contours there also look less “Gaussian” than the ones at $z=1$, so this is not very surprising. We note that similar subtleties have been seen before in Fisher and MCMC comparison studies, and we refer the interested reader to Wolz et al. (2012); Hawken et al. (2012). The final, marginalised $f$ constraints from all methods we have used are summarised in Table 2.

5 Summary and conclusions

In this paper we have assessed the performance of the commonly used TNS model in the context of Stage IV galaxy surveys, considering the multipoles of the redshift space halo power spectrum as our observable. We considered two redshifts, $z=0.5$ and $z=1$, and made use of PICOLA simulations to perform maximum likelihood, MCMC, and Fisher matrix analyses. Here we summarise our main results and conclude.

We first determined the TNS model’s ranges of validity using the monopole and quadrupole power spectra, which contain most of the cosmological information. We note that while this approach is appropriate for the purposes and limitations of this work, in an actual data analysis the bias on $f$ needs to be evaluated directly as a function of $k_{\rm max}$. In Bose et al. (2019b), this point is investigated thoroughly by running a large number of MCMC analyses. In the MCMC analysis that followed, we varied the model’s four nuisance parameters and the logarithmic growth rate, $f$.

The analysis using the $k_{\rm max}$ from Table 1 shows a significant degeneracy between $f$ and $\sigma_{v}$ that has also been found previously (Zheng et al., 2017; Bose et al., 2017). The improvement on the TNS constraints at $z=1$ is mainly due to the much higher $k_{\rm max}$ at $z=1$ compared to that at $z=0.5$.

To investigate the possible improvement by adding information from the hexadecapole, we performed two distinct MCMC analyses at $z=0.5$ and $1$. One excludes the hexadecapole using the range of validity found in section 2, and one includes it. It has been shown that adding the hexadecapole up to the same $k_{\rm max}$ as the monopole and quadrupole biases the estimation of $f$ (see also Figure 4). We followed what was done in the BOSS analysis in Beutler et al. (2017), and restricted its range to a maximum $k_{\rm max,4}=0.129\,h/{\rm Mpc}$ for $z=0.5$ and $k_{\rm max,4}=0.05\,h/{\rm Mpc}$ for $z=1$, so that the estimation of $f$ remains unbiased at the required $2\sigma$ level. Our results are summarised in Table 2.

The addition of the hexadecapole with a restricted $k_{\rm max}$ improves the $f$ constraints, but at a much lower level than using the full $P(k,\mu)$ method, which is equivalent to taking $P_{0},P_{2}$ and $P_{4}$ up to the same $k_{\rm max}$ and leads to biased estimates of $f$.

Finally, we performed a comprehensive Fisher matrix analysis to quickly explore the parameter space and test whether we can reproduce our MCMC analysis results. Using the multipole expansion formalism for the Fisher matrix we reached very good agreement with the MCMC, as shown in Table 2. We also compared the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case with the full $P(k,\mu)$ case, showing that the former gives much more conservative error estimates (as well as avoiding bias in the estimate of $f$ itself).

It is useful to try and assess how our forecasted constraints on the growth rate, given the $k_{\rm max}$ cuts, compare to current requirements for Stage IV surveys. We will refer to the Euclid-like forecasts for the growth rate $f$ presented in Amendola et al. (2018) (see Table 4). At $z=1$ with the same survey volume as we use, they find a fractional error on f of $1.4\%$ in their pessimistic scenario, which has assumed a number density $16\%$ lower than ours. Their adopted model is linear Kaiser multiplied by the BAO smearing function proposed in (Eisenstein et al., 2007). They use the full $P(k,\mu)$ with a $k_{\rm max}\sim 0.2\,{\rm h/Mpc}$. Their set of redshift dependent parameters includes the growth rate $f(z)$, the linear bias $b(z)$, the residual shot noise $P_{\rm s}$, the angular diameter distance $D_{A}(z)$, and the Hubble rate $H(z)$, but the latter two are not marginalised over. Instead, they are projected onto the set of parameters they depend on, which helps with the information on the background (shape) parameters they also vary. Since the models, parameter sets, forecasting method, and $k_{\rm max}$ choices are different, it is very difficult to compare directly, but note that we find the same fractional error of $1.4\%$ using the full $P(k,\mu)$, and a factor of 2 larger error using the hexadecapole cut needed to keep $f(z)$ unbiased at the $2\sigma$ level.

In summary, our study reinforces the need for accurate modelling of the non-linear redshift space power spectrum in light of upcoming Stage IV galaxy surveys. It also shows how reliable forecasts can be obtained if the forecast procedure closely follows what is being done in real data analyses. To our knowledge, a forecasting galaxy clustering analysis for Stage-IV surveys using an accurate theoretical template that can be readily applied to data, and the multipole expansion formalism with different, realistic $k_{\rm max}$ limits informed by simulations, has not been performed in the literature. In addition, we stress that even if a fit to simulations is not available, one should still perform conservative forecasts using the multipole expansion formalism with different ranges for the different multipoles.

That being said, the model considered here is partly phenomenological, and ideally one would want a fully perturbative model that consistently models bias and redshift space distortions, and under which one can obtain a handle on the modelling uncertainties in very general cosmological settings. Several studies have been performed in this direction recently (see, for example, Desjacques et al., 2018b, a; Ivanov and Sibiryakov, 2018; Ding et al., 2018; Hand et al., 2017). In the next paper of this series (Bose et al., 2019a) we also push in this direction and perform an extensive investigation using other proposed models like the Effective Theory of Large Scale Structure (see, for example, Perko et al., 2016; de la Bella et al., 2018; Lewandowski et al., 2018), and compare their performance to the TNS model.

Acknowledgments

We would like to thank the two anonymous referees for their very useful comments and suggestions that improved the quality of this manuscript. We are grateful to Hans Winther for providing the simulation data. We thank Florian Beutler for useful discussions. We acknowledge use of open source software (Jones et al., 01; Hunter, 2007; Lewis et al., 2000; McKinney, 2010; Van Der Walt et al., 2011). KM acknowledges support from the UK Science & Technology Facilities Council through grant ST/N000668/1 and from the UK Space Agency through grant ST/N00180X/1. BB acknowledges support from the Swiss National Science Foundation (SNSF) Professorship grant No.170547. AP is a UK Research and Innovation Future Leaders Fellow, grant MR/S016066/1, and also acknowledges support from the UK Science & Technology Facilities Council through grant ST/S000437/1. Numerical computations for this research were done on the Sciama High Performance Compute (HPC) cluster which is supported by the ICG, SEPNet, and the University of Portsmouth.

Appendix A Fisher and MCMC comparison for the $P_{0}+P_{2}+P_{4}|_{\rm restricted}$ case